# Pushing the Limits of the Periodic Table -A Review on Atomic Relativistic Electronic Structure Theory and Calculations for the Superheavy Elements *

Odile R. Smits, Paul Indelicato, Witold Nazarewicz, Morten Piibeleht, Peter Schwerdtfeger

<figure><img src="image_1.png" /></figure>

## To cite this version:

Odile R. Smits, Paul Indelicato, Witold Nazarewicz, Morten Piibeleht, Peter Schwerdtfeger. Pushing the Limits of the Periodic Table -A Review on Atomic Relativistic Electronic Structure Theory and Calculations for the Superheavy Elements *. Physics Reports, 2023, 1035, pp.1-57. 10.1016/j.physrep.2023.09.004. hal-03884148

HAL Id: hal-03884148

https://hal.science/hal-03884148

Submitted on 6 Dec 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Pushing the Limits of the Periodic Table – A Review on Atomic Relativistic Electronic Structure Theory and Calculations for the Superheavy Elements*

O. R. Smits<sup>a,1</sup>, P. Indelicato<sup>b,2</sup>, W.

Nazarewicz<sup>c,3</sup>, M. Piibeleht<sup>1</sup>, P.

Schwerdtfeger<sup>d,1</sup>

¹ Centre for Theoretical Chemistry and Physics, The New Zealand Institute for Advanced Study, Massey University Auckland, 0632 Auckland, New Zealand

²Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-PSL Research University, Collège de France, Case 74; 4, place Jussieu, F-75005 Paris, France

³Facility for Rare Isotope Beams and Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, USA

*In memoriam to two of the pioneers in this field, Jean-Paul Desclaux (Grenoble) and Sigurd Hofmann (Darmstadt)*

December 6, 2022

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>The Dirac Equation in Strong Coulomb Fields</td><td>4</td></tr><tr><td>2.1</td><td>The QED Lagrangian</td><td>4</td></tr><tr><td>2.2</td><td>The Many-Electron Dirac-Coulomb-Breit Hamiltonian</td><td>5</td></tr><tr><td>2.3</td><td>The one-particle Dirac equation</td><td>7</td></tr><tr><td>2.3.1</td><td>Point nucleus and self-adjointness</td><td>7</td></tr><tr><td>2.3.2</td><td>Nuclear Recoil and Uehling terms</td><td>8</td></tr><tr><td>2.3.3</td><td>Finite nuclear charge distributions</td><td>10</td></tr><tr><td>2.3.4</td><td>1s energy level reaching the negative energy continuum</td><td>14</td></tr><tr><td>2.4</td><td>Electron states in the super-critical region</td><td>15</td></tr><tr><td>2.4.1</td><td>Energy-projected Dirac equation</td><td>15</td></tr><tr><td>2.4.2</td><td>Hartree-Fock-Bogoliubov equation analogy</td><td>17</td></tr><tr><td>2.4.3</td><td>Perturbative approach</td><td>20</td></tr><tr><td>2.4.4</td><td>Analytical continuation</td><td>21</td></tr><tr><td>2.4.5</td><td>Gamow states</td><td>22</td></tr><tr><td>2.5</td><td>Positron production in the super-critical regime</td><td>25</td></tr><tr><td>2.6</td><td>Experimental perspective: heavy-ion collisions</td><td>27</td></tr><tr><td>3</td><td>Multi-Configuration Dirac-Hartree-Fock</td><td>29</td></tr><tr><td>4</td><td>Quantum Electrodynamic Effects</td><td>32</td></tr><tr><td>4.1</td><td>S-matrix formalism</td><td>34</td></tr><tr><td>4.2</td><td>Two-times Green's function method</td><td>38</td></tr><tr><td>4.3</td><td>Covariant evolution-operator procedure</td><td>38</td></tr><tr><td>4.4</td><td>Calculation of QED corrections</td><td>39</td></tr><tr><td>4.4.1</td><td>One-electron radiative corrections</td><td>39</td></tr><tr><td>4.4.2</td><td>Two-electron radiative and non-radiative corrections</td><td>43</td></tr></table>

<footnote><sup>a</sup>e-mail: smits.odile.rosette@gmail.com</footnote>

<footnote><sup>b</sup>e-mail: paul.indelicato@lkb.upmc.fr</footnote>

<footnote><sup>c</sup>e-mail: witek@frib.msu.edu</footnote>

<footnote><sup>d</sup>e-mail: peter.schwerdtfeger@gmail.com</footnote>
<table><tr><td>4.4.3</td><td>Effective QED Hamiltonians</td><td>44</td></tr><tr><td>5</td><td>Electron Correlation</td><td>47</td></tr><tr><td>6</td><td>Atomic Structure Calculations of the Superheavy Elements</td><td>50</td></tr><tr><td>6.1</td><td>Dominant Electron Configurations</td><td>50</td></tr><tr><td>6.2</td><td>Ionization potentials</td><td>51</td></tr><tr><td>6.3</td><td>QED effects</td><td>53</td></tr><tr><td>6.4</td><td>Atomic static dipole polarizabilities</td><td>55</td></tr><tr><td>6.5</td><td>Electron Localization Function</td><td>56</td></tr><tr><td>6.6</td><td>Examples of relativistic effects on the chemistry of SHE</td><td>57</td></tr><tr><td>7</td><td>General Considerations</td><td>58</td></tr><tr><td>7.1</td><td>Placing new elements on the periodic table</td><td>58</td></tr><tr><td>7.2</td><td>Dominant ground state configurations predicted by electronic structure calculations</td><td>60</td></tr><tr><td>7.3</td><td>Periodic Table - How far can we go?</td><td>61</td></tr><tr><td>8</td><td>Conclusions</td><td>61</td></tr><tr><td>9</td><td>Appendix A: The Self-Adjointness of the Dirac-Coulomb Hamiltonian</td><td>63</td></tr><tr><td>10</td><td>Appendix B: The Rigged Hilbert Space Formalism</td><td>64</td></tr></table>

**Abstract** We review the progress in atomic structure theory with a focus on superheavy elements and the aim to predict their ground state configuration and element's placement in the periodic table. To understand the electronic structure and correlations in the regime of large atomic numbers, it is important to correctly solve the Dirac equation in strong Coulomb fields, and also to take into account quantum electrodynamic effects. We specifically focus on the fundamental difficulties encountered when dealing with the many-particle Dirac equation. We further discuss the possibility for future many-electron atomic structure calculations going beyond the critical nuclear charge $Z_{crit} \approx 170$, where levels such as the 1s shell dive into the negative energy continuum ($E_{n\kappa} < -m_ec^2$). The nature of the resulting Gamow states within a rigged Hilbert space formalism is highlighted.

# 1 Introduction

The periodic table (PT) of the elements, introduced by Dmitri Mendeleev and Lothar Meyer, is based on the Pauli and Aufbau (building-up) principle [1]. Arguably, the PT is the most important and useful tool concerning the electronic structure of atoms and molecules [2–5]. Chemical and physical similarities between the elements within a group or period obtained from their measurable properties is often hailed as a building block of the PT, but these patterns also follow from the underlying electronic shell structure of the atoms. Despite many controversies concerning the PT, for example, the starting and ending points of the f-block elements, the placement of the lightest elements hydrogen and helium, observed anomalies in chemical behavior or even the shape and visual representation [4, 6–8], it is still going strong after 150 years. Furthermore, with the nuclear synthesis of the 7p block elements up to oganesson with nuclear charge $Z = 118$ [9, 10], the full 7th period of the PT is now complete. Hence, what remains to be solved is how the PT can successfully be extended both theoretically and experimentally into the superheavy element region beyond $Z = 118$ [11–15]. A progress in this direction has been made by placing the unknown elements up to nuclear charge $Z = 172$ into the Periodic Table [16, 17], see for example Fig. 1.
The existence and properties of new superheavy elements beyond oganesson depends on both nuclear and electronic structure properties [18]. There are, however, a number of open questions and major challenges to both electronic and nuclear structure theory concerning the accurate prediction of physical and chemical properties of the superheavy elements.* For example, to correctly place an element into the PT and predict its basic properties, one should gain knowledge of its atomic shell structure, such as ground and excited electronic states and underlying dominant configurations [11, 12]. In the case of dense spectra, which are prominent in open-shell systems as well as in the superheavy element region where high principal quantum number and angular momentum states are occupied, detailed knowledge of low-lying excited electronic states are required within a window of a few eV. This is often a very challenging task as both relativistic and electron correlation effects play a major role requiring sophisticated multi-reference methods at the relativistic Dirac-Coulomb-Breit level of theory. Currently, the heaviest element for which it is possible to compare theory and experiment is lawrencium ($Z = 103$) [19, 20].

Moreover, the Dirac-Coulomb Hamiltonian has its limits in strong Coulomb fields as beyond the critical nuclear charge of $Z_{crit} \approx 170$ for finite-size nuclei, the

<figure><img src="image_3.png" /><figcaption>Fig. 1 Pyykkö's periodic table extended to $Z = 172$ (with permission from PCCP [17]).</figcaption></figure>

<footnote>*Here we define the starting point of the superheavy element region at the transactinides, $Z \ge 103$</footnote>
1s electron level dives into the negative energy continuum below $E = -m_ec^2$ [21–32]. At the single-particle level of theory, the correct description and interpretation of the resulting resonances can be given in terms of Gamow states [33–37], but how such diving states can correctly and accurately be described within a multi-electron framework, and how the PT can be extended beyond the critical nuclear charge, are open questions.

At high nuclear charge, the PT is ultimately limited by the nuclear stability, not by its electronic shell structure [18, 38]. For nuclear structure theory and corresponding predictions of nuclear stability of isotopes see for example Refs. [18, 38, 39] and references therein. Here we focus solely on the discussion of relativistic electronic structure theory in the superheavy element region [14, 40–42].

The outline of this Review is as follows. We first discuss the Dirac equation and its peculiarities compared to the non-relativistic Schrödinger equation, specifically for electrons in strong Coulomb fields. We discuss the critical nuclear charge in detail to clarify the region of validity of the Dirac-Coulomb Hamiltonian and discuss how states embedded in the negative energy continuum should be interpreted. The process of spontaneous pair creation in a supercritical field is analyzed including most recent references. The importance of quantum electrodynamics (QED) effects and how these can be treated in strong Coulomb fields is outlined. The major problem of correctly describing electron correlation for the accurate prediction of electronic spectra in the superheavy element region is addressed. We review the current status of electronic structure calculations for the transactinides and discuss the placement of the elements beyond oganesson into the PT based on quantum theoretical predictions. The literature on this topic is vast [28, 43–45], including a rigorous mathematical treatment of the Dirac equation and its generalizations [29, 46–50].

# 2 The Dirac Equation in Strong Coulomb Fields

## 2.1 The QED Lagrangian

Electronic structure theory is based on the QED sector of the Standard Model of particle physics. Within the Standard Model, electrons are spin-1/2 Dirac fermions, and their dynamics is described by the QED Lagrangian density

$$ \mathcal{L}_{\text{QED}} = i\hbar c \bar{\psi}(x) \gamma^{\mu} \partial_{\mu} \psi(x) - m_e c^2 \bar{\psi}(x) \psi(x) \\ - \frac{1}{4} F_{\mu\nu} F^{\mu\nu} - e \bar{\psi}(x) \gamma^{\mu} A_{\mu}(x) \psi(x), \qquad (1) $$

where $\psi(x)$ is the field operator and $\gamma^\mu$ are the Dirac matrices. The first two terms in (1) are the kinetic and mass terms describing the free electrons with mass $m_e$, whereas the third term describes the photon field $A^\mu = (\phi, A)$, corresponding to the electromagnetic scalar and vector potentials (with $F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu$). The last term corresponds to the interaction between electrons and photons, with the elementary charge $e$ acting as the coupling constant. The interaction picture represented by Eq. (1) has been extensively used in quantum field theory and it has been demonstrated to work to astonishingly high accuracy.
It would be highly desirable to treat the QED Lagrangian for a many-electron system in an external Coulomb field to avoid divergencies that appear in perturbative treatments [51]. Such a direct treatment could in principle be performed through lattice gauge theory which is mathematically well defined [52, 53]. However, the long-range nature of the Coulomb potential, related to the zero rest-mass of the photon, currently prevents any accurate computational treatment using lattice gauge theory in finite boxes [54]. Treating the required large boxes is currently computationally too demanding. However, progress in this field has recently been made on the nuclear length scale. For instance, a combined lattice QCD+QED approach has been used to successfully calculate hadron and meson mass differences, such as the proton-to-neutron mass splitting, and its dependence on both the strong and electromagnetic coupling constants [55, 56].

## 2.2 The Many-Electron Dirac-Coulomb-Breit Hamiltonian

Atomic physics calculations are performed in the Hamiltonian formalism derived from the Langrangian (1) by a Legendre transformation [57–60]. The resulting first-quantized N-particle Hamiltonian can be written in atomic units (i.e., $\hbar = 1, e = 1, m_e = 1$) as [45, 61, 62]:

$$ 
\begin{aligned}
H_D &= \sum_{k=1}^{N} h_k + \sum_{k<l}^{N} V^{ee}(r_{kl}) + H_{QED} + H_{\text{other}}, \\
h_k &= -ic\alpha_k \cdot \nabla_k + \beta_k m_e c^2 + V(r_k),
\end{aligned} 
\quad (2) $$

where $\alpha = \gamma^0\gamma$, $\beta = \gamma^0$, $r_{kl} = |r_k - r_l|$ is the inter-electronic distance, and $h_k$ is the single-particle Dirac Hamiltonian with an external potential $V(r)$, which can be the physical nuclear potential (accounting for the finite extent of atomic nuclei), or an effective potential also including electron screening, providing a better starting point for perturbative calculations [63]. The full electron-electron interaction $V^{ee}(r_{kl})$, as derived from QED, will be discussed in Sec. 4.4.2. The electron-electron interaction is often approximated by

$$ V^{ee}(r_{kl}) = \frac{1}{r_{kl}} - \frac{1}{2r_{kl}} \left[ \alpha_k \cdot \alpha_l + \frac{(\alpha_k \cdot r_{kl})(\alpha_l \cdot r_{kl})}{r_{kl}^2} \right], \quad (3) $$

were the first term is the classical Coulomb interaction, which is the dominant contribution. The frequency-independent Breit interaction (second term) contains magnetic interactions and retardation effects up to order $1/c^2$ and is an important correction to the fine structure in atoms. Together, Eqs. (2) and (3) form the Dirac-Coulomb-Breit Hamiltonian, the starting point of most applications in relativistic electronic structure theory. The importance of the effect of the Breit contribution to the 1s shell energy of superheavy elements has been pointed out quite early [64]. QED effects, represented by $H_{QED}$, which are the focus of Sec. 4, are often included using effective Hamiltonians [65–69]. The Hamiltonian may also include additional terms, represented by $H_{\text{other}}$, such as the ones arising for example from the hyperfine structure [62, 70], the nucleus-electron Breit term [71], or from weak interactions [72].
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Schematic spectrum of the one-particle Dirac operator with potential (in SI units) $V(r) = -\frac{e^2}{4\pi\epsilon_0} \frac{Z}{r}$ showing the discrete {$\phi^d$} and positive {$\phi_+$} and negative {$\phi_-$} energy continuum states. (a) negatively charged particle of mass $m$ in a Coulomb potential with $eZ > 0$, (b) free particle ($eZ=0$), and (c) the charge conjugated case of an antiparticle of charge $+e$ and mass $m$ in a Coulomb potential with $eZ < 0$.</figcaption></figure>

It is worth mentioning that the Hamiltonian (2) is not Lorentz invariant, as the Breit operator accounts for magnetic interactions and retardation effects only to order $1/c^2$ [73]. However, the corresponding deviations are supposed to be small compared to other sources of errors, such as from the approximate treatment of electron correlation [74, 75]. For inner shells, the all-order retardation contribution may not be negligible. Including the Breit operator in a self-consistent process to obtain its contribution to all-orders can also have a strong effect [76].

In order to describe electrons in the field of high nuclear charges, one first requires a detailed understanding of the spectrum of the Dirac or Dirac-Coulomb-Breit Hamiltonian in strong Coulomb fields [22, 24, 26, 28]. One of the major differences between the (many-particle) Dirac operator and its non-relativistic counterpart, is that the Dirac operator is not bounded from below and features a continuum of negative-energy states, as shown in Fig. 2. This gives rise to difficulties with variational approaches that have plagued the atomic physics and quantum chemistry communities for a long time [77–85]. This is now seen, however, as more of a technical problem than a fundamental one† discussed in more detail in Sec. 3.

<footnote>†We distinguish between problems of fundamental nature as those where knowledge to solve a particular problem is not yet available (such as problems involving physics beyond the standard model, the foundation of quantum field theory and Haag’s theorem [86], etc.) and those where knowledge is in principle available but the solution of the problem can be very hard to obtain (such as electron correlation and QED to all orders) or can be solved based on existing theory (such as resonant states embedded in the scattering continuum).</footnote>
## 2.3 The one-particle Dirac equation

In order to solve the many-electron problem, one must first understand the single-particle case. Thus, in the following, we consider the stationary Dirac equation for a single particle.

### 2.3.1 Point nucleus and self-adjointness

In strong Coulomb fields, a difficulty arises for the Dirac equation modelled with a point nuclear charge (PNC). To illustrate this, it suffices to consider the radial form of the one-particle Dirac-Coulomb equation:

$$ \begin{pmatrix} m_e c^2 + V(r) - E_{n\kappa} & c \left(-\frac{d}{dr} + \frac{\kappa}{r}\right) \\ c \left(\frac{d}{dr} + \frac{\kappa}{r}\right) & -m_e c^2 + V(r) - E_{n\kappa} \end{pmatrix} \begin{pmatrix} P_{n\kappa}(r) \\ Q_{n\kappa}(r) \end{pmatrix} = 0, \quad (4) $$

with the corresponding four-component orbital spinor

$$ \psi_{n\kappa\mu}(r) = \frac{1}{r} \left[ \begin{array}{c} P_{n\kappa}(r)\chi_{\kappa\mu}(\theta, \phi) \\ iQ_{n\kappa}(r)\chi_{-\kappa\mu}(\theta, \phi) \end{array} \right], \quad (5) $$

where $\kappa = \pm(j + \frac{1}{2})$ for $j = \ell \mp \frac{1}{2}$. The bound state eigenvalues for the point nuclear charge, $V(r) = -Z/r$ are [87–89]

$$ E_{n\kappa}(Z) = m_e c^2 \left( 1 + \frac{(Z\alpha)^2}{[n - |\kappa| + \sqrt{\kappa^2 - (Z\alpha)^2}]^2} \right)^{-1/2}, \quad (6) $$

where $\alpha$ is the fine-structure constant ($\alpha^{-1} = 137.035999206(11)$ [90]). The solution (6) is known as the Sommerfeld fine-structure formula [91]. A historical overview is given in Weinberg’s book on the quantum theory of fields [92].

It is apparent that a problem occurs when $Z > Z_{c_p} = |\kappa|/\alpha$, as $E_{n\kappa}(Z)$ becomes imaginary [93]. The range of such large Z-values is usually referred to as the critical nuclear charge region. At the onset of the imaginary solutions, Eq. (6) simplifies to

$$ E_{n,\kappa}(Z_{c_p}) = m_e c^2 (n - |\kappa|) (n^2 - 2|\kappa|n + 2\kappa^2)^{-1/2} \ge 0, \quad (7) $$

and one obtains $E_{1,-1} = 0$ for 1s, $E_{2,-1} = E_{2,1} = m_e c^2/\sqrt{2}$ for 2s and $2p_{1/2}$ at $Z_{c_p} = 1/\alpha \simeq 137.036$, and $E_{2,-2} = 0$ for $2p_{3/2}$ at $Z_{c_p} = 2/\alpha \simeq 274.072$. The difference between the behaviour of the nonrelativistic and relativistic 1s energies with increasing nuclear charge is shown in Fig. 3.

The presence of the critical charge distinguishes the Dirac equation from the standard Schrödinger equation with a Coulomb potential of a point nuclear charge, where all values $Z \ge 1$ are allowed, although one would run into similar problems with the Schrödinger equation for potentials of the form $V(r) = -Z/r^n$ with $n \ge 2$ [95].

To treat atoms with nuclear charges beyond a certain critical charge, $Z > Z_{nsa}$, where the Dirac operator becomes non-self-adjoint (nsa), one has to carefully choose an appropriate self-adjoint extension to the basic Dirac-Coulomb operator together with the correct operator domain [29, 30, 46, 96–100]. For example, this
<figure><img src="image_5.png" /><figcaption>Fig. 3 Nuclear charge dependence of the 1s energy levels for hydrogen-like atoms at various levels of theory using the Dirac equation. If not otherwise stated the results are from Ref. [31]. The models considered are: PNC - point nuclear charge; FNC - finite nuclear charge distribution; recoil - nuclear recoil effects according to Eq. (8) [94]; recoil+VP - includes the Uehling vacuum polarization term [94]; VP+SE - includes major QED corrections from vacuum polarization and self-energy; NR - nonrelativistic results. Results are also shown for the Ar-like system with FNC.</figcaption></figure>

can be done by adding additional operators such as the nuclear recoil and Uehling terms, discussed in section 2.3.2, or by removing the problematic singularity in the Coulomb term at zero by working with a realistic finite-size nuclear charge distribution to regularize the Coulomb interaction. The mathematical problem arises due to the singularity of the Coulomb operator $-Z/r$ at the origin. As a result, the Dirac operator is not (essentially) self-adjoint anymore in the critical nuclear charge region. In fact, $H_D$ becomes non-self-adjoint [98] for a j-state at $Z \ge Z_{nsa} = \sqrt{j(j+1)/\alpha}$. For the 1s level this corresponds to $Z \ge \sqrt{3}/(2\alpha) \approx 118.677$ [96, 97, 101], which lies just above the nuclear charge of oganesson ($Z = 118$). This was pointed out as early as in 1928 by Gordon [89]. For a more rigorous mathematical analysis on the self-adjointness of the point-charge Dirac-Coulomb operator we refer the reader to Sec. 9 and the literature cited therein.

On a historical note, the onset of imaginary solutions for the Dirac equation with the bare Coulomb operator led Feynman to the conclusion that elements above $Z = 137$ should not exist. Hence, the element with nuclear charge 137 is sometimes (jokingly) called Feynmanium [102].

### 2.3.2 Nuclear Recoil and Uehling terms

For a point-like nucleus, the nuclear recoil operator can be approximated by [94]

$$H_{\text{NRB}} = -\frac{1}{2M}\Delta + i\frac{(Z\alpha)}{2Mr}\left[\boldsymbol{\alpha} \cdot \boldsymbol{\nabla} + \frac{1}{r^2}(\boldsymbol{\alpha} \cdot \mathbf{r})(\mathbf{r} \cdot \boldsymbol{\nabla})\right], \quad (8)$$
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Different contributions to the 1s energy level for hydrogen-like atoms, evaluated using the MCDFGME code [106]. Higher-order VP includes the Wichmann and Kroll (WK) correction (order $\alpha(Z\alpha)^3$) as well as approximation to the $\alpha(Z\alpha)^5$ and $\alpha(Z\alpha)^7$ potential contributions. Two-loop self-energy corrections SE-SE, SE-VP and S(VP)E are from Refs. [107–113]. Loop-after-loop VP is approximated by solving the Dirac equation including the Uehling potential. Finite nuclear size correction and uncertainties on nuclear size are from Ref. [114]. See also [115, 116] and references therein.</figcaption></figure>

where *M* is the mass of the nucleus (for a more concise QED treatment see [103]). This recoil operator can be added to the one-particle Dirac-Coulomb operator. For a more detailed discussion of nuclear recoil effects see Refs. [104, 105].

In Ref. [94], both the recoil correction and the Uehling potential $V_U$ for a point nucleus were included in the Dirac equation to see how that would change the $Z_{c_p}$. A value of $Z_{c_p}(1s) = 144$ is then obtained. These additional operators do not necessarily secure the self-adjointness of $H_D + H_{NRB}$ in the critical charge region, hence, a careful analysis of the 1s eigenfunction at the origin is still required. More details on vacuum polarization (VP) and on the Uehling potential can be found in Sec. 4.

Nevertheless, numerical calculations show that the critical charge for the 1s state before reaching the onset of the negative energy continuum is $Z_c(1s) \approx 144.75$ due to the nuclear recoil, and $Z_{c_p}(1s) \approx 143.95$ due to the nuclear recoil-plus-Uehling term as shown in Fig. 3 [94]. Moreover, the diving of the $2p_{1/2}$, $2s$ and $3s$ levels comes at nuclear charges of $Z_{c_p} \approx 146$, 165, and 193 respectively. The lifting of the $2s/2p_{1/2}$ level degeneracy due to the nuclear recoil becomes thus quite sizable at high-Z values. For $Z\alpha < 1$, the results including nuclear recoil and Uehling terms are close to the point nuclear charge (PNC) case, as is the steep descent of the energy levels towards the critical nuclear charge. On the other hand, Fig. 4 demonstrates that around $Z = 120$ the finite nuclear size correction becomes more important than that originating from the nuclear recoil and Uehling terms. This is addressed in the following section.
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Radial proton (left) and neutron (right) densities of doubly-magic nuclei $^{48}$Ca, $^{208}$Pb, $^{302}$Og, and $^{472}$164 obtained in nuclear DFT with three different energy density functionals. The shaded areas indicate the spread of DFT predictions. (Modified from [125].)</figcaption></figure>

### 2.3.3 Finite nuclear charge distributions

By considering a finite nuclear charge distribution, $\rho_N(r)$, the problematic singularity at zero is removed. As a result $H_D$ becomes self-adjoint for $Z > Z_{cp}$ with real eigenvalues and real radial functions for the discrete spectrum, and thus represents the most natural self-adjoint extension to the PNC Dirac Hamiltonian. This was already realized by Schiff, Snyder and Weinberg as early as in 1939 [93]: In all these cases where the energy cannot be brought to diagonal form, one must take into account either existing deviations from the assumed potential, such as the breakdown of the Coulomb law at small distances, or the reaction of the pair field itself on the external field.

The potential for an electron interacting with a nuclear charge distribution is given by

$$V(r) = - \int d\mathbf{R} \frac{\rho_N(\mathbf{R})}{|\mathbf{r} - \mathbf{R}|} \qquad (9)$$

The nuclear charge densities should in principle be obtained using nuclear density functional theory (DFT) based on realistic energy density functionals, see Sec. 2.4.2. To obtain the nuclear charge density from computed proton and neutron density distributions, several corrections have to be considered [117–119]. The nucleon structure is taken into account by folding with the intrinsic form factor of the free nucleons expressed in terms of the Sachs form factors [120]. The spurious center-of-mass motion can be corrected by an unfolding with the width of the centre-of-mass vibrations. Finally, one should include the contribution from the spin-orbit currents [121]. Note that, for the deformed nuclei, the spin-orbit contributions change gradually as the single-particle spin-orbit strength becomes highly fragmented by deformation and nucleonic pairing (nucleonic superconductivity) [119]. Precise nuclear charge densities are essential for interpreting atomic experiments searching for new physics [122, 123] or for studying effects related to fundamental symmetry violations [124].
Realistic nuclear modeling of charge densities is particularly important for the superheavy nuclei, the existence of which depends on the interplay between the short-ranged attractive nuclear force and long-ranged electrostatic repulsion, which rapidly grows with Z. Since the Coulomb repulsion minimizes the total binding energy of the nucleus by increasing the average distance between protons, the total energy is significantly lowered by pushing protons toward the nuclear surface. This mismatch between interaction ranges in superheavy nuclei results in Coulomb frustration effects [18, 38], which are expected to produce exotic topologies of nucleonic densities, such as voids (bubbles) or tori. Figure 5 shows the proton and neutron density distributions of several nuclei predicted by nuclear DFT [125]. The superheavy nuclei such as $^{302}$Og, and $^{472}$164 show a clear central depression in the proton density distributions resulting in a semi-bubble structure. The properties of Coulomb-frustrated superheavy nuclei, including their characteristic density distributions and shell structure, have been investigated in numerous studies, see Refs. [125–127] and references cited therein.

In the absence of predictions based on realistic nuclear models, schematic approximations for $\rho_N$ are often applied. These are sufficient for most applications in heavy element research. There is a range of nuclear charge models in use and, for several of these models, analytical expressions for the integral (9) in terms of standard functions can be found in Ref. [128]. Most implementations in numerical atomic structure programs apply the (spherical) Fermi two-parameter model [129, 130]

$$ \rho_N(R) = \frac{\rho_0}{1 + e^{(R - R_0)/a}}, \qquad (10) $$

where $R_0$ is the half-density radius, a is the diffuseness parameter, and $\rho_0$ is a normalization constant such that $\int \rho_N(\mathbf{R})d\mathbf{R} = Z$. For many nuclei, this model reasonably agrees with nuclear DFT calculations. It is to be noted, however, that a simple model like (10) is bound to fail for superheavy nuclei that exhibit appreciable Coulomb frustration effects, see Fig. 5. Nevertheless, for the valence shell this nuclear charge model should perform reasonably well even for the superheavy elements. For example, the Fermi charge distribution have been used for electronic structure calculations of Ref. [14] in the superheavy element region up to $Z = 173$.

For the homogeneous nuclear charge distribution analytical expressions for the radial Dirac components of the wave function exist. In that category of nuclear models, the simplest one is the uniformly charged spherical shell or top slice (TS) model, with the nuclear charge being smeared out over a spherical nuclear surface at radius $R_0$ [128, 131, 132]

$$ \rho(r) = \frac{Z}{4\pi r^2} \delta(r - R_0). \qquad (11) $$

It results in a potential of the form $V(r) = -Z/R_0$ for $0 \le r \le R_0$ together with the usual Coulomb term $V(r) = -Z/r$ at $r > R_0$. This approximation cuts off the problematic singularity of the Coulomb potential at nuclear radius $R_0$ and therefore secures the self-adjointness in the region $|E| < m_e c^2$ of the discrete spectrum [21]. To express the radial Dirac components analytically, one divides the solution of the Dirac equation into the two regions $[0, R_0]$ and $[R_0, \infty)$ with an additional boundary condition at $R_0$ to match the two wave functions (see also Sec. 2.4.4) [28].
The potential $V(r)$ for the TS model is, however, discontinuous in its first derivative and is therefore often extended to the homogeneously charged sphere (HCS) model of the form [133]

$$ \rho(r) = \rho_0 \Theta(1 - r/R_0), \qquad (12) $$

where $\rho_0 = 3Z/4\pi R_0^3$ and $\Theta(x)$ is the Heaviside step function [128, 132]. The resulting HCS potential is of the form $V(r) = V_0 + V_2r^2$, where $V_0 = -3Z/2R_0$ and $V_2 = V_0/2R_0^2$. This potential is discontinuous in its second derivative.‡ It is clear that by choosing $V_0 = -Z/R_0$ and $V_2 = 0$ the TS model is recovered. This results in a special case of a Fuchs-type differential equation for which analytical solutions to the Dirac equation can be formulated similarly to the procedure used for the TS model [135]. The HCS model was recently used to study isotope shifts using a modified nuclear parameter $\delta\langle r^2 \rangle \to \delta\langle r^2\gamma \rangle$, such that the electronic structure factor $\tilde{F}_i$ becomes isotope independent [136–138]. Note that this approximate expression for the energy shift is valid only when $\alpha ZR_0 \ll 1$ [139], where $R_0$ is expressed in atomic units. For nuclear charge $Z > 137$ this expression is manifestly wrong as $\gamma$ becomes imaginary.

The HCS model can be further generalized by using a Taylor expansion for the nuclear density around the origin [128]

$$ \rho(x) = \Theta(1-x) \sum_{i=0}^{n} a_i x^i, \qquad (13) $$

with $x = r/R_0$ resulting in a power series for $V(r)$. Breit introduced the simple potential $V(r) = V_0 + V_2r^n$, where $V_0 = -(n+1)Z/nR_0$ and $V_2 = Z/nR_0^{(n+1)}$ [133]. For these nuclear models one can derive the radial Dirac wave function from a polynomial expansion. [128, 135]. In the region $r < R_0$ for $\kappa > 0$ the radial wave function is expressed as [140]

$$ \begin{aligned} P_{n\kappa}(r) &= N_{n\kappa} r^{\kappa} \left\{ r - \left[ \frac{(E_{n\kappa} - V_0)(E_{n\kappa} + 2c^2 - V_0)}{2c^2(3 + 2\kappa)} + \frac{V_2(1 + 2\kappa)}{(E_{n\kappa} + 2c^2 - V_0)(3 + 2\kappa)} \right] r^3 + \dots \right\} \\ Q_{n\kappa}(r) &= N_{n\kappa} r^{\kappa} \left\{ \frac{c(1 + 2\kappa)}{(E_{n\kappa} + 2c^2 - V_0)} - \frac{E_{n\kappa} - V_0}{2c} r^2 + \dots \right\}, \end{aligned} \qquad (14) $$

and for $\kappa < 0$

$$ \begin{aligned} P_{n\kappa}(r) &= N_{n\kappa} r^{|\kappa|} \left\{ 1 - \frac{(E_{n\kappa} - V_0)(E_{n\kappa} + 2c^2 + V_0)}{2c^2(1 + 2|\kappa|)} r^2 + \dots \right\} \\ Q_{n\kappa}(r) &= N_{n\kappa} r^{|\kappa|} \left\{ -\frac{E_{n\kappa} - V_0}{c(1 + 2|\kappa|)} r + \left[ \frac{(E_{n\kappa} - V_0)^2(E_{n\kappa} + 2c^2 + V_0)}{2c^3(1 + 2|\kappa|)(3 + 2|\kappa|)} + \frac{V_2}{c(3 + 2|\kappa|)} \right] r^3 + \dots \right\}. \end{aligned} \qquad (15) $$

Since the exponents of the $r^{|\kappa|+i}$ terms in (14) and (15) are integers, there is no problem at the origin and the derivative norm exists. Furthermore, the wave function is locally absolutely continuous, unlike for the PNC case. To show this more rigorously, one applies the Weyl-Weidmann limit point - limit circle theorem [141]

<footnote>‡Because of the discontinuity in the potential at $R_0$, one has to set one of the grid points in numerical program packages at the nuclear boundary to avoid numerical instabilities [134].</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 6** The nuclear-size contributions to the ground-state energies of the Li-like uranium isotopes using a deformed Fermi model for $\rho_N$, a fitted Fermi model, and a uniform charge distribution. (From [145].)</figcaption></figure>

and shows that the Dirac operator is self-adjoint in the range $E_{n\kappa} \in [-m_e c^2, m_e c^2]$, with the Sobolev space $W_{1,2}(\mathbb{R}_+)^2$ as the natural domain of the Dirac operator. Hence, for the Dirac equation with a finite-size nuclear charge distribution, the only critical charge is at the onset of the negative energy continuum at $E = -m_e c^2$. Full analytic expressions for the Dirac wave functions for TS and uniformly charged nucleus have been derived for $s_{1/2}$ and $p_{1/2}$ orbitals and used in the evaluation of the self-energy with finite size contribution [142, 143].

Shifting from a point nucleus to a model that accounts for the finite nuclear charge distribution leads to a noticeable contribution to the total electronic energy. The difference in energy originating from the use of different nuclear charge models is far smaller [144]. For example, Fig. 6 shows the calculated ground state energy shift in Li-like uranium due to the finite nuclear charge distribution for the Fermi and uniform charge distributions [145].

When introducing a finite nuclear charge into the Dirac equation, the degeneracy between the states of the same $(nj)$ but with different $\kappa$ quantum numbers is lifted. This is most prominently seen between the $2s_{1/2}$ and $2p_{1/2}$ levels. This lifting of degeneracy already appears at the nonrelativistic level between levels of same $n$ but different $\ell$ quantum numbers, but to a much smaller extent compared to the relativistic case [128].

Figure 7 shows the energy difference $\Delta E$ between the $2p_{1/2}$ and $2p_{3/2}$ orbitals and the $2s_{1/2}$ orbital for the hydrogen-like and the Be-like state [40]. The lifting of degeneracy by the finite size of the nuclear charge for the hydrogen-like system can be qualitatively explained by perturbation theory. However, in the small region inside the nucleus, the perturbing potential is so large that a first-order calculation for high nuclear charges is insufficient [146]. In contrast to the hydrogen-like energy difference, in multi-electron systems the $2s$ shell lies below the $2p$ shell for nuclear charges up to about $Z = 120$. This comes from the different effective screening of the nucleus for these two shells, which gave rise in the early history of quantum theory to the Slater rules [147]. For nuclear charges beyond
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Orbital energy difference $\Delta E$ (in a.u.) of the $2p_{1/2}$ (purple) and $2p_{3/2}$ (orange) states relative to the $2s_{1/2}$ state (in a.u.). The dashed lines are hydrogenic energy differences. The solid lines are multi-reference energies for Be-like $J = 0$ states systems involving the major configurations $1s^2 2s^2, 1s^2 2p_{1/2}^2, 1s^2 2p_{3/2}^2$.</figcaption></figure>

$Z = 120$, the $1s^2 2p_{1/2}^2$ $J = 0$ configuration lies below the $1s^2 2s^2$ configuration, as demonstrated for the Be-like $J = 0$ state in Fig. 7 and in Ref. [40]. This is because, in strong Coulomb fields, the Coulomb operator starts to dominate over the electron-electron repulsion and the atom behaves more hydrogen-like. As a result of this effect, the $2p_{1/2}$ level dives into the negative energy continuum at a far earlier stage at $Z_c(2p_{1/2}) \approx 218$ compared to the $2s$ level at $Z_c(2s) \approx 247$ [31], see discussion in Sec. 2.4.4 for more details.

### 2.3.4 $1s$ energy level reaching the negative energy continuum

Figure 3 shows the $1s$ energy level as a function of nuclear charge for hydrogen-like systems in the FNC variant, computed using the relativistic atomic program package GRASP [148]. The calculations predict a critical charge of $Z_c(1s) = 170.161$ (170.017 including QED effects) before diving into the negative energy continuum [31].

Using different models of nuclear charge distribution, the predictions for the critical nuclear charge can vary widely between $Z_c = 164$ to 174 for the $1s$ level [149, 150], but more realistically between 168 to 172 using the uniform nuclear charge distribution and neutron numbers varying between $N = Z$ and $N = 3Z$. This is demonstrated in Fig. 8, which shows the relation between the rms nuclear charge radius $R_{ch}$ and the proton number $Z$. The critical charge as a function of $R_{ch}$ has been computed using the analytical expressions of Ref. [151]. Filled squares mark the experimental charge radii [114]. The lines denote the relation between $R_{ch}$ and $Z$ for three different neutron to proton ratios [152, 153] and the semi-empirical relation [128]. From the intercept between the dashed and dash-dotted lines, an estimate for the critical charge is $Z_c(1s) = 170.26$, with a nuclear charge radius of $R_{ch} = 7.19$ fm.
<figure><img src="image_10.png" /><figcaption>**Fig. 8** The nuclear charge radius $R_{ch}$ as a function of Z, using phenomenological expressions with different neutron/proton ratios (solid lines) and the expression by Andrae [128] (dashed line). Experimentally known charge radii [114] are marked by orange squares. The critical charge as a function of nuclear radius obtained with the analytical expression of Ref. [151] is shown by a dash-dotted line. The $Z_c (R_{ch})$ Numerical values (black square) have been obtained using the MDFGME code [14, 154] with a Fermi nuclear charge model.</figcaption></figure>

In the context of the above discussion, it is interesting to notice that because of the mass scaling $\sqrt{m}$ of the Dirac equation (4) the critical charge for muonic atoms ($m_\mu/m_e = 206.7682830(46)$ [155]) for a point nucleus is more than an order of magnitude larger $Z_{C_p}^\mu(1s) \approx 1966$ compared to the electronic case. Taking into consideration the finite nuclear radius, the critical value shifts to $Z_c^\mu(1s) \approx 2200$ [156]. As in the free-particle case, the small component becomes large and takes over for $E \to -m_e c^2$.

## 2.4 Electron states in the super-critical region

In 1969, Pieper and Greiner [135] analyzed in detail the analytical solutions for FNC models as the limit $E_{n\kappa} = -m_e c^2$ is approached for different $(n\kappa)$ states. The coefficients in the r-expansion in (14) and (15) do not exhibit any pathological behavior, but the radial functions and eigenvalues become complex in the critical region $E_{n\kappa} < -m_e c^2$ and thus lie outside the natural domain of the self-adjoint Dirac operator. As a result, the Dirac-Hamiltonian eigenstates embedded in the continuum cannot readily be reached by standard atomic structure theory. In the following, we discuss some of the approaches to deal with this problem.

### 2.4.1 Energy-projected Dirac equation

The relation between the absence of self-adjointness and the appearance of the negative energy continuum in the spectrum of the Dirac operator was studied by restricting the Hilbert space to the subspace defined by the positive energy
continuum states. This can be effectively achieved by means of the projection technique, analogous to the Feshbach projection technique [157, 158] used in the context of open quantum systems. Effectively, in this method, the negative-energy continuum space is removed [71]. The resulting single-particle Dirac Hamiltonian, the so-called no-pair external field Dirac Hamiltonian, becomes:

$$ \hat{h}^{+} = \Lambda^{+}(h_{\mathrm{D}} + V_{\mathrm{ext}})\Lambda^{+} \qquad (16) $$

where $\Lambda^+$ is the projection operator onto the free-particle positive energy subspace of the free-particle Dirac Hamiltonian $H_D^{FP}$. As long as $H_D^{FP}$ has no zero eigenvalues, the operator $\Lambda^+$ can be written as

$$ \Lambda^{+} = \frac{1}{2} \left( 1 + \frac{H_{\mathrm{D}}^{\mathrm{FP}}}{|H_{\mathrm{D}}^{\mathrm{FP}}|} \right) = \frac{\alpha \cdot \mathbf{p} + \beta mc}{\sqrt{\mathbf{p}^2 + m^2 c^2}} \qquad (17) $$

where the quotient $H_D^{FP}/|H_D^{FP}|$ is called the sign operator.§ The eigenvalues of the free-particle Dirac Hamiltonian $h_D$ are $|E| \ge m_e c^2$ [29]. The projected Dirac Hamiltonian (16) can be traced back to Bethe and Salpeter [43, 160], and is therefore sometimes referred to as the Bethe-Salpeter operator [161]. As discussed in [162], the projection operator effectively removes the pair creation and annihilation terms from the Dirac Hamiltonian, i.e., removes the coupling to the pair creation/annihilation channel.

Intuitively one would expect that various mathematical problems with the Dirac equation might disappear if the negative-energy continuum states are projected out. However, if the external field $V_{\mathrm{ext}}(r)$ is the simple $1/r$ potential corresponding to a point nucleus, $\hat{h}^+$ also has a critical charge at which it becomes non-self-adjoint, just like the standard Dirac operator. In fact, the critical charge of $\hat{h}^+$,

$$ Z_c = \left( \frac{2}{\pi} + \frac{\pi}{2} \right) \alpha^{-1} \approx 124.16, \qquad (18) $$

is lower than $\alpha^{-1} \approx 137$ [71, 161].¶ The no-pair approach based on the free Dirac Hamiltonian has therefore been criticized in Ref. [82], where it is shown that it does not prevent continuum dissolution and that projection operators from the bound Dirac Hamiltonian must be used instead. The necessity to use projection operators for correlation orbitals is shown in Ref. [76].

Unlike the Dirac equation, the no-pair operator has a lower bound in the subcritical region. This result was further refined in Refs. [163, 164], which demonstrated that the operator's eigenvalues are strictly positive [163, 164], in contrast to the point nucleus Dirac equation, for which the eigenvalues go to zero for increasing nuclear charge up to $Z\alpha = 1$.

<footnote>§The Hamiltonian $H_D^{FP}$, while similar to, is not the same as the no-pair Hamiltonian often used in relativistic quantum chemistry to avoid the continuum dissolution. In that case, the projection operator is usually constructed from the positive energy eigenstates of the full external-field Dirac Hamiltonian, and does not span quite the same space as that of free-particle states. Furthermore, the corresponding projection operators depends on the nuclear charge distribution [59, 76, 82, 159].</footnote>

<footnote>¶As discussed in Sec. 9, the Dirac equation with a $1/r$ potential has another critical nuclear charge at $Z_c = (\sqrt{3}/2)\alpha^{-1} \approx 118.68$, when the condition $||H_D\phi||_2 < \infty$ is imposed to guarantee self-adjointness. The projected equation exhibits a similar critical charge at $Z_c = (3/4)\alpha^{-1} \approx 102.78$ [71, Eq. (2.9)].</footnote>
<table><caption>**Table 1** Comparison of Dirac-Coulomb calculations with experimental values for the 1s<sub>1/2</sub> ionization potential and transition energies of <sup>238</sup>U<sup>91+</sup>. All energies in eV. The rows correspond to the standard hydrogenic Dirac-Coulomb (DC) equation with point nucleus (PNC), finite nucleus (FNC), and the free-particle projected Dirac-Coulomb (PDC) equation with a FNC approximation. The experimental values are taken from Refs. [116, 165] by picking the values with the lowest uncertainty. The homogeneous uniformly charged sphere model was used. The difference with experiment and calculation for the FNC value is due to QED corrections which are not included here.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Ionization potential</th><th colspan="2">1s<sub>1/2</sub> → 2p<sub>1/2</sub></th><th colspan="2">1s<sub>1/2</sub> → 2p<sub>3/2</sub></th></tr><tr><th><i>E</i></th><th>Δ<i>E</i><sub>exp</sub></th><th><i>E</i></th><th>Δ<i>E</i><sub>exp</sub></th><th><i>E</i></th><th>Δ<i>E</i><sub>exp</sub></th></tr></thead><tbody><tr><td>DC / PNC</td><td>132, 279.93</td><td>454.83</td><td>98, 064.45</td><td>458.84</td><td>102, 630.10</td><td>451.98</td></tr><tr><td>DC / FNC</td><td>132, 083.55</td><td>258.45</td><td>97, 872.42</td><td>266.81</td><td>102, 433.71</td><td>255.59</td></tr><tr><td>PDC / FNC</td><td>140, 474.30</td><td>8, 649.20</td><td>105, 686.10</td><td>8, 080.49</td><td>110, 767.33</td><td>8, 589.21</td></tr><tr><td>Exp.</td><td>131, 825.10 ± 4.20</td><td></td><td>97, 605.61 ± 16.00</td><td></td><td>102, 178.12 ± 4.33</td><td></td></tr></tbody></table>

The projection equation with a finite nuclear potential was initially thought to remove all the problems with the negative energy continuum. Table 2.4.1 benchmarks the no-pair approximation against Dirac-Coulomb calculations for the 1s<sub>1/2</sub> ionization potential and transition energies of <sup>238</sup>U<sup>91+</sup>. Such highly charged atoms are important for precision tests of QED [165], and QED results agree with experiments to a few eV [116]. Unlike in the Dirac-Coulomb variant, the results of the free-particle projected approach shown in Table 2.4.1 compare poorly with experiment. This indicates that the projected Dirac Hamiltonian appears to be a far worse starting point than the standard Dirac equation for further QED refinements. The reasonable choice of projection operators for the whole range of nuclear charges Z remains a challenging problem. At this stage, keeping the physically relevant negative-energy continuum and dealing with directly it seems to be a better solution. However, this requires to correctly describe resonance states with $E \le -m_e c^2$ as discussed in Secs. 2.4.2-2.4.5.

### 2.4.2 Hartree-Fock-Bogoliubov equation analogy

It is instructive to make an analogy between the one-particle Dirac-Coulomb Eq. (4) and one-quasiparticle Hartree-Fock-Bogoliubov (HFB; or Bogoliubov-de Gennes) equation used in the density functional theory (DFT) of superconductors and atomic nuclei.

The HFB equation in the coordinate representation [166, 167] can be written as:

$$ \begin{bmatrix} h - \lambda & \Delta \\ -\Delta^* & -h^* + \lambda \end{bmatrix} \begin{bmatrix} u_i \\ v_i \end{bmatrix} = E_i \begin{bmatrix} u_i \\ v_i \end{bmatrix}, \qquad (19) $$

where h is the single-particle Hamiltonian; Δ is the pairing mean-field; λ is the chemical potential (or Fermi level); $E_i$ is the quasi-particle energy; and $u_i(r, \sigma)$ and $v_i(r, \sigma)$ are the upper and lower components of quasi-particle wave functions, respectively, that depend on the spatial coordinates r and spin σ. The main DFT ingredient is the energy density functional (EDF) that depends on the particle and pair densities and currents. The mean-fields h and Δ are determined self-consistently from the one-body densities and the assumed EDF.
The quasiparticle vectors are two-component wave functions $u_i(\mathbf{r}, \sigma)$ and $v_i(\mathbf{r}, \sigma)$, which acquire specific asymptotic properties [166–169] determining the asymptotic behavior of local densities. As shown in Fig. 9, the quasiparticle energy spectrum $E_i$ of HFB consists of discrete bound states, resonances, and non-resonant continuum states. The bound HFB solutions exist only in the energy region $|E| < -\lambda$. The quasiparticle continuum with $|E| > -\lambda$ consists of non-resonant (scattering) continuum and quasiparticle resonances.

The HFB equation (19) possesses the *quasiparticle-quasihole symmetry*. Namely, for each quasiparticle state $(u_i, v_i)$ and energy $E_i$ there exists a conjugate quasihole state $(v_i^*, u_i^*)$ of opposite energy $-E_i$. That is, the spectrum is composed of pairs of states with opposite energies, see Fig. 9. The conjugate states can be related through a discrete symmetry, such as time reversal [170]. In the HFB vacuum, corresponding to even number of fermions, all negative-energy eigenstates are occupied by quasiparticles. This set of quasihole states is referred to as the Bogoliubov sea [171, 172]. It follows from the projection property of the generalized HFB density matrix that if a positive-energy one-quasiparticle state is occupied, its conjugated negative-energy partner is empty [173], and vice-versa.

The Bogoliubov sea is infinitely deep, in a full analogy with the sea of negative-energy states of the Dirac equation. In practice, since infinite sums over the Bogoliubov sea cannot be carried out when computing local HFB densities, the number of HFB-active states must be truncated. Two different ways of achieving this goal are most often implemented, namely, solution of the HFB equations in a finite Hartree-Fock space [174] and truncation of the quasiparticle space. The second method corresponds to truncating directly the quasiparticle space and using a renormalization or regularization technique to account for the truncated states [167, 169, 175–179].

The proper treatment of nuclear quasi-particle HFB continuum is important for accurate description of ground-state properties and excitations [169, 171, 177, 180, 181]. Within the real-energy HFB framework, the HFB equations must be solved by imposing the scattering boundary conditions on the quasiparticle vectors. If the outgoing boundary conditions are imposed, the unbound HFB eigenstates have complex energies; within such Gamow HFB (GHFB) approach [182] the imaginary energies are related to the particle decay width.

The quasi-particle HFB continuum can also be treated in an approximate way by means of a discretization method. The commonly used approach is to impose the box boundary conditions [169, 177, 183, 184], in which HFB eigenvectors $(u_i, v_i)$ are spanned by a basis of $\mathcal{L}^2$-integrable orthonormal functions defined on a lattice in coordinate space and vanish at box boundaries. In this approach, referred to as the $\mathcal{L}^2$ discretization, quasi-particle continuum of HFB is represented by a finite number of box states. The structure of the discretized continuum depends on the size and geometry of the box [185]. In the context of the Dirac equation, scalar confinements at the level of strong Coulomb fields need to be explored, for example within a finite element approach [186, 187].

There are two kinds of quasiparticle HFB resonances. The *particle resonances* represent metastable states that have large particle (upper) component, i.e., the

<footnote>‡Confinement potentials need to be introduced in scalar form, i.e. added to the mass term. Adding a confinement to the potential term causes the spectrum to become completely continuous [28, 29, 188].</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 9** One-quasiparticle HFB spectrum. The bound states exist in the energy region $|E| < -\lambda$, where $\lambda$ is the chemical potential (negative for a particle-bound system).</figcaption></figure>

normalization of $u_i$ is much larger than that of $v_i$. The deep-hole resonances are associated with excitations of low-lying hole states of the s.p. Hamiltonian $h$. For those states, the lower component $v_i$ dominates. The deep-hole resonances acquire decay width through the coupling to the pairing channel [168, 169].

Quasiparticle resonances can be directly calculated using coordinate-space Green's function technique [189, 190] and GHFB [182]. For approaches based on the $L^2$-discretization, approximate methods have been developed to deal with HFB resonances. Since the HFB quasiparticle resonances are highly-localized states whose energies are weakly affected by the box size, the stabilization method based on box solutions with different box sizes [177, 191] can be used to obtain the resonance energies and widths. Besides the stabilization method, a straightforward smoothing and fitting technique that utilizes the smoothed occupation numbers obtained from the dense spectrum of box states has been successfully used [177].

Summarizing this section, there are many similarities between the single-particle Dirac problem and one-quasiparticle HFB problem:

- The corresponding equations have a similar two-component form.

- In both cases, the energy spectra are symmetric with respect to zero energy. In the Dirac case, this is related to charge conjugation. In the HFB case, this is due to the quasiparticle-quasihole symmetry. For a recent discussion of particle-hole symmetries of multi-fermion systems (such as band insulators or superconductors) and the charge-conjugation symmetry of relativistic Dirac fermions, see Ref. [192].
- In both cases, the resonances can be divided into particle resonances with the upper component dominating over the lower component and the hole resonances, for which the lower component dominates. At $Z \approx Z_c$, the diving states resemble hole resonances of HFB.

- In both cases, one deals with spectra that are partly discrete and partly continuous. The continuum space contains metastable states (resonances) that are embedded in the non-resonant background.

- The Dirac and HFB spectra are bound neither from above nor from below. This leads to a variational collapse (Dirac) and difficulties with the use of the imaginary time method (for both Dirac and HFB), see, e.g., Ref. [193] for a remedy.

- In both cases, one has to deal with continuum-space truncations.

Those analogies can be helpful when tackling similar problems or interpreting similar phenomena with the Dirac equation. See also Refs. [192, 194] for relevant examples.

### 2.4.3 Perturbative approach

For narrow resonances with energies close to $E = -m_e c^2$, the energy eigenstates can be obtained perturbatively. To this end, one can employ the two-potential approach [195] to the decay of a metastable state [196, 197]. Within this method, the potential describing the decaying system can be decomposed into $V = V_0 + V'$, where $V_0$ represents the bound-state potential of a closed quantum system and $V'$ is the closing potential. When applied to the diving states, one can assume $V_0$ in the form of the Coulomb potential of the finite nuclear charge distribution with $Z_0 < Z_c$ and $V' = (Z'/Z_0)V_0$, where $Z > Z_c$ and $Z' = Z - Z_0$ [24, 28]. This decomposes the overcritical Dirac Hamiltonian into $H_D = H_{D0} + V'$. Seeking for an expression of the discrete 1s state as a solution to the overcritical Hamiltonian, the approximate eigenvector is chosen to be

$$ \psi_E(\mathbf{x}) = a(E)\psi_{1s}^0(\mathbf{x}) + \int_{-\infty}^{-m_e c^2} dE' b(E', E)\psi_{E'}^0(\mathbf{x}), \quad (20) $$

where $\psi_{1s}^0(\mathbf{x})$, the 1s bound state, and $\psi_{E'}^0(\mathbf{x})$, a continuum state with energy $E'$, are the solutions of the total Dirac equation just before diving. $a(E)$ and $b(E', E)$ are coefficients to be determined [198]. This leads to the perturbative expression for the 1s state energy embedded in the continuum

$$ E_{1s}^{\text{cr}} = E_{1s}^{0} + \Delta E_{1s} + F_{1s}(E), \quad (21) $$

where

$$ \Delta E_{1s} = \langle \psi_{1s}^{0}(\mathbf{x}) | V'(\mathbf{x}) | \psi_{1s}^{0}(\mathbf{x}) \rangle \propto Z' \quad (22) $$

and

$$ F_{1s}(E) = \int dE' \frac{|\langle \psi_{E'}^0(\mathbf{x}) | V'(\mathbf{x}) | \psi_{1s}^0(\mathbf{x}) \rangle|^2}{E - E'} \propto -Z'^2. \quad (23) $$

The dash in the integral (23) indicates the Cauchy principal value. The function $F_{1s}(E)$ in (21) introduces an energy distribution to $E_{1s}^0 + \Delta E_{1s}$ with a width of

$$ \Gamma_E = 2\pi|V_E|^2 \propto \gamma Z'^2. \quad (24) $$

The width can be interpreted in terms of the positron escape width [24, 28].
### 2.4.4 Analytical continuation

One-particle resonances embedded in the negative energy continuum can be found by extending the Dirac eigenvalue problem into the complex domain. One approach is based on solving the Dirac-equation eigenproblem with the incoming boundary condition. The resulting discrete resonant (Gamow) states have complex energies $E = E_0 + i\Gamma/2$ with a positive imaginary part $\Gamma$, termed supercritical in the remainder of this review. This interpretation differs from the usual complex-energy description of decaying Gamow states for which $E = E_0 - i\Gamma/2$. Indeed, the supercritical negative-energy electron resonances can be interpreted in terms of resonances in scattering of positive-energy positron propagating backwards in time according to the Feynman-Stückelberg interpretation [199, 200].

Complex-energy solutions for the differential equation are obtained using the appropriate boundary conditions, analogous to states in the discrete region. This has, for example, been studied for the spectrum of the Dirac equation with a spherical well potential [201–203] and for a Coulomb cut-off potential [28, 151, 194, 203, 204], for which the solutions can be analytically expressed. Alternatively, solutions to the Dirac equation can be analytically continued into the complex plane by complex scaling or by introducing a complex absorbing potential [205–209]. In the following, we discuss the complex-energy solutions following the analytical continuation approach of Ref. [204], in which the nuclear potential is assumed to be constant inside the sphere of radius $R_0$.

At distances up to a cut-off radius $R_0$, the solution to the radial Dirac equation is given by the Bessel functions:

$$ \begin{pmatrix} P(r) \\ Q(r) \end{pmatrix} = C \sqrt{\beta r} \left( J_{\pm(1/2-\kappa)}(\beta r) \frac{\mp J_{\mp(1/2+\kappa)}(\beta r)}{\beta \frac{E+m_e c^2 + \frac{Z^\alpha}{R_0}}{E+m_e c^2 + \frac{Z^\alpha}{R_0}}} \right), \quad (25) $$

where $\beta = \sqrt{(E + Z\alpha/R_0)^2 - m^2c^4}$ and upper (lower) signs correspond to $\kappa < 0$ ($\kappa > 0$). For $r > R_0$, the solutions are given by the Dirac equation with a Coulomb potential. A combination of exponential and confluent hypergeometric functions satisfy the boundary conditions [210]:

$$ \begin{pmatrix} P(E, r) \\ Q(E, r) \end{pmatrix} = \left( \frac{\sqrt{m_e c^2 + E}}{-\sqrt{m_e c^2 - E}} \right) e^{ikr} \rho^{i\tau} \begin{pmatrix} f_1(E, r) \\ f_2(E, r) \end{pmatrix} \quad (26) $$

Here, $\tau = \sqrt{(Z\alpha)^2 - \kappa^2}$, $\rho = -2ikr$, $-ik = \sqrt{(m_e c^2 - E)(m_e c^2 + E)}$, and the functions $f_i$ contain Kummer’s confluent hypergeometric functions. The full analytical form can be found in Ref. [204]. ** The poles of the S-matrix correspond to the resonant states; these are found by matching the $P/Q$ ratio of (25) and (26) at $R_0$. This results in real eigenenergies for solutions in the domain $E_0 \in [m_e c^2, -m_e c^2]$. Solutions with $E_0 \le -m_e c^2$ are embedded in the negative energy continuum, and are of the form $E = E_0 + \frac{i}{2}\Gamma$ with real energies $E_0 < -m_e c^2$ and widths $\Gamma > 0$. The states in the continuum diverge as $r \to \infty$ and are identified as Gamow wave functions, see Sec. 2.4.5 below for a detailed description. Note that at the critical energy $E = -m_e c^2$ the upper Dirac component $P$ in Eq. (26) vanishes at large distances. This means that close to $Z_c$ the diving resonances resemble the deep-hole HFB states discussed in Sec. 2.4.2.

<footnote>**A linear combination of the two-parameter Tricomia function and the exponential terms $e^{ikr}$ and $e^{-ikr}$, is given in Ref. [151].</footnote>
<figure><img src="image_12.png" /><figcaption>Fig. 11 Individual unnormalized P and Q components (a), (b) and the unnormalized density of the 1s wave function (c), (d) for an atom with Z = 185 and Rcut = 0.031 in units of h/(mc).</figcaption></figure>

Asymptotically, the resonant states $u_n(E_n, r)$ obey the outgoing (or incoming) boundary condition

$$u_n(E_n, r) \underset{r \to \infty}{\sim} O_l(k_n r) \sim e^{ik_n r} \quad (27)$$

where $k_n = \gamma_n - i\kappa_n$ (for details see Ref. [212]). As shown in Fig. 13, the bound states with $k_n = i\kappa_n$ ($\kappa_n > 0$) lie on the positive imaginary k-axis while the antibound (or virtual) states with $\kappa_n < 0$ lie on the negative imaginary k-axis. The decaying resonant states with ($\kappa_n, \gamma_n > 0$) lie in the fourth quadrant of the complex k-plane while the capturing resonant states with ($\kappa_n > 0, \gamma_n < 0$) lie in the third quadrant. The resonant-state trajectories in complex k-plane near the continuum thresholds $E = \pm m_ec^2$ have been analysed in Ref. [203].

The single-particle resonant states, augmented by complex-energy scattering continuum states $u(k, r)$ lying on the contour $\mathcal{L}$ obey the Berggren completeness relation [212]:

$$\sum_n |u_n\rangle\langle u_n| + \int_{\mathcal{L}} |u(k)\rangle\langle u(k)| dk = 1. \quad (28)$$

Since complex-energy continuum states belong to the RHS, the metric has to be generalized by introducing a biorthogonal basis for the radial wave functions. In particular, contrary to the Hilbert space situation, no complex conjugation appears in the radial wave functions of bra vectors [212, 221]. That is why the radial densities of 1s states shown in Figs. 11 and 12 are defined through squared upper and lower Dirac components [182]. Moreover, the radial integrals must be regularized as the Gamow states with $\mathrm{Im}(k) < 0$ exponentially diverge as $r \to \infty$, see Figs. 11 and 12. This can be done by various techniques [224–226], including
<figure><img src="image_13.png" /><figcaption>Fig. 12 Similar as in Fig. 11 but for an atom with Z = 300.</figcaption></figure>

the external complex scaling method [227]. The very reason for the asymptotic
growth of the Gamow state wave function at large *r* is the fact that such a state
represents the stationary approach to the intrinsically time-dependent problem of
decay. Indeed, the exponential temporal decrease of the wave function amplitude
must be complemented by its exponential spatial increase, and this assures that
the particle number is conserved [228].

It is important to note, that due to the charge conjugation property of the
Dirac equation, the appearance of the electron Gamow state in the negative-energy
continuum results in the presence of a positron resonant state in the positive-energy
continuum [203, 204]. This suggest an interpretation of diving electron states in
terms of positron scattering resonances, see Sec. 2.5.

As discussed in Sec. 2.4.2, resonances can also be described within the real-
energy framework of standard quantum mechanics. The commonly used approach
is based on the dense continuum discretization, elimination of the smooth non-
resonant background, and fitting the resonance peaks [177]. Another approach
is the stabilization method, in which resonances are extracted from phase shifts
obtained from box solutions obtained by assuming different box sizes [177, 191,
229, 230]. For very narrow resonances, perturbative methods, such as the two-
potential method of Sec. 2.4.3 can also be used.

Despite some work on resonances embedded in the continuum [28, 151, 203, 204], a direct utilization of Dirac Gamow states in atomic many-body calculations is practically nonexistent. The basic mathematical formulation rests on the rigged Hilbert space structure which comes with its own challenges. To be of use in atomic structure calculations of the superheavy elements, Dirac Gamow states need be studied within a multi-electron framework. Computing Dirac Berggren ensemble
<figure><img src="image_14.png" /><figcaption>**Fig. 13** Location of resonant states in the complex momentum plane. The Berggren completeness relation, Eq. (28), used in the decay context involves the bound states (b) lying on the imaginary k-axis, scattering states on the $\mathcal{L}$ contour (solid thick line), and resonant decaying states (d) in the fourth quadrant of the complex k-plane lying between the real axis and $\mathcal{L}$. For problems involving capture, the capturing resonant states (c) need to be considered and the scattering contour needs to be moved to the third quadrant. The antibound states (a) can be included in the generalized completeness relation, see For general expansions of the resolvent, see Refs. [214, 215]. The antibound (virtual) states (a) can be included in the generalized completeness relation; in this case the scattering contour has to be slightly deformed [222, 223].</figcaption></figure>

defined in Eq. (28), which can be used in a numerical atomic structure program packages, will offer many exciting avenues.

## 2.5 Positron production in the super-critical regime

The QED vacuum is unstable in the presence of a strong electromagnetic field above the Schwinger field limit, $E_S = m_e^2 c^3 / e\hbar = 1.32 \times 10^{-18} \text{ V m}^{-1}$ (or the equivalent intensity of $I_S = 2.3 \times 10^{29} \text{ Wm}^{-2}$) [231], and decays by emitting electron-positron pairs [232–234]. In the case of a potential barrier, it results in the much discussed and debated Klein’s paradox.

As pointed out in Ref. [235], pair production cannot be described within a one-body Dirac theory: it requires quantum field theoretical treatment within a time-dependent rigged Fock-space formalism that dynamically couples particles (electrons) and holes (positrons) in the Dirac continuum. A close non-relativistic analogy to this problem is a two-nucleon nuclear decay of a Gamow resonance [236]. A concise mathematical treatment in terms of incoming and outgoing electron-positron states is given by Rumpf, where the outgoing basis may be connected with the ingoing one by a unitary Bogoliubov transformation [237–239]. For further details see Ref. [240].
As discussed in Sec. 2.4.5, the resonance states of the supercritical Dirac equa-
tion are the Gamow states. The physical interpretation of an electron state em-
bedded in the continuum was extensively studied by the Frankfurt group [28].
According to these works, if an empty level is embedded in the negative energy
continuum, the initially neutral vacuum can spontaneously decay into a positron
and a bound electron with a supercritical energy. In such a case, an empty level in
the Dirac sea is interpreted as a positronic state, with the positron escaping the
supercritical field. After two positrons are emitted, the supercritical K-shell has
been successively filled with two electrons, and the Pauli principle prevents further
decay [23, 28, 135, 241–245].<sup>‡‡</sup> The resonance’s width has been interpreted as the
positron escape width with the characteristic time $\tau_E = \hbar/\Gamma$ for the pair creation
process.

This picture of pair creation was debated [151, 203, 246] on the basis of the
unitarity of the *S*-matrix. Indeed, the unitarity of the partial scattering matrix is
equivalent to the absence of inelastic channels, in particular, the absence of spon-
taneous electron-positron creation, in which it has been proven that for a *static*
external field the probability of pair creation is exactly zero [29, p. 298]. However,
the probability for pair creation does not go exactly to zero as the time derivative
of the external field approaches zero. Instead, in the adiabatic limit one observes a
sudden jump in the probability of adiabatic pair creation for critical fields which
may be defined as spontaneous pair creation [29]. That is, one requires only a
weak time dependence to trigger pair creation. Consequently, rather than to talk
about “spontaneous pair creation”, it has been recommended to use “adiabatic
pair creation” [247, 248]. Recently, the vacuum polarization energy decline and
spontaneous positron emission in QED under Coulomb supercriticality were ex-
plored within the Dirac-Coulomb problem with an external static or adiabatically
slowly varying spherically symmetric Coulomb potential created by a uniformly
charged sphere [249]. It was found that in the supercritical region the vacuum
polarization energy is a decreasing function of the Coulomb charge, resulting in a
decay, with a vacuum polarization energy $\varepsilon_{\text{VP}}^{\text{ren}} \sim -Z^4/R(Z)$, which provides the
required energy for positron emission ([249] Eq. (104)). Here $R(Z)$ is the nuclear
radius. The vacuum polarization and its effect on the value of supercritical $Z$ are
also studied in [250]. This debate could, however, have been avoided by referring
to Thaller’s work on scattering operators [29].

In principle, one could initiate pair creation using intensive laser fields above
the Schwinger limit $I_S$ (see Ref. [251] for a recent review). One proposal is by using
multiple<sup>§§</sup> focused beams from x-ray free electron lasers [252]. Repeated cycles of
particle creation and annihilation can take place in tune with the laser frequency
and the production of a few hundred particle pairs per laser period can occur.
As an analogous approach, Ref. [253] proposed a model of the quantum Dirac
field realized by ultra-cold fermionic atoms in an optical lattice. Here, numerical
simulations demonstrate the effect of spontaneous pair creation in the optical
analogue system. Yet another possibility is to use a strong laser beam coupled to
an atomic or molecular system with a strong Coulomb field as found for example in

<footnote>‡‡We could, in principle, excite an electron from the filled Gamow $1s_{1/2}$ state in the con-
tinuum into one of the discrete states above $-m_e c^2$. This creates another hole in the state
embedded in the negative energy continuum and the possibility for yet another pair creation.</footnote>

<footnote>§§An electromagnetic plane wave that fulfills $E^2 = B^2$ and $\mathbf{E} \cdot \mathbf{B} = 0$ cannot produce electron-positron pairs.</footnote>
graphene [246, 254–256]. A Schwinger-like production of hot electron-hole plasma in semi-metallic graphene has been claimed to be observed for the first time only very recently [257].

## 2.6 Experimental perspective: heavy-ion collisions

It was proposed, that pair creation should occur in the collision between two bare nuclei with total charge number exceeding the critical value, such as the case for two U$^{92+}$ ions with a combined nuclear charge of $Z = 184$ [28, 244]. The collision system will have a supercritical regime time for ~2.3 × 10$^{-21}$ s for the U$^{92+}$+U$^{92+}$ collision at center-of-mass energy of $E_{cm} = 740$ MeV [245] as shown in Fig. 14. The expected lifetime of the supercritical resonance state is ~392 × 10$^{-21}$ s, which is two orders of magnitude shorter than the time required for vacuum decay. The probability of pair production is therefore estimated to be around 1% for the 1s level [209]. Early attempt to observe this effect [258] using ions without a $1s_{1/2}$ hole failed. The use of cooled U$^{92+}$ ions in the ESR storage ring of GSI/FAIR [259] could allow to observe this effect for the first time. A test experiment using collisions of a Xe$^{54+}$ beam on a Xe gas jet target is underway [260].

<figure><img src="image_15.png" /><figcaption>**Fig. 14** The low-lying energy levels formed by the collision of two uranium nuclei as functions of time. The arrows *a*, *b*, and *c* denote different dynamical pair-creation mechanisms and the arrow *d* indicates the spontaneous pair creation. The 1s state dives into the negative-energy continuum for about 1 × 10$^{-21}$ s. Figure taken from [209]; see also [202].</figcaption></figure>

The spontaneous emission is not the only process that can occur during the collision.¶¶ It is generally masked by a dynamical positron emission, which is induced by the time-dependent potential of the colliding nuclei above the Coulomb

<footnote>¶¶One should not forget possible weak decay processes such as the electron capture that is a common decay mode of proton rich nuclei, albeit the time frame for weak decays is much longer than for nuclear or electronic transitions [261]. Take for example the work on relativistic quantum dynamic calculations of the probability of K-vacancy production in the Xe-Xe$^{54+}$ collision at 30 MeV [262].</footnote>
**Fig. 15** Schematic Feynman diagram for the dynamical pair creation for the (inelastic) collision of two heavy nuclei with mass numbers and nuclear charges ($A_1, Z_1$) and ($A_2, Z_2$) respectively, where the outgoing nucleus binds an electron ($Z'_2 + e^-$). Two colliding nuclei create a strong electromagnetic field, strong enough to generate an electron positron pair [267]. Colliding nuclei are represented by normal lines while wavy lines refer to virtual photons and the lines with arrows correspond to leptons (electrons and positron). The double line represents a bound electron.

barrier [209, 263–266]. In this mechanism, the two colliding nuclei create a strong electromagnetic field, strong enough to generate electron-positron pairs. The pair creation in heavy atom collisions is visualized by the Feynman diagrams in Fig. 15 [202].

While the spontaneous pair creation works only in the supercritical regime, the dynamical pair creation takes place in both subcritical [264] and supercritical modes if the collision energy is high enough [268, 269]. Experimental verification of spontaneous pair creation and the distinction from the dynamical process is however challenging as the energy-differential spectra of emitted positrons by spontaneous vacuum decay are indistinguishable from the spectra of positrons emitted by the dynamical process. There are however a range of different approaches that should make vacuum decay observable. One example is by collisions with nuclear sticking, in which nuclei are bound to each other for some period of time by nuclear forces allowing for few nuclear rotations. In this very short time frame, typically of the order of $1 \times 10^{-21}$ s to $1 \times 10^{-20}$ s [270, 271], there is an increase in pair creation probability that can only be explained with the spontaneous pair creation mechanism [27, 244, 266]. Additionally, it has been shown that the pair-production probability varies as a function of nuclear collision velocities in the supercritical and subcritical region, allowing for the detection of vacuum decay experimentally [209, 272]. The impact of the vacuum polarization on the value of $Z_c$ in the case of heavy ions collisions is considered in Ref. [250]. Moreover, it has been argued [209] that the positron spectra for symmetric collisions of heavy ions with $83 \le Z \le 96$ as a function of the collision energy should show a signature of the transition to the supercritical regime.
### 3 Multi-Configuration Dirac-Hartree-Fock

With very few exceptions [273, 274], one treats the multi-electron Dirac equation within mean-field theory, that is either at the D-HF (Dirac-Hartree-Fock) level or by using D-DFT (Dirac density functional theory) [275, 276], with the latter method being more popular in molecular calculations. It is fair to say that the accuracy of current density functional approximations cannot compete with wave-function-based methods (for a recent critical analysis on DFT see Ref. [277]), especially when QED effects need to be included. At an early stage of atomic structure calculations, however, DFT in the form of D-HF-Slater theory did play an important role as electron correlation is approximately included in such a scheme [278]. Here, we focus on modern multi-reference D-HF electronic structure theory for static correlation describing correctly the states of a given $J^\pi$ symmetry, with $J$ being the total angular momentum and $\pi$ the parity. Dynamic electron correlation and its effects on atomic structure is described in Sec. 5 below.

Like in the nonrelativistic HF case, to obtain the correct ground state symmetry and low-lying electronic transitions in open-shell cases, one requires the correct description of static correlations. In finite basis-set calculations this requires a set of Slater determinants in a multi-reference treatment within a nonrelativistic or relativistic coupling scheme. In relativistic atomic numerical program packages such as GRASP [60, 148, 279–281] or MDFGME [14, 154], this is done through linear combinations of multi-shell configurational state functions (CSF’s) within a jj-coupling scheme [45]:

$$ \Psi_i(J^\pi, M_J) = \sum_r c_{ri} \Phi_r (\gamma_\nu J^\pi, M_J), \qquad (29) $$

where the $\Phi_r$ wavefunctions share the same overall total angular momentum $J$, corresponding $M_J$, and parity $\pi$. The quantity $\gamma_\nu$ stands for all other values such as angular momentum recoupling and seniority numbers [45]. Each CSF $\Phi_r$ is a linear combination of Slater determinants

$$ \Phi_r (\gamma_\nu J^\pi, M_J) = \sum_i d_i \left| \begin{array}{c} \phi_1^i(r_1) \dots \phi_N^i(r_1) \\ \vdots \quad \ddots \quad \vdots \\ \phi_1^i(r_N) \dots \phi_N^i(r_N) \end{array} \right|, \qquad (30) $$

where $\phi$ are the Dirac four-component orbital spinors defined in Eq. (5), and the coefficients $d_i$ are determined such that the CSF is an eigenstate to both $J^2$ and $J_z$. The eigenvalues and eigenvectors (configuration mixing coefficients $c_{ri}$) are then obtained by diagonalizing the Hamiltonian matrix $H_{ij} = \langle \Psi_i(J^\pi, M_J) | H_D | \Psi_j(J^\pi, M_J) \rangle$.

Multi-reference methods (including complete active space SCF) used in the quantum chemistry community have been reviewed extensively [282–284]. A comprehensive account on MCSCF theory in relativistic atomic structure calculations (usually termed MCDHF) has been provided in a textbook [45] and several publications [285, 286]. The construction of these multi-reference functions can be a formidable task if many high angular momentum open-shell $j$-states are involved [45]. The multi-reference treatment, therefore, provides a challenge for superheavy element calculations where the electronic spectrum becomes very dense and, as
a result, the multi-reference space becomes huge.*** In addition, SCF conver-
gence problems can arise for nearly-degenerate states. Nonetheless, for few-electron
systems, high-accuracy in excitation energies can be achieved if both QED and
dynamic correlation effects are included, see Secs. 4.4 and 5, respectively. High-
accuracy atomic structure calculations are also required, for example, in the search
of physics beyond the standard model (BSM) [116, 122, 287–296].

Numerical program packages, such as MCHF for the nonrelativistic [297] case
or GRASP and MDFGME for the relativistic case, apply the finite difference
method (FDM) [298]. Alternatively, the finite element method (FEM) employ-
ing, for example, B-splines (piecewise polynomials) [186, 299–301] can be used,
as implemented for example in the program AMBiT [302]. The use of B-splines
has certain advantages in relativistic atomic structure calculations [300]. As the
radial wave functions are restricted to an interval [0, R<sub>c</sub>], the atoms are spherically
confined within a radius R<sub>c</sub> set large enough (usually around 40 au) to achieve
accurate numerical results. This discretizes the positive and negative real-energy
continuum. It thus allows for an easy implementation of projection operators [76].
This method could therefore be well suited to approximately describe diving oc-
cupied levels with E < -m<sub>e</sub>c<sup>2</sup> at charges Z > Z<sub>c</sub>.†† The virtual space created
can be used for a successive electron correlation procedure, such as configura-
tion interaction or coupled cluster or MCDF. In all these numerical procedures
one usually chooses exponentially spaced grid points (called knots in FEM) with
r = r<sub>0</sub>e<sup>t</sup>, t > 0, to describe the radial wave function φ(r) accurately in the near nu-
clear region. We note that the correct description of the wave function in the inner
core region is mandatory for the accurate treatment of relativistic effects [303, 304].
B-splines have also been used to create basis sets to perform many-body perturba-
tion theory [300, 305, 305, 306] or to do MCDF calculations, as they can be used
to implement projection operators with the nucleus and electrons average poten-
tial, and obtain correlation orbitals [76]. More recently an improved method, the
dual kinetic balance [186], has been proposed to obtain basis sets free of spurious
states.

The systems of coupled integro-differential equations obtained in multi-configuration
methods are intrinsically very non-linear. In particular exchange potentials for
correlation orbitals are inversely proportional to the square of the configuration
weight, and can then be huge. Initial configuration state functions for an SCF cal-
culation are usually obtained from either the Thomas-Fermi model or from single-
particle Dirac-Coulomb solutions using screened nuclear charges [45]. However, se-
vere convergence problems can be experienced when, for example, diffuse orbitals
are involved such as for high angular momentum functions or negatively charged
atoms, or when doing correlation calculations with highly-excited configurations.
In such cases, choosing the right initial guess becomes important. Convergence
issues within the MCDF procedure have been discussed in Refs. [40, 45, 307, 308].
In some cases the problem occurs due to the relativistic nature of the atom or
ion being studied. When going to very high-Z the angular coupling goes from
LSJ coupling to almost pure JJ coupling. In that case the weight of some of the

<footnote>***This is similar to the strong correlation problem in solid state physics to describe, for example, metallic systems.</footnote>

<footnote>††The accuracy of such a discretization procedure has been shown to be poor, when it comes to the description of narrow resonance states [177]. In such a case, the preferred method to deal with these resonances is the Gamow-state framework.</footnote>
configurations contributing to a given *LSJ* level becomes very small and severe
convergence problems are observed [40].

When the four components of the spinor in relativistic methods are each al-
lowed to vary independently, the matrix representations of the Dirac operator will
fail to give the right formal nonrelativistic limit, resulting in an energy below the
true numerical value, known as variational collapse or finite basis set disease [79].
It arises whenever one wants to expand wave functions in a given basis replac-
ing operators by their matrix representation. To prevent such an unwanted effect,
certain boundary conditions such as the kinetic balance (which is automatically
considered in numerical calculations) have to be imposed which ensures the cor-
rect relation between the large and small component [45, 79, 85]. In finite basis
set treatments of the D-HF equations, using for example Slater or Gaussian type
basis sets, small errors may nevertheless occur due to variational problems (pro-
lapse). Since the kinetic balance condition implicitly projects onto the positive
energy states, it is possible that, due to the incompleteness of the basis set, the
total energy lies below the one obtained from numerical DHF calculations [159].
This can be avoided by freezing the inner core functions such that core orbitals
are sufficiently well described, or by restricting the size of the *s* and *p* basis sets, or
by making use of specifically derived prolapse-free Gaussian basis sets [309–312].

Upon inclusion of the Breit operator in Eq. (2), coupling of the positive and
negative continuum states occurs due to electron-electron interaction, leading to
the non-existence of a discrete spectrum. This is known as continuum dissolu-
tion or the Brown-Ravenhall disease [77], and can be avoided by removing all
Slater determinants containing negative-energy orbitals using a projection opera-
tor, effectively eliminating electron-positron pair contributions [59]. The projection
operator is usually constructed from the positive energy eigenstates of the full ex-
ternal field Dirac Hamiltonian, leading to the no-pair Hamiltonian of Sec. 2.4.1.
(For a recent discussion on this topic see [313].) For the case where photon-matter
field interactions are removed, a single Slater determinant (D-HF solution) auto-
matically includes the HF projection operators on positive energy states [314], i.e.,
the low-frequency Breit interaction has been shown to cause no variational failure
when included in the iterative solutions of the D-HF equations [315, 316]. Thus,
the Breit interaction has been successfully applied perturbatively [74] as well as
in variational treatments [45, 76, 316–321], where the solutions of the Dirac-Breit-
HF equations serve as a starting point for further electron correlation and QED
treatments.

An other issue with Dirac-Fock codes is the fact that for levels originating from
the same *LS* level, they may give wrong values. It was shown in [322] that the
$2p_{1/2} - 2p_{3/2}$ fine structure energy in B-like ions and the $2p^5 J = 3/2 - 2p^5 J = 1/2$
one in F-like ions did not provide the right value for light elements. The non-
relativistic limit obtained by setting the speed of light to a high value was not
zero as it should have been. At the time the proposed solution was to remove
the energy splitting obtained for $c \to \infty$ from the relativistic value. More recently
it was shown that this effect could be handled by doing large scale correlation
calculations to obtain those level energies, including all single excitations, even
the ones obeying the Brillouin theorem [323]. The same issue was also identified
in the evaluation of forbidden transitions probabilities [324].
# 4 Quantum Electodynamic Effects

Besides the corrections stemming from relativistic electron correlation described in Sec. 3, corrections issued from bound-state quantum electrodynamics must be added to get accurate predictions. The need for such corrections was demonstrated by two famous experimental discoveries. The first discovery, made by Lamb and Retherford, was the non-degeneracy between the $2p_{1/2}$ and $2s_{1/2}$ states, in contradiction to the Dirac equation, which gives degenerate levels [325]. The second discovery, made by Kusch and Foley, was that the electron Landé *g*-factor is not exactly equal to 2 in Na and Ga [326], later understood to be due to the anomalous magnetic moment of the electron. The experimental discoveries were followed by the theoretical work of Bethe [327], Feynman [200, 328], Schwinger [234, 329–331] and Tomonaga [332], which lead to the foundation of QED, the principle of which remains unchallenged up to now [116].

The derivation of the different QED contributions starts from the QED Lagrangian (1). Several methods have been proposed to calculate all-order QED corrections which are necessary for applications to high-Z elements. However, it is not trivial to define physical particle states in the presence of an external gauge field within the framework of gauge invariant quantum field theory [240]. Pioneering works on all-order vacuum polarization [333–335] have led to the modern calculations. The first accurate all-order calculation of the $1s$ self-energy [336, 337] showed that $Z\alpha$ expansions used up to that time were non-convergent at medium- and high-Z. A first attempt to evaluate the $1s_{1/2}$ state self-energy in superheavy elements was done in Ref. [338]. It was followed by the calculation of the self-energy contribution of the $1s_{1/2}$ level for finite nuclei up to $Z = 170$ [339]. This evaluation has recently been extended to all states up to $n = 5$ and $J = 5/2$ [340]. The method described in Ref. [336] is based on the *S*-Matrix formalism, which allow a full treatment of QED corrections in one-electron systems and to calculate corrections to the electron-electron interaction in few-electron systems beyond the no-pair approximation [341, 342], provided there is a well isolated reference system. A review of QED corrections in low-Z one-electron systems can be found in Ref. [343].

The Bethe-Salpeter equation [43] is a real two-body equation that has been used to derive, for example, higher-order recoil corrections in hydrogen [344] beyond what can be obtained with the Breit equation (see, e.g., Ref. [343] and references therein). The Bethe-Salpeter equation has, however, some fundamental problems [345, 346] and it becomes soon intractable for many-electron systems. For the efficient treatment of many-electron systems one requires a Hamiltonian approach (e.g., the Dirac-Coulomb-Breit Hamiltonian as a starting point) with additional effective QED perturbation terms that describe the multi-electron system to the required accuracy. The Bethe-Salpeter equation can in principle be transformed into two independent equations that match the equations of Hamiltonian relativistic quantum mechanics [347]; there is also the quasipotential approach [348]).

Three methods have been developed to deal with bound state QED (BSQED) calculations, in particular in heavy-elements. The original one is based on the *S*-matrix formalism. A detailed description of the *S*-matrix formalism and review of QED calculations based on it can be found in Refs. [143, 349]. An overview of this method is given in subsection 4.1. Approaches capable of dealing with quasi-
degenerate reference states have been proposed by using (i) a method based on the two-times Green function [350–352] (Subsec. 4.2) and (ii) a covariant version of RMBPT based on the time-evolution operator, which allow to treat more easily degenerate and quasi-degenerate states [353–356] (Subsec. 4.3). In practice, the complexity of the involved calculation is the main limitation to the use of any of these approaches, and approximate methods had to be devised.

BSQED is usually based on the Furry bound picture [357]. The unperturbed Dirac Hamiltonian $H_D$ contains the Coulomb field of the nucleus, such that the Coulomb potential is included to all orders. The electron-electron interaction is treated as a perturbation given by the potential

$$V_{\epsilon,g} = gH_Ie^{-\epsilon|t|}, \qquad (31)$$

where g is a formal expansion parameter and the interaction Hamiltonian is

$$H_I = j^\mu A_\mu - \delta M(x), \qquad (32)$$

which contains a mass renormalization term. As the electromagnetic interactions can act at an infinite distance, the term $e^{-\epsilon|t|}$ is added to turn off adiabatically the interaction at $t = \pm\infty$ to recover the unperturbed states before and after the interaction.

The electron-positron field operators defined on an appropriate Fock-space are expanded in terms of electron and positron annihilation and creation operators,

$$\psi(x) = \sum_{E_n > -m_e c^2} a_n \phi_n(x) + \sum_{E_m < -m_e c^2} b_m^\dagger \phi_m(x), \qquad (33)$$

while the BSQED Hamiltonian is given as [240],

$$H_0 = \sum_{E_n > -m_e c^2} E_n a_n^\dagger a_n \phi_n(x) - \sum_{E_m < -m_e c^2} E_m b_m^\dagger b_m \phi_m(x), \qquad (34)$$

where $a_n$ an electron annihilation operator for an electron in state n, with energy $E_n > -m_e c^2$ and $b^\dagger$ a positron creation operator for a positron in state m with energy $E_m < -m_e c^2$. For the Gamow states that dive into the negative energy continuum and have complex energies, the formalism has to be further extended. It should be noted that the formalism in Eqs.(33) and (34) is the proper quantum-field theory replacement for the Hamiltonian with projection operators given in Eq. (16), which is based on the Dirac sea definition of the positrons. Yet, for practical applications in many-electron systems, the BSQED formalism is too difficult to use, and has not been used beyond second-order corrections.

The expressions (33) and (34) are usually formulated in terms of the positive ($E_n > 0$) and negative ($E_n < 0$) spectrum of the Dirac operator, loosely termed electronic and positronic states [143]. Such terminology originates from a free-particle QED formalism [358, 359], and was later adopted for Coulomb fields describing a point nuclear charge where the lower part of the discrete spectrum terminates at $E_n = 0$ at $Z\alpha = 1$. As already pointed out, for the general case of a finite nucleus the energy can become negative and eventually the state can dive below $E = -m_e c^2$ for $Z \approx 170$. Hence the terminology of positive and negative energy states makes only sense if one shifts the spectrum up by $m_e c^2$ where the lower continuum starts then at $E < 0$.
## 4.1 S-matrix formalism

The evaluation of the energy shift in QED for an isolated q-electron state with no real photons $|N_q; 0\rangle = |n_1, \dots, n_q; 0\rangle$, is made through the Gell-Mann and Low theorem [360, 361], symmetrized by Sucher [362]

$$ \Delta E_{N_q} = \lim_{g \to 0} \frac{i\epsilon g}{2} \frac{\partial}{\partial g} \log \langle N_q; 0 | S_{\epsilon,g} | N_q; 0 \rangle, \quad (35) $$

where the adiabatic S-matrix is given by

$$ S_{\epsilon,g} = \lim_{t \to \infty} U_{\epsilon,g}(-t, t), \quad (36) $$

and $U_{\epsilon,g}$ is the adiabatic evolution operator defined as

$$ U_{\epsilon,g}(t_1, t_2) = T e^{-i \int_{t_1}^{t_2} dt V_{\epsilon,g}(t)}, \quad (37) $$

where T is the time ordering operator.

The next step is to expand the connected adiabatic S-matrix in power of the coupling constant g as done in Refs. [143, 349, 359, 363, 364]

$$ \begin{aligned} \left. g \frac{\partial}{\partial g} \log \langle S_{\epsilon,g} \rangle_C \right|_{g=1} &= \frac{\langle S_{\epsilon,1}^{(1)} \rangle_C + 2 \langle S_{\epsilon,1}^{(2)} \rangle_C + 3 \langle S_{\epsilon,1}^{(3)} \rangle_C + \dots}{1 + \langle S_{\epsilon,1}^{(1)} \rangle_C + \langle S_{\epsilon,1}^{(2)} \rangle_C + \langle S_{\epsilon,1}^{(3)} \rangle_C + \dots} \\ &= \langle S_{\epsilon,1}^{(1)} \rangle_C + 2 \langle S_{\epsilon,1}^{(2)} \rangle_C - \langle S_{\epsilon,1}^{(1)} \rangle_C^2 \\ &\quad + 3 \langle S_{\epsilon,1}^{(3)} \rangle_C - 3 \langle S_{\epsilon,1}^{(1)} \rangle_C \langle S_{\epsilon,1}^{(2)} \rangle_C + \langle S_{\epsilon,1}^{(1)} \rangle_C^3 \\ &\quad + 4 \langle S_{\epsilon,1}^{(4)} \rangle_C - 4 \langle S_{\epsilon,1}^{(1)} \rangle_C \langle S_{\epsilon,1}^{(3)} \rangle_C - 2 \langle S_{\epsilon,1}^{(2)} \rangle_C^2 \\ &\quad + 4 \langle S_{\epsilon,1}^{(1)} \rangle_C^2 \langle S_{\epsilon,1}^{(2)} \rangle_C - \langle S_{\epsilon,1}^{(1)} \rangle_C^4, \end{aligned} \quad (38) $$

where the connected S-matrix is defined by

$$ \langle S_{\epsilon,g} \rangle_C = \langle N_q; 0 | S_{\epsilon,g} | N_q; 0 \rangle_C = \sum_j \langle S_{\epsilon,1}^{(j)} \rangle_C, \quad (39) $$

with

$$ \langle S_{\epsilon,1}^{(j)} \rangle_C = \langle N_q; 0 | S_{\epsilon,1}^{(j)} | N_q; 0 \rangle_C. \quad (40) $$

Connected diagrams are diagrams with external legs, which are bound-state wave functions like the ones in Figs. 16, 18 and 19. The disconnected diagrams, which have only closed loops, only contribute to the energy of the vacuum. Examples of disconnected diagrams for one- and two-electron systems are shown in Fig. 17. Each order in Eq. (38) has poles at $\epsilon^j$, which cancel out only if all terms of a given order are calculated simultaneously. For example, for $j = 2$ both terms of $2 \langle S_{\epsilon,1}^{(2)} \rangle_C - \langle S_{\epsilon,1}^{(1)} \rangle_C^2$ must be calculated together to cancel the $1/\epsilon^2$ pole.

The S-matrix formalism is not limited to the evaluation of QED energy shifts. It can also be used for the evaluation of radiative corrections to one- [365, 366] and two-photon [367] emission probability for example, and line shapes [368].
From the definition of the S-matrix (36) and the evolution operator (37) one obtains for the matrix element of order j:

$$S_{\epsilon,g}^{(j)} = \frac{(-ig)^j}{j!} \int d^4 x_j \dots \int d^4 x_1 e^{-\epsilon|t_j|} \dots e^{-\epsilon|t_1|} T [H_I(x_j) \dots H_I(x_1)]. \quad (41)$$

These matrix elements can be expressed in terms of the electron propagator and photon propagator. The electron propagator is connected to the Dirac bound electron Green's function by

$$ \begin{aligned} S_F (x, y) &= \langle 0 | T [\psi (x) \bar{\psi} (y)] | 0 \rangle \\ &= \left\{ \begin{array}{ll} \sum_{E_n > 0} \phi_n (x) \bar{\phi}_n (y) & t_x > t_y \\ -\sum_{E_n < 0} \phi_n (x) \bar{\phi}_n (y) & t_x < t_y \end{array} \right. \\ &= \frac{-i}{2\pi} \int_{C_F} dz G (x_2, x_1, z(1+i\delta)) \gamma^0 e^{-iz(t_2-t_1)}. \end{aligned} \quad (42) $$

The electron Green's function in (42) is the solution of [143, 349]

$$(-i\alpha \cdot \nabla_2 + V(|\mathbf{x}_2|) + \beta m - z) G(\mathbf{x}_2, \mathbf{x}_1, z(1+i\delta)) = \delta(\mathbf{x}_2 - \mathbf{x}_1). \quad (43)$$

The energies of the bound states are given by the poles of the Green's function along the real axis.

The contraction of the two photon field operators in (41) gives

$$\langle 0 | A_{\mu} (x_2) A_{\nu} (x_1) | 0 \rangle = g_{\mu\nu} D_F (x_2 - x_1) \quad (44)$$

where

$$ \begin{aligned} D_F (x_2 - x_1) &= -\frac{i}{(2\pi)^4} \int d^4 q \frac{e^{-iq \cdot (x_2 - x_1)}}{q^2 + i\delta} \\ &= \frac{1}{(2\pi i)} \int_{-\infty}^{+\infty} dq_0 H (x_2 - x_1, q_0) e^{-iq_0 (t_2 - t_1)} \end{aligned} \quad (45) $$

In Eq. (45), $H(x_2 - x_1, q_0)$ is the photon Green's function, given by

$$ \begin{aligned} H(x_2 - x_1, q_0) &= -\frac{e^{-bx_21}}{4\pi x_21} \\ x_{21} &= |x_2 - x_1|, \quad b = -i(q_0^2 + i\delta)^{\frac{1}{2}}, \quad \Re(b) > 0. \end{aligned} \quad (46) $$

As noted by Dyson, the expansion in power of $\alpha$ of Eq. (35) has a radius of convergence equal to zero [52]. The series in $\alpha$ is thus only an asymptotic series that diverges for $n \ge 1/\alpha$. Thanks to the small value of $\alpha$, this is not an issue unlike for the strong interactions.

The first-order contribution in Eq. (38), the mass renormalization term, can be written as:

$$\Delta E_n^{(1)} = \lim_{\epsilon \to 0} \frac{1}{2} i \epsilon \langle S_{\epsilon, 1}^{(1)} \rangle_c = -\delta m \int d\mathbf{x} \phi_n^\dagger(\mathbf{x}) \gamma^0 \phi_n(\mathbf{x}). \quad (47)$$

The three possible second-order connected diagrams are shown in Figs. 16 (first order, one-electron QED corrections) and 18 (electron-electron interaction). They
(a)

(b)

**Fig. 16** Bound state QED corrections of lowest order with the usual labelling of Feynman diagrams. (a) vacuum polarisation; (b) one-electron self-energy. Elementary charge *e* is included for clarity. *D*<sub>*F*</sub> and *S*<sub>*F*</sub> are the Dyson (photon) and Feynman (electron/positron) propagators respectively. The double line represents a propagator in the field of the nucleus. Ψ<sub>*nℓj*</sub> represents a bound electron wave function.

<figure><img src="image_16.png" /></figure>

**Fig. 17** Example of disconnected diagrams of order $\alpha$, which only contribute to the vacuum energy. (a): one electron case, (b): two-electron case.

originate from the second-order term in Eq. (38) which can be explicitly written
(a)

**Fig. 18** Similar as in Fig. 16 but for the electron-electron interaction Feynman diagram.

(a)

(b)

**Fig. 19** Similar as in Fig. 16 but for the second-order electron-electron interaction Feynman diagrams: (a) ladder diagram; (b) crossed diagram.

(in natural units) as

$$ \langle S_{\epsilon, 1}^{(2a)} \rangle_c = \frac{1}{(4\pi i)} \int_{-\infty}^{+\infty} dq_0 \int d^4 x_2 \int d^4 x_1 e^{-\epsilon(|t_2|+|t_1|)} e^{-iq_0(t_2-t_1)} \frac{e^{-bx_21}}{4\pi x_21} \\ \left\{ \sum_{m_2 n_2 m_1 n_1} e^{i(E_{n_2}-E_{m_2})t_2} e^{i(E_{n_1}-E_{m_1})t_1} \right. \\ \times \phi_{n_2}^{\dagger}(\mathbf{x}_2) \gamma^0 \gamma^{\mu} \phi_{m_2}(\mathbf{x}_2) \phi_{n_1}^{\dagger}(\mathbf{x}_1) \gamma^0 \gamma^{\nu} \phi_{m_1}(\mathbf{x}_1) \\ \times \langle N_q; 0 | : a_{n_2}^{\dagger} a_{m_2} a_{n_1}^{\dagger} a_{m_1} : | N_q; 0 \rangle \\ \left. - 2Tr \left[ \gamma^{\mu} \frac{-i}{2\pi} \int_{-\infty}^{+\infty} dz G(\mathbf{x}_2, \mathbf{x}_2, z(1+i\delta)) \gamma^0 \right] \right. \\ \times \sum_{nm} e^{i(E_n-E_m)t_1} \phi_n^{\dagger}(\mathbf{x}_1) \gamma^0 \gamma^{\nu} \phi_m(\mathbf{x}_1) \langle N_q; 0 | a_n^{\dagger} a_m | N_q; 0 \rangle \\ \left. + \frac{-i}{\pi} \int_{-\infty}^{+\infty} dz \sum_{nm} e^{i(E_n t_2 - E_m t_1 - iz(t_2-t_1))} \right. \\ \times \phi_n^{\dagger}(\mathbf{x}_2) \gamma^0 \gamma^{\mu} G(\mathbf{x}_2, \mathbf{x}_1, z(1+i\delta)) \gamma^0 \gamma^{\nu} \phi_m(\mathbf{x}_1) \\ \left. \times \langle N_q; 0 | a_n^{\dagger} a_m | N_q; 0 \rangle \right\}. \qquad (48) $$

Those diagrams are of the order of $\alpha/\pi$ since they have two vertices.
## 4.2 Two-times Green's function method

This method is based on the generalisation of the Green's function (42) to a system of N electrons [350, 352]. The 2N-times Green's function is defined as

$$G(x'_1 \cdots x'_N; x_1 \cdots x_N) = \langle 0 | T [\psi(x'_1) \cdots \psi(x'_N) \bar{\psi}(x_1) \cdots \bar{\psi}(x_N)] | 0 \rangle. \quad (49)$$

It can be expressed as

$$ 
\begin{aligned}
& G(x'_1, \dots, x'_N; x_1, \dots, x_N) && (50) \\
&= \frac{\langle 0 | T [\psi_{\text{in}}(x'_1) \cdots \psi_{\text{in}}(x'_N) \bar{\psi}_{\text{in}}(x_N) \cdots \bar{\psi}_{\text{in}}(x_1)] \exp\{-i \int d^4 z \, H_I(z)\} | 0 \rangle}{\langle 0 | T [\exp\{-i \int d^4 z \, H_I(z)\}] | 0 \rangle} \\
&= \left\{ \sum_{m=0}^{\infty} \frac{(-i)^m}{m!} \int d^4 y_1 \cdots d^4 y_m \langle 0 | [T \psi_{\text{in}}(x'_1) \cdots \psi_{\text{in}}(x'_N) \bar{\psi}_{\text{in}}(x_N) \cdots \bar{\psi}_{\text{in}}(x_1) \right. \\
& \qquad \left. \times H_I(y_1) \cdots H_I(y_m)] | 0 \rangle \right\} \left\{ \sum_{l=0}^{\infty} \frac{(-i)^l}{l!} \int d^4 z_1 \cdots d^4 z_l \langle 0 | T [H_I(z_1) \cdots H_I(z_l)] | 0 \rangle \right\}^{-1}.
\end{aligned}
$$

Two-times Green's function method starts by keeping only two times in Eq. (50), setting $t_1 \equiv t_2 \cdots \equiv t_N \equiv t$ and $t'_1 \equiv t'_2 \cdots \equiv t'_N \equiv t'$. This operation does not lead to any loss of information. For an isolated level a of an N-electron atom, with an unperturbed energy $E_a^{(0)}$, the energy shift is given by

$$ \Delta E_a = \frac{\frac{1}{2\pi i} \oint_{\Gamma} dE (E - E_a^{(0)}) \Delta \mathcal{G}_{aa}(E)}{1 + \frac{1}{2\pi i} \oint_{\Gamma} dE \Delta \mathcal{G}_{aa}(E)}, \quad (51) $$

where $\mathcal{G}_{aa}(E)$ is the mean value for state a of the Fourier transform of the two-times Green's function (50). A perturbation expansion of $\mathcal{G}_{aa}(E)$ in powers of the fine structure constant $\alpha$ generates results similar to those shown in section 4.1.

In the case of two quasi-degenerate levels, one can use the 4-times Green's function in a similar manner.

$$ G(x'_1 x'_2; x_1 x_2) = \frac{\langle 0 | T [\psi_{\text{in}}(x'_1) \psi_{\text{in}}(x'_2) \bar{\psi}_{\text{in}}(x_2) \bar{\psi}_{\text{in}}(x_1)] \exp\{-i \int d^4 z \, H_I(z)\} | 0 \rangle}{\langle 0 | T [\exp\{-i \int d^4 z \, H_I(z)\}] | 0 \rangle}. \quad (52) $$

In this case, the perturbation expansion is realized on the subspace containing the two quasi-degenerate levels. This procedure can formally be extended to any number of quasi-degenerate states.

## 4.3 Covariant evolution-operator procedure

The covariant evolution-operator method has been developed in Refs. [356, 369–372]. The method also starts from the evolution operator, and applies Relativistic Many-Body Perturbation Theory methods (RMBPT). It is closely related to the two-times Green's function method discussed in the previous subsection. In contrast to the S-matrix formalism, which cannot handle quasi-degenerate states due to the energy-conservation condition caused by the integration over all times, the covariant evolution operator procedure has been successfully applied to quasi-degenerate states. The method is based on the fact that at $t = 0$, the Green's
(a)

(b)

**Fig. 20** Self-energy screening diagrams of order $(\alpha/\pi)^2$. The other notations are defined in the legend of Fig. 16.

operator is equivalent to the RMBPT wave operator, which is obtained as solution of a generalized Bloch equation. Additionally, in the standard evolution operator, time runs only in the forward direction and is therefore not relativistically covariant. By allowing the time to evolve forwards as well as backward, the relativistic covariance is restored. This method allows to evaluate non-QED many body effects to high-order and to take into account first and second order QED diagrams as well, at least in simple systems.

## 4.4 Calculation of QED corrections

There are several kinds of QED corrections that need to be evaluated for computing transition energies in superheavy elements. In the first category, there are one-electron corrections like the self-energy and the vacuum polarization shown in Fig. 16. These corrections concern all atoms. The Feynman diagram for the electron-electron interaction is shown in Fig. 18. It contains the Breit interaction and all-order retardation corrections. This diagram can be iterated to provide higher-order corrections to the electron-electron interaction as shown in Fig. 19. The latter diagrams and similar ones with more photons provide QED corrections to the correlation energy. The full ladder diagram in Fig. 19(a) contains both the correlation contribution present in many-body theories like RMBPT, MCDF or RCI, and pure QED corrections, which involve positrons. The crossed diagrams in Fig. 19(b) provides pure QED corrections not included in many-body calculations. A last category of Feynman diagrams contains electron-electron interaction corrections to one-electron correction, like self-energy screening presented in Fig. 20. These two last categories concern atoms with at least two electrons.

### 4.4.1 One-electron radiative corrections

The energy shift due to one-electron radiative corrections of order $i$, corresponding to an ensemble of diagrams with $2i$ vertices, is formally of order $(\alpha/\pi)^{imec^2}$, but
after renormalization it can be written as:

$$ \Delta E_{(n,\kappa)}^{(i)} = \left(\frac{\alpha}{\pi}\right)^i \frac{(Z\alpha)^4}{n^3} F_{n,\kappa}^{(i)}(Z\alpha) m_e c^2, \qquad (53) $$

where $F_{n,\kappa}^{(i)}(Z\alpha)$ is a slowly varying function of Z for a level of quantum numbers $(n, \kappa)$. For low Z, one can write an expansion of $F_{n,\kappa}^{(i)}(Z\alpha)$ as an expansion in powers of $Z\alpha$ and $\log((Z\alpha)^{-2})$. The lower-order coefficients of this expansion can be found in Refs. [373, 374]. As shown in [336], this expansion is not convergent at medium to high-Z. For superheavy elements, we will thus only consider results evaluated to all orders in $Z\alpha$.

We now discuss the evaluation of the self-energy diagram 16(b). The $1s_{1/2}$ self-energy has been evaluated to all-orders for $5 \le Z \le 120$ for point nucleus [68, 337, 338, 375, 376]. The finite nuclear-size correction is important for high-Z and small values of the principal quantum number and for s and p states. It is negligible for larger values of $|\kappa|$. The self-energy for the $1s_{1/2}$ state has been evaluated in [142, 143, 338, 377]. In Ref. [338], it was evaluated up to $Z = 160$ and in Ref. [339] up to $Z = 170$. An extension of the evaluation of $F_{1s}^{(1)}(Z\alpha)$ for $i = 1$ for point nuclei up to $Z = 137$ has been performed recently [378] and for uniformly charged nuclei up to $Z = 135$ [379] with radii in the range 1.5 fm to 7.3 fm. The work from Ref. [68] has been extended recently to $Z = 170$ [340]. In both works, the self-energy for a given Z value is calculated for a specific nuclear size, using the Fermi model. A comparison between the different values of $F_{1s}^{(1)}(Z\alpha)$ for $i = 1$ and finite nuclear size is shown in Fig. 21. All calculations are in good agreement with each other, that is within the differences of the nuclear model and size applied.

In the case of excited states, the self-energy has been evaluated for $ns$, $np_{1/2}$, $np_{3/2}$ and $nd_{3/2}$ states for $5 \le Z \le 110$ and $2 \le n \le 5$ in Refs. [336, 375, 376, 380, 381]. The point-nucleus self-energy of $ns$, $np$ and $nd$ states up to $n = 5$ can also be found in Ref. [68] for $10 \le Z \le 120$. Reference [382] contains the values of $F_{n,\kappa}^{(1)}(Z\alpha)$ for $nd_{3/2}$ to $ng_{9/2}$ up to $n = 5$ for a point nucleus. The finite size of correction for 2s states and $2p_{1/2}$ states can be found in [142, 143, 377] for $26 \le Z \le 100$, for $10 \le Z \le 120$ in [68] and for $100 \le Z \le 170$ in [340]. The comparison between different theoretical values with and without finite size correction for 2s, $2p_{1/2}$ and $2p_{3/2}$ states is shown in Fig. 21. The divergence of the point-nucleus values when $Z \to \alpha^{-1} \approx 137$ for states with $|\kappa| = 1$ is visible for both 2s and $2p_{1/2}$ states. It is in fact even more pronounced than for the $1s$ state.

For larger values of $n$ and $Z > 120$ there are no published results for $F_{n,\kappa}^{(1)}(Z\alpha)$. A large scale effort has been recently undertaken to provide values of $F_{n,\kappa}^{(1)}(Z\alpha)$ to cover the range of interest for superheavy elements. The values of $F_{n,\kappa}^{(1)}(Z\alpha)$ with all possible $\kappa$ for all $1 \le n \le 10$ and $Z$ up to 137 have been evaluated for point nuclei [378]. The point nucleus values for all possible $|\kappa| > 1$ for $n = 5$ are plotted in Fig. 22, together with values from Refs. [68, 340] for p and d states. The functions $F_{n,\kappa}^{(1)}(Z\alpha, R)$ including finite nuclear size correction for $3 \le n \le 6$ and $|\kappa| = 1$ ($s$ and $p_{1/2}$ states) have been evaluated for $Z$ up to 135, with values of $R$, the mean spherical radius in the 1.5 fm to 7.3 fm [379].

The vacuum polarization correction of order one in $(\alpha/\pi)$, presented in Fig. 16(a) can be evaluated with good enough accuracy by using an expansion in power of $Z\alpha$. The potential of order $\alpha(Z\alpha)$, with only one interaction with the nucleus is
<figure><img src="image_17.png" /><figcaption>Fig. 21 Values of the $F^{(1)}(Z\alpha)$ function in the high-Z and supercritical region. Top: comparison between finite size values for $1s_{1/2}$. Bottom: comparison between finite-size and point nucleus values for the $n=2$ shell. References: [a]=[379], [a1]=[378], [b]=[338], [c]=[339], [d]=[68], [e]=[340].</figcaption></figure>

called the Uehling potential $V_U$. The next term in the expansion, of order $\alpha(Z\alpha)^3$ is called the Wichmann-Kroll potential $V_{WK}$. All orders calculations of the vacuum polarization have been performed in [383, 384]. The Uehling potential [385] is evaluated as

$$V_U(r) = -\frac{2\alpha Z}{3\pi r} \int_{[1,\infty)} \exp(-2\alpha^{-1}\xi r) \left(1 + \frac{1}{2\xi^2}\right) \frac{\sqrt{\xi^2 - 1}}{\xi^2} d\xi \quad (54)$$

if one treats the nucleus as a point charge. An analytical formula for the Uehling potential in terms of modified Bessel functions has been provided in Ref. [386]. The expression in (54) can be extended to a finite nuclear charge distribution [387, 388]. The Uehling potential can be added to the Dirac equation potential, providing an easy way to include the loop-after-loop vacuum polarization correction to all orders [81].
<figure><img src="image_18.png" /><figcaption>Fig. 22 Similar as in Fig. 21 but for $F_{5p_{3/2}}^{(1)}(Z\alpha)$, $F_{5d_{3/2}}^{(1)}(Z\alpha)$, $F_{5d_{5/2}}^{(1)}(Z\alpha)$, $F_{5f_{5/2}}^{(1)}(Z\alpha)$, $F_{5f_{7/2}}^{(1)}(Z\alpha)$, $F_{5g_{7/2}}^{(1)}(Z\alpha)$, and $F_{5g_{9/2}}^{(1)}(Z\alpha)$ functions in the high-Z and supercritical region.</figcaption></figure>

The Wichmann-Kroll contribution $V_{WK}$ to the vacuum polarization is of order $\alpha(Z\alpha)^3$; it can be written approximately for $r \to 0$ as [333],

$$V_{WK}(r) \approx \frac{\alpha(Z\alpha)^3}{\pi} \left[ \left( -\frac{3}{2}\zeta(3) + \frac{\pi^2}{6} - \frac{7}{9} \right) \frac{1}{r} + 2\pi\zeta(3) - \frac{\pi^3}{4} + \left( -6\zeta(3) + \frac{\pi^4}{16} - \frac{\pi^2}{6} \right) r + O(r^2) \right] \quad (55)$$

where $\zeta(n)$ is the Riemann zeta function. For more details see for example Refs. [384, 389]. An efficient numerical method to evaluate $V_{WK}$ without low-$r$ expansion is given in [390]. Additional terms on this expansion of order $\alpha(Z\alpha)^5$ and $\alpha(Z\alpha)^7$, corresponding to 5 and 7 interactions with the nucleus in the vacuum-polarization loop, are approximately known. They have been used in muonic atoms for many years [333, 391–393]. Numerical methods to evaluate them can be found in Ref. [390].

One should also add next-order terms with $i = 2$ in (53). These terms are of the order of $(\alpha/\pi) \approx 2 \times 10^{-3}$ compared to the two-vertex terms (see Fig. 4) but the leading coefficients in their $Z\alpha$ expansion can be large. These terms represent, e.g., two-loop self-energy, two-loop vacuum-polarization corrections, and mixed self-energy vacuum-polarization terms. The corresponding Feynman diagrams are shown for example in Refs. [110, 321, 349]. Some of these terms can be easily calculated, such as the Källen-Sabry contribution to the vacuum polarization [394] for which a potential is known [387]. The two-loop self-energy terms have been evaluated for one- [107, 108, 111–113, 395, 396] and three-electron atoms [397] for $30 \le Z \le 100$, but only for $n = 1$ and $n = 2$ states. Mixed self-energy vacuum polarization diagrams have been evaluated in [110, 398, 399]. Whilst these diagrams are important for inner shell electron energies, they are not expected to contribute significantly to the outer-shell energies of superheavy elements compared to correlation effects. The evaluation of the diagrams of order $(\alpha/\pi)^2$, which are easier
<figure><img src="image_19.png" /><figcaption>**Fig. 23** Dirac energy, self-energy, and vacuum polarization near $Z_c = 1/\alpha$ for a point nucleus. The Uehling and Wichmann and Kroll vacuum polarization contributions have been evaluated with the MDFGME code and the self-energy is from [378] and evaluated to all order in $Z\alpha$.</figcaption></figure>

to calculate, like the Källen-Sabry term or the loop-after-loop Uehling contribution, can provide the needed order of magnitude to assess the importance of the uncalculated ones on specific cases.

In Fig. 23 we show the evolution of the QED contributions of order $\alpha$ and of the Dirac energy for the $1s$ state as a function of $1/\alpha - Z$ with non-integer values of $Z$, to show what happens near the critical $Z_c = 1/\alpha$ in the case of a point nucleus. It shows that the self-energy becomes nearly independent of $Z$ and that the vacuum polarization becomes the dominant contribution almost one order of magnitude larger than the self-energy. The total energy to that order becomes close to $-m_ec^2$. Although the Wichmann-Kroll contribution is very small, it remains to be checked how the all-order vacuum polarization and the sum of higher-order contributions would behave. The same comparison for finite size nucleus does not show the same effect: the self-energy and vacuum polarization remain of the same size and their values are strongly reduced.

### 4.4.2 Two-electron radiative and non-radiative corrections

Concerning the calculation of atomic spectra of heavy and superheavy elements, the bottleneck in terms of accuracy in such many-electron systems still lies in the treatment of electron correlation (see discussion in Sec. 5). There are, however, mixed terms between radiative corrections and the electron-electron interaction that need to be considered. The main one is known as self-energy screening (see Fig. 20). It has been evaluated by direct calculation of the Feynman diagrams only for $n=1$ and $n=2$ states [400, 401]. These terms containing self-energy loops cannot be put into the form of an exact potential, and thus cannot be easily generalized to arbitrary atoms. It would therefore be useful to formulate an approximate QED potential that could describe the screened Lamb-shift and other atomic properties sufficiently accurately and could be successfully used in molecular calculations.
One could then introduce a (model) perturbation Hamiltonian to represent the radiative part of QED corrections of the form

$$ \Delta\tilde{H}^{\text{QED}} = V_{\text{U}} + V_{\text{WK}} + h_{\text{SE}} + h_{\text{h.o.t.}}, \qquad (56) $$

where $V_U$ is the Uehling potential, $V_{WK}$ is the Wichmann and Kroll potential, $h_{SE}$ the self energy model potential and $h_{h.o.t.}$ represents two-loops contributions. The aim of this operator is to include approximate QED corrections to the electron-electron operators, with negligible errors compared to the electron correlation treatment. In addition, the matrix elements of this QED Hamiltonian can be added to the CI matrix or to the Hamiltonian matrix and differential equation in the MCDF procedure.

The non-radiative part in Eq. (2) is dominated by the electron-electron interaction, which is obtained by evaluating the Feynman diagram of Fig. 18. It is given in atomic units and in the Coulomb gauge by

$$ V(r_{ij}, \omega_{ij}) = \frac{1}{r_{ij}} - \frac{\alpha_i \cdot \alpha_j}{r_{ij}} - \frac{\alpha_i \cdot \alpha_j}{r_{ij}} (\cos(\alpha\omega_{ij}r_{ij}) - 1) + (\alpha_i \cdot \nabla_i)(\alpha_j \cdot \nabla_j) \frac{\cos(\alpha\omega_{ij}r_{ij}) - 1}{(\alpha\omega_{ij})^2 r_{ij}}, \qquad (57) $$

where $\omega_{ij} = E_i - E_j$ is the energy of the photon exchanged between the two electrons. The $\nabla$ operator acts only on $r_{ij}$ and not on the following wave function. The Breit operator [402–404] in Eq. (3) corresponds to the expansion of Eq. (57) in powers of $\alpha = 1/c$ up to the second order. It is then independent of $\omega_{ij}$. The frequency-dependent part is called higher-order retardation. This finite frequency dependent Breit interaction has been explored in detail in many works [66, 76, 314, 405–407]. The main difficulty lies in the definition of $\omega_{ij}$ in CI or MCDF calculations, where the energy of an individual orbitals is not physical and can also reach very negative values, much lower than $-m_e c^2$ [66, 76, 407]. The gauge dependence of the resulting energy shift has been discussed in detail in Refs. [74, 408].

The second-order diagrams of Fig. 19 have been evaluated in Refs. [341, 342, 409, 410] for the ground state and $n = 2$ excited states of two-electron atoms. They contains specific QED corrections beyond what can be obtained by many-body treatment of the interaction in Eq. (57) with the necessary projection operators. These corrections contains the positron part of the ladder diagram 19 (a) and the contribution from the cross-ladder diagram 19 (b).

### 4.4.3 Effective QED Hamiltonians

To bring the self-energy term into a useful effective Hamiltonian form, $h_{SE}$, is the most challenging part as this operator is inherently non-local. Nevertheless, many attempts were made to estimate the self-energy shift in atomic spectra by approximations. Earlier ones were summarized in [411]. Approximations based on effective Z values to account for electron screening were introduced in the early versions of GRASP [412]. The Welton approximation [413] was introduced in the MDFGME code in 1987 [65] for s-states and generalized to $\ell \ge 0$ in [66]. Effective operators directly based on BSQED have been introduced more recently [68, 340].
In Ref. [411], a local self-energy potential was introduced in a simple Gaussian form

$$h_{\text{SE}}(r) = (b_0 + b_1 Z + b_2 Z^2) \exp \left\{ -(\beta_0 + \beta_1 Z + \beta_2 Z^2) r^2 \right\}, \quad (58)$$

which serves as a rough estimate. Here $b_i$ and $\beta_i$ are adjustable parameters and the Gaussian is located close to the nucleus.

A far more accurate expression for an effective self-energy Hamiltonian has been proposed in Ref. [67] to be

$$h_{\text{SE}}(r) = \Phi_{\text{mag}}(r) + \Phi_{\text{el}}(r) + \Phi_{\text{low}}(r), \quad (59)$$

with the magnetic form factor

$$\Phi_{\text{mag}}(r) = \frac{\alpha}{4\pi m} i\gamma \cdot \nabla \left[ \phi(r) \left( \int_{1}^{\infty} dt \frac{e^{2\text{trm}}}{t^2\sqrt{t^2-1}} - 1 \right) \right], \quad (60)$$

where $\phi(r)$ is the electric potential of the nucleus. The last two terms are contributions from the electric form factor decomposed into a high- and a low-frequency part,

$$\Phi_{\text{el}}(r) = A(Z) \frac{\alpha}{\pi} \phi(r) \int_{1}^{\infty} dt \frac{e^{-2\text{trm}}}{\sqrt{t^2-1}} \left[ \left(1 - \frac{1}{2t^2}\right) \left(\log(t^2-1) + 4\log\left(\frac{1}{Z\alpha} + \frac{1}{2}\right)\right) - \frac{3}{2} + \frac{1}{t^2} \right] \quad (61)$$

The (long-range) low-frequency contribution is given by

$$\Phi_{\text{low}}(r) = -B(Z)Z^4\alpha^5me^{-Zr/a_B}, \quad (62)$$

where $B(Z) = 0.074 + 0.35Z\alpha$ is a coefficient adjusted to reproduce the radiative shifts for the high Coulomb p-levels [67], and $a_B$ is the Bohr radius. This expression of the self-energy operator was implemented into the program GRASP [280] by simply replacing the Coulomb potential $-Z/r$ by its extension to the finite nucleus case [321]. Later, this has been correctly folded into the self-energy potential leading to more complicated expressions, which slightly improves the self-energy shifts [414]. To improve the SE corrections for the s-levels, and especially for the $1s_{1/2}$ level, in multi-electron systems the prefactor $A(Z)$ in (61) was chosen to be dependent on the principal quantum number $n$ [321]. More recently, both coefficients $A(Z)$ and $B(Z)$ were refitted and made dependent on the angular momentum $\ell$ [414].

To go beyond models with adjustable parameters, one can in principle go back to first principle QED and use a spectral decomposition of the self-energy operator

$$h_{\text{SE}}^{\text{nl}} = \sum_{i,j} |\psi_i\rangle D_{i,j}^{\text{SE}} \langle\psi_j| \quad (63)$$

where $\{\psi_i\}$ represents a complete set of hydrogenic wave functions (including both continua), and the matrix elements $D_{i,j}^{\text{SE}}$ need to come from accurate self-energy calculations. This has been explored [415] with limited success because of the basis set restrictions imposed and the underlying slow convergence of this sum, which is well known from direct exact QED evaluations [143, 336, 337]. Furthermore, the off-diagonal elements $D_{i,j}^{\text{SE}}$ with $i \neq j$ are crucial and cannot be neglected. To
<figure><img src="image_20.png" /><figcaption>Fig. 24 Energy levels for the dominant configurations of the Group 10 elements Ni, Pd, Pt, and Ds. The values for Ni, Pd, and Pt are from the NIST database [419]. The Ds levels are from Ref. [42]. Different colors are used to distinguish between the three different configurations: green [(n-1)d⁸ ns²], blue [(n-1)d⁹ ns] and black [3d¹⁰]. For Pd, there are intruder states (not shown here) arising from the [(n-1)d⁹ np] configuration (for Pt from the [(n-1)d⁹ np] and [(n-1)d⁸ ns np] configurations), which mix with several of the low energy states shown. Thus, some configuration assignments (especially for the ³P₀ level) are approximate at best. For Ds, a dense spectrum arising from the [6d⁷ 7s² 7p] configuration intrudes into the normal spectrum and only few predicted lines of even parity are shown [42]. Adapted from Ref. [1].</figcaption></figure>

this end, an additional exponential type semi-local operator has been added to reduce the matrix elements $D_{ij}$ in size for the subsequent spectral decomposition [68, 416],

$$h_{SE} = h_{SE}^{sl} + h_{SE}^{nl}, \qquad (64)$$

with

$$h_{SE}^{sl}(r) = \sum_{\kappa} V_{\kappa}(r) P_{\kappa}. \qquad (65)$$

The semi-local operator

$$V_{\kappa}(r) = A_{\kappa} e^{-r/\lambda_{c}} \qquad (66)$$

differentiates between the different $\kappa$ states through the projection operator $P_\kappa$ (for the definition see Ref. [416]). Here $\lambda_c$ is the Compton wavelength. For details see Refs.[68, 416]. A recent extension to superheavy elements up to nuclear charge Z = 170 has been carried out in Ref. [340]. This scheme gives very accurate results for the self-energy. One wonders if a semilocal ansatz in the same form of a pseudopotential applied commonly in electronic structure theory could be efficiently used as well [417, 418] for all-electron QED treatments in molecules for example. It would certainly be an improvement to the original local ansatz [411] and possibly of sufficient accuracy in molecular calculations.
## 5 Electron Correlation

The accurate computational treatment of both static and dynamic electron correlation in atomic open-shell multi-electron systems is a daunting task. Even for the lightest elements such as nickel, a correct description of the many low-lying states arising from the $3d^8 4s^2$, $3d^9 4s$ and $3d^{10}$ configurations is currently not available. For example, using Gaussian type basis sets (GTOs), Ref. [420] applied large-scale complete active space second-order perturbation theory (CASPT2) calculations including relativistic corrections for nickel correlating 18 electrons within 14 orbitals. These calculations resulted in the following excitation energies with respect to the $^3D(d^9s)$ ground state (j-averaged experimental values set in parentheses): $^3F -0.08$ eV (0.03 eV), $^1D 0.32$ eV (0.33 eV), $^1S 1.77$ eV (1.74 eV). Figure 24 shows the energy levels for the Group 10 elements. From this it is clear that the correct prediction of the ground state symmetry is difficult for atoms with dense spectra. This problem will become worse when degenerate high angular momentum states are involved such as in the lanthanides and actinides and for the superheavy elements.

One of the main workhorses in relativistic atomic structure theory is the configuration interaction (CI) method with a predetermined set of CSFs where the radial shapes of the one-electron orbital spinors remain unchanged. In a typical CI procedure, the active virtual and core space are systematically increased and higher angular momentum functions added to test convergence against the final value. The resulting CI wave functions are then used for calculating QED effects, albeit QED matrix elements can be directly added to the CI matrix resulting in correlated QED calculations (this still needs to be explored for the SHEs). These CI techniques are invaluable for obtaining accurate properties to, for example, test the standard model [421]. However, as the size of a CI calculation scales exponentially with the excitation level (the number of determinants is $N_{det} \sim n^m N_v^m / (m!)^2$ with n being the number of electrons, $N_v$ the number of virtual orbitals and m the excitation level), the CI method is often combined with many-body perturbation theory for electron correlation (CI+MBPT) to allow for an efficient treatment of important core excitations [422]. Again, because of the large computer time involved one rarely goes beyond second-order MBPT, although calculations for atoms with one valence electron (Cs and Tl for example) have been performed up to third order [423, 424]. A mix of MBPT and CC methods has also been used for evaluating electron affinities for Ca and Sr [425].

A very popular electron correlation method within the quantum chemistry community is coupled cluster (CC) theory originally proposed by Coester and Kuemmel [426] for nuclear interactions, and subsequently brought into electronic structure theory [427]. There are several excellent papers, books and reviews on CC applications [428–435]. In CC theory, the ground state wave function is related to the DHF ground state configuration by an exponential operator containing the cluster operator $T$,

$$ \Psi_0 = e^T \Psi_0^{\text{DHF}}, \qquad (67) $$

with $T = T_1 + T_2 + ...$ and $T_n$ are the n-particle excitation operators. For example, if one restricts to double excitations only ($T = T_2$, CCSD), the cluster operator is

$$ T_2 = \sum_{i<j,r<s} t_{ij}^{rs} c_r^\dagger c_s^\dagger c_i c_j, \qquad (68) $$
where $t_{ij}^{rs}$ are the coupled-cluster amplitudes, determined through a variational procedure, and $c_r^\dagger$ and $c_i$ are electron creation and annihilation operators for single-particle states $r$ (virtual) and $i$ (occupied), respectively.

This scheme can be extended from the Hilbert space to the Fock space formalism (FSCC) where in addition electrons are removed (ionization) or added (attachment) [436–440]. FSCC theory has been successfully applied for atomic properties of heavy and superheavy elements [19, 441–445]. One may like to chose a multireference wave function instead of $\Psi_0^{\text{DHF}}$ for the starting point in coupled-cluster theory, termed multi-reference coupled-cluster (MRCC) theory [446–453], which, however, has its computational challenges [454]. Because of the steep computational scaling, single-reference coupled-cluster theory is usually restricted to CCSD(T) (the golden standard of quantum chemistry), where the single and double contributions to the coupled-cluster amplitudes are determined variationally, and the triples are obtained by perturbation theory. Most coupled-cluster calculations are using GTOs as the underlying basis set, but the coupled-cluster scheme can easily be adapted to numerical procedures such as FEM using B-splines [452, 455, 456].

In atomic structure theory a popular approximation to CC theory is the so-called all-order CI method (AOCI) where one keeps only the linear terms in the expansion, i.e., $\Psi_0 = (1 + T_1 + T_2 + ... )\Psi_0^{\text{DHF}}$ [457–460]. For core excitations, one can add MBPT (AOCI+MBPT) [457]. This method has successfully been used for predicting spectra of multi-electron systems with large number of electrons [460, 461]. For an alternative method combining configuration interaction with perturbation theory with a large number of valence electrons see Ref. [462]. However, high accuracy electronic structure calculations of spectra with uncertainties below $1 \times 10^{-3}$ eV are usually limited to few electron systems [66, 463–466]. Such calculations are important to test QED and subsequently the standard model as well [116, 287, 467–470]. Despite the computational limitations and challenges to treat relativistic many-electron systems with a high number of electrons, the ionization potential and electron affinity have recently been obtained for gold to meV accuracy compared to experiment [75]. Gold is a special case that can still be treated accurately as the ionization and electron attachment involves mainly the valence 6s shell. The computational cost comes, however, from a very soft polarizable 5d core giving rise to large core and core-valence correlation effects. For these calculations, relativistic CC theory up to pentuple excitations were required, that is single reference DHF-CCSDTQ(P) theory including Breit and lowest-order QED interactions, to achieve almost (but not quite) experimental accuracy [75].

For systems with two or more electrons in open-shells, as this is the case for most actinide and transactinide elements, the electronic many-body treatment remains a major challenge for the foreseeable future. For example, using MRCI plus relativistic corrections for the ionization potential of neutral uranium, $U[5f^36d^17s^2](^5L_6) \rightarrow U^+[5f^37s^2](^4I_{9/2})$, gave a value of 6.062 eV (at the DHF+Breit level one gets 5.540 eV) [471] compared to the experimental value of 6.19405(6) eV [472] achieving practically an accuracy in the 0.1 eV region. Such calculations are computationally expensive and often may not be sufficient to predict the correct sequence of states within a window of about 1 eV. For comparison, pseudopotentials which replace the core by an effective Hamiltonian give results accurate to about 0.1 eV for atomic spectra [473, 474], and are therefore not al-
ways suitable for applications in atomic physics if high accuracy is required. An all-order correlation potential method based on Green's functions has been developed [475, 476] and applied to many-electron systems into the superheavy element region. Future developments may include density matrix renormalization group (DRMG) methods [477], or even machine learning algorithms to estimate the electron correlation error compared to experiment. It is clear that for dense spectra, efficient electron correlation methods and algorithms need to be further developed to efficiently study atoms with a high number of electrons (such as the superheavy elements) including QED effects [478, 479]. An additional difficulty comes from the very large size reached by the calculation of configuration with two-large angular momenta orbitals like $6d^n 5g^m$ for example. The $6d^9 5g^9 J = 2$ LSJ configuration has 1895 JJ configurations and $\approx 30.300$ determinants, leading to $\approx 1.3 \times 10^7$ Coulomb integrals, $\approx 3.8 \times 10^7$ magnetic integrals and $\approx 3.1 \times 10^7$ retardation integrals. The $6d^{65} 5g^{12} J = 0$ LSJ configuration leads to 3360 JJ configurations and $\approx 256.000$ determinants and the number of angular integrals cannot be indexed by a 32 bit integer. Such unusual configurations, including for example the $5g^{18} J = 0$ configuration, may become competing candidates for the ground state of the Og isoelectronic sequence at very large Z.

Finally, the question arises if one should include the negative energy continuum (NEC) in the electron correlation procedure. For the lighter elements, the gap from the lowest bound state to the negative energy continuum is almost $E = 2m_e c^2$, larger than the gap to the positive energy continuum. The correlation term estimated perturbatively turned out to be of the order of $E_{cor}(\text{NEC}) \sim (Z\alpha)^3$ [300]. For the $Z = 50$ He-like ion, $E_{cor}(\text{NEC})= 0.00471$ eV [300], which implies that this term needs to be included for high precision tests on few electron systems with high nuclear charge. We need to be reminded that in the $Z \to \infty$ nonrelativistic limit the electron correlation contribution for a He-like atom is -0.046663254 a.u. [480]. For the relativistic case, this limit is not accurately known [313], see also Ref. [481] for a detailed discussion. However, perturbation theory will eventually break down if the occupied level comes close to the negative energy continuum, or even dives into it. Ref. [482] studied the effect of removing the no-virtual-pair approximation on the correlation energy of the He isoelectronic sequence. They showed that for He-like Ds ($Z = 110$) the correlation energy changes from -2.019 eV to -1.391 eV due to the NEC inclusion. The use of projection operators using a B-spline basis set build with direct Dirac-Fock potentials and their effect on the correlation energy were studied in [76] for the ground state of He-like ions. The larger impact on the relativistic correlation energy was shown to be due to the magnetic and retardation interaction. To illustrate the effect near $Z_c$ we have used the MDFGME code (2022 version) to calculate more precisely the correlation energy for the $1s^2 \, 1S_0$, with a fully relaxed wave function including 44 configurations, from $1s^2$ to $7i^2$ and projection operators. The Breit interaction was included in the self-consistent process. The relativistic Coulomb contribution changes sign around $Z = 90$. The magnetic contribution is largely dominant, with a value of -80 eV at $Z = 170$, while the retardation contribution going to 25 eV and the Coulomb part to 15 eV as can be seen in Fig. 25. It should be noted that a phenomenological inclusion of the negative energy continuum is not really consistent as the crossed diagram of Fig. 18, not included in the correlation contribution, is expected to contribute.
<figure><img src="image_21.png" /><figcaption>**Fig. 25** He-like ions contribution to the correlation energy obtained with a fully MCDF calculations using B-spline basis set and the full operator from Eq. (48). Orbitals up to 7i have been included.</figcaption></figure>

# 6 Atomic Structure Calculations of the Superheavy Elements

## 6.1 Dominant Electron Configurations

Electron configurations are needed for placing the elements into their correct place in the PT and to discuss their chemical behaviour [1]. Systematic D-HF calculations of total atomic energies of ground state configurations and symmetries (within the *jj*-coupling scheme) for the Li ($Z = 3$) to Db ($Z = 105$) isoelectronic series up to Og, including QED effects, have been provided in Ref. [483]. Predicted dominant configurations, ionization potentials, electron affinities and dipole polarizabilities for the transactinides, $Z = 102-122$, are collected in Table 2. Concerning the elements with nuclear charge $Z = 105-110$, the Table does not distinguish between the $j = 3/2$ and $5/2$ occupations for the 6d level. Moreover, whenever spin-orbit coupling becomes large, the assigned (nonrelativistic) LS symmetry has to be taken with some care. For example, the ground-state electron configuration of Fl has a $J = 0$ spin and positive parity, and one expects strong mixing between the $^3P$ and $^1S$ states due to spin-orbit coupling. Furthermore, for a variety of elements, the electronic ground and low lying excited configuration states mix, making it difficult to unambiguously assign a ground-state configuration.

Table 3 lists the ground state linear combination of the CSFs for a selection of elements, obtained within a multi-reference treatment using GRASP [148]. For Rf, there is a single dominant configuration, $6d_{3/2}^2$ in contrast to Db, for which a strong mixing is predicted between the $6d_{3/2}$ and $6d_{5/2}$ levels (the two CSFs with identical configurations have different seniority numbers $\nu=0$ and 2 for the $6d_{5/2}^2$ occupation [45]). Also for Ds (2 holes in the 6d shell) a substantial mixing between two configurations is expected. Less mixing is predicted between the strongly spin-orbit separated $7p_{1/2}$ and $7p_{3/2}$ levels in Fl, Mc, and Lv. Hence, the ground states
<table><caption>**Table 2** Atomic number *Z* and element symbol *E* according to IUPAC and predicted ground-state properties: dominant valence shell configuration and term symbol for the ground electronic state, ionization potential *I*<sub>p</sub> and electron affinity *E*<sub>A</sub> in eV, and dipole polarizability *α*<sub>D</sub> in atomic units. If possible, the most accurate value was taken from the literature. For error estimates and methods used see the cited references.</caption><thead><tr><th>Z</th><th>E</th><th>configuration</th><th>I<sub>p</sub></th><th>E<sub>A</sub></th><th>α<sub>D</sub></th><th>Refs.</th></tr></thead><tbody><tr><td>102</td><td>No</td><td>[Rn]5f<sup>14</sup>7s<sup>2</sup></td><td>6.626</td><td>-</td><td>111</td><td>[484–486]</td></tr><tr><td>103</td><td>Lr</td><td>[No]6d<sup>3</sup><sub>3/2</sub>(<sup>2</sup>D<sub>3/2</sub>)</td><td>4.96</td><td>-</td><td>323</td><td>[19, 20, 487]</td></tr><tr><td>104</td><td>Rf</td><td>[No]6d<sup>5</sup><sub>3/2</sub>(<sup>3</sup>F<sub>2</sub>)</td><td>6.01</td><td>-</td><td>115</td><td>[487, 488]</td></tr><tr><td>105</td><td>Db</td><td>[No]6d<sup>2</sup><sub>2</sub>(<sup>3</sup>F<sub>3/2</sub>)</td><td>6.814</td><td>1.189</td><td>42.5</td><td>[489–491]</td></tr><tr><td>106</td><td>Sg</td><td>[No]6d<sup>4</sup>(<sup>5</sup>D<sub>0</sub>)</td><td>7.7</td><td>-</td><td>40.7</td><td>[137, 489]</td></tr><tr><td>107</td><td>Bh</td><td>[No]6d<sup>5</sup>(<sup>6</sup>S<sub>5/2</sub>)</td><td>8.6</td><td>-</td><td>38.4</td><td>[137, 489]</td></tr><tr><td>108</td><td>Hs</td><td>[No]6d<sup>6</sup>(<sup>5</sup>D<sub>4</sub>)</td><td>9.5</td><td>-</td><td>36.2</td><td>[137, 489]</td></tr><tr><td>109</td><td>Mt</td><td>[No]6d<sup>7</sup>(<sup>4</sup>F<sub>9/2</sub>)</td><td>10.4</td><td>-</td><td>34.2</td><td>[137, 489]</td></tr><tr><td>110</td><td>Ds</td><td>[No]6d<sup>8</sup>(<sup>3</sup>F<sub>4</sub>)</td><td>9.562</td><td>0.830</td><td>32.3</td><td>[491]</td></tr><tr><td>111</td><td>Rg</td><td>[No]6d<sup>4</sup><sub>3/2</sub>d<sup>5</sup><sub>5/2</sub>(<sup>2</sup>D<sub>5/2</sub>)</td><td>11.03</td><td>1.97</td><td>30.6</td><td>[488, 489, 492]</td></tr><tr><td>112</td><td>Cn</td><td>[No]6d<sup>10</sup>(<sup>1</sup>S<sub>0</sub>)</td><td>12.02</td><td>0</td><td>27.64</td><td>[488, 492, 493]</td></tr><tr><td>113</td><td>Nh</td><td>[Cn]7p<sup>1</sup><sub>1/2</sub>(<sup>2</sup>P<sub>1/2</sub>)</td><td>7.49</td><td>0.73</td><td>29.85</td><td>[492, 494–496]</td></tr><tr><td>114</td><td>Fl</td><td>[Cn]7p<sup>2</sup><sub>1/2</sub>(<sup>1</sup>S<sub>0</sub>)</td><td>8.65</td><td>0</td><td>30.59</td><td>[488, 492, 493, 497]</td></tr><tr><td>115</td><td>Mc</td><td>[Cn]7p<sup>2</sup><sub>1/2</sub>p<sup>1</sup><sub>3/2</sub>(<sup>2</sup>P<sub>3/2</sub>)</td><td>5.574</td><td>0.313</td><td>70.5</td><td>[444, 498]</td></tr><tr><td>116</td><td>Lv</td><td>[Cn]7p<sup>2</sup><sub>1/2</sub>p<sup>2</sup><sub>3/2</sub>(<sup>3</sup>P<sub>2</sub>)</td><td>6.855</td><td>0.776</td><td>-</td><td>[444]</td></tr><tr><td>117</td><td>Ts</td><td>[Cn]7p<sup>2</sup><sub>1/2</sub>p<sup>3</sup><sub>3/2</sub>(<sup>2</sup>P<sub>3/2</sub>)</td><td>7.654</td><td>1.602</td><td>76.3</td><td>[444, 499]</td></tr><tr><td>118</td><td>Og</td><td>[Cn]7p<sup>6</sup>(<sup>1</sup>S<sub>0</sub>)</td><td>8.888</td><td>0.076</td><td>57.98</td><td>[500–503]</td></tr><tr><td>119</td><td>Uue</td><td>[Og]8s<sup>1</sup>(<sup>2</sup>S<sub>1/2</sub>)</td><td>4.783</td><td>0.663</td><td>169.7</td><td>[443, 504, 505]</td></tr><tr><td>120</td><td>Ubn</td><td>[Og]8s<sup>2</sup>(<sup>1</sup>S<sub>0</sub>)</td><td>5.851</td><td>0.021</td><td>162.6</td><td>[445]</td></tr><tr><td>121</td><td>Ubu</td><td>[Ubn]8p<sup>1</sup><sub>1/2</sub>(<sup>2</sup>P<sub>1/2</sub>)</td><td>4.447</td><td>-</td><td>-</td><td>[488]</td></tr><tr><td>122</td><td>Ubb</td><td>[Ubn]7d<sup>1</sup><sub>3/2</sub>8p<sup>1</sup><sub>1/2</sub>(<sup>1</sup>D<sub>2</sub>)</td><td>5.651</td><td>-</td><td>-</td><td>[488]</td></tr></tbody></table>

<table><caption>**Table 3** Selected ground state CSF's obtained within a multi-reference DHF treatment using GRASP.</caption><thead><tr><th>Z</th><th>E</th><th>CSF</th></tr></thead><tbody><tr><td>104</td><td>Rf</td><td>ψ<sub>J=2</sub> = 0.9300φ(6d<sup>2</sup><sub>3/2</sub>) + 0.0350φ(6d<sup>1</sup><sub>3/2</sub>d<sup>2</sup><sub>5/2</sub>) - 0.1230φ(6d<sup>2</sup><sub>5/2</sub>)</td></tr><tr><td>105</td><td>Db</td><td>ψ<sub>J=3/2</sub> = 0.7869φ(6d<sup>3</sup><sub>3/2</sub>) - 0.5083φ(6d<sup>2</sup><sub>3/2</sub>d<sup>1</sup><sub>5/2</sub>) -<br>0.2510φ(6d<sup>1</sup><sub>3/2</sub>d<sup>2</sup><sub>5/2,0</sub>) + 0.2327φ<sub>2</sub>(6d<sup>1</sup><sub>3/2</sub>d<sup>2</sup><sub>5/2,2</sub>) - 0.0727φ(6d<sup>3</sup><sub>5/2</sub>)</td></tr><tr><td>110</td><td>Ds</td><td>ψ<sub>J=4</sub> = 0.9775φ(6d<sup>4</sup><sub>5/2</sub>) + 0.2110φ(6d<sup>3</sup><sub>3/2</sub>d<sup>5</sup><sub>5/2</sub>)</td></tr><tr><td>114</td><td>Fl</td><td>ψ<sub>J=0</sub> = 0.9955φ(7p<sup>2</sup><sub>1/2</sub>) - 0.095φ(7p<sup>2</sup><sub>3/2</sub>)</td></tr><tr><td>115</td><td>Mc</td><td>ψ<sub>J=3/2</sub> = 0.9932φ(7p<sup>2</sup><sub>1/2</sub>p<sup>1</sup><sub>3/2</sub>) + 0.099φ(7p<sup>1</sup><sub>1/2</sub>p<sup>2</sup><sub>3/2</sub>) -<br>0.061φ(7p<sup>3</sup><sub>3/2</sub>)</td></tr><tr><td>116</td><td>Lv</td><td>ψ<sub>J=0</sub> = 0.9984φ(7p<sup>2</sup><sub>1/2</sub>p<sup>2</sup><sub>3/2</sub>) + 0.057φ(7p<sup>1</sup><sub>3/2</sub>p<sup>3</sup><sub>3/2</sub>)</td></tr></tbody></table>

of the 7p block elements can be unambiguously assigned to a single dominant configuration.

## 6.2 Ionization potentials

Ionization potentials ($I_p$) and electron affinities ($E_A$) are important for discussing the chemical behaviour of an element. For example, Mulliken's (empirical) defini-
<figure><img src="image_22.png" /><figcaption>**Fig. 26** Ionization potentials (in eV) for the s-, d-, and p-block elements along period 4-7 and the four first elements of period 8. Experimental values are from NIST [419] and from various authors (see text). For the transactinides, see Table 2. The Period 7 solid lines correspond to ionization potentials obtained by a linear fit, which are explained in the main text.</figcaption></figure>

tion of the (dimensionless) electronegativity $\chi$ [506] of an element takes the sum of both properties in units of eV, $\chi = (0.187/eV)(I_p + E_A) + 0.17$.

The predicted ionization potentials for the transactinides are summarized in Fig. 26. Only for the elements up to the actinides ionization potentials have been determined experimentally. The experimental value of $4.96_{-0.04}^{+0.05}$ eV for Lr is in good agreement with the value of 4.963(15) eV obtained from relativistic coupled cluster calculations, CCSD(T), which include Breit and QED contributions [19, 20].

Concerning the accuracy of the data listed in Table 2, the decrease in the ionization potential of Ds compared to Mt and Rg is a reminder for the limited accuracy of the electron correlation treatment of open-shell systems, as this decrease is most likely due to the fact that the CIPT (configuration interaction with perturbation theory) method of Refs. [137, 489] overestimates ionization potentials of the last 6d-block elements (they give error bars of about 1 eV). Their CIPT values show a steady increase in the ionization potential from Mt (10.3 eV) to Ds (11.2 eV) to Rg (12.2 eV) and Cn (13.1 eV). We therefore took their lowest listed values (fitting parameter $a=0$ in Table III of Ref. [489]) to show a more likely trend for the transactinides in Fig. 26, i.e., 5.91 eV (Db), 6.83 eV (Sg), 7.70 eV (Bh), 8.57 eV (Hs), 9.43 eV (Mt), and 10.3 eV (Ds). The most accurate values for the ionization potentials, electron affinities and polarizabilities (including excited states to confirm the correct ground state symmetry) of a transactinide element come from Fock-space coupled cluster theory [436, 439, 440, 442, 443].
Figure 26 shows a smooth increase in the ionization potentials for the 6d (Periodic 7) elements due to an increase in nuclear charge causing the 6d shell to become more compact. This is consistent with the lighter d-block elements. They do, however, start at a lower ionization potential compared to the lighter systems. This can be explained by the more diffuse and relativistically-expanded 6d orbitals. However, starting at Bh, where occupation of the relativistically-destabilized $5d_{5/2}$ shell begins, we observe a higher ionization potential compared to the lower d-block elements in the PT, most likely due to a less sufficient screening of the nucleus by the 6d orbitals compared to the 5d orbitals. A drop in the ionization potential from the Group 12 to the Group 13 is expected as the underlying ns-shell and fully occupied ($n-1$d-shell are more compact compared to the np-shell. Moreover, the s-shell undergoes a strong relativistic stabilization.

## 6.3 QED effects

Calculated vacuum polarization $\Delta E_{n\ell j}^{\text{VP}}$ and self-energy $\Delta E_{n\ell j}^{\text{SE}}$ (absolute value) contributions to the total electronic energy for element Ubn ($Z = 120$) are shown in Fig. 27. As both VP and SE operators act in the close vicinity of the nucleus [31], the major QED contributions come from the $1s_{1/2}$ and $p_{1/2}$ shells. We also observe a decrease in QED contributions with increasing principal quantum number *n* for fixed ($\ell j$), a decrease with increasing angular quantum number $\ell$ for fixed *n* and increasing total quantum number *j* for fixed $n\ell$. This can all be explained by the changing density around the nucleus with changing combinations of ($n\ell j$).

As we are interested in valence shell properties, one needs to know if QED-related changes in shell occupations stem from the valence shell or from relaxation effects of deeper-lying core shells. As shown in Eq. (53), QED effects of the valence shell vary as $1/n^3$. Figures 21 and 22 show the strong decrease of $F_{(n,\kappa)}^{(1)}(Z\alpha)$ with increasing $|\kappa|$ values for calculations with finite-size correction. This is illustrated in Fig. 27 for Ubn ($Z = 120$). We therefore list in Table 4 the QED per-shell D-HF contributions to the ionization potentials for some selected atoms, Cn, Nh, and Ubn ($Z = 112, 113$ and 120), as well as the QED contributions to the electron affinities of Nh and Uue ($Z = 119$).

Regarding the valence shell ionization potential, for Cn most of the QED effect comes from the relaxation of the $7s$ shell. The large contribution from $6p_{3/2}$ may come as a surprise, whereas for the $6p_{1/2}$ shell we have cancellation effects between the VP and SE contributions. For the Nh ionization, we see a slightly different picture with a small QED contribution due to an almost perfect cancellation of the $7p_{1/2}$ and $7s^2$ shell contributions. For the Ubn ionization potential, we obtain the largest contribution from the $8s$ shell, but the $7p_{3/2}$ and $7s$ shells also have substantial contributions (our result here deviates from the value given in Ref. [321]). Regarding the individual shell contributions to the electron affinities (defined here as positive values) of Nh and Uuh, we observe for both elements that the main contribution comes from the shell where the electron occupation is altered, however there are also large contributions of shells where the occupations stay the same. These results show that relaxation effects are important. For Cn, the dominant contribution does not come from the shell where the electron occupation is altered as QED effects are much larger for relaxing the *s* shell than removing an electron from the underlying *d*-shell.
<figure><img src="image_23.png" /><figcaption>**Fig. 27** Absolute values for the vacuum polarization $|\Delta E_{n\ell j}^{\text{VP}}|$ (solid lines) and self-energy $|\Delta E_{n\ell j}^{\text{SE}}|$ (dashed lines) for different $(n\ell j)$ shells (per electron) of Ubn ($Z=120$). The formalism of Ref. [67] for the self-energy and the Uehling potential with a finite nuclear charge distribution was used [321]. All shell contributions to the VP have a negative sign, all shell contributions to the SE have a positive sign except for the $4f_{5/2}$ and $5f_{5/2}$ shells.</figcaption></figure>

We note that Og was predicted to be the first rare gas element with a non-zero electron affinity of 0.080 eV [41, 503, 507], where the Breit interaction contributes with $-3 \times 10^{-4}$ eV and QED with $-3 \times 10^{-3}$ eV, the latter already comparable in size to the quadruple contributions in a CC treatment [503].

As shown in Fig. 28, the VP and SE contributions for the valence shell ionization potentials approximately obey a simple power law [321]

$$E(Z) = CZ^{\gamma}. \qquad (69)$$

The exponent $\gamma$ of the VP is roughly 40 to 50% larger than the one of the SE, and we observe a lower scaling for the valence-$p$ shell compared to the $s$-states. An logarithmic-scale extrapolation to high $Z$ shows that the VP and SE curves cross at $Z = 160$ for the Group 11 and $Z = 139$ for the Group 1 elements.

For the superheavy elements beyond nuclear charge $Z = 120$, the accurate treatment of electron correlation methods still remains the major bottleneck. It could therefore be sufficient to include QED and Breit interactions at a perturbative level, especially if higher-$\ell$ states are involved and shell relaxation effects are small. For example, Ref. [14] performed calculations on a series of transactinides with $Z \ge 140$. For the elements with $Z = 171, 172$, and $173$ QED effects contribute to the ionization potential by 1.7%, 0.1%, and 1.2%, respectively (the very low value for $Z = 172$ comes from an almost exact screening of the SE and an exact cancellation of VP in the neutral and singly ionized case [14]).
<table><caption>**Table 4** QED contributions to the ionization potentials of Cn, Nh and Ubn and electron affinities of Nh and Uue obtained from D-HF calculations using GRASP [321]. For Nh the QED contributions are shown in eV rather than in percentage as cancellation effects in the valence shell lead to a very small total QED contribution to the ionization potential.</caption><thead><tr><th>Element</th><th colspan="3">ionization potentials</th><th colspan="2">electron affinities</th></tr><tr><th>Transition</th><th><sup>1</sup>S<sub>0</sub>(6d<sup>10</sup>7s<sup>2</sup>)<br/>→<sup>2</sup>D<sub>5/2</sub>(6d<sup>4</sup><sub>3/2</sub>d<sup>5</sup><sub>5/2</sub>7s<sup>2</sup>)</th><th><sup>2</sup>P<sub>1/2</sub>(7s<sup>2</sup>7p<sup>1</sup><sub>1/2</sub>)<br/>→<sup>1</sup>S<sub>0</sub>(7s<sup>2</sup>)</th><th><sup>1</sup>S<sub>0</sub>(8s<sup>2</sup>)<br/>→<sup>2</sup>S<sub>1/2</sub>(8s<sup>1</sup>)</th><th><sup>2</sup>P<sub>1/2</sub>(7p<sup>1</sup><sub>1/2</sub>)<br/>→<sup>1</sup>S<sub>0</sub>(7p<sup>1</sup><sub>1/2</sub>)</th><th><sup>2</sup>S<sub>1/2</sub>(8s<sup>1</sup>)<br/>→<sup>1</sup>S<sub>0</sub>(8s<sup>2</sup>)</th></tr></thead><tbody><tr><td>QED tot (eV)</td><td>0.0227</td><td>-0.0003</td><td>0.0110</td><td>0.0021</td><td>-0.0020</td></tr><tr><td>QED VP (eV)</td><td>-0.0150</td><td>0.0029</td><td>-0.0060</td><td>0.0005</td><td>0.0031</td></tr><tr><td>QED SE (eV)</td><td>0.0377</td><td>-0.0032</td><td>0.0170</td><td>0.0026</td><td>-0.0051</td></tr><tr><td>Contributions<br/>to total</td><td>7s (79.9%)<br/>6d<sub>5/2</sub> (-17.6%)<br/>6p<sub>3/2</sub> (22.2%)<br/>6s (14.5%)</td><td>7p<sub>1/2</sub> (-0.0109 eV)<br/>7s<sup>2</sup> (0.0104 eV)</td><td>8s (78.1%)<br/>7p<sub>3/2</sub> (35.3%)<br/>7s (-16.7%)<br/>6p<sub>3/2</sub> (4.9%)</td><td>7s (356%)<br/>7p<sub>1/2</sub> (-227%)<br/>6d<sub>5/2</sub> (18.9%)<br/>6d<sub>3/2</sub> (11.9%)<br/>6p<sub>3/2</sub> (-20.6%)<br/>6p<sub>1/2</sub> (-29.7%)</td><td>8s (123 %)<br/>7p<sub>3/2</sub> (-68.9%)<br/>7s (41.6 %)<br/>6p<sub>3/2</sub> (-7.3%)<br/>6s (6.2 %)</td></tr><tr><td>QED effect</td><td>6d<sub>5/2</sub> (7.9%)<br/>6p<sub>1/2</sub> (2.0%)<br/>1s (-0.7%)</td><td></td><td>6s (-3.8%)</td><td></td><td></td></tr></tbody></table>

The accurate treatment of QED effects is more important when it comes to the prediction of inner shell ionization potentials [508, 509], especially for high-Z few electron systems [116], or for K-capture rates for neutron deficient nuclei [510]. For the latter, little theoretical work in this direction has been done so far [511].

## 6.4 Atomic static dipole polarizabilities

Atomic static dipole polarizabilities $\alpha_D$ are very useful quantities for chemical reactions and are, for example, used in the simulation of temperature dependent atom-at-a-time experiments of transactinides absorbed on surfaces such as gold or quartz [493, 495, 498, 512-517]. Dipole polarizabilities are proportional to the inverse of the ionization potential, $\alpha_D \sim I_p^{-1}$ [518], which can be argued from the sum-over-states formula. Empirically one approximately finds $\alpha_D = 6.67 I_p^{-2}$ a.u. (with RMSD of 0.93 a.u.).

Calculated dipole polarizabilities are listed in Table 2. They have recently been reviewed for the known elements in the PT [519], where periodic trends were also discussed. Note that we only discuss here the scalar component of the polarizability tensor as shown in Table 2. If strong spin-orbit coupling is involved, one requires knowledge of the $M_J$-resolved components for the Stark effect in open-shell atoms [520].

The start of a new shell occupation usually increases the polarizability; this is seen for Lr, Nh, Mc and element 119. The rather large DHF+CI+Breit+QED polarizability for Lr comes with large estimated uncertainties [487, 519]. Nevertheless, we see a decreasing trend in dipole polarizabilities across the transactinide series, and perhaps the spin-orbit splitting is not large enough to explain an increase in $\alpha_D$ from Sg to Bh upon occupation or the $6d_{5/2}$ shell.
<figure><img src="image_24.png" /><figcaption>Fig. 28 (a) Self-energy (SE) and (b) vacuum polarization (VP) energy contributions to the ionization potential against the atomic number Z for different atoms within a Group of the Periodic Table. Data are taken from Ref. [321]</figcaption></figure>

## 6.5 Electron Localization Function

The concept of the electron localization function (ELF) was introduced in Ref. [521], and later applied to atoms, molecules and the solid state [522, 523]. More recently, the concept was extended to nucleon localization functions (NLF) in nuclear structure theory [524–526]. The ELF provides a measure of finding an electron in the vicinity of another same-spin electron located at a given position $r$:

$$D_{\sigma}(\mathbf{r}) = \left[ 1 + \left( \frac{\tau_{\sigma}(\mathbf{r})\rho_{\sigma}(\mathbf{r}) - \frac{1}{4}|\nabla\rho_{\sigma}(\mathbf{r})|^2}{\rho_{\sigma}(\mathbf{r})\tau_{\sigma}^{\text{TF}}(\mathbf{r})} \right)^2 \right]^{-1}, \quad (70)$$

where spin $\sigma$ is ↑ or ↓, $\rho_{\sigma}$ is the electron spin density, $\tau_{\sigma}$ is the kinetic energy density, $\nabla\rho_{\sigma}$ is the electron density gradient, and $\tau_{\sigma}^{\text{TF}}$ is the Thomas-Fermi kinetic energy (of the uniform electron gas). The ELF can be seen as a tool to identify regions where electrons are localized or delocalized (for a review see for example Ref. [527]).
<figure><img src="image_25.png" /><figcaption>**Fig. 29** ELF's $D_σ(r)$ from nonrelativistic (NR, top), scalar relativistic (SR, middle) and fully relativistic (R, bottom) DFT calculations for Hg and Cn (left), Pb and Fl (center), and Rn and Og (right). Adopted from Ref. [528].</figcaption></figure>

The ELF assumes values between 0 and 1, where a value close to 1 indicates that the probability of finding two same-spin electrons close to each other is very low (high level of electron localization) and the ELF value of 0.5 corresponding to the limit of a hypothetical uniform Fermi gas of the same density (high level of electron delocalization). ELFs are therefore ideal to show the changes in electron localization due to the influence by other atoms or by relativistic effects. The discussion of ELF and NFL in superheavy nuclei can be found in Ref. [502].

To demonstrate the importance of relativistic effects for the superheavy elements we show in Fig. 29 the ELF's for Cn, Fl, and Og in comparison with their lighter congeners. The apparent delocalization is due to scalar relativistic effects in the case of Cn, and due to spin-orbit effects for both Fl and Og. The appreciable relativistic effects are already notable for Hg. It has been argued that the large spin-orbit splitting in the 7p shell of Og with 10.13 eV is responsible for a uniform-gas-like behavior in the valence region [502]; see also discussion in Ref. [500].

## 6.6 Examples of relativistic effects on the chemistry of SHE

Relativistic effects have a profound influence on the chemistry of the transactinides. For example, due to the large relativistic 7s stabilization, Cn is expected to behave like a rare gas with a low predicted melting and boiling point of −18(30) °C and 77(10) °C, respectively (at the nonrelativistic level ca. 370 °C and 970 °C respectively) [529, 530]. In contrast, relativistic effects can cause the Group-18 element Og to be a semi-conductor and a solid at room temperature (RT) with an electron localization function resembling that of a Thomas-Fermi gas [502, 531, 532]. The trends in melting points for the Group 12, 14 and 18 of the periodic table are summarized in Fig. 30. Each of the predicted melting points for the three super-heavy elements are very close to room temperature and are highly influenced by relativistic effects.
<figure><img src="image_26.png" /><figcaption>**Fig. 30** Melting points for the Group 12, 14 and 18 elements of the periodic table. Experimental values are from [533]. For Cn, Fl and Og the melting points are taken from computational simulations [528–530].</figcaption></figure>

As another example, we take the closed $7p_{1/2}^2$ shell of the Fl atom (Period 7, Group 14), for which the relativistic effects are predicted to cause a large $7p_{1/2} - 7p_{3/2}$ spin-orbit splitting. This explains the maximum in the ionization potential in Fig. 26 and its expected low melting temperature [530]. Furthermore, the d-block transactinides are expected to always ionize out of the 6d shell as the fully occupied 7s shell undergoes a very strong relativistic energetic stabilization and contraction. This was pointed out early on for Rg, which adopts the $d^9s^2$ configuration rather than the usual $d^{10}s$ arrangement [441, 534]. For further reading, see Refs. [514, 535–537].

# 7 General Considerations

## 7.1 Placing new elements on the periodic table

The correct placement of the elements in the PT (see Fig. 1) has been a matter of intense debate [2, 13, 15]. Besides the evident ordering of elements according to their atomic number Z, the additional two principles for the placement of atoms into the PT are the Pauli principle, derived from the spin-statistics theorem, and the Aufbau principle, derived from mean-field theory [1, 5]. The Pauli principle dictates that only one electron can be put into a (spin-)orbital, or in the relativistic case into an orbital with quantum numbers (n, l, j, mj). The electrons then fill the orbitals in the order of increasing orbital energy, leading to the Aufbau principle, where the leading configuration for each element can be obtained from either Kohn-Sham or from MCSCF theory [1]. By formatting the PT in two dimension, the configurations of valence electrons are repeated within a row with increased
<figure><img src="image_27.png" /><figcaption>**Fig. 31** Predicted ground state configurations for the elements with atomic numbers Z= 1 Empty orbitals in green, filled orbitals in blue and the partly filled orbitals in yellow. 19-172 [16]. This diagram visualizes the difficulty with placing the SHE elements in the PT. Whereas some elements have a definite place, a handful of elements, such as Z = 121 – 125, cannot be uniquely placed. The elements Z = 133 and Z = 134 share the same place, just as Z = 148, 149 and 150. Also the 5g and 6f elements contain higher-order occupations of 6f and 7d state.</figcaption></figure>

principal quantum numbers. Since, in the nonrelativistic approach, the valence electrons dictate the angular distribution of the wave functions, there is a periodicity of the chemical properties [1, 538, 539]. Relativistic effects, on the other hand, can substantially alter the behavior of an element. This is especially observed for the main group and late transition elements [1, 3, 514, 540–543]. These relativistic effects as well as electron correlation contributions lead to exceptions to the Aufbau principle, which cause difficulties with the element placement. Electron configurations can alter within a group of the PT as well, and the Group 10 elements serve as a perfect example. Here one has the dominant electron configuration 3d⁸4s² for Ni, 4d¹⁰ for Pd, 5d⁹6s¹ for Pt, and 6d⁸7s² for Ds. Yet, inspection
of the excited states reveals that apart from the principal quantum number, all three configurations $(n-1)d^8ns^2$, $(n-1)d^9ns^1$ and $(n-1)d^{10}$ are close in energy [1]. Another difficulty is the correct placement (starting and ending point) of the *f*-block elements [544, 545], see Fig. 31. This makes the PT a somehow fuzzy concept in the superheavy region.

## 7.2 Dominant ground state configurations predicted by electronic structure calculations

A number of authors have tried to predict the dominant electron configuration(s) beyond the element with nuclear charge $Z = 118$ [11, 12, 16, 17, 534, 535, 546, 547]. Early attempts were made by using D-HF-Slater calculations for the range $Z = 118-131$ [11, 548, 549]. However, these mean-field studies did not explicitly include static and dynamic electron correlations. Since the total energies of different configurations differ little from one another in many cases, inclusion of configuration interaction can lead to a change of the ground-state symmetry and configuration. To improve the early results, multi-configuration Dirac-Fock calculations have been carried out [16] with inclusion of a total angular momentum coupling scheme (see Sec. 3) and the Breit correction. These results were extended to chemically plausible ions [17], to differentiate between free ions and ions in chemical compounds. The resulting PT is shown in Fig. 1. In this work, the starting point of the 5g elements has been placed at $Z = 121$ (in contrast to Ref. [16] that places the starting point at $Z =126$). Because of the high density of states in the region beyond $Z = 120$, it becomes quite difficult to pick the dominant configuration and ground state symmetry. For example, Ref. [14] showed from average level (AL) calculations that for the element with $Z = 140$ the ground state configuration is $8s^28p^27d6f^{35}g^{14}$ in agreement with Ref. [16], but this could change if the state of specific symmetry is optimized and more configurations are included. This shows that more work is required to determine the ground state configurations and symmetries, as well as associated chemical behaviour [17]. Possible candidates for configurations for the elements up to $Z = 172$, taken from Ref. [16], are visualized in Fig. 31. Needless to say that for the superheavy elements relativistic effects are crucial and need to be correctly accounted for as they can lead to a noticeable violation of simple regularities such as for the Group 11 elements where Rg adopts a configuration with a hole in the $6d_{5/2}$ shell instead of the 7s shell [441].

The third pillar (beside the atomic number and electron configuration) for building the PT comes from chemical properties, similar as how the plausible chemical ions were taken into consideration for the placement of atoms in Ref. [17]. We have learned in the past few decades that relativistic effects can alter chemical properties substantially [540, 541], especially in the superheavy element region [1, 18, 514]. Regarding the chemical properties of the SHE beyond $Z = 120$, there is no information available except from the computational studies. For example, it has been shown [550] that the 5g-electrons for the elements with atomic numbers $Z = 125-129$ are core-like and only act as spectators, similar to the 4f electrons in the lanthanides. It is clear that more work needs to be done to study the chemical properties in the $Z > 120$ region, which will help designing future atom-at-a-time chemistry experiments [551].
## 7.3 Periodic Table - How far can we go?

Oganesson was the heaviest element and nihonium the last element to be added officially into the PT [552]. Thus the 7th period of elements is now complete. The question arises if one can go much further in the atomic number. From an electronic point of view, there is no limitation to Z. While the correct description of multi-electron systems with $Z > Z_c \approx 170$ is still a difficult problem, and the inclusion of Gamow states for multi-electron systems needs to be addressed, the real limitation to the PT comes from the nuclear stability [18, 553]. In the transactinide region, $Z > 103$, in early days known as the sea of instability [554], the half-life of the elements varies between hours ($^{266}_{103}\text{Lr}$) and seconds ($^{294}_{118}\text{Og}$) or below. Although there is a small predicted region of increased stability between Z=114-126 and N = 184 with predicted lifetimes of hours or even days [18, 555-557], it is currently not clear how far the PT can be extended from the nuclear point of view. Indeed, the IUPAC defines an element to exist if its lifetime is longer than $T_{el} \approx 1 \times 10^{-14}$ s, which is the time it takes for electron cloud to form around the nucleus. This means that for atomic nuclei living shorter than $T_{el}$ it makes no sense to talk about atoms and chemistry [18, 553].

The lifetimes of most known superheavy nuclei are governed by the competition between α-decay and spontaneous fission (SF). The corresponding lifetimes predicted by a particular DFT model [558] are shown in Fig. 32. For a survey of various predictions of α-decay and SF lifetimes, see Refs. [261] and [559], respectively. The shortest SF half-lives, reaching down to $10^{-10}$ s, are predicted for nuclei from a narrow corridor formed by $^{280}\text{Hs}$, $^{284}\text{Fl}$, and $^{284}\text{Og}$. This corridor of fission instability separates the regions of superheavy nuclei synthesized in hot- and cold-fusion reactions. Moving towards more neutron-rich nuclei beyond N = 184, dramatic decrease of SF lifetimes, below $10^{-15}$ s is expected, see Fig. 32(b) and Ref. [560].

It is to be noted that predictions of nuclear models in the region of superheavy nuclei are very sensitive to both input (forces, functionals) and theoretical framework used. Consequently, theoretical lifetime estimates, especially for SF, often differ by many orders of magnitude [559]. The heaviest nuclei synthesized so far are all proton-rich; hence, they can in principle decay by means of electron capture or β$^+$/EC process. So far, no such decay modes have been observed in the upper superheavy region, and this indicates that they cannot compete with the dominant α-decay and SF modes. Indeed, according to theory β$^+$/EC lifetimes shorter than 1 s are expected in nuclei that lie rather far from the current superheavy region [261].

Experiments to synthesize new superheavy nuclei and elements beyond Og are underway [561, 562]. If discovered, these systems will be crucial for testing many-body nuclear structure theories [18]. To explore their chemistry will be very challenging, however.

# 8 Conclusions

Atomic structure theory developed enormously over the past decades to the extend where QED can be tested to high precision for few-electron systems. However, for
<figure><img src="image_28.png" /><figcaption>**Fig. 32** Summary of theoretical predictions [558] for decay modes of superheavy nuclei obtained with nuclear DFT: (a) SF half-lives; (b) α-decay half-lives.</figcaption></figure>

many-electron systems the accurate description of both QED and electron correlation effects remains a major challenge [479]. But even here progress has been made for elements with large atomic numbers [14]. While the negative energy continuum is required for QED, it creates a formidable conceptual and computational challenge, especially when bound states approach the negative energy continuum threshold. The correct description of diving (Gamow) states within a multi-electron formalism including QED effects still needs to be explored. It is clear that effects associated with the negative energy continuum distinguishes the Dirac from the Schrödinger equation in both mathematical and in physical terms [29]. Once these problems are solved, there is no limitation to the treatment of atoms beyond the critical nuclear charge. The PT, seen as the foundation for chemistry, is therefore *not* limited to a certain nuclear charge region, but limited by nuclear stability [18]. The future looks bright for superheavy element synthesis and associated chemistry experiments, which require the support of accurate electronic and nuclear structure theory.
<table><caption>**Table 5** Norm and expectation value existences for different ranges of nuclear charges Z and parameter $\pm\gamma(Z) = \pm\sqrt{1 - (Z\alpha)^2}$ appearing in the radial 1s function of hydrogenic atoms with potential $V(r) = -Z/r$. The symbol S stands for the Sobolev norm ($S = 1$ for Dirac and 2 for Schrödinger) which requires the gradient norm to exist for the Dirac-Coulomb operator (first and second derivatives for the Schrödinger case).</caption><thead><tr><th>System</th><th>range</th><th>||<i>&phi;</i>||<sub>2</sub></th><th>&lang;<i>&phi;</i>|<i>H</i>|<i>&phi;</i>&rang;</th><th>||<i>&phi;</i>||<sub>2</sub><sup>(S)</sup></th><th>||<i>H</i><i>&phi;</i>||<sub>2</sub></th></tr></thead><tbody><tr><td>Schrödinger</td><td>Z &gt; 0</td><td>yes</td><td>yes</td><td>yes</td><td>yes</td></tr><tr><td>Dirac +&gamma;(Z)</td><td>0 &lt; Z&alpha; &lt; &radic;3/2</td><td>yes</td><td>yes</td><td>yes</td><td>yes</td></tr><tr><td>Dirac +&gamma;(Z)</td><td>&radic;3/2 &le; Z&alpha; &lt; 1</td><td>yes</td><td>yes</td><td>no</td><td>no</td></tr><tr><td>Dirac -&gamma;(Z)</td><td>0 &lt; Z&alpha; &lt; &radic;3/2</td><td>no</td><td>no</td><td>no</td><td>no</td></tr><tr><td>Dirac -&gamma;(Z)</td><td>&radic;3/2 &le; Z&alpha; &lt; 1</td><td>yes</td><td>no</td><td>no</td><td>no</td></tr></tbody></table>

# 9 Appendix A: The Self-Adjointness of the Dirac-Coulomb Hamiltonian

In the two appendices we address some of the more mathematical features of the Dirac equation, the self-adjointness problem and (in the next section) the rigged Hilbert space formalism for Gamow states. Both aspects often lead to some misunderstandings in the community and are therefore discussed briefly here.

To explain the non-self-adjointness at critical charge in correct mathematical terms, we require the $L_2$-norm of the derivative $||\frac{d}{dr}\phi||_2$ to exist (or the gradient norm for the three-dimensional case) for the eigensolutions as the Dirac equation is a first-order differential equation (see discussion of Sobolev spaces further below). The radial solutions for a point charge nucleus $\phi(r)$ have the general form $\phi(r) = a_{n\kappa}(2Zr)^{\gamma}e^{-Zr}f_{n\kappa,Z}^{P,Q}(r)$ with the exponent $\gamma = \pm\sqrt{\kappa^2 - (Z\alpha)^2}$ and $f_{n\kappa,Z}^{P,Q}(r)$ containing expressions of $P_{n\kappa}(r)$ and $Q_{n\kappa}(r)$ in terms of confluent hypergeometric functions [89]. Unlike for the Schrödinger equation, $\gamma$ is a non-integer and derivatives lead to negative r-exponents. As a result, both $||\frac{d}{dr}\phi||_2$ and $||H_D\phi||_2$ become infinite if $\gamma \le \frac{1}{2}$ leading to $Z_{c_1\alpha} = \sqrt{3}/2$ as discussed in Sec. 2.4.1. In contrast, for the nonrelativistic radial Schrödinger equation, all derivative norms exist, i.e., $||\frac{d^n}{dr^n}\phi_{NR}||_2 < \infty$, as only integers appear in the r-exponent. This leads often to misunderstandings as the $L_2$ Hilbert space only requires the norm $||\phi||_2$ to exist, but as soon as we introduce an unbound differential operator such as $H_D$ we have to deal with the domain of such an operator and the existence of certain expectation values and derivative norms. For a more detailed analysis using the Weyl's limit point - limit circle theorem, generalized to the Dirac equation by Weidmann [141], the reader is referred to a recent paper by Gallone [99]. The various norm existences for different nuclear charge regions are summarized in Table 5.

To rephrase the self-adjointness condition in different terms, a self-adjoint extension of the radial Dirac operator should have the following domain [29, 141]: $\text{dom}(H_D) = \{\phi \in L_2(\mathbb{R}_+)^2 | \text{ each component of } \phi \text{ is locally absolutely continuous;}$ $H_D\phi \in L_2(\mathbb{R}_+)^2$; deficiency indices $d\{\phi(r = 0)\} = (0,0)\}$ (see also Hogreve [98]), where $L_2(\mathbb{F})^n \equiv L_2(\mathbb{F}) \otimes \mathbb{C}^n$ over a field $\mathbb{F}$ ($\mathbb{F} \equiv \mathbb{R}_+$ and $n = 2$ for the radial Dirac equation) [29]. This allows us to select the Sobolev space $W_{1,2}(\mathbb{R}_+)^2$ [29, 46] as the natural domain for the (unbound) Dirac operator (or similarly for a four-component wave function in the three-dimensional case $W_{1,2}(\mathbb{R}^3)^2$) lying dense in $L_2(\mathbb{R}_+)^2$, such that $\text{dom}(H_D) \subseteq W_{1,2}(\mathbb{R}^3)^2 \subseteq L_2(\mathbb{R}^3)^2$. The domain problem of the Dirac-Coulomb operator has been very recently discussed and critically analyzed
by Estaban [563]. In the subcritical nuclear charge region ($\sqrt{3}/2 \le Z\alpha \le 1$), $\phi$ is an eigenfunction to a non-self-adjoint Dirac-Coulomb Hamiltonian with real eigenvalues and norm $||\phi||_2 < \infty$, but does not belong to $\text{dom}(H_D)$ ($H_D$ self-adjoint)! More generally, one looks for the largest subdomain of the Hilbert space that remains invariant under the action of certain powers of required operators (observables) including the Hamiltonian of the system, which is known as the maximal invariant subspace of the algebra generated by these operators [564].

A note of caution should be added here. If we eliminate the small component and focus on the resulting second-order differential equation, the underlying Sobolev space is now $W_{2,2}(\mathbb{R}^3)^2$, which makes the conditions more stringent for the norm existence. In any case, we seek for an appropriate self-adjoint extension of $H_D$ [48, 101, 565] as physics does not restrict atoms to a maximum critical charge (except for nuclear instability which is an entirely different matter [39]). The necessary boundary condition for securing the self-adjointness of the Dirac operator at the origin has been discussed in detail by Kuleshov [151] and Gitman [48, 566]. The physical realization of this boundary condition is the introduction of a finite nuclear charge density [21].

Last we mention that mathematically, there are many self-adjoint extensions which can be realized for the Dirac operator. For example, Kato showed that a potential energy matrix of the form $V_{ik} = a/2r + b$ with $a < 1$ and $b > 0$ makes the Dirac operator essentially self-adjoint on certain domains [567]. For a more recent discussion on possible self-adjoint extensions we refer to [565].

# 10 Appendix B: The Rigged Hilbert Space Formalism

If we consider the spectrum of the Dirac-Coulomb operator, $\sigma(H_D)$, we need to include the discrete (d) and both the positive (+) and negative (−) continuum states (c) (cf. Figure 2), i.e., $\sigma(H_D) = \{\phi^d\} \cup \{\phi_+^c\} \cup \{\phi_-^c\}$. The continuum states are important in scattering (resonance) theory and are essential for the quantum mechanical completeness relation. It is well known, however, that unlike the discrete states, the continuum states lead to domain problems for unbound operators, i.e., they are not normalizable and therefore do not belong to the quantum mechanically relevant $L_2$ space. As a result, von Neumann’s original Hilbert space formalism requires an extension to include such (generalized) functions. Continuum states are properly defined within an extended or rigged Hilbert space ($\mathcal{G}$) formalism [568–574] originally proposed for the quantum mechanical framework by Roberts, Antoine and Bohm [575–579].‡‡‡ In fact, rigged Hilbert spaces are the structures required for both the discrete orthonormal and continuous bases to coexist [580]. Their existence for self-adjoint operators on separable Hilbert spaces is guaranteed by the Gelfand-Maurin nuclear spectral theorem.

In strict mathematical terms, the rigged Hilbert space (RHS) $\mathcal{G}$ is defined as a triple of topological vector spaces $\mathcal{G} = (\Phi, \mathcal{H}, \Phi^\times)$, called the Gelfand triple [581], generated by an infinite dimensional separable Hilbert space $\mathcal{H}$ such that the denseness relation is $\Phi \subseteq \mathcal{H} \subseteq \Phi^\times$. Here, $\Phi$ is a (complete) nuclear Fréchet space, also called test-function space (not necessarily a Hilbert space, which enables to

<footnote>‡‡‡The term rigged Hilbert space is misleading as $\mathcal{G}$ is not a Hilbert space per se, but is generated from a Hilbert space $\mathcal{H}$.</footnote>
use the nuclear spectral theorem of Gelfand and Maurin [568, 581]), and $\Phi^\times$ is the topological dual (or topological conjugate) of $\Phi$, that is the complete space of continuous anti-linear functionals on $\Phi$ (also called distribution or Schwartz space). The RHS structure includes an inductive limit of a sequence of topological spaces $\Phi^{(n)}$ in which the topologies get rapidly coarser with increasing $n$ [574, 582], e.g., we might think of a series of Sobolev spaces $W_{k+1,2} \subseteq W_{k,2}$, with $W_{1,2} \subseteq L_2$. It is clear from this example that these subspaces have different norms (or semi-norms for the more general nuclear spaces [574]).

The RHS formalism provides a correct mathematical foundation to Dirac's original bra and ket notation [574, 583], used extensively in quantum theory. Needless to say that the Dirac delta "function" is a distribution required for the properties of continuum states belonging to $\Phi^\times$ rather than $\Phi$. To cite Bohm [569]: *The difference between [the rigged] Hilbert space formulation and the usual [von Neumann] Hilbert space formulation appears to be minor from a physicists point of view, but is essential from a mathematical point of view and leads far to tremendous mathematical simplification; in fact it justifies the mathematically undefined operations that the physicists have been accustomed to in their calculations.* The RHS formalism can easily be generalized to a rigged Fock space formalism required in quantum field theory [584–586].§§§

To be more specific, the Gelfand triple for the spectrum of the Dirac operator is chosen as $\Phi \subseteq \text{dom}(H_D) \subseteq L_2 \subseteq \Phi^\times$. Vectors in $\Phi$ will be complex linear, and the vectors in $\Phi^\times$ are complex antilinear continuous functionals compatible with the scalar product in the underlying Hilbert space, e.g., $F: \psi \in \Phi \to \mathbb{C}$. For example, we define the action $F \in \Phi^\times$ on $\psi \in \Phi$ as an extension to the Hilbert space product $F(\psi) = \langle\psi|F\rangle = \langle F|\psi\rangle^*$. This action is linear to the right and antilinear to the left. Continuity is defined such that if $\psi_n \to \psi$ for $n \to \infty$, $\psi_n \in \Phi$, $\psi \in L_2$ then $F(\psi_n) \to F(\psi)$ in $\mathbb{C}$.

States with given quantum numbers $n, \kappa$, when diving into the negative energy continuum, have complex eigenenergies [151, 204] and therefore move out of the natural domain of the self-adjoint operator $H_D$. We need to give such states belonging to a subset of distributions in the space $\Phi^\times$ a physical interpretation, however, we need first to interpret such generalized eigenfunctions (or eigenfunctions) from a mathematical point of view.

Let $A : \Phi \to \Phi$ a (closed) linear operator and $A^\times : \Phi^\times \to \Phi^\times$ its natural extension of the usual adjoint operator $A^\dagger$ such that $F(A\psi) = A^\times F(\psi) = \langle A\psi|F \rangle = \langle \psi|A^\times F \rangle$ for all $\psi \in \Phi$, $F \in \Phi^\times$. According to the Riesz representation theorem, for every $F \in \Phi^\times$ and linear operator $A$ there exist a unique complex function $\phi$ with complex eigenvalue $\lambda$, $A\phi = \lambda\phi$ such that for all $\psi \in \Phi$ we have $F(A\psi) = \langle A\psi|F \rangle = \langle \psi|A^\times\phi \rangle = \lambda^*\langle \psi|\phi \rangle$ [569]. We call $\phi$ the generalized eigenfunction of $A$ with respect to $F$. These extended eigenstates can be used in the normal way using Dirac's notation keeping in mind that these may not be normalizable and, in general, have complex eigenvalues. Of course, the scalar product $\langle\psi|\phi\rangle$ always needs to be finite which may require a specific (semi)norm definition or regularization of integrals.

<footnote>§§§As quantum operators (in first quantization) act on Hilbert spaces, quantum field operators act on Fock spaces.</footnote>

<footnote>There is always another rigged Hilbert space $\Phi \subseteq \mathcal{H} \subseteq \Phi^*$, where $\Phi^*$ is the dual space of $\Phi$ containing the continuous, linear functionals over $\Phi$ [575, 587]. Dirac's bras and kets belong to $\Phi^*$ and $\Phi^\times$, respectively, and both spaces are isometrically isomorph.</footnote>
**Acknowledgements** We acknowledge financial support by the Program Hubert Curien Du-
mont d'Urville New Zealand - France Science & Technology Support Program number 43245QC,
and the Marsden Fund of the Royal Society of New Zealand. We thank Profs. Trond Saue
(Toulouse), Ephraim Eliav (Tel Aviv), W. H. Eugen Schwarz (Siegen), James Avery (Copen-
hagen) and Vladimir Shabaev (St. Petersburg) for many stimulating and critical discussions.
This material is also based upon work supported by the U.S. Department of Energy, Office of
Science, Office of Nuclear Physics under award number DE-SC0013365. P.I. is a member of
the Allianz Program of the Helmholtz Association, contract no EMMI HA-216 "Extremes of
Density and Temperature: Cosmic Matter in the Laboratory".

**References**

1. P. Schwerdtfeger, O. R. Smits, and P. Pyykkö, Nature Rev. Chem. **4**, 359 (2020), URL https://doi.org/10.1038/s41570-020-0195-y.

2. E. R. Scerri, in *Philosophy of Chemistry*, edited by A. I. Woody, R. F. Hendry, and P. Needham (North-Holland, Amsterdam, 2012), vol. 6 of *Handbook of the Philosophy of Science*, pp. 329–338, URL http://www.sciencedirect.com/science/article/pii/B9780444516756500244.

3. P. Pyykkö, Chemical Reviews **112**, 371 (2012), URL https://doi.org/10.1021/cr200042e.

4. C. Cao, R. E. Vernon, W. H. E. Schwarz, and J. Li, Frontiers in Chemistry **8**, 813 (2021), ISSN 2296-2646, URL https://www.frontiersin.org/article/10.3389/fchem.2020.00813.

5. W. H. E. Schwarz, U. Müller, and F. Kraus, Zeitschrift für anorganische und allgemeine Chemie **648**, e202200008 (2022), https://onlinelibrary.wiley.com/doi/pdf/10.1002/zaac.202200008, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/zaac.202200008.

6. G. Restrepo, Chemistry - A European Journal **25**, 15430 (2019), https://onlinelibrary.wiley.com/doi/pdf/10.1002/chem.201902802, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/chem.201902802.

7. S. Shaik, E. Cremades, and S. Alvarez, Angew. Chem. Int. Ed. **58**, 13194 (2019), https://onlinelibrary.wiley.com/doi/pdf/10.1002/anie.201904584, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/anie.201904584.

8. G. Steinhauser, Nachr. Chemie **67**, 8 (2019), https://onlinelibrary.wiley.com/doi/pdf/10.1002/nadc.20194086536, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nadc.20194086536.

9. Y. Oganessian, J. Phys. G: Nucl. Part. Phys. **34**, R165 (2007), URL https://doi.org/10.1088/2F0954-3899/2F34/2F4/2Fr01.

10. Y. T. Oganessian, Radiochim. Acta **99**, 429 (2011), URL https://doi.org/10.1524/ract.2011.1860.

11. B. Fricke, W. Greiner, and J. T. Waber, Theoret. Chim. Acta **21**, 235 (1971), URL 10.1007/BF01172015.

12. B. Fricke and J. McMinn, Naturwiss. **63**, 162 (1976).

13. E. Scerri, Scientific American **308**, 68 (2013).
14. P. Indelicato, J. Bieroń, and P. Jönsson, Theor. Chem. Acc. **129**, 495 (2011), URL https://doi.org/10.1007/s00214-010-0887-3.

15. P. Pyykkö, Pure and Applied Chemistry **91**, 1959 (2019), URL https://doi.org/10.1515/pac-2019-0801.

16. V. I. Nefedov, M. B. Trzhaskovskaya, and V. G. Yarzhemskii, Doklady Physical Chemistry **408**, 149 (2006), ISSN 1608-3121, URL https://doi.org/10.1134/S0012501606060029.

17. P. Pyykkö, Physical Chemistry Chemical Physics **13**, 161 (2011), URL https://doi.org/10.1039/C0CP01575J.

18. S. A. Giuliani, Z. Matheson, W. Nazarewicz, E. Olsen, P.-G. Reinhard, J. Sadhukhan, B. Schuetrumpf, N. Schunck, and P. Schwerdtfeger, Rev. Mod. Phys. **91**, 011001 (2019), URL https://link.aps.org/doi/10.1103/RevModPhys.91.011001.

19. T. K. Sato, M. Asai, A. Borschevsky, T. Stora, N. Sato, Y. Kaneya, K. Tsukada, C. E. Düllmann, K. Eberhardt, E. Eliav, et al., Nature **520**, 209 (2015), URL https://doi.org/10.1038/nature14342.

20. T. K. Sato, M. Asai, A. Borschevsky, R. Beerwerth, Y. Kaneya, H. Makii, A. Mitsukai, Y. Nagame, A. Osa, A. Toyoshima, et al., J. Am. Chem. Soc. **140**, 14609 (2018), pMID: 30358998, https://doi.org/10.1021/jacs.8b09068, URL https://doi.org/10.1021/jacs.8b09068.

21. I. Y. Pomeranchuk and Y. A. Smorodinsky, J. Fiz. USSR **9**, 97 (1945).

22. P.-G. Reinhard, W. Greiner, and H. Arenhövel, Nucl. Phys. A **166**, 173 (1971), ISSN 0375-9474, URL http://www.sciencedirect.com/science/article/pii/0375947471904210.

23. Y. B. Zeldovich and V. S. Popov, Soviet Physics Uspekhi **14**, 673 (1972), URL https://doi.org/10.1070/pu1972v014n06abeh004735.

24. B. Müller, H. Peitz, J. Rafelski, and W. Greiner, Phys. Rev. Lett. **28**, 1235 (1972), URL https://link.aps.org/doi/10.1103/PhysRevLett.28.1235.

25. V. S. Popov and V. D. Mur, Sov. J. Nucl. Phys. **18** (1974).

26. J. Reinhardt and W. Greiner, Rep. Prog. Phys. **40**, 219 (1977), URL https://doi.org/10.1088/2F0034-4885/2F40/2F3/2F001.

27. J. Reinhardt, B. Müller, and W. Greiner, Phys. Rev. A **24**, 103 (1981), URL https://link.aps.org/doi/10.1103/PhysRevA.24.103.

28. W. Greiner, B. Müller, and J. Rafelski, *Quantum Electrodynamics of Strong Fields* (Springer, Berlin, 1985).

29. B. Thaller, *The Dirac equation* (Springer, Berlin, 1992), URL https://doi.org/10.1007/978-3-662-02753-0.

30. D. M. Gitman, A. D. Levin, I. V. Tyutin, and B. L. Voronov, Phys. Scr. **87**, 038104 (2013), URL https://doi.org/10.1088/0031-8949/87/03/038104.

31. P. Schwerdtfeger, L. F. Pašteka, A. Punnett, and P. O. Bowman, Nucl. Phys. A **944**, 551 (2015), ISSN 0375-9474, URL 10.1016/j.nuclphysa.2015.02.005.

32. V. Shabaev, A. Bondarev, D. Glazov, Y. Kozhedub, I. Maltsev, A. Malyshev, R. Popov, D. Tumakov, and I. Tupitsyn, arXiv preprint arXiv:1910.01373 (2019).

33. G. Gamow, Z. Phys. **51**, 204 (1928), URL https://link.springer.com/article/10.1007/BF01343196.

34. G. Gamow, Z. Phys. **52**, 510 (1929), URL https://doi.org/10.1007/BF01339451.
35. A. J. F. Siegert, Phys. Rev. **56**, 750 (1939), URL https://link.aps.org/doi/10.1103/PhysRev.56.750.

36. A. Bohm, M. Gadella, and J. Dollard, *Dirac Kets, Gamow Vectors and Gel'fand Triplets: The Rigged Hilbert Space Formulation of Quantum Mechanics. Lectures in Mathematical Physics at the University of Texas at Austin* (Springer, berlin, 1989).

37. O. Civitarese and M. Gadella, Physics Reports **396**, 41 (2004), ISSN 0370-1573, URL https://www.sciencedirect.com/science/article/pii/S0370157304001085.

38. W. Nazarewicz, Nature Phys. **14**, 537 (2018), URL https://doi.org/10.1038/s41567-018-0163-3.

39. W. Nazarewicz, Journal of Physics G: Nuclear and Particle Physics **43**, 044002 (2016), URL https://doi.org/10.1088/2F0954-3899/2F43/2F4/2F044002.

40. P. Indelicato, J. Santos, S. Boucard, and J.-P. Desclaux, Eur. Phys. J. D **45**, 155 (2007), URL http://link.springer.com/article/10.1140/2Fepjd/2Fe2007-00229-y.

41. B. G. C. Lackenby, V. A. Dzuba, and V. V. Flambaum, Phys. Rev. A **98**, 042512 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.042512.

42. B. G. C. Lackenby, V. A. Dzuba, and V. V. Flambaum, Phys. Rev. A **101**, 012514 (2020), URL https://link.aps.org/doi/10.1103/PhysRevA.101.012514.

43. E. E. Salpeter and H. A. Bethe, Phys. Rev. **84**, 1232 (1951), URL https://link.aps.org/doi/10.1103/PhysRev.84.1232.

44. W. Greiner et al., *Relativistic quantum mechanics*, vol. 2 (Springer, 2000).

45. I. P. Grant, *Relativistic quantum theory of atoms and molecules: theory and computation*, vol. 40 (Springer Science & Business Media, 2007).

46. R. D. Richtmyer and C. Burdorf, *Principles of advanced mathematical physics*, vol. 1 (Springer, 1978).

47. V. G. Bagrov and D. Gitman, *Exact solutions of relativistic wave equations*, vol. 39 (Springer Science & Business Media, 1990).

48. D. M. Gitman, I. V. Tyutin, and B. L. Voronov, *Self-adjoint extensions in quantum mechanics: general theory and applications to Schrödinger and Dirac equations with singular potentials*, vol. 62 (Springer Science & Business Media, 2012).

49. I. Sargsjan et al., *Sturm-Liouville and Dirac Operators*, vol. 59 (Springer Science & Business Media, 2012).

50. V. G. Bagrov and D. Gitman, *The Dirac equation and its solutions*, vol. 4 (Walter de Gruyter GmbH & Co KG, 2014).

51. G. Magnifico, T. Felser, P. Silvi, and S. Montangero, Nature Communications **12**, 3600 (2021), URL https://doi.org/10.1038/s41467-021-23646-3.

52. F. J. Dyson, Phys. Rev. **85**, 631 (1952), URL http://journals.aps.org/pr/abstract/10.1103/PhysRev.85.631.

53. T. Heinzl, A. Ilderton, and B. King, Phys. Rev. Lett. **127**, 061601 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.127.061601.

54. J. B. Kogut and E. Dagotto, Phys. Rev. Lett. **59**, 617 (1987), URL https://link.aps.org/doi/10.1103/PhysRevLett.59.617.
55. S. Borsanyi, S. Durr, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, L. Lellouch, T. Lippert, A. Portelli, K. K. Szabo, et al., Science **347**, 1452 (2015), ISSN 0036-8075, https://science.sciencemag.org/content/347/6229/1452.full.pdf, URL https://science.sciencemag.org/content/347/6229/1452.

56. D. K. Sinclair and J. B. Kogut, *Lattice QED in external electromagnetic fields* (2021), URL https://arxiv.org/abs/2111.01990.

57. B. Fricke, J. Desclaux, and J. Waber, Phys. Rev. Lett. **28**, 714 (1972), URL https://doi.org/10.1103/PhysRevLett.28.714.

58. J. P. Desclaux, At. Data, Nucl. Data Tab. **12**, 311 (1973), URL https://doi.org/10.1016/0092-640X(73)90020-X.

59. J. Sucher, Phys. Rev. A **22**, 348 (1980), URL https://link.aps.org/doi/10.1103/PhysRevA.22.348.

60. C. F. Fischer, M. Godefroid, T. Brage, P. Jönsson, and G. Gaigalas, Journal of Physics B: Atomic, Molecular and Optical Physics **49**, 182004 (2016), URL https://doi.org/10.1088/0953-4075/49/18/182004.

61. I. P. Grant, in *Relativistic Effects in Atoms, Molecules, and Solids*, edited by G. L. Malli (Plenum Press, New York, 1983), pp. 73–88.

62. W. R. Johnson, *Atomic structure theory* (Springer, 2007).

63. K. T. Cheng, M. H. Chen, W. R. Johnson, and J. Sapirstein, Canadian Journal of Physics **86**, 33 (2008), ISSN 0008-4204, 1208-6045, URL http://www.nrcresearchpress.com/doi/10.1139/p07-106.

64. P. Indelicato, J. Phys. B: At. Mol. Opt. Phys. **19**, 1719 (1986), URL https://doi.org/10.1088/0022-3700/19/12/012.

65. P. Indelicato, O. Gorceix, and J. P. Desclaux, J. Phys. B: At. Mol. Opt. Phys. **20**, 651 (1987), URL http://dx.doi.org/10.1088/0022-3700/20/4/007.

66. P. Indelicato and J. P. Desclaux, Phys. Rev. A **42**, 5139 (1990), URL http://journals.aps.org/pra/abstract/10.1103/PhysRevA.42.5139.

67. V. V. Flambaum and J. S. M. Ginges, Phys. Rev. A **72**, 052115 (2005), URL https://link.aps.org/doi/10.1103/PhysRevA.72.052115.

68. V. M. Shabaev, I. I. Tupitsyn, and V. A. Yerokhin, Phys. Rev. A **88**, 012513 (2013), URL https://link.aps.org/doi/10.1103/PhysRevA.88.012513.

69. P. Pyykkö and L.-B. Zhao, Journal of Physics B: Atomic, Molecular and Optical Physics **36**, 1469 (2003), ISSN 0953-4075, 1361-6455, URL http://stacks.iop.org/0953-4075/36/i=8/a=302?key=crossref.9362092d92024faa2ceaa5016292b8ed.

70. H. M. Pilkuhn, *Relativistic Quantum Mechanics: Theory of Atomic Bound States* (Kluwer Academic Publishers Group., 2008).

71. G. Hardekopf and J. Sucher, Phys. Rev. A **31**, 2020 (1985), URL https://link.aps.org/doi/10.1103/PhysRevA.31.2020.

72. W. Greiner, B. Müller, et al., *Gauge theory of weak interactions*, vol. 5 (Springer, 1996).

73. J. Mourad and H. Sazdjian, Journal of Physics G: Nuclear and Particle Physics **21**, 267 (1995), URL https://doi.org/10.1088/0954-3899/21/3/004.

74. O. Gorceix and P. Indelicato, Phys. Rev. A **37**, 1087 (1988), URL https://link.aps.org/doi/10.1103/PhysRevA.37.1087.

75. L. F. Pašteka, E. Eliav, A. Borschevsky, U. Kaldor, and P. Schwerdtfeger, Phys. Rev. Lett. **118**, 023002 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.118.023002.