8 Hartree Method

8.1 Born–Oppenheimer Approximation: Concluding Remarks

In the framework of first-order perturbation theory, the off-diagonal nuclear kinetic energy terms will not contribute to the total energy; contributions from these terms arise only in second-order perturbation theory, and are generally quite small. Thus, in the Born–Oppenheimer approximation, we neglect the off-diagonal transitions. In other words, we assume that nuclear motion cannot cause electronic transitions and the electrons always remain in their ground state.

We see that in Equation 2.8, $E^n_E(\vec{R}) + V_{II}(\vec{R})$ is the potential in which the ions move. It defines a surface in a space spanned by nuclear coordinates, which is called Born–Oppenheimer potential energy surface (BOPES). For each $n$, there is a BOPES. Normally, we assume $n = 0$, that is, the electrons are in the ground state. Thus, if the dynamics are adiabatic, the electrons will stay on the BOPES given by $E^0_E(\vec{R}) + V_{II}(\vec{R})$ and will not make any transitions to other surfaces corresponding to $n \neq 0$.

The Born–Oppenheimer approximation is an excellent approximation for semiconductors and insulators as there is an energy gap in the electronic spectrum. This gap is of the order of electron volt, whereas the typical phonon energies are of the order of 0.01 eV. Thus, in these materials, the nuclear (ionic) motion cannot excite electrons to higher energies and consequently the off-diagonal terms are zero. On the other hand, in metallic systems, the energy gap between the highest occupied and lowest unoccupied states is zero and ionic motion can cause transitions of vanishing energy. In that case, the Born–Oppenheimer approximation can break down. Thus, in first-principles simulations of the metallic system, one has to be very watchful.

8.2 Hartree Method

In the previous section, we saw that the electronic motion can be separated from the ionic motion using the Born–Oppenheimer approximation. Now, we shall focus on solving the electronic problem and return to the ionic problem later. To do this, we shall assume that the ions are fixed at their positions and try to solve Equation 2.4, which can be written as

\[ \begin{aligned} \Bigg( -\frac{\hbar^2}{2m_e}\sum_i \nabla_i^2 + V_{IE} + V_{EE} \Bigg) \psi(&\vec{r}_1, \sigma_1,\, \vec{r}_2, \sigma_2,\, \ldots, \\ &\vec{r}_N, \sigma_N) \\ &= E\psi( \vec{r}_1, \sigma_1,\, \vec{r}_2, \sigma_2,\, \ldots,\, \vec{r}_N, \sigma_N) \end{aligned} \tag{8.1}\]

For simplicity, we have dropped the subscript $E$ of energy $E$ of Equation 2.4. Equation Eq. 8.1 is a complicated equation for many-electron wave function $\psi$, which is a function of all electronic coordinates (position as well as spin). The simplest approximation is due to Hartree (1928) who wrote it as the product of single-particle functions as

\[ \psi_H = \psi_1(\vec{r}_1, \sigma_1)\,\psi_2(\vec{r}_2, \sigma_2) \cdots \psi_N(\vec{r}_N, \sigma_N) \tag{8.2}\]

where $\psi_i(\vec{r}, \sigma)$ denotes the wave function of a single electron $i$ with spin $\sigma$ and is also called orbital. It has the usual interpretation that $|\psi_i|^2\,d^3r$ gives the probability of finding the electron with spin $\sigma$ in volume element $d^3r$.

In this section, we shall be neglecting the spin of the electrons as it plays no role in the Hartree method. Now, we shall use the variational principle for ground-state energy to find $E$ and $\psi_i$. The variational principle states that the expectation value of the energy

\[ E = \langle \psi_H \,|\, H \,|\, \psi_H \rangle \geq E_0 \tag{8.3}\]

where $E_0$ is the ground-state energy of the system. We shall assume that $\psi_i$ are orthonormal. Now, substituting Equation Eq. 8.2 and

\[ H = -\sum_i \frac{\hbar^2}{2m_e} \nabla_i^2 + V_{IE} + V_{EE} \]

in Equation Eq. 8.3, where

\[ V_{IE} = -\frac{1}{4\pi\varepsilon_0} \sum_i \sum_\ell \frac{Ze^2}{\left|\vec{r}_i - \vec{R}_\ell\right|} \]

and

\[ V_{EE} = \frac{1}{8\pi\varepsilon_0} \sum_{\substack{ij \\ i \neq j}} \frac{e^2}{\left|\vec{r}_i - \vec{r}_j\right|} \]

we get

