8.4: Perturbative Approaches

The situation becomes much more difficult if we need to account for direct interactions between the particles. Let us assume that the interaction may be reduced to that between their pairs (as it is the case at their Coulomb interaction and most other interactions ${ }^{30}$ ), so that it may be described by the following "pair-interaction" Hamiltonian $$\hat{U}_{\text {int }}=\frac{1}{2} \sum_{k, k^{\prime}=1 \atop k \neq k^{\prime}}^{N} \hat{u}_{\text {int }}\left(\mathbf{r}_{k}, \mathbf{r}_{k^{\prime}}\right),$$ with the front factor of $1 / 2$ compensating the double-counting of each particle pair. The translation of this Pairoperator to the second-quantization form may be done absolutely similarly to the derivation of Eq. (87), Hamiltonian and gives a similar (though naturally more involved) result $$\hat{U}_{\text {int }}=\frac{1}{2} \sum_{j, j^{\prime}, l, l^{\prime}} u_{j j^{\prime \prime} l^{\prime}} \hat{a}_{j}^{\dagger} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{l^{\prime}} \hat{a}_{l},$$ where the two-particle matrix elements are defined similarly to Eq. (82): $$u_{i j^{\prime \prime} l^{\prime}} \equiv\left\langle\beta_{j} \beta_{j^{\prime}}\left|\hat{u}_{\mathrm{int}}\right| \beta_{l} \beta_{l^{\prime}}\right\rangle .$$ The only new feature of Eq. (114) is a specific order of the indices of the creation operators. Note the mnemonic rule of writing this expression, similar to that for Eq. (87): each term corresponds to moving a pair of particles from states $l$ and $l^{\prime}$ to states $j$ ’ and $j$ (in this order!) factored with the corresponding two-particle matrix element (115).

However, with the account of such term, the resulting Heisenberg equations of the time evolution of the creation/annihilation operators are nonlinear, so that solving them and calculating observables from the results is usually impossible, at least analytically. The only case when some general results may be obtained is the weak interaction limit. In this case, the unperturbed Hamiltonian contains only single-particle terms such as (79), and we can always (at least conceptually :-) find such a basis of orthonormal single-particle states $\beta_{j}$ in which that Hamiltonian is diagonal in the Dirac representation:

$$\hat{H}^{(0)}=\sum_{j} \varepsilon_{j}^{(0)} \hat{a}_{j}^{\dagger} \hat{a}_{j} .$$ Now we can use Eq. (6.14), in this basis, to calculate the interaction energy as a first-order perturbation: $$\begin{aligned} E_{\mathrm{int}}^{(1)} &=\left\langle N_{1}, N_{2}, \ldots\left|\hat{U}_{\mathrm{int}}\right| N_{1}, N_{2}, \ldots\right\rangle=\frac{1}{2}\left\langle N_{1}, N_{2}, \ldots\left|\sum_{j, j^{\prime}, l, l^{\prime}} u_{j \prime \prime l^{\prime}} \hat{a}_{j}^{\dagger} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{l^{\prime}} \hat{a}_{l}\right| N_{1}, N_{2}, \ldots\right\rangle \\ & \equiv \frac{1}{2} \sum_{j, j^{\prime}, l, l^{\prime}} u_{i j l l^{\prime}}\left\langle N_{1}, N_{2}, \ldots\left|\hat{a}_{j}^{\dagger} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{l^{\prime}} \hat{a}_{l}\right| N_{1}, N_{2}, \ldots\right\rangle . \end{aligned}$$ Since, according to Eq. (63), the Dirac states with different occupancies are orthogonal, the last long bracket is different from zero only for three particular subsets of its indices:

