Photoionization of Atomic Systems Using the Random-Phase Approximation Including Relativistic Interactions

Abstract

Approximation methods are unavoidable in solving a many-electron problem. One of the most successful approximations is the random-phase approximation (RPA). Miron Amusia showed that it can be used successfully to describe atomic photoionization processes of many-electron atomic systems. In this article, the historical reasons behind the term “random-phase approximation” are revisited. A brief introduction to the relativistic RPA (RRPA) developed by Walter Johnson and colleagues is provided and some of its illustrative applications are presented.

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 relativistic-random-phase approximation (RRPA) which is based on the Dirac equation rather than the Schrödinger equation [2,3,4,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

$$\begin{array}{c}\hfill H^{\left(N\right)}\left(q_1,q_2,\dots ,q_N\right)=\left\{\sum _{i=1}^N\left(-\frac{{\hslash }^2}{2m}{\nabla }_i^2-\frac{Z{e}^2}{r_i}\right)\right\}\\ \hfill +\left\{\sum _{i<j=1}^N\frac{e^2}{r_{ij}}\right\}=\left\{\sum _{i=1}^{N}f\left({\overrightarrow{r}}_i\right)\right\}+\left\{\frac{1}{2}\sum _{i=1;i\ne j}^N\sum _{j=1}^N\frac{e^2}{r_{ij}}\right\},\end{array}$$

(1)

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,7,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,17,18,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,3,4,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

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\Psi \left(t\right)=[\sum _i\sum _j{c}_i^{†}\langle i\left|f\right|j\rangle c_j\\ \hfill +\frac{1}{2}\sum _i\sum _j\sum _k\sum _l{c}_i^{†}{c}_j^{†}\langle ij\left|v\right|kl\rangle c_l{c}_k]|\Psi (t)\rangle ,\end{array}$$

(2a)

where

$$\langle ij\left|v\right|kl\rangle =\int d{q}_1\int d{q}_2{\psi }_i^{\ast }\left(q_1\right){\psi }_j^{\ast }\left(q_2\right)v(q_1,q_2){\psi }_k\left(q_1\right){\psi }_l\left(q_2\right),$$

(2b)

with the subscripts i, j, k, l denoting the set of four one-electron quantum numbers and the arguments ($q_r$) denoting the four coordinates of the $r^{th}$ 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\left(q\right)={\psi }_i\left(\overrightarrow{r}\right){\chi }_i\left(\zeta \right),$$

(2c)

wherein the spin part is either

$$\begin{array}{c}\hfill {\chi }_i\left(\zeta \right)=\alpha =\left[\begin{array}{c}1\\ 0\end{array}\right],\end{array}$$

(2d)

for $m_{s_i}=+\frac{1}{2}$ i.e., ↑, or

$$\begin{array}{c}\hfill {\chi }_i\left(\zeta \right)=\beta =\left[\begin{array}{c}0\\ 1\end{array}\right],\end{array}$$

(2e)

for $m_{s_i}=-\frac{1}{2}$ i.e., ↓.

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

$$\begin{array}{c}\hfill H=[\sum _{i,\alpha }\sum _{j,\beta }{c}_{i\alpha }^{†}\left(\int {\psi }_{i\alpha }^{\ast }\left(q\right)f\left(q\right){\psi }_{j\beta }^{\ast }\left(q\right)dq\right)c_{j\beta }\\ \hfill +\frac{1}{2}\sum _{i\alpha }\sum _{j\beta }\sum _{k\gamma }\sum _{l\delta }{c}_{i\alpha }^{†}{c}_{j\beta }^{†}\\ \hfill \times \left(\int dqd{q}^{\prime }{\psi }_{i\alpha }^{\ast }\left(q\right){\psi }_{j\beta }^{\ast }\left(q\right)v(q,q^{\prime }){\psi }_{l\delta }\left(q\right){\psi }_{k\gamma }\left(q\right)\right)c_{k\gamma }{c}_{l\delta }],\end{array}$$

(3a)

or more compactly using the Dirac notation for the integrals as

$$\begin{array}{c}\hfill H=\sum _{i,\alpha }\sum _{j,\beta }{c}_{i\alpha }^{†}\langle i\alpha \left|f\right|j\beta \rangle c_{j\beta }+\frac{1}{2}\sum _{i\alpha }\sum _{j\beta }\sum _{k\gamma }\sum _{l\delta }{c}_{i\alpha }^{†}{c}_{j\beta }^{†}\\ \hfill \times \langle i\alpha ,j\beta \left|v\right|l\delta ,k\gamma \rangle c_{k\gamma }{c}_{l\delta }.\end{array}$$

(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\ne j}^N\sum _{j=1}^{N}v(\overrightarrow{r_i},\overrightarrow{r_j})\right\}$$

in Equation (1) in a one-electron operator

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

by writing

$$\begin{array}{c}\hfill H^{\left(N\right)}\left(q_1,q_2,\dots ,q_N\right)=\sum _{i=1}^{N}f\left({\overrightarrow{r}}_i\right)+\sum _{i=1}^{N}F\left({\overrightarrow{r}}_i\right)\\ \hfill +\frac{1}{2}\sum _{i=1;i\ne j}^N\sum _{j=1}^{N}v(\overrightarrow{r_i},\overrightarrow{r_j})-\sum _{i=1}^{N}F\left({\overrightarrow{r}}_i\right),\end{array}$$

(4a)

or briefly as

$$\begin{array}{c}\hfill H^{\left(N\right)}\left(q_1,q_2,\dots ,q_N\right)=\\ \hfill (f+F)+(H_2-F)=O_1+O^{\prime },\end{array}$$

(4b)

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

$$\begin{array}{c}\hfill \left\{\frac{1}{2}\sum _{i=1;i\ne j}^N\sum _{j=1}^{N}v(\overrightarrow{r_i},\overrightarrow{r_j})-\sum _{i=1}^{N}F\left({\overrightarrow{r}}_i\right)\right\}\end{array}$$

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

$$\begin{array}{c}\hfill E^{\left(N\right)}=\langle {\Psi }^{\left(N\right)}|H^{\left(N\right)}|{\Psi }^{\left(N\right)}\rangle ,\end{array}$$

(5)

is minimized.

When $O^{\prime }$ is neglected, the unperturbed ground state wavefunction of the N-electron system is expressible as a determinant:

$$\begin{array}{c}\hfill {\mathsf{\Phi }}^{\left(N\right)}=\frac{1}{\sqrt{N!}}{\left[\begin{array}{cccc}{\psi }_{1\uparrow }\left(1\right)& \dots & \dots & {\psi }_{1\uparrow }\left(N\right)\\ {\psi }_{1\downarrow }\left(1\right)& \dots & \dots & {\psi }_{1\downarrow }\left(N\right)\\ \dots & \dots & \dots & \dots \\ {\psi }_{\frac{N}{2}\uparrow }\left(1\right)& \dots & \dots & {\psi }_{\frac{N}{2}\uparrow }\left(N\right)\\ {\psi }_{\frac{N}{2}\downarrow }\left(1\right)& \dots & \dots & {\psi }_{\frac{N}{2}\downarrow }\left(N\right)\end{array}\right]}_{N\times N}=\\ \hfill \frac{1}{\sqrt{N!}}{\left[\begin{array}{cccc}{\psi }_1\left(1\right)& \dots & \dots & {\psi }_1\left(N\right)\\ {\psi }_1\left(1\right)& \dots & \dots & {\psi }_1\left(N\right)\\ \dots & \dots & \dots & \dots \\ {\psi }_{N-1}\left(1\right)& \dots & \dots & {\psi }_{N-1}\left(N\right)\\ {\psi }_N\left(1\right)& \dots & \dots & {\psi }_N\left(N\right)\end{array}\right]}_{N\times N},\end{array}$$

(6)

wherein

$$\begin{array}{c}\hfill \left[f\left(\overrightarrow{r}\right)+F\left(\overrightarrow{r}\right)\right]{\psi }_{i\sigma }\left(\overrightarrow{r}\right)={\epsilon }_i{\psi }_{i\sigma }\left(\overrightarrow{r}\right),\end{array}$$

(7)

with ${\psi }_{i\sigma }\left(\overrightarrow{r}\right)={\psi }_{i\downarrow }\left(\overrightarrow{r}\right)$ or ${\psi }_{i\uparrow }\left(\overrightarrow{r}\right)$. Observe that the one-particle eigenvalue ${\epsilon }_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.

$$\begin{array}{c}\hfill {\mathsf{\Phi }}^{\left(N\right)}=\frac{1}{\sqrt{N!}}{\left[\begin{array}{ccccc}{\psi }_1\left(1\right)& \dots & \dots & \dots & {\psi }_1\left(N\right)\\ {\psi }_1\left(1\right)& \dots & \dots & \dots & {\psi }_1\left(N\right)\\ \dots & \dots & \dots & \langle j|i\rangle ={\psi }_i\left(q_j\right)& \dots \\ {\psi }_{N-1}\left(1\right)& \dots & \dots & \dots & {\psi }_{N-1}\left(N\right)\\ {\psi }_N\left(1\right)& \dots & \dots & \dots & {\psi }_N\left(N\right)\end{array}\right]}_{N\times N}.\end{array}$$

(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 ${\epsilon }_i>{\epsilon }_{N/2}$. In an ordered set of spin-orbitals, let us denote $p\le 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\ne p$. The choice of the operator $F\left(\overrightarrow{r}\right)$ which makes the energy functional (Equation (5)) a minimal is the one for which the matrix element

$$\begin{array}{c}\hfill \langle q\left|F\right|p\rangle =\sum _{i=1}^N\left[\langle iq\left|v\right|ip\rangle -\langle qi\left|v\right|ip\rangle \right],\end{array}$$

(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 $p^{th}$ spin-orbital replaced by the excited $q^{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{array}{c}\hfill [\left(-\frac{{\hslash }^2}{2m}{\nabla }^2-\frac{Z{e}^2}{r}\right){\psi }_p\left(\xi \right)\\ \hfill +\sum _{i=1}^N\left[\int d^4{V}^{\prime }\frac{{\psi }_i^{\ast }\left({\xi }^{\prime }\right){\psi }_i\left({\xi }^{\prime }\right){\psi }_p\left(\xi \right)e^2}{\overrightarrow{r}-\overrightarrow{r^{\prime }}}\right]\\ \hfill -\sum _{i=1}^N\delta (m_{s_p},m_{s_i})\left[\int d^4{V}^{\prime }\frac{{\psi }_i^{\ast }\left({\xi }^{\prime }\right)\left({\psi }_p\left({\xi }^{\prime }\right){\psi }_i\left(\xi \right)\right)e^2}{\overrightarrow{r}-\overrightarrow{r^{\prime }}}\right]]=\\ \hfill {\epsilon }_p{\psi }_p\left(\xi \right).\end{array}$$

(10a)

The four-dimensional integration $\int d^4{V}^{\prime }$ 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 non-ferromagnetic systems so that the number of spin-up and -down terms is equal. The energies ${\epsilon }_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{array}{c}\hfill [\left(-\frac{{\hslash }^2}{2m}{\nabla }^2-\frac{Z{e}^2}{r}\right){\psi }_p\left(\overrightarrow{r}\right)\\ \hfill +2\sum _{i=1}^{N/2}\left[\int d{V}^{\prime }{\left|{\psi }_i\left(\overrightarrow{r^{\prime }}\right)\right|}^{2}v(\overrightarrow{r},\overrightarrow{r^{\prime }})\right]{\psi }_p\left(\overrightarrow{r}\right)\\ \hfill -\sum _{i=1}^{N/2}{\psi }_i\left(\overrightarrow{r}\right)\left[\int d{V}^{\prime }{\psi }_i^{\ast }\left(\overrightarrow{r^{\prime }}\right){\psi }_p\left(\overrightarrow{r^{\prime }}\right)v(\overrightarrow{r},\overrightarrow{r^{\prime }})\right]]=\\ \hfill {\epsilon }_p{\psi }_p\left(\overrightarrow{r}\right),\end{array}$$

(10b)

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

$$\begin{array}{c}\hfill \frac{e^2}{\overrightarrow{r}-{\overrightarrow{r}}^{\prime }}=v(\overrightarrow{r},\overrightarrow{r^{\prime }}).\end{array}$$

(11)

Now,

$$\begin{array}{c}\hfill F{\psi }_p\left(\overrightarrow{r_2}\right)=\\ \hfill 2\sum _{i=1}^{N/2}\int d^3\overrightarrow{r_1}{\left|{\psi }_i\left(\overrightarrow{r_1}\right)\right|}^{2}v(\overrightarrow{r_1},\overrightarrow{r_2}){\psi }_p\left(\overrightarrow{r_2}\right)\\ \hfill -\sum _{i=1}^{N/2}\int d^3\overrightarrow{r_1}{\psi }_i^{\ast }\left(\overrightarrow{r_1}\right)v(\overrightarrow{r_1},\overrightarrow{r_2}){\psi }_p\left(\overrightarrow{r_1}\right){\psi }_i\left(\overrightarrow{r_2}\right),\end{array}$$

hence the HF SCF equation becomes

$$\begin{array}{c}\hfill [\left(-\frac{{\hslash }^2}{2m}{\nabla }^2-\frac{Z{e}^2}{r}\right){\psi }_p\left(\overrightarrow{r}\right)+F{\psi }_p\left(\overrightarrow{r}\right)=\\ \hfill (f+F){\psi }_p\left(\overrightarrow{r}\right)={\epsilon }_p{\psi }_p\left(\overrightarrow{r}\right).\end{array}$$

(12)

We write the momentum-dependent energies as $\epsilon \left(\overrightarrow{k}\right)$ or equivalently as $\epsilon \left(\overrightarrow{k}\right)$, since $\overrightarrow{p}=\hslash \overrightarrow{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.

Adding the jellium potential $V\left(\overrightarrow{r}\right)$ in the HF SCF equation we get

$$\begin{array}{c}\hfill [-\frac{{\hslash }^2}{2m}{\nabla }^2{\psi }_p\left(\overrightarrow{r}\right)+\cancel{V\left(\overrightarrow{r}\right){\psi }_p\left(\overrightarrow{r}\right)}\\ \hfill +\cancel{2\sum _{i=1}^{N/2}\left[\int d{V}^{\prime }{\left|{\psi }_i\left(\overrightarrow{r^{\prime }}\right)\right|}^{2}v(\overrightarrow{r},\overrightarrow{r^{\prime }})\right]{\psi }_p\left(\overrightarrow{r}\right)}\\ \hfill -\sum _{i=1}^{N/2}{\psi }_i\left(\overrightarrow{r}\right)\left[\int d{V}^{\prime }{\psi }_i^{\ast }\left(\overrightarrow{r^{\prime }}\right){\psi }_p\left(\overrightarrow{r^{\prime }}\right)v(\overrightarrow{r},\overrightarrow{r^{\prime }})\right]]=\\ \hfill {\epsilon }_p{\psi }_p\left(\overrightarrow{r}\right).\end{array}$$

(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{array}{c}\hfill \left[-\frac{{\hslash }^2}{2m}{\nabla }^2-V^{exchange}\left(q\right)\right]{\psi }_p\left(q\right)=\\ \hfill {\epsilon }_p{\psi }_p\left(\overrightarrow{r}\right),\end{array}$$

(14)

where

$$\begin{array}{c}\hfill V^{exchange}\left(q\right){\psi }_p\left(q\right)=\sum _{i=1}^N{V}_i^{exchange}\left(q\right){\psi }_p\left(q\right)=\\ \hfill \sum _{i=1}^N{\psi }_i\left(q\right)\left[\int d{q}^{\prime }\frac{{\psi }_i^{\ast }\left(q^{\prime }\right){\psi }_p\left(q^{\prime }\right)}{|\overrightarrow{r}-\overrightarrow{r^{\prime }}|}\right],\end{array}$$

(15a)

i.e.,

$$\begin{array}{c}\hfill V^{exchange}\left(q\right){\psi }_p\left(q\right)=\\ \hfill \sum _{i=1}^{N/2}{\psi }_i\left(\overrightarrow{r}\right)\left[{\psi }_i^{\ast }\left(\overrightarrow{r^{\prime }}\right){\psi }_p\left(\overrightarrow{r^{\prime }}\right)v(\overrightarrow{r},\overrightarrow{r^{\prime }})\right].\end{array}$$

(15b)

Using (i) box-normalized wavefunctions

$$\begin{array}{c}\hfill {\psi }_{\overrightarrow{k}\sigma }\left(\overrightarrow{r}\right)=L^{-3/2}{e}^{i\overrightarrow{k}·\overrightarrow{r}}{\chi }_{\sigma }\left(\zeta \right),\end{array}$$

(16)

and (ii)

$$\begin{array}{c}\hfill \varphi \left(\overrightarrow{r_1}\right)=\int d^3\overrightarrow{r_2}\frac{e^{i(\overrightarrow{k}-\overrightarrow{k^{\prime }})·\overrightarrow{r_2}}}{r_{12}}=\frac{4\pi e^{i(\overrightarrow{k}-\overrightarrow{k^{\prime }})·\overrightarrow{r_1}}}{|\overrightarrow{k}-\overrightarrow{k^{\prime }}{|}^2},\end{array}$$

(17)

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

$$\begin{array}{c}\hfill -\frac{{\hslash }^2}{2m}{\nabla }^2{\psi }_{\overrightarrow{k}}\left(\overrightarrow{r_1}\right)+{\epsilon }_{\overrightarrow{k}}{\psi }_{\overrightarrow{k}}\left(\overrightarrow{r_1}\right)=\epsilon \left(\overrightarrow{k}\right){\psi }_{\overrightarrow{k}}\left(\overrightarrow{r_1}\right),\end{array}$$

(18)

with the exchange term given by

$$\begin{array}{c}\hfill {\epsilon }_{\overrightarrow{k}}=\frac{4\pi e^2}{L^3}\sum _{\overrightarrow{k^{\prime }}}\frac{1}{|\overrightarrow{k}-\overrightarrow{k^{\prime }}{|}^2}.\end{array}$$

(19)

It follows that

$$\begin{array}{c}\hfill \epsilon \left(\overrightarrow{p}\right)=\frac{p^2}{2m}{\epsilon }_{exchange}\left(\overrightarrow{p}\right),\end{array}$$

(20)

where

$$\begin{array}{c}\hfill {\epsilon }_{exchange}\left(\overrightarrow{p}\right)=\frac{-e^2{p}_f}{\hslash \pi }\left[1+\frac{p_f^2-p^2}{2p_{f}p}ln\left|\frac{p+p_f}{p-p_f}\right|\right],\end{array}$$

(21)

$p_f$ 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

$$\begin{array}{c}\hfill E_{HF}^{atom}=\langle {\psi }^{\left(N\right)}\left|H\right|{\psi }^{\left(N\right)}\rangle =\\ \hfill \sum _{i=1}^N\langle i\left|f\right|i\rangle +\frac{1}{2}\sum _{j=1}^N\sum _{i=1}^N\left[\langle ij\left|g\right|ij\rangle -\langle ij\left|g\right|ji\rangle \right],\end{array}$$

(22)

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

$$\begin{array}{c}\hfill E_{HF}^{electrongasinjelliumpotential}=E_{KE}+E_{exchangecorrelation},\end{array}$$

(23)

where

$$\begin{array}{c}\hfill E_{KE}=2\frac{L^3}{\left(2\pi {\hslash }^2\right)}{\int }_{p=0}^{p=p_f}{p}^{2}dp{\int }_{\theta =0}^{\theta =\pi }sin\theta d\theta \\ \hfill \times {\int }_{\varphi =0}^{\varphi =2\pi }d\varphi \left[\frac{\overrightarrow{p}·\overrightarrow{p}}{2m}\right]=\frac{{\hslash }^2{L}^3}{10{\pi }^{2}m}{k}_f^5,\end{array}$$

(24)

with

$$\begin{array}{c}\hfill k_f={\left(\frac{3{\pi }^{2}N}{V}\right)}^{\frac{1}{3}},\end{array}$$

(25)

and

$$\begin{array}{c}\hfill E_{exchangecorrelation}=2\frac{L^3}{\left(2\pi {\hslash }^2\right)}\\ \hfill \times {\int }_{p=0}^{p=p_f}{p}^{2}dp{\int }_{\theta =0}^{\theta =\pi }sin\theta d\theta {\int }_{\varphi =0}^{\varphi =2\pi }d\varphi \left[\frac{1}{2}{\epsilon }_{exchange}\left(\overrightarrow{p}\right)\right]\\ \hfill =\frac{{\hslash }^2{L}^3}{10{\pi }^{2}m}{k}_f^5,\end{array}$$

(26a)

$$\begin{array}{c}\hfill E_{exchangecorrelation}=-\frac{V{e}^2}{4{\pi }^3}{\int }_{k=0}^{k=k_f}dk[2k_f{k}^2\\ \hfill +k(k_f^2-k^2)ln\left(\frac{k_f+k}{k_f-k}\right)].\end{array}$$

(26b)

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

$$\begin{array}{c}\hfill N\times \left(\frac{4}{3}\pi r_s^3\right)=V=\frac{3{\pi }^{2}N}{k_f^3},\end{array}$$

(27)

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

$$\begin{array}{c}\hfill \frac{E_{KE}}{N}=\frac{3{\hslash }^2}{10m}{\left(\frac{9\pi }{4}\right)}^{\frac{2}{3}}\frac{1}{r_s^2}=\frac{2.21}{r_s^2}Ryd,\end{array}$$

(28a)

and the exchange correlation energy per electron is

$$\begin{array}{c}\hfill \frac{E_{exchangecorrelation}}{N}=-\frac{0.916}{r_s}Ryd.\end{array}$$

(28b)

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

$$\begin{array}{c}\hfill \frac{E_{HF}}{N}=\left(\frac{2.21}{r_s^2}-\frac{0.916}{r_s}\right)Ryd.\end{array}$$

(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

$$\begin{array}{c}\hfill H_0=H_{el}+H_b+H_{el-b},\end{array}$$

(30)

where

$$\begin{array}{c}\hfill H_{EL}=\sum _{i=1}^N\frac{p_i^2}{2m}+\frac{1}{2}{e}^2\sum _{j=1;j\ne i}^N\sum _{i=1}^N\frac{exp(-\mu |\overrightarrow{r_i}-\overrightarrow{r_j}\left|\right)}{|\overrightarrow{r_i}-\overrightarrow{r_j}|},\end{array}$$

(31)

represents the many-electron part,

$$\begin{array}{c}\hfill H_b=\frac{1}{2}{e}^2\int \int \int d^3\overrightarrow{x}\int \int \int d^3\overrightarrow{x^{\prime }}\frac{{\rho }_{\overrightarrow{x}}^{+}{\rho }_{\overrightarrow{x^{\prime }}}^{+}exp(-\mu |\overrightarrow{x}-\overrightarrow{x^{\prime }}\left|\right)}{|\overrightarrow{x}-\overrightarrow{x^{\prime }}|},\end{array}$$

(32)

represents the jellium background and

$$\begin{array}{c}\hfill H_{el-b}=-e^2\sum _{i=1}^N\int \int \int d^3\overrightarrow{x}\frac{{\rho }_{\overrightarrow{x}}^{+}exp(-\mu |\overrightarrow{x}-\overrightarrow{r_i}\left|\right)}{|\overrightarrow{x}-\overrightarrow{r_i}|},\end{array}$$

(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 $\mu$ in the exponential terms in Equations (31)–(33) is introduced to avoid some divergences; solutions are finally sought in the limit $\mu \to 0$. As a result of carrying out the integrals in Equations (32) and (33), the Hamiltonian (Equation (30)) turns out to be

$$\begin{array}{c}\hfill H=H_{el}-\frac{1}{2}{e}^2\frac{N^2}{V}\frac{4\pi }{{\mu }^2},\end{array}$$

(34)

which manifestly diverges in the limit $\mu \to 0$. This is commonly referred to as the ${\mu }^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_{\overrightarrow{k_i}{\sigma }_i}^{†}$ and the annihilation operator $c_{\overrightarrow{k_i}{\sigma }_i}$, $\overrightarrow{k_i}{\sigma }_i$ being the momentum (in units of ℏ) 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^3=V)\to \infty$ (specifically in this order, with $L^{-1}<<\mu$). Using the three steps described, after some tedious algebra, one finds that terms corresponding to momentum transfer $\overrightarrow{q}$ in the two-electron interactions term for which $\overrightarrow{q}=\overrightarrow{0}$ in the $H_{el}$ part of the Hamiltonian cancels the abovementioned ${\mu }^2$-divergence, and along with the limits sought as per (iii), the Hamiltonian in the second quantized notation is

$$\begin{array}{c}\hfill H_0=\sum _{\overrightarrow{k}}\sum _{\sigma }\frac{{\hslash }^2\overrightarrow{k^2}}{2m}{c}_{\overrightarrow{k}\sigma }^{†}{c}_{\overrightarrow{k}\sigma }\\ \hfill +\frac{1}{2}\frac{e^2}{V}\sum _{\overrightarrow{k},{\sigma }_1}\sum _{\overrightarrow{p},{\sigma }_2}\sum _{\overrightarrow{q}\ne \overrightarrow{0}}\left(\frac{4\pi }{q^2}{c}_{\overrightarrow{k}+\overrightarrow{q}{\sigma }_1}^{†}{c}_{\overrightarrow{p}-\overrightarrow{q}{\sigma }_2}^{†}{c}_{\overrightarrow{p}{\sigma }_2}{c}_{\overrightarrow{k}{\sigma }_1}\right).\end{array}$$

(35)

Scaling

$$\begin{array}{c}\hfill \tilde{\overrightarrow{k}}=r_s\overrightarrow{k},\end{array}$$

(36a)

$$\begin{array}{c}\hfill \tilde{\overrightarrow{p}}=r_s\overrightarrow{p},\end{array}$$

(36b)

$$\begin{array}{c}\hfill \tilde{V}=\frac{V}{r_s^3},\end{array}$$

(36c)

and

$$\begin{array}{c}\hfill \tilde{q}=r_{s}q.\end{array}$$

(36d)

allows us to write the Hamiltonian as

$$\begin{array}{c}\hfill H_0=\left(\frac{e^2}{a_0{r}_0^2}\right)[\sum _{\tilde{\overrightarrow{k}}}\sum _{\sigma }\frac{\widetilde{\overrightarrow{k^2}}}{2}{c}_{\tilde{\overrightarrow{k}}\sigma }^{†}{c}_{\tilde{\overrightarrow{k}}\sigma }\\ \hfill +\frac{1}{2}\frac{r_0}{\tilde{V}}\sum _{\tilde{\overrightarrow{k}}}\sum _{\tilde{\overrightarrow{p}}}\sum _{\tilde{\overrightarrow{q}}\ne \overrightarrow{0}}\sum _{{\sigma }_1}\sum _{{\sigma }_2}\left(\frac{4\pi }{{\tilde{q}}^2}{c}_{\tilde{\overrightarrow{k}}+\tilde{\overrightarrow{q}}{\sigma }_1}^{†}{c}_{\tilde{\overrightarrow{p}}-\tilde{\overrightarrow{q}}{\sigma }_2}^{†}{c}_{\tilde{\overrightarrow{p}}{\sigma }_2}{c}_{\tilde{\overrightarrow{k}}{\sigma }_1}\right)],\end{array}$$

(37)

where we have introduced a dimensionless variable $r_0=\frac{r_s}{a_0}$, $a_0$ being the Bohr radius. In the high density limit $r_0\to \infty$, 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.

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 $\rho$ (Figure 2) sets in oscillations of the electron gas described by the classical equation of motion:

$$\begin{array}{c}\hfill m\frac{d^2\xi }{d{t}^2}=\left(\frac{1}{{\epsilon }_0}e\overline{\rho }\xi \right)(-e),\end{array}$$

(38)

wherein $e\overline{\rho }\xi$ denotes the surface charge density $\sigma$, the static average volume charge density being written as

$$\begin{array}{c}\hfill \overline{\rho }=\frac{N}{N\frac{4}{3}\pi r_s^3}=\frac{3}{4\pi r_s^3}.\end{array}$$

(39)

The zero-point energy of the plasma oscillations is $\frac{1}{2}\hslash {\omega }_p$ wherein the natural frequency of the plasma oscillations is

$$\begin{array}{c}\hfill {\omega }_p=\sqrt{\frac{\overline{\rho }{e}^2}{m{\epsilon }_0}}=\sqrt{\frac{3e^2}{m{r}_s^3}}.\end{array}$$

(40)

Using the Fourier decomposition of the inter-electron interaction,

$${\displaystyle \frac{e^2}{r_{ij}}=\frac{1}{V}\sum _{\overrightarrow{k}}{c}_{\overrightarrow{k}}{e}^{i\overrightarrow{k}·(\overrightarrow{r_i}-\overrightarrow{r_i})},}$$

the potential energy of the ith electron due to one electron charge uniformly smeared throughout the box is $e^2\int \int \int d^3{\overrightarrow{r}}_j$, i.e., ${\displaystyle e^2\sum _{\overrightarrow{k}}{c}_{\overrightarrow{k}}\left\{\frac{1}{V}\int \int \int d^3{\overrightarrow{r}}_j{e}^{i\overrightarrow{k}·(\overrightarrow{r_i}-\overrightarrow{r_i})}\right\}}$.

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

$${\displaystyle P\left(\overrightarrow{r_i}\right)=\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N\frac{e^2}{r_{ij}}=\frac{1}{V}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N\sum _{\overrightarrow{k}}{c}_{\overrightarrow{k}}{e}^{i\overrightarrow{k}·(\overrightarrow{r_i}-\overrightarrow{r_i})},}$$

(41a)

where $c_{\overrightarrow{k}}=\frac{4\pi e^2}{\left|{\overrightarrow{k}}^2\right|}$, except for $\overrightarrow{k}=\overrightarrow{0}$. The term corresponding to $\overrightarrow{k}=\overrightarrow{0}$ cancels the positive jellium; hence the potential energy of the $i^{th}$ electron due to all the electrons and the positive background is

$${\displaystyle U\left(\overrightarrow{r_i}\right)=\frac{1}{V}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N\sum _{\begin{array}{c}\overrightarrow{k}\\ \overrightarrow{k}\ne \overrightarrow{0}\end{array}}\frac{4\pi e^2}{k^2}{e}^{i\overrightarrow{k}·(\overrightarrow{r_i}-\overrightarrow{r_i})}.}$$

(41b)

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

$${\displaystyle N=\sum _{\zeta }\int \int \int d^3\overrightarrow{r}{\psi }^{†}\left(q\right)\psi \left(q\right)=\sum _{\zeta }\int \int \int d^3\overrightarrow{r}\rho \left(q\right)=\int \int \int d^3\overrightarrow{r}\rho \left(\overrightarrow{r}\right)=\sum _{i=1}^N\int \int \int d^3\overrightarrow{r}\delta (\overrightarrow{r}-\overrightarrow{r_i}),}$$

(42)

and the electron density is

$${\displaystyle \rho \left(\overrightarrow{r}\right)=\sum _{i=1}^N\delta (\overrightarrow{r}-\overrightarrow{r_i})=\frac{1}{V}\sum _{k=1}^N{\rho }_{\overrightarrow{k}}{e}^{i\overrightarrow{k}·\overrightarrow{r}},}$$

(43)

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

$${\displaystyle {\rho }_{\overrightarrow{k}}=\sum _{i=1}^N{e}^{-i\overrightarrow{k}·\overrightarrow{r_i}}.}$$

(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{\overrightarrow{r_i}}=m\dot{\overrightarrow{v_i}}=-\overrightarrow{{\nabla }_i}U\left(\overrightarrow{r_i}\right),$$

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

$${\displaystyle {\ddot{\rho }}_{\overrightarrow{k}}=-\sum _{i=1}^N{(\overrightarrow{k}·\dot{\overrightarrow{r_i}})}^2{e}^{-i\overrightarrow{k}·\overrightarrow{r_i}}-\frac{1}{V}\frac{4\pi N{e}^2}{m}{\rho }_{\overrightarrow{k}}-\frac{1}{V}\frac{4\pi e^2}{m}\sum _{\begin{array}{c}{\overrightarrow{k}}^{\prime }\ne \overrightarrow{k}\\ {\overrightarrow{k}}^{\prime }\ne \overrightarrow{0}\end{array}}{\rho }_{{\overrightarrow{k}}^{\prime }}({\rho }_{\overrightarrow{k}}-{\rho }_{{\overrightarrow{k}}^{\prime }}).}$$

(45)

The first term in Equation (45) is quadratic in k. It can be ignored if ${〈{\left(\overrightarrow{k}·\dot{\overrightarrow{r_i}}\right)}^2〉}_{average}\ll {\omega }_p^2$; i.e., if one limits k to be small. The ‘upper bound’ on the wave number, denoted by $k_c$, of the plasma oscillations is

$$k_{max}=k_c\approx \frac{{\omega }_p}{v_f}.$$

(46)

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

$$\begin{array}{c}\hfill {\displaystyle \int \int \int d^3\overrightarrow{r}\rho \left(\overrightarrow{r}\right)=\sum _{i=1}^N\int \int \int d^3\overrightarrow{r}\delta (\overrightarrow{r}-\overrightarrow{r_i})=N,}\end{array}$$

(47a)

corresponding to

$$\begin{array}{c}\hfill {\displaystyle \rho \left(\overrightarrow{r}\right)=\sum _{i=1}^N\delta (\overrightarrow{r}-\overrightarrow{r_i}).}\end{array}$$

(47b)

Using the fact that the Fourier expansion of the charge density

$${\displaystyle \rho \left(\overrightarrow{r}\right)=\frac{1}{V}\sum _{\overrightarrow{k}=1}^N{\rho }_{\overrightarrow{k}}{e}^{i\overrightarrow{k}·\overrightarrow{r}},}$$

(48)

with ${\displaystyle {\rho }_{\overrightarrow{k}}=\sum _{i=1}^N{e}^{-i\overrightarrow{k}·\overrightarrow{r_i}}}$ and ${\displaystyle {\rho }_{\overrightarrow{k}}^{\ast }=\sum _{j=1}^N{e}^{+i\overrightarrow{k}·\overrightarrow{r_j}}.}$

In the third term on the right-hand side of Equation (45), we have ${\displaystyle {\rho }_{\overrightarrow{k}}=\sum _{i=1}^N{e}^{-i\overrightarrow{k}·\overrightarrow{r_i}}}$ and ${\displaystyle {\rho }_{\overrightarrow{k}-\overrightarrow{k^{\prime }}}=\sum _{i=1}^N{e}^{i\overrightarrow{k^{\prime }}-\overrightarrow{k}·\overrightarrow{r_i}}}$, 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

$${\ddot{\rho }}_{\overrightarrow{k}}=-\frac{4\pi \overline{\rho }{e}^2}{m}{\rho }_{\overrightarrow{k}}={\omega }_p^2{\rho }_k,$$

(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

$$\begin{array}{c}\hfill {\omega }_p=\sqrt{\frac{4\pi \overline{\rho }{e}^2}{m}}=\sqrt{\frac{4\pi \left(\frac{3}{4\pi r_s^3}\right)e^2}{m}}=\sqrt{\frac{3e^2}{m{r}_s^3}}.\end{array}$$

(50a)

The zero-point energy of the plasma oscillations is $\frac{1}{2}\hslash {\omega }_p$, where

$$\begin{array}{c}\hfill \hslash {\omega }_p=\frac{2\sqrt{3}}{r_s^{\frac{3}{2}}}\mathrm{Ryd}.\end{array}$$

(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

$${\displaystyle H_0=\sum _{i=1}^N\frac{p_i^2}{2m}+\frac{1}{2}\frac{e^2}{V}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N\sum _{i=1}^N\sum _{\begin{array}{c}\overrightarrow{k}\\ \overrightarrow{k}\ne \overrightarrow{0}\end{array}}\frac{4\pi e^2}{k^2}{e}^{i\overrightarrow{k}·(\overrightarrow{r_i}-\overrightarrow{r_i})},}$$

(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 _{\begin{array}{c}\overrightarrow{k}\\ \overrightarrow{k}\ne \overrightarrow{0}\end{array}}\frac{1}{k^2}\sum _{i=1}^N{e}^{i\overrightarrow{k}·\overrightarrow{r_i}}\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{e}^{-i\overrightarrow{k}·\overrightarrow{r_j}}=\sum _{i=1}^N\frac{p_i^2}{2m}+\frac{2\pi e^2}{V}\sum _{\begin{array}{c}\overrightarrow{k}\\ \overrightarrow{k}\ne \overrightarrow{0}\end{array}}\frac{1}{k^2}\left({\rho }_{\overrightarrow{k}}^{\ast }{\rho }_{\overrightarrow{k}}-N\right),$$

(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 .$$

(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 $\overrightarrow{k}$ values limited by Equation (46). They therefore proposed canonical transformations $(q,p)\to (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

$$\begin{array}{c}\hfill H_k=\frac{P_k^{†}{P}_k}{2}+\frac{1}{2}{\omega }^2{Q}_k^{†}{Q}_k,\end{array}$$

(53a)

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

$$\begin{array}{c}\hfill h_{SHO}=\frac{p^2}{2}+\frac{1}{2}{\omega }^2{q}^2.\end{array}$$

(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

$$\begin{array}{c}\hfill H_1=\sum _{\begin{array}{c}\overrightarrow{k}\\ \overrightarrow{k}<\overrightarrow{k_c}\end{array}}\frac{1}{2}{P}_{\overrightarrow{k}}^{†}{P}_{\overrightarrow{k}}-M_k{P}_{\overrightarrow{k}}^{†}{\rho }_{\overrightarrow{k}},\end{array}$$

(54a)

with

$$\begin{array}{c}\hfill M_k=\sqrt{\frac{4\pi e^2}{V{k}^2}}.\end{array}$$

(54b)

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

$$\begin{array}{c}\hfill P_{\overrightarrow{k}}^{†}=P_{-\overrightarrow{k}},\end{array}$$

(55a)

and

$$\begin{array}{c}\hfill Q_{\overrightarrow{k}}^{†}=Q_{-\overrightarrow{k}},\end{array}$$

(55b)

we see on recognizing that the summation over $\overrightarrow{k}$ and that over $-\overrightarrow{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

$$\begin{array}{c}\hfill \frac{\partial \psi }{\partial Q_{\overrightarrow{k}}}=0,\end{array}$$

(56a)

however, limited by $k<k_c$. The derivative operator is the momentum,

$$\begin{array}{c}\hfill P_k=-i\hslash \frac{\partial }{\partial Q_{\overrightarrow{k}}},\end{array}$$

(56b)

hence

$$\begin{array}{c}\hfill P_k\psi =0\mathrm{for}k<k_c,\end{array}$$

(56c)

and

$$[Q_k,P_{k^{\prime }}]=i\hslash {\delta }_{k,k^{\prime }},$$

(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 ,$$

(58)

and the Hamiltonian $(H_0+H_1)$ describes the same quantum system. We seek a transformation affected by an operator

$$\begin{array}{c}\hfill U=e^{\frac{i}{\hslash }s},\end{array}$$

(59a)

with

$$\begin{array}{c}\hfill S=\sum _{\overrightarrow{k};\overrightarrow{k}<\overrightarrow{k_c}}{M}_k{Q}_{\overrightarrow{k}}{\rho }_{\overrightarrow{k}},\end{array}$$

(59b)

and

$$\begin{array}{c}\hfill S^{†}=\sum _{\overrightarrow{k};\overrightarrow{k}<\overrightarrow{k_c}}{M}_k{Q}_{\overrightarrow{k}}^{†}{\rho }_{\overrightarrow{k}}^{\ast }=\sum _{\overrightarrow{k};\overrightarrow{k}<\overrightarrow{k_c}}{M}_k{Q}_{-\overrightarrow{k}}{\rho }_{-\overrightarrow{k}}=S,\end{array}$$

(59c)

which gives

$$\begin{array}{c}\hfill U^{†}=e^{-\frac{i}{\hslash }s}=U^{-1},\end{array}$$

(59d)

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

$$\begin{array}{c}\hfill {\left(P_k\right)}_{new}=U^{-1}\left(-i\hslash \frac{\partial U}{\partial Q_{\overrightarrow{k}}}+U{P}_k\right)=P_k-i\hslash U^{-1}\frac{\partial U}{\partial Q_{\overrightarrow{k}}}=P_k+U^{-1}{[P_k,U]}_{-}=P_k+M_k{\rho }_{\overrightarrow{k}},\end{array}$$

(60a)

and the $i^{th}$ component of the operator is

$$\begin{array}{c}\hfill {\left(\overrightarrow{p_i}\right)}_{new}=\overrightarrow{p_i}-i\sum _{\overrightarrow{k};\overrightarrow{k}<\overrightarrow{k_c}}{M}_{\overrightarrow{k}}{Q}_{\overrightarrow{k}}\overrightarrow{k}{e}^{i\overrightarrow{k}·\overrightarrow{r_i}}.\end{array}$$

(60b)

Essentially, under the transformation under consideration, the operators $\overrightarrow{r_i},Q_{\overrightarrow{k}},{\rho }_{\overrightarrow{k}}$ remain invariant, but ${\left(\overrightarrow{p}\right)}_{new}$ and ${\left(P_k\right)}_{new}$ are different. We now ask what the transformed Hamiltonian,

$$H_{new}=U^{-1}(H_0+H_1)U,$$

(61)

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

$$\mathfrak{H}=H_{new}=\sum _{i=1}^N\frac{p_i^2}{2m}+\sum _{\begin{array}{c}\overrightarrow{k}\\ \overrightarrow{k}<\overrightarrow{k_c}\end{array}}\frac{1}{2}\left(P_{\overrightarrow{k}}^{†}{P}_{\overrightarrow{k}}+{\omega }_p^2{Q}_{\overrightarrow{k}}^{†}{Q}_{\overrightarrow{k}}\right)-\sum _{\overrightarrow{k};\overrightarrow{k}\ne \overrightarrow{0}}\frac{2\pi e^2}{V{k}^2}N+H_{s.r.}+H_{int}+K,$$

(62)

where

$$H_{s.r.}=\frac{1}{2}\sum _{\overrightarrow{k};\overrightarrow{k}\ne \overrightarrow{0},\overrightarrow{k}>\overrightarrow{k_c}}\left({\rho }_{\overrightarrow{k}}^{\ast }{\rho }_{\overrightarrow{k}}-N\right),$$

(63)

$$H_{int}=-\frac{i}{2m}\sum _j\sum _{\begin{array}{c}\overrightarrow{k}\\ \overrightarrow{k}\ne \overrightarrow{0}\end{array}}{M}_{\overrightarrow{k}}{Q}_{\overrightarrow{k}}\overrightarrow{k}·(2\overrightarrow{p_j}+\hslash \overrightarrow{k})e^{-i\overrightarrow{k}·\overrightarrow{r_j}},$$

(64)

and

$$K=\frac{1}{2m}\sum _{\begin{array}{c}-\overrightarrow{k}\\ \overrightarrow{k}<\overrightarrow{k_c}\end{array}}^{k\ne l}\sum _{\begin{array}{c}\overrightarrow{l}\\ \overrightarrow{l}<\overrightarrow{k_c}\end{array}}{M}_{-\overrightarrow{k}}{M}_{\overrightarrow{k}}(\overrightarrow{k}·\overrightarrow{l})\left\{\sum _j\left(Q_{-\overrightarrow{k}}{e}^{+i\overrightarrow{k}·\overrightarrow{r_j}}\times Q_{\overrightarrow{k}}{e}^{-i\overrightarrow{k}·\overrightarrow{r_j}}\right)\right\},$$

(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 $\mathfrak{H}=H_{new}$ makes it possible to drop the operator K and justifies the term random-phase approximation. The rest of the Hamiltonian is

$$\mathfrak{H}=H_{new}=\sum _{i=1}^N\frac{p_i^2}{2m}+\sum _{\begin{array}{c}\overrightarrow{k}\\ \overrightarrow{k}<\overrightarrow{k_c}\end{array}}\frac{1}{2}\left(P_{\overrightarrow{k}}^{†}{P}_{\overrightarrow{k}}+{\omega }_p^2{Q}_{\overrightarrow{k}}^{†}{Q}_{\overrightarrow{k}}\right)-\frac{N}{V}\sum _{\overrightarrow{k};\overrightarrow{k}\ne \overrightarrow{0}}\frac{2\pi e^2}{k^2}+H_{s.r.}+H_{int},$$

(66)

in which $H_{s.r.}$ (Equation (63)) represents a set of quasi-particles interacting via short-range screened-Coulomb potential and given by

$$H_{s.r.}=\frac{1}{V}\sum _{\overrightarrow{k};\overrightarrow{k}\ne \overrightarrow{0},\overrightarrow{k}>\overrightarrow{k_c}}\frac{2\pi e^2}{k^2}\left({\rho }_{\overrightarrow{k}}^{\ast }{\rho }_{\overrightarrow{k}}-N\right),$$

(67)

and

$$-\frac{N}{V}\sum _{\overrightarrow{k};\overrightarrow{k}\ne \overrightarrow{0}}\frac{2\pi e^2}{k^2}$$

(68)

is the self-energy of the electron gas.

$H_{int}$ 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 $H_{int}$ term gets dropped, but in the process, the first two terms get somewhat modified, and the new approximate Hamiltonian becomes

$$\mathfrak{H}=H_{new}=\{\sum _{\overrightarrow{k};\overrightarrow{k}<\overrightarrow{k_c}}\frac{1}{2}(P_{\overrightarrow{k}}^{†}{P}_{\overrightarrow{k}}+{\omega }_p^2{Q}_{\overrightarrow{k}}^{†}{Q}_{\overrightarrow{k}})\}+\{\sum _{i=1}^N\frac{p_i^2}{2m}(1-\frac{{\beta }^2}{6})+H_{s.r.}\}-\{\sum _{\overrightarrow{k};\overrightarrow{k}\ne \overrightarrow{0}}\frac{2\pi e^2}{V{k}^2}N\},$$

(69)

where

$$\beta =\frac{k_c}{k_F},$$

(70)

and

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

(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 }_{\overrightarrow{k}}){\psi }_{new}=0\mathbf{for}k<k_c.$$

(72)

The kinetic energy part in the newer Hamiltonian is diminished by the factor $(1-\frac{{\beta }^2}{6})$; actual calculations show that $\beta \simeq 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). $H_{s.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 frozen-orbital 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})u_i={\epsilon }_i{u}_i\dots i=1,2,\dots ,N,$$

(73)

where $u_i$ represent the four-component bispinor, $h_0$ represents the Dirac Hamiltonian,

$$h_0=\overrightarrow{\alpha }·\overrightarrow{p}+\beta m-\frac{Z{e}^2}{r}(\hslash =1,c=1),$$

(74)

${\epsilon }_i$ represent the DHF eigenvalue and $V_{DHF}\left(\overrightarrow{r}\right)$ represents the DHF potential, given by

$$V_{DHF}u\left(\overrightarrow{r}\right)=\sum _{j=1}^N{e}^2\int d^3{r}^{\prime }\frac{\left\{{\left(u_j^{†}{u}_j\right)}^{\prime }u-{\left(u_j^{†}u\right)}^{\prime }{u}_j\right\}}{|\overrightarrow{r}-\overrightarrow{r^{\prime }}|},$$

(75)

and the prime denotes the argument over which integration is carried out. Solutions to the DHF equations are best represented by a Slater determinant

$$\begin{array}{c}\hfill {\psi }^{\left(N\right)}=\frac{1}{\sqrt{N!}}\left|\begin{array}{ccccc}{u}_1\left(1\right)& \dots & \dots & \dots & u_1\left(N\right)\\ u_1\left(1\right)& \dots & \dots & \dots & u_1\left(N\right)\\ \dots & \dots & \dots & \dots & \langle N|i\rangle \\ \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots \\ u_N\left(1\right)& \dots & \dots & \dots & u_N\left(N\right)\end{array}\right|,\end{array}$$

(76a)

where

$$u_{n\kappa m}=\frac{1}{r}\left(\begin{array}{c}i{G}_{n\kappa }\left(r\right){\Omega }_{\kappa m}\left(\widehat{r}\right)\\ F_{n\kappa }\left(r\right){\Omega }_{-\kappa m}\left(\widehat{r}\right)\end{array}\right)=\left(\begin{array}{c}{u}_{+}\\ u_{-}\end{array}\right),$$

(76b)

with

$$\begin{array}{c}\hfill \mathrm{for}j=\ell +\frac{1}{2},{\Omega }_{\kappa m}=\left(\begin{array}{c}\sqrt{(}\frac{j+m}{2j})Y_{\left(\ell =j-\frac{1}{2}\right)\left(m_{{\ell }^{\prime }=m-\frac{1}{2}}\right)\left(\widehat{r}\right)}\\ \sqrt{(}\frac{j-m}{2j})Y_{\left(\ell =j-\frac{1}{2}\right)\left(m_{{\ell }^{\prime }=m+\frac{1}{2}}\right)\left(\widehat{r}\right)}\end{array}\right)\end{array}$$

$$\begin{array}{c}\hfill \mathrm{for}j=\ell -\frac{1}{2},{\Omega }_{\kappa m}=\left(\begin{array}{c}-\sqrt{(}\frac{j-m+1}{2j+2})Y_{\left(\ell =j+\frac{1}{2}\right)\left(m_{{\ell }^{\prime }=m-\frac{1}{2}}\right)\left(\widehat{r}\right)}\\ \sqrt{(}\frac{j+m+1}{2j+2})Y_{\left(\ell =j+\frac{1}{2}\right)\left(m_{{\ell }^{\prime }=m+\frac{1}{2}}\right)\left(\widehat{r}\right)}\end{array}\right).\end{array}$$

(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 }_{+}{e}^{-i\omega t}+{\nu }_{-}{e}^{+i\omega t},$$

(77)

with the positive and negative frequency driving terms respectively denoted by

$${\nu }_{+}=\overrightarrow{\alpha }·\overrightarrow{A};{\nu }_{-}={\nu }_{+}^{†}.$$

(78)

$\overrightarrow{A}$ being the vector potential of the electromagnetic field, the DHF orbitals must be represented by time-dependent functions described by

$$u_i\left(\overrightarrow{r}\right)\to u_i\left(\overrightarrow{r}\right)+w_{i-}\left(\overrightarrow{r}\right)e^{i\omega t}+w_{i+}\left(\overrightarrow{r}\right)e^{-i\omega t}+\cdots$$

(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,3,4,5] introduced the very same linearization in the TD-DHF system of coupled integro-differential equations for the orbitals $w_{i\pm }$.The orbitals $w_{i+}$ represent perturbation of the DHF orbitals due to the positive frequency part of the perturbation, and the orbitals $w_{i-}$ 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}-{\epsilon }_i\mp \omega )w_{i\pm }=\left({\nu }_{\pm }-V_{\pm }^{\left(1\right)}\right)u_i+\sum _j{\lambda }_{ij\pm }{u}_j;i=1,2,\dots ,N,$$

(80)

with

$$V_{\pm }^{\left(1\right)}{u}_i\left(\overrightarrow{r}\right)=\sum _{j=1}^N\int d^3{r}^{\prime }\frac{\left[{\left(u_j^{†}{w}_{j\pm }\right)}^{\prime }{u}_i+{\left(w_{j\mp }^{†}{u}_j\right)}^{\prime }{u}_i-{\left(w_{j\mp }^{†}{u}_i\right)}^{\prime }{u}_j-{\left(u_j^{†}{u}_i\right)}^{\prime }{u}_j{w}_{j\pm }\right]}{|\overrightarrow{r}-\overrightarrow{r^{\prime }}|},$$

(81)

which includes the Coulomb correlations that are omitted in the DHF method. The factors ${\lambda }_{ij\pm }$ in Equation (80) are the Lagrange’s variational multipliers, introduced in the algebraic equations to ensure orthogonality of the perturbed orbitals $w_{j\mp }$ with respect to the unperturbed ones $u_j$. Omission of the driving terms ${\nu }_{\pm }$ gives us the fundamental RRPA equations:

$$\pm (h_0+V-{\epsilon }_i)w_{i\pm }+V_{\pm }^{\left(1\right)}{u}_i\mp \sum _j{\lambda }_{ij\pm }{u}_j=\omega w_{i\pm };i=1,2,\dots ,N,$$

(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 $w_{i\pm }$ 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 $w_{i\pm }\left(\overrightarrow{r}\right)$ corresponding to the frequency $\omega$ brought about by the interaction (Equation (76)) is

$$\begin{array}{c}\hfill T=\sum _{j=1}^N\int d^{3}r\left(w_{i+}^{†}{v}_{+}{u}_i+w_{i-}^{†}{v}_{-}{u}_i\right),\end{array}$$

(83a)

i.e.,

$$\begin{array}{c}\hfill T=\sum _{j=1}^N\int d^{3}r\left(w_{i+}^{†}\overrightarrow{\alpha }·\overrightarrow{A}{u}_i+w_{i-}^{†}\overrightarrow{\alpha }·\overrightarrow{A}{u}_i\right).\end{array}$$

(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.

Time-forward diagrams represent correlations in the final state in which configuration-interaction 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.

3. Illustrative Examples

As a first example, consider the dipole photoelectron angular distribution resulting from the photoionization of the 5s subshell of the closed-shell Xe atom. The general form of the dipole angular distribution for linearly polarized incident radiation for subshell i is given as the differential photoionization cross section by [25]

$$\frac{d{\sigma }_i}{d\Omega }=\frac{{\sigma }_i}{4\pi }[1+{\beta }_i{P}_2(cos\theta )],$$

(84)

where ${\sigma }_i$ is the total subshell cross section, $P_2\left(x\right)=(3x^2-1)2$, $\theta$ is the angle between photon polarization and photoelectron momentum directions, and ${\beta }_i$ is the dipole angular distribution asymmetry parameter. Nonrelativistically, for a closed shell atom, ${\beta }_i$ = 2 and is independent of energy for an initial ns subshell [24]. From a physical point of view, this occurs because, nonrelativistically, there is only a single final state partial wave, characterized by ns $\to \epsilon$p, so that there is nothing to interfere with, and the angular distribution is just determined by the symmetry of the $\epsilon$p-wave. When relativistic interactions are included, the situation is changed in that there are then two possible photoionizing transitions ns $\to \epsilon$p${}_{1/2}$ and ns $\to \epsilon$p${}_{3/2}$, and these can interfere with one another, giving rise to both a deviation of ${\beta }_i$ from the value two, and an energy dependence. This behavior will be most evident near Cooper minima [26] owing to the fact that the two relativistic channels exhibit the minima at slightly different energies [27], so that in this region the matrix elements for the two relativistic channels can be vastly different, both in magnitude and phase.

The Xe 5s photoionization cross section and $\beta$ parameter calculated using RRPA, including correlation in the form of interchannel coupling with 5p and 4d relativistic photoionization channels, are depicted in Figure 4. The cross section exhibits a deep minimum, indicative of the Cooper minima in the 5s $\to \epsilon$p${}_{1/2}$ and 5s$\to \epsilon$p${}_{3/2}$ channels. Correspondingly, the $\beta$ parameter shows a deep minimum around that energy, and over a significant energy range, $\beta$ is energy-dependent. It is to be emphasized, that, while the nonrelativistic RPAE does pretty well on the cross section [1], the $\beta$ parameter predicted by the nonrelativistic calculation is constant and equal to two at all energies. Thus, it is evident that the addition of relativistic effects to the original nonrelativistic RPAE brings out additional physical effects. It should also be noted that, although the RRPA result for the cross section is rather good in the Cooper minimum region, the predictions for $\beta$ are not, owing to the mission of satellite photoionization channels in the RRPA [28].

The next example is the phenomenon of spin-orbit-activated interaction interchannel coupling (SOIAIC). Basically, this results from the interchannel coupling among the photoionization channels emanating from the two members of a spin-orbit doublet, nl, with $j=l\pm \frac{1}{2}$. The phenomenon was discovered experimentally in the photoionization of the 3d subshell of Xe [30] and subsequently explained theoretically [31]. The explanation given in [31] was based upon the nonrelativistic RPAE that was artificially made to include the spin-orbit splitting of the Xe 3d spin-orbit split levels. Subsequently, RRPA was applied and verified the explanation [32]. The results, both experimental and theoretical, are shown in Figure 5.

The key point here is that both the 3d${}_{5/2}$ and 3d${}_{3/2}$ cross sections exhibit shape resonances [33] just above their respective thresholds. As a result, the 3d${}_{3/2}$ cross section is significantly larger than the 3d${}_{5/2}$ cross section just above the 3d${}_{3/2}$ threshold. As a general rule, when a large photoionization cross section is degenerate with a small one, the small cross section is altered owing to interchannel coupling (configuration interaction in the continuum) [34]. This is exactly what is seen here in the form of the structure in the 3d${}_{5/2}$ cross section just above the 3d${}_{3/2}$ threshold. It is evident that the relativistic spin-orbit splitting of the 3d thresholds is crucial to the existence of this SOIAIC phenomenon, and it is also clear that SOIAIC is likely to be in evidence for the photoionization of inner subshells that have a significant spin-orbit splitting and exhibit a near-threshold shape resonance. Furthermore, as a result of the fact that the photoionization cross sections of inner shells of confined atoms often have significant near-threshold maxima [35], owing to the phenomenon of confinement resonances [36], the SOIAIC should be much more generally exhibited [37].

The asymptotic branching ratios of spin-orbit doublets have been a topic of interest since the 1960s. Earlier, it was expected that far above thresholds, the ratio of cross sections for the members of a spin-orbit doublet, nl, with $j=l\pm \frac{1}{2}$, respectively, should approach $(l+1)/l$, known as the statistical value, at asymptotically high energies [38]. However, it was later shown that relativistic effects on the initial state wave functions would cause the ratio to drop below the statistical [39,40]. Recently, experimental technology has improved to the point that this prediction has been verified experimentally [41]. In addition, it had been found that, in the vicinity of inner shell thresholds, there can be very large swings of the branching ratio over relatively small energy ranges [41,42]. However, all of this phenomenology cannot occur within a nonrelativistic framework; relativistic interactions are required.

As a particular example, the calculated branching ratio for Xe 5p using the RRPA with coupling among relativistic channels from all subshells except 1s [42] is depicted in Figure 6 over a very large energy range, from threshold to 500 a.u. (approximately 13.6 keV), where it is seen that, at the highest energies, that branching ratio is about 1.6, well below the statistical value for an initial np doublet, in keeping with the earlier predictions. In addition, there are seen to be significant excursions from monotone decreasing results in the vicinity of the n = 4, n = 3 and n = 2 thresholds. These excursions are the result of interchannel coupling between the relativistic photoionization channels resulting from the 5p${}_j$ initial states and the channels from the inner shells. The fact that there are significant excursions of the branching ratios in these energy ranges means that the actual interchannel coupling matrix elements are themselves dependent upon relativistic interactions.

Figure 6 also shows the results of RRPA calculation including only coupling among the five 5p relativistic photoionization channels. It is evident that the truncated RRPA results agree quite well with the fully coupled branching ratios away from the inner shell thresholds but do not reproduce the significant excursions from the smooth curve in the vicinity of the inner shell thresholds, thereby showing conclusively that these excursions are the result of interchannel coupling with the inner shell photoionization channels.

The Cooper minimum [26] mentioned above in connection with Xe 5s photoionization, is a ubiquitous phenomenon that pervades the photoionization of outer and/or near-outer subshells of all of the elements of the periodic table [43,44]. Among the interesting facets of the influence of relativistic interactions which cause a single nonrelativistic Cooper minimum to be split into several relativistic Cooper minima dependent upon the total angular momentum, j, is of the initial and final states of the relativistic photoionizing transition [45,46]. In addition, the locations of Cooper minima depend very sensitively on many-body correlations, in addition to relativistic interactions; as a matter of fact, the Xe 5s Cooper minimum is below threshold in the discrete region at the level of single particle calculations [47] but appear in the continuum in calculations that include significant many-body effects [5]. Thus, RRPA is an ideal formalism to study these relativistic effects in Cooper minima. For completeness, it should be pointed out here that in the neighborhood of the Xe 5s Cooper minimum, quadrupole effects become important in the angular distribution, although the total subshell cross section is virtually unaffected. The interference between dipole and quadrupole photoionizing transitions leads to extra terms in the expression for the differential cross section, but $\beta$ remains unaffected to the first order [29]. In addition, calculations of the effects of quadrupole transitions on the differential cross section (photoelectron angular distribution) using the RRPA methodology has been shown to be in good agreement with the experiment [48].

Recent work on the splitting of Cooper minima (CM) for heavy and superheavy atoms has illustrated the importance of including correlation in the calculations [49]. As an example, given in

Table 1. Positions of the Cooper minima (CM) in Dirac–Fock (DF), two-channel RRPA and RRPA with full coupling in photoelectron energy in (a.u.) [5].

Atoms	Photoelectron Energy (a.u.)			Photoelectron Energy (a.u.)
	of the 6s $\to \mathit{\epsilon }$p${}_{3/2}$ CM			of the 6s $\to \mathit{\epsilon }$p${}_{1/2}$ CM
	DF	RRPA	RRPA (Full	DF	RRPA	RRPA (Full
		(2 Channel)	Coupling)		(2 Channel)	Coupling)
Hg (Z = 80)	2.73	2.17	4.17	0.55	0.67	3.67
Rn (Z = 86)	3.63	1.43	5.93			4.43
Ra (Z = 88)	4.01	1.38	6.38			3.88
No (Z = 102)	10.91	5.70	11.70			6.70
Cn (Z = 112)	32.53	29.32	24.82			4.82
Og (Z = 118)	62.12	63.02	171.02			3.52

are the positions of the various 6s CM for six elements obtained at three levels of calculation: independent particle Dirac–Fock (DF), two-channel RRPA coupling on the two relativistic channels arising from 6s photoionization, and RRPA with full coupling of all the channels that might affect the result. To begin with, there are huge differences between the positions of the 6s $\to \epsilon$p${}_{3/2}$ CM and the 6s $\to \epsilon$p${}_{1/2}$ CM, and the differences increase with Z to an astounding degree in Og—a splitting of more than 4.5 keV. This comes about owing to the spin-orbit force which is attractive for the 6s $\to \epsilon$p${}_{1/2}$ final state but repulsive for the 6s $\to \epsilon$p${}_{3/2}$ final state. In addition, it is clear from Table 1 that correlation in the form of interchannel coupling induces rather significant changes in the location of 6s $\to \epsilon$p${}_{3/2}$ CM compared to the two-channel and DF results, changes that generally increase with Z, indicating the crucial importance of the interchannel coupling in the determination of the position of the CM in these heavy and superheavy elements.

4. Concluding Remarks

It is clear from the above exposition and examples that RRPA, which includes significant many-body correlation including initial state two-particle two-hole terms and final state interchannel coupling, can be suitably applied to a number of aspects of atomic photoionization and give physical insight into what makes the results for these processes what they are. Furthermore, it must be noted that only a few of the many examples that have been studied over the years are presented above. Thus, it is evident that the RRPA methodology has contributed greatly to our understanding of atomic photoionization as well as many other atomic processes [50]. However, it must be recalled that all of this was made possible by the pioneering work on the RPAE method by Miron Amusia and his collaborators.