Calculations of QED Effects with the Dirac Green Function

Abstract

Modern spectroscopic experiments in few-electron atoms reached the level of precision at which an accurate description of quantum electrodynamics (QED) effects is mandatory. In many cases, theoretical treatment of QED effects need to be performed without any expansion in the nuclear binding strength parameter $Z\alpha$ (where Z is the nuclear charge number and $\alpha$ is the fine-structure constant). Such calculations involve multiple summations over the whole spectrum of the Dirac equation in the presence of the binding nuclear field, which can be evaluated in terms of the Dirac Green function. In this paper we describe the technique of numerical calculations of QED corrections with the Dirac Green function, developed in numerous investigations during the last two decades.

1. Introduction

Few-electron highly-charged ions are widely considered as important tools in testing quantum electrodynamics (QED) theory in the presence of the binding nuclear field [1,2,3]. Since the nuclear field in highly-charged ions is strong, its binding strength cannot be used as an expansion parameter and theoretical investigations of QED effects should be carried out to all orders in $Z\alpha$, where Z is the nuclear charge number and $\alpha$ is the fine structure constant. This is achieved by working in the so-called Furry picture, where the classical binding field of the nucleus is included into the zeroth-order approximation.

Interaction of the electron(s) bound in the field of the nucleus with the quantized radiation field gives rise to the QED effects, which are accounted for by an expansion in powers of $\alpha$. General expressions for individual QED corrections are derived within the dedicated methods, most notably, the adiabatic S-matrix formalism by Gell-Mann, Low and Sucher [4,5] and by the two-time Green function method by Shabaev [6].

The major difficulty encountered in calculations of QED corrections comes from the presence of infinite summations over the whole spectrum of the Dirac equation with the binding nuclear potential. These sums can be interpreted in terms of the so-called bound electron propagators, or the Dirac Green function.

Calculations of QED effects with the Dirac Green functions started in 1970th with computations of the one-loop self-energy [7,8,9] and vacuum-polarization [10,11]. Over the past years, the number and the complexity of QED calculations performed to all orders in the binding field has been increasing rapidly. These calculations have been successful not only in improving the achievable precision but also in extending the range of the studied effects, from the classical Lamb shift to the QED corrections to the hyperfine structure, the g factor, the transition amplitudes, the nuclear magnetic shielding, etc. This progress was due to not merely the increased computing speed and the availability of parallel computer resources, but also due to the development of new computational algorithms and methods.

With the present work we summarize the computational technique developed for calculations of various QED corrections with the bound-electron propagators, paying particular attention to the notoriously problematic diagrams with several propagators inside the radiative photon loop. This technique was developed in numerous calculations performed over the last two decades, notably, in refs. [12,13,14,15,16].

The relativistic units ($\hslash =c=m=1$) and the Heaviside charge units ($\alpha =e^2/4\pi$, $e<0$) will be used throughout this paper.

2. Dirac Green Function

The electron propagator $S(x_2,x_1)$ is standardly defined as the vacuum expectation value of the time-ordered product of the electron-positron field operators,

$$S(x_2,x_1)=-i\langle 0|T\mathsf{\Psi }(x_2)\overline{\mathsf{\Psi }}(x_1)|0\rangle ,$$

(1)

where T denotes the time-ordered product, $\overline{\mathsf{\Psi }}={\mathsf{\Psi }}^{†}{\gamma }^0$, and $\mathsf{\Psi }$ is the electron-positron field operator (see, e.g., ref. [1]),

$$\mathsf{\Psi }(x)=\sum _k{\hat{a}}_k{\phi }_k^{(+)}+\sum _k{\hat{b}}_k^{†}{\phi }_k^{(-)}.$$

(2)

Here, ${\hat{a}}^{†}$ (${\hat{b}}^{†}$) and $\hat{a}$ ($\hat{b}$) are the electron (positron) creation and annihilation operators, respectively; ${\phi }_k^{(\pm )}(x)={\psi }_k^{(\pm )}(\mathit{x})exp(-i{\epsilon }_k^{(\pm )}t)$ are single-particle electron (positron) states in the external field $A(x)$, and ${\psi }_k^{(\pm )}$ are the positive- and negative-energy eigenfunctions of the time-independent Dirac Hamiltonian ${\mathcal{H}}_D$,

$${\mathcal{H}}_D{\psi }_k(\mathit{x})\equiv \left[\mathit{\alpha }·(\mathit{p}-e\mathbf{A})+\beta m+e{A}^0\right]{\psi }_k(\mathit{x})={\epsilon }_k{\psi }_k(\mathit{x}),$$

(3)

where $\beta ={\gamma }^0$, $\mathit{\alpha }=\beta \mathit{\gamma }$, and $x=(t,\mathit{x})$ is a four-vector. Substituting Equation (2) into Equation (1), we get

$$\begin{array}{cc}\hfill S(x_2,x_1)=& -i\theta (t_2-t_1)\sum _k{\phi }_k^{(+)}(x_2)\overline{{\phi }_k^{(+)}}(x_1)\hfill \\ & +i\theta (t_1-t_2)\sum _k{\phi }_k^{(-)}(x_2),\overline{{\phi }_k^{(-)}}(x_1),\hfill \end{array}$$

(4)

where $\theta (t)$ is the Heaviside step function. This expression can be conveniently rewritten in an equivalent form

$$S(x_2,x_1)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }d\omega e^{-i\omega (t_2-t_1)}\sum _n\frac{{\psi }_n({\mathit{x}}_2){\overline{\psi }}_n({\mathit{x}}_1)}{\omega -{\epsilon }_n(1-i0)},$$

(5)

where the summation is carried out over both positive and negative energy states. Equivalence of these two representations for the electron propagator can be checked by performing the $\omega$ integration in Equation (5) by Cauchy’s theorem.

It can be easily shown (see, e.g., ref. [17]) that the electron propagator satisfies the differential equation

$$\left[i{\cancel{\partial }}_2-e\cancel{A}(x_2)-m\right]S(x_2,x_1)={\delta }^4(x_2-x_1),$$

(6)

where slashed symbols denote contractions with $\gamma$ matrices, $\cancel{\partial }={\gamma }_{\mu }{\partial }^{\mu }$ and $\cancel{A}={\gamma }_{\mu }{A}^{\mu }$. In the absence of the external field, this equation can be solved in a closed form. The result is the free-electron propagator,

$$S^{(0)}(x_2-x_1)=\int \frac{d^{4}p}{{(2\pi )}^4}{e}^{-ip·(x_2-x_1)}\frac{\cancel{p}+m}{p^2-m^2+i0},$$

(7)

where $p=(p^0,\mathit{p})$ is a four-vector.

Within the Feynman-diagram technique (see, e.g., ref. [6]), the integration over the time components of the arguments of the electron propagator is usually carried out in the general form, so that in practical calculations one deals with the Fourier transform of $S(x_2,x_1)$ with respect to the time variable $\tau =t_2-t_1$. The result is referred to as the Dirac Green function,

$$\begin{array}{cc}\hfill G(E,{\mathit{x}}_2,{\mathit{x}}_1)& ={\int }_{-\infty }^{\infty }d\tau e^{iE\tau }S(x_2,x_1){\gamma }^0=\sum _n\frac{{\psi }_n({\mathit{x}}_2){\psi }_n^{†}({\mathit{x}}_1)}{E-{\epsilon }_n(1-i0)}.\hfill \end{array}$$

(8)

From Equation (6), we deduce that $G(E,{\mathit{x}}_2,{\mathit{x}}_1)$ satisfies the differential equation

$$(E-{\mathcal{H}}_D)G(E,{\mathit{x}}_2,{\mathit{x}}_1)={\delta }^3({\mathit{x}}_2-{\mathit{x}}_1),$$

(9)

where ${\mathcal{H}}_D$ is the Dirac Hamiltonian, see Equation (3).

In this work we are interested in the Dirac Green function in a central field. In this case the angular structure of $G(E,{\mathit{x}}_2,{\mathit{x}}_1)$ follows from Equation (8) and from the angular dependence of the Dirac solutions [18],

$${\psi }_n(\mathit{x})=\left(\genfrac{}{}{0pt}{}{g_n(x){\chi }_{{\kappa }_n{\mu }_n}(\hat{\mathit{x}})}{i{f}_n(x){\chi }_{-{\kappa }_n{\mu }_n}(\hat{\mathit{x}})}\right),$$

(10)

where $g_n(r)$ and $f_n(r)$ are the upper and the lower radial components of the wave function, respectively, ${\chi }_{\kappa \mu }(\hat{\mathit{x}})$ is the spin-angular spinor, $\kappa$ is the relativistic angular-momentum quantum number, $\mu$ is the angular momentum projection, $x=\left|\mathit{x}\right|$, and $\hat{\mathit{x}}=\mathit{x}/\left|\mathit{x}\right|$. We thus obtain the standard partial-wave representation of the Dirac Green function [19,20]

$$G(E,{\mathit{x}}_2,{\mathit{x}}_1)=\sum _{\kappa }\left(\begin{array}{cc}{G}_{\kappa }^{11}(E,x_2,x_1){\pi }_{\kappa }^{++}({\hat{\mathit{x}}}_1,{\hat{\mathit{x}}}_2)& -i{G}_{\kappa }^{12}(E,x_2,x_1){\pi }_{\kappa }^{+-}({\hat{\mathit{x}}}_1,{\hat{\mathit{x}}}_2)\\ i{G}_{\kappa }^{21}(E,x_2,x_1){\pi }_{\kappa }^{-+}({\hat{\mathit{x}}}_1,{\hat{\mathit{x}}}_2)& G_{\kappa }^{22}(E,x_2,x_1){\pi }_{\kappa }^{--}({\hat{\mathit{x}}}_1,{\hat{\mathit{x}}}_2)\end{array}\right),$$

