Is Hartree-Fock (HF) a ground state theory?

Question

To calculate the orbitals in Hartree-Fock (HF) theory we imply the variations principle, so we try to find the wavefunction which minimises the energy $\left<\psi|H|\psi\right>$. This variational principle does only apply to the ground state of a system, so it is said that also HF theory is a ground state theory and can not describe excites states.

But is this really true? If we do a HF calculation we usually get a set of virtual unoccupied orbitals. We could take an electron from an occupied orbital and put it in a virtual orbital. Wouldn’t this correspond to an excited state (at least approximately)? Something like this is often done in introductory courses when discussing the H2-molecule: One says that if both electrons are in the lowest orbital (sigma bonding orbital) it is the ground state. But one could put one electron in the higher orbital (sigma antibodies orbital) and this is often referred to as excited state.

My question in short:

Can HF describe excited states i.e. do the virtual unoccupied orbitals correspond to excited states? If yes, why is this so (the variational principle does not apply to excited states)? If no, what is the interpretation of these virtual orbitals?

Comment: Is there some connection to Koopmans Theorem which says that the energy of a Hartree Fock molecular orbital corresponds (approximately) to the ionisation energy /electron affinity.

Answer 1

Hartree-Fock is a groundstate method. The equations are constructed to optimize the groundstate energy and you will not obtain reasonable excited state energies by replacing an occupied orbital with an unoccupied one in your Slater determinant. This yields a very poor approximation for an excited state. Virtual orbitals do not correspond to excited states in such a simple manner.

But the results from a Hartree-Fock calculation can be used in Post-Hartree-Fock Methods to obtain excited states.

The virtual orbital eigenvalues can be interpreted as electron attachment energy/electron affinity. By looking at the energy of $E(N) - E(N+1)=-\epsilon_v$ where N stands for the number of electrons. $E(N)$ stands for the energy of your groundstate calculated within HF with $N$ occupied orbitals. $E(N+1)$ stands for the energy of a determinant that you obtain by occupying one additional virtual orbital with index $v$ while keeping the former $N$ occupied orbitals as they were. The difference of the energies of these two determinants corresponds to the electron affinity of the unoccupied orbital and is given by its negative orbital energy eigenvalue $\epsilon_v$. This is basically the same as Koopmans theorem but now for electron affinity instead of ionization potential.

But the energies calculated using Koopmans theorem are usually not good approximations. The determinant energies with $N\pm 1$ are bad approximations, when we keep using orbitals that were optimized for a $N$ electron system. If you did a HF calculation starting witn $N+1$ electrons, you would get different orbitals since your effective potential that is part of the Fock Operator changes with the number of electrons. The orbitals of a $N$ and a $N\pm1$ calculation are often similar but not the same. This difference is also often called orbital relaxation when we take the $N$ to $N-1$ case. Since the electrons feel the repulsion of one electron less and thus can "relax" compared to the system with one electron more. And to get these "relaxed" orbitals you would need to do new HF calculation with just $N-1$ electrons, instead of using the orbitals from the $N$ calculation.

Answer 2

From a practical point of view, HF is (almost) always used for ground state calculations, as was nicely explained by Hans Wurst. However, within certain limits, HF calculations can be performed for excited states as well, as I will show below.