(i) $j \neq j^{\prime}, l=j$, and $l^{\prime}=j^{\prime} .$ In this case, the four-operator product in Eq. (117) is equal to $\hat{a}_{j}^{\dagger} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{j^{\prime}} \hat{a}_{j}$, and applying the commutation rules twice, we can bring it to the so-called normal ordering, with each creation operator standing to the right of the corresponding annihilation operator, thus forming the particle number operator $(72)$ : $$\hat{a}_{j}^{\dagger} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{j^{\prime}} \hat{a}_{j}=\pm \hat{a}_{j}^{\dagger} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{j} \hat{a}_{j^{\prime}}=\pm \hat{a}_{j}^{\dagger}\left(\pm \hat{a}_{j} \hat{a}_{j^{\prime}}^{\dagger}\right) \hat{a}_{j^{\prime}}=\hat{a}_{j}^{\dagger} \hat{a}_{j} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{j^{\prime}}=\hat{N}_{j} \hat{N}_{j^{\prime}},$$ with a similar sign of the final result for bosons and fermions.

(ii) $j \neq j^{\prime}, l=j^{\prime}$, and $l^{\prime}=j$. In this case, the four-operator product is equal to $\hat{a}_{j}^{\dagger} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{j} \hat{a}_{j^{\prime}}$, and bringing it to the form $\hat{N}_{j} \hat{N}_{j^{\prime}}$, requires only one commutation: $$\hat{a}_{j}^{\dagger} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{j} \hat{a}_{j^{\prime}}=\hat{a}_{j}^{\dagger}\left(\pm \hat{a}_{j} \hat{a}_{j^{\prime}}^{\dagger}\right) \hat{a}_{j^{\prime}}=\pm \hat{a}_{j}^{\dagger} \hat{a}_{j} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{j^{\prime}}=\pm \hat{N}_{j} \hat{N}_{j^{\prime}},$$ with the upper sign for bosons and the lower sign for fermions.

