**1. Introduction**

About a half-century ago, Miron Amusia showed that the photoionization of the noble gas atoms could be described quite accurately using the random-phase approximation with exchange (RPAE); see [1] and references therein. RPAE includes two-particle two-hole correlations in the initial state of the photoionization process and coupling among the ionization channels in the final state, essentially configuration interaction in the continuum. The methodology did not include relativistic interactions. Building upon Amusia's work, Walter Johnson and his co-workers developed a relativistic version, the relativisticrandom-phase approximation (RRPA) which is based on the Dirac equation rather than the Schrödinger equation [2–5]. This advance allowed us to study heavier systems, where relativistic effects are large, and processes that are only possible owing to relativistic interactions. In this paper, an exposition of the random-phase approximation is presented along with some of the important advances using the relativistic version of Amusia's methodology.

## **2. The Random-Phase Approximation**

The random-phase approximation can be formulated in a few alternative ways. An extensive review of these methods and their equivalence is not attempted in this article. Instead, we limit ourselves to comment on the expression "random-phase approximation" and the "linearization" process that it entails. These expressions are widely used in the literature, but the historical reasons behind this terminology are seldom discussed. The scope of this article is further limited to only illustrating some of the applications of the relativistic RPA, RRPA, for which the earlier nonrelativistic work of [1] set the stage.

#### *2.1. Beyond the Hartree–Fock Method*

The simplest N-electron atomic problem, even if we leave out relativistic effects such as the screen-orbit interaction, is described by the Hamiltonian

**Citation:** Deshmukh, P.C.; Manson, S.T. Photoionization of Atomic Systems Using the Random-Phase Approximation Including Relativistic Interactions. *Atoms* **2022**, *10*, 71. https://doi.org/10.3390/ atoms10030071

Academic Editor: Eugene T. Kennedy

Received: 31 May 2022 Accepted: 1 July 2022 Published: 11 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

$$H^{(N)}(q\_1, q\_2, \dots, q\_N) = \left\{ \sum\_{i=1}^N \left( -\frac{\hbar^2}{2m} \nabla\_i^2 - \frac{Ze^2}{r\_i} \right) \right\}$$

$$+ \left\{ \sum\_{i$$

consisting of one- and two- electron operators. The electron–electron Coulomb interaction must be described formally in terms of charge densities which require the electron wavefunctions, which is to say that one needs the solution of the N-electron Schrödinger equation even before the differential equation is formulated. The single electron atomic problem has analytical solutions at both non-relativistic and relativistic levels, but approximation methods are necessary to solve the *N*-electron problem for *N* ≥ 2. The many-electron problem is a vexing one on account of two types of correlations between the electrons: (a) statistical (also called Fermi–Dirac or exchange) correlations and (b) Coulomb correlations. The Hartree–Fock self-consistent field (HF-SCF) provides an effective strategy to obtain solutions that are excellent approximations [6–8] to the required solutions using an *antisymmetrized* product of *N* single-electron wave wavefunctions.

The two-electron terms render the electron–electron potential non-local. The Hartree– Fock method obtains atomic wavefunctions employing numerical procedures, using variational calculus to obtain a self-consistent field in the *frozen orbital approximation*.

Since the wavefunction employed in the HF method is antisymmetrized, statistical correlations are accounted for, but not the *Coulomb correlations*. In fact, the Coulomb correlations are *defined* to be just those that are left out of the HF method. There is no method available to account for the Coulomb correlations exactly.

One must employ approximation methods. Various approximation methods have been developed to address the Coulomb correlations, such as the Multi-Configuration Hartree–Fock (MCHF), also called the Configuration Interaction (CI) method [9], diagrammatic perturbation theory [10], Greens function method [11], etc. A very successful approximation method to treat the electron correlations is the random-phase approximation (RPA). There are various routes to the RPA such as the method of canonical transformations [12], the equation of motion method [13], and the linearized Time-Dependent Hartree–Fock method—*commonly known as the random-phase approximation with exchange*, or RPAE, for short [1,14,15]. All of these routes to RPA are equivalent; they depend on employing a linear approximation to the electron correlations.

We shall first briefly visit salient features of the method of canonical transformation of the Hamiltonian employed by Bohm and Pines [12,16–19] since: (i) their method lucidly illustrates the *linear* approximation to electron correlations and (ii) explains the RPA which involves cancellation of terms associated with the term *random-phase*, arrived at using a *linearization process. This linearization process is the heart of the RPA*. We shall then review the linearization of the Time-Dependent Hartree–Fock system (TDHF) of equations developed by Dalgarno and Victor (1968). This approach is equivalent to that of Bohm–Pines on account of the linearization process that drives it. It is especially insightful toward appreciating the treatment of the many-electron system *beyond* the Hartree–Fock model. The *linearized* TDHF thus provides the essential platform toward methods in the analysis of atomic collisions and photoionization processes [1] employing the technique known as the random-phase-approximation with exchange (RPAE). Finally, we shall summarize the relativistic improvisation of the RPAE, called the relativistic-random-phase approximation (RRPA) developed by Johnson, Lin, and Dalgarno [2–5,14] which is arrived at by linearizing the Time-Dependent Dirac–Hartree–Fock (TDDHF, often abbreviated as TDDF) and illustrate some of its applications.

#### *2.2. Linear approximation to Coulomb correlations*

Strong electron–electron correlations in a free electron gas were dealt with by Bohm and Pines by subjecting the many-electron Hamiltonian to a series of canonical transformations. These transformations result in weakly interacting *elementary excitations* (plasmons) which represent collective elementary excitations of the electron gas.

Using the second quantized notation [20] for the electron creation (*c*†) and annihilation operator (*c*), the Schrödinger equation for an N-electron free electron gas is

$$\frac{\partial}{\partial t}|\Psi(t)\rangle = \left[\sum\_{i}\sum\_{j}c\_{i}^{\dagger}\langle i|f|j\rangle c\_{j}\right]$$

$$\left[ +\frac{1}{2}\sum\_{i}\sum\_{j}\sum\_{k}\sum\_{l}c\_{i}^{\dagger}c\_{j}^{\dagger}\langle ij|v|kl\rangle c\_{l}c\_{k}\right]|\Psi(t)\rangle,\tag{2a}$$

where

$$<\langle ij|v|kl\rangle = \int dq\_1 \int dq\_2 \psi\_i^\*(q\_1)\psi\_j^\*(q\_2)v(q\_1,q\_2)\psi\_k(q\_1)\psi\_l(q\_2),\tag{2b}$$

with the subscripts *i*, *j*, *k*, *l* denoting the set of four one-electron quantum numbers and the arguments (*qr*) denoting the four coordinates of the *rth* electron, three of which being the space coordinates and the fourth being the spin coordinate. The operator f is a *single-electron* operator, similar to that in Equation (1), but consisting of only the kinetic energy terms, the electrons being *free*. In Equation (2b), a typical spin-orbital is represented by

$$
\psi\_i(q) = \psi\_i(\vec{r})\chi\_i(\mathbb{Q}),\tag{2c}
$$

wherein the spin part is either

$$\chi\_i(\mathbb{J}) = \mathfrak{a} = \begin{bmatrix} 1 \\ 0 \end{bmatrix} \, \tag{2d}$$

for *msi*= +12 i.e., ↑, or

$$\chi\_i(\zeta) = \beta = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \,\,\,\,\,\tag{2e}$$

for *msi* = −12 i.e., ↓.

The second quantized Hamiltonian in Equation (2a) is equivalently written as

$$H = \left(\sum\_{i,\alpha} \sum\_{j,\beta} c\_{i\alpha}^{\dagger} \left(\int \psi\_{i\alpha}^{\*}(q) f(q) \psi\_{j\beta}^{\*}(q) dq\right) c\_{j\beta}\right)$$

$$+ \frac{1}{2} \sum\_{i\alpha} \sum\_{j\beta} \sum\_{k\gamma} \sum\_{l\delta} c\_{i\alpha}^{\dagger} c\_{j\beta}^{\dagger}$$

$$\times \left(\int dq dq' \psi\_{i\alpha}^{\*}(q) \psi\_{j\beta}^{\*}(q) \upsilon(q, q') \psi\_{l\delta}(q) \psi\_{k\gamma}(q)\right) c\_{k\gamma} c\_{l\delta}\Big),\tag{3a}$$

or more compactly using the Dirac notation for the integrals as

