I doubt any answer will be satisfactory. My opinion is that we are still very far from a mathematical justification. If we accept the mathematical foundations of quantum mechanics, and if we make the approximation that the nucleus of the atom is just one heavy thing with $N$ positive charges, then the motion of the $N$ electrons is governed by a linear equation (Schrödinger) in ${\mathbb R}^{3N}$. The unknown is a function $\psi(r^1,\ldots,r^N,t)$ with the property (Pauli exclusion) that it has full skew-symmetry. For instance, $$\psi(r^2,r^1,\ldots,r^N,t)=-\psi(r^1,r^2,\ldots,r^N,t).$$ In practice, we look for steady states $e^{i\omega t}\phi(r^1,r^2,\ldots,r^N)$. Then $\omega$ is the energy level.
Because of the very large space dimension, one cannot perform reliable calculations on computer, when $N$ is larger than a few units. One attempt to simplify the problem has been to postulate that $\phi$ is a Slatter determinant, which means that $$\phi(r^1,r^2,\ldots,r^N)=\|a_i(r^j)\|_{1\le i,j\le N}.$$ The unknown is then an $N$-tuple of functions $a_i$ over ${\mathbb R}^3$. Of course, we do not expect that steady states be really Slater determinants; after all, the Schrödinger equation does not preserve the class of Slater determinants. Thus there is a price to pay, which is to replace the Schrödinger equation by an other one, obtained by an averaging process (Hartree--Fock model). The drawback is that the new equation is non-linear. Such approximate states have been studied by P.-L. Lions & I. Catto in the 90's.
Update. Suppose $N=2$ only. If we think to $\phi$ as a finite-dimensional object instead of an $L^2$-function, then it is nothing but a skew-symmetric matrix $A$. Approximation à la Slater consists in writing $A\sim XY^T-YX^T$, where $X$ and $Y$ are vectors. In other words, one approximate $A$ by a rank-two skew-symmetric matrix. The approximation must be in terms of the Hilbert-Schmidt norm (also named Frobenius, Schur): this norm is natural because of the requirement $\|\phi\|_{L^2}=N$. If $\pm a_1,\ldots,\pm a_m$ are the pairs of eigenvalues of $A$, with $0\le a_1\le\ldots\le a_m$, then the best Slater approximation $B$ satisfies $\|B\|^2=2a_m^2$, $\|A-B\|^2=2(a_1^2+\cdots+a_{m-1}^2)$. Not that good. Imagine how much worse it can be if $N$ is larger than $2$.
I doubt any answer will be satisfactory. My opinion is that we are still very far from a mathematical justification. If we accept the mathematical foundations of quantum mechanics, and if we make the approximation that the nucleus of the atom is just one heavy thing with $N$ positive charges, then the motion of the $N$ electrons is governed by a linear equation (Schrödinger) in ${\mathbb R}^{3N}$. The unknown is a function $\psi(r^1,\ldots,r^N,t)$ with the property (Pauli exclusion) that it has full skew-symmetry. For instance, $$\psi(r^2,r^1,\ldots,r^N,t)=-\psi(r^1,r^2\ldots,r^N,t).$$ $\psi(r^2,r^1,\ldots,r^N,t)=-\psi(r^1,r^2,\ldots,r^N,t).$$In practice, we look for steady states e^{i\omega t}\phi(r^1,r^2\ldots,r^N)t}\phi(r^1,r^2,\ldots,r^N). Then \omega is the energy level. Because of the very large space dimension, one cannot perform reliable calculations on computer, when N is larger than a few units. One attempt to simplify the problem has been to postulate that \phi is a Slatter determinant, which means that$$\phi(r^1,r^2\ldots,r^N)=\|a_i(r^j)\|_{1\le$\phi(r^1,r^2,\ldots,r^N)=\|a_i(r^j)\|_{1\le i,j\le N}.$$The unknown is then an N-tuple of functions a_i over {\mathbb R}^3. Of course, we do not expect that steady states be really Slatter Slater determinants; after all, the Schrödinger equation does not preserve the class of Slatter Slater determinants. Thus there is a price to pay, which is to replace the Schrödinger equation by an other one, obtained by an averaging process (Hartree--Fock model). The drawback is that the new equation is non-linear. Such approximate states have been studied by P.-L. Lions & I. Catto in the 90's. 1 I doubt any answer will be satisfactory. My opinion is that we are still very far from a mathematical justification. If we accept the mathematical foundations of quantum mechanics, and if we make the approximation that the nucleus of the atom is just one heavy thing with N positive charges, then the motion of the N electrons is governed by a linear equation (Schrödinger) in {\mathbb R}^{3N}. The unknown is a function \psi(r^1,\ldots,r^N,t) with the property (Pauli exclusion) that it has full skew-symmetry. For instance,$$\psi(r^2,r^1,\ldots,r^N,t)=-\psi(r^1,r^2\ldots,r^N,t).$$In practice, we look for steady states e^{i\omega t}\phi(r^1,r^2\ldots,r^N). Then \omega is the energy level. Because of the very large space dimension, one cannot perform reliable calculations on computer, when N is larger than a few units. One attempt to simplify the problem has been to postulate that \phi is a Slatter determinant, which means that$$\phi(r^1,r^2\ldots,r^N)=\|a_i(r^j)\|_{1\le i,j\le N}. The unknown is then an $N$-tuple of functions $a_i$ over ${\mathbb R}^3$. Of course, we do not expect that steady states be really Slatter determinants; after all, the Schrödinger equation does not preserve the class of Slatter determinants. Thus there is a price to pay, which is to replace the Schrödinger equation by an other one, obtained by an averaging process. The drawback is that the new equation is non-linear. Such approximate states have been studied by P.-L. Lions & I. Catto in the 90's.