Hartree-Fock methods

Why Hartree-Fock? Derivation of Hartree-Fock equations in coordinate space

Hartree-Fock (HF) theory is an algorithm for finding an approximative expression for the ground state of a given Hamiltonian. The basic ingredients are

• Define a single-particle basis \( \{\psi_{\alpha}\} \) so that

• The term \( \hat{u}^{\mathrm{HF}} \) is a single-particle potential to be determined by the HF algorithm.
• The HF algorithm means to choose \( \hat{u}^{\mathrm{HF}} \) in order to have

• The variational principle ensures that \( E^{\mathrm{HF}} \ge E_0 \), with \( E_0 \) the exact ground state energy.

Variational Calculus and Lagrangian Multipliers

The calculus of variations involves problems where the quantity to be minimized or maximized is an integral.

In the general case we have an integral of the type $$ E[\Phi]= \int_a^b f(\Phi(x),\frac{\partial \Phi}{\partial x},x)dx, $$ where \( E \) is the quantity which is sought minimized or maximized. The problem is that although \( f \) is a function of the variables \( \Phi \), \( \partial \Phi/\partial x \) and \( x \), the exact dependence of \( \Phi \) on \( x \) is not known. This means again that even though the integral has fixed limits \( a \) and \( b \), the path of integration is not known. In our case the unknown quantities are the single-particle wave functions and we wish to choose an integration path which makes the functional \( E[\Phi] \) stationary. This means that we want to find minima, or maxima or saddle points. In physics we search normally for minima. Our task is therefore to find the minimum of \( E[\Phi] \) so that its variation \( \delta E \) is zero subject to specific constraints. In our case the constraints appear as the integral which expresses the orthogonality of the single-particle wave functions. The constraints can be treated via the technique of Lagrangian multipliers

Let us specialize to the expectation value of the energy for one particle in three-dimensions. This expectation value reads $$ E=\int dxdydz \psi^*(x,y,z) \hat{H} \psi(x,y,z), $$ with the constraint $$ \int dxdydz \psi^*(x,y,z) \psi(x,y,z)=1, $$ and a Hamiltonian $$ \hat{H}=-\frac{1}{2}\nabla^2+V(x,y,z). $$ We will, for the sake of notational convenience, skip the variables \( x,y,z \) below, and write for example \( V(x,y,z)=V \).

The integral involving the kinetic energy can be written as, with the function \( \psi \) vanishing strongly for large values of \( x,y,z \) (given here by the limits \( a \) and \( b \)), $$ \int_a^b dxdydz \psi^* \left(-\frac{1}{2}\nabla^2\right) \psi dxdydz = \psi^*\nabla\psi|_a^b+\int_a^b dxdydz\frac{1}{2}\nabla\psi^*\nabla\psi. $$ We will drop the limits \( a \) and \( b \) in the remaining discussion. Inserting this expression into the expectation value for the energy and taking the variational minimum we obtain $$ \delta E = \delta \left\{\int dxdydz\left( \frac{1}{2}\nabla\psi^*\nabla\psi+V\psi^*\psi\right)\right\} = 0. $$

The constraint appears in integral form as $$ \int dxdydz \psi^* \psi=\mathrm{constant}, $$ and multiplying with a Lagrangian multiplier \( \lambda \) and taking the variational minimum we obtain the final variational equation $$ \delta \left\{\int dxdydz\left( \frac{1}{2}\nabla\psi^*\nabla\psi+V\psi^*\psi-\lambda\psi^*\psi\right)\right\} = 0. $$ We introduce the function \( f \) $$ f = \frac{1}{2}\nabla\psi^*\nabla\psi+V\psi^*\psi-\lambda\psi^*\psi= \frac{1}{2}(\psi^*_x\psi_x+\psi^*_y\psi_y+\psi^*_z\psi_z)+V\psi^*\psi-\lambda\psi^*\psi, $$ where we have skipped the dependence on \( x,y,z \) and introduced the shorthand \( \psi_x \), \( \psi_y \) and \( \psi_z \) for the various derivatives.

For \( \psi^* \) the Euler-Lagrange equations yield $$ \frac{\partial f}{\partial \psi^*}- \frac{\partial }{\partial x}\frac{\partial f}{\partial \psi^*_x}-\frac{\partial }{\partial y}\frac{\partial f}{\partial \psi^*_y}-\frac{\partial }{\partial z}\frac{\partial f}{\partial \psi^*_z}=0, $$ which results in $$ -\frac{1}{2}(\psi_{xx}+\psi_{yy}+\psi_{zz})+V\psi=\lambda \psi. $$ We can then identify the Lagrangian multiplier as the energy of the system. The last equation is nothing but the standard Schroedinger equation and the variational approach discussed here provides a powerful method for obtaining approximate solutions of the wave function.

Definitions and notations

Before we proceed we need some definitions. We will assume that the interacting part of the Hamiltonian can be approximated by a two-body interaction. This means that our Hamiltonian is written as the sum of some onebody part and a twobody part $$ \begin{equation} \hat{H} = \hat{H}_0 + \hat{H}_I = \sum_{i=1}^A \hat{h}_0(x_i) + \sum_{i < j}^A \hat{v}(r_{ij}), \tag{1} \end{equation} $$ with $$ \begin{equation} H_0=\sum_{i=1}^A \hat{h}_0(x_i). \tag{2} \end{equation} $$ The onebody part \( u_{\mathrm{ext}}(x_i) \) is normally approximated by a harmonic oscillator potential or the Coulomb interaction an electron feels from the nucleus. However, other potentials are fully possible, such as one derived from the self-consistent solution of the Hartree-Fock equations to be discussed here.

Our Hamiltonian is invariant under the permutation (interchange) of two particles. Since we deal with fermions however, the total wave function is antisymmetric. Let \( \hat{P} \) be an operator which interchanges two particles. Due to the symmetries we have ascribed to our Hamiltonian, this operator commutes with the total Hamiltonian, $$ [\hat{H},\hat{P}] = 0, $$ meaning that \( \Psi_{\lambda}(x_1, x_2, \dots , x_A) \) is an eigenfunction of \( \hat{P} \) as well, that is $$ \hat{P}_{ij}\Psi_{\lambda}(x_1, x_2, \dots,x_i,\dots,x_j,\dots,x_A)= \beta\Psi_{\lambda}(x_1, x_2, \dots,x_i,\dots,x_j,\dots,x_A), $$ where \( \beta \) is the eigenvalue of \( \hat{P} \). We have introduced the suffix \( ij \) in order to indicate that we permute particles \( i \) and \( j \). The Pauli principle tells us that the total wave function for a system of fermions has to be antisymmetric, resulting in the eigenvalue \( \beta = -1 \).

In our case we assume that we can approximate the exact eigenfunction with a Slater determinant $$ \begin{equation} \Phi(x_1, x_2,\dots ,x_A,\alpha,\beta,\dots, \sigma)=\frac{1}{\sqrt{A!}} \left| \begin{array}{ccccc} \psi_{\alpha}(x_1)& \psi_{\alpha}(x_2)& \dots & \dots & \psi_{\alpha}(x_A)\\ \psi_{\beta}(x_1)&\psi_{\beta}(x_2)& \dots & \dots & \psi_{\beta}(x_A)\\ \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots \\ \psi_{\sigma}(x_1)&\psi_{\sigma}(x_2)& \dots & \dots & \psi_{\sigma}(x_A)\end{array} \right|, \tag{3} \end{equation} $$ where \( x_i \) stand for the coordinates and spin values of a particle \( i \) and \( \alpha,\beta,\dots, \gamma \) are quantum numbers needed to describe remaining quantum numbers.

The single-particle function \( \psi_{\alpha}(x_i) \) are eigenfunctions of the onebody Hamiltonian \( h_i \), that is $$ \hat{h}_0(x_i)=\hat{t}(x_i) + \hat{u}_{\mathrm{ext}}(x_i), $$ with eigenvalues $$ \hat{h}_0(x_i) \psi_{\alpha}(x_i)=\left(\hat{t}(x_i) + \hat{u}_{\mathrm{ext}}(x_i)\right)\psi_{\alpha}(x_i)=\varepsilon_{\alpha}\psi_{\alpha}(x_i). $$ The energies \( \varepsilon_{\alpha} \) are the so-called non-interacting single-particle energies, or unperturbed energies. The total energy is in this case the sum over all single-particle energies, if no two-body or more complicated many-body interactions are present.

