He's taking the many-electron wavefunction, squaring it, and integrating over all but one of the coordinates.
I.e.,
[tex]
\rho(\vec r) \equiv N \int d^3r_2 d^3r_3 d^3r_4\ldots d^3r_N |\Psi(\vec r,\vec r_2,\vec r_3,\vec r_4,\ldots,\vec r_N)|^2
[/tex]
it desn't matter which one is singled out because the wavefunction is either symmetric or antisymmetric in all it's coordinates and thus the squared wavefunction is symmetric in it's coordinates. multiplication by N is because we want rho normalized such that
[tex]\int d^3r \rho = N[/tex]
rather than normalized to 1.