(11)

where ${\pi }_{\kappa }^{\pm \pm }({\hat{\mathit{x}}}_1,{\hat{\mathit{x}}}_2)={\sum }_{\mu }{\chi }_{\pm \kappa \mu }({\hat{\mathit{x}}}_1){\chi }_{\pm \kappa \mu }^{†}({\hat{\mathit{x}}}_2)$, and $G_{\kappa }^{ij}(E,x_2,x_1)$ are the radial components of the Dirac Green function.

For a static potential [$e{A}^0(\mathit{x})=V(x)$, $\mathbf{A}(x)=0$], Equation (9) in the matrix form reads

$$\left(\begin{array}{cc}E-m-V(x)& -(\mathit{\sigma }·\mathit{p})\\ -(\mathit{\sigma }·\mathit{p})& E+m-V(x)\end{array}\right)G(E,\mathit{x},{\mathit{x}}^{\prime })=\delta (\mathit{x}-{\mathit{x}}^{\prime })I,$$

(12)

where I is the 2 × 2 identity matrix. Substituting Equation (11) and using the identities

$$(\mathit{\sigma }·\mathit{p})f(x){\chi }_{\kappa \mu }(\hat{\mathit{x}})=i(\frac{\partial }{\partial x}+\frac{1+\kappa }{x})f(x){\chi }_{-\kappa \mu }(\hat{\mathit{x}}),$$

(13)

and

$$\delta (\mathit{x}-{\mathit{x}}^{\prime })=\frac{1}{x{x}^{\prime }}\delta (x-x^{\prime })\sum _{\kappa \mu }{\chi }_{\kappa \mu }(\hat{\mathit{x}}){\chi }_{\kappa \mu }^{†}({\hat{\mathit{x}}}^{\prime }),$$

(14)

we obtain the equation for the radial Dirac Green function,

$$\left(EI-h_{D,\kappa }\right)G_{\kappa }(E,x,x^{\prime })\equiv \left(\begin{array}{cc}{\displaystyle E-m-V(x)}& {\displaystyle \frac{d}{dx}-\frac{\kappa -1}{x}}\\ {\displaystyle -\frac{d}{dx}-\frac{\kappa +1}{x}}& {\displaystyle E+m-V(x)}\end{array}\right)G_{\kappa }(E,x,x^{\prime })=\frac{1}{x{x}^{\prime }}\delta (x-x^{\prime })I,$$

(15)

where $h_{D,\kappa }$ is the radial Dirac Hamiltonian and $G_{\kappa }$ is the 2 × 2 matrix of radial components of the Green function $G_{\kappa }^{ij}$, defined by Equation (11).

2.1. Representation in Terms of Regular and Irregular Solutions

The solution of an inhomogeneous differential Equation (15) can be constructed from the solutions of the corresponding homogeneous equation bounded at infinity (${\varphi }_{\kappa }^{\infty }$) and at origin (${\varphi }_{\kappa }^0$),

$$\begin{array}{cc}\hfill G_{\kappa }(E,x,x^{\prime })=& \frac{1}{{\mathsf{\Delta }}_{\kappa }(E)}\left[{\varphi }_{\kappa }^{\infty }(x){\varphi }_{\kappa }^{0^T}(x^{\prime })\theta (x-x^{\prime })+{\varphi }_{\kappa }^0(x){\varphi }_{\kappa }^{{\infty }^T}(x^{\prime })\theta (x^{\prime }-x)\right],\hfill \end{array}$$

(16)

where the superscript T denotes the transposition, ${\varphi }_{\kappa }^0$ and ${\varphi }_{\kappa }^{\infty }$ are the two-component solutions of the homogeneous radial Dirac equation, and ${\mathsf{\Delta }}_{\kappa }(E)$ is their Wronskian,

$$\begin{array}{c}\hfill {\mathsf{\Delta }}_{\kappa }(E)=x^2{\varphi }_{\kappa }^{0^T}(x)\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right){\varphi }_{\kappa }^{\infty }(x),\end{array}$$

(17)

which is independent of x. When the energy parameter E of the Green function is an eigenvalue of the Dirac Hamiltonian, the two solutions ${\varphi }_{\kappa }^0$ and ${\varphi }_{\kappa }^{\infty }$ coincide (up to a constant factor) and their Wronskian vanishes, ${\mathsf{\Delta }}_{\kappa }(E_n)=0$. This gives rise to poles of the Green function. The Green function has also branch points at $E=\pm m$, with cuts along the real axis for $\left|E\right|>m$, as will be discussed in more details below.

For the point-nucleus Coulomb potential [$V(x)=-Z\alpha /x$] the Equation (15) can be solved analytically [19] in terms of the Whittaker functions. The result is commonly referred to as the Dirac-Coulomb Green function. The radial Dirac-Coulomb Green function is represented by the form (16), with the functions ${\varphi }^0$ and ${\varphi }^{\infty }$ given by [8]

$${\varphi }_C^0(x)=\left(\genfrac{}{}{0pt}{}{{\varphi }_{C,+}^0(x)}{{\varphi }_{C,-}^0(x)}\right),{\varphi }_C^{\infty }(x)=\left(\genfrac{}{}{0pt}{}{{\varphi }_{C,+}^{\infty }(x)}{{\varphi }_{C,-}^{\infty }(x)}\right),$$

(18)

$$\begin{array}{ccc}\hfill {\varphi }_{C,\pm }^0(x)& =& \frac{\sqrt{1\pm \epsilon }}{x^{3/2}}\left[(\lambda -\nu )M_{\nu -(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.),\lambda }(2cx)\mp \left(\kappa -\frac{\alpha Z}{c}\right)M_{\nu +(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.),\lambda }(2cx)\right],\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\varphi }_{C,\pm }^{\infty }(x)& =& \frac{\sqrt{1\pm \epsilon }}{x^{3/2}}\left[\left(\kappa +\frac{\alpha Z}{c}\right)W_{\nu -(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.),\lambda }(2cx)\pm W_{\nu +(\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.),\lambda }(2cx)\right],\hfill \end{array}$$

(20)

and

$${\mathsf{\Delta }}_{C,\kappa }(E)=-4c^2\frac{\mathsf{\Gamma }(1+2\lambda )}{\mathsf{\Gamma }(\lambda -\nu )},$$

(21)

where $\epsilon =E/m$, $c=\sqrt{1-{\epsilon }^2}$, $\lambda =\sqrt{{\kappa }^2-{(\alpha Z)}^2}$, $\nu =Z\alpha \epsilon /c$, and $M_{\alpha ,\beta }$ and $W_{\alpha ,\beta }$ are the Whittaker functions of the first and the second kind [21], respectively. We mention the opposite sign of the present definition of the Green function as compared to the definition of refs. [1,8].

Zeros of the Wronskian (21) correspond to the bound-state energy levels, $\lambda -\nu =-n_r$ ($n_r=0,1,\dots$ is the radial quantum number), which yields the well-known formula for the Dirac bound energies,

$$E_{\kappa ,n_r}=m{\left[1+{\left(\frac{\alpha Z}{\lambda +n_r}\right)}^2\right]}^{-1/2}.$$

(22)

The cut structure of the Dirac-Coulomb Green function is defined by that of the square root $\sqrt{m^2-E^2}$. The square root function is defined to be positive in the gap $-m<E<m$ on the real E-axis. Outside of the gap, the sign of the square root is fixed by the condition $\mathrm{Re}(\sqrt{m^2-E^2})>0$. Special care should be taken evaluating the Green function for real energies $\left|E\right|>m$. Behaviour of the Green function on the real E axis is defined by the sign of the infinitesimal addition in the energy denominator of Equations (5) and (8). In our case the addition is negative and, therefore, the cut $E>m$ should be approached from the upper half of the E plane, and the cut $E<-m$ from the lower half. So, e.g., starting from the gap $-m<E<m$ and approaching the branch cut $E>m$ from the upper half-plane, we have the following prescription [22] for the analytical continuation of the square root: $\sqrt{m^2-E^2}\to -i\sqrt{E^2-m^2}$.

In the limit of $Z\to 0$, the Dirac-Coulomb Green function is reduced to the free Dirac Green function. The corresponding radial solutions are given by [8]

$$\begin{array}{ccc}\hfill {\varphi }_{F,\pm }^0(x)& =& \left(\begin{array}{c}1\\ i\frac{\kappa }{\left|\kappa \right|}\end{array}\right)\sqrt{1\pm \epsilon }{j}_{l_{\pm \kappa }}(icx),\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill {\varphi }_{F,\pm }^{\infty }(x)& =& \left(\begin{array}{c}1\\ i\frac{\kappa }{\left|\kappa \right|}\end{array}\right)\sqrt{1\pm \epsilon }{h}_{l_{\pm \kappa }}^{(1)}(icx),\hfill \end{array}$$

(24)

where $l_{\kappa }=|\kappa +1/2|-1/2$, $j(z)$ and $h^{(1)}(z)$ are spherical Bessel functions, and in $\left(\cdots \right)$ the upper value is chosen for the “+” component and the lower, for the “-” component. The Wronskian of the above solutions is ${\mathsf{\Delta }}_{F,\kappa }(E)=1/c$.