(iii) All indices are equal to each other, giving $\hat{a}_{j}^{\dagger} \hat{a}_{j^{\prime}}^{\dagger} \hat{a}_{l^{\prime}} \hat{a}_{l}=\hat{a}_{j}^{\dagger} \hat{a}_{j}^{\dagger} \hat{a}_{j} \hat{a}_{j} .$ For fermions, such an operator (that "tries" to create or to kill two particles in a row, in the same state) immediately gives the null-vector. In the case of bosons, we may use Eq. (74) to commute the internal pair of operators, getting $$\hat{a}_{j}^{\dagger} \hat{a}_{j}^{\dagger} \hat{a}_{j} \hat{a}_{j}=\hat{a}_{j}^{\dagger}\left(\hat{a}_{j} \hat{a}_{j}^{\dagger}-\hat{I}\right) \hat{a}_{j}=\hat{N}_{j}\left(\hat{N}_{j}-\hat{I}\right) .$$ Note, however, that this expression formally covers the fermion case as well (always giving zero). As a result, Eq. (117) may be rewritten in the following universal form: $$E_{\mathrm{int}}^{(1)}=\frac{1}{2} \sum_{j, j^{\prime} \atop j \neq j^{\prime}} N_{j} N_{j^{\prime}}\left(u_{j j i j^{\prime} \prime^{\prime}} \pm u_{j j j^{\prime} j}\right)+\frac{1}{2} \sum_{j} N_{j}\left(N_{j}-1\right) u_{j i j i}$$ The corollaries of this important result are very different for bosons and fermions. In the former case, the last term usually dominates, because the matrix elements (115) are typically the largest when all basis functions coincide. Note that this term allows a very simple interpretation: the number of the diagonal matrix elements it sums up for each state $(j)$ is just the number of interacting particle pairs residing in that state. In contrast, for fermions the last term is zero, and the interaction energy is proportional to the difference of the two terms inside the first parentheses. To spell them out, let us consider the case when there is no direct spin-orbit interaction. Then the vectors $|\beta\rangle_{j}$ of the single-particle state basis may be represented as direct products $\left|o_{j}\right\rangle \otimes\left|m_{j}\right\rangle$ of their orbital and spin-orientation parts. (Here, for the brevity of notation, I am using $m$ instead of $m_{s}$.) For spin- $1 / 2$ particles, including electrons, $m_{j}$ may equal only either $+1 / 2$ or $-1 / 2$; in this case the spin part of the first matrix element, proportional to $u_{j j} j j j$, equals $$\left\langle m\left|\otimes\left\langle m^{\prime} \| m\right\rangle \otimes\right| m^{\prime}\right\rangle,$$ where, as in the general Eq. (115), the position of a particular state vector in each direct product is encoding the particle’s number. Since the spins of different particles are defined in different Hilbert spaces, we may move their state vectors around to get $$\left\langle m\left|\otimes\left\langle m^{\prime} \| m\right\rangle \otimes\right| m^{\prime}\right\rangle=(\langle m \mid m\rangle)_{1} \times\left(\left\langle m^{\prime} \mid m^{\prime}\right\rangle\right)_{2}=1,$$ for any pair of $j$ and $j$ ’. On the other hand, the second matrix element, $u_{j j^{\prime} j^{\prime} j}$, is factored by $$\left\langle m\left|\otimes\left\langle m^{\prime} \| m^{\prime}\right\rangle \otimes\right| m\right\rangle=\left(\left\langle m \mid m^{\prime}\right\rangle\right)_{1} \times\left(\left\langle m^{\prime} \mid m\right\rangle\right)_{2}=\delta_{m m^{\prime}} .$$ In this case, it is convenient to rewrite Eq. (121) in the coordinate representation, using singleparticle wavefunctions called spin-orbitals $$\psi_{j}(\mathbf{r}) \equiv\left\langle\mathbf{r} \mid \beta_{j}\right\rangle=(\langle\mathbf{r} \mid o\rangle \otimes|m\rangle)_{j} .$$ They differ from the spatial parts of the usual orbital wavefunctions of the type (4.233) only in that their index $j$ should be understood as the set of the orbital-state and the spin-orientation indices. ${ }^{31}$ Also, due to the Pauli-principle restriction of numbers $N_{j}$ to either 0 or 1, Eq. (121) may be also rewritten without the explicit occupancy numbers, with the understanding that the summation is extended only over the pairs of occupied states. As a result, it becomes $$E_{\mathrm{int}}^{(\mathrm{l})}=\frac{1}{2} \sum_{j, j^{\prime} \atop j \neq j^{\prime}} \int d^{3} r \int d^{3} r^{\prime}\left[\begin{array}{l} \psi_{j}^{*}(\mathbf{r}) \psi_{j^{\prime}}^{*}\left(\mathbf{r}^{\prime}\right) u_{\mathrm{int}}\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \psi_{j}(\mathbf{r}) \psi_{j^{\prime}}\left(\mathbf{r}^{\prime}\right) \\ -\psi_{j}^{*}(\mathbf{r}) \psi_{j^{\prime}}^{*}\left(\mathbf{r}^{\prime}\right) u_{\mathrm{int}}\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \psi_{j^{\prime}}(\mathbf{r}) \psi_{j}\left(\mathbf{r}^{\prime}\right) \end{array}\right]$$ In particular, for a system of two electrons, we may limit the summation to just two states $\left(j, j^{\prime}=\right.$ $1,2)$. As a result, we return to Eqs. (39)-(41), with the bottom (minus) sign in Eq. (39), corresponding to the triplet spin states. Hence, Eq. (126) may be considered as the generalization of the direct and exchange interaction balance picture to an arbitrary number of orbitals and an arbitrary total number $N$ of electrons. Note, however, that this formula cannot correctly describe the energy of the singlet spin states, corresponding to the plus sign in Eq. (39), and also of the entangled triplet states. ${ }^{32}$ The reason is that the description of entangled spin states, given in particular by Eqs. (18) and (20), requires linear superpositions of different Dirac states. (A proof of this fact is left for the reader’s exercise.)