Let us denote the ground state energy by \( E_0 \). According to the variational principle we have $$ E_0 \le E[\Phi] = \int \Phi^*\hat{H}\Phi d\mathbf{\tau} $$ where \( \Phi \) is a trial function which we assume to be normalized $$ \int \Phi^*\Phi d\mathbf{\tau} = 1, $$ where we have used the shorthand \( d\mathbf{\tau}=dx_1dr_2\dots dr_A \).

In the Hartree-Fock method the trial function is the Slater determinant of Eq. (3) which can be rewritten as $$ \Phi(x_1,x_2,\dots,x_A,\alpha,\beta,\dots,\nu) = \frac{1}{\sqrt{A!}}\sum_{P} (-)^P\hat{P}\psi_{\alpha}(x_1) \psi_{\beta}(x_2)\dots\psi_{\nu}(x_A)=\sqrt{A!}\hat{A}\Phi_H, $$ where we have introduced the antisymmetrization operator \( \hat{A} \) defined by the summation over all possible permutations of two particles.

It is defined as $$ \begin{equation} \hat{A} = \frac{1}{A!}\sum_{p} (-)^p\hat{P}, \tag{4} \end{equation} $$ with \( p \) standing for the number of permutations. We have introduced for later use the so-called Hartree-function, defined by the simple product of all possible single-particle functions $$ \Phi_H(x_1,x_2,\dots,x_A,\alpha,\beta,\dots,\nu) = \psi_{\alpha}(x_1) \psi_{\beta}(x_2)\dots\psi_{\nu}(x_A). $$

Both \( \hat{H}_0 \) and \( \hat{H}_I \) are invariant under all possible permutations of any two particles and hence commute with \( \hat{A} \) $$ \begin{equation} [H_0,\hat{A}] = [H_I,\hat{A}] = 0. \tag{5} \end{equation} $$ Furthermore, \( \hat{A} \) satisfies $$ \begin{equation} \hat{A}^2 = \hat{A}, \tag{6} \end{equation} $$ since every permutation of the Slater determinant reproduces it.

The expectation value of \( \hat{H}_0 \) $$ \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau} = A! \int \Phi_H^*\hat{A}\hat{H}_0\hat{A}\Phi_H d\mathbf{\tau} $$ is readily reduced to $$ \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau} = A! \int \Phi_H^*\hat{H}_0\hat{A}\Phi_H d\mathbf{\tau}, $$ where we have used Eqs. (5) and (6). The next step is to replace the antisymmetrization operator by its definition and to replace \( \hat{H}_0 \) with the sum of one-body operators $$ \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau} = \sum_{i=1}^A \sum_{p} (-)^p\int \Phi_H^*\hat{h}_0\hat{P}\Phi_H d\mathbf{\tau}. $$

The integral vanishes if two or more particles are permuted in only one of the Hartree-functions \( \Phi_H \) because the individual single-particle wave functions are orthogonal. We obtain then $$ \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau}= \sum_{i=1}^A \int \Phi_H^*\hat{h}_0\Phi_H d\mathbf{\tau}. $$ Orthogonality of the single-particle functions allows us to further simplify the integral, and we arrive at the following expression for the expectation values of the sum of one-body Hamiltonians $$ \begin{equation} \int \Phi^*\hat{H}_0\Phi d\mathbf{\tau} = \sum_{\mu=1}^A \int \psi_{\mu}^*(x)\hat{h}_0\psi_{\mu}(x)dx d\mathbf{r}. \tag{7} \end{equation} $$

We introduce the following shorthand for the above integral $$ \langle \mu | \hat{h}_0 | \mu \rangle = \int \psi_{\mu}^*(x)\hat{h}_0\psi_{\mu}(x)dx, $$ and rewrite Eq. (7) as $$ \begin{equation} \int \Phi^*\hat{H}_0\Phi d\tau = \sum_{\mu=1}^A \langle \mu | \hat{h}_0 | \mu \rangle. \tag{8} \end{equation} $$

The expectation value of the two-body part of the Hamiltonian is obtained in a similar manner. We have $$ \int \Phi^*\hat{H}_I\Phi d\mathbf{\tau} = A! \int \Phi_H^*\hat{A}\hat{H}_I\hat{A}\Phi_H d\mathbf{\tau}, $$ which reduces to $$ \int \Phi^*\hat{H}_I\Phi d\mathbf{\tau} = \sum_{i\le j=1}^A \sum_{p} (-)^p\int \Phi_H^*\hat{v}(r_{ij})\hat{P}\Phi_H d\mathbf{\tau}, $$ by following the same arguments as for the one-body Hamiltonian.

Because of the dependence on the inter-particle distance \( r_{ij} \), permutations of any two particles no longer vanish, and we get $$ \int \Phi^*\hat{H}_I\Phi d\mathbf{\tau} = \sum_{i < j=1}^A \int \Phi_H^*\hat{v}(r_{ij})(1-P_{ij})\Phi_H d\mathbf{\tau}. $$ where \( P_{ij} \) is the permutation operator that interchanges particle \( i \) and particle \( j \). Again we use the assumption that the single-particle wave functions are orthogonal.

We obtain $$ \begin{align} \int \Phi^*\hat{H}_I\Phi d\mathbf{\tau} = \frac{1}{2}\sum_{\mu=1}^A\sum_{\nu=1}^A &\left[ \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\mu}(x_i)\psi_{\nu}(x_j) dx_idx_j \right.\\ &\left. - \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j) \hat{v}(r_{ij})\psi_{\nu}(x_i)\psi_{\mu}(x_j) dx_idx_j \right]. \tag{9} \end{align} $$ The first term is the so-called direct term. It is frequently also called the Hartree term, while the second is due to the Pauli principle and is called the exchange term or just the Fock term. The factor \( 1/2 \) is introduced because we now run over all pairs twice.

The last equation allows us to introduce some further definitions. The single-particle wave functions \( \psi_{\mu}(x) \), defined by the quantum numbers \( \mu \) and \( x \) are defined as the overlap $$ \psi_{\alpha}(x) = \langle x | \alpha \rangle . $$

We introduce the following shorthands for the above two integrals $$ \langle \mu\nu|\hat{v}|\mu\nu\rangle = \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\mu}(x_i)\psi_{\nu}(x_j) dx_idx_j, $$ and $$ \langle \mu\nu|\hat{v}|\nu\mu\rangle = \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j) \hat{v}(r_{ij})\psi_{\nu}(x_i)\psi_{\mu}(x_j) dx_idx_j. $$

Derivation of Hartree-Fock equations in coordinate space

Let us denote the ground state energy by \( E_0 \). According to the variational principle we have $$ E_0 \le E[\Phi] = \int \Phi^*\hat{H}\Phi d\mathbf{\tau} $$ where \( \Phi \) is a trial function which we assume to be normalized $$ \int \Phi^*\Phi d\mathbf{\tau} = 1, $$ where we have used the shorthand \( d\mathbf{\tau}=dx_1dx_2\dots dx_A \).

In the Hartree-Fock method the trial function is a Slater determinant which can be rewritten as $$ \Psi(x_1,x_2,\dots,x_A,\alpha,\beta,\dots,\nu) = \frac{1}{\sqrt{A!}}\sum_{P} (-)^PP\psi_{\alpha}(x_1) \psi_{\beta}(x_2)\dots\psi_{\nu}(x_A)=\sqrt{A!}\hat{A}\Phi_H, $$ where we have introduced the anti-symmetrization operator \( \hat{A} \) defined by the summation over all possible permutations p of two fermions. It is defined as $$ \hat{A} = \frac{1}{A!}\sum_{p} (-)^p\hat{P}, $$ with the the Hartree-function given by the simple product of all possible single-particle function $$ \Phi_H(x_1,x_2,\dots,x_A,\alpha,\beta,\dots,\nu) = \psi_{\alpha}(x_1) \psi_{\beta}(x_2)\dots\psi_{\nu}(x_A). $$