The numerical computation of the Whittaker functions required in calculations of QED corrections with the Dirac-Coulomb Green function was first tackled by Mohr in refs. [8,9] (see also the review [1]). His approach enabled an accurate computation of the Whittaker functions in a wide range of arguments, including high values of the relativistic angular parameter $\kappa$. A disadvantage of this numerical approach was that it required the extended-precision arithmetic to be used in a certain range of the arguments. A more economical variation of this approach was reported in ref. [12]. It allowed a computation of Whittaker functions within the standard double-precision arithmetics, for not very high partial waves ($\left|\kappa \right|<40$), which turned out to be sufficient for many practical applications.

The representation (16) can also be used in computations of the Dirac Green function for potentials other than the point-Coulomb potential. In particular, ref. [10] presented a numerical approach for computing the Dirac Green function for the potential induced by the nuclear charge distribution given by the shell nuclear model $\rho (r)\propto \delta (r-R)$. The computation of the Dirac Green function for the homogeneously charged nuclear model $\rho (r)\propto \theta (r-R)$ was reported in ref. [1].

In practical calculations, more realistic models of the nuclear-charge distribution are often required, first of all, the two-parameter Fermi distribution. A numerical approach for the computation of the Dirac Green function with the spherically-symmetric Fermi nuclear model was described in ref. [23]. This approach can be easily generated for the case of an arbitrary central potential approaching the Coulomb potential in the limit of $r\to \infty$, in particular, for a wide class of screened nuclear potentials.

2.2. Finite Basis Set Representations

Using the spectral representation of the Green function (8), we can express the radial Dirac Green function as

$$G_{\kappa }(E,x,x^{\prime })=\sum _n\frac{{\varphi }_{\kappa ,n}(x){\varphi }_{\kappa ,n}^T(x^{\prime })}{E-{\epsilon }_{\kappa ,n}},$$

(25)

where ${\varphi }_{\kappa ,n}$ are the two-component radial Dirac functions with energies ${\epsilon }_{\kappa ,n}$ satisfying the radial Dirac equation

$$h_{D,\kappa }{\varphi }_{\kappa ,n}(x)={\epsilon }_{\kappa ,n}{\varphi }_{\kappa ,n}(x).$$

(26)

The sum over n in Equation (25) should be understood as a summation over the discrete part of the spectrum and the integration over the positive and negative continuum parts of the spectrum.

A very useful approach to the numerical evaluation of the Dirac Green function is provided by the finite basis set method. In this method, the radial Dirac solutions are approximately represented by linear combinations of (a finite set of) two-component basis functions $u_i(x)$,

$$\begin{array}{c}\hfill {\varphi }_n(x)=\sum _{i=1}^N{c}_i{u}_i(x).\end{array}$$

(27)

within this representation, the solution of the radial Dirac Equation (26) is reduced to a generalized eigenvalue problem for the coefficients $c_i$,

$$\begin{array}{c}\hfill \frac{1}{2}\left[\langle u_i|h_D|u_k\rangle +\langle u_k|h_D|u_i\rangle \right]c_k=E\langle u_i|u_k\rangle c_k,\end{array}$$

(28)

where the summation over repeated indices is implied and $i,k=1\dots N$. This equation can be solved numerically by the standard methods of linear algebra, which yields the set of N eigenvectors and eigenvalues of the radial Dirac equation. After that, by using Equation (25) one obtains a finite basis-set representation of the radial Dirac Green function.

The choice of the basis functions $u_k$ can vary. One of the most successful implementations is delivered by the dual kinetic balance (DKB) basis [24] constructed with B-splines [25]. Within this method, the radial Dirac solutions are represented as