Now comes a very important fact: the approximate result (126), added to the sum of unperturbed energies $\varepsilon_{j}^{(0)}$, equals the sum of exact eigenenergies of the so-called Hartree-Fock equation: $^{33}$ $$\begin{aligned} &\left(-\frac{\hbar^{2}}{2 m} \nabla^{2}+u(\mathbf{r})\right) \psi_{j}(\mathbf{r}) \\ &+\sum_{j^{\prime} \neq j} \int\left[\psi_{j^{\prime}}^{*}\left(\mathbf{r}^{\prime}\right) u_{\mathrm{int}}\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \psi_{j}(\mathbf{r}) \psi_{j^{\prime}}\left(\mathbf{r}^{\prime}\right)-\psi_{j^{\prime}}^{*}\left(\mathbf{r}^{\prime}\right) u_{\mathrm{int}}\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \psi_{j^{\prime}}(\mathbf{r}) \psi_{j}(\mathbf{r})\right] d^{3} r^{\prime}=\varepsilon_{j} \psi_{j}(\mathbf{r}), \end{aligned}$$ where $u(\mathbf{r})$ is the external-field potential acting on each particle separately - see the second of Eqs. (79). An advantage of this equation in comparison with Eq. (126) is that it allows the (approximate) calculation of not only the energy spectrum of the system, but also the corresponding spin-orbitals, taking into account their electron-electron interaction. Of course, Eq. (127) is an integro-differential rather than just a differential equation. There are, however, efficient methods of numerical solution of such equations, typically based on iterative methods. One more important practical trick is the exclusion of the filled internal electron shells (see Sec. 3.7) from the explicit calculations, because the shell states are virtually unperturbed by the valence electron effects involved in typical atomic phenomena and chemical reactions. In this approach, the Coulomb field of the shells, described by fixed, pre-calculated, and tabulated pseudo-potentials, is added to that of the nuclei. This approach dramatically cuts the computing resources necessary for systems of relatively heavy atoms, enabling a pretty accurate simulation of electronic and chemical properties of rather complex molecules, with thousands of electrons. ${ }^{34}$ As a result, the Hartree-Fock approximation has become the de-facto baseline of all socalled $a b$-initio ("first-principle") calculations in the very important field of quantum chemistry. ${ }^{35}$

In departures from this baseline, there are two opposite trends. For larger accuracy (and typically smaller systems), several "post-Hartree-Fock methods", notably including the configuration interaction method, ${ }^{36}$ that are more complex but may provide higher accuracy, have been developed.

There is also a strong opposite trend of extending such $a b$ initio ("first-principle") methods to larger systems while sacrificing some of the results’ accuracy and reliability. The ultimate limit of this trend is applicable when the single-particle wavefunction overlaps are small and hence the exchange interaction is negligible. In this limit, the last term in the square brackets in Eq. (127) may be ignored, and the multiplier $\psi_{i}(\mathbf{r})$ taken out of the integral, which is thus reduced to a differential equation formally just the Schrödinger equation for a single particle in the following self-consistent effective potential: $$u_{\mathrm{ef}}(\mathbf{r})=u(\mathbf{r})+u_{\mathrm{dir}}(\mathbf{r}), \quad u_{\mathrm{dir}}(\mathbf{r})=\sum_{j^{\prime} \neq j} \int \psi_{j^{\prime}}^{*}\left(\mathbf{r}^{\prime}\right) u_{\mathrm{int}}\left(\mathbf{r}, \mathbf{r}^{\prime}\right) \psi_{j^{\prime}}\left(\mathbf{r}^{\prime}\right) d^{3} r^{\prime}$$ This is the so-called Hartree approximation - that gives reasonable results for some systems, ${ }^{37}$ especially those with low electron density.

