Why does Hartree-Fock (HF) theory even work?

Question

Let’s say we have $N$ electrons and we want to derive the Hartree-Fock (HF) equations. The first step would be to define a Slater determinant of $N$ electrons:

$$\psi(x_1,x_2,… x_N) = \frac{1}{\sqrt{N!}}\begin{vmatrix}\phi_{1}(x_1) & \phi_{2}(x_1) & … & \phi_{N}(x_1)\\ .\\ .\\ .\\ \phi_{1}(x_N) & \phi_{2}(x_N) & … & \phi_{N}(x_N)\end{vmatrix} $$

then we would use the Lagrange minimisation principle to get our HF equations, which are:

$$f\phi_k = \varepsilon_k\phi_k \qquad \forall k=1,\dots,N.$$

$$f := h + \sum_{n=1}^N J_n - K_n,$$

We note the following: We have $N$ electrons in our system so we get a $N\times N$ slater determinant and $N$ molecular wave functions $ \phi_k$.

If we would try to solve this HF equation we could simply put $N$ trial orbitals into the HF equation and solve it iteratively. In practice however, one approximates the molecular orbitals by a linear combination of e.g. $M>N$ basis functions: $$\phi_k=\sum_{m=1}^M{C_{mk} \xi_m}$$

If we put this into the HF equation above, we eventually get a matrix equation (Roothan equations).$$FC=\epsilon SC $$ which can be solved on the computer.

My question is:

There are $M>N$ expansion coefficients $C_m$. By using these $M$ coefficients we eventually get $2M$ molecular orbitals $\phi_k$ with $k=1,\cdots,2M$. Or put in words: we get two molecular orbitals for every basis function we use (due to spin). The lowest $N$ orbitals are the occupied orbitals the highest $2M–N$ are the virtual orbitals. These $2M$ molecular orbitals now correspond to a slater determinant with $2M$ electrons

$$\psi(x_1,x_2,\cdots, x_{2M}) = \frac{1}{\sqrt{2M!}}\begin{vmatrix}\phi_{1}(x_1) & \phi_{2}(x_1) & … & \phi_{2M}(x_1)\\ .\\ .\\ .\\ \phi_{1}(x_{2M}) & \phi_{2}(x_{2M}) & … & \phi_{2M}(x_{2M})\end{vmatrix} $$

But now we have a problem: We wanted to describe a $N$-electron system and stared out by assuming a $N$-electron slater determinant. Now we used this expansion into basis functions and ended up with a slater determinant describing $ 2M>N$ electrons. Isn’t this somehow unphysical? How do we know that our result does even describe a $N$-electron system in the right way? Since our Slater determinant has a size of $2M$ how is it possible to just fill the lowest orbitals and ignore all other?

Answer 1

I think the confusion here is about the number of orbitals vs. the number of electrons in the Slater determinant. In the first equation the indices of $x_i$ enumerate the electrons, but the indices of $\phi_j$ are actually shortcuts for the orbitals: $N$ electrons occupy $N$ different states, but it doesn't mean that there are only $N$ states - we simply use only these $N$ state for constructing the particular $N$-electron wave function.

A clearer but more cumbersome notation would be $$\psi_{j_1, j_2,..., j_N}(x_1,x_2,… x_N) = \frac{1}{\sqrt{N!}}\begin{vmatrix} \phi_{j_1}(x_1) & \phi_{j_2}(x_1) & … & \phi_{j_N}(x_1)\\ .\\ .\\ .\\ \phi_{j_1}(x_N) & \phi_{j_2}(x_N) & … & \phi_{j_N}(x_N)\end{vmatrix}, $$ where $j_1,...,j_N$ is a selection of $N$ orbitals out of $M$ single-particle eigenstates of the one-particle Hamiltonian (i.e., any of the $j_i$ can take values in the range $1...M$, but of course they should be all different to give a non-zero Slater determinant.)

Answer 2

Your expansion coefficient matrix $C$ tells you how to transform your original set of unoptimized basis functions $\xi$ to a set of optimized orbitals $\phi$. The sets have the same size, i.e. if you start out with $M$ basis functions you will also end up with $M$ orbitals.

I think the point that you are missing is that the number $M$ may be larger than the actual number of orbitals that we need to construct our Slater determinant, let that be $N$. We will only use the $N$ lowest optimized orbitals (lowest with respect to the orbital energy i.e. the associated eigenvalue) to construct the Slater determinant. These orbitals will be the $N$ occupied orbitals, which leaves you with $M-N$ unoccupied orbitals that are not used. The number $M$ can be anything as long as it is equal or greater than $N$. The size of the determinant is independent of $M$.

Answer to the remark in comments:

There are in principle $N$ 1-particle functions that minimize the energy of the $N$-particle wavefunction that has the form of a Slaterdeterminant. To find these $N$ 1-particle functions we expand the Hartree-Fock equations into a basis set. This leads us to the Roothan equations, which we can solve by iteration. The Roothan equations system can have size $M \times M$. And it will yield $M$ 1-particle functions as solutions. The lowest $N$ solutions are the best approximations to the "exact" $N$ 1-particles functions that the Hartree-Fock equations define, to construct the best Slater determinant.

So there are exactly $N$ 1-particle functions which we do not know and I would consider these $N$ functions to be the solutions of the Hartree Fock equations. We can find good approximations to these $N$ 1-particle functions by using a basis set expansion which yields the Roothan equations. And once we have these we can construct the $N$-particle wavefunction in Slater determinant form.