My textbook "Density Functional Theory of Atoms and Molecules" by Parr and Yang says that any N-representable density is derivable from a single determinantal wavefunction. A density $\rho$ is N-representable if it satisfies $$ \tag{1} \rho(\mathbf r) = \int d\sigma_1\int |\psi(\mathbf x_1, \mathbf x_2, ... \mathbf x_N)|^2 d^3\mathbf x_2 ... d^3\mathbf x_N $$ where $\psi(\mathbf x_1, \mathbf x_2, ... \mathbf x_N)$ is antisymmetric in exchanging any pair of its arguments, and $x \equiv (\mathbf r, \sigma)$. An antisymmetric wavefunction is in general a linear combination of Slater determinants. What I want to prove is that the above $\rho$ can also be written as $$ \tag{2} \rho(\mathbf r) = \sum_{\sigma=1}^2 \sum_{i=1}^{N_\sigma} |\phi_{i\sigma}(\mathbf r)|^2 $$ for some set of wavefunctions $\{\phi_i\}$. But I don't know where to go beyond this point.
Can someone help me to prove this statement?
EDIT: The Kohn-Sham DFT apparently unconditionally assumes that there exists a non-interacting system having ground state density that is identical to the exact ground state density of the original interacting system. Therefore, stated in a different way, my question can also be understood to ask about the existence of this reference system, does it really always exist for an arbitrary real electronic system? If it doesn't always exist, can one also specify which conditions does the original interacting system have to have in order for its non-interacting "twin" system to exist.