$$H = \sum\_{i,\alpha} \sum\_{j,\emptyset} c\_{i\alpha}^{\dagger} \langle i\alpha | f | j\beta \rangle c\_{j\beta} + \frac{1}{2} \sum\_{i\alpha} \sum\_{j\emptyset} \sum\_{k\gamma} \sum\_{l\delta} c\_{i\alpha}^{\dagger} c\_{j\beta}^{\dagger}$$

$$\times \langle i\alpha, j\beta | \upsilon | l\delta, k\gamma \rangle c\_{k\gamma} c\_{l\delta}.\tag{3b}$$

Without compromising the above Hamiltonian in any way, we can place the *most part* of the two-electron interactions

$$\left\{ \frac{1}{2} \sum\_{i=1: i \neq j}^{N} \sum\_{j=1}^{N} v\left(\vec{r}\_{i\prime}^{\*} \,\_{j} \vec{r}\_{j}^{\*}\right) \right\}$$

in Equation (1) in a one-electron operator

$$\left\{ \sum\_{i=1}^{N} F(\overrightarrow{r}\_i) \right\}$$

by writing

$$\begin{aligned} H^{(N)}(q\_1, q\_2, \dots, q\_N) &= \sum\_{i=1}^N f(\vec{r}\_i) + \sum\_{i=1}^N F(\vec{r}\_i) \\ &+ \frac{1}{2} \sum\_{i=1; i \neq j}^N \sum\_{j=1}^N v(\vec{r}\_i, \vec{r}\_j) - \sum\_{i=1}^N F(\vec{r}\_i), \end{aligned} \tag{4a}$$

or briefly as

$$H^{(N)}(q\_{1'}q\_{2'}\dots\_{\prime}q\_N) = $$

$$(f+F) + (H\_2 - F) = O\_1 + O^{\prime},\tag{4b}$$

where *O*1 is a one-electron operator and *O* is only a small part of the full Hamiltonian. The notations employed in Equations (4a) and (4b) are self-explanatory. Essentially, the Nelectron Hamiltonian is re-written such that it can be approximated by (*f* + *F*), since

$$\left\{ \frac{1}{2} \sum\_{i=1; i\neq j}^{N} \sum\_{j=1}^{N} \upsilon(\vec{r}\_i, \vec{r}\_j) - \sum\_{i=1}^{N} F(\vec{r}\_i) \right\}$$

is small; we shall treat it perturbatively. The choice of the operator *F* is so made that the total energy functional

$$E^{(N)} = \langle \Psi^{(N)} | H^{(N)} | \Psi^{(N)} \rangle,\tag{5}$$

is minimized.

When *O* is neglected, the unperturbed ground state wavefunction of the N-electron system is expressible as a determinant:

$$\Phi^{(N)} = \frac{1}{\sqrt{N!}} \begin{bmatrix} \Psi\_{\mathbb{T}\uparrow}(1) & \dots & \dots & \Psi\_{\mathbb{T}\uparrow}(N) \\ \Psi\_{\mathbb{T}\downarrow}(1) & \dots & \dots & \Psi\_{\mathbb{T}\downarrow}(N) \\ \dots & \dots & \dots & \dots \\ \Psi\_{\frac{N}{2}\mathbb{T}}(1) & \dots & \dots & \Psi\_{\frac{N}{2}\mathbb{T}}(N) \\ \Psi\_{\frac{N}{2}\mathbb{L}}(1) & \dots & \dots & \Psi\_{\frac{N}{2}\mathbb{L}}(N) \end{bmatrix}\_{N\times N} = $$
 
$$\frac{1}{\sqrt{N!}} \begin{bmatrix} \Psi\_{1}(1) & \dots & \dots & \Psi\_{1}(N) \\ \Psi\_{1}(1) & \dots & \dots & \Psi\_{1}(N) \\ \dots & \dots & \dots & \dots \\ \Psi\_{N-1}(1) & \dots & \dots & \Psi\_{N-1}(N) \\ \Psi\_{N}(1) & \dots & \dots & \Psi\_{N}(N) \end{bmatrix}\_{N\times N} \tag{6}$$

wherein

$$\left[f(\vec{r}) + F(\vec{r})\right]\psi\_{i\sigma}(\vec{r}) = \varepsilon\_i \psi\_{i\sigma}(\vec{r}),\tag{7}$$

with *ψ<sup>i</sup>σ*(*r*) = *ψi*↓(*r*) or *ψi*<sup>↑</sup>(*r*). Observe that the one-particle eigenvalue *εi* is doubly degenerate with the spin-orbitals for spin up and down being linearly independent. In the second determinant in Equation (6), we have only re-designated the single-electron spin-orbitals.

$$\boldsymbol{\Phi}^{(N)} = \frac{1}{\sqrt{N!}} \begin{bmatrix} \psi\_1(1) & \dots & \dots & \dots & \psi\_1(N) \\ \psi\_1(1) & \dots & \dots & \dots & \psi\_1(N) \\ \dots & \dots & \dots & \zeta(|i|) = \psi\_i(q\_j) & \dots \\ \psi\_{N-1}(1) & \dots & \dots & \dots & \psi\_{N-1}(N) \\ \psi\_N(1) & \dots & \dots & \dots & \psi\_N(N) \end{bmatrix}\_{N \times N} \tag{8}$$

Re-designation of the single-electron spin-orbitals in Equation (6) is appropriate for closed shell atoms for which the random-phase approximation is applicable. Wave functions of the *excited* unperturbed states are also Nth order determinants made up of eigenfunctions of Equation (7), but with one or more *εi* > *<sup>ε</sup>N*/2. In an ordered set of spin-orbitals, let us denote *p* ≤ *N* and *q* > *N*. Thus, we denote a typical spin-orbital in the ground state Slater determinant by *p* and an excited state spin-orbital corresponding to a single excitation by *q*. The operator (*f* + *F*) is diagonal with respect to one-electron functions, and *q* = *p*. The choice of the operator *<sup>F</sup>*(*r*) which makes the energy functional (Equation (5)) a minimal is the one for which the matrix element

$$<\langle q|F|p\rangle = \sum\_{i=1}^{N} \left[ \langle iq|v|ip\rangle - \langle qi|v|ip\rangle \right]\_{\prime} \tag{9}$$

as can be shown using a variational method under the *frozen orbital approximation*; i.e., for the excited state we use a Slater determinant only the *pth* spin-orbital replaced by the excited *<sup>q</sup>th*—*all other spin-orbitals in the Slater determinant (8) retain their occupancies*.

The one-electron Hartree–Fock equation satisfied by the SCF ground-state spin-orbitals is

$$\begin{aligned} \left[ \left( -\frac{\hbar^2}{2m}\nabla^2 - \frac{Ze^2}{r} \right) \psi\_p(\vec{\xi}) \right] \\ + \sum\_{i=1}^N \left[ \int d^4 V' \frac{\psi\_i^\* \left( \vec{\xi}' \right) \psi\_i(\vec{\xi}') \psi\_p(\vec{\xi}) e^2}{\vec{r} - \vec{r}'} \right] \\ - \sum\_{i=1}^N \delta(m\_{x\_{p'}}, m\_{s\_i}) \left[ \int d^4 V' \frac{\psi\_i^\* \left( \vec{\xi}' \right) \left( \psi\_p(\vec{\xi}') \psi\_i(\vec{\xi}) \right) e^2}{\vec{r} - \vec{r}'} \right] \\ + \varepsilon\_p \psi\_p(\vec{\xi}). \end{aligned} \tag{10a}$$

The four-dimensional integration *d*4*V* in Equation (10a) includes integration over the three (continuous) space coordinates and the discrete summation over the spin-coordinate. The second and the third terms on the left-hand side of Equation (10a) involve two-electron terms; the second is the Coulomb term and the third is the exchange term. We consider nonferromagnetic systems so that the number of spin-up and -down terms is equal. The energies *ε p* are doubly degenerate with two linearly independent functions corresponding to up- and down-spins. Carrying out the summation over the spin variable, Equation (10a) simplifies to

$$\begin{aligned} \left[ \left( -\frac{\hbar^2}{2m}\nabla^2 - \frac{Ze^2}{r} \right) \psi\_p(\vec{r}) \right] \\ + 2\sum\_{i=1}^{N/2} \left[ \int dV' |\psi\_i(\vec{r'})|^2 \nu(\vec{r}, \vec{r'}) \right] \psi\_p(\vec{r}) \\ - \sum\_{i=1}^{N/2} \psi\_i(\vec{r}) \left[ \int dV' \psi\_i^\*(\vec{r'}) \psi\_p(\vec{r'}) \nu(\vec{r}, \vec{r'}) \right] \right] = \\ \varepsilon\_p \psi\_p(\vec{r}), \end{aligned} \tag{10b}$$

in which we have written the inter-electron Coulomb interaction as

$$\frac{\epsilon^2}{\vec{r} - \vec{r}'} = v(\vec{r}, \vec{r}'). \tag{11}$$

Now,

$$\begin{aligned} F\psi\_p(\vec{r\_2}) &= 0\\ 2\sum\_{i=1}^{N/2} \int d^3 \vec{r\_1} |\psi\_i(\vec{r\_1})|^2 \upsilon(\vec{r\_1}, \vec{r\_2}) \psi\_p(\vec{r\_2})\\ -\sum\_{i=1}^{N/2} \int d^3 \vec{r\_1} \psi\_i^\*(\vec{r\_1}) \upsilon(\vec{r\_1}, \vec{r\_2}) \psi\_p(\vec{r\_1}) \psi\_i(\vec{r\_2}), \end{aligned}$$

hence the HF SCF equation becomes

$$\left[\left(-\frac{\hbar^2}{2m}\nabla^2 - \frac{Ze^2}{r}\right)\psi\_p(\vec{r}) + F\psi\_p(\vec{r}) = \right.$$

$$(f+F)\psi\_p(\vec{r}) = \varepsilon\_p \psi\_p(\vec{r}).\tag{12}$$

We write the momentum-dependent energies as *ε*( *k*) or equivalently as *ε*( *k*), since *p* = *h*¯ *k*. The free electron gas is the only many-electron system for which the HF SCF equation can be obtained analytically. The linearization of the Coulomb approximation was developed by Bohm and Pines [17] in which the positive charges of the nuclei were considered to be spread out uniformly over a volume *V* as jellium (Figure 1). The electron gas is also spread out over the volume *V* in which the electron wavefunction is box-normalized.

**Figure 1.** The uniform spread of the potential generated by the positive nuclei as a jellium.

Adding the jellium potential *<sup>V</sup>*(*r*) in the HF SCF equation we ge<sup>t</sup>

$$\begin{aligned} \left[ -\frac{\hbar^2}{2m}\nabla^2\psi\_p(\vec{r}) + \underbrace{V(\vec{r})\psi\_{\vec{p}}(\vec{r})} \right] \\ + 2\sum\_{i=1}^{N/2} \left[ \int dV' |\psi\_{\vec{r}}(\vec{r'})|^2 \upsilon(\vec{r}, \vec{r'}) \right] \overline{\psi\_p(\vec{r})} \\ - \sum\_{i=1}^{N/2} \psi\_{\vec{r}}(\vec{r}) \left[ \int dV' \psi\_i^\*(\vec{r'}) \psi\_p(\vec{r'}) \upsilon(\vec{r}, \vec{r'}) \right] \right] = \\ \varepsilon\_p \psi\_p(\vec{r}). \end{aligned} \tag{13}$$

The attractive jellium potential (second term on the left-hand side of Equation (13)) exactly cancels the electron–electron Coulomb repulsion term (third term). We are then left with

$$\begin{aligned} \left[ -\frac{\hbar^2}{2m}\nabla^2 - V^{\text{exchange}}(q) \right] \psi\_p(q) &= \\ \varepsilon\_p \psi\_p(\vec{r}), \end{aligned} \tag{14}$$

where

$$V^{\text{exchang}\wp}(q)\psi\_p(q) = \sum\_{i=1}^{N} V\_i^{\text{exchang}\wp}(q)\psi\_p(q) = $$
 
$$\sum\_{i=1}^{N} \psi\_i(q) \left[ \int dq' \frac{\psi\_i^\*(q')\psi\_p(q')}{|\vec{r} - \vec{r'}|} \right], \tag{15a}$$

i.e.,

$$V^{\text{exclungen}}(q)\psi\_p(q) = $$

$$\sum\_{i=1}^{N/2} \psi\_i(\vec{r}) \left[\psi\_i^\*(\vec{r}^\flat)\psi\_p(\vec{r}^\flat)\upsilon(\vec{r},\vec{r}^\flat)\right].\tag{15b}$$

Using (i) box-normalized wavefunctions

$$
\psi\_{\vec{k}\vec{r}}(\vec{r}) = L^{-3/2} \epsilon^{i\vec{k}\cdot\vec{r}} \chi\_{\vec{r}}(\zeta),
\tag{16}
$$

and (ii)

$$\phi(\vec{r\_1}) = \int d^3 \vec{r\_2} \frac{e^{i(\vec{k} - \vec{k'}) \cdot \vec{r\_2}}}{r\_{12}} = \frac{4\pi e^{i(\vec{k} - \vec{k'}) \cdot \vec{r\_1}}}{|\vec{k} - \vec{k'}|^2},\tag{17}$$

the Hartree–Fock equation for the free-electron gas (with exchange) in the positive jellium potential becomes

$$-\frac{\hbar^2}{2m}\nabla^2\psi\_{\vec{k}}(\vec{r}\_1) + \varepsilon\_{\vec{k}}\psi\_{\vec{k}}(\vec{r}\_1) = \varepsilon(\vec{k})\psi\_{\vec{k}}(\vec{r}\_1),\tag{18}$$

with the exchange term given by

$$
\varepsilon\_{\vec{k}} = \frac{4\pi\varepsilon^2}{L^3} \sum\_{\vec{k}'} \frac{1}{|\vec{k} - \vec{k}'|^2}. \tag{19}
$$

It follows that

$$
\varepsilon(\vec{p}) = \frac{p^2}{2m} \varepsilon\_{\text{exchange}}(\vec{p}) , \tag{20}
$$

where

$$\varepsilon\_{\text{exchange}}(\vec{p}) = \frac{-\varepsilon^2 p\_f}{\hbar \pi \tau} \left[ 1 + \frac{p\_f^2 - p^2}{2p\_f p} \ln \left| \frac{p + p\_f}{p - p\_f} \right| \right],\tag{21}$$

*pf* being the electron momentum at the highest occupied energy level; viz., Fermi level. Whereas the HF SCF energy of an atom described by the Hamiltonian (Equation (1)) is

$$E\_{HF}^{atom} = \langle \Psi^{(N)} | H | \Psi^{(N)} \rangle = $$

$$\sum\_{i=1}^{N} \langle i | f | i \rangle + \frac{1}{2} \sum\_{j=1}^{N} \sum\_{i=1}^{N} \left[ \langle ij | \mathbf{g} | ij \rangle - \langle ij | \mathbf{g} | j i \rangle \right], \tag{22}$$

that of an electron gas in the jellium potential of the positive charges is

$$E\_{HF}^{\text{electron gas in jellium potential}} = E\_{KE} + E\_{\text{exchange correlation}} \tag{23}$$

where

$$E\_{KE} = 2\frac{L^3}{\left(2\pi\hbar^2\right)} \int\_{p=0}^{p=\mu\_f} p^2 dp \int\_{\theta=0}^{\theta=\pi} \sin\theta d\theta$$

$$\times \int\_{\phi=0}^{\phi=2\pi} d\phi \left[\frac{\vec{p}\cdot\vec{p}}{2m}\right] = \frac{\hbar^2 L^3}{10\pi^2 m} k\_{f'}^5 \tag{24}$$

with

> *kf* = <sup>3</sup>*π*<sup>2</sup>*NV* 13 , (25)

and

$$E\_{\text{exchange correlation}} = 2\frac{L^3}{(2\pi\hbar^2)}$$

$$\times \int\_{p=0}^{p=p\_f} p^2 dp \int\_{\theta=0}^{\theta=\pi} \sin\theta d\theta \int\_{\phi=0}^{\phi=2\pi} d\phi \left[\frac{1}{2}\varepsilon\_{\text{exchange}}(\vec{p})\right]$$

$$= \frac{\hbar^2 L^3}{10\pi^2 m} k\_{f'}^5 \tag{26a}$$

$$E\_{\text{exchange correlation}} = -\frac{V\varepsilon^2}{4\pi^3} \int\_{k=0}^{k=k\_f} dk \left[2k\_f k^2\right]$$

$$+k(k\_f^2 - k^2) \ln\left(\frac{k\_f + k}{k\_f - k}\right) \,\text{.}\tag{26b}$$

If we now consider the entire physical volume under consideration to consist of spheres of radius *rs* (Seitz parameter in Bohr units), each having *one* unit of electron charge, then

$$N \times \left(\frac{4}{3}\pi r\_s^3\right) = V = \frac{3\pi^2 N}{k\_f^3},\tag{27}$$

then the K.E. contribution to the average HF ground state energy *per electron* in a freeelectron gas is

$$\frac{E\_{KE}}{N} = \frac{3\hbar^2}{10m} \left(\frac{9\pi}{4}\right)^{\frac{2}{3}} \frac{1}{r\_\text{s}^2} = \frac{2.21}{r\_\text{s}^2} \text{Ryd} \tag{28a}$$

and the exchange correlation energy per electron is

$$\frac{E\_{\text{exchange correlation}}}{N} = -\frac{0.916}{r\_s} \text{Ryd.} \tag{28b}$$

Thus, for electron gas in the SCF jellium potential, the average Hartree–Fock energy per electron is

$$\frac{E\_{HF}}{N} = \left(\frac{2.21}{r\_s^2} - \frac{0.916}{r\_s}\right) \text{Ryd.}\tag{29}$$

A first order perturbative treatment gives essentially the same result as above. electron– electron exchange interactions reduce the energy *below* that of the Sommerfeld gas in a positive jellium potential; exchange energy is negative.

In the mid-1950s, Bohm and Pines improvised on the above model by considering a random mutual displacement of the centers of the positive and negative charge densities (Figure 2). In the jellium potential, these are coincident; their mutual displacement can be considered to have been triggered by a random event, but once displaced, the positive and negative charges are set in oscillations of the plasma as the system seeks its original configuration. Bohm and Pines modeled these oscillations using a harmonic oscillator potential, inclusive of a dispersive wavelength-dependent frequency of the plasma oscillations.

The Hamiltonian for *N* electrons in a volume *V* together with a uniform positive charge background jellium distribution is

$$H\_0 = H\_{cl} + H\_b + H\_{cl-b\star} \tag{30}$$

where

$$H\_{EL} = \sum\_{i=1}^{N} \frac{p\_i^2}{2m} + \frac{1}{2} \varepsilon^2 \sum\_{\substack{j=1 \ j \neq i}}^{N} \sum\_{i=1}^{N} \frac{\exp\left(-\mu \left| \vec{r}\_i - \vec{r}\_j^\* \right| \right)}{\left| \vec{r}\_i - \vec{r}\_j^\* \right|},\tag{31}$$

represents the many-electron part,

$$H\_b = \frac{1}{2}e^2 \int \int \int d^3 \vec{x} \int \int \int d^3 \vec{x'} \frac{\rho\_{\vec{x}}^+ \rho\_{\vec{x'}}^+ \exp(-\mu|\vec{x} - \vec{x'}|)}{|\vec{x} - \vec{x'}|},\tag{32}$$

represents the jellium background and

$$H\_{cl-b} = -\varepsilon^2 \sum\_{i=1}^{N} \int \int \int d^3 \vec{x}' \frac{\rho\_{\vec{x}}^+ \exp\left(-\mu|\vec{x} - \vec{r}\_i|\right)}{|\vec{x} - \vec{r}\_i|},\tag{33}$$

represents the interaction between the electrons and the jellium background.The *N*-electron system in the jellium background potential constitutes an electrically neutral system, but the relative displacements of the positive and negative charges allow for plasma oscillations of the electron gas. A mathematical device using the coefficient *μ* in the exponential terms in Equations (31)–(33) is introduced to avoid some divergences; solutions are finally sought in the limit *μ* → 0. As a result of carrying out the integrals in Equations (32) and (33), the Hamiltonian (Equation (30)) turns out to be

$$H = H\_{cl} - \frac{1}{2}e^2 \frac{N^2}{V} \frac{4\pi}{\mu^2} \tag{34}$$

which manifestly diverges in the limit *μ* → 0. This is commonly referred to as the *μ*2- divergence.

It is most convenient to: (i) use the second quantized representations of the Hamiltonian for the bulk electron gas in a uniform positive background jellium potential (Equation (30)) using electron creation operator *c*†*kiσi* and the annihilation operator *<sup>c</sup>ki<sup>σ</sup>i* , *ki<sup>σ</sup>i* being the momentum (in units of *h*¯) and spin quantum numbers; (ii) employ the Fourier representation of the screened Coulomb interaction that appears in Equations (31)–(33); and finally (iii) seek the limits (*L*<sup>3</sup> = *V*) → ∞ (specifically in this order, with *L*−<sup>1</sup> << *μ*). Using the three steps described, after some tedious algebra, one finds that terms corresponding to momentum transfer*q* in the two-electron interactions term for which*q* =0 in the *Hel* part of the Hamiltonian cancels the abovementioned *μ*2-divergence, and along with the limits sought as per (iii), the Hamiltonian in the second quantized notation is

$$H\_0 = \sum\_{\vec{k}} \sum\_{\sigma} \frac{\hbar^2 \vec{k}^2}{2m} c\_{\vec{k}r}^{\dagger} c\_{\vec{k}r}^{\dagger}$$

$$= \frac{1}{2} \frac{e^2}{V} \sum\_{\vec{k}, \sigma\_1} \sum\_{\vec{p}, \sigma\_2} \sum\_{\vec{q} \neq \vec{0}} \left( \frac{4\pi}{q^2} c\_{\vec{k} + \vec{q}\sigma\_1}^{\dagger} c\_{\vec{p} - \vec{q}\sigma\_2}^{\dagger} c\_{\vec{p}\sigma\_2} c\_{\vec{k}r\_1} \right). \tag{35}$$

Scaling

$$=r\_s \tilde{k}\_\prime \tag{36a}$$

$$
\vec{\vec{p}} = r\_s \vec{p}\_\prime \tag{36b}
$$

$$\mathcal{V} = \frac{V}{r\_s^{\overline{3}}},\tag{36c}$$

and

$$
\vec{q} = r\_{\text{s}} \eta. \tag{36d}
$$

˜

allows us to write the Hamiltonian as

$$H\_0 = \left(\frac{e^2}{a\_0 r\_0^2}\right) \left[\sum\_{\hat{k}} \sum\_{\sigma} \frac{\vec{k}^2}{2} c\_{\hat{k}\sigma}^{\dagger} c\_{\hat{k}\sigma}\right.$$

$$+ \frac{1}{2} \frac{r\_0}{\vec{V}} \sum\_{\hat{k}} \sum\_{\hat{\mathcal{P}}} \sum\_{\hat{\mathcal{P}}} \sum\_{\sigma\_1} \sum\_{\sigma\_2} \left(\frac{4\pi}{\vec{q}^2} c\_{\hat{k}+\hat{q}\sigma\_1}^{\dagger} c\_{\hat{p}-\hat{q}\sigma\_2}^{\dagger} c\_{\hat{p}\sigma\_2} c\_{\hat{k}\sigma\_1}\right)\right].\tag{37}$$

where we have introduced a dimensionless variable *r*0 = *rs a*0 , *a*0 being the Bohr radius. In the high density limit *r*0 → <sup>∞</sup>, the second term in Equation (37) can be treated using first order perturbation theory even if the electron–electron interactions in the second term are quite strong. The result in the first order turns out to be essentially the same as in Equation (29), but higher order perturbation theory does not converge. Therefore, Bohm and Pines developed a non-perturbative approximation by carrying out canonical transformation of the Hamiltonian to represent pseudoparticles (elementary excitations of the many-electron gas) called *plasmons* which represent collective oscillations of the electron gas. The approximation involves linearization of the Hamiltonian concomitant with the neglect of certain terms whose phases are random and hence cancelable. Prior to discussing the canonical transformation of the Hamiltonian, we briefly visit their earlier semi-classical model which helps build insight in the *linearization* process and also in the *approximation* involved in the concomitant cancellation of terms having *random phases*.

˜

*k*

In the semi-classical model, both the electron gas and the positive charge in the bulk medium are considered to be uniformly spread over the entire volume with their collective centers coincident. An incidental movement of the electron density *ρ* (Figure 2) sets in oscillations of the electron gas described by the classical equation of motion:

$$m\frac{d^2\xi^\imath}{dt^2} = \left(\frac{1}{\varepsilon\_0}\varepsilon\rho\xi^\imath\right)(-\varepsilon),\tag{38}$$

wherein *eρξ*¯ denotes the surface charge density *σ*, the static average volume charge density being written as

$$\rho = \frac{N}{N\frac{4}{3}\pi r\_s^3} = \frac{3}{4\pi r\_s^3}.\tag{39}$$

The zero-point energy of the plasma oscillations is 12 *h*¯ *<sup>ω</sup>p* wherein the natural frequency of the plasma oscillations is

$$
\omega\_p = \sqrt{\frac{\rho \epsilon^2}{m \varepsilon\_0}} = \sqrt{\frac{3 \epsilon^2}{m r\_s^3}}.\tag{40}
$$

**Figure 2.** The 'jellium model' considers the collection of positive nuclear charges and negative electron charges smeared out uniformly across the region of a metal; the centers of positive and negative charges being coincident. A slight relative displacement of the centers of positive and negative parts of the total charge density sets in plasma oscillations.

Using the Fourier decomposition of the inter-electron interaction,

$$\frac{e^2}{r\_{i\bar{j}}} = \frac{1}{V} \sum\_{\vec{k}} c\_{\vec{k}} e^{i\vec{k} \cdot (\vec{r\_i} - \vec{r\_i})}.$$

the potential energy of the ith electron due to one electron charge *uniformly smeared* throughout the box is *e*2 *d*3 *rj*, i.e., *e*2 ∑ *kck* 0 1 *V d*3 *rje i k*·(*ri*−*ri*) 1 .

 Hence, the potential energy due to *all* the electrons is:

$$P(\vec{r\_i}) = \sum\_{j=1 \atop j \neq i}^{N} \frac{e^2}{r\_{ij}} = \frac{1}{V} \sum\_{j=1 \atop j \neq i}^{N} \sum\_{\vec{k}} c\_{\vec{k}} e^{i\vec{k} \cdot (\vec{r\_i} - \vec{r\_i})} {}\_{\prime} \tag{41a}$$

,

where *ck* = 4*πe*<sup>2</sup> *k*2 , except for*k* =0. The term corresponding to*k* =0 cancels the positive jellium; hence the potential energy of the *ith* electron due to all the electrons *and* the positive background is

$$\mathcal{U}(\vec{r\_i}) = \frac{1}{V} \sum\_{j=1 \atop j \neq i}^{N} \sum\_{\vec{k} \atop \vec{k} \neq \vec{0}} \frac{4\pi e^2}{\vec{k}^2} e^{i\vec{k} \cdot (\vec{r}\_i - \vec{r\_i})}.\tag{41b}$$

Now, in terms of the electron field operators, the total number of electrons is

$$N = \sum\_{\vec{\gamma}} \iiint d^3 \vec{r} \Psi^\dagger(q) \Psi(q) = \sum\_{\vec{\gamma}} \iiint d^3 \vec{r} \rho(q) = \iiint d^3 \vec{r} \rho(\vec{r}) = \sum\_{i=1}^N \iiint d^3 \vec{r} \delta(\vec{r} - \vec{r}\_i) \tag{42}$$

and the electron density is

$$\rho(\vec{r}) = \sum\_{i=1}^{N} \delta(\vec{r} - \vec{r}\_i) = \frac{1}{V} \sum\_{k=1}^{N} \rho\_{\vec{k}} e^{i\vec{k}\cdot\vec{r}},\tag{43}$$

wherein we have used the Fourier expansion of the charge density with the Fourier components being given by

$$\rho\_{\vec{k}} = \sum\_{i=1}^{N} e^{-i\vec{k}\cdot\vec{r\_i}}.\tag{44}$$

Identifying the force on the electron force as the negative gradient of the potential in Equation (41b), we arrive at the semi-classical equation of motion for the harmonic oscillator

$$m\ddot{\vec{r}}\_i = m\dot{\vec{v}}\_i = -\nabla\_i \mathcal{U}(\vec{r\_i})\_{\prime\prime}$$

which translates to the equation of motion for density fluctuations of the Fourier components in Equation (43):

$$\vec{\rho}\_{\vec{k}} = \left( -\sum\_{i=1}^{N} \left( \vec{k} \cdot \dot{\vec{r}}\_{i}^{\prime} \right)^{2} e^{-i\vec{k} \cdot \vec{r}\_{i}^{\prime}} - \frac{1}{V} \frac{4\pi Ne^{2}}{m} \rho\_{\vec{k}} - \frac{1}{V} \frac{4\pi e^{2}}{m} \sum\_{\substack{\vec{k}^{\prime} \neq \vec{k} \\ \vec{k}^{\prime} \neq \vec{0}}} \rho\_{\vec{k}^{\prime}} (\rho\_{\vec{k}} - \rho\_{\vec{k}^{\prime}}) \right. \tag{45}$$

The first term in Equation (45) is quadratic in *k*. It can be ignored if *k* · ˙*ri* <sup>2</sup>*average ω*<sup>2</sup> *p*; i.e., if one limits *k* to be small. The 'upper bound' on the wave number, denoted by *kc*, of the plasma oscillations is

$$k\_{\max} = k\_{\mathfrak{c}} \approx \frac{\omega\_p}{v\_f}.\tag{46}$$

Now, the integral of the charge density over the entire space adds up to the total number of electrons N, i.e.,

$$\int\iint d^3 \vec{r} \rho(\vec{r}) = \sum\_{i=1}^{N} \int\iint d^3 \vec{r} \delta(\vec{r} - \vec{r}\_i) = N,\tag{47a}$$

corresponding to

$$\rho(\vec{r}) = \sum\_{i=1}^{N} \delta(\vec{r} - \vec{r}\_i^\*). \tag{47b}$$

Using the fact that the Fourier expansion of the charge density

$$\rho(\vec{r}) = \frac{1}{V} \sum\_{\vec{k}=1}^{N} \rho\_{\vec{k}} e^{i\vec{k}\cdot\vec{r}},\tag{48}$$

$$\text{with } \rho\_{\vec{k}} = \sum\_{i=1}^{N} e^{-i\vec{k} \cdot \vec{r\_i}} \text{ and } \rho\_{\vec{k}}^{\*} = \sum\_{j=1}^{N} e^{+i\vec{k} \cdot \vec{r\_j}}.$$

In the third term on the right-hand side of Equation (45), we have *ρk* = *N* ∑ *i*=1 *e* −*i k*·*ri* and *N* 

*ρk*−*k* = ∑ *i*=1 *e i k* −*k*·*ri* , which involve oscillatory terms consisting of phase factors of modulus unity. It is like carrying out a sum of vectors in a complex plane whose directions are *random*, and one expects this to be a zero-sum addition. Thus: (i) neglecting the first term (enabled by placing an upper limit on *k*) and (ii) *linearizing* Equation (45) (i.e., *neglecting the quadratic terms*, taking advantage of the *random-phases*), we obtain

$$
\psi\_{\vec{k}} = -\frac{4\pi\mathfrak{p}c^2}{m}\rho\_{\vec{k}} = \omega\_p^2 \rho\_{k\prime} \tag{49}
$$

which essentially is an equation of motion for a simple harmonic oscillator. Quantized collective excitations of this system are elementary excitations. They are pseudo-particles called plasmons. The frequency of plasma oscillations is

$$
\omega\_p = \sqrt{\frac{4\pi\bar{\rho}e^2}{m}} = \sqrt{\frac{4\pi\left(\frac{3}{4\pi r\_s^3}\right)e^2}{m}} = \sqrt{\frac{3e^2}{mr\_s^3}}.\tag{50a}
$$

The zero-point energy of the plasma oscillations is 1 2*h*¯ *<sup>ω</sup>p* , where

$$
\hbar\omega\_p = \frac{2\sqrt{3}}{r\_s^{\frac{3}{2}}} \text{Ryd.}\tag{50b}
$$

In the Bohm–Pines method of canonical transformation of the Hamiltonian discussed below, the significance of the approximation involving *linearization* of the Hamiltonian concomitant with the neglect of terms having *random phases* gets *further* accentuated.

We have seen that the Hamiltonian for a bulk electron gas in a uniform positive background jellium potential is

$$H\_0 = \sum\_{i=1}^{N} \frac{p\_i^2}{2m} + \frac{1}{2} \frac{e^2}{V} \sum\_{j=1}^{N} \sum\_{i=1}^{N} \sum\_{\vec{k} \atop \vec{k} \neq \vec{0}} \frac{4\pi e^2}{k^2} e^{i\vec{k}\cdot(\vec{r}\_i - \vec{r}\_i)},\tag{51a}$$

and noting that the *j* = *i* term adds up to *N*, the total number of electrons, we arrive at

$$H\_0 = \sum\_{i=1}^{N} \frac{p\_i^2}{2m} + \frac{2\pi e^2}{V} \sum\_{\substack{\vec{k} \\ \vec{k}\neq \vec{0}}} \frac{1}{k^2} \sum\_{i=1}^{N} e^{i\vec{k}\cdot\vec{r\_i}} \sum\_{\substack{j=1 \\ j\neq i}}^{N} e^{-i\vec{k}\cdot\vec{r\_j}} = \sum\_{i=1}^{N} \frac{p\_i^2}{2m} + \frac{2\pi e^2}{V} \sum\_{\substack{\vec{k} \\ \vec{k}\neq \vec{0}}} \frac{1}{k^2} (\rho\_{\vec{k}}^\* \rho\_{\vec{k}} - N), \tag{51b}$$

where the last form is obtained by adding and subtracting the term corresponding to the *j* = *i*.

The quantum problem to be solved for the above Hamiltonian is

$$H\_0 \psi = E \psi. \tag{52}$$

Bohm and Pines recognized that the classical model which yielded plasma oscillations described by Equation (49) would be an approximation to a quantum model. One ought to seek a transformation of the Hamiltonian (Equations (51a) and (51b)) such that plasma oscillations appear explicitly as a set of *Hamiltonians for simple harmonic oscillators* for various *k* values limited by Equation (46). They therefore proposed canonical transformations (*q*, *p*) → (*Q*, *P*) of the Hamiltonian in Equations (51a) and (51b) to a *new* set of generalized coordinates Q and momenta P such that the new quantum Hamiltonian would have the *form*

$$H\_k = \frac{P\_k^\dagger P\_k}{2} + \frac{1}{2}\omega^2 Q\_k^\dagger Q\_{k\prime} \tag{53a}$$

which is characteristic of the Hamiltonian for a simple harmonic oscillator represented by the Hamiltonian (in units of *m* = 1)

$$
\hbar\_{SHO} = \frac{p^2}{2} + \frac{1}{2}\omega^2 q^2. \tag{53b}
$$

The transformation we seek is not inspired by actual measurements of the new coordinates and momenta; it is motivated only by seeking the form in Equation (53a). Hence, the operators *Q* & *P* need not necessarily be Hermitian. The Bohm–Pines strategy consists of starting with an *auxiliary* Hamiltonian

$$H\_1 = \sum\_{\substack{\vec{k} \\ \vec{k} < \vec{k}\_{\vec{k}}}} \frac{1}{2} P\_{\vec{k}}^{\dagger} P\_{\vec{k}} - M\_k P\_{\vec{k}}^{\dagger} \rho\_{\vec{k}'} \tag{54a}$$

with

$$M\_k = \sqrt{\frac{4\pi c^2}{Vk^2}}.\tag{54b}$$

We do not demand the operators *Q* & *P* to be Hermitian. Instead, by choosing

$$P\_{\vec{k}}^{\dagger} = P\_{-\vec{k}'} \tag{55a}$$

and

$$\mathcal{Q}\_{\vec{k}}^{\dagger} = \mathcal{Q}\_{-\vec{k}'} \tag{55b}$$

we see on recognizing that the summation over *k* and that over − *k* is equivalent considering the symmetry in the momentum space that, *H*1 is Hermitian, even if *Q* and *P* are not. The wavefunction depends only on the original set of electron coordinates; it cannot depend on any additional degrees of freedom. It is therefore judicious to employ subsidiary conditions

$$\frac{\partial \psi}{\partial Q\_{\vec{k}}} = 0,\tag{56a}$$

however, limited by *k* < *kc*. The derivative operator is the momentum,

$$P\_k = -i\hbar \frac{\partial}{\partial Q\_{\vec{k}}} \,, \tag{56b}$$

hence

> *Pkψ* = 0 for *k* < *kc*, (56c)

and

$$[Q\_{k'}P\_{k'}] = i\hbar \delta\_{k,k'} \tag{57}$$

which is just the uncertainty relation for canonically conjugate coordinates and momenta. Use of Equation (56c) in Equation (54a) ensures that

$$(H\_0 + H\_1)\psi = E\psi,\tag{58}$$

and the Hamiltonian (*<sup>H</sup>*0 + *<sup>H</sup>*1) describes the same quantum system. We seek a transformation affected by an operator

$$
\mathcal{U} = e^{\frac{\hbar}{\hbar}s},\tag{59a}
$$

with

$$S = \sum\_{\vec{k}\vec{\mu}\vec{k} < \vec{k}\_c} M\_k Q\_{\vec{k}} \rho\_{\vec{k}'} \tag{59b}$$

and

> *S*† = ∑ *k*; *k*<*k c MkQ*†*kρ*<sup>∗</sup>*k* = ∑ *k*; *k*<*k c MkQ*−*kρ*−*<sup>k</sup>* = *S*, (59c)

which gives

$$
\mathcal{U}^{\dagger} = \varepsilon^{-\dot{\mathbb{F}}s} = \mathcal{U}^{-1},\tag{59d}
$$

and we see that the transformation is unitary. It follows that

$$(P\_k)\_{new} = \mathcal{U}^{-1} \left( -i\hbar \frac{\partial \mathcal{U}}{\partial Q\_{\tilde{k}}} + \mathcal{U}P\_k \right) = P\_k - i\hbar \mathcal{U}^{-1} \frac{\partial \mathcal{U}}{\partial Q\_{\tilde{k}}} = P\_k + \mathcal{U}^{-1}[P\_k, \mathcal{U}]\_- = P\_k + M\_k \rho\_{\tilde{k}'} \tag{60a}$$

and the *ith* component of the operator is

$$(\vec{p}\_i)\_{new} = \vec{p}\_i - i \sum\_{\vec{k} \neq \vec{k} \, < \vec{k}\_c} M\_{\vec{k}} Q\_{\vec{k}} \vec{k} x^{i\vec{k} \cdot \vec{r}\_i}. \tag{60b}$$

Essentially, under the transformation under consideration, the operators *ri*, *Qk*, *ρk* remain invariant, but (*p*)*new* and (*Pk*)*new* are different. We now ask what the transformed Hamiltonian,

$$H\_{ncvw} = \mathcal{U}^{-1}(H\_0 + H\_1)\mathcal{U},\tag{61}$$

is. After some tedious algebra, it turns out to be

$$\mathfrak{H} = H\_{nuc} = \sum\_{i=1}^{N} \frac{p\_i^2}{2m} + \sum\_{\substack{\vec{k} \\ \vec{k} < \vec{k}\_c}} \frac{1}{2} \left( P\_{\vec{k}}^{\dagger} P\_{\vec{k}} + \omega\_p^2 Q\_{\vec{k}}^{\dagger} Q\_{\vec{k}} \right) - \sum\_{\substack{\vec{k} : \vec{k} \neq \vec{0}}} \frac{2\pi e^2}{Vk^2} N + H\_{s.r.} + H\_{int} + K\_{\prime} \tag{62}$$

where

$$H\_{s.r.} = \frac{1}{2} \sum\_{\vec{k}; \vec{k} \neq \vec{0}, \vec{k} > \vec{k}\_c} (\rho\_{\vec{k}}^{\*} \rho\_{\vec{k}} - N)\_{\prime} \tag{63}$$

$$H\_{\rm int} = -\frac{i}{2m} \sum\_{\vec{j}} \sum\_{\substack{\vec{k} \\ \vec{k} \neq \vec{0}}} M\_{\vec{k}} Q\_{\vec{k}} \vec{k} \cdot (2\vec{p}\_{\vec{j}} + \hbar \vec{k}) e^{-i\vec{k} \cdot \vec{r\_{j}}},\tag{64}$$

and

$$K = \frac{1}{2m} \sum\_{\substack{-\vec{k} \\ \vec{k} < \vec{k}\_c \\ \vec{k} < \vec{k}\_c \\ \vec{k} < \vec{k}\_c}}^{k \neq l} \sum\_{\substack{\vec{l} \\ \vec{k} < \vec{k}\_c \\ \vec{k} < \vec{k}\_c}} M\_{-\vec{k}} M\_{\vec{k}}(\vec{k} \cdot \vec{l}) \left\{ \sum\_j \left( Q\_{-\vec{k}} e^{+i\vec{k} \cdot \vec{r\_j}} \times Q\_{\vec{k}} e^{-i\vec{k} \cdot \vec{r\_j}} \right) \right\},\tag{65}$$

The new Hamiltonian has a manifestly complicated form. The term *K* (Equation (65)) is quadratic in the new coordinates and has *random phases* which would cancel out in a *linearization* process, as explained earlier in the context of the classical model and arrived at Equation (49). Linearization of the H = *Hnew* makes it possible to drop the operator K and justifies the term *random-phase approximation*. The rest of the Hamiltonian is

$$\mathfrak{H} = H\_{\text{new}} = \sum\_{i=1}^{N} \frac{p\_i^2}{2m} + \sum\_{\substack{\vec{k} \\ \vec{k} < \vec{k}\_{\text{c}}}} \frac{1}{2} \left( P\_{\vec{k}}^{\dagger} P\_{\vec{k}} + \omega\_p^2 Q\_{\vec{k}}^{\dagger} Q\_{\vec{k}} \right) - \frac{N}{V} \sum\_{\substack{\vec{k} : \vec{k} \neq \vec{0}}} \frac{2\pi e^2}{k^2} + H\_{\text{sf.}} + H\_{\text{int.}} \tag{66}$$

in which *Hs*.*r*. (Equation (63)) represents a set of *quasi-particles* interacting via *short-range* screened-Coulomb potential and given by

$$H\_{\rm s.r.} = \frac{1}{V} \sum\_{\vec{k}: \vec{k} \neq \vec{0}, \vec{k} > \vec{k}\_c} \frac{2\pi e^2}{k^2} (\rho\_{\vec{k}}^\* \rho\_{\vec{k}} - N)\_{\prime} \tag{67}$$

and

$$-\frac{N}{V} \sum\_{\vec{k}:\vec{k}\neq\vec{0}} \frac{2\pi e^2}{k^2} \tag{68}$$

is the self-energy of the electron gas.

*Hint* is accounted for by a further canonical transformation of the Hamiltonian (in which *K* is ignored) written in terms of transformed coordinates and momenta. Using the *random-phase approximation* concomitant with *linearization of the transformed Hamiltonian* (i.e., neglect of quadratic terms in the newer set of coordinates), the *Hint* term gets dropped, but in the process, the first two terms ge<sup>t</sup> somewhat modified, and the new approximate Hamiltonian becomes

$$\mathfrak{H} = H\_{\text{ncw}} = \{ \sum\_{\vec{k}, \vec{k} < \vec{k}\_c} \frac{1}{2} (P\_{\vec{k}}^{\dagger} P\_{\vec{k}} + \omega\_p^2 Q\_{\vec{k}}^{\dagger} Q\_{\vec{k}}) \} + \{ \sum\_{i=1}^N \frac{p\_i^2}{2m} (1 - \frac{\beta^2}{6}) + H\_{\text{s.r.}} \} - \{ \sum\_{\vec{k}, \vec{k} \ne \vec{0}} \frac{2\pi e^2}{V k^2} N \},\tag{69}$$

where

$$
\beta = \frac{k\_c}{k\_F} \,\tag{70}
$$

and

$$
\omega\_k^2 = \omega\_p^2 + \frac{2}{m} E\_F k^2 \tag{71}
$$

expresses a weak *k*-dependent dispersion of the plasma frequency. We have another subsidiary condition, similar to Equation (56c):

$$(P\_k + M\_k \rho\_{\tilde{k}}) \psi\_{new} = 0 \quad \text{for} \quad k < k\_c. \tag{72}$$

The kinetic energy part in the newer Hamiltonian is diminished by the factor (1 − *β*26 ); actual calculations show that *β* 0.7 and hence the kinetic energy part is reduced by about 8%. The long-range part of the interaction is what leads to the plasma oscillations corresponding to the first curly bracket in Equation (69). *Hs*.*r*. denotes the short-range screened-Coulomb interaction between the new pseudo-particles, which are elementary excitations called plasmons. The Hartree–Fock approximation accounted for only the static part of the density fluctuations of the collective behavior of an electron gas. The frozenorbital approximation that leads us to the Koopmans theorem highlights this approximation which limits it to the neglect of the Coulomb correlations. The method of canonical transformation of the Hamiltonian enables us address the Coulomb correlations albeit in an approximate manner by systematic and straightforward interpretation of Equation (69) in which the first curly bracket represents the collective oscillations of the electron gas resulting from the long-range part of the Coulomb interaction. A quantum of these oscillations is the plasmon. The second curly bracket represents the Hamiltonian for the screened Coulomb interaction, and the third represents the self-energy of the electron gas.

$$E\_{BP} = \frac{2.21}{r\_s^2} - \frac{0.916}{r\_s} + \frac{\sqrt{3}}{2r\_s^{3/2}} \beta^2 - \frac{0.916}{r\_s} \left(\frac{\beta^2}{2} - \frac{\beta^4}{48}\right).$$

The Bohm–Pines method elucidates the physical content of the *random-phase approximation (RPA) and the linearization process* it involves. There are other methods of arriving at the RPA, such as the Equation of Motion method [13] and the Greens function method [11]. The approximation is equivalent to summing over all the ring diagrams (along with the diagrams for the exchange interaction corresponding to each Coulomb term) in Feynman diagrammatic perturbation theory. Another equivalent approach to the RPA(E) results from the linearization of the Time-Dependent Hartree–Fock (TDHF) method developed by Dalgarno and Victor [14] and Amusia [15], and its relativistic version, namely the linearized Time Dependent Dirac–Hartree–Fock (TDDHF, often briefly denoted as TDDF) developed by Johnson and Lin [4]. In the next section, we summarize the linearized TDHF/TDDF approximations.

#### *2.3. Linearization of TDHF and that of TDDF Formalism*

The Hartree–Fock self-consistent field (HF-SCF) method accounts for correlations in many-electron dynamics that result by demanding that a many-fermion wavefunction must be anti-symmetric with respect to every exchange of pairs of the elementary particles. These correlations are therefore equivalently referred to as exchange correlations or as statistical (Fermi–Dirac) correlations. The Pauli Exclusion Principle automatically follows from it; hence, they are also sometimes called the Pauli correlations. The HF-SCF, however, only accounts for a static average of the density fluctuations of the many-electron system and thus leaves out what are known as *Coulomb correlations*. In the previous section, we discussed the RPA which employs a *linearization* technique and provides for a very successful methodology to account for the Coulomb correlations. We now proceed to discuss the RRPA [4], which is essentially based on *linearization* of the Time-Dependent Dirac–Hartree–Fock (TD-DHF, or just TDDF) family of coupled integro-differential equations. The linearized TDHF (RPAE) approximation [Amusia] was the precursor to the RRPA; it employs the *linearization* of the TDHF family of equations. The RPAE employs two-component spin-orbitals obtained as SCF solutions to the non-relativistic Schrödinger equation for the many-electron system. These spin-orbitals go into the construct of the Slater determinantal single-configuration wavefunctions. The RRPA is the relativistic extension of the RPAE. It employs four-component bi-spinor SCF solutions to the Dirac equation for the many-electron system which appear in the Slater determinant. The bi-spinors (spin-orbitals) are, however, admittedly time-dependent allowing for density fluctuations of the many-electron system. The resulting TD-DHF (TDHF) equations are non-linear; making a *linear approximation* to the TD-DHF (TDHF) equations result in the RRPA (RPAE).

The time-independent DHF equations for an *N*-electron closed-shell atomic system are

$$(h\_0 + V\_{DHF})\mu\_i = \varepsilon\_i \mu\_i \qquad \dots \text{i} = 1, 2, \dots, N\_\prime \tag{73}$$

where *ui* represent the four-component bispinor, *h*0 represents the Dirac Hamiltonian,

$$
\hbar\_0 = \vec{\pi} \cdot \vec{p} + \beta m - \frac{Ze^2}{r} \quad (\hbar = 1, c = 1), \tag{74}
$$

*εi* represent the DHF eigenvalue and *VDHF*(*r*) represents the DHF potential, given by

$$V\_{DHF}u(\vec{r}) = \sum\_{j=1}^{N} e^2 \int d^3r' \frac{\left\{ \left( u\_j^\dagger u\_j \right)' u - \left( u\_j^\dagger u \right)' u\_j \right\}}{|\vec{r} - \vec{r}'|},\tag{75}$$

and the prime denotes the argumen<sup>t</sup> over which integration is carried out. Solutions to the DHF equations are best represented by a Slater determinant

$$\psi^{(N)} = \frac{1}{\sqrt{N!}} \begin{vmatrix} u\_1(1) & \dots & \dots & \dots & u\_1(N) \\ u\_1(1) & \dots & \dots & \dots & u\_1(N) \\ \dots & \dots & \dots & \dots & \langle N|i \rangle \\ \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots \\ u\_N(1) & \dots & \dots & \dots & u\_N(N) \end{vmatrix} \tag{76a}$$

where

$$u\_{n\times n} = \frac{1}{r} \begin{pmatrix} iG\_{n\times}(r)\Omega\_{\text{km}}(\boldsymbol{\theta})\\ F\_{\text{m}\times}(r)\Omega\_{-\kappa m}(\boldsymbol{\theta}) \end{pmatrix} = \begin{pmatrix} \mu\_{+}\\ \mu\_{-} \end{pmatrix},\tag{76b}$$

with

$$\text{for}\quad j=\ell+\frac{1}{2},\quad \Omega\_{\text{km}}=\begin{pmatrix} \sqrt{\binom{j+m}{2j}} Y\_{\binom{\ell-j-\frac{1}{2}}{2}\binom{m}{\ell-m-\frac{1}{2}}}(^{\ell})\\ \sqrt{\binom{j-m}{2j}} Y\_{\binom{\ell-j-\frac{1}{2}}{2}\binom{m}{\ell-m+\frac{1}{2}}}(^{\ell}) \end{pmatrix}$$

$$\text{for}\quad j=\ell-\frac{1}{2},\quad\Omega\_{\text{km}}=\begin{pmatrix} -\sqrt{\left(\frac{j-m+1}{2j+2}\right)}Y\_{\left(\ell-j+\frac{1}{2}\right)\left(m\_{\ell'=m-\frac{1}{2}}\right)}(\hat{r})\\ \sqrt{\left(\frac{j+m+1}{2j+2}\right)}Y\_{\left(\ell-j+\frac{1}{2}\right)\left(m\_{\ell'=m+\frac{1}{2}}\right)}(\hat{r}) \end{pmatrix}.\tag{76c}$$

The electron densities of the DHF one-electron bispinors (spin-orbitals) represent only a time-average since the DHF model ignores electron correlations. Due to the electron correlations in the initial and the final state of a transition affected by what may be represented by an interaction operator

$$
\Omega = \nu\_{+} \boldsymbol{e}^{-i\omega t} + \nu\_{-} \boldsymbol{e}^{+i\omega t},
\tag{77}
$$

with the positive and negative frequency driving terms respectively denoted by

$$\nu\_{+} = \vec{\mathfrak{a}} \cdot \vec{A}; \quad \nu\_{-} = \nu\_{+}^{\dagger}. \tag{78}$$

*A* being the vector potential of the electromagnetic field, the DHF orbitals must be represented by time-dependent functions described by

$$u\_{i}(\vec{r}) \to u\_{i}(\vec{r}) + w\_{i-}(\vec{r})e^{\imath \omega \dagger} + w\_{i+}(\vec{r})e^{-\imath \omega \dagger} + \cdots \tag{79}$$

The dots ... at the end in Equation (79) represent higher harmonics. If we *rebuild* the (Time-Dependent) Dirac–Hartree–Fock scheme with all the higher harmonics, we obtain Non-Linear Time-Dependent Dirac–Hartree–Fock equations. Dalgarno and Victor [14] proposed the RPA linearization of the (non-relativistic TD-HF equations by dropping the higher harmonics. Following a similar logic, Johnson, Lin and Dalgarno [2–5] introduced the very same linearization in the TD-DHF system of coupled integro-differential equations for the orbitals *wi*±.The orbitals *wi*+ represent perturbation of the DHF orbitals due to the positive frequency part of the perturbation, and the orbitals *wi*− represent perturbation of the DHF orbitals due to the negative frequency part. The linearized Time-Dependent Dirac–Hartree–Fock (L-TslatdetD-DHF) equations are

$$(h\_0 + V\_{DHF} - \varepsilon\_i \mp \omega) w\_{i\pm} = (\nu\_{\pm} - V\_{\pm}^{(1)}) u\_i + \sum\_j \lambda\_{i j \pm} u\_{j \prime}; \quad i = 1, 2, \dots, N\_{\prime} \tag{80}$$

with

$$V\_{\pm}^{(1)}u\_i(\vec{r}) = \sum\_{j=1}^{N} \int d^3r' \frac{\left[\left(u\_j^{\dagger}w\_{j\pm}\right)'u\_i + \left(w\_{j\mp}^{\dagger}u\_j\right)'u\_i - \left(w\_{j\mp}^{\dagger}u\_i\right)'u\_j - \left(u\_j^{\dagger}u\_i\right)'u\_jw\_{j\pm}\right]}{|\vec{r}-\vec{r}'|},\tag{81}$$

which includes the Coulomb correlations that are omitted in the DHF method. The factors *<sup>λ</sup>ij*± in Equation (80) are the Lagrange's variational multipliers, introduced in the algebraic equations to ensure orthogonality of the perturbed orbitals *wj*∓ with respect to the unperturbed ones *uj*. Omission of the driving terms *ν*± gives us the fundamental RRPA equations:

$$\pm \left( h\_0 + V - \varepsilon\_i \right) w\_{i\pm} + V\_{\pm}^{(1)} u\_i \mp \sum\_j \lambda\_{i j \pm} u\_j = \omega w\_{i\pm}; \quad i = 1, 2, \dots, N\_\prime \tag{82}$$

The eigenvalues of the RRPA equations provide the linear approximation to the excitation spectra in both the discrete and the continuum. The positive and negative components *wi*± of the eigenfunctions describe the correlations which are omitted in the DHF formalism, respectively, in the excited state (both discrete and continuum) and in the initial states. The amplitude of transition from an initial state to the excited state described by the RRPA function *wi*±(*r*) corresponding to the frequency *ω* brought about by the interaction (Equation (76)) is

$$T = \sum\_{j=1}^{N} \int d^3r \left( w\_{i+}^{\dagger} v\_+ u\_i + w\_{i-}^{\dagger} v\_- u\_i \right),\tag{83a}$$

i.e.,

$$T = \sum\_{j=1}^{N} \int d^3r \left( w\_{i+}^{\dagger} \vec{\pi} \cdot \vec{A} u\_i + w\_{i-}^{\dagger} \vec{\pi} \cdot \vec{A} u\_i \right). \tag{83b}$$

It is a very important property of the (R)RPA equations that the transition matrix element is invariant under gauge transformations of the electromagnetic potentials. In actual calculations, one often employs truncated RRPA in which only the most important interchannel coupling is used. This leads to a slight disagreement in the estimation of the transition matrix elements in the length gauge and in the velocity gauge.

It is very convenient to have a pictorial representation of the correlations that are addressed in a many-electron theory. The diagrammatic representation developed by Feynman, first presented in the Spring of 1948 at the Pocono Conference, is suitable for our purpose. In AMO sciences, we (normally) do not work with positrons, but there are 'hole' states which are vacant states normally occupied by electrons. Thus we represent evolution of atomic states by vertical solid lines with reference to a time-axis going from the bottom to the top (left to right is an alternative convention). The atomic state lines are sometimes referred to as the 'trunk' of the diagram. A vertex in the diagram represents an intersection of a photon wavy line and the trunk. Particle lines point upwards and hole lines point downwards. Summing over only the ring graphs, as shown by Gell-Mann and Brueckner [21] has precisely the same effect as the linearization approximation that results in the random-phase approximation introduced by Bohm and Pines. Electron correlations are interpreted by recognizing that electrons exchange virtual photons which mediate the interaction between the electrons. The electromagnetic interaction is treated at the level of quantum theory. A positron is anelectron propagating backward in time. Figure 3 shows some of the lowest order diagrams [22] which contribute to the RRPA matrix elements.

**Figure 3.** Lowest order Feynman diagrams which contribute to the RRPA transition amplitude for the transition *a* → *a*¯. Time axis is from the bottom to the top of the page. The dashed lines represent electron–electron correlations, and the wiggly lines correspond to the photon operator. Arrows pointing upward (downward) are the electron (hole): (**a**) represents uncorrelated transition matrix element; (**b**,**<sup>c</sup>**) represent, respectively, the first order time-forward (i.e., final state) Coulomb and exchange terms; and (**d**,**<sup>e</sup>**) represent, respectively, the first order time-backward (i.e., initial state) Coulomb and exchange terms. (**f**–**h**) represent higher order ring diagrams.

Time-forward diagrams represent correlations in the final state in which configurationinteraction in the continuum is taken care of. In the RRPA, this corresponds to interchannel coupling. The time-backward diagrams represent correlations in the initial state. Variants of the RRPA that offer some advantages include the R-MCTD (Relativistic Multi-Configuration Tamm–Dancoff) method [23] and the RRPA-with-relaxation (RRPA-R) method [24]. The particle and the hole creation and annihilation operators that are referenced in the ring diagrams (Figure 3) are defined with respect to a vacuum consisting of a closed-shell Hartree–Fock (or Dirac–Hartree–Fock) fermion system, such as that described by Equation (8). *Two-particle two-hole* correlations are included, accommodating both creation and annihilation of a particle-hole pair. The focus of this article is to discuss the linearization approximation involved in the RPA and illustrate a few applications of the RRPA; hence, we omit an elaboration of the RRPA-R and R-MCTD techniques.