Our functional is written as $$ E[\Phi] = \sum_{\mu=1}^A \int \psi_{\mu}^*(x_i)\hat{h}_0(x_i)\psi_{\mu}(x_i) dx_i + \frac{1}{2}\sum_{\mu=1}^A\sum_{\nu=1}^A \left[ \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\mu}(x_i)\psi_{\nu}(x_j)dx_idx_j- \int \psi_{\mu}^*(x_i)\psi_{\nu}^*(x_j) \hat{v}(r_{ij})\psi_{\nu}(x_i)\psi_{\mu}(x_j)dx_idx_j\right] $$ The more compact version reads $$ E[\Phi] = \sum_{\mu}^A \langle \mu | \hat{h}_0 | \mu\rangle+ \frac{1}{2}\sum_{\mu\nu}^A\left[\langle \mu\nu |\hat{v}|\mu\nu\rangle-\langle \nu\mu |\hat{v}|\mu\nu\rangle\right]. $$

Since the interaction is invariant under the interchange of two particles it means for example that we have $$ \langle \mu\nu|\hat{v}|\mu\nu\rangle = \langle \nu\mu|\hat{v}|\nu\mu\rangle, $$ or in the more general case $$ \langle \mu\nu|\hat{v}|\sigma\tau\rangle = \langle \nu\mu|\hat{v}|\tau\sigma\rangle. $$

The direct and exchange matrix elements can be brought together if we define the antisymmetrized matrix element $$ \langle \mu\nu|\hat{v}|\mu\nu\rangle_{AS}= \langle \mu\nu|\hat{v}|\mu\nu\rangle-\langle \mu\nu|\hat{v}|\nu\mu\rangle, $$ or for a general matrix element $$ \langle \mu\nu|\hat{v}|\sigma\tau\rangle_{AS}= \langle \mu\nu|\hat{v}|\sigma\tau\rangle-\langle \mu\nu|\hat{v}|\tau\sigma\rangle. $$ It has the symmetry property $$ \langle \mu\nu|\hat{v}|\sigma\tau\rangle_{AS}= -\langle \mu\nu|\hat{v}|\tau\sigma\rangle_{AS}=-\langle \nu\mu|\hat{v}|\sigma\tau\rangle_{AS}. $$ The antisymmetric matrix element is also hermitian, implying $$ \langle \mu\nu|\hat{v}|\sigma\tau\rangle_{AS}= \langle \sigma\tau|\hat{v}|\mu\nu\rangle_{AS}. $$

With these notations we rewrite the Hartree-Fock functional as $$ \begin{equation} \int \Phi^*\hat{H_I}\Phi d\mathbf{\tau} = \frac{1}{2}\sum_{\mu=1}^A\sum_{\nu=1}^A \langle \mu\nu|\hat{v}|\mu\nu\rangle_{AS}. \tag{10} \end{equation} $$

Adding the contribution from the one-body operator \( \hat{H}_0 \) to (10) we obtain the energy functional $$ \begin{equation} E[\Phi] = \sum_{\mu=1}^A \langle \mu | h | \mu \rangle + \frac{1}{2}\sum_{{\mu}=1}^A\sum_{{\nu}=1}^A \langle \mu\nu|\hat{v}|\mu\nu\rangle_{AS}. \tag{11} \end{equation} $$ In our coordinate space derivations below we will spell out the Hartree-Fock equations in terms of their integrals.

If we generalize the Euler-Lagrange equations to more variables and introduce \( N^2 \) Lagrange multipliers which we denote by \( \epsilon_{\mu\nu} \), we can write the variational equation for the functional of \( E \) $$ \delta E - \sum_{\mu\nu}^A \epsilon_{\mu\nu} \delta \int \psi_{\mu}^* \psi_{\nu} = 0. $$ For the orthogonal wave functions \( \psi_{i} \) this reduces to $$ \delta E - \sum_{\mu=1}^A \epsilon_{\mu} \delta \int \psi_{\mu}^* \psi_{\mu} = 0. $$

Variation with respect to the single-particle wave functions \( \psi_{\mu} \) yields then $$ \sum_{\mu=1}^A \int \delta\psi_{\mu}^*\hat{h_0}(x_i)\psi_{\mu} dx_i + \frac{1}{2}\sum_{{\mu}=1}^A\sum_{{\nu}=1}^A \left[ \int \delta\psi_{\mu}^*\psi_{\nu}^*\hat{v}(r_{ij})\psi_{\mu}\psi_{\nu} dx_idx_j- \int \delta\psi_{\mu}^*\psi_{\nu}^*\hat{v}(r_{ij})\psi_{\nu}\psi_{\mu} dx_idx_j \right]+ $$ $$ \sum_{\mu=1}^A \int \psi_{\mu}^*\hat{h_0}(x_i)\delta\psi_{\mu} dx_i + \frac{1}{2}\sum_{{\mu}=1}^A\sum_{{\nu}=1}^A \left[ \int \psi_{\mu}^*\psi_{\nu}^*\hat{v}(r_{ij})\delta\psi_{\mu}\psi_{\nu} dx_idx_j- \int \psi_{\mu}^*\psi_{\nu}^*\hat{v}(r_{ij})\psi_{\nu}\delta\psi_{\mu} dx_idx_j \right]- \sum_{{\mu}=1}^A E_{\mu} \int \delta\psi_{\mu}^* \psi_{\mu}dx_i - \sum_{{\mu}=1}^A E_{\mu} \int \psi_{\mu}^* \delta\psi_{\mu}dx_i = 0. $$

Although the variations \( \delta\psi \) and \( \delta\psi^* \) are not independent, they may in fact be treated as such, so that the terms dependent on either \( \delta\psi \) and \( \delta\psi^* \) individually may be set equal to zero. To see this, simply replace the arbitrary variation \( \delta\psi \) by \( i\delta\psi \), so that \( \delta\psi^* \) is replaced by \( -i\delta\psi^* \), and combine the two equations. We thus arrive at the Hartree-Fock equations $$ \begin{equation} \left[ -\frac{1}{2}\nabla_i^2+ \sum_{\nu=1}^A\int \psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\nu}(x_j)dx_j \right]\psi_{\mu}(x_i) - \left[ \sum_{{\nu}=1}^A \int\psi_{\nu}^*(x_j)\hat{v}(r_{ij})\psi_{\mu}(x_j) dx_j\right] \psi_{\nu}(x_i) = \epsilon_{\mu} \psi_{\mu}(x_i). \tag{12} \end{equation} $$ Notice that the integration \( \int dx_j \) implies an integration over the spatial coordinates \( \mathbf{r_j} \) and a summation over the spin-coordinate of fermion \( j \). We note that the factor of \( 1/2 \) in front of the sum involving the two-body interaction, has been removed. This is due to the fact that we need to vary both \( \delta\psi_{\mu}^* \) and \( \delta\psi_{\nu}^* \). Using the symmetry properties of the two-body interaction and interchanging \( \mu \) and \( \nu \) as summation indices, we obtain two identical terms.

The two first terms in the last equation are the one-body kinetic energy and the electron-nucleus potential. The third or direct term is the averaged electronic repulsion of the other electrons. As written, the term includes the self-interaction of electrons when \( \mu=\nu \). The self-interaction is cancelled in the fourth term, or the exchange term. The exchange term results from our inclusion of the Pauli principle and the assumed determinantal form of the wave-function. Equation (12), in addition to the kinetic energy and the attraction from the atomic nucleus that confines the motion of a single electron, represents now the motion of a single-particle modified by the two-body interaction. The additional contribution to the Schroedinger equation due to the two-body interaction, represents a mean field set up by all the other bystanding electrons, the latter given by the sum over all single-particle states occupied by \( N \) electrons.