\[ \begin{aligned} E ={}& -\sum_i \int \psi_i^*\, \frac{\hbar^2}{2m_e} \nabla_i^2\, \psi_i\, d^3r_i \\[6pt] &- \frac{1}{4\pi\varepsilon_0} \sum_i \sum_\ell \int \psi_i^*\, \frac{Ze^2} {\left|\vec{r}_i-\vec{R}_\ell\right|} \, \psi_i\, d^3r_i \\[6pt] &+ \frac{1}{8\pi\varepsilon_0} \sum_i \sum_{\substack{j \\ i\neq j}} \int \psi_i^*\psi_j^*\, \frac{e^2} {\left|\vec{r}_i-\vec{r}_j\right|} \, \psi_i\psi_j\, d^3r_i\, d^3r_j \end{aligned} \tag{8.4}\] where we have assumed that $\psi_i$ are orthonormal. We note that $E$ is a functional of $\psi$ (see Appendix D).

Now, we vary the wave functions $\psi_i$ so as to minimize $E$ subject to constraint

\[ \int \psi_i^*\, \psi_i\, d^3r_i = 1 \tag{8.5}\]

A more general treatment would be to demand orthonormality during the variational procedure, but this would add more constraint equations and complicate the algebra. However, the above treatment using Equation Eq. 8.5 also gives the same result. This is done, most easily, by using the Lagrange multiplier method. Thus, we minimize

\[ E - \sum_i \varepsilon_i \left[ \int \psi_i^*\, \psi_i\, d^3r_i - 1 \right] \tag{8.6}\]

without any constraint. Here, $\varepsilon_i$ are the Lagrange multipliers that are yet to be determined. Taking the functional derivative of Equation Eq. 8.6 (see Appendix D), we get

\[ \frac{\delta E}{\delta \psi_i^*} - \varepsilon_i \psi_i = 0 \]

where $\delta E / \delta \psi_i^*$ is the functional derivative of $E$.

Thus, we get \[ \begin{aligned} & -\frac{\hbar^2}{2m_e} \nabla_i^2\, \psi_i \\[6pt] &- \frac{1}{4\pi\varepsilon_0} \sum_\ell \frac{Ze^2} {\left|\vec{r}_i-\vec{R}_\ell\right|} \,\psi_i \\[6pt] &+ \frac{1}{4\pi\varepsilon_0} \sum_{j\neq i} \int \frac{e^2\,|\psi_j|^2} {\left|\vec{r}_i-\vec{r}_j\right|} \, d^3r_j\; \psi_i \\[6pt] &= \varepsilon_i \psi_i \end{aligned} \tag{8.7}\] Equation Eq. 8.7 is known as the Hartree equation and by solving this equation, we obtain all $\varepsilon_i$ and $\psi_i$. Now, we populate the lowest $N$ levels in accordance with the Pauli exclusion principle at absolute zero temperature. The treatment can be generalized to finite temperature by using the Fermi–Dirac statistics.

Now, the number density $n(\vec{r})$ can be written as

\[ n(\vec{r}) = \sum_{\text{occupied}} \psi_i^*\psi_i = \sum_{\text{occupied}} |\psi_i|^2 \tag{8.8}\]

Also, we define Hartree potential as

\[ v_H(\vec{r}_i) = \frac{1}{4\pi\varepsilon_0} \sum_{j \neq i} \int \frac{e^2\,|\psi_j|^2}{\left|\vec{r}_i - \vec{r}_j\right|}\, d^3r_j \]

and

\[ v_I(\vec{r}_i) = -\frac{1}{4\pi\varepsilon_0} \sum_\ell \frac{Ze^2}{\left|\vec{r}_i - \vec{R}_\ell\right|} \tag{8.9}\]

The Hartree potential can be rewritten as

\[ \begin{aligned} v_H(\vec{r}_i) &= \frac{1}{4\pi\varepsilon_0} \iint \sum_{j \neq i} \frac{e^2\,|\psi_j|^2}{\left|\vec{r}_i - \vec{r}\right|} \,\delta(\vec{r} - \vec{r}_j)\, d^3r\, d^3r_j \\[6pt] &= \frac{1}{4\pi\varepsilon_0} \int \frac{d^3r}{\left|\vec{r}_i - \vec{r}\right|} \, e^2 \sum_{j \neq i} |\psi_j(\vec{r})|^2 \end{aligned} \tag{8.10}\]

We see that the Hartree potential is the average electrostatic potential created by the rest of the electrons. For a system containing a large number of electrons such as a solid, the $j = i$ term in Equation Eq. 8.10 will have negligible contribution in comparison to the total contribution. For such systems, $v_H(\vec{r})$ can be approximated by