However, in dense electron systems (such as typical atoms, molecules, and condensed matter) the exchange interaction, described by the second term in the square brackets of Eqs. (126)-(127), may be as high as $\sim 30 \%$ of the direct interaction, and frequently cannot be ignored. The tendency of taking this interaction in the simplest possible form is currently dominated by the Density Functional Theory, ${ }^{38}$ universally known by its acronym DFT. In this approach, the equation solved for each eigenfunction $\psi_{j}(\mathbf{r})$ is a differential, Schrödinger-like Kohn-Sham equation $$\left[-\frac{\hbar^{2}}{2 m} \nabla^{2}+u(\mathbf{r})+u_{\mathrm{dir}}^{\mathrm{KS}}(\mathbf{r})+u_{\mathrm{xc}}(\mathbf{r})\right] \psi_{j}(\mathbf{r})=\varepsilon_{j} \psi_{j}(\mathbf{r}),$$ where $$u_{\mathrm{dir}}^{\mathrm{KS}}(\mathbf{r})=-e \phi(\mathbf{r}), \quad \phi(\mathbf{r})=\frac{1}{4 \pi \varepsilon_{0}} \int d^{3} r^{\prime} \frac{\rho\left(\mathbf{r}^{\prime}\right)}{\left|\mathbf{r}-\mathbf{r}^{\prime}\right|}, \quad \rho(\mathbf{r})=-e n(\mathbf{r})$$ and $n(\mathbf{r})$ is the total electron density in a particular point, calculated as $$n(\mathbf{r}) \equiv \sum_{j} \psi_{j}^{*}(\mathbf{r}) \psi_{j}(\mathbf{r}) .$$ The most important feature of the Kohn-Sham Hamiltonian is the simplified description of the exchange and correlation effects by the effective exchange-correlation potential $u_{\mathrm{xc}}(\mathbf{r})$. This potential is calculated in various approximations, most of them valid only in the limit when the number of electrons in the system is very high. The simplest of them (proposed by Kohn et al. in the 1960s) is the Local Density Approximation (LDA) in which the effective exchange potential at each point is a function only of the electron density $n$ at the same point $\mathbf{r}$, taken from the theory of a uniform gas of free electrons. ${ }^{39}$ However, for many tasks of quantum chemistry, the accuracy given by the LDA is insufficient, because inside molecules the density $n$ typically changes very fast. As a result, DFT has become widely accepted in that field only after the introduction, in the $1980 \mathrm{~s}$, of more accurate, though more cumbersome models for $u_{\mathrm{xc}}(\mathbf{r})$, notably the so-called Generalized Gradient Approximations (GGAs). Due to its relative simplicity, DFT enables the calculation, with the same computing resources and reasonable precision, some properties of much larger systems than the methods based on the Hartree-Fock theory. As the result, is has become a very popular tool of ab initio calculations. This popularity is enhanced by the availability of several advanced DFT software packages, some of them in the public domain.

Please note, however, that despite this undisputable success, this approach has its problems. From my personal point of view, the most offensive of them is the implicit assumption of the unphysical Coulomb interaction of an electron with itself (by dropping, on the way from Eq. (128) to Eq. (130), the condition $j^{\prime} \neq j$ at the calculation of $u_{\text {dir }}^{\mathrm{KS}}$ ). As a result, for a reasonable description of some effects, the available DFT packages are either inapplicable at all or require substantial artificial tinkering. ${ }^{40}$

Unfortunately, because of lack of time/space, for details I have to refer the interested reader to specialized literature. ${ }^{41}$

---

${ }^{30}$ A simple but important example from the condensed matter theory is the so-called Hubbard model, in which particle repulsion limits their number on each of localized sites to either 0 , or 1 , or 2 , with negligible interaction of the particles on different sites - though the next-neighbor sites are still connected by tunneling - as in Eq. (112).

${ }^{31}$ The spin-orbitals (125) are also close to spinors (13), besides that the former definition takes into account that the spin $s$ of a single particle is fixed, so that the spin-orbital may be indexed by the spin’s orientation $m \equiv m_{s}$ only. Also, if an orbital index is used, it should be clearly distinguished from $j$, i.e. the set of the orbital and spin indices. This is why I believe that the frequently met notation of spin-orbitals as $\psi_{j, s}(\mathbf{r})$ may lead to confusion.