The Hartree-Fock equation is an example of an integro-differential equation. These equations involve repeated calculations of integrals, in addition to the solution of a set of coupled differential equations. The Hartree-Fock equations can also be rewritten in terms of an eigenvalue problem. The solution of an eigenvalue problem represents often a more practical algorithm and the solution of coupled integro-differential equations. This alternative derivation of the Hartree-Fock equations is given below.

Analysis of Hartree-Fock equations in coordinate space

A theoretically convenient form of the Hartree-Fock equation is to regard the direct and exchange operator defined through $$ \begin{equation*} V_{\mu}^{d}(x_i) = \int \psi_{\mu}^*(x_j) \hat{v}(r_{ij})\psi_{\mu}(x_j) dx_j \end{equation*} $$ and $$ \begin{equation*} V_{\mu}^{ex}(x_i) g(x_i) = \left(\int \psi_{\mu}^*(x_j) \hat{v}(r_{ij})g(x_j) dx_j \right)\psi_{\mu}(x_i), \end{equation*} $$ respectively.

The function \( g(x_i) \) is an arbitrary function, and by the substitution \( g(x_i) = \psi_{\nu}(x_i) \) we get $$ \begin{equation*} V_{\mu}^{ex}(x_i) \psi_{\nu}(x_i) = \left(\int \psi_{\mu}^*(x_j) \hat{v}(r_{ij})\psi_{\nu}(x_j) dx_j\right)\psi_{\mu}(x_i). \end{equation*} $$ We may then rewrite the Hartree-Fock equations as $$ \hat{h}^{HF}(x_i) \psi_{\nu}(x_i) = \epsilon_{\nu}\psi_{\nu}(x_i), $$ with $$ \hat{h}^{HF}(x_i)= \hat{h}_0(x_i) + \sum_{\mu=1}^AV_{\mu}^{d}(x_i) - \sum_{\mu=1}^AV_{\mu}^{ex}(x_i), $$ and where \( \hat{h}_0(i) \) is the one-body part. The latter is normally chosen as a part which yields solutions in closed form. The harmonic oscilltor is a classical problem thereof. We normally rewrite the last equation as $$ \hat{h}^{HF}(x_i)= \hat{h}_0(x_i) + \hat{u}^{HF}(x_i). $$

Hartree-Fock by varying the coefficients of a wave function expansion

Another possibility is to expand the single-particle functions in a known basis and vary the coefficients, that is, the new single-particle wave function is written as a linear expansion in terms of a fixed chosen orthogonal basis (for example the well-known harmonic oscillator functions or the hydrogen-like functions etc). We define our new Hartree-Fock single-particle basis by performing a unitary transformation on our previous basis (labelled with greek indices) as $$ \begin{equation} \psi_p^{HF} = \sum_{\lambda} C_{p\lambda}\phi_{\lambda}. \tag{13} \end{equation} $$ In this case we vary the coefficients \( C_{p\lambda} \). If the basis has infinitely many solutions, we need to truncate the above sum. We assume that the basis \( \phi_{\lambda} \) is orthogonal. A unitary transformation keeps the orthogonality, as discussed in exercise 1 below.

It is normal to choose a single-particle basis defined as the eigenfunctions of parts of the full Hamiltonian. The typical situation consists of the solutions of the one-body part of the Hamiltonian, that is we have $$ \hat{h}_0\phi_{\lambda}=\epsilon_{\lambda}\phi_{\lambda}. $$ The single-particle wave functions \( \phi_{\lambda}({\bf r}) \), defined by the quantum numbers \( \lambda \) and \( {\bf r} \) are defined as the overlap $$ \phi_{\lambda}({\bf r}) = \langle {\bf r} | \lambda \rangle . $$

In our discussions hereafter we will use our definitions of single-particle states above and below the Fermi (\( F \)) level given by the labels \( ijkl\dots \le F \) for so-called single-hole states and \( abcd\dots > F \) for so-called particle states. For general single-particle states we employ the labels \( pqrs\dots \).

In Eq. (11), restated here $$ E[\Phi] = \sum_{\mu=1}^A \langle \mu | h | \mu \rangle + \frac{1}{2}\sum_{{\mu}=1}^A\sum_{{\nu}=1}^A \langle \mu\nu|\hat{v}|\mu\nu\rangle_{AS}, $$ we found the expression for the energy functional in terms of the basis function \( \phi_{\lambda}({\bf r}) \). We then varied the above energy functional with respect to the basis functions \( |\mu \rangle \). Now we are interested in defining a new basis defined in terms of a chosen basis as defined in Eq. (13). We can then rewrite the energy functional as $$ \begin{equation} E[\Phi^{HF}] = \sum_{i=1}^A \langle i | h | i \rangle + \frac{1}{2}\sum_{ij=1}^A\langle ij|\hat{v}|ij\rangle_{AS}, \tag{14} \end{equation} $$ where \( \Phi^{HF} \) is the new Slater determinant defined by the new basis of Eq. (13).

Using Eq. (13) we can rewrite Eq. (14) as $$ \begin{equation} E[\Psi] = \sum_{i=1}^A \sum_{\alpha\beta} C^*_{i\alpha}C_{i\beta}\langle \alpha | h | \beta \rangle + \frac{1}{2}\sum_{ij=1}^A\sum_{{\alpha\beta\gamma\delta}} C^*_{i\alpha}C^*_{j\beta}C_{i\gamma}C_{j\delta}\langle \alpha\beta|\hat{v}|\gamma\delta\rangle_{AS}. \tag{15} \end{equation} $$

We wish now to minimize the above functional. We introduce again a set of Lagrange multipliers, noting that since \( \langle i | j \rangle = \delta_{i,j} \) and \( \langle \alpha | \beta \rangle = \delta_{\alpha,\beta} \), the coefficients \( C_{i\gamma} \) obey the relation $$ \langle i | j \rangle=\delta_{i,j}=\sum_{\alpha\beta} C^*_{i\alpha}C_{i\beta}\langle \alpha | \beta \rangle= \sum_{\alpha} C^*_{i\alpha}C_{i\alpha}, $$ which allows us to define a functional to be minimized that reads $$ \begin{equation} F[\Phi^{HF}]=E[\Phi^{HF}] - \sum_{i=1}^A\epsilon_i\sum_{\alpha} C^*_{i\alpha}C_{i\alpha}. \end{equation} $$

Minimizing with respect to \( C^*_{i\alpha} \), remembering that the equations for \( C^*_{i\alpha} \) and \( C_{i\alpha} \) can be written as two independent equations, we obtain $$ \frac{d}{dC^*_{i\alpha}}\left[ E[\Phi^{HF}] - \sum_{j}\epsilon_j\sum_{\alpha} C^*_{j\alpha}C_{j\alpha}\right]=0, $$ which yields for every single-particle state \( i \) and index \( \alpha \) (recalling that the coefficients \( C_{i\alpha} \) are matrix elements of a unitary (or orthogonal for a real symmetric matrix) matrix) the following Hartree-Fock equations $$ \sum_{\beta} C_{i\beta}\langle \alpha | h | \beta \rangle+ \sum_{j=1}^A\sum_{\beta\gamma\delta} C^*_{j\beta}C_{j\delta}C_{i\gamma}\langle \alpha\beta|\hat{v}|\gamma\delta\rangle_{AS}=\epsilon_i^{HF}C_{i\alpha}. $$

We can rewrite this equation as (changing dummy variables) $$ \sum_{\beta} \left\{\langle \alpha | h | \beta \rangle+ \sum_{j}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}\right\}C_{i\beta}=\epsilon_i^{HF}C_{i\alpha}. $$ Note that the sums over greek indices run over the number of basis set functions (in principle an infinite number).

Defining $$ h_{\alpha\beta}^{HF}=\langle \alpha | h | \beta \rangle+ \sum_{j=1}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}, $$ we can rewrite the new equations as $$ \begin{equation} \sum_{\gamma}h_{\alpha\beta}^{HF}C_{i\beta}=\epsilon_i^{HF}C_{i\alpha}. \tag{16} \end{equation} $$ The latter is nothing but a standard eigenvalue problem. Compared with Eq. (12), we see that we do not need to compute any integrals in an iterative procedure for solving the equations. It suffices to tabulate the matrix elements \( \langle \alpha | h | \beta \rangle \) and \( \langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS} \) once and for all. Successive iterations require thus only a look-up in tables over one-body and two-body matrix elements. These details will be discussed below when we solve the Hartree-Fock equations numerical.

