GRASP: The Future?

Abstract

The theoretical foundations of relativistic electronic structure theory within quantum electrodynamics (QED) and the computational basis of the atomic structure code GRASP are briefly surveyed. A class of four-component basis set is introduced, which we denote the CKG-spinor set, that enforces the charge-conjugation symmetry of the Dirac equation. This formalism has been implemented using the Gaussian function technology that is routinely used in computational quantum chemistry, including in our relativistic molecular structure code, BERTHA. We demonstrate that, unlike the kinetically matched two-component basis sets that are widely employed in relativistic quantum chemistry, the CKG-spinor basis is able to reproduce the well-known eigenvalue spectrum of point-nuclear hydrogenic systems to high accuracy for all atomic symmetry types. Calculations are reported of third- and higher-order vacuum polarization effects in hydrogenic systems using the CKG-spinor set. These results reveal that Gaussian basis set expansions are able to calculate accurately these QED effects without recourse to the apparatus of regularization and in agreement with existing methods. An approach to the evaluation of the electron self-energy is outlined that extends our earlier work using partial-wave expansions in QED. Combined with the treatment of vacuum polarization effects described in this article, these basis set methods suggest the development of a comprehensive ab initio approach to the calculation of radiative and QED effects in future versions of the GRASP code.

1. Introduction

Electronic structure theory has driven the development of atomic, molecular and solid state physics for almost a century. The nonrelativistic self-consistent field model pioneered by Douglas Hartree in the pre-computer era formed the foundations of atomic structure theories and a theoretical understanding of the structure of the Periodic Table of the elements. This model was incorporated into atomic structure computer codes which have remained under constant development for more than half a century. These codes include the nonrelativistic ATSP and relativistic GRASP programs that are based on the multiconfigurational Hartree–Fock formalism described in the well-known 1957 lectures by Hartree [1].

The development of relativistic electronic structure theories has been motivated by the need to incorporate physical effects not described by the Schrödinger equation, which are usually most evident in systems containing heavy elements. In this article, we restrict our attention to electronic structure models based on the four-component Dirac equation for mean-field potentials [2,3], with subsequent treatment of many-body effects using perturbative, configurational or coupled-cluster expansions, or the two-particle reduced density matrix approach. The formalisms for these computational approaches are equally applicable to the determination of atomic and molecular structures and are widely employed in existing computational codes.

The developers of the GRASP codes, of which GRASP2018 [4,5] is the most recent release, have been motivated both by the experimental demands of atomic spectroscopy and by the desire to formulate computational methods that are consistent with the underlying theory of quantum electrodynamics (QED), achieving both accuracy and completeness. The structure of the code is sketched in Section 2. It has become a widely-used tool in atomic spectroscopy, with many applications in solar, stellar and plasma physics and as a platform for the exploration of relativistic, many-body and QED effects, especially in systems containing heavy elements. It has also been used to gain insights into molecular and solid state physics, in investigations of particle physics beyond the Standard Model, and in recent searches for the experimental signatures of dark matter.

The relativistic molecular electronic structure code BERTHA [6,7] has been developed within the same theoretical framework as GRASP. These common foundations are described in detail in [8] and are summarized in Section 3. Potential synergies between GRASP and BERTHA that exploit their common structures and shared formulation form the basis of this article. Both codes generate numerical four-component solutions of the Dirac equation which define distinct non-interacting fields of electrons or positrons. Furry provided a covariant derivation of the bound-state interaction picture of QED within the Feynman-Dyson formalism [9]. Charged particles in the Furry picture interact only through coupling to the Maxwell field, emitting and absorbing virtual photons. Both GRASP and BERTHA, which assume that the exchange of energy involves single photons, represented by Coulomb, Breit or other interactions, are computational implementations of the Furry picture. Radiative corrections, electron self-energy and vacuum polarization effects, involve one-photon exchanges between the uncoupled electron and positron fields.