$$\begin{array}{cc}\hfill {\varphi }_{\kappa ,n}(x)=& \sum _{i=1}^{\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_i\left(\begin{array}{c}{B}_i(x)\\ {\displaystyle \frac{1}{2m}(\frac{d}{dx}+\frac{\kappa }{x})B_i(x)}\end{array}\right)+\sum _{i=1}^{\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{i+\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(\begin{array}{c}{\displaystyle \frac{1}{2m}(\frac{d}{dx}-\frac{\kappa }{x})B_i(x)}\\ B_i(x)\end{array}\right),\hfill \end{array}$$

(29)

where ${\left\{B_i(x)\right\}}_{i=1}^{\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ is the set of B-splines [25] on the interval $(0,R)$, where R is the cavity radius, chosen to be sufficiently large in order to have no influence on the calculated properties of the atom. The B-splines are chosen to vanish at $x=0$ and $x=R$, thus yielding the zero boundary conditions for the wave functions, $\varphi (0)=\varphi (R)=0$.

It needs to be stressed that the DKB ansatz (29) assumes that the potential in the Dirac equation is regular at $r\to 0$. This means that it can be used for solving the Dirac equation for an extended-nucleus potential, but not for the point-nucleus Coulomb potential. The advantages of the DKB basis is the absence of the so-called spurious states, the correct behaviour of the upper and lower radial components at $r\to 0$ and, as a consequence, an improved convergence of the calculated atomic properties with increase of the size of the basis set.

For the point nuclear model, one often uses the simpler ansatz of ref. [26],

$$\begin{array}{cc}\hfill {\varphi }_{\kappa ,n}(x)=& \sum _{i=1}^{\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_i\left(\begin{array}{c}{B}_i(x)\\ 0\end{array}\right)+\sum _{i=1}^{\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{c}_{i+\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\left(\begin{array}{c}0\\ B_i(x)\end{array}\right).\hfill \end{array}$$

(30)

It should be noted that the ansatz (30) leads to appearance of spurious eigenstates (highly oscillating eigenvectors with unphysical energies), as analytically proved in ref. [24]. In practical calculations, these spurious states do not cause significant problems (since their contributions to integrals is very small due to rapid oscillations), but their presence is manifested in a slower convergence of the calculated results as $N\to \infty$.

We also mention the space-discretization method for the solution of the Dirac equation [27,28], which can be regarded as a variant of the finite basis-set method with the basis constructed with $\delta$-functions. For practical calculations, the B-spline basis has the advantages of being more compact and consisting of continuous functions, while eigenvectors in the space-discretization method are defined on a grid only. However, the space-discretization method was successfully used in many calculations of QED effects by the Göteburg group, notably, in refs. [29,30,31]. Moreover, this method apparently yields a better convergence than the B-spline approach in calculations of the Wichmann-Kroll vacuum-polarization corrections [32].

2.3. Discussion

In Section 2.1 and Section 2.2 we described the two main representations of the bound electron propagator used in modern calculations of QED effects in atomic spectra. The first one is the representation in terms of the regular and irregular Dirac solutions (in what follows, the Green’s function approach) and the second is the finite basis set method. The other known representations of the Dirac-Coulomb Green function are not discussed in the present work, since they have not been proved useful in the calculations we consider here. In particular, the Sturmian expansion of the Dirac-Green function, widely used in the literature for the description of multiphoton processes (see, e.g., refs. [33,34,35,36]), does not seem to be useful for the calculations considered here. The main reasons are the numerical character of calculations and the lack of convergence of the Sturmian expansion when the energy argument is in the complex plane.

We now give a comparative discussion of the two main approaches. The basis-set method has several attractive features. The corresponding numerical routine is relatively simple, flexible and can easily incorporate any spherical-symmetric potential. Moreover, this method allows one to perform summations over a part of the Dirac spectrum (e.g., over the positive or the negative part only) and evaluate sums over spectrum with energy denominators different from the one in Equation (25).

Another attractive feature of the basis-set method is that it provides an approximation to the Green function which is a continuous function of the radial arguments at $x\approx x^{\prime }$. This is not the case for the exact Green function (16), whose components contain the discontinuous step function $\theta (x-x^{\prime })$ (which yields a $\delta$-function in Equation (15) after differentiation). This feature is often referred to as the radial ordering, since the exact Green function depends on $x_{<}$ and $x_{>}$, rather than just on x and $x^{\prime }$. This feature complicates the numerical evaluation of matrix elements, especially for higher-order diagrams with multiple radial integrations.

The basis-set method has also some important draw-backs as compared to the Green-function approach. It has an additional parameter, the number of basis functions N, and the final result should be investigated for stability when N is increased. In practice, the dependence on the basis size often sets a limitation on the accuracy of calculations. In addition, the number of partial waves (i.e., the maximal value of $\left|\kappa \right|$) included in the numerical evaluation is rather limited in the basis-set method. The typical number of partial waves employed in actual calculations with the basis-set method is ∼20, while in the Green-function approach it can be of order ${10}^4$ and more.

We conclude that the basis-set method has computational advantages for a restricted (but sufficiently broad) class of problems, where the partial-wave expansion is well converging and (or) the required numerical accuracy is not very high. The Green-function approach is preferable for problems where (i) high numerical accuracy is needed, (ii) large numerical cancellations occur, (iii) the partial-wave expansion does not converge rapidly, (iv) the contribution of high-energy intermediate electron states is enhanced, leading to slow convergence of the basis-set calculations with respect to N.

We now mention some of the calculations of QED corrections which used the above-mentioned methods for computing the bound electron propagators. Historically, the first was the Green-function approach elaborated, most notably, by Mohr in refs. [8,9]. This method was developed further in calculations of the one-loop self-energy [12,37,38,39,40], the self-energy correction to the hyperfine splitting and g factor [13,41,42], the screened QED corrections [43,44], and the QED corrections to the magnetic shielding [45]. The B-spline basis-set method was used in calculations of the two-photon exchange diagrams [46,47,48,49], the one-loop self-energy [50], the nuclear recoil [51,52,53], the screened QED corrections [54,55]. The space-discretization method was extensively applied by the Göteburg group in calculations of the first-order self-energy and vacuum polarization [32], two-photon exchange [29,31], the self-energy corrections to the bound-electron g factor [30], to the hyperfine structure [56], to the electron-electron interaction [57].

3. General Formulas

In the present work we will consider actual calculations of three contributions originating from the electron self-energy, specifically, the matrix elements of the self-energy operator, the magnetic vertex operator, and the double vertex operator, graphically represented in Figure 1, Figure 2 and Figure 3. The corresponding diagrams involve one, two, and three bound electron propagators in the radiative photon loop, respectively. Calculations of self-energy diagrams with one electron propagator started already in 1970th [7,8,9]. First calculations of the vertex diagrams with two electron propagators were performed in 1990th [43,58,59,60,61], whereas the double vertex diagrams have been tackled only relatively recently [45,62]. There have been no calculations of diagrams with more than three bound electron propagators in the radiative loop performed so far.

The matrix element of the one-loop self-energy operator depicted on Figure 1 yields the dominant contribution to the Lamb shift of the energy levels. It is given by

$$\begin{array}{c}\hfill \langle a|\Sigma ({\epsilon }_a)|a\rangle =\frac{i}{2\pi }{\int }_{-\infty }^{\infty }d\omega \sum _n\frac{\langle an|I(\omega )|na\rangle }{{\epsilon }_a-\omega -u{\epsilon }_n},\end{array}$$

(31)

where the summation over n is extended over the complete spectrum of the Dirac equation and $u\equiv 1-i0$ ensures the positions of the singularities of the Green function with respect to the integration contour. $I(\omega )$ is the operator of the electron-electron interaction, defined as

$$I(\omega ,{\mathit{r}}_1,{\mathit{r}}_2)=e^2{\alpha }_1^{\mu }{\alpha }_2^{\nu }{D}_{\mu \nu }(\omega ,{\mathit{r}}_{12}),$$

(32)

where ${\alpha }^{\mu }=(1,\mathit{\alpha })$ are the Dirac matrices, ${\mathit{r}}_{12}={\mathit{r}}_1-{\mathit{r}}_2$, and $D_{\mu \nu }(\omega ,{\mathit{r}}_{12})$ is the photon propagator. The photon propagator takes the simplest form in the Feynman gauge, where it is given by

$$D_{\mu \nu }(\omega ,{\mathit{r}}_{12})=g_{\mu \nu }\frac{e^{i\sqrt{{\omega }^2+i0}{r}_{12}}}{4\pi r_{12}},$$

(33)

with $r_{12}=\left|{\mathit{r}}_{12}\right|$.

The matrix element of the magnetic vertex operator depicted in Figure 2 is the most problematic part of the self-energy correction to the g factor [14]. The magnetic vertex operator, accompanied by the corresponding reducible part, is defined by its matrix elements as

$$\begin{array}{cc}\hfill \langle a\left|{\mathsf{\Lambda }}_{\mathrm{vr}}({\epsilon }_a)\right|a\rangle =& \frac{i}{2\pi }{\int }_{-\infty }^{\infty }d\omega \sum _{n_1{n}_2}\frac{\langle a{n}_2\left|I(\omega )\right|n_{1}a\rangle \left[\langle n_1|V_g|n_2\rangle -\langle n_1|n_2\rangle \langle a|V_g|a\rangle \right]}{({\epsilon }_a-\omega -u{\epsilon }_{n_1})({\epsilon }_a-\omega -u{\epsilon }_{n_2})},\hfill \end{array}$$

(34)

where $V_g$ is the effective magnetic operator responsible for the g factor [14], $V_g=(1/{\mu }_a){[\mathit{r}\times \mathit{\alpha }]}_z$, with ${\mu }_a$ being the angular momentum projection of the reference state a. We note that the scalar product $\langle n_1|n_2\rangle$ in Equation (34) can be trivially performed due to the orthogonality of the wave functions, $\langle n_1|n_2\rangle ={\delta }_{n_1{n}_2}$, but we find it convenient to keep it in the integral form.

The double-vertex operator matrix element shown in Figure 3 is the most problematic part of the self-energy correction to the nuclear shielding [45,63]. It is defined, together with the corresponding reducible parts, as

$$\begin{array}{cc}\hfill 〈a\left|{\mathsf{\Lambda }}_{\mathrm{dvr}}({\epsilon }_a)\right|a〉=2\frac{i}{2\pi }{\int }_{-\infty }^{\infty }d\omega & \{\sum _{n_1{n}_2{n}_3}\frac{\langle a{n}_3\left|I(\omega )\right|n_{1}a\rangle }{({\epsilon }_a-\omega -u{\epsilon }_{n_1})({\epsilon }_a-\omega -u{\epsilon }_{n_2})({\epsilon }_a-\omega -u{\epsilon }_{n_3})}\hfill \\ \hfill & \times [\langle n_1|V_g|n_2\rangle \langle n_2|V_{\mathrm{hfs}}|n_3\rangle -\langle n_1|n_2\rangle \langle n_2|V_{\mathrm{hfs}}|n_3\rangle \langle a|V_g|a\rangle \hfill \\ \hfill & -\langle n_1|V_g|n_2\rangle \langle n_2|n_3\rangle \langle a|V_{\mathrm{hfs}}|a\rangle +\langle n_1|n_2\rangle \langle n_2|n_3\rangle \langle a|V_g|a\rangle \langle a|V_{\mathrm{hfs}}|a\rangle ]\hfill \\ \hfill & -\sum _{{\mu }_{a^{\prime }}{n}_2}\frac{\langle a{a}^{\prime }\left|I(\omega )\right|a^{\prime }a\rangle }{{(-\omega +i0)}^2}\langle a|V_g|n_2\rangle \frac{1}{{\epsilon }_a-{\epsilon }_{n_2}}\langle n_2|V_{\mathrm{hfs}}|a\rangle \},\hfill \end{array}$$

(35)

where $V_{\mathrm{hfs}}$ is the effective magnetic operator responsible for the hyperfine interaction [15], $V_{\mathrm{hfs}}=(1/{\mu }_a){[\mathit{r}\times \mathit{\alpha }]}_z/r^3$, $a^{\prime }$ denotes the reference state a with a different momentum projection (${\mu }_{a^{\prime }}$), and the factor of 2 in the front accounts for two equivalent diagrams.

The above general formulas for the self-energy and magnetic-vertex matrix elements contain ultraviolet (UV) divergences. The standard approach to handle them [64] is to separate out one or two first terms of the expansion of the electron propagators in terms of the interaction with the binding nuclear field. In order to get UV-finite results, the self-energy operator needs a subtraction of the two first terms of the potential expansion,

$$\begin{array}{c}\hfill \Sigma ({\epsilon }_a)\to {\Sigma }^{(2+)}({\epsilon }_a)=\Sigma ({\epsilon }_a)-{\Sigma }^{(0)}({\epsilon }_a)-{\Sigma }^{(1)}({\epsilon }_a),\end{array}$$

(36)

whereas the vertex operator needs the subtraction of the first term only,

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_{\mathrm{vr}}({\epsilon }_a)\to {\mathsf{\Lambda }}_{\mathrm{vr}}^{(1+)}({\epsilon }_a)={\mathsf{\Lambda }}_{\mathrm{vr}}({\epsilon }_a)-{\mathsf{\Lambda }}_{\mathrm{vr}}^{(0)}({\epsilon }_a),\end{array}$$

(37)

where the superscript indicates the number of interactions with the binding field in the electron propagator(s) inside the radiative photon loop. The double-vertex operator ${\mathsf{\Lambda }}_{\mathrm{dvr}}$ contains three electron propagators inside the loop and thus is UV finite.

The separated terms containing zero and one interaction with the binding field (${\Sigma }^{(0)}$, ${\Sigma }^{(1)}$, ${\mathsf{\Lambda }}_{\mathrm{vr}}^{(0)}$) are regularized by using the dimensional regularization and calculated in momentum space. Their calculation does not involve bound electron propagators and thus is beyond the scope of the present paper; we refer the reader for the original investigations, ref. [12] for the self-energy matrix element, and ref. [14] for the magnetic vertex matrix element.

4. Angular Integration

The integration over the angular variables in the above formulas is conveniently carried out with help of the following representation of the matrix elements of the electron-electron interaction operator,

$$\langle ab\left|I(\omega )\right|cd\rangle =\alpha \sum _{L=L_{\min }}^{L_{\max }}{J}_L(abcd)R_L(\omega ,abcd),$$

(38)

where $J_L$ contains all the dependence on the angular momenta projections, $R_L$ are the radial integrals defined in Appendix A, and the summation over L goes from $L_{\min }=max(|j_a-j_c|,|j_b-j_d|)$ to $L_{\max }=min(j_a+j_c,j_b+j_d)$, with $j_n$ being the total angular momentum of the electron state n. The function $J_L$ is given by

$$\begin{array}{ccc}\hfill J_L(abcd)& =& \sum _{m_L}\frac{{(-1)}^{L-m_L+j_c-{\mu }_c+j_d-{\mu }_d}}{2L+1}{C}_{j_a{\mu }_a,j_c-{\mu }_c}^{L{m}_L}{C}_{j_d{\mu }_d,j_b-{\mu }_b}^{L{m}_L},\hfill \end{array}$$

(39)

where $C_{j_1{\mu }_1,j_2{\mu }_2}^{j\mu }$ denotes the Clebsch-Gordan coefficient and ${\mu }_n$ is the angular momentum projection of the electron state n.

Substituting Equation (38) into Equation (31) and performing the sum of two Clebsch-Gordan coefficients, we immediately obtain the result for the matrix element of the self-energy operator,

$$\begin{array}{cc}\hfill \langle a|\Sigma ({\epsilon }_a)|a\rangle =& \frac{i\alpha }{2\pi }{\int }_{-\infty }^{\infty }d\omega \sum _{n,L}\frac{{(-1)}^{j_a-j_n+L}}{2j_a+1}\frac{R_L(\omega ,anna)}{{\epsilon }_a-\omega -u{\epsilon }_n}.\hfill \end{array}$$

(40)

In order to perform the angular integrations in the magnetic vertex operator, we first apply the Wigner-Eckart theorem to the matrix element of the magnetic interaction $V_g$ (which is the rank-1 spherical tensor),

$$\begin{array}{c}\hfill \langle n_1|V_g|n_2\rangle =\frac{{(-1)}^{j_1-{\mu }_1}}{\sqrt{3}}{C}_{j_2-{\mu }_2,j_1{\mu }_1}^{10}(n_1\left|\right|V_g\left|\right|n_2),\end{array}$$

(41)

where $(.||.||.)$ denotes the reduced matrix element. Now we can perform the angular integration in the magnetic vertex matrix element as

$$\begin{array}{cc}\hfill \sum _{{\mu }_1{\mu }_2}\langle a{n}_2\left|I(\omega )\right|n_{1}a\rangle & \langle n_1|V_g|n_2\rangle =\sum _L{X}_L{R}_L(\omega ,a{n}_2{n}_{1}a)(n_1\left|\right|V_g\left|\right|n_2),\hfill \end{array}$$

(42)

where ${\mu }_1$ and ${\mu }_2$ are the angular momentum projections of the electron states $n_1$ and $n_2$, respectively, and $X_L$ are the angular coefficients defined by

$$\begin{array}{c}\hfill X_L=\sum _{{\mu }_1{\mu }_2}\frac{{(-1)}^{j_1-{\mu }_1}}{\sqrt{3}}{C}_{j_2-{\mu }_2,j_1{\mu }_1}^{10}{J}_L(a{n}_2{n}_{1}a).\end{array}$$

(43)

Performing the summation of three Clebsch-Gordan coefficients with help of formulas from ref. [65], we obtain

$$\begin{array}{c}\hfill X_L=\frac{{(-1)}^{j_1-j_2}{\mu }_a}{\sqrt{j_a(j_a+1)(2j_a+1)}}\left\{\begin{array}{ccc}{j}_1& j_2& 1\\ j_a& j_a& L\end{array}\right\},\end{array}$$

(44)

where $\{\dots \}$ denotes the $6j$-symbol.

Analogously, evaluating summations of four Clebsch-Gordan coefficients with help of formulas from ref. [65], we perform the angular integration for the double vertex matrix elements,

$$\begin{array}{cc}\hfill \sum _{{\mu }_1{\mu }_2{\mu }_3}& \langle a{n}_3\left|I(\omega )\right|n_{1}a\rangle \langle n_1|V_g|n_2\rangle \langle n_2|V_{\mathrm{hfs}}|n_3\rangle =\hfill \\ & \sum _L{Z}_L{R}_L(\omega ,a{n}_3{n}_{1}a)(n_1\left|\right|V_g\left|\right|n_2)(n_2\left|\right|V_{\mathrm{hfs}}\left|\right|n_3),\hfill \end{array}$$

(45)

where the angular coefficients $Z_L$ are

$$\begin{array}{c}\hfill Z_L=\sum _{j_k}{(-1)}^{j_k+j_a}{\left[C_{j_a{\mu }_a,10}^{j_k{\mu }_a}\right]}^2\left\{\begin{array}{ccc}1& j_1& j_2\\ j_a& L& j_3\\ j_k& j_a& 1\end{array}\right\},\end{array}$$

(46)

where $\{\dots \}$ denotes the $9j$-symbol. In practical calculations, summations over the angular momentum projections can be just as well carried out numerically.

The formulas above are written in terms of explicit summations over the Dirac spectrum, assuming the spectral representation of the radial Green function. In order to use the analytical representation of the Green function in terms of regular and irregular solutions, we would need to rewrite these formulas, identifying the components of the radial Green function, as

$$\begin{array}{c}\hfill \sum _n\frac{g_{\kappa ,n}(x)g_{\kappa ,n}(x^{\prime })}{E-{\epsilon }_n}\to G_{\kappa }^{11}(E,x,x^{\prime }),etc.\end{array}$$

(47)

This is possible but often leads to long and unnecessary cumbersome expressions, especially for complicated diagrams with multiple radial integrations. One can avoid this tedious work by introducing [66] the following formal representation for the radial Green function

$$G_{\kappa }(E,x,x^{\prime })={\psi }_{\kappa }(E,x){\tau }_{\kappa }^T(E,x^{\prime }),$$

(48)

where the two-component functions ${\psi }_{\kappa }$ and ${\tau }_{\kappa }$ depend on one radial argument only. The price to pay is that ${\psi }_{\kappa }$ and ${\tau }_{\kappa }$ have different forms depending on the ordering of the radial arguments x and $x^{\prime }$,

$${\psi }_{\kappa }(E,x)=\frac{1}{{\mathsf{\Delta }}_{\kappa }^{1/2}}\times \left\{\genfrac{}{}{0pt}{}{{\varphi }_{\kappa }^0(E,x),\mathrm{when}x<x^{\prime },}{{\varphi }_{\kappa }^{\infty }(E,x),\mathrm{when}x>x^{\prime },}\right.$$

(49)

and

$${\tau }_{\kappa }(E,x^{\prime })=\frac{1}{{\mathsf{\Delta }}_{\kappa }^{1/2}}\times \left\{\genfrac{}{}{0pt}{}{{\varphi }_{\kappa }^{\infty }(E,x^{\prime }),\mathrm{when}x<x^{\prime },}{{\varphi }_{\kappa }^0(E,x^{\prime }),\mathrm{when}x>x^{\prime }.}\right.$$

(50)

Here,${\varphi }_{\kappa }^0$ and ${\varphi }_{\kappa }^{\infty }$ are the two-component solutions of the radial Dirac equation bounded at the origin and infinity, respectively, and ${\mathsf{\Delta }}_{\kappa }$ is their Wronskian, see Equations (16) and (17). Employing the representation (48), we can immediately use formulas written via summations over the Dirac spectrum for calculations with the analytical representation of the Green function. The only complication is that the Green function is discontinuous when the two radial arguments are equal, $x=x^{\prime }$. This implies that radial integrations in different matrix elements cannot be performed independently. Their computation requires a special procedure, described in Section 7.

5. Choice of the Integration Contour

The formulas presented so far contained the integration over the virtual-photon energy $\omega$ performed along the real axis. This choice of the integration contour, however, is not favorable for numerical calculations, since the Dirac Green function is a highly oscillating function for large and real energy arguments and $x,x^{\prime }\to \infty$. It is advantageous to deform the integration contour to the region of large imaginary $\omega$ since the Dirac Green function acquires an exponentially damping factor in this case. Deforming the contour of integration, one should take care about poles and branch cuts of the integrand, however.

The analytical structure of the Dirac Green function is outlined in Section 2. The branch cuts of the photon propagator (15) are defined by the square root function, which should be understood as a limit of the regularized expression with a photon mass $\mu >0$,

$$\begin{array}{cc}\hfill \sqrt{{\omega }^2+i0}& \to \underset{\mu \to 0}{lim}\sqrt{{\omega }^2-{\mu }^2+i0}=\underset{\mu \to 0}{lim}\sqrt{\omega -\mu +i0}\sqrt{\omega +\mu -i0}.\hfill \end{array}$$

(51)

The photon propagator thus has two branch cuts starting from $\omega =\mu -i0$ and $\omega =-\mu +i0$. The analytical structure of the integrand of the self-energy matrix element is shown in Figure 4.

Figure 4 also presents the deformed contour of the $\omega$ integration, which we found to be optimal for most practical calculations. Specifically, the contour $C_{LH}$ consists of the low-energy part $C_L$ and the high-energy part $C_H$.

The low-energy part of the integration contour $C_L$ consists of two parts, $C_{L+}$ and $C_{L-}$, the first of which runs on the upper bank of the cut of the photon propagator and the second, on the lower bank and in the opposite direction. On the upper bank $\sqrt{{\omega }^2}=\omega$, whereas on the lower bank $\sqrt{{\omega }^2}=-\omega$. The integrands for $C_{L+}$ and $C_{L-}$ differ only by the sign of $\omega$ in the photon propagator (and the overall sign due to the opposite directions of the integration), thus allowing the following simplification,

$$e^{i\omega x_{12}}\to e^{i\omega x_{12}}-e^{-i\omega x_{12}}=2isin(\omega x_{12}).$$

The high-energy part $C_H$ is parallel to the imaginary axis and consists of two parts $C_{H-}=(\mathsf{\Delta }-i\infty ,\mathsf{\Delta }-iϵ)$ and $C_{H+}=(\mathsf{\Delta }+iϵ,\mathsf{\Delta }+i\infty )$. The integrands for $C_{H+}$ and $C_{H-}$ are typically complex conjugated, so that one can perform the integration over $C_{H+}$ only, take the real part of the result and multiply by two.

In the general case of an excited reference state, the low-energy part $C_L$ is bent in the complex plane, in order to avoid singularities coming from virtual bound states with energies ${\epsilon }_n<{\epsilon }_a$ in the electron propagator. Specifically, the contours $C_{L+}$ and $C_{L-}$ consist of 3 sections: $(0,{\delta }_{x,1}-i{\delta }_y)$, $({\delta }_{x,1}-i{\delta }_y,{\delta }_{x,2})$, and $({\delta }_{x,2},\mathsf{\Delta })$, as shown on Figure 4. The parameters of the contour ${\delta }_{x,1}$, ${\delta }_{x,2}$, ${\delta }_y$, and $\mathsf{\Delta }$ may be chosen differently. In our calculations, we used the following choice (assuming the reference state a to be an excited state): $\mathsf{\Delta }=Z\alpha {\epsilon }_a$; ${\delta }_{x,1}={\epsilon }_a-{\epsilon }_{1s}$; ${\delta }_{x,2}=2{\delta }_{x,1}$; ${\delta }_y={\delta }_{x,1}/2$. If the reference state a is the ground state, there is no need to bend the low-energy part of the contour in the complex plane (as there are no intermediate states with energy $0<{\epsilon }_n<{\epsilon }_a$); so we just integrate along the real axis (setting ${\delta }_y=0$).

We note that the described contour $C_{LH}$ resembles the contour used by P. Mohr in his calculations [8]. The difference is that he did not bend the low-energy part in the complex plane and used a different choice of the parameter $\mathsf{\Delta }$, $\mathsf{\Delta }={\epsilon }_a$.

Another choice of the $\omega$ integration contour frequently encountered in the literature (e.g., in refs. [50,67,68]) is the standard Wick rotation from the real into the imaginary axis, $\omega \to i\omega$. In this case the intermediate states with energy $0<{\epsilon }_n\le {\epsilon }_a$ lead to appearance of the pole terms, which need a special treatment. Apart for the pole contributions, small energy differences ${\epsilon }_a-{\epsilon }_n$ appear in the denominators of the electron propagators, leading to a rapidly varying structure of the integrand for small $\omega$ in this choice of the contour, which may lead to numerical difficulties.

6. Infrared Divergencies

In this section we address the infrared (IR) reference-state divergencies which appear in the QED corrections involving bound-electron propagators.

We start with pointing out that the bound-state QED corrections do not possess the standard free-QED IR divergences which arise when the four-momentum p of the intermediate electron states approaches the mass shell, $\rho =(m^2-p^2)/m^2=(m^2-p_0^2+{\mathit{p}}^2)/m^2\to 0$. For the bound-state QED corrections, the intermediate electron states are always off mass shell ($p_0={\epsilon }_a<m\Rightarrow \rho >0$) and the would-be IR divergences are cut off by the binding energy of the reference state. However, the bound-state QED corrections often contain IR divergences of a different kind, also known as the reference-state divergences. They appear when two or more denominators in the electron propagators vanish at $\omega \to 0$. Specifically, IR divergences arise in the magnetic vertex operator (34) when $n_1=n_2=a$ and in the double vertex operator (35) when $n_1=n_3=a$.

The general approach to the treatment of the IR divergencies is to separate out the divergent contributions, regularize them by introducing a finite photon mass $\mu$ in the photon propagator, evaluate the integral over $\omega$ analytically, and separate out the $\mu$-dependent divergent terms. The divergent terms should of course cancel out when all relevant contributions are summed together. The evaluation of the IR-divergent integrals with the finite photon mass is illustrated in Appendix B.

Using formulas from Appendix B, the magnetic vertex matrix element (34) can be transformed to a form that is explicitly free from any IR divergences,

$$\begin{array}{cc}\hfill \langle a\left|{\mathsf{\Lambda }}_{\mathrm{vr}}({\epsilon }_a)\right|a\rangle =& \frac{i}{2\pi }{\int }_{-\infty }^{\infty }d\omega [\sum _{n_1{n}_2}\frac{\langle a{n}_2\left|I(\omega )\right|n_{1}a\rangle \left[\langle n_1|V_g|n_2\rangle -{\delta }_{n_1{n}_2}\langle a|V_g|a\rangle \right]}{({\epsilon }_a-\omega -u{\epsilon }_{n_1})({\epsilon }_a-\omega -u{\epsilon }_{n_2})}\hfill \\ & -\sum _{{\mu }_{a^{\prime }}{\mu }_{a^{″}}}\frac{\langle a{a}^{″}\left|I(\omega )\right|a^{\prime }a\rangle \left[\langle a^{\prime }|V_g|a^{″}\rangle -{\delta }_{a^{\prime }{a}^{″}}\langle a|V_g|a\rangle \right]}{{(-\omega +i0)}^2}]\hfill \\ & +\frac{\alpha }{\pi }\sum _{{\mu }_{a^{\prime }}{\mu }_{a^{″}}}\langle a{a}^{″}|{\alpha }_{\mu }{\alpha }^{\mu }ln{x}_{12}|a^{\prime }a\rangle \left[\langle a^{\prime }|V_g|a^{″}\rangle -{\delta }_{a^{\prime }{a}^{″}}\langle a|V_g|a\rangle \right],\hfill \end{array}$$

(52)

where $a^{\prime }$ and $a^{″}$ denote the reference state a with a different angular momentum projection (${\mu }_{a^{\prime }}$ and ${\mu }_{a^{″}}$, respectively).

For the double-vertex matrix element (35) the situation is somewhat more complicated because there are two types of divergences, the one $\propto 1/\mu$ coming from three vanishing denominators ($n_1=n_2=n_3=a$) and the other $\propto ln\mu$ coming from two vanishing denominators ($n_1=n_3=a\ne n_2$). Still, the divergences in the double-vertex matrix can be handled with help of formulas from Appendix B analogously to that for the magnetic vertex case.

There exists also a more economic method of handling IR divergences in actual calculations. It relies on the fact that the matrix elements (34) and (35) are defined so that they are overall IR finite, i.e., have a well-defined limit at $\mu \to 0$. This means that they can be numerically evaluated with the zero photon mass. As long as the $\omega$ integration is performed after all parts of the integrand are combined together, the integrand will have a smooth small-$\omega$ behaviour when integrated along the low-energy part of the integration contour $C_{LH}$. The would-be IR divergences will be cancelled numerically at a given $\omega$ between different parts of the integrand. It is easy to check that for the magnetic vertex matrix element, both methods give the identical numerical results. For the double vertex matrix element, the numerical treatment of IR divergences was used in the calculation of the diamagnetic shielding in refs. [45,63].

It should be mentioned that vanishing denominators in the electron propagators could arise not only from the intermediate states $n=a$, but also from the intermediate states having the same energy but the opposite parity as the reference state (e.g., $n=2p_{1/2}$ for $a=2s$ for the point nuclear model). Such intermediate states do not cause IR divergences, since the radial matrix element in the numerator vanishes due to the orthogonality of the wave functions, as can be seen from formulas in Appendix B.

A different approach for handling the IR divergencies was used by the Notre-Dame group [67,69,70]. In their works, numerical calculations were performed with an explicit regularization parameter $\delta$ shifting the position of the reference-state energy, ${\epsilon }_a\to {\epsilon }_a(1-\delta )$; the numerical limit of $\delta \to 0$ was then performed in the end of the calculations.

7. Radial Integration

In actual calculations it is very important to find an efficient way to perform multiple radial integrations. The number of radial integrations is two for the self-energy matrix element, three for the magnetic vertex, and four for the double-vertex matrix element. In what follows we will assume that the analytical representation of the Dirac Green function is used, since in the basis-set representation, the radial integrations do not cause particular difficulties.

We now formulate a general numerical approach suitable for carrying out multiple radial integrations, first introduced in the context of the two-loop self-energy in ref. [66]. We start with the simplest case of the self-energy matrix element, in which the radial integral is two-dimensional. The two-dimensional radial integrals can be schematically represented to be a linear combination of terms of the following structure

$${\int }_0^{\infty }d{x}_1{\int }_0^{\infty }d{x}_{2}H(x_1)I(x_{<})L(x_{>})M(x_2),$$

(53)

where $x_{>}=max(x_1,x_2)$ and $x_{<}=min(x_1,x_2)$ and H, I, L, M are some functions of the specified arguments. It is important that the integrand can be represented as a product of functions that depend on one radial argument only, some of those being $x_{<}$ and $x_{>}$. In particular, $I(x_{<})$ involves ${\varphi }^0(x_{<})$ from the Dirac Green function and $j_l(\omega x_{<})$ from the photon propagator, and $L(x_{>})$ involves ${\varphi }^{\infty }(x_{>})$ and $h_l^{(1)}(\omega x_{>})$. It is clear that if we store all functions on a suitably chosen radial grid, it should be possible to compute the integral (53) just by summing up the pre-stored data.

In order to do this, we introduce a 3-dimensional radial grid $\left\{r_{i,j,k}\right\}$ on the interval $(0,r_{\max })$ as follows. First, we fill the elements of the first layer, $r_{i,0,0}$ with $i=0,\dots ,N_i$, which coarsely span the whole interval, e.g., as

$$r_{i,0,0}=r_0\frac{1-t^2}{t^2},$$

(54)

where t is uniformly distributed on the interval $(t_{\min },1)$ and $t_{\min }\approx 0$ is defined by the cavity radius $r_{\max }$. Next, we introduce a finer grid of the second layer. Specifically, on each interval $(r_{i,0,0},r_{i+1,0,0})$ we introduce the set of the Gauss-Legendre abscissae ${\left\{r_{i,j,0}\right\}}_{j=1}^{N_j}$. We see that in order to perform the outer radial integration, it is sufficient to know the integrand on the grid $\left\{r_{i,j,0}\right\}$.

In order to perform the inner radial integral, we have to split the integration interval at the point $x_2=x_1$, since it is the discontinuity point of the integrand. We achieve this by introducing a yet finer grid of the third layer, $\left\{r_{i,j,k}\right\}$. Specifically, for fixed values of i and j, the set ${\left\{r_{i,j,k}\right\}}_{k=1}^{N_k}$ represents the Gauss-Legendre abscissae on the interval $(r_{i,j,0},r_{i,j+1,0})$ if $j<N_j$ and on the interval $(r_{i,j,0},r_{i+1,0,0})$ if $j=N_j$. Now, for each point $r_{i,j,0}$ of the outer radial integration, we can perform the inner integral splitted into parts, $(0,r_{i,j,0})$ and $(r_{i,j,0},r_{\max })$.

We conclude that when all functions in the integrand of Equation (53) are stored on the radial grid $\left\{r_{i,j,k}\right\}$, the two-dimensional numerical integration can be carried out by just summing up the pre-stored numerical values. The described procedure can be easily generalized for integrals of higher dimensions. So, for a computation of a four-dimensional radial integral, we need a 5-dimensional grid $\left\{r_{i,j,k,l,m}\right\}$ introduced in the same way as $\left\{r_{i,j,k}\right\}$.

An additional complication arises from the fact that the regular Dirac solution ${\varphi }_{\kappa }^0(\epsilon ,r)$ in the Dirac Green function has typically the exponentially-growing behaviour for large values of the radial argument and complex values of $\epsilon$, whereas the irregular solutions ${\varphi }_{\kappa }^{\infty }(\epsilon ,r)$ is exponentially decreasing in this region. In order to avoid numerical overflow and underflow, we store the “normalized” solutions ${\tilde{\varphi }}^0$ and ${\tilde{\varphi }}^{\infty }$, with the approximate large-r and small-r behaviour pulled out,

$$\begin{array}{cc}\hfill {\varphi }_{\kappa }^0(\epsilon ,r)& =r^{\left|\kappa \right|}{e}^{cr}{\tilde{\varphi }}_{\kappa }^0(\epsilon ,r),\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill {\varphi }_{\kappa }^{\infty }(\epsilon ,r)& =r^{-\left|\kappa \right|}{e}^{-cr}{\tilde{\varphi }}_{\kappa }^{\infty }(\epsilon ,r),\hfill \end{array}$$

(56)

where $c=\sqrt{1-{(\epsilon /m)}^2}$. When ${\varphi }_{\kappa }^0(\epsilon ,r)$ and ${\varphi }_{\kappa }^{\infty }(\epsilon ,r)$ multiply together in the Dirac Green function, the result is usually in the range accessible in the standard double-precision (8-byte) arithmetics. A similar normalization is required also for the regular $j_l$ and irregular $h_l^{(1)}$ Bessel solutions, originating from the photon propagator. With these precautions, we are able to perform calculations completely within the standard double-precision arithmetics typically for $\kappa \le 50$. For $\kappa \le 100$, it is usually possible to use the quadruple-precision arithmetics for computation of Dirac (${\varphi }_{\kappa }^0$, ${\varphi }_{\kappa }^{\infty }$) and Bessel ($j_l$, $h_l^{(1)}$) solutions but the double-precision arithmetics for the radial integrations. For even higher values of $\kappa$, use of the extended-precision arithmetics becomes unavoidable [71].

8. Magnetically-Perturbed Green Function

The computation of radial integrations in diagrams with various kind of potentials can be significantly accelerated by introducing the first-order perturbations of the Green function by this potentials. Such an approach was used long ago by Gyulassy in his evaluation of the vacuum polarization [72]. More recently, similar algorithms were used in calculations of various self-energy corrections (in particular, in refs. [23,73,74]).

In this section we describe the computation of the Dirac Green function perturbed by a magnetic potential $V_g$, which will be referred to as the magnetically-perturbed Green function. Specifically, we are interested in the radial part of the magnetically-perturbed Green function, defined as

$$\begin{array}{c}\hfill {\mathcal{G}}_{{\kappa }_1{\kappa }_2}(x_1,x_3)={\int }_0^{\infty }d{x}_2{x}_2^2{G}_{{\kappa }_1}(x_1,x_2)V_g(x_2)G_{{\kappa }_2}(x_2,x_3),\end{array}$$

(57)

where $V_g(x)=x{\sigma }_x$ is the radial part of the magnetic potential $V_g(\mathit{x})$. Using the representation of the radial Green functions in terms of the regular and irregular Dirac solutions, see Equation (16), we obtain the following expressions for the magnetically-perturbed Green function. For $x_1<x_3$, we get

$$\begin{array}{cc}\hfill {\mathcal{G}}_{{\kappa }_1{\kappa }_2}& (x_1,x_3)={\varphi }_{{\kappa }_1}^{\infty }(x_1){\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{00}(x_1){\varphi }_{{\kappa }_2}^{{\infty }^T}(x_3)\hfill \\ & +{\varphi }_{{\kappa }_1}^0(x_1){\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{\infty \infty }(x_3){\varphi }_{{\kappa }_2}^{0^T}(x_3),\hfill \\ & +{\varphi }_{{\kappa }_1}^0(x_1)\left[{\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{\infty 0}(x_3)-{\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{\infty 0}(x_1)\right]{\varphi }_{{\kappa }_2}^{{\infty }^T}(x_3),\hfill \end{array}$$

(58)

whereas for $x_1>x_3$,

$$\begin{array}{cc}\hfill {\mathcal{G}}_{{\kappa }_1{\kappa }_2}& (x_1,x_3)={\varphi }_{{\kappa }_1}^{\infty }(x_1){\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{00}(x_3){\varphi }_{{\kappa }_2}^{{\infty }^T}(x_3)\hfill \\ & +{\varphi }_{{\kappa }_1}^0(x_1){\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{\infty \infty }(x_1){\varphi }_{{\kappa }_2}^{0^T}(x_3),\hfill \\ & +{\varphi }_{{\kappa }_1}^{\infty }(x_1)\left[{\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{0\infty }(x_1)-{\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{0,\infty }(x_3)\right]{\varphi }_{{\kappa }_2}^{0^T}(x_3),\hfill \end{array}$$

(59)

where for simplicity we assumed that ${\varphi }_{\kappa }^0$ and ${\varphi }_{\kappa }^{\infty }$ are normalized so that their Wronskian is unity and the functions ${\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}$ are defined by the integrals

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{00}(x)=& {\int }_0^{x}d{x}_2{x}_2^2{\varphi }_{{\kappa }_1}^{0^T}(x_2)V_g(x_2){\varphi }_{{\kappa }_2}^0(x_2),\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{0\infty }(x)=& {\int }_0^{x}d{x}_2{x}_2^2{\varphi }_{{\kappa }_1}^{0^T}(x_2)V_g(x_2){\varphi }_{{\kappa }_2}^{\infty }(x_2),\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{\infty 0}(x)=& {\int }_0^{x}d{x}_2{x}_2^2{\varphi }_{{\kappa }_1}^{{\infty }^T}(x_2)V_g(x_2){\varphi }_{{\kappa }_2}^0(x_2),\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}^{\infty \infty }(x)=& {\int }_x^{\infty }d{x}_2{x}_2^2{\varphi }_{{\kappa }_1}^{{\infty }^T}(x_2)V_g(x_2){\varphi }_{{\kappa }_2}^{\infty }(x_2).\hfill \end{array}$$

(63)

We observe that after storing the functions ${\varphi }_{\kappa }(x)$ and ${\mathsf{\Phi }}_{{\kappa }_1,{\kappa }_2}(x)$ on a radial grid, we are able to construct the magnetically-perturbed Green function ${\mathcal{G}}_{{\kappa }_1{\kappa }_2}(x_1,x_2)$ for any radial arguments needed in our computation. The integral functions ${\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}(x)$ are evaluated by numerical integration with help of Gauss-Legendre quadratures. It is important that only one integral over $(0,\infty )$ needs to be evaluated (for a given value of the energy argument) in order to store ${\mathsf{\Phi }}_{{\kappa }_1{\kappa }_2}(x)$ on the whole radial grid. Analogously to the case of the plain Green function, all manipulations with the regular and irregular solutions need to be carried out after normalizing them according to Equations (55) and (56) in order to prevent numerical overflow.

9. Numerical Calculations

In this section we demonstrate the technique described in previous sections with three examples of actual calculations. The first one is the calculation of the one-loop self-energy correction to the Lamb shift of a hydrogen-like ion. In

Table 1. Numerical results for the one-loop self-energy correction for the $2s$ state of hydrogen-like calcium ($Z=20$), for the point nucleus, in terms of the standard scaled function $F(Z\alpha )=\delta E/[(\alpha /\pi ){(Z\alpha )}^4/n^3]$, where $\delta E$ is a contribution to the energy in relativistic units. $S_l$ denotes the sum of partial-wave expansion ${\sum }_{\left|\kappa \right|=1}^l$; $\delta S_l$ is the increment with respect to the previous line.

	l	${\mathit{S}}_{\mathit{l}}$	$\mathit{\delta }{\mathit{S}}_{\mathit{l}}$
$\langle {\Sigma }^{(2+)}\rangle$	1	82.268 19
	2	85.541 56	3.273 37
	3	86.515 41	0.973 85
	4	86.967 68	0.452 26
	5	87.223 70	0.256 02
	10	87.675 13	0.451 44
	15	87.790 89	0.115 75
	20	87.836 31	0.045 42
	30	87.869 99	0.033 68
	40	87.881 51	0.011 52
	50	87.886 57	0.005 06
	60	87.889 17	0.002 61
	∞	87.894 34 (26)	0.005 17 (26)
$\langle {\Sigma }^{(0+1)}\rangle$		−84.387 704
Total		3.506 64 (26)
P. J. Mohr [75]		3.506 648 (2)
Refs. [76,77]		3.506 647 (5)

we present numerical results for the one-loop self-energy correction to the Lamb shift of the $2s$ state of hydrogen-like calcium ($Z=20$), for the point nuclear model.

The many-potential part $\langle {\Sigma }^{(2+)}\rangle$ defined by Equation (36) is calculated in coordinate space by the method described in the present work. Specifically, the one-potential Green function was calculated by the method described in Section 8. (Alternatively, it can also be calculated as a derivative over the nuclear charge Z, as described in ref. [12].) For the radial integration, we used a four-dimensional grid $\left\{r_{i,j,k,l}\right\}$ constructed as discussed in Section 7 with the number of integration points $(N_i,N_j,N_k,N_l)=(15,10,6,6)$. The $\omega$ integration is carried out along the contour $C_{LH}$ introduced in Section 5 using the Gauss-Legendre quadratures, after mapping of the integration intervals to the range $(0,1)$. The summation over the partial waves was extended up to $\left|\kappa \right|=60$, with the remaining tail of the expansion estimated by a least-square fitting in polynomials in $1/\left|\kappa \right|$. The remaining zero- and one-potential part $\langle {\Sigma }^{(0+1)}\rangle$ is calculated in momentum space. Their computation is relatively simple and can be performed up to essentially arbitrary accuracy. This part is not discussed here; we refer the reader to the original work [12].

As follows from Table 1, the uncertainty of the final numerical result for the self-energy correction comes exclusively from the truncation of the partial-wave expansion. It can be seen that despite the inclusion of 60 partial waves, the resulting accuracy is significantly lower than that of the best literature values. There are two ways described in the literature that allow to achieve a better numerical precision. One method was developed originally by Mohr [8,9,75] and extended by Jentschura and Mohr [38,39]. This method involves a summation of many thousands of partial waves and usage of extended-precision arithmetics in order to obtain very accurate numerical results. Another method was developed in ref. [76]. It involves an additional subtraction in ${\Sigma }^{(2+)}$ which greatly accelerates the convergence of the partial-wave expansion and allows one to obtain accurate numerical results with just 20–30 partial waves. Both these methods are difficult to extend for computations of more complicated diagrams, unfortunately.

Table 2. Numerical results for the self-energy correction to the g factor of the $2s$ state of hydrogen-like calcium ($Z=20$), for the point nucleus (in units of ${10}^{-6}$).

	l	${\mathit{S}}_{\mathit{l}}$	$\mathit{\delta }{\mathit{S}}_{\mathit{l}}$
$\langle {\mathsf{\Lambda }}_{\mathrm{vr}}^{(2+)}\rangle$	1	36.130 52
	2	17.563 35	−18.567 17
	3	14.605 25	−2.958 10
	4	13.586 86	−1.018 39
	5	13.115 22	−0.471 64
	10	12.489 16	−0.626 06
	15	12.379 38	−0.109 78
	20	12.343 92	−0.035 46
	25	12.328 78	−0.015 15
	30	12.321 11	−0.007 67
	35	12.316 76	−0.004 35
	∞	12.306 43 (50)	−0.010 33 (50)
$\langle {\mathsf{\Lambda }}_{\mathrm{vr}}^{(0+1)}\rangle$		2237.914 11
$\langle {\mathsf{\Lambda }}_{\mathrm{ir}}\rangle$		75.453 02
Total		2325.673 56 (50)
Ref. [14]		2325.674 (5)

presents our numerical results for the self-energy correction to the g factor of the $2s$ state of hydrogen-like calcium ($Z=20$), for the point nuclear model. The many-potential part $\langle {\mathsf{\Lambda }}_{\mathrm{vr}}^{(2+)}\rangle$ is calculated in coordinate space by the method described in the present work. It is important that we calculate the magnetic vertex after subtracting two first terms of its potential expansion, not just one as in Equation (37). This is done in order to accelerate the convergence of the partial-wave expansion of the matrix element, following refs. [14,30]. The subtracted part $\langle {\mathsf{\Lambda }}_{\mathrm{vr}}^{(0+1)}\rangle$ is calculated in momentum space as described in ref. [14]. The irreducible part $\langle {\mathsf{\Lambda }}_{\mathrm{ir}}\rangle$ is expressed as a non-diagonal matrix element of the self-energy operator; its numerical values were taken from ref. [14]. The total result presented in Table 2 is in good agreement with the previous value obtained in ref. [14]. Its numerical uncertainty comes exclusively from the truncation of the partial-wave expansion. Even more accurate results can be achieved if one extends the partial-wave expansion further, as was done for the $1s$ state in ref. [71], but it requires significant efforts and intensive usage of extended-precision arithmetics.

In

Table 3. Numerical results for the self-energy correction to the nuclear magnetic shielding constant $\sigma$ of the $1s$ state of hydrogen-like calcium ($Z=20$), for the point nucleus, in units of the scaled function $D(Z\alpha )=\delta \sigma /\left[{\alpha }^2{(Z\alpha )}^3\right]$, where $\delta \sigma$ is a contribution to the shielding constant.

	l	${\mathit{S}}_{\mathit{l}}$	$\mathit{\delta }{\mathit{S}}_{\mathit{l}}$
$\langle {\mathsf{\Lambda }}_{\mathrm{dvr}}\rangle$	1	−3.409 2
	2	−5.550 9	−2.141 7
	3	−6.559 8	−1.008 9
	4	−7.111 6	−0.551 7
	5	−7.438 5	−0.327 0
	10	−7.941 8	−0.503 2
	15	−7.986 1	−0.044 3
	20	−7.968 3	0.017 7
	25	−7.945 0	0.023 3
	30	−7.925 5	0.019 5
	35	−7.910 4	0.015 1
	∞	−7.846 7 (32)	0.063 7 (32)
$\langle {\mathsf{\Lambda }}_{\mathrm{der}}\rangle$		7.782 4
$\langle {\mathsf{\Lambda }}_{\mathrm{vr},\mathrm{Zee}}\rangle$		1.760 7
$\langle {\mathsf{\Lambda }}_{\mathrm{vr},\mathrm{hfs}}\rangle$		−0.404 9
$\langle {\mathsf{\Lambda }}_{\mathrm{po}}\rangle$		−2.217 0
Total		−0.925 5 (32)

we present numerical results for the self-energy correction to the diamagnetic shielding constant of the $1s$ state of hydrogen-like calcium ($Z=20$), for the point nuclear model. The many-potential part $\langle {\mathsf{\Lambda }}_{\mathrm{dvr}}\rangle$ is calculated in coordinate space by the method described in the present work. We observe a slow convergence of the partial-wave expansion of the results presented in the table. It can probably be accelerated by separating out the leading term of the potential expansion (i.e., the contribution of the free propagators) and calculating it in the momentum space, but this has not been accomplished so far. The other contributions to the shielding constant are defined in ref. [63]; the corresponding numerical results are taken from that work. Again, we observe that the dominant uncertainty of the final result comes from the truncation of the partial-wave expansion.

10. Summary

In this paper we described the technique used in modern calculations of QED corrections with the bound-electron propagators, including the notoriously problematic diagrams with several propagators inside the radiative photon loop. The bound-electron propagators are described by the Green function of the Dirac equation with the binding nuclear potential. We considered two most widely used ways to represent the Dirac Green function, the representation via the regular and irregular Dirac solutions and the finite basis set representation. These representations are applicable for a wide range of binding potentials, including the case of the nuclear field modified by a spherically-symmetric screening potential caused by the presence of other electrons in the atom.

We demonstrated that the dominant uncertainty of the obtained results usually comes from the truncation of the partial-wave expansion. Further extension of the partial-wave expansion is possible but often associated with large technical difficulties. In view of this, it is important to look for ways to accelerate convergence of the partial-wave expansion. This was accomplished for the one-loop self-energy in ref. [76], for the self-energy correction to the g factor in refs. [14,41], and for the self-energy correction to the hyperfine splitting in refs. [15,42]. Unfortunately, all these methods turned out to be problem-specific, i.e., they do not allow straightforward extensions to more complicated corrections. It would be thus of great importance to find a more universal approach to improve the convergence of the partial-wave expansion in such calculations.