Hartree-Fock algorithm

Our Hartree-Fock matrix is thus $$ \hat{h}_{\alpha\beta}^{HF}=\langle \alpha | \hat{h}_0 | \beta \rangle+ \sum_{j=1}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}. $$ The Hartree-Fock equations are solved in an iterative waym starting with a guess for the coefficients \( C_{j\gamma}=\delta_{j,\gamma} \) and solving the equations by diagonalization till the new single-particle energies \( \epsilon_i^{\mathrm{HF}} \) do not change anymore by a prefixed quantity.

Normally we assume that the single-particle basis \( |\beta\rangle \) forms an eigenbasis for the operator \( \hat{h}_0 \), meaning that the Hartree-Fock matrix becomes $$ \hat{h}_{\alpha\beta}^{HF}=\epsilon_{\alpha}\delta_{\alpha,\beta}+ \sum_{j=1}^A\sum_{\gamma\delta} C^*_{j\gamma}C_{j\delta}\langle \alpha\gamma|\hat{v}|\beta\delta\rangle_{AS}. $$ The Hartree-Fock eigenvalue problem $$ \sum_{\beta}\hat{h}_{\alpha\beta}^{HF}C_{i\beta}=\epsilon_i^{\mathrm{HF}}C_{i\alpha}, $$ can be written out in a more compact form as $$ \hat{h}^{HF}\hat{C}=\epsilon^{\mathrm{HF}}\hat{C}. $$

The Hartree-Fock equations are, in their simplest form, solved in an iterative way, starting with a guess for the coefficients \( C_{i\alpha} \). We label the coefficients as \( C_{i\alpha}^{(n)} \), where the subscript \( n \) stands for iteration \( n \). To set up the algorithm we can proceed as follows:

• We start with a guess \( C_{i\alpha}^{(0)}=\delta_{i,\alpha} \). Alternatively, we could have used random starting values as long as the vectors are normalized. Another possibility is to give states below the Fermi level a larger weight.
• The Hartree-Fock matrix simplifies then to (assuming that the coefficients \( C_{i\alpha} \) are real)

Solving the Hartree-Fock eigenvalue problem yields then new eigenvectors \( C_{i\alpha}^{(1)} \) and eigenvalues \( \epsilon_i^{HF(1)} \).

• With the new eigenvalues we can set up a new Hartree-Fock potential

• TODO: add code with Hartree-Fock for nuclear system

Analysis of Hartree-Fock equations and Koopman's theorem

We can rewrite the ground state energy by adding and subtracting \( \hat{u}^{HF}(x_i) \) $$ E_0^{HF} =\langle \Phi_0 | \hat{H} | \Phi_0\rangle = \sum_{i\le F}^A \langle i | \hat{h}_0 +\hat{u}^{HF}| j\rangle+ \frac{1}{2}\sum_{i\le F}^A\sum_{j \le F}^A\left[\langle ij |\hat{v}|ij \rangle-\langle ij|\hat{v}|ji\rangle\right]-\sum_{i\le F}^A \langle i |\hat{u}^{HF}| i\rangle, $$ which results in $$ E_0^{HF} = \sum_{i\le F}^A \varepsilon_i^{HF} + \frac{1}{2}\sum_{i\le F}^A\sum_{j \le F}^A\left[\langle ij |\hat{v}|ij \rangle-\langle ij|\hat{v}|ji\rangle\right]-\sum_{i\le F}^A \langle i |\hat{u}^{HF}| i\rangle. $$ Our single-particle states \( ijk\dots \) are now single-particle states obtained from the solution of the Hartree-Fock equations.

Using our definition of the Hartree-Fock single-particle energies we obtain then the following expression for the total ground-state energy $$ E_0^{HF} = \sum_{i\le F}^A \varepsilon_i - \frac{1}{2}\sum_{i\le F}^A\sum_{j \le F}^A\left[\langle ij |\hat{v}|ij \rangle-\langle ij|\hat{v}|ji\rangle\right]. $$ This form will be used in our discussion of Koopman's theorem.

In the atomic physics case we have $$ E[\Phi^{\mathrm{HF}}(N)] = \sum_{i=1}^H \langle i | \hat{h}_0 | i \rangle + \frac{1}{2}\sum_{ij=1}^N\langle ij|\hat{v}|ij\rangle_{AS}, $$ where \( \Phi^{\mathrm{HF}}(N) \) is the new Slater determinant defined by the new basis of Eq. (13) for \( N \) electrons (same \( Z \)). If we assume that the single-particle wave functions in the new basis do not change when we remove one electron or add one electron, we can then define the corresponding energy for the \( N-1 \) systems as $$ E[\Phi^{\mathrm{HF}}(N-1)] = \sum_{i=1; i\ne k}^N \langle i | \hat{h}_0 | i \rangle + \frac{1}{2}\sum_{ij=1;i,j\ne k}^N\langle ij|\hat{v}|ij\rangle_{AS}, $$ where we have removed a single-particle state \( k\le F \), that is a state below the Fermi level.

Calculating the difference $$ E[\Phi^{\mathrm{HF}}(N)]- E[\Phi^{\mathrm{HF}}(N-1)] = \langle k | \hat{h}_0 | k \rangle + \frac{1}{2}\sum_{i=1;i\ne k}^N\langle ik|\hat{v}|ik\rangle_{AS} + \frac{1}{2}\sum_{j=1;j\ne k}^N\langle kj|\hat{v}|kj\rangle_{AS}, $$ we obtain $$ E[\Phi^{\mathrm{HF}}(N)]- E[\Phi^{\mathrm{HF}}(N-1)] = \langle k | \hat{h}_0 | k \rangle +\sum_{j=1}^N\langle kj|\hat{v}|kj\rangle_{AS} $$ which is just our definition of the Hartree-Fock single-particle energy $$ E[\Phi^{\mathrm{HF}}(N)]- E[\Phi^{\mathrm{HF}}(N-1)] = \epsilon_k^{\mathrm{HF}} $$

Similarly, we can now compute the difference (we label the single-particle states above the Fermi level as \( abcd > F \)) $$ E[\Phi^{\mathrm{HF}}(N+1)]- E[\Phi^{\mathrm{HF}}(N)]= \epsilon_a^{\mathrm{HF}}. $$ These two equations can thus be used to the electron affinity or ionization energies, respectively. Koopman's theorem states that for example the ionization energy of a closed-shell system is given by the energy of the highest occupied single-particle state. If we assume that changing the number of electrons from \( N \) to \( N+1 \) does not change the Hartree-Fock single-particle energies and eigenfunctions, then Koopman's theorem simply states that the ionization energy of an atom is given by the single-particle energy of the last bound state. In a similar way, we can also define the electron affinities.

As an example, consider a simple model for atomic sodium, Na. Neutral sodium has eleven electrons, with the weakest bound one being confined the \( 3s \) single-particle quantum numbers. The energy needed to remove an electron from neutral sodium is rather small, 5.1391 eV, a feature which pertains to all alkali metals. Having performed a Hartree-Fock calculation for neutral sodium would then allows us to compute the ionization energy by using the single-particle energy for the \( 3s \) states, namely \( \epsilon_{3s}^{\mathrm{HF}} \).

From these considerations, we see that Hartree-Fock theory allows us to make a connection between experimental observables (here ionization and affinity energies) and the underlying interactions between particles. In this sense, we are now linking the dynamics and structure of a many-body system with the laws of motion which govern the system. Our approach is a reductionistic one, meaning that we want to understand the laws of motion in terms of the particles or degrees of freedom which we believe are the fundamental ones. Our Slater determinant, being constructed as the product of various single-particle functions, follows this philosophy.