In common with most research in relativistic electronic structure theory (see, for example [10,11,12,13,14,15,16], the development of GRASP and BERTHA has focussed on the construction of positive-energy bound-states and many-body corrections that exclude explicit reference to the negative-energy sector of the spectrum. This approach does not exclude consideration of QED corrections, but restricts their calculation to the evaluation of bound-state expectation values of effective interaction operators, examples of which have been implemented in GRASP. In Section 3.3.3, we consider a new class of four-component basis set that is consistent with the Furry picture, the charge-conjugation symmetries of the free-particle Dirac equation and the relativistic Gaussian basis set technologies that we have developed for molecular physics and quantum chemistry. We demonstrate in Section 4 that this basis set is able to calculate the vacuum polarization effects attributed to the Wichmann-Kroll interaction from a first-principles implementation of QED. This is achieved without recourse to the use of effective potentials or the machinery of regularization that is usually invoked to handle divergences in QED [17]. This formulation echoes our approach to the calculation of the electron self-energy [18,19,20] and the work of Lindgren and collaborators [21,22], but is based on quite different computational methods that have evolved through our work on molecular structure. We conclude by suggesting that these basis set methods may be usefully employed in future developments of GRASP.

2. Building Atomic Structure Programs

The software repository [5] of the CompAS (Computational Atomic Structure) collaboration [23] supplies the nonrelativistic ATSP and the relativistic GRASP programs, both based on multiconfiguration Hartree–Fock ideas expounded in Hartree’s 1957 lectures [1]. Bound electron orbitals, the building blocks of nonrelativistic configurational states (CSF) in ATSP have the familiar form:

$$\begin{array}{cc}\hfill {\psi }_{nlm\sigma }\left(\mathbf{r}\right)& =P_{nl}\left(r\right)Y_{lm}(\theta ,\phi ){\phi }_{\sigma }\hfill \\ \hfill l=0,1,2,\dots ,& m=l,l-1,\dots ,-l,\sigma =\pm 1/2,\hfill \\ \hfill {\phi }_{1/2}=& \left(\begin{array}{c}1\\ 0\end{array}\right),{\phi }_{-1/2}=\left(\begin{array}{c}0\\ 1\end{array}\right)\hfill \end{array}$$

(1)

in spherical polar coordinates, where $l,m$ are orbital angular momentum quantum numbers and $\sigma =\pm 1/2$ the spin projection. The label n enumerates successive bound orbitals of this symmetry, replaced by the energy, $\epsilon$, above threshold for ionization to the continuum. We shall refer to the $2\times (2l+1)$ orbitals labelled $nlm\sigma$ as the $nl$-subshell, with the common radial factor $P_{nl}\left(r\right)$.

In relativistic programs like GRASP these are replaced by 4-component spinors [2,8]

$$\begin{array}{c}\hfill {\psi }_{n\kappa m}\left(\mathbf{r}\right)={\displaystyle \frac{1}{r}}\left(\begin{array}{c}{P}_{n\kappa }\left(r\right){\chi }_{\kappa ,m}(\theta ,\phi )\\ i{Q}_{n\kappa }\left(r\right){\chi }_{-\kappa ,m}(\theta ,\phi )\end{array}\right),\end{array}$$

(2)

where the 2-component coupled spin–orbit function is

$$\begin{array}{c}\hfill {\chi }_{\kappa m}(\theta ,\phi )=\sum _{\sigma =\pm 1/2}\langle l,1/2,j,m|l,m-\sigma ,1/2,\sigma \rangle Y_{l,m-\sigma }(\theta ,\phi ){\varphi }_{\sigma },\\ \hfill l=j+\eta /2,\kappa =\eta (j+1/2),(\eta =\pm 1),m=-j,\dots ,j.\end{array}$$

(3)

In this case each $n\kappa$-subshell consists of $2j+1=2\left|\kappa \right|$ orbitals, with common radial functions $P_{n\kappa }\left(r\right),Q_{n\kappa }\left(r\right)$. It is often convenient to replace the $n\kappa$ label with the more familiar set $nljm$ with $j=l\pm 1/2$, as $P_{n\kappa }\to P_{nl}\left(r\right)$, and $Q_{n\kappa }\sim O(1/c)$ in the limit $c\to \infty$.

These orbital functions are the building blocks for many-electron wave functions. The simplest is the Slater determinant

$$\begin{array}{c}\hfill \langle x_1,\dots ,x_n|j_1,j_2,\dots ,j_n\rangle ={\displaystyle \frac{1}{\sqrt{n!}}}\left|\begin{array}{cccc}{\psi }_{j_1}\left(x_1\right)& {\psi }_{j_2}\left(x_1\right)& \dots & {\psi }_{j_n}\left(x_1\right)\\ {\psi }_{j_1}\left(x_2\right)& {\psi }_{j_2}\left(x_2\right)& \dots & {\psi }_{j_n}\left(x_2\right)\\ ⋮& ⋮& ⋮& \\ {\psi }_{j_1}\left(x_n\right)& {\psi }_{j_2}\left(x_n\right)& \dots & {\psi }_{j_n}\left(x_n\right)\end{array}\right|,\end{array}$$

(4)

that vanishes if more than one electron orbital is assigned the same state label $j_r$ or the same location $x_s$, thus ensuring that no orbital has more than one occupant in accordance with Pauli’s exclusion principle. GRASP programs classify configurational states (CSF)—notionally a linear combination of Slater determinants—in $jj$-coupling: $\Phi (\gamma ,J^{\pi },M)$ denotes an n-electron CSF with total angular momentum quantum numbers $J,M$, and parity $\pi$; $\gamma$ includes all additional information, in particular the angular momentum coupling scheme, needed to specify $\Phi (\gamma ,J^{\pi },M)$ uniquely using Racah’s theory of angular momentum ([8], Sections 6.6–6.8). CSF in ATSP, $\Phi (\gamma ,L,M_L,S,M_S)$, are formed in much the same way. The atomic state functions, ASFs, are linear combinations of CSFs of the same symmmetry, $\Psi (\Gamma ,J^{\pi },M)={\sum }_{\gamma }{c}_{\Gamma \gamma }\Phi (\gamma ,J^{\pi },M)$.

Given estimates of the orbitals, we can construct matrix elements ${\mathit{H}}_{\gamma {\gamma }^{\prime }}$ of the Hamiltonian of an N-electron atom with respect to the CSFs $\Phi (\gamma ,J^{\pi },M)$. Diagonalization of $\mathit{H}$ generates ASFs $\Psi (\Gamma ,J^{\pi },M)$ and ASF energies $E(\Gamma ,J^{\pi })$ with the same spectroscopic classification.

The GRASP orbital equations are described in ([8], Chapter 7). Here, as in [8], we adopt Hartree atomic units for which $e=\hslash =1/4\pi {ϵ}_0=m_e=1$, and $c\simeq 137$. We also adhere to the notational conventions of [8], in that the $4\times 4$ Dirac matrices $\mathit{\alpha }$ and $\beta$ are employed when discussing computational implementations of relativistic electronic structure theory, and the $4\times 4$${\gamma }_{\mu }$ matrices appear when discussing the formulation of the theory within quantum electrodynamics. The correspondence between these notations and their properties are discussed in ([8], Appendix A.2) and summarized by the relations ${\gamma }_0=\beta$ and ${\gamma }_{\mu }=\beta {\alpha }_{\mu }\mathrm{for}\mu =1,2,3$. If the electron mass, $m_e$, and charge, e, appear in formulae it is to emphasise dimensionality, but it is implicitly understood that they should be assigned unit value. The (unquantized) multi-configurational Dirac-Hartree–Fock–Breit (MCDHFB) Hamiltonian for an N-electron atom with point charge nucleus Z has the form

$$\begin{array}{c}\hfill H=\sum _{i=1}^N(h_i-Z/r_i)+\sum _{i<j}{V}_{ij},\end{array}$$

(5)

where

$$\begin{array}{c}\hfill h_i=c{\mathit{\alpha }}_i·{\mathit{p}}_i+({\beta }_i-1)c^2\end{array}$$

(6)

is the Dirac Hamiltonian for electron i. The electron-electron interaction potential $V_{ij}$ is just the classical Coulomb potential, augmented by the Breit interaction. A standard variational argument generates coupled sets of integro-differential equations, one for each participating subshell A, of the form ([8], Section 7.3)

$$\begin{array}{c}\hfill \left[\begin{array}{cc}-{\displaystyle \frac{Z_A\left(r\right)}{r}}-{\epsilon }_{AA}& c\left(-{\displaystyle \frac{d}{dr}}+{\displaystyle \frac{{\kappa }_A}{r}}\right)\\ c\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{{\kappa }_A}{r}}\right)& -2c^2-{\epsilon }_{AA}-{\displaystyle \frac{Z_A\left(r\right)}{r}}\end{array}\right]\left[\begin{array}{c}{P}_A\left(r\right)\\ Q_A\left(r\right)\end{array}\right]=-{\displaystyle \frac{1}{r}}\left[\begin{array}{c}{X}_A^{+}\left(r\right)\\ X_A^{-}\left(r\right)\end{array}\right].\end{array}$$

(7)

The potential energy $-Z_A\left(r\right)/r$ is the sum of the nuclear Coulomb potential energy $-Z/r$ and, the ‘direct’ (or classical) Coulomb repulsion from subshells $B,C,\dots$. The exchange repulsion energy from other electrons in subshell A is also included in $Z_A\left(r\right)$. The quantities $X_A^{\pm }\left(r\right)/r$ on the right are ‘exchange’ potentials together with any Lagrange multipliers, ${\epsilon }_{AB},\dots$, required to preserve orbital orthonormality. Each such equation, given estimates of the interaction with other subshells, provides an estimate of $P_A\left(r\right),Q_A\left(r\right)$. They can be used to modify the interaction potentials for the next iteration of the equations until all orbitals have stabilized to agreed precision.

Equation (7) recognizes only subshells containing at least one electron. Other one-electron states of that symmetry are ignored, in particular the negative-energy solutions representing positrons in the atomic field. These are never calculated in MCDHFB models but, though ignored, they are implicitly part of the model.

Numerical methods to solve Equation (7) in GRASP and similar programs in the Hartree tradition are mostly based on finite differences although a number of recent programs have used B-spline expansions. Relativistic molecular calculations require basis set methods proposed by Roothaan [24] and Hall [25]. This approach has rarely been used for atomic problems, and there were many failed attempts to construct 4-spinor basis sets before the need for kinetic matching of the two radial components was established ([8], Section 5.7). The BERTHA relativistic molecular package [6,7] relies on the atomic kinetically-matched spinor basis sets described in Section 3.3.1. Four-spinor basis set algorithms automatically generate both electron and positron wave functions; the positron solutions play no part in atomic and molecular structure. Calculation of additional contributions suggested by QED require a computational scheme with further symmetry constraints on the 4-spinor basis functions to balance the representation of electrons and positrons.

3. QED of Atoms and Molecules

3.1. Relativistic Wave Equations