\[ v_H(\vec{r}) \simeq \frac{1}{4\pi\varepsilon_0} \int \frac{e^2\, n(\vec{r}')}{\left|\vec{r}' - \vec{r}\right|}\, d^3r' \tag{8.11}\]

Thus, the Hartree equation can be written as

\[ \left[ -\frac{\hbar^2}{2m_e} \nabla^2 + v_I(\vec{r}) + v_H(\vec{r}) \right] \psi_i(\vec{r}) = \varepsilon_i\, \psi_i(\vec{r}) \tag{8.12}\]

Equation Eq. 8.12 has a very familiar form and is similar to Schrödinger’s equation for a single electron moving in a potential $v_I(\vec{r}) + v_H(\vec{r})$. However, there is a big difference: in the Hartree equation, the potential depends on the solution itself (see Equation Eq. 8.10). The way this equation is solved is as follows. First, a guess is made for $v_H(\vec{r})$, and Equation Eq. 8.12 is solved. Then, using the solutions, $v_H(\vec{r})$ is computed and the equation is solved again and the potential is computed. This procedure is repeated until the input and the output potentials are very close, within some tolerance. This procedure is called a self-consistent procedure.

Under assumption (Eq. 8.2), the total wave function can be written as a simple product of single-particle wave functions and leads to Equation Eq. 8.12, which is an equation for a single electron moving in an effective potential $v_I(\vec{r}) + v_H(\vec{r})$. Thus, we see that the Hartree approximation reduces the many-electron problem to a one-electron problem. This kind of approximation is also known as one-electron approximation or independent particle approximation. We shall see that this is the common philosophy followed in the Hartree, H–F, and the DFT. In all such theories, the problem is reduced to the equation of the form given in Equation Eq. 8.12. One-electron approximation has been extremely successful in understanding and predicting a materials’ properties.

Although one-electron approximation is very successful, it neglects correlations between electrons that arise due to the Coulomb repulsion between electrons. Because of the Coulomb repulsion, electrons move in such a way as to avoid each other as far as possible. For example, if one electron moves, it affects the motion of all other electrons. Thus, owing to the Coulomb repulsion, the motion of electrons is coupled or correlated. Such correlations between electronic motions are called Coulomb or electron correlations. But in one-electron theories, an electron moves like an independent particle and thus, such theories miss out the effect of correlated motion or electron correlations. Thus, electron correlation is a many-body effect, and to understand it, one has to go beyond one-electron theories.

From Equation Eq. 8.4, total energy can be written as

\[ E = \sum_i \varepsilon_i - \frac{e^2}{2} \sum_{\substack{ij \\ i \neq j}} \iint \frac{1}{4\pi\varepsilon_0} \frac{|\psi_i|^2\,|\psi_j|^2}{\left|\vec{r}_i - \vec{r}_j\right|} \, d^3r_i\, d^3r_j \tag{8.13}\]

Here, the sum is only over occupied levels ($T = 0$). Using Equation Eq. 8.9, the second term in Equation Eq. 8.4 can be written as

\[ \begin{aligned} \int \sum_i v_I(\vec{r}_i)\,|\psi_i(\vec{r}_i)|^2\, d^3r_i &= \iint \sum_i v_I(\vec{r}_i)\,\delta(\vec{r} - \vec{r}_i)\,|\psi_i(\vec{r}_i)|^2\, d^3r_i\, d^3r \\[6pt] &= \int \sum_i v_I(\vec{r})\,|\psi_i(\vec{r})|^2\, d^3r \end{aligned} \tag{8.14}\]

\[ = \int v_I(\vec{r})\, n(\vec{r})\, d^3r \tag{8.15}\]

Although the Hartree approximation introduces great simplification and reduces the many-body problem to a one-body problem, it has two major shortcomings. First, it includes the electron–electron interaction only in an average manner. Second, the many-body Hartree wave function violates the antisymmetry requirement.

8.2.1 Interpretation of $\varepsilon_i$

We had introduced $\varepsilon_i$ in Equation Eq. 8.6 as a Lagrange multiplier. To interpret it physically, let us remove the $k$th electron from the system and calculate the total energy. We assume that by removing one electron from the system, the orbitals $\psi_i$ remain unchanged. The assumption is valid if the system is very large such as a solid.

Let us call the energy of $(N - 1)$ electrons as $E'$. Then, from Equation Eq. 8.4, $E'$ can be written as

\[ \begin{aligned} E' &= -\sum_{i \neq k} \int \psi_i^*\, \frac{\hbar^2}{2m_e} \nabla_i^2\, \psi_i\, d^3r_i - \sum_{i \neq k} \sum_\ell \int \psi_i^*\, \frac{1}{4\pi\varepsilon_0} \frac{Ze^2}{\left|\vec{r}_i - \vec{R}_\ell\right|} \,\psi_i\, d^3r_i \\[6pt] &\quad + \frac{1}{2} \sum_{\substack{i \neq k \\ j \neq k \\ i \neq j}} \iint \psi_i^*\psi_j^*\, \frac{1}{4\pi\varepsilon_0} \frac{e^2}{\left|\vec{r}_i - \vec{r}_j\right|} \,\psi_i\psi_j\, d^3r_i\, d^3r_j \end{aligned} \tag{8.16}\]

\[ \begin{aligned} E - E' &= -\frac{\hbar^2}{2m_e} \int \psi_k^*\, \nabla_k^2\, \psi_k\, d^3r_k - \sum_\ell \int \psi_k^*\, \frac{1}{4\pi\varepsilon_0} \frac{Ze^2}{\left|\vec{r}_k - \vec{R}_\ell\right|} \,\psi_k\, d^3r_k \\[6pt] &\quad + \frac{1}{2} \sum_{j \neq k} \iint \psi_k^*\psi_j^*\, \frac{1}{4\pi\varepsilon_0} \frac{e^2}{\left|\vec{r}_k - \vec{r}_j\right|} \,\psi_k\psi_j\, d^3r_k\, d^3r_j \\[6pt] &\quad + \frac{1}{2} \sum_{i \neq k} \iint \psi_i^*\psi_k^*\, \frac{1}{4\pi\varepsilon_0} \frac{e^2}{\left|\vec{r}_i - \vec{r}_k\right|} \,\psi_k\psi_i\, d^3r_k\, d^3r_j \\[6pt] &= -\frac{\hbar^2}{2m_e} \int \psi_k^*\, \nabla_k^2\, \psi_k\, d^3r_k - \sum_\ell \int |\psi_k|^2\, \frac{1}{4\pi\varepsilon_0} \frac{Ze^2}{\left|\vec{r}_k - \vec{R}_\ell\right|} \,d^3r_k \\[6pt] &\quad + \sum_{j \neq k} \iint |\psi_k|^2\,|\psi_j|^2\, \frac{1}{4\pi\varepsilon_0} \frac{e^2}{\left|\vec{r}_k - \vec{r}_j\right|} \,d^3r_k\, d^3r_j \end{aligned} \tag{8.17}\] \[ = \varepsilon_k \tag{8.18}\]

(from Equation Eq. 8.7)

\[ E' - E = -\varepsilon_k \tag{8.19}\]

that is, $-\varepsilon_k$ is the energy required to remove the $k$th electron from the system. Thus, $\varepsilon_k$ can be interpreted as single-particle energy.

8.3 Hartree–Fock (H–F) Method

There are two main problems with the Hartree method:

Since electrons are fermions, the many-body wave function should be antisymmetric with respect to exchange of two electrons. Hartree’s assumption (Equation Eq. 8.2) does not satisfy this requirement.
The electrostatic interaction between electrons is treated in an average way.

In the H–F method (Fock 1930), the first problem is addressed. We want the many-body wave function to satisfy the requirement

\[ \psi(\vec{r}_1\sigma_1, \ldots, \vec{r}_i\sigma_i, \ldots, \vec{r}_j\sigma_j, \ldots) = -\psi(\vec{r}_1\sigma_1, \ldots, \vec{r}_j\sigma_j, \ldots, \vec{r}_i\sigma_i, \ldots) \tag{8.20}\]

If we have our trial wave function as the determinant, we can easily satisfy this requirement

\[ \psi(\vec{r}_1\sigma_1, \ldots, \vec{r}_N\sigma_N) = \frac{1}{\sqrt{N!}} \begin{vmatrix} \psi_1(\vec{r}_1\sigma_1) & \psi_1(\vec{r}_2\sigma_2) & \cdots & \psi_1(\vec{r}_N\sigma_N) \\ \psi_2(\vec{r}_1\sigma_1) & \psi_2(\vec{r}_2\sigma_2) & \cdots & \\ \vdots & \vdots & \ddots & \vdots \\ \psi_N(\vec{r}_1\sigma_1) & \psi_N(\vec{r}_2\sigma_2) & \cdots & \psi_N(\vec{r}_N\sigma_N) \end{vmatrix} \tag{8.21}\]

where $N$ is the number of electrons. The determinantal form of the many-body wave function was first proposed by Slater (1929) and is known as the Slater determinant.

--- title: "Hartree Method" format: html: toc: true toc-depth: 3 toc-title: "Contents" number-sections: true theme: light: cosmo dark: darkly theme-toggle: true highlight-style: github code-fold: false html-math-method: mathjax smooth-scroll: true include-in-header: text: | <style> body { font-size: 1.05rem; line-height: 1.75; } h2 { margin-top: 2.5rem; padding-bottom: 0.4rem; border-bottom: 2px solid var(--bs-border-color); } h3 { margin-top: 1.8rem; color: var(--bs-body-color); } h4 { color: var(--bs-secondary-color); font-style: italic; } .boxed-eq { background: var(--bs-tertiary-bg); border-left: 4px solid #4c6ef5; padding: 0.8rem 1.2rem; border-radius: 4px; margin: 1rem 0; } </style> --- ## Born–Oppenheimer Approximation: Concluding Remarks {#sec-10-bo-remarks} In the framework of first-order perturbation theory, the off-diagonal nuclear kinetic energy terms will not contribute to the total energy; contributions from these terms arise only in second-order perturbation theory, and are generally quite small. Thus, in the Born–Oppenheimer approximation, we neglect the off-diagonal transitions. In other words, we assume that nuclear motion cannot cause electronic transitions and the electrons always remain in their ground state. We see that in Equation 2.8, $E^n_E(\vec{R}) + V_{II}(\vec{R})$ is the potential in which the ions move. It defines a surface in a space spanned by nuclear coordinates, which is called **Born–Oppenheimer potential energy surface (BOPES)**. For each $n$, there is a BOPES. Normally, we assume $n = 0$, that is, the electrons are in the ground state. Thus, if the dynamics are adiabatic, the electrons will stay on the BOPES given by $E^0_E(\vec{R}) + V_{II}(\vec{R})$ and will not make any transitions to other surfaces corresponding to $n \neq 0$. The Born–Oppenheimer approximation is an excellent approximation for semiconductors and insulators as there is an energy gap in the electronic spectrum. This gap is of the order of electron volt, whereas the typical phonon energies are of the order of 0.01 eV. Thus, in these materials, the nuclear (ionic) motion cannot excite electrons to higher energies and consequently the off-diagonal terms are zero. On the other hand, in metallic systems, the energy gap between the highest occupied and lowest unoccupied states is zero and ionic motion can cause transitions of vanishing energy. In that case, the Born–Oppenheimer approximation can break down. Thus, in first-principles simulations of the metallic system, one has to be very watchful. --- ## Hartree Method {#sec-10-hartree} In the previous section, we saw that the electronic motion can be separated from the ionic motion using the Born–Oppenheimer approximation. Now, we shall focus on solving the electronic problem and return to the ionic problem later. To do this, we shall assume that the ions are fixed at their positions and try to solve Equation 2.4, which can be written as $$ \begin{aligned} \Bigg( -\frac{\hbar^2}{2m_e}\sum_i \nabla_i^2 + V_{IE} + V_{EE} \Bigg) \psi(&\vec{r}_1, \sigma_1,\, \vec{r}_2, \sigma_2,\, \ldots, \\ &\vec{r}_N, \sigma_N) \\ &= E\psi( \vec{r}_1, \sigma_1,\, \vec{r}_2, \sigma_2,\, \ldots,\, \vec{r}_N, \sigma_N) \end{aligned} $$ {#eq-10-many-electron-schrodinger} For simplicity, we have dropped the subscript $E$ of energy $E$ of Equation 2.4. Equation @eq-10-many-electron-schrodinger is a complicated equation for many-electron wave function $\psi$, which is a function of all electronic coordinates (position as well as spin). The simplest approximation is due to Hartree (1928) who wrote it as the product of single-particle functions as $$ \psi_H = \psi_1(\vec{r}_1, \sigma_1)\,\psi_2(\vec{r}_2, \sigma_2) \cdots \psi_N(\vec{r}_N, \sigma_N) $$ {#eq-10-hartree-product} where $\psi_i(\vec{r}, \sigma)$ denotes the wave function of a single electron $i$ with spin $\sigma$ and is also called **orbital**. It has the usual interpretation that $|\psi_i|^2\,d^3r$ gives the probability of finding the electron with spin $\sigma$ in volume element $d^3r$. In this section, we shall be neglecting the spin of the electrons as it plays no role in the Hartree method. Now, we shall use the variational principle for ground-state energy to find $E$ and $\psi_i$. The variational principle states that the expectation value of the energy $$ E = \langle \psi_H \,|\, H \,|\, \psi_H \rangle \geq E_0 $$ {#eq-10-variational-principle} where $E_0$ is the ground-state energy of the system. We shall assume that $\psi_i$ are orthonormal. Now, substituting Equation @eq-10-hartree-product and $$ H = -\sum_i \frac{\hbar^2}{2m_e} \nabla_i^2 + V_{IE} + V_{EE} $$ in Equation @eq-10-variational-principle, where $$ V_{IE} = -\frac{1}{4\pi\varepsilon_0} \sum_i \sum_\ell \frac{Ze^2}{\left|\vec{r}_i - \vec{R}_\ell\right|} $$ and $$ V_{EE} = \frac{1}{8\pi\varepsilon_0} \sum_{\substack{ij \\ i \neq j}} \frac{e^2}{\left|\vec{r}_i - \vec{r}_j\right|} $$ we get $$ \begin{aligned} E ={}& -\sum_i \int \psi_i^*\, \frac{\hbar^2}{2m_e} \nabla_i^2\, \psi_i\, d^3r_i \\[6pt] &- \frac{1}{4\pi\varepsilon_0} \sum_i \sum_\ell \int \psi_i^*\, \frac{Ze^2} {\left|\vec{r}_i-\vec{R}_\ell\right|} \, \psi_i\, d^3r_i \\[6pt] &+ \frac{1}{8\pi\varepsilon_0} \sum_i \sum_{\substack{j \\ i\neq j}} \int \psi_i^*\psi_j^*\, \frac{e^2} {\left|\vec{r}_i-\vec{r}_j\right|} \, \psi_i\psi_j\, d^3r_i\, d^3r_j \end{aligned} $$ {#eq-10-energy-functional} where we have assumed that $\psi_i$ are orthonormal. We note that $E$ is a functional of $\psi$ (see Appendix D). Now, we vary the wave functions $\psi_i$ so as to minimize $E$ subject to constraint $$ \int \psi_i^*\, \psi_i\, d^3r_i = 1 $$ {#eq-10-normalization-constraint} A more general treatment would be to demand orthonormality during the variational procedure, but this would add more constraint equations and complicate the algebra. However, the above treatment using Equation @eq-10-normalization-constraint also gives the same result. This is done, most easily, by using the **Lagrange multiplier method**. Thus, we minimize $$ E - \sum_i \varepsilon_i \left[ \int \psi_i^*\, \psi_i\, d^3r_i - 1 \right] $$ {#eq-10-lagrangian} without any constraint. Here, $\varepsilon_i$ are the Lagrange multipliers that are yet to be determined. Taking the functional derivative of Equation @eq-10-lagrangian (see Appendix D), we get $$ \frac{\delta E}{\delta \psi_i^*} - \varepsilon_i \psi_i = 0 $$ where $\delta E / \delta \psi_i^*$ is the functional derivative of $E$. Thus, we get $$ \begin{aligned} & -\frac{\hbar^2}{2m_e} \nabla_i^2\, \psi_i \\[6pt] &- \frac{1}{4\pi\varepsilon_0} \sum_\ell \frac{Ze^2} {\left|\vec{r}_i-\vec{R}_\ell\right|} \,\psi_i \\[6pt] &+ \frac{1}{4\pi\varepsilon_0} \sum_{j\neq i} \int \frac{e^2\,|\psi_j|^2} {\left|\vec{r}_i-\vec{r}_j\right|} \, d^3r_j\; \psi_i \\[6pt] &= \varepsilon_i \psi_i \end{aligned} $$ {#eq-10-hartree-equation} Equation @eq-10-hartree-equation is known as the **Hartree equation** and by solving this equation, we obtain all $\varepsilon_i$ and $\psi_i$. Now, we populate the lowest $N$ levels in accordance with the Pauli exclusion principle at absolute zero temperature. The treatment can be generalized to finite temperature by using the Fermi–Dirac statistics. Now, the number density $n(\vec{r})$ can be written as $$ n(\vec{r}) = \sum_{\text{occupied}} \psi_i^*\psi_i = \sum_{\text{occupied}} |\psi_i|^2 $$ {#eq-10-number-density} Also, we define **Hartree potential** as $$ v_H(\vec{r}_i) = \frac{1}{4\pi\varepsilon_0} \sum_{j \neq i} \int \frac{e^2\,|\psi_j|^2}{\left|\vec{r}_i - \vec{r}_j\right|}\, d^3r_j $$ and $$ v_I(\vec{r}_i) = -\frac{1}{4\pi\varepsilon_0} \sum_\ell \frac{Ze^2}{\left|\vec{r}_i - \vec{R}_\ell\right|} $$ {#eq-10-ion-potential} The Hartree potential can be rewritten as $$ \begin{aligned} v_H(\vec{r}_i) &= \frac{1}{4\pi\varepsilon_0} \iint \sum_{j \neq i} \frac{e^2\,|\psi_j|^2}{\left|\vec{r}_i - \vec{r}\right|} \,\delta(\vec{r} - \vec{r}_j)\, d^3r\, d^3r_j \\[6pt] &= \frac{1}{4\pi\varepsilon_0} \int \frac{d^3r}{\left|\vec{r}_i - \vec{r}\right|} \, e^2 \sum_{j \neq i} |\psi_j(\vec{r})|^2 \end{aligned} $$ {#eq-10-hartree-potential-rewritten} We see that the Hartree potential is the average electrostatic potential created by the rest of the electrons. For a system containing a large number of electrons such as a solid, the $j = i$ term in Equation @eq-10-hartree-potential-rewritten will have negligible contribution in comparison to the total contribution. For such systems, $v_H(\vec{r})$ can be approximated by $$ v_H(\vec{r}) \simeq \frac{1}{4\pi\varepsilon_0} \int \frac{e^2\, n(\vec{r}')}{\left|\vec{r}' - \vec{r}\right|}\, d^3r' $$ {#eq-10-hartree-potential-approx} Thus, the Hartree equation can be written as $$ \left[ -\frac{\hbar^2}{2m_e} \nabla^2 + v_I(\vec{r}) + v_H(\vec{r}) \right] \psi_i(\vec{r}) = \varepsilon_i\, \psi_i(\vec{r}) $$ {#eq-10-hartree-compact} Equation @eq-10-hartree-compact has a very familiar form and is similar to Schrödinger's equation for a single electron moving in a potential $v_I(\vec{r}) + v_H(\vec{r})$. However, there is a big difference: in the Hartree equation, the potential depends on the solution itself (see Equation @eq-10-hartree-potential-rewritten). The way this equation is solved is as follows. First, a guess is made for $v_H(\vec{r})$, and Equation @eq-10-hartree-compact is solved. Then, using the solutions, $v_H(\vec{r})$ is computed and the equation is solved again and the potential is computed. This procedure is repeated until the input and the output potentials are very close, within some tolerance. This procedure is called a **self-consistent procedure**. Under assumption (@eq-10-hartree-product), the total wave function can be written as a simple product of single-particle wave functions and leads to Equation @eq-10-hartree-compact, which is an equation for a single electron moving in an effective potential $v_I(\vec{r}) + v_H(\vec{r})$. Thus, we see that the Hartree approximation reduces the many-electron problem to a one-electron problem. This kind of approximation is also known as **one-electron approximation** or **independent particle approximation**. We shall see that this is the common philosophy followed in the Hartree, H–F, and the DFT. In all such theories, the problem is reduced to the equation of the form given in Equation @eq-10-hartree-compact. One-electron approximation has been extremely successful in understanding and predicting a materials' properties. Although one-electron approximation is very successful, it neglects correlations between electrons that arise due to the Coulomb repulsion between electrons. Because of the Coulomb repulsion, electrons move in such a way as to avoid each other as far as possible. For example, if one electron moves, it affects the motion of all other electrons. Thus, owing to the Coulomb repulsion, the motion of electrons is coupled or correlated. Such correlations between electronic motions are called **Coulomb or electron correlations**. But in one-electron theories, an electron moves like an independent particle and thus, such theories miss out the effect of correlated motion or electron correlations. Thus, electron correlation is a many-body effect, and to understand it, one has to go beyond one-electron theories. From Equation @eq-10-energy-functional, total energy can be written as $$ E = \sum_i \varepsilon_i - \frac{e^2}{2} \sum_{\substack{ij \\ i \neq j}} \iint \frac{1}{4\pi\varepsilon_0} \frac{|\psi_i|^2\,|\psi_j|^2}{\left|\vec{r}_i - \vec{r}_j\right|} \, d^3r_i\, d^3r_j $$ {#eq-10-total-energy} Here, the sum is only over occupied levels ($T = 0$). Using Equation @eq-10-ion-potential, the second term in Equation @eq-10-energy-functional can be written as $$ \begin{aligned} \int \sum_i v_I(\vec{r}_i)\,|\psi_i(\vec{r}_i)|^2\, d^3r_i &= \iint \sum_i v_I(\vec{r}_i)\,\delta(\vec{r} - \vec{r}_i)\,|\psi_i(\vec{r}_i)|^2\, d^3r_i\, d^3r \\[6pt] &= \int \sum_i v_I(\vec{r})\,|\psi_i(\vec{r})|^2\, d^3r \end{aligned} $$ {#eq-10-second-term-transform} $$ = \int v_I(\vec{r})\, n(\vec{r})\, d^3r $$ {#eq-10-second-term-final} Although the Hartree approximation introduces great simplification and reduces the many-body problem to a one-body problem, it has two major shortcomings. First, it includes the electron–electron interaction only in an average manner. Second, the many-body Hartree wave function violates the antisymmetry requirement. --- ### Interpretation of $\varepsilon_i$ {#sec-10-epsilon-interpretation} We had introduced $\varepsilon_i$ in Equation @eq-10-lagrangian as a Lagrange multiplier. To interpret it physically, let us remove the $k$th electron from the system and calculate the total energy. We assume that by removing one electron from the system, the orbitals $\psi_i$ remain unchanged. The assumption is valid if the system is very large such as a solid. Let us call the energy of $(N - 1)$ electrons as $E'$. Then, from Equation @eq-10-energy-functional, $E'$ can be written as $$ \begin{aligned} E' &= -\sum_{i \neq k} \int \psi_i^*\, \frac{\hbar^2}{2m_e} \nabla_i^2\, \psi_i\, d^3r_i - \sum_{i \neq k} \sum_\ell \int \psi_i^*\, \frac{1}{4\pi\varepsilon_0} \frac{Ze^2}{\left|\vec{r}_i - \vec{R}_\ell\right|} \,\psi_i\, d^3r_i \\[6pt] &\quad + \frac{1}{2} \sum_{\substack{i \neq k \\ j \neq k \\ i \neq j}} \iint \psi_i^*\psi_j^*\, \frac{1}{4\pi\varepsilon_0} \frac{e^2}{\left|\vec{r}_i - \vec{r}_j\right|} \,\psi_i\psi_j\, d^3r_i\, d^3r_j \end{aligned} $$ {#eq-10-energy-N-minus-1} $$ \begin{aligned} E - E' &= -\frac{\hbar^2}{2m_e} \int \psi_k^*\, \nabla_k^2\, \psi_k\, d^3r_k - \sum_\ell \int \psi_k^*\, \frac{1}{4\pi\varepsilon_0} \frac{Ze^2}{\left|\vec{r}_k - \vec{R}_\ell\right|} \,\psi_k\, d^3r_k \\[6pt] &\quad + \frac{1}{2} \sum_{j \neq k} \iint \psi_k^*\psi_j^*\, \frac{1}{4\pi\varepsilon_0} \frac{e^2}{\left|\vec{r}_k - \vec{r}_j\right|} \,\psi_k\psi_j\, d^3r_k\, d^3r_j \\[6pt] &\quad + \frac{1}{2} \sum_{i \neq k} \iint \psi_i^*\psi_k^*\, \frac{1}{4\pi\varepsilon_0} \frac{e^2}{\left|\vec{r}_i - \vec{r}_k\right|} \,\psi_k\psi_i\, d^3r_k\, d^3r_j \\[6pt] &= -\frac{\hbar^2}{2m_e} \int \psi_k^*\, \nabla_k^2\, \psi_k\, d^3r_k - \sum_\ell \int |\psi_k|^2\, \frac{1}{4\pi\varepsilon_0} \frac{Ze^2}{\left|\vec{r}_k - \vec{R}_\ell\right|} \,d^3r_k \\[6pt] &\quad + \sum_{j \neq k} \iint |\psi_k|^2\,|\psi_j|^2\, \frac{1}{4\pi\varepsilon_0} \frac{e^2}{\left|\vec{r}_k - \vec{r}_j\right|} \,d^3r_k\, d^3r_j \end{aligned} $$ {#eq-10-E-minus-Eprime} $$ = \varepsilon_k $$ {#eq-10-equals-epsilon-k} (from Equation @eq-10-hartree-equation) or $$ E' - E = -\varepsilon_k $$ {#eq-10-koopmans} that is, $-\varepsilon_k$ is the energy required to remove the $k$th electron from the system. Thus, $\varepsilon_k$ can be interpreted as **single-particle energy**. --- ## Hartree–Fock (H–F) Method {#sec-10-hf} There are two main problems with the Hartree method: 1. Since electrons are fermions, the many-body wave function should be antisymmetric with respect to exchange of two electrons. Hartree's assumption (Equation @eq-10-hartree-product) does not satisfy this requirement. 2. The electrostatic interaction between electrons is treated in an average way. In the H–F method (Fock 1930), the first problem is addressed. We want the many-body wave function to satisfy the requirement $$ \psi(\vec{r}_1\sigma_1, \ldots, \vec{r}_i\sigma_i, \ldots, \vec{r}_j\sigma_j, \ldots) = -\psi(\vec{r}_1\sigma_1, \ldots, \vec{r}_j\sigma_j, \ldots, \vec{r}_i\sigma_i, \ldots) $$ {#eq-10-antisymmetry} If we have our trial wave function as the determinant, we can easily satisfy this requirement $$ \psi(\vec{r}_1\sigma_1, \ldots, \vec{r}_N\sigma_N) = \frac{1}{\sqrt{N!}} \begin{vmatrix} \psi_1(\vec{r}_1\sigma_1) & \psi_1(\vec{r}_2\sigma_2) & \cdots & \psi_1(\vec{r}_N\sigma_N) \\ \psi_2(\vec{r}_1\sigma_1) & \psi_2(\vec{r}_2\sigma_2) & \cdots & \\ \vdots & \vdots & \ddots & \vdots \\ \psi_N(\vec{r}_1\sigma_1) & \psi_N(\vec{r}_2\sigma_2) & \cdots & \psi_N(\vec{r}_N\sigma_N) \end{vmatrix} $$ {#eq-10-slater-determinant} where $N$ is the number of electrons. The determinantal form of the many-body wave function was first proposed by Slater (1929) and is known as the **Slater determinant**.

8.1 Born–Oppenheimer Approximation: Concluding Remarks

8.2 Hartree Method

8.2.1 Interpretation of \(\varepsilon_i\)

8.3 Hartree–Fock (H–F) Method