With similar arguments as in atomic physics, we can now use Hartree-Fock theory to make a link between nuclear forces and separation energies. Changing to nuclear system, we define $$ E[\Phi^{\mathrm{HF}}(A)] = \sum_{i=1}^A \langle i | \hat{h}_0 | i \rangle + \frac{1}{2}\sum_{ij=1}^A\langle ij|\hat{v}|ij\rangle_{AS}, $$ where \( \Phi^{\mathrm{HF}}(A) \) is the new Slater determinant defined by the new basis of Eq. (13) for \( A \) nucleons, where \( A=N+Z \), with \( N \) now being the number of neutrons and \( Z \) th enumber of protons. If we assume again that the single-particle wave functions in the new basis do not change from a nucleus with \( A \) nucleons to a nucleus with \( A-1 \) nucleons, we can then define the corresponding energy for the \( A-1 \) systems as $$ E[\Phi^{\mathrm{HF}}(A-1)] = \sum_{i=1; i\ne k}^A \langle i | \hat{h}_0 | i \rangle + \frac{1}{2}\sum_{ij=1;i,j\ne k}^A\langle ij|\hat{v}|ij\rangle_{AS}, $$ where we have removed a single-particle state \( k\le F \), that is a state below the Fermi level.

Calculating the difference $$ E[\Phi^{\mathrm{HF}}(A)]- E[\Phi^{\mathrm{HF}}(A-1)] = \langle k | \hat{h}_0 | k \rangle + \frac{1}{2}\sum_{i=1;i\ne k}^A\langle ik|\hat{v}|ik\rangle_{AS} + \frac{1}{2}\sum_{j=1;j\ne k}^A\langle kj|\hat{v}|kj\rangle_{AS}, $$ which becomes $$ E[\Phi^{\mathrm{HF}}(A)]- E[\Phi^{\mathrm{HF}}(A-1)] = \langle k | \hat{h}_0 | k \rangle +\sum_{j=1}^A\langle kj|\hat{v}|kj\rangle_{AS} $$ which is just our definition of the Hartree-Fock single-particle energy $$ E[\Phi^{\mathrm{HF}}(A)]- E[\Phi^{\mathrm{HF}}(A-1)] = \epsilon_k^{\mathrm{HF}} $$

Similarly, we can now compute the difference (recall that the single-particle states \( abcd > F \)) $$ E[\Phi^{\mathrm{HF}}(A+1)]- E[\Phi^{\mathrm{HF}}(A)]= \epsilon_a^{\mathrm{HF}}. $$ If we then recall that the binding energy differences $$ BE(A)-BE(A-1) \hspace{0.5cm} \mathrm{and} \hspace{0.5cm} BE(A+1)-BE(A), $$ define the separation energies, we see that the Hartree-Fock single-particle energies can be used to define separation energies. We have thus our first link between nuclear forces (included in the potential energy term) and an observable quantity defined by differences in binding energies.

We have thus the following interpretations (if the single-particle fields do not change) $$ BE(A)-BE(A-1)\approx E[\Phi^{\mathrm{HF}}(A)]- E[\Phi^{\mathrm{HF}}(A-1)] = \epsilon_k^{\mathrm{HF}}, $$ and $$ BE(A+1)-BE(A)\approx E[\Phi^{\mathrm{HF}}(A+1)]- E[\Phi^{\mathrm{HF}}(A)] = \epsilon_a^{\mathrm{HF}}. $$ If we use \( {}^{16}\mbox{O} \) as our closed-shell nucleus, we could then interpret the separation energy $$ BE(^{16}\mathrm{O})-BE(^{15}\mathrm{O})\approx \epsilon_{0p^{\nu}_{1/2}}^{\mathrm{HF}}, $$ and $$ BE(^{16}\mathrm{O})-BE(^{15}\mathrm{N})\approx \epsilon_{0p^{\pi}_{1/2}}^{\mathrm{HF}}. $$

Similalry, we could interpret $$ BE(^{17}\mathrm{O})-BE(^{16}\mathrm{O})\approx \epsilon_{0d^{\nu}_{5/2}}^{\mathrm{HF}}, $$ and $$ BE(^{17}\mathrm{F})-BE(^{16}\mathrm{O})\approx\epsilon_{0d^{\pi}_{5/2}}^{\mathrm{HF}}. $$ We can continue like this for all \( A\pm 1 \) nuclei where \( A \) is a good closed-shell (or subshell closure) nucleus. Examples are \( {}^{22}\mbox{O} \), \( {}^{24}\mbox{O} \), \( {}^{40}\mbox{Ca} \), \( {}^{48}\mbox{Ca} \), \( {}^{52}\mbox{Ca} \), \( {}^{54}\mbox{Ca} \), \( {}^{56}\mbox{Ni} \), \( {}^{68}\mbox{Ni} \), \( {}^{78}\mbox{Ni} \), \( {}^{90}\mbox{Zr} \), \( {}^{88}\mbox{Sr} \), \( {}^{100}\mbox{Sn} \), \( {}^{132}\mbox{Sn} \) and \( {}^{208}\mbox{Pb} \), to mention some possile cases.

We can thus make our first interpretation of the separation energies in terms of the simplest possible many-body theory. If we also recall that the so-called energy gap for neutrons (or protons) is defined as $$ \Delta S_n= 2BE(N,Z)-BE(N-1,Z)-BE(N+1,Z), $$ for neutrons and the corresponding gap for protons $$ \Delta S_p= 2BE(N,Z)-BE(N,Z-1)-BE(N,Z+1), $$ we can define the neutron and proton energy gaps for \( {}^{16}\mbox{O} \) as $$ \Delta S_{\nu}=\epsilon_{0d^{\nu}_{5/2}}^{\mathrm{HF}}-\epsilon_{0p^{\nu}_{1/2}}^{\mathrm{HF}}, $$ and $$ \Delta S_{\pi}=\epsilon_{0d^{\pi}_{5/2}}^{\mathrm{HF}}-\epsilon_{0p^{\pi}_{1/2}}^{\mathrm{HF}}. $$

Exercise 1: Derivation of Hartree-Fock equations

Consider a Slater determinant built up of single-particle orbitals \( \psi_{\lambda} \), with \( \lambda = 1,2,\dots,N \).

The unitary transformation $$ \psi_a = \sum_{\lambda} C_{a\lambda}\phi_{\lambda}, $$ brings us into the new basis. The new basis has quantum numbers \( a=1,2,\dots,N \).

a) Show that the new basis is orthonormal.

b) Show that the new Slater determinant constructed from the new single-particle wave functions can be written as the determinant based on the previous basis and the determinant of the matrix \( C \).

c) Show that the old and the new Slater determinants are equal up to a complex constant with absolute value unity.

Hint.

Exercise 2: Derivation of Hartree-Fock equations

Consider the Slater determinant $$ \Phi_{0}=\frac{1}{\sqrt{n!}}\sum_{p}(-)^{p}P \prod_{i=1}^{n}\psi_{\alpha_{i}}(x_{i}). $$ A small variation in this function is given by $$ \delta\Phi_{0}=\frac{1}{\sqrt{n!}}\sum_{p}(-)^{p}P \psi_{\alpha_{1}}(x_{1})\psi_{\alpha_{2}}(x_{2})\dots \psi_{\alpha_{i-1}}(x_{i-1})(\delta\psi_{\alpha_{i}}(x_{i})) \psi_{\alpha_{i+1}}(x_{i+1})\dots\psi_{\alpha_{n}}(x_{n}). $$

a) Show that $$ \langle \delta\Phi_{0}|\sum_{i=1}^{n}\left\{t(x_{i})+u(x_{i}) \right\}+\frac{1}{2} \sum_{i\neq j=1}^{n}v(x_{i},x_{j})|\Phi_{0}\rangle=\sum_{i=1}^{n}\langle \delta\psi_{\alpha_{i}}|\hat{t}+\hat{u} |\phi_{\alpha_{i}}\rangle +\sum_{i\neq j=1}^{n}\left\{\langle\delta\psi_{\alpha_{i}} \psi_{\alpha_{j}}|\hat{v}|\psi_{\alpha_{i}}\psi_{\alpha_{j}}\rangle- \langle\delta\psi_{\alpha_{i}}\psi_{\alpha_{j}}|\hat{v} |\psi_{\alpha_{j}}\psi_{\alpha_{i}}\rangle\right\} $$