${ }_{32}$ Indeed, due to the condition $j^{\prime} \neq j$, and Eq. (124), the calculated negative exchange interaction is limited to electron state pairs with the same spin direction - such as the factorable triplet states (\uparrow\uparrow and $\downarrow \downarrow$ ) of a twoelectron system, in which the contribution of $E_{\mathrm{ex}}$, given by Eq. (41), to the total energy is also negative.

${ }^{33}$ This equation was suggested in 1929 by Douglas Hartree for the direct interaction and extended to the exchange interaction by Vladimir Fock in 1930. To verify its compliance with Eq. (126), it is sufficient to multiply all terms of Eq. (127) by $\psi^{*}(\mathbf{r})$, integrate them over all $\mathbf{r}$-space (so that the right-hand side would give $\left.\varepsilon_{j}\right)$, and then sum these single-particle energies over all occupied states $j$.

${ }^{34}$ For condensed-matter systems, this and other computational methods are applied to single elementary spatial cells, with a limited number of electrons in them, using cyclic boundary conditions.

${ }^{35}$ See, e.g., A. Szabo and N. Ostlund, Modern Quantum Chemistry, Revised ed., Dover, $1996 .$

${ }^{36}$ That method, in particular, allows the calculation of proper linear superpositions of the Dirac states (such as the entangled states for $N=2$, discussed above) which are missing in the generic Hartree-Fock approach - see, e.g., the just-cited monograph by Szabo and Ostlund.

${ }^{37}$ An extreme example of the Hartree approximation is the Thomas-Fermi model of heavy atoms (with $Z>>1$ ), in which atomic electrons, at each distance $r$ from the nucleus, are treated as an ideal, uniform Fermi gas, with a certain density $n(r)$ corresponding to the local value $u_{\mathrm{ef}}(r)$, but a global value of their highest full single-particle energy, $\varepsilon=0$, to ensure the equilibrium. (The analysis of this model is left for the reader’s exercise.)

${ }^{38}$ It had been developed by Walter Kohn and his associates (notably Pierre Hohenberg) in 1965-66, and eventually (in 1998) was marked with a Nobel Prize in Chemistry for W. Kohn.

${ }_{39}$ Just for the reader’s reference: for a uniform, degenerate Fermi-gas of electrons (with the Fermi energy $\varepsilon_{\mathrm{F}}$ >> $\left.k_{\mathrm{B}} T\right)$, the most important, exchange part $u_{\mathrm{x}}$ of $u_{\mathrm{xc}}$ may be calculated analytically: $u_{\mathrm{x}}=-(3 / 4 \pi) e^{2} k_{\mathrm{F}} / 4 \pi \varepsilon_{0}$, where the Fermi momentum $k_{\mathrm{F}}=\left(2 m_{\mathrm{e}} \varepsilon_{\mathrm{F}}\right)^{1 / 2} / \hbar$ is defined by the electron density: $n=2(4 \pi / 3) k_{\mathrm{F}}^{3} /(2 \pi)^{3} \equiv k_{\mathrm{F}}^{3} / 3 \pi^{2}$.

${ }^{40}$ As just a few examples, see N. Simonian et al., J. Appl. Phys. 113, 044504 (2013); M. Medvedev et al., Science 335,49 (2017); A. Hutama et al., J. Phys. Chem. C 121, 14888 (2017).

${ }^{41}$ See, e.g., either the monograph by R. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford U. Press, 1994, or the later textbook J. A. Steckel and D. Sholl, Density Functional Theory: Practical Introduction, Wiley, 2009. For a popular review and references to more recent work in this still-developing field, see A. Zangwill, Phys. Today $\mathbf{6 8}, 34$ (July 2015).