The most important point is that the variational principle is applicable not only for ground, but excited states as well! We just have to make sure our approximate $n$-th excited state does not collapse to a lower energy one; this is ensured by the parametrization $$ |\Phi_n\rangle=\hat{P}_n|\Theta_n(a_1,...,a_p)\rangle \ , $$ where $a_1,...,a_p$ are the variational parameters in the trial function $|\Theta_n\rangle$, and $\hat{P}_n$ is a projection operator orthogonal to the lowest $n-1$ exact eigenstates of $\hat{H}$ (wiping out the low-energy components of $|\Theta_n\rangle$ and ensuring $\langle\Phi_n|\hat{H}|\Phi_n\rangle\geq E_n$). If we knew these first $n-1$ eigenstates, then the construction of $\hat{P}_n$ would be trivial: $$ \hat{P}_n=\hat{I}-\sum_{k=1}^{n-1}|\Psi_k\rangle\langle\Psi_k| \ , $$ but of course we generally do not know these states. However, if the desired $n$-th state has different symmetry than the lower ones (i.e. they belong to different irreducible representations of some group $G$ satisfying $[\hat{H},\hat{R}]=0$ for $\forall {\hat{R}}\in G$), then one can easily enforce the orthogonality by projecting $|\Theta_n\rangle$ onto the basis of the appropriate irrep. For example, if we used a function of strictly $P$ symmetry for a helium calculation (no other restrictions imposed), then the variational optimization would tend to the lowest lying $P$ state, which happens to be $^3P$ (since the Ansatz does not have $S$ components, it cannot fall to any of the lower energy states). If we wanted the lowest lying $^1P$ state instead ($E(^1P)>E(^3P)$ according to Hund's rule), then we would have to enforce the singlet spin symmetry too, which in a two-electron case amounts to making the (spatial) Ansatz symmetric under exchange of coordinates. And so on. The main point is that single-state variational calculations can be used to approximate the lowest-lying state of any symmetry sector in the spectrum of $\hat{H}$.

This means that $-$ under the restrictions discussed above $-$ HF calculations can also be performed for excited states, if the state of the given symmetry can be represented by a Slater determinant (alternatively, a symmetry-adapted configuration can be used, if you are willing to relax the definition of HF a bit).

Just for fun, let us calculate the HF energy of the lowest-lying unnatural parity $^3P_g$ state of a helium-like system (a triplet $P$ state, which is invariant under $\boldsymbol{r}_1,\boldsymbol{r}_2\rightarrow-\boldsymbol{r}_1,-\boldsymbol{r}_2$). This is a relatively simple example because the calculation boils down to determining a single (radial) function. In the case $M=+1$, $M_S=+1$, the state can be represented by the Slater determinant $$ |\phi_{2p_0}\uparrow,\phi_{2p_{+1}}\uparrow\rangle =\frac{1}{\sqrt{2}}\left(|\phi_{2p_0}\phi_{2p_{+1}}\rangle-|\phi_{2p_{+1}}\phi_{2p_{0}}\rangle\right)\otimes|\uparrow\uparrow\rangle \ , $$ where the orbitals are written as $$ \phi_{2p_{m}}(\boldsymbol{r})=\frac{1}{r}\chi(r)Y_1^{m}(\theta,\phi) $$ for $m=0,+1$. The energy expression $$ E=\langle\phi_{2p_0}\uparrow,\phi_{2p_{+1}}\uparrow|\hat{H}|\phi_{2p_0}\uparrow,\phi_{2p_{+1}}\uparrow\rangle $$ is minimized with respect to the radial function $\chi$. After substituting and carrying out the angular integrations (using the Laplace expansion of $1/r_{12}$), we find $$ E[\chi]=2\int_0^{\infty}\mathrm{d}r\chi^{*}(r)\left[-\frac{1}{2}\frac{\mathrm{d}^2}{\mathrm{d}r^2}+\frac{1}{r^2}-\frac{Z}{r}\right]\chi(r) \\ + \int_0^{\infty}\mathrm{d}r\int_0^{\infty}\mathrm{d}r' |\chi(r)|^2|\chi(r')|^2\left[\frac{1}{r_{>}}-\frac{1}{5}\frac{r_{<}^2}{r_{>}^3}\right] \ , $$ where $r_{>}=\text{max}(r,r')$ and $r_{<}=\text{min}(r,r')$. The centrifugal term for the $p$ orbitals is $1(1+1)/(2r^2)=1/r^2$. The variation of the functional with normalization condition $$ \int_0^{\infty}\mathrm{d}r|\chi(r)|^2=1 $$ leads to the HF equation: $$ \left[ -\frac{1}{2}\frac{\mathrm{d}^2}{\mathrm{d}r^2}+\frac{1}{r^2}-\frac{Z_{\text{eff}}(r)}{r}\right]\chi(r)=\varepsilon_{2p}\chi(r) \ , $$ where $$ Z_{\text{eff}}(r)=Z -r \int_0^{\infty}\mathrm{d}r' |\chi(r')|^2\left[\frac{1}{r_{>}}-\frac{1}{5}\frac{r_{<}^2}{r_{>}^3}\right] \ . $$ This equation can be solved numerically (virtual orbitals could also be generated, if needed), and for $Z=2$, we find $\varepsilon_{2p}\approx-0.216 \, E_h$ and $E_{\text{HF}}\approx-0.701 \, E_h$, in acceptable agreement with the exact energy $E_{\text{exact}}\approx-0.711 \, E_h$.

\begin{array}{l|c|c} & \varepsilon_{2p} \, / \, E_h & E \, / \, E_h \\ \hline \text{PT[1]} & -0.500 & -0.672 \\ \text{coordinate scaling} & -0.349 & -0.699 \\ \text{HF} & -0.216 & -0.701 \\ \text{HF (g.s)} & +0.348 & -0.553 \\ \text{exact} & \text{n/a} & -0.711 \end{array}

For comparison, we included in the table the first-order PT energy (built on the non-interacting system) $$ E^{[1]} = \left[-\frac{1}{4}Z^2+ \frac{21}{128}Z\right]E_h \ , $$ as well as the coordinate-scaled energy (found with optimized exponent hydrogenic orbitals): $$ E_{\text{opt}}= 2\left(-\frac{1}{2}\frac{Z_{\text{opt}}^2}{2^2} \, E_h\right)= -\frac{1}{4}Z_{\text{opt}}^2 \, E_h \ \ \ , \ \ \ Z_{\text{opt}}=Z-\frac{21}{64} \ . $$

Note that the HF energy is a bit better than $E_{\text{opt}}$, owing to the greater variational flexibility. As a warning, we also showcase the energies calculated with the $2p$ orbitals of the ground state HF calculation (denoted as HF (g.s.)). As one can see, one obtains a completely different quasiparticle energy and a wrong total energy, signalling the fact that virtual orbitals of the HF ground state cannot be directly used for excited states.

Of course, the applicability of this technique is limited, and the above considerations do not help if you need higher energy states from a given symmetry species (e.g. excited $^1S$ states of helium). To find such excited states while retaining some of the simplicity of HF, you should use Multi-Configurational SCF (MCSCF) with appropriate weight factors for the balanced optimization of multiple states.