Exercise 3: Hartree-Fock calculations in atomic physics

In this exercise we will develop two simple models for studying the helium atom (with two electrons) and the beryllium atom with four electrons. The reason why we chose these systems is simply due to the fact that we can use hydrogen-like single-particle state functions. These functions allow for simple evaluations of the two-body integrals needed in a Hartree-Fock calculation. In most cases, these integrals have a simple closed form solution, as indicated in the table below.

After having introduced the Born-Oppenheimer approximation which effectively freezes out the nucleonic degrees of freedom, the Hamiltonian for \( N \) electrons takes the following form $$ \hat{H} = \sum_{i=1}^{N} t(x_i) - \sum_{i=1}^{N} k\frac{Ze^2}{r_i} + \sum_{i < j}^{N} \frac{ke^2}{r_{ij}}, $$ with \( k=1.44 \) eVnm. Througout this work we will use atomic units, this means that \( \hbar=c=e=m_e=1 \). The constant \( k \) becomes also equal 1. The resulting energies have to be multiplied by \( 2\times 13.6 \) eV in order to obtain energies in eletronvolts.

We can rewrite our Hamiltonians as $$ \hat{H} = \hat{H_0} + \hat{H_I} = \sum_{i=1}^{N}\hat{h}_0(x_i) + \sum_{i < j}^{N}\frac{1}{r_{ij}}, $$ where we have defined \( r_{ij}=| {\bf r}_i-{\bf r}_j| \) and \( \hat{h}_0(x_i) = \hat{t}(x_i) - \frac{Z}{r_i} \) The variable \( x \) contains both the spatial coordinates and the spin values. The first term of the previous equation, \( H_0 \), is the sum of the \( N \) one-body Hamiltonians \( \hat{h}_0 \). Each individual Hamiltonian \( \hat{h}_0 \) contains the kinetic energy operator of an electron and its potential energy due to the attraction of the nucleus. The second term, \( H_I \), is the sum of the \( N(N-1)/2 \) two-body interactions between each pair of electrons. Note that the double sum carries a restriction \( i < j \).

As basis functions for our calculations we will use hydrogen-like single-particle functions. This means the onebody operator is diagonal in this basis for states \( i,j \) with quantum numbers \( nlm_lsm_s \) with energies $$\langle i|\hat{h}_0| j\rangle = -Z^2/2n^2\delta_{i,j}. $$ The quantum number \( n \) refers to the number of nodes of the wave function. Observe that this expectation value is independent of spin.

We will in all calculations here restrict ourselves to only so-called \( s \) -waves, that is the orbital momentum \( l \) is zero. We will also limit the quantum number \( n \) to \( n\le 3 \). It means that every \( ns \) state can accomodate two electrons due to the spin degeneracy. In the calculations you will need the Coulomb interaction with matrix elements involving single-particle wave functions with \( l=0 \) only, the so-called \( s \)-waves. We need only the radial part since the spherical harmonics for the \( s \)-waves are rather simple. We omit single-particle states with \( l> 0 \). Our radial part of the wave functions are $$ R_{n0}(r)=\left(\frac{2Z}{n}\right)^{3/2}\sqrt{\frac{(n-1)!}{2n\times n!}}L_{n-1}^1(\frac{2Zr}{n})\exp{(-\frac{Zr}{n})}, $$ where \( L_{n-1}^1(r) \) are the so-called Laguerre polynomials. These wave functions can then be used to compute the direct part of the Coulomb interaction $$ \langle \alpha\beta| V| \gamma\delta\rangle = \int r_1^2dr_1 \int r_2^2dr_2R_{n_{\alpha}0}^*(r_1) R_{n_{\beta}0}^*(r_2) \frac{1}{| {\bf r}_1-{\bf r}_2|}R_{n_{\gamma}0}(r_1)R_{n_{\delta}0}(r_2) $$ Observe that this is only the radial integral and that the labels \( \alpha\beta\gamma\delta \) refer only to the quantum numbers \( nlm_l \), with \( m_l \) the projection of the orbital momentum \( l \). A similar expression can be found for the exchange part. Since we have restricted ourselves to only \( s \)-waves, these integrals are straightforward but tedious to calculate. As an addendum to this project we list all closed-form expressions for the relevant matrix elements. Note well that these matrix elements do not include spin. When setting up the final antisymmetrized matrix elements you need to consider the spin degrees of freedom as well. Please pay in particular special attention to the exchange part and the pertinent spin values of the single-particle states.

We will also, for both helium and beryllium assume that the many-particle states we construct have always the same total spin projection \( M_S=0 \). This means that if we excite one or two particles from the ground state, the spins of the various single-particle states should always sum up to zero.

a) We start with the helium atom and define our single-particle Hilbert space to consist of the single-particle orbits \( 1s \), \( 2s \) and \( 3s \), with their corresponding spin degeneracies.

Set up the ansatz for the ground state \( |c\rangle = |\Phi_0\rangle \) in second quantization and define a table of single-particle states. Construct thereafter all possible one-particle-one-hole excitations \( |\Phi_i^a\rangle \) where \( i \) refer to levels below the Fermi level (define this level) and \( a \) refers to particle states. Define particles and holes. The Slater determinants have to be written in terms of the respective creation and annihilation operators. The states you construct should all have total spin projection \( M_S=0 \). Construct also all possible two-particle-two-hole states \( |\Phi_{ij}^{ab}\rangle \) in a second quantization representation.

b) Define the Hamiltonian in a second-quantized form and use this to compute the expectation value of the ground state (defining the so-called reference energy and later our Hartree-Fock functional) of the helium atom. Show that it is given by $$ E[\Phi_0] = \langle c | \hat{H}| c \rangle = \sum_{i} \langle i | \hat{h}_0 | i\rangle+ \frac{1}{2}\sum_{ij}\left[\langle ij |\frac{1}{r}|ij\rangle-\langle ij |\frac{1}{r}|ji\rangle\right]. $$ Define properly the sums keeping in mind that the states \( ij \) refer to all quantum numbers \( nlm_lsm_s \). Use the values for the various matrix elements listed at the end of the project to find the value of \( E \) as function of \( Z \) and compute \( E \) as function of \( Z \).

c) Hereafter we will limit ourselves to a system which now contains only one-particle-one-hole excitations beyond the chosen state \( |c\rangle \). Using the possible Slater determinants from exercise a) for the helium atom, find the expressions (without inserting the explicit values for the matrix elements first) for $$ \langle c | \hat{H}| \Phi_i^a \rangle, $$ and $$ \langle \Phi_i^a | \hat{H}| \Phi_j^b \rangle. $$ Represent these expressions in a diagrammatic form, both for the onebody part and the two-body part of the Hamiltonian.

Insert then the explicit values for the various matrix elements and set up the final Hamiltonian matrix and diagonalize it using for example Octave, Matlab, Python, C++ or Fortran as programming tools.

Compare your results from those of exercise b) and comment your results. The exact energy with our Hamiltonian is \( -2.9037 \) atomic units for helium. This value is also close to the experimental energy.

d) We repeat exercises b) and c) but now for the beryllium atom. Define the ansatz for \( |c\rangle \) and limit yourself again to one-particle-one-hole excitations. Compute the reference energy \( \langle c | \hat{H}| c \rangle \) as function of \( Z \). Thereafter you will need to set up the appropriate Hamiltonian matrix which involves also one-particle-one-hole excitations. Diagonalize this matrix and compare your eigenvalues with \( \langle c | \hat{H}| c \rangle \) as function of \( Z \) and comment your results. The exact energy with our Hamiltonian is \( -14.6674 \) atomic units for beryllium. This value is again close to the experimental energy.\newline\newline\newline With a given energy functional, we can perform at least two types of variational strategies

1. Vary the Slater determinant by changing the spatial part of the single-particle

1. Expand the single-particle functions in a known basis and vary the coefficients,

The unitary transformation \[ \psi_p = \sum_{\lambda} C_{p\lambda}\phi_{\lambda}, \] brings us into the new basis \( \psi \). The new basis is orthonormal and \( C \) is a unitary matrix.

e) Minimizing with respect to \( C^*_{p\alpha} \), remembering that \( C^*_{p\alpha} \) and \( C_{p\alpha} \) (and that the indices contain all single-particle quantum numbers including spin) are independent and defining $$ h_{\alpha\gamma}^{HF}=\langle \alpha | h | \gamma \rangle+ \sum_{p}\sum_{\beta\delta} C^*_{p\beta}C_{p\delta}\langle \alpha\beta|V|\gamma\delta\rangle_{AS}, $$ show that you can write the Hartree-Fock equations as $$ \sum_{\gamma}h_{\alpha\gamma}^{HF}C_{p\gamma}=\epsilon_p^{\mathrm{HF}}C_{p\alpha}. $$ Explain the meaning of the different terms and define the Hartree-Fock operator in second quantization. Write down its diagrammatic representation as well. The greek letters refer to the wave functions in the original basis (in our case the hydrogen-like wave functions) while roman letters refer to the new basis.

f) The Hartree-Fock equations with a variation of the coefficients \( C_{p\alpha} \) lead to an eigenvalue problem whose eigenvectors are the coefficients \( C_{p\alpha} \) and eigenvalues are the new single-particle energies. Use the single-particle orbits \( 1s-3s \) and set up the Hartree-Fock matrix for both the helium atom and the beryllium atom. Find after the first diagonalization the new single-particle energies and the new ground state energy. Compare these results with those you obtained under the minimization of the ground states as functions of \( Z \) and the full diagonalization. When setting up the Hartree-Fock matrix in the first iteration, our guess for the coefficients \( C_{p\beta} \) etc. is \( C_{p\beta}=1 \) for \( p=\beta \) and zero else. The final stage is to set up an iterative scheme where you use the new wave functions determined via the coefficients \( C_{p\alpha} \) to solve iteratively the Hartree-Fock equations till a given self-consistency is reached. A typical way of doing this is to compare the single-particle energies from the previous iteration with those obtained from the new diagonalization. If the total difference is smaller than a prefixed value, the iterative process is stopped.

Compare these results with the those you obtained under the minimization of the ground states as functions of \( Z \) and the full diagonalization. Discuss your results.

Solution.

"""
Program for solving Hartree Fock iteratively.
1. Import radial integral results
2. Set up HF matrix using C = I as initial guess
3. Solve eigenvalue problem using numpy.linalg
4. Use resulting eigenvectors to assemble new HF matrix
5. Repeat steps 3 and 4 until convergence is met.
6. Calculate ground state energy from resulting C.
"""

from __future__ import division
from sympy import *
from numpy import *
import pickle, sys

Z = Symbol("Z")

class HF:
    """
    Class for solving the Hartree Fock equations iteratively.
    """

    def __init__(self, N, basis, Z_value, first_C='identity'):
        """
        N is the number of particles in the system and basis
        is the single-particle basis for the system.
        """
        # Read in the radial_integrals from pickled object
        with open("radial_integrals.p", "rb") as infile:
            self.radial_integrals = pickle.load(infile)
        
        self.N = N # number of particles in the system
        self.basis = basis # single-particle basis for the system
        self.Z_value = Z_value # atomic number of the atom
        self.n = len(basis) # number of single-particle basis states
        self.ek = array((0,0,0,0,0,0)) # new single-particle energies
        self.E = 0 # energy
        # Set up the first coefficient matrix to be used
        if first_C == 'identity':
            self.C = identity(self.n)
        elif first_C == 'zero':
            self.C = zeros((self.n, self.n))
        elif first_C == 'rand':
            self.C = random.rand(self.n, self.n)
        else:
            print "first_C argument not understood"
            print "Legal values are: 'identity', 'zero', 'rand'"
            sys.exit(1)

        self.h_HF = zeros((self.n, self.n))
        self.assemble_HF_matrix() # set up HF matrix for C = I

    def h_0(self, p, q):
        """
        Takes the integer values of states and returns the 
        asymmetrized twobody matrix element <pq||rs>.
        """
        n1, s1 = self.basis[p]
        n2, s2 = self.basis[q]

        if n1 != n2 or s1 != s2:
            return 0
        else:
            return -Z**2/(2*n1**2)

    def rad(self, n1, n2, n3, n4):
        """
        Returns the radial integral <n1, n2|v|n3, n4>.
        """
        return self.radial_integrals[n1-1, n2-1, n3-1, n4-1]

    def h_1(self, p, q, r, s):
        """
        Takes the integer values of four basis-states and returns
        the asymmetrized twobody matrix element <pq||rs>.
        """
        n1, s1 = self.basis[p]
        n2, s2 = self.basis[q]
        n3, s3 = self.basis[r]
        n4, s4 = self.basis[s]

        if s1 == s2 == s3 == s4:
            return self.rad(n1, n2, n3, n4) - self.rad(n1, n2, n4, n3)
        if s1 == s3 and s2 == s4:
            return self.rad(n1, n2, n3, n4)
        if s1 == s4 and s2 == s3:
            return -self.rad(n1, n2, n4, n3)
        else:
            return 0

    def assemble_HF_matrix(self):
        """
        Assemble the HF matric from the coefficient matrix.
        """
        n, N = self.n, self.N
        C = self.C

        for a in range(n):
            for g in range(n):
                s = self.h_0(a,g)
                for p in range(N):
                    for b in range(n):
                        for d in range(n):
                            s += C[p,b]*C[p,d]*self.h_1(a,b,g,d)

                self.h_HF[a,g] = s.subs(Z, self.Z_value)

    def reorder_coefficients(self):
        ek, C = self.ek, self.C

        # Sort eigenvalues and coefficient matrix using numpy.argsort
        indices = argsort(ek)        
        ek = ek[indices]
        C = C[:, indices]

        self.ek, self.C = ek, C.T

    def calc_energy(self):
        """
        Calculates the ground state energy from the 
        current coefficient matric.
        """
        n, N = self.n, self.N
        C = self.C

        e = 0
        for p in range(N):
         for a in range(n):
           for b in range(n):
             e += C[p,a]*C[p,b]*self.h_0(a,b)
             for q in range(N):
              for c in range(n):
               for d in range(n):
                e += 0.5*C[p,a]*C[q,b]*C[p,c]*C[q,d]*self.h_1(a,b,c,d)

        self.E = e.subs(Z, Z_value).evalf()
        return self.E

    def solve(self, tol=1e-6, max_iters=40):
        iterations = 0
        n, N = self.n, self.N
        Ep = 0
        ekp = array((0,0,0,0,0,0))

        while iterations < max_iters:
            iterations +=1

            # Find eigenvalues and eigenvector of HF matrix
            self.ek, self.C = linalg.eig(self.h_HF)
            
            # Reorder eigenvalues and eigenvector
            self.reorder_coefficients()

            # Assemble the new HF matrix
            self.assemble_HF_matrix()
            
            # Test tolerance of lowest eigenvalue
            print self.calc_energy()
            error = sum(abs(ekp - self.ek[0]))
            if error < tol:
                print "Solver converged after %d iterations." % (iterations)
                return

            Ep = self.E
            ekp = self.ek[0]

        print "Solver failed to converge in %d iterations." % (iterations)


N = 4
Z_value = 4
basis = [(1,1), (1,-1), (2,1), (2,-1), (3,1), (3,-1)]

print "Solving with initial guess C=I."
solver = HF(N, basis, Z_value, first_C='identity')
solver.solve(max_iters=100)

print "\n\n\n Solving with initial guess C=0."
solver = HF(N, basis, Z_value, first_C='zero')
solver.solve(max_iters=100)

print "\n\n\n Solving with random initial guess."
solver = HF(N, basis, Z_value, first_C='rand')
solver.solve(max_iters=100)

We conclude by listing the matrix elements for the radial integrals to be used for the direct part and the exchange part. Note again that these integrals do not include spin. The nomenclature is \( 1=1s \), \( 2=2s \), and \( 3=3s \), with no spin degrees of freedom.

• 1