In order to solve the equation of an electron in a spherical potential, Hartree first introduced atomic units to eliminate physical constants. Then he converted the Laplacian from Cartesian to spherical coordinates to show that the solution was a product of a radial function ${\displaystyle P(r)/r}$  and a spherical harmonic with an angular quantum number ${\displaystyle \ell }$ , namely ${\displaystyle \psi =(1/r)P(r)S_{\ell }(\theta ,\phi )}$ . The equation for the radial function was^[2][3][4]

${\displaystyle {\frac {\mathrm {d} ^{2}P(r)}{\mathrm {d} r^{2}}}+\left\{2[E-v(r)]-{\frac {\ell (\ell +1)}{r^{2}}}\right\}P(r)=0.}$ 

In mathematics, the Hartree equation, named after Douglas Hartree, is

${\displaystyle i\,\partial _{t}u+\nabla ^{2}u=V(u)u}$ 

in ${\displaystyle \mathbb {R} ^{d+1}}$  where

${\displaystyle V(u)=\pm |x|^{-n}*|u|^{2}}$ 

and

${\displaystyle 0<n<d}$ 

The non-linear Schrödinger equation is in some sense a limiting case.

The wavefunction which describes all of the electrons, ${\displaystyle \Psi }$ , is almost always too complex to calculate directly. Hartree's original method was to first calculate the solutions to Schrödinger's equation for individual electrons 1, 2, 3, ${\displaystyle ...}$ , p, in the states ${\displaystyle \alpha ,\beta ,\gamma ,...,\pi }$ , which yields individual solutions: ${\displaystyle \psi _{\alpha }(\mathbf {x} _{1}),\psi _{\beta }(\mathbf {x} _{2}),\psi _{\gamma }(\mathbf {x} _{3}),...,\psi _{\pi }(\mathbf {x} _{p})}$ . Since each ${\displaystyle \psi }$  is a solution to the Schrödinger equation by itself, their product should at least approximate a solution. This simple method of combining the wavefunctions of the individual electrons is known as the Hartree product:^[5]

${\displaystyle \Psi (\mathbf {x} _{1},\mathbf {x} _{2},\mathbf {x} _{3},...,\mathbf {x} _{p})=\psi _{\alpha }(\mathbf {x} _{1})\psi _{\beta }(\mathbf {x} _{2})\psi _{\gamma }(\mathbf {x} _{3})...\psi _{\pi }(\mathbf {x} _{p})}$ 

This Hartree product gives us the wavefunction of a system (many-particle) as a combination of wavefunctions of the individual particles. It is inherently mean-field (assumes the particles are independent) and is the unsymmetrized version of the Slater determinant ansatz in the Hartree–Fock method. Although it has the advantage of simplicity, the Hartree product is not satisfactory for fermions, such as electrons, because the resulting wave function is not antisymmetric. An antisymmetric wave function can be mathematically described using the Slater determinant.

Let's start from a Hamiltonian of one atom with Z electrons. The same method with some modifications can be expanded to a monoatomic crystal using the Born–von Karman boundary condition and to a crystal with a basis.

${\displaystyle {\hat {H}}=-{\frac {\hbar ^{2}}{2m}}\sum _{i}{\nabla ^{2}}_{\mathbf {r} _{i}}-\sum _{i}{\frac {Ze^{2}}{4\pi \epsilon _{0}|\mathbf {r} _{i}|}}+{\frac {1}{2}}\sum _{i\neq j}{\frac {e^{2}}{4\pi \epsilon _{0}|\mathbf {r} _{i}-\mathbf {r} _{j}|}}}$ 

The expectation value is given by

${\displaystyle \langle \psi |{\hat {H}}|\psi \rangle =\int \psi ^{*}(\mathbf {r} _{1},s_{1},...,\mathbf {r} _{Z},s_{Z}){\hat {H}}\psi (\mathbf {r} _{1},s_{1},...,\mathbf {r} _{Z},s_{Z})\prod _{i}d\mathbf {r} _{i}}$ 

Where the ${\displaystyle s_{i}}$  are the spins of the different particles. In general we approximate this potential with a mean field which is also unknown and needs to be found together with the eigenfunctions of the problem. We will also neglect all relativistic effects like spin-orbit and spin-spin interactions.

At the time of Hartree the full Pauli exclusion principle was not yet invented, it was only clear the exclusion principle in terms of quantum numbers but it was not clear that the wave function of electrons shall be anti-symmetric. If we start from the assumption that the wave functions of each electron are independent we can assume that the total wave function is the product of the single wave functions and that the total charge density at position ${\displaystyle \mathbf {r} }$  due to all electrons except i is

${\displaystyle \rho (\mathbf {r} )=-e\sum _{i\neq j}|\phi _{n_{j}}(\mathbf {r} )|^{2}}$ 

Where we neglected the spin here for simplicity.

This charge density creates an extra mean potential:

${\displaystyle \nabla ^{2}V(\mathbf {r} )=-{\frac {\rho (\mathbf {r} )}{\epsilon _{0}}}}$ 

The solution can be written as the Coulomb integral

${\displaystyle V(\mathbf {r} )={\frac {1}{4\pi \epsilon _{0}}}\int {\frac {\rho (\mathbf {r'} )}{|\mathbf {r} -\mathbf {r'} |}}d\mathbf {r'} =-{\frac {e}{4\pi \epsilon _{0}}}\sum _{i\neq j}\int {\frac {|\phi _{n_{j}}(\mathbf {r'} )|^{2}}{|\mathbf {r} -\mathbf {r'} |}}d\mathbf {r'} }$ 

If we now consider the electron i this will also satisfy the time independent Schrödinger equation

${\displaystyle \left[-{\frac {\hbar \nabla ^{2}}{2m}}-{\frac {Ze^{2}}{4\pi \epsilon _{0}|\mathbf {r} |}}-eV(\mathbf {r} )\right]\phi _{n_{i}}=\mathrm {E} _{i}\phi _{n_{i}}}$ 

This is interesting on its own because it can be compared with a single particle problem in a continuous medium where the dielectric constant is given by:

${\displaystyle \varepsilon (\mathbf {r} )={\frac {\epsilon _{0}}{1+{\frac {4\pi \epsilon _{0}}{Ze}}|\mathbf {r} |V(\mathbf {r} )}}}$ 

Where ${\displaystyle V(\mathbf {r} )<0}$  and ${\displaystyle \varepsilon (\mathbf {r} )>\epsilon _{0}}$ 

Finally, we have the system of Hartree equations

${\displaystyle \left[-{\frac {\hbar \nabla ^{2}}{2m}}-{\frac {Ze^{2}}{4\pi \epsilon _{0}|\mathbf {r} |}}+{\frac {e^{2}}{4\pi \epsilon _{0}}}\sum _{i\neq j}\int {\frac {|\phi _{n_{j}}(\mathbf {r'} )|^{2}}{|\mathbf {r} -\mathbf {r'} |}}d\mathbf {r'} \right]\phi _{n_{i}}=\mathrm {E} _{i}\phi _{n_{i}}}$ 

This is a non linear system of integro-differential equations, but it is interesting in a computational setting because we can solve them iteratively.

Namely, we start from a set of known eigenfunctions (which in this simplified mono-atomic example can be the ones of the hydrogen atom) and starting initially from the potential ${\displaystyle V(\mathbf {r} )=0}$  computing at each iteration a new version of the potential from the charge density above and then a new version of the eigen-functions, ideally these iterations converge.

From the convergence of the potential we can say that we have a "self consistent" mean field, i.e. a continuous variation from a known potential with known solutions to an averaged mean field potential. In that sense the potential is consistent and not so different from the originally used one as ansatz.

In 1928 J. C. Slater and J. A. Gaunt independently showed that given the Hartree product approximation:

${\displaystyle \psi (\mathbf {r} _{1},s_{1},...,\mathbf {r} _{Z},s_{Z})=\prod _{i}^{Z}\phi _{n_{i}}(\mathbf {r} _{i},s_{i})}$ 

They started from the following variational condition

${\displaystyle \delta \left(\langle \prod _{i}\phi _{n_{i}}(\mathbf {r} _{i},s_{i})|{\hat {H}}|\prod _{i}\phi _{n_{i}}(\mathbf {r} _{i},s_{i})\rangle -\sum _{i}\epsilon _{i}\langle \phi _{n_{i}}(\mathbf {r} _{i},s_{i})|\phi _{n_{i}}(\mathbf {r} _{i},s_{i})\rangle \right)=0}$ 

where the ${\displaystyle \epsilon _{i}}$  are the Lagrange multipliers needed in order to minimize the functional of the mean energy ${\displaystyle \langle \psi |{\hat {H}}|\psi \rangle }$ . The orthogonal conditions acts as constraints in the scope of the lagrange multipliers. From here they managed to derive the Hartree equations.

In 1930 Fock and Slater independently then used the slater determinant instead of the Hartree product for the wave function

${\displaystyle \psi (\mathbf {r} _{1},s_{1},...,\mathbf {r} _{Z},s_{Z})={\frac {1}{\sqrt {Z!}}}det{\begin{bmatrix}\phi _{n_{1}}(\mathbf {r} _{1},s_{1})&\phi _{n_{1}}(\mathbf {r} _{2},s_{2})&...&\phi _{n_{1}}(\mathbf {r} _{Z},s_{Z})\\\phi _{n_{2}}(\mathbf {r} _{1},s_{1})&\phi _{n_{2}}(\mathbf {r} _{2},s_{2})&...&\phi _{n_{2}}(\mathbf {r} _{Z},s_{Z})\\...&...&...&...\\\phi _{n_{Z}}(\mathbf {r} _{1},s_{1})&\phi _{n_{Z}}(\mathbf {r} _{2},s_{2})&...&\phi _{n_{Z}}(\mathbf {r} _{Z},s_{Z})\end{bmatrix}}}$ 

This determinant guarantees the exchange symmetry (i.e. if the two columns are swapped the determinant change sign) and the Pauli principle if two electronic states are identical there are two identical rows and therefore the determinant is zero.

They then applied the same variational condition as above

${\displaystyle \delta \left(\langle \psi (\mathbf {r} _{i},s_{i})|{\hat {H}}|\psi (\mathbf {r} _{i},s_{i})\rangle -\sum _{i}\epsilon _{i}\langle \phi _{n_{i}}(\mathbf {r} _{i},s_{i})|\phi _{n_{i}}(\mathbf {r} _{i},s_{i})\rangle \right)=0}$ 

Where now the ${\displaystyle \phi _{n_{i}}}$  are a generic orthogonal set of eigen-functions ${\displaystyle \langle \phi _{n_{i}}(\mathbf {r} ,s_{i})|\phi _{n_{j}}(\mathbf {r} ,s_{j})\rangle =\delta _{ij}}$  from which the wave function is built. The orthogonal conditions acts as constraints in the scope of the lagrange multipliers. From this they derived the Hartree–Fock method.