The early history of relativistic wave equations is well documented by, for example, the books [26,27,28]. Dirac’s ‘negative energy states’ and their stop-gap ‘hole theory’ [29,30,31,32,33] made it impossible to regard their relativistic wave equation as describing motion of a single particle in the manner of the Schrödinger equation, inspiring the development of quantum field theory. The unquestioned acceptance of the physical reality of Dirac’s ‘hole theory’ created major problems in the design of mathematical and computational algorithms for relativistic atomic and molecular structure. Variational methods of nonrelativistic quantum mechanics are routinely justified by the existence of a finite lower bound to the atomic or molecular Hamiltonian, so that iterative decreasing sequences of level energies may converge to a limit according to Cauchy’s general principle of convergence. The spectra of Dirac Hamiltonians are infinite above and below leading to an extensive literature regarding problems associated with the application of the linear variation method to relativistic electronic structure theory. A review of this literature and a discussion of both the linear variation method and the application of min-max principles to relativistic electronic structure problems may be found in ([8], Section 5.5). This is where relativistic spinor adaptations of the ideas of Roothaan [24] and Hall [25] and an understanding of their symmetry structure came to the rescue.

3.2. Quantized Electron and Positron Fields

This section follows ([28], Section 14.1). The Dirac equation for free particles (in atomic units) is

$$\begin{array}{c}\hfill (i{\gamma }^{\mu }{\partial }_{\mu }-m_{e}c)\psi \left(x\right)=0,x=\{x^0,x^1,x^2,x^3\}\equiv (ct,\mathit{x}).\end{array}$$

(8)

We can think of this as a field equation, interpreting $\psi \left(x\right)$ as an operator on the states of the Dirac field rather than providing particle wave functions. Label the states of the Dirac field ${\Phi }_n$, where the index n runs over all possible field modes, with the vacuum state ${\Phi }_0$. We define the wave functions$U_r\left(x\right),V_s\left(x\right)$ such that

$$\begin{array}{c}\hfill U_r\left(x\right)=\langle {\Phi }_0,\psi \left(x\right){\Phi }_r\rangle ,V_s\left(x\right)=\langle {\Phi }_s,\psi \left(x\right){\Phi }_0\rangle \end{array}$$

(9)

as solutions of the equations

$$(i{\gamma }^{\mu }{\partial }_{\mu }-m_{e}c)U_r\left(x\right)=0,(i{\gamma }^{\mu }{\partial }_{\mu }-m_{e}c)V_s\left(x\right)=0.$$

Invariance with respect to time translation allows us to write

$$\begin{array}{c}\hfill U_r\left(x\right)=e^{-i{E}_{r}t}{u}_r\left(\mathit{x}\right),(E_r>0),V_s\left(x\right)=e^{+i{E}_{s}t}{v}_s\left(\mathit{x}\right),(E_s>0)\end{array}$$

(10)

for positive and negative energy modes, respectively, so that the $\mathit{x}$-dependent amplitudes satisfy the time-independent Dirac equations

$$\begin{array}{c}\hfill (c\mathit{\alpha }·\mathit{p}+\beta m_e{c}^2)u_r\left(\mathit{x}\right)=E_r{u}_r\left(\mathit{x}\right),(c\mathit{\alpha }·\mathit{p}+\beta m_e{c}^2)v_s\left(\mathit{x}\right)=-E_s{v}_s\left(\mathit{x}\right).\end{array}$$

(11)

Expand the field operators in the form

$$\begin{array}{c}\hfill \psi \left(x\right)=\sum _r{a}_r{U}_r\left(x\right)+\sum _s{b}_s^{†}{V}_s\left(x\right),\\ \hfill {\psi }^{†}\left(x\right)=\sum _r{a}_r^{†}{U}_r^{†}\left(x\right)+\sum _s{b}_s{V}_s^{†}\left(x\right)\end{array}$$

(12)

with anticommuting creation and annihilation operators such that

$$\begin{array}{c}\hfill \{a_r^{†},a_{r^{\prime }}\}={\delta }_{r,r^{\prime }},\{b_s^{†},b_{s^{\prime }}\}={\delta }_{s,s^{\prime }};\end{array}$$

(13)

All other anticommutators vanish. If we set the equal-time anticommutator of the field operators ${\psi }^{†}\left(x\right),\psi \left(x\right)$ to be

$$\begin{array}{c}\hfill \{\psi (ct,\mathit{x}),\tilde{\psi }(ct,\mathit{y})\}={\gamma }^0{\delta }^3(\mathit{x}-\mathit{y})\end{array}$$

(14)

where $\tilde{\psi }\left(x\right)={\psi }^{†}\left(x\right){\gamma }^0$; then

$$\begin{array}{c}\hfill \sum _r{u}_r(ct,\mathit{x})u_r^{†}(ct,\mathit{y})+\sum _s{v}_s(ct,\mathit{x})v_s^{†}(ct,\mathit{y})={\delta }^3(\mathit{x}-\mathit{y})\end{array}$$

(15)

so that the complete spectrum of the Dirac free particle field operator has both positive- and negative-energy solutions. The formalism implies that electrons and positrons obey Fermi statistics: there can be only one electron in state r or positron in state s. In all physically realistic models, the electron and positron spectra are quite distinct, although their wave functions satisfy the same Dirac equation. This formalism makes it possible to think of separate electron and positron fields. Then

$$\begin{array}{c}\hfill N^e=\sum _r{n}_r,E^e=\sum _r{n}_r{E}_r,N^p=\sum _s{n}_s,E^p=\sum _s{n}_s{E}_s,\end{array}$$

(16)

are, respectively, the total number and total energy of the quanta of the electron and positron fields. The combined totals are

$$N=N^e+N^p,E=E^e+E^p.$$

This construction can also be used with minimal coupling to an external electromagnetic field, $a_{\mu }\left(x\right)$, so that

$$\begin{array}{c}\hfill ({\gamma }^{\mu }(i{\partial }_{\mu }-e{a}_{\mu }/c)+\beta m_e{c}^2)\psi \left(x\right)=0,a_{\mu }=(c\varphi \left(\mathit{x}\right),\mathit{a}\left(\mathit{x}\right)).\end{array}$$

(17)

The electron field and the positron fields are again in separate domains. The charge conjugation symmetry is broken in an atomic environment as shown schematically in Figure 1.

Quantization of the electron and positron fields is independent of the particular method of solution of the Dirac equation. Figure 1 illustrates schematically the symmetry of the free fermion levels $n^{+}$ and $n^{-}$ about the relativistic energy zero. The charge conjugation operator, C of Section 3.3.2, sets up a bijective mapping between the $n^{+}$ and $n^{-}$ energy states of free fermions. This symmetry is broken in the atomic or molecular environment; all realistic physical models identify separate electron and positron domains, which we can think of as defining separate quantum fields although the wave functions are solutions of the same Dirac equation.

3.3. Basis Set Spinor Expansions

The formalism in current use for solving relativistic atomic and molecular structure problems using a linear combination of spinor basis functions is surveyed in ([8], Sections 5.6–5.11). An orbital a in subshell A having the structure (2) is approximated by

$$\begin{array}{c}\hfill {\psi }_a\left(\mathit{x}\right)=\left(\begin{array}{c}{\displaystyle \sum _{\mu =1}^N{c}_{\mu ,a}^{+}M[+,\mu ;\mathit{x}]}\\ {\displaystyle \sum _{\mu =1}^N{c}_{\mu ,a}^{-}M[-,\mu ;\mathit{x}]}\end{array}\right)\end{array}$$

(18)

where $\mu$ labels the large- and small-component basis spinors, $M[\pm ,\mu ;\mathit{x}]$, which in a molecule may be centred on more than one nucleus. The time-independent Dirac Hamiltonian for an electron can be written as in (11)

$$\begin{array}{c}\hfill {\tilde{h}}_D=\mathcal{K}+\mathcal{M}+\mathcal{V},\mathcal{K}=c\mathit{\alpha }·\mathit{p},\mathcal{M}=\beta c^2,\mathcal{V}=V\left(\mathit{x}\right)\end{array}$$

(19)

where the potential energy $V\left(\mathit{x}\right)$ may include contributions from interactions with other electrons and nuclei. We can now form a Rayleigh quotient by taking the expectation of ${\tilde{h}}_D$ with respect to ${\psi }_a\left(\mathit{x}\right)$:

$$\begin{array}{c}\hfill \mathcal{R}=\langle {\mathit{c}}_a^{†}\mathit{H}{\mathit{c}}_a\rangle /\langle {\mathit{c}}_a^{†}\mathit{S}{\mathit{c}}_a\rangle ,{\mathit{c}}_a=\left(\begin{array}{c}{\mathit{c}}_a^{+}\\ {\mathit{c}}_a^{-}\end{array}\right),\end{array}$$

(20)

where the elements of the N-vectors ${\mathit{c}}_a^{\pm }$ are the expansion coefficients of (18),

$$\mathit{H}=\left(\begin{array}{cc}{c}^2{\mathit{S}}^{++}+{\mathit{V}}^{++}& c{\mathit{K}}^{+-}\\ c{\mathit{K}}^{-+}& -c^2{\mathit{S}}^{--}+{\mathit{V}}^{--}\end{array}\right),\mathit{S}=\left(\begin{array}{cc}{\mathit{S}}^{++}& 0\\ 0& {\mathit{S}}^{--}\end{array}\right)$$

The $N\times N$ matrices ${\mathit{K}}^{\beta ,{\beta }^{\prime }}={{\mathit{K}}^{{\beta }^{\prime },\beta }}^{†},{\mathit{V}}^{\beta ,\beta },{\mathit{S}}^{\beta ,\beta }$ have elements

$$\begin{array}{cc}\hfill S_{\mu \nu }^{\beta {\beta }^{\prime }}& ={\delta }_{\beta {\beta }^{\prime }}\int M{[\beta ,\mu ;\mathit{x}]}^{†}M[{\beta }^{\prime },\nu ;\mathit{x}]d\mathit{x}\hfill \\ \hfill V_{\mu \nu }^{\beta {\beta }^{\prime }}& ={\delta }_{\beta {\beta }^{\prime }}\int M{[\beta ,\mu ;\mathit{x}]}^{†}V\left(\mathit{x}\right)M[{\beta }^{\prime },\nu ;\mathit{x}]d\mathit{x}\hfill \\ \hfill K_{\mu \nu }^{\beta {\beta }^{\prime }}& ={\delta }_{\beta ,-{\beta }^{\prime }}\int M{[\beta ,\mu ;\mathit{x}]}^{†}\mathbf{\sigma }·\mathit{p}M[-\beta ,\nu ;\mathit{x}]d\mathit{x}\hfill \end{array}$$

so that the matrix $\mathit{H}$ is Hermitian. The Rayleigh-Ritz variational method gives Galerkin equations

$$\begin{array}{c}\hfill \left(\begin{array}{cc}(c^2-E){\mathit{S}}^{++}+{\mathit{V}}^{++}& c{\mathit{K}}^{+-}\\ c{\mathit{K}}^{-+}& (-c^2-E){\mathit{S}}^{--}+{\mathit{V}}^{--}\end{array}\right)\left(\begin{array}{c}{\mathit{c}}^{+}\\ {\mathit{c}}^{-}\end{array}\right)=0.\end{array}$$

(21)

This generalized eigenvalue equation for $2N$ eigenvalues and eigenfunctions requires basis 4-spinors to have specific relations between their components to give physically useful results. Following (2), we assume

$$\begin{array}{c}\hfill M[+,\mu ;\mathit{x}]=[f_{\kappa }^{+}\left(r\right)/r].{\chi }_{\kappa ,m}(\theta ,\phi ),M[-,\mu ;\mathit{x}]=i[f_{\kappa }^{-}\left(r\right)/r].{\chi }_{-\kappa ,m}(\theta ,\phi ).\end{array}$$

(22)

The kinetic matching restriction ([8], Section 5.7) of the radial components

$$\begin{array}{c}\hfill f_{\kappa }^{-}\left(r\right)={\displaystyle \frac{1}{2c}}\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{\kappa }{r}}\right)f_{\kappa }^{+}\left(r\right)\end{array}$$

(23)

ensures that the matrix $\mathit{H}$ reduces to the corresponding $N\times N$ Schrödinger matrix with respect to the radial functions $f_{\kappa }^{+}\left(r\right)$ for positive-energy states in the nonrelativistic limit. More specifically, this construction—we refer to ([8], Chapter 5) for the details—gives the nonrelativistic kinetic energy matrix ([8], Equation (5.7.13)).

$$\begin{array}{c}\hfill {\mathit{T}}_{nr}^{++}={\mathit{K}}^{+-}{\left({\mathit{S}}^{--}\right)}^{-1}{\mathit{K}}^{-+}={\mathit{S}}^{--},\end{array}$$

(24)

so that (21) reduces to the matrix Schrödinger equation

$$\begin{array}{c}\hfill [{\mathit{T}}_{nr}^{++}+{\mathit{V}}^{++}-ϵ{\mathit{S}}^{++}]{\mathit{c}}^{+}=0.\end{array}$$

(25)

where $ϵ=E-m{c}^2$ is the eigenvalue relative to the nonrelativistic zero of energy. Note that kinetic matching is not a direct relativistic effect. It is a consequence of the internal structure of the Dirac operator, essential for ensuring correct behaviour in the nonrelativistic limit. A similar argument reveals that (23) holds with a change of sign for negative energy solutions, leading to the Schrödinger equation $[{\mathit{T}}_{nr}^{--}-{\mathit{V}}^{--}-{ϵ}^{\prime }{\mathit{S}}^{--}]{\mathit{c}}^{-}=0$ with the repulsive potential $[-V^{--}]$ for low energy positrons with energy ${\epsilon }^{\prime }>0$. Successful types of kinetically matched spinor basis functions are listed in ([8], Sections 5.8–5.10) along with numerical examples.

3.3.1. KG-Spinors

The fact that the upper (“large”) radial component of a 4-spinor may be of any form used in nonrelativistic atomic models defines the class of K-spinors (K denotes ‘kinetically matched’) in which the radial components are related by (23). We consider only the Gaussian (KG-spinor) sets (G denotes ‘Gaussian function’), from which the BERTHA relativistic molecular program [6] is built:

$$\begin{array}{c}\hfill f_{\kappa ,i}^{+}\left(r\right)=\left\{\begin{array}{cc}\hfill r^K{e}^{-{\lambda }_{K,i}{r}^2},\hfill & \kappa =-K\\ r^{K+1}{e}^{-{\lambda }_{K,i}{r}^2},& \kappa =+K\end{array}\right.\end{array}$$

(26)

where the exponents are ${\lambda }_{K,i}=\alpha {\beta }^{i-1},i=1,2,\dots ,N$ with suitable chosen parameters $\alpha ,\beta ,N$ that may depend on $K=\left|\kappa \right|$. From (23), the lower components are

$$\begin{array}{c}\hfill f_{\kappa ,i}^{-}\left(r\right)=\left\{\begin{array}{cc}\hfill C_{\kappa ,i}{r}^{K+1}{\lambda }_{K,i}{e}^{-{\lambda }_{K,i}{r}^2},\hfill & \kappa =-K\\ C_{\kappa ,i}{r}^K{e}^{-{\lambda }_{K,i}{r}^2}\left((2K+1)-2{\lambda }_{K,i}{r}^2\right),& \kappa =+K\end{array}\right.\hfill \end{array}$$

(27)

where $C_{\kappa ,i}$ are arbitrary constants, chosen to enforce normalization. KG-spinors have been shown to give excellent results for atomic and molecular structure and properties but they lack the symmetry needed to give a balanced representation of negative-energy states.

3.3.2. Charge Conjugation

The charge conjugation matrix C sets up a bijective mapping between the elements of the domains of free electrons and positrons. In the Dirac representation, this matrix is given by ([8], Appendix A.2):

$$\begin{array}{c}\hfill \psi \to {\psi }_c=C{\overline{\psi }}^t,C=i{\gamma }^2{\gamma }^0=\left(\begin{array}{cc}0& -i{\sigma }^2\\ -i{\sigma }^2& 0\end{array}\right),C^t=C^{†}=-C,\end{array}$$

(28)

where ${\overline{\psi }}^t={\psi }^{*}{\gamma }^0$ is Dirac conjugation. The asterisk denotes complex conjugation, and superscript t denotes vector transposition. Minimal coupling of a particle with charge q to an external electromagnetic field $a_{\mu }$ replaces (8) by

$$\begin{array}{c}\hfill (i{\gamma }^{\mu }{\partial }_{\mu }-m_{e}c)\psi \left(x\right)={\displaystyle \frac{q}{c}}{\gamma }^{\mu }{a}_{\mu }\psi \left(x\right).\end{array}$$

(29)

A simple calculation shows that

$$\begin{array}{c}\hfill (i{\gamma }^{\mu }{\partial }_{\mu }-m_{e}c){\psi }_c\left(x\right)=-{\displaystyle \frac{q}{c}}{\gamma }^{\mu }{a}_{\mu }{\psi }_c\left(x\right).\end{array}$$

(30)

so that charge conjugation reverses the sign of the particle’s charge. The expectation value of an operator $\mathcal{O}$ is given by

$$\begin{array}{c}\hfill \langle \mathcal{O}\rangle =\int \overline{\psi }\mathcal{O}\psi d^{3}x,{\langle \mathcal{O}\rangle }_c=\int {\overline{\psi }}_c\mathcal{O}{\psi }_c{d}^{3}x,\end{array}$$

(31)

so that for free particles we have

$$\begin{array}{cc}\hfill {\langle x^{\mu }\rangle }_c=\langle x^{\mu }\rangle ,{\langle p^{\mu }\rangle }_c\hfill & \hfill =-\langle p^{\mu }\rangle ,{\langle j^{\mu }\rangle }_c=-\langle j^{\mu }\rangle ,\hfill \\ \hfill {\langle \mathbf{\Sigma }\rangle }_c=-\langle \mathbf{\Sigma }\rangle ,{\langle \mathit{L}\rangle }_c\hfill & \hfill =-\langle \mathit{L}\rangle ,{\langle \mathit{J}\rangle }_c=-\langle \mathit{J}\rangle ,\hfill \end{array}$$

These electron/positron symmetry relations are broken by interaction with other particles in atoms and molecules. Figure 1 compares the spectrum of the Dirac operator for free particles on the left, with a typical atomic potential on the right, generated using a basis set which respects free particle charge conjugation. The free spectrum is symmetric about energy $E=0$, while the atomic potential lowers all eigenvalues, introducing bound states and making $E_{n^{+}}\ne -E_{n^{-}}$.

Atomic and molecular structure calculations write the Dirac Hamiltonian in the form

$$\begin{array}{c}\hfill H_D\left(q\right):=c\mathit{\alpha }·(\mathit{p}-q\mathit{A})+\beta m_e{c}^2+q\varphi \end{array}$$

(32)

so that for particles of charge q

$$\begin{array}{c}\hfill {\langle H_D\left(q\right)\rangle }_c=-\langle H_D(-q)\rangle .\end{array}$$

(33)

The standard 4-spinor in spherical coordinates with real radial amplitudes $f_{\kappa }^{\beta }\left(r\right)$ has the form

$$\begin{array}{c}\hfill {\psi }_{\kappa ,m}\left(\mathit{x}\right)={\displaystyle \frac{1}{r}}\left(\begin{array}{c}{f}_{\kappa }^{+}\left(r\right){\chi }_{\kappa m}(\theta ,\phi )\\ i{f}_{\kappa }^{-}\left(r\right){\chi }_{-\kappa m}(\theta ,\phi )\end{array}\right);\end{array}$$

(34)

applying (28) gives the charge conjugate spinor

$$\begin{array}{c}\hfill {\psi }_{-\kappa ,-m}{\left(\mathit{x}\right)}_c=-i{(-1)}^{m+1/2}{\displaystyle \frac{1}{r}}\left(\begin{array}{c}{f}_{\kappa }^{-}\left(r\right){\chi }_{-\kappa ,-m}(\theta ,\phi )\\ i{f}_{\kappa }^{+}\left(r\right){\chi }_{\kappa ,-m}(\theta ,\phi )\hfill \end{array}\right).\end{array}$$

(35)

3.3.3. CKG-Spinors

The charge conjugation operator discussed in Section 3.3.2 imposes a symmetry on four-component Dirac spinors under the transformation $e\to -e$. A complete set that satisfies charge conjugation symmetry must consist of all spin-angular symmetry types, and contain elements classified by positive or negative energy. In a Gaussian representation, we denote such a basis the CKG set (C denotes ’charge conjugate’), which expands an arbitrary Dirac spinor $\psi \left(x\right)$ through the relation:

$$\begin{array}{c}\hfill \psi \left(\mathit{x}\right)=\sum _{\alpha }\sum _{\mu }{c}_{\mu ,\alpha }{G}^{\alpha }[\mu ;\mathit{x}]\end{array}$$

(36)

where $\alpha =\pm 1$ denotes positive- and negative-energy, respectively, and setting $\mu =\{\kappa ,m,i\}$,

$$\begin{array}{c}\hfill G^{\pm }[\kappa ,m,i;\mathit{x}]=N_{\kappa ,i}^{\pm }\left(\begin{array}{c}{P}_{\kappa ,i}^{\pm }\left(r\right){\chi }_{\kappa ,m}(\theta ,\phi )\\ i{Q}_{\kappa ,i}^{\pm }\left(r\right){\chi }_{-\kappa ,m}(\theta ,\phi )\end{array}\right)\end{array}$$

(37)

with normalization constants $N_{\kappa ,i}^{\pm }$. For a given symmetry-type, $\kappa$, we choose a Gaussian basis set for the large components of positive energy free-particle states, from which the corresponding small component functions are derived by application of the free-particle Dirac equation, which fixes the relative amplitudes of the components. For positive energies ($\alpha =+1$), the explicit forms of the radial components are

$$\begin{array}{cc}{P}_{-K,i}^{+}\left(r\right)\hfill & =r^{K}exp[-{\lambda }_{K,i}{r}^2]\hfill \\ Q_{-K,i}^{+}\left(r\right)\hfill & ={(c+E_{-K,i}/c)}^{-1}(-2{\lambda }_{K,i}{r}^{K+1})exp[-{\lambda }_{K,i}{r}^2]\hfill \end{array}$$

(38)

$$\begin{array}{cc}{P}_{+K,i}^{+}\left(r\right)\hfill & =r^{K+1}exp[-{\lambda }_{|K,i}{r}^2]\hfill \\ Q_{+K,i}^{+}\left(r\right)\hfill & ={(c+E_{+K,i}/c)}^{-1}(2K+1-2{\lambda }_{K,i}{r}^2)r^{K}exp[-{\lambda }_{K,i}{r}^2].\hfill \end{array}$$

(39)

The behaviour of these functions at small r resembles that of the analytic solutions of the radial free-particle functions, which are proportional to r times a spherical Bessel function, $j_{\ell }\left(pr\right)$ ([8], Equation (3.2.24)). Coincidently, this is also consistent with the small r behaviour of Dirac spinors for finite nuclear models, for which the potential approaches a constant value.

With these choices for positive-energy radial basis spinors, charge conjugation symmetry dictates the construction of negative-energy radial basis functions to be:

$$\begin{array}{c}\hfill P_{\kappa ,i}^{\pm }\left(r\right)=Q_{-\kappa ,i}^{\mp }\left(r\right),N_{\kappa ,i}^{\pm }=N_{-\kappa ,i}^{\mp }.\end{array}$$

(40)

The parameter $E_{\pm K,i}$, corresponds to the magnitude of the average free-particle energy. In the so-called “dual kinetic balance” basis [12], this parameter is simply set to $c^2$, on the grounds that we are mainly interested in positive-energy bound-states. As high-energy states play a significant role in QED correction calculations, we choose

$$E_{\pm K,i}=+c\sqrt{\langle p_{\pm K,i}^2\rangle +c^2}$$

(41)

where

$$\langle p_{-K,i}^2\rangle =(2K+3){\lambda }_{K,i},\langle p_{+K,i}^2\rangle =(2K+1){\lambda }_{K,i}$$

The normalization constants, $N_{\kappa ,i}^{\pm }$, are chosen so that

$${\int }_0^{2\pi }d\phi {\int }_0^{\pi }sin\vartheta d\vartheta {\int }_0^{\infty }{r}^{2}dr{G}^{\pm †}[\kappa ,m,i;\mathit{x}]G^{\pm }[\kappa ,m,i;\mathit{x}]=1$$

for all valid $\kappa$ and m.

Figure 2 vividly illustrates the necessity of charge conjugation symmetry in a calculation of the $K=1$ contribution to the vacuum charge density. The rapid oscillations of the KG-spinor radial density integrate to zero, whilst the CKG-spinor calculation radial densities cancel within rounding errors for all radii by construction. This ensures charge and parity conservation in interaction matrix elements.

Table 1. CKG-spinor calculation of eigenvalues (a.u.), relative to the nonrelativistic zero of energy, of Hg${}^{79+}$ with a point nucleus. Exact Bohr-Sommerfeld eigenvalues in the second column, followed by computed eigenvalues (units $E_h$) for the three lowest bound states of each symmetry $\kappa$ in order.

$\mathit{\kappa }$	${\mathit{E}}_{\mathit{Analytic}}$	${\mathit{E}}_{{\mathit{n}}_1,\mathit{\kappa }}$	${\mathit{E}}_{{\mathit{n}}_2,\mathit{\kappa }}$	${\mathit{E}}_{{\mathit{n}}_3,\mathit{\kappa }}$
−1	−3532.1921489294	−3532.1921489289	−904.8478012876	−392.0836928780
1	−904.8478012882	−904.8478012878	−392.0836928781	−216.4247478039
−2	−817.8074977480	−817.8074977480	−366.1427114567	−205.5771277760
2	−366.1427114567	−366.1427114567	−205.5771277760	−131.1010555098
−3	−358.9868485160	−358.9868485160	−202.5363034958	−129.6328330777
3	−202.5363034958	−202.5363034958	−129.6328330776	−89.9621990162
−4	−201.0765233582	−201.0765233582	−128.8823613985	−89.5273365309
4	−128.8823613985	−128.8823613985	−89.5273365309	−65.7655876909
−5	−128.4392341889	−128.4392341889	−89.2702733629	−65.6035374887
5	−89.2702733629	−89.2702733629	−65.6035374887	−50.2279904000
−6	−89.1002663743	−89.1002663743	−65.4963124418	−50.1561023779
6	−65.4963124418	−65.4963124419	−50.1561023780	−39.6314776684
−7	−65.4200746697	−65.4200746697	−50.1049761176	−39.5955492450
7	−50.1049761176	−50.1049761176	−39.5955492450	−32.0742810655
−8	−50.0667420260	−50.0667420260	−39.5686766900	−32.0546823612
8	−39.5686766900	−39.5686766900	−32.0546823613	−26.4929690886
−9	−39.5478161969	−39.5478161969	−32.0394670086	−26.4815336906
9	−32.0394670087	−32.0394670087	−26.4815336906	−22.2529707815
−10	−32.0273112566	−32.0273112566	−26.4723972662	−22.2459315332
10	−26.4723972662	−26.4723972662	−22.2459315332	−18.9559506743

reports eigenvalues of a test calculation on the hydrogenic ion Hg${}^{+79}$ obtained with a large CKG-spinor basis. The exponents for all values of $K=\left|\kappa \right|$ were ${\lambda }_{K,i}=\alpha {\beta }^{i-1}\{i=1,2,\dots ,100\}$ with $\alpha =0.001,\beta =1.40$. The table reports only the first three numerical eigenvalues of each symmetry $\kappa$ together with the lowest exact analytic eigenvalue. The first two columns agree to 10 decimal places for all $\kappa$ with the minor exception of the first two lines $\left|\kappa \right|=1$. The atomic symmetry requires that point-nuclear hydrogenic eigenvalues are degenerate for states that share the same $n-$ and $j-$quantum numbers. We find 12 significant figure agreement between the $2s_{1/2}\left(E_{n_2,-1}\right)$ and $2p_{1/2}\left(E_{n_1,1}\right)$ eigenvalues. This pairwise numerical agreement persists across all symmetry-types considered in Table 1 up to $\kappa =\pm 10$; the numerical values of the degenerate pairs of energies coincide with the corresponding Dirac-Sommerfeld eigenvalues.

4. QED Corrections

Vacuum Polarization

A long-standing ambition has been to evaluate numerically on the fly the contribution of vacuum polarization to orbital energies for each atomic state. CKG-spinors defined in Section 3.3.3 make this possible. This is because we can diagonalize the Dirac Hamiltonian

$$h_Z=h_0-Z/r$$

for all Z in terms of a common basis, simplifying evaluation of the matrix elements appearing in Equations (42) and (43) below. The principle of this calculation was demonstrated by Persson et al. [21] who used a similar technology based on partial wave expansion to evaluate the Wichmann-Kroll potential [34,35]. They started from the Feynman diagram, Figure 3, expanding the photon propagator in partial waves ([21], Equation (2.6)).

$$\begin{array}{cc}\hfill U_{VP}\left({\mathit{x}}_1\right)=-\frac{2i}{\pi }& \int d^3{\mathit{x}}_2{\int }_0^{\infty }dk{\int }_{-\infty }^{\infty }\frac{dz}{2\pi }\sum _{l=0}^{\infty }(2l+1)j_l\left(k{r}_1\right){\mathit{C}}^l\left(1\right){\alpha }^{\mu }\left(1\right)\hfill \\ \hfill & \times \sum _t{\displaystyle \frac{{\psi }_t^{†}\left({\mathit{x}}_2\right)j_l\left(k{r}_2\right){\mathit{c}}^l\left(2\right){\alpha }_{\mu }\left(2\right){\psi }_t\left({\mathit{x}}_2\right)}{z-E_t(1-i\eta )}},\hfill \end{array}$$

(42)

where ${\mathit{C}}^l\left(i\right)$ denotes the spherical tensor of order l associated with the coordinate labelled i. Any vacuum polarization acting on the electron in atomic orbital a is due to the closed loop in which electrons circulate in one direction and positrons in the other. In Figure 3B the currents cancel by charge conjugation symmetry so that the photon has zero energy [36]; the external field breaks the charge conjugation symmetry of electrons and positrons as illustrated in Figure 1 so that the exchanged photon in Figure 3A has a finite energy. (Schweber [37] makes this point in a review of the Furry Picture). Equation (42) can be simplified giving

$$\begin{array}{cc}\hfill U_{VP}\left(r\right)=-\frac{1}{\pi }& {\int }_0^{\infty }dk{j}_0\left(kr\right)\sum _{\kappa }(2j_{\kappa }+1)\hfill \\ \hfill & \times \sum _n\mathrm{sgn}\left(E_{n\kappa ,Z}\right)\langle n\kappa ,Z|j_0\left(k{r}_2\right)|n\kappa ,Z\rangle .\hfill \end{array}$$

(43)

where the index t in (42) has been replaced by n. The summation over n includes the complete set of radial functions for each $\kappa$. The case $Z=0$ applies to freely moving electrons, and it is easy to show that pairs of matrix elements related by charge conjugation symmetry cancel. This symmetry is broken when $Z\ne 0$.

Figure 2 has already demonstrated that the CKG-spinor result for the zero potential term, $Z=0$, is consistent with Furry’s theorem. We can also use the CKG-spinors to evaluate the one-potential term:

$$\begin{array}{cc}\hfill V_0^{\left(1\right)}\left({\mathit{r}}_{\mathbf{1}}\right)=-\frac{4}{\pi }& {\int }_0^{\infty }dk{j}_0\left(k{r}_1\right)\sum _{\kappa }(2j_{\kappa }+1)\hfill \\ \hfill & \times \sum _p^{+}\sum _q^{-}{\displaystyle \frac{\langle p\kappa ,0|j_0\left(kr\right)|q\kappa ,0\rangle \langle q\kappa ,0|Z|p\kappa ,0\rangle }{E_{p\kappa ,0}-E_{q\kappa ,0}}},\hfill \end{array}$$

(44)

using the free particle Dirac Hamiltonian, $Z=0$, to evaluate the kets $|p\kappa ,0\rangle$ and free particle eigenvalues $E_{p\kappa ,0}$. CKG-spinor evaluation of the matrix elements $\langle p\kappa ,0\left|j_0\left(kr\right)\right|q\kappa ,0\rangle$ is straightforward. Alternatively, we can use free electron partial waves whose radial amplitudes are proportional to spherical Bessel functions, giving matrix elements proportional to the vertex integrals [19].

Table 2. Partial wave contributions to $\langle 1s\left|U_{VP}\right|1s\rangle$ in Hg${}^{79+}$. The CKG-spinor exponents are ${\lambda }_i=\alpha {\beta }^{i-1}\{i=1,2,\dots N\}$, with $\alpha =0.1,\beta =1.9,N=90$.

$\mathit{K}=\left|\mathit{\kappa }\right|$	$\langle {\mathit{V}}_{\mathit{Z},\mathit{K}}\rangle$	$\langle {\mathit{V}}_{0,\mathit{K}}^{\left(1\right)}\rangle$	$\langle {\mathit{V}}_{\mathit{Z},\mathit{K}}\rangle -\langle {\mathit{V}}_{0,\mathit{K}}^{\left(1\right)}\rangle$
1	3.29166528	3.22434547	6.73198 (−2)
2	2.81964973	2.81325408	6.39565 (−3)
3	2.40679616	2.40522784	1.56831 (−3)
4	2.07156051	2.07095146	6.09053 (−4)
5	1.79881093	1.79850086	3.10045 (−4)
6	1.57229887	1.57211452	1.84348 (−4)
7	1.38083155	1.38071232	1.29222 (−4)
8	1.21719228	1.21711175	8.05279 (−5)
Sum	16.5588053	16.42822183	7.65869 (−2)

compares partial wave contributions to the $1s$ vacuum polarization in Hg${}^{79+}$ with a point nucleus. The first column is the CKG-spinor contribution from (43) using atomic spinors with $Z=80$. Column 2 shows the corresponding unrenormalized one-potential contribution from (44) with the $|p\kappa ,0\rangle$ generated in the same CKG-spinor basis. Following [21], we interpret the difference in the last column as the contribution of Wichmann-Kroll and higher order contributions.

Table 3. Wichmann-Kroll contribution to vacuum polarization of hydrogenic $1s_{1/2}$ electrons. $\langle V^{WK}\rangle$, column 2, is the expectation of the WK-potential from [35] integrated over the analytic $1s_{1/2}$ charge density. Columns 2 and 3 are CKG-spinor results for comparison.

Z	$\langle {\mathit{V}}^{\mathbf{WK}}\rangle$	$\langle {\mathit{V}}_0^{\left(3\right)}\rangle$	$\langle {\mathit{V}}_0^{(3+)}\rangle$
10	3.232201 (−7)	3.740969 (−7)	3.714822 (−7)
20	1.862744 (−5)	2.036149 (−5)	2.042291 (−5)
30	1.964703 (−4)	2.093781 (−4)	2.114400 (−4)
40	1.046683 (−3)	1.098815 (−3)	1.120773 (−3)
50	3.869888 (−3)	4.023198 (−3)	4.158461 (−3)
60	1.144670 (−2)	1.181736 (−2)	1.242842 (−2)
70	2.926072 (−2)	3.005133 (−2)	3.230899 (−2)
80	6.784298 (−2)	6.935930 (−2)	7.667064 (−2)

compares the values of the Wichmann-Kroll contribution for a range of values of Z with the corresponding CKG-spinor calculation. $\langle V^{WK}\rangle$, column 2, is the expectation of the WK-potential defined in [35] integrated over the analytic $1s_{1/2}$ charge density. The value of $\langle V_0^{(3+)}\rangle$ for $Z=80$ in the last column is extrapolated from the sum in Table 2.

We have not yet been able to make comparisons with extensive calculations such as those of Sapirstein and Cheng [38], who used a coordinate space formalism to calculate vacuum polarization energies for hydrogenic ions with finite size nuclei and principal quantum numbers $n=1,\dots ,5$. They also calculated results for lithium-, sodium- and copper-like ions using Kohn-Sham potentials and first-order screening corrections. These are brute force computations. They found it necessary to use an extremely fine radial grid of up to 50,000 points when calculating electron wave functions to control the accuracy of the photon energy integration, where they encounter numerical instabilities giving unphysical behaviour at very high energies $\omega \approx {10}^{6}m{c}^2$ or greater.

5. Electron Self-Energy

There seems no reason why charge symmetric CKG-spinors should not be equally useful in calculating the electron self-energy correction [39,40] to atoms and molecules. The partial wave renormalization approach was developed in [18,19] and, with rather different technology, in [21]. The partial wave expressions to be evaluated for the self-energy of a bound orbital a can be reduced to the general form

$$\begin{array}{c}\hfill E_{\kappa }\left(a\right)=2\pi \alpha \int {\displaystyle \frac{d^{3}k}{{\left(2\pi \right)}^3}}{\displaystyle \frac{1}{k}}\sum _n{\displaystyle \frac{\langle a|{\alpha }_{\mu }{e}^{-i\mathit{k}·\mathit{x}}|n\kappa ,Z\rangle \langle n\kappa ,Z|{\alpha }^{\mu }{e}^{i\mathit{k}·\mathit{x}}|a\rangle }{E_a-E_{n\kappa ,Z}-ck\mathrm{sgn}\left(E_{n\kappa ,Z}\right)}}\end{array}$$

(45)

where the sum over intermediate virtual states $|n\kappa ,Z\rangle$ runs over the complete Dirac spectrum. The integrand is finite for each $|n\kappa ,Z\rangle$, but the partial wave sum over $\kappa$ diverges logarithmically whether for free ($Z=0$) or bound ($Z>0$) sums ([18], Table 1). We write $E_{\kappa }\left(a\right)=M_{\kappa }\left(a\right)$ when the intermediate states are free, $Z=0$, and $E_{\kappa }\left(a\right)=B_{\kappa }\left(a\right)$ for $Z>0$. The divergent sum $M\left(a\right)={\sum }_{\kappa }{M}_{\kappa }\left(a\right)$, the renormalization counter term, is usually interpreted as already forming part of the observed mass of the electron. Partial wave renormalization amounts to calculating the convergent sum $E\left(a\right)={\sum }_{\kappa }{R}_{\kappa }\left(a\right)$, where $R_{\kappa }\left(a\right)=B_{\kappa }\left(a\right)-M_{\kappa }\left(a\right)$ ([18], Table 2).

The use of CKG-spinors for self-energy has still to be investigated. The results of [18,19] for $Z>0$ used a Wick rotation, $k_0\to i\omega$, for the configuration space propagator in (45) similar to that used by Sapirstein and Cheng [38] in their work on vacuum polarization. The counter terms $M_{\kappa }\left(a\right)$ treated the bound state as a wave packet superposition of free states of the same symmetry ${\kappa }_a$: $|a\rangle ={\sum }_n|n{\kappa }_a,0\rangle \langle n{\kappa }_a,0|a\rangle$, so that the matrix elements in the numerator of (45) are proportional to the free electron vertex integrals [20]. The different algorithms all give similar accuracy.

6. Conclusions

QED was formulated [26,27,28,37], when computers were in their infancy and the most advanced computational aids were hand or electric calculators. Relativistic atomic structure programs like GRASP initially adapted algorithms in the 1960s ([3], Table 3) on lines proposed in Hartree’s 1957 lectures [1]. The symmetry structure of Dirac spinors on which GRASP is built [2,8], not yet fully understood in 1957, would not have been possible without the fundamental work of Racah [41,42,43] on quantum theory of angular momentum. Deeper understanding of the symmetry structure of spinor basis sets, developed over the last 30 years, along with advances in computer power and computational software, have motivated the fusion of computational methods from atomic and molecular physics and quantum electrodynamics discussed in this article.

Those familiar with quantum field theory will no doubt question why it is possible to perform calculations of radiative corrections in atoms without encountering the endemic infinities of perturbative QED. Weinberg ([28], Section 1.3) discusses the historic ‘problem of infinities’, noting that major progress on QED in the late 1940s required confidence in renormalization techniques. Similarly ‘Renormalization theory is a complicated conceptual system. …The renormalization procedure can be viewed as a technical device for circumventing—and discarding—the infinite results that occur in perturbative calculations in quantum field theory.’ ([27], p. 595) The reliance on regularization and renormalization has largely restricted serious attempts to calculate radiative corrections without expansion in powers of $Z\alpha$ to one- or two-electron systems. Spectroscopic applications need reliable data for many-electron atoms and ions beyond the first rows of the Periodic Table. GRASP and other programs rely on approximate formulae for self-energy and vacuum polarization based on hydrogenic calculations or effective potential approaches.

The complete reconciliation of conventional QED with the approach we have presented in this article and in our earlier work on the self energy [18,19] poses some open questions. There are clues: it is natural to use central field orbitals in atomic and molecular calculations, so that partial wave amplitudes at the vertices of Feynman diagrams are constrained by conservation of 3-momentum, angular momentum and parity as discussed in [20]. In the calculations of vacuum polarization reported here, and of the self-energy in [18,19], the Feynman diagrams are evaluated as partial wave expansions. Each term in these expansions is finite, and calculated without the use of any form of regularization. In the case of the self-energy, a partial wave mass-renormalization counter-term was introduced that renders the calculation of the physical effect finite, and in agreement with the results of conventional QED and experiment. We have demonstrated that the calculation of the higher-order terms in the vacuum polarization may also be achieved without regularization using the CK-basis, by subtracting corresponding terms in the partial wave expansions defined by Equations (43) and (44). This computational procedure is successful because the zero-potential term, and all higher-order terms involving an even number of external field interactions, vanish. This is a formal requirement of Furry’s theorem, but is not satisfied for all choices of basis set. The kinetically-balanced KG-basis set fails this fundamental test, as shown in Figure 2 and in results presented in [44]. The CKG-basis set passes this test by construction, so only physically meaningful contributions participate in the calculation of physical effects. Strong numerical cancellations are observed between the partial wave contributions in Table 2, which would be overwhelmed by any numerical failure of Furry’s theorem in the zero-potential term.

Our calculations of the leading-order self energy in [18,19] include external field contributions to all orders in $Z\alpha$, but the treatment of vacuum polarization discussed here includes only contributions of order $\alpha {\left(Z\alpha \right)}^3$ and higher. The missing part of the vacuum polarization of order $\alpha \left(Z\alpha \right)$ is represented by the expectation value of the Uehling potential. This well-understood effective potential incorporates charge renormalization and is conventionally derived using some form of regularization to tame the divergences that arise in its derivation. Our treatments of the self-energy and vacuum polarization would appear to be more symmetric and mutually complementary if the formally divergent part of the one-potential term, Equation (44), could be identified and eliminated within the partial wave expansion. It is an implicit assumption of our analysis, and of earlier work [21], that Equations (43) and (44) each contain the finite parts of the Uehling interaction and the divergent parts requiring renormalization. By the same reasoning, the sum of the differences for each partial wave results in a convergent series of finite, higher-order terms, but excludes the effects of the Uehling potential. A detailed analysis of the formulation of the vacuum polarization effect in configuration-space is presented in [17]. We are yet to develop this analysis to achieve an all-order treatment of vacuum polarization within the CKG-basis representation, though the calculation of the higher-order terms represents significant progress in that direction.

Since the calculation of bound-state matrix elements of the Uehling potential is readily achieved to high accuracy using any of the conventional methods of atomic physics, the combination of our methods for higher-order effects with a conventional calculation of the leading-order Uehling effect already offers a new computational treatment of QED effects. The unified approach that we have explored in this article points, however, to the development of a first-principles treatment mean-field, many-body and radiative effects within a common basis set and without the use of the QED machinery of regularization.

The adoption of a pragmatic point of view might suggest that the use of a patchwork of computational approaches in relativistic electronic structure theory currently achieves the accuracy required by experimental data in most practical applications. One should never forget, however, that a change of viewpoint, motivated by aesthetics, physical transparency or appeals to fundamental symmetries, may suggest new and better ways of doing things.