We consider a load balancing model where a Poisson stream of jobs arrive at a system of many servers whose service time distribution possesses a finite second moment. A small fraction of arrivals pass through the so called power-of-choice algorithm, which assigns a job to the shortest among $\ell$, $\ell\ge 2$, randomly chosen queues, and the remaining jobs are assigned to queues chosen uniformly at random. The system is analyzed at critical load in an asymptotic regime where both the number of servers and the usual heavy traffic parameter associated with individual queue lengths grow to infinity. The first main result is a hydrodynamic limit, where the empirical measure of the diffusively normalized queue lengths is shown to converge to a path in measure space whose density is given by the unique solution of a parabolic PDE with nonlocal coefficients. Further, two forms of an invariance principle are proved, corresponding to two different assumptions on the initial distribution, where individual normalized queue lengths converge weakly to solutions of SDE. In one of these results, the limit is given by a McKean-Vlasov SDE, and propagation of chaos holds. The McKean-Vlasov limit is closely related to limit results for Brownian particles on $\mathbb{R}_+$ interacting through their rank (with a specific interaction). However, an entirely different set of tools is required, as the collection of $n$ prelimit particles does not obey a Markovian evolution on $\mathbb{R}_+^n$.