Calculations of QED effects with the Dirac Green function
aa r X i v : . [ phy s i c s . a t o m - ph ] M a y Calculations of QED effects with the Dirac Green function
Vladimir A. Yerokhin and Anna V. Maiorova
Center for Advanced Studies, Peter the Great St. Petersburg Polytechnic University,Polytekhnicheskaya 29, 195251 St. Petersburg, Russia
Modern spectroscopic experiments in few-electron atoms reached the level of precision at whichan 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 nuclearbinding strength parameter Zα (where Z is the nuclear charge number and α is the fine-structureconstant). Such calculations involve multiple summations over the whole spectrum of the Diracequation in the presence of the binding nuclear field, which can be evaluated in terms of the DiracGreen function. In this paper we describe the technique of numerical calculations of QED correctionswith the Dirac Green function, developed in numerous investigations during the last two decades. I. INTRODUCTION
Few-electron highly-charged ions are widely consideredas important tools in testing quantum electrodynamics(QED) theory in the presence of the binding nuclear field[1–3]. Since the nuclear field in highly-charged ions isstrong, its binding strength cannot be used as an expan-sion parameter and theoretical investigations of QED ef-fects should be carried out to all orders in Zα , where Z is the nuclear charge number and α is the fine structureconstant. This is achieved by working in the so-calledFurry picture, where the classical binding field of the nu-cleus is included into the zeroth-order approximation.Interaction of the electron(s) bound in the field of thenucleus with the quantized radiation field gives rise tothe QED effects, which are accounted for by an expan-sion in powers of α . General expressions for individualQED corrections are derived within the dedicated meth-ods, most notably, the adiabatic S -matrix formalism byGell-Mann, Low and Sucher [4, 5] and by the two-timeGreen function method by Shabaev [6].The major difficulty encountered in calculations ofQED corrections comes from the presence of infinite sum-mations over the whole spectrum of the Dirac equationwith the binding nuclear potential. These sums can beinterpreted in terms of the so-called bound electron prop-agators, or the Dirac Green function.Calculations of QED effects with the Dirac Green func-tions started in 1970th with computations of the one-loop self-energy [7–9] and vacuum-polarization [10, 11].Over the past years, the number and the complexity ofQED calculations performed to all orders in the bind-ing field has been increasing rapidly. These calculationshave been successful not only in improving the achievableprecision but also in extending the range of the studiedeffects, from the classical Lamb shift to the QED correc-tions to the hyperfine structure, the g factor, the tran-sition amplitudes, the nuclear magnetic shielding, etc.This progress was due to not merely the increased com-puting speed and the availability of parallel computerresources, but also due to the development of new com-putational algorithms and methods.With the present work we summarize the computa- tional technique developed for calculations of variousQED corrections with the bound-electron propagators,paying particular attention to the notoriously problem-atic diagrams with several propagators inside the radia-tive photon loop. This technique was developed in nu-merous calculations performed over the last two decades,notably, in Refs. [12–16].The relativistic units (¯ h = c = m = 1) and the Heav-iside charge units ( α = e / π , e <
0) will be usedthroughout this paper.
II. DIRAC GREEN FUNCTION
The electron propagator S ( x , x ) is standardly de-fined as the vacuum expectation value of the time-orderedproduct of the electron-positron field operators, S ( x , x ) = − i h | T Ψ( x ) Ψ( x ) | i , (1)where T denotes the time-ordered product, Ψ = Ψ † γ ,and Ψ is the electron-positron field operator (see, e.g.,Ref. [1]), Ψ( x ) = X k ˆ a k ϕ (+) k + X k ˆ b † k ϕ ( − ) k . (2)Here, ˆ a † (ˆ b † ) and ˆ a (ˆ b ) are the electron (positron)creation and annihilation operators, respectively; ϕ ( ± ) k ( x ) = ψ ( ± ) k ( x ) exp( − iε ( ± ) k t ) are single-particleelectron (positron) states in the external field A ( x ), and ψ ( ± ) k are the positive- and negative-energy eigenfunctionsof the time-independent Dirac Hamiltonian H D , H D ψ k ( x ) ≡ (cid:2) α · ( p − e A ) + βm + eA (cid:3) ψ k ( x ) = ε k ψ k ( x ) , (3)where β = γ , α = β γ , and x = ( t, x ) is a four-vector.Substituting Eq. (2) into Eq. (1), we get S ( x , x ) = − iθ ( t − t ) X k ϕ (+) k ( x ) ϕ (+) k ( x )+ iθ ( t − t ) X k ϕ ( − ) k ( x ) , ϕ ( − ) k ( x ) , (4)where θ ( t ) is the Heaviside step function. This expressioncan be conveniently rewritten an equivalent form S ( x , x ) = 12 π Z ∞−∞ dω e − iω ( t − t ) X n ψ n ( x ) ψ n ( x ) ω − ε n (1 − i , (5)where the summation is carried out over both positiveand negative energy states. Equivalence of these two rep-resentations for the electron propagator can be checkedby performing the ω integration in Eq. (5) by Cauchy’stheorem.It can be easily shown (see, e.g. , Ref. [17]) that theelectron propagator satisfies the differential equation (cid:2) i /∂ − e /A ( x ) − m (cid:3) S ( x , x ) = δ ( x − x ) , (6)where slashed symbols denote contractions with γ ma-trices, /∂ = γ µ ∂ µ and /A = γ µ A µ . In the absence of theexternal field, this equation can be solved in a closedform. The result is the free-electron propagator, S (0) ( x − x ) = Z d p (2 π ) e − ip · ( x − x ) /p + mp − m + i , (7)where p = ( p , p ) is a four-vector.Within the Feynman-diagram technique (see, e.g.,Ref. [6]), the integration over the time components of thearguments of the electron propagator is usually carriedout in the general form, so that in practical calculations one deals with the Fourier transform of S ( x , x ) withrespect to the time variable τ = t − t . The result isreferred to as the Dirac Green function, G ( E, x , x ) = Z ∞−∞ dτ e iEτ S ( x , x ) γ = X n ψ n ( x ) ψ † n ( x ) E − ε n (1 − i . (8)From Eq. (6), we deduce that G ( E, x , x ) satisfies thedifferential equation( E − H D ) G ( E, x , x ) = δ ( x − x ) , (9)where H D is the Dirac Hamiltonian, see Eq. (3).In this work we are interested in the Dirac Green func-tion in a central field. In this case the angular structureof G ( E, x , x ) follows from Eq. (8) and from the angulardependence of the Dirac solutions [18], ψ n ( x ) = (cid:18) g n ( x ) χ κ n µ n ( ˆ x ) if n ( x ) χ − κ n µ n ( ˆ x ) (cid:19) , (10)where g n ( r ) and f n ( r ) are the upper and the lower radialcomponents of the wave function, respectively, χ κµ ( ˆ x )is the spin-angular spinor, κ is the relativistic angular-momentum quantum number, µ is the angular momen-tum projection, x = | x | , and ˆ x = x / | x | . We thusobtain the standard partial-wave representation of theDirac Green function [19, 20] G ( E, x , x ) = X κ (cid:18) G κ ( E, x , x ) π ++ κ ( ˆ x , ˆ x ) − i G κ ( E, x , x ) π + − κ ( ˆ x , ˆ x ) i G κ ( E, x , x ) π − + κ ( ˆ x , ˆ x ) G κ ( E, x , x ) π −− κ ( ˆ x , ˆ x ) (cid:19) , (11)where π ±± κ ( ˆ x , ˆ x ) = P µ χ ± κµ ( ˆ x ) χ †± κµ ( ˆ x ), and G ijκ ( E, x , x ) are the radial components of the Dirac Greenfunction.For a static potential [ eA ( x ) = V ( x ), A ( x ) = 0], Eq. (9) in the matrix form reads (cid:18) E − m − V ( x ) − ( σ · p ) − ( σ · p ) E + m − V ( x ) (cid:19) G ( E, x , x ′ ) = δ ( x − x ′ ) I , (12)where I is the 2 × σ · p ) f ( x ) χ κµ ( ˆ x ) = i (cid:16) ∂∂x + 1 + κx (cid:17) f ( x ) χ − κµ ( ˆ x ) , (13)and δ ( x − x ′ ) = 1 xx ′ δ ( x − x ′ ) X κµ χ κµ ( ˆ x ) χ † κµ ( ˆ x ′ ) , (14)we obtain the equation for the radial Dirac Green function, (cid:0) E I − h D,κ (cid:1) G κ ( E, x, x ′ ) ≡ E − m − V ( x ) ddx − κ − x − ddx − κ + 1 x E + m − V ( x ) G κ ( E, x, x ′ ) = 1 xx ′ δ ( x − x ′ ) I , (15)where h D,κ is the radial Dirac Hamiltonian and G κ is the 2 × G ijκ , defined by Eq. (11). A. Representation in terms of regular andirregular solutions
The solution of an inhomogeneous differential equa-tion (15) can be constructed from the solutions of the corresponding homogeneous equation bounded at infin-ity ( φ ∞ κ ) and at origin ( φ κ ), G κ ( E, x, x ′ ) = 1∆ κ ( E ) (cid:20) φ ∞ κ ( x ) φ T κ ( x ′ ) θ ( x − x ′ )+ φ κ ( x ) φ ∞ T κ ( x ′ ) θ ( x ′ − x ) (cid:21) , (16)where the subscript T denotes the transposition, φ κ and φ ∞ κ are the two-component solutions of the homogeneousradial Dirac equation, and ∆ κ ( E ) is their Wronskian,∆ κ ( E ) = x φ T κ ( x ) (cid:18) − (cid:19) φ ∞ κ ( x ) , (17)which is independent on x . When the energy param-eter E of the Green function is an eigenvalue of theDirac Hamiltonian, the two solutions φ κ and φ ∞ κ coin-cide (up to a constant factor) and their Wronskian van-ishes, ∆ κ ( E n ) = 0. This gives rise to poles of the Greenfunction. The Green function has also branch points at E = ± m , with cuts along the real axis for | E | > m , aswill be discussed in more details below.For the point-nucleus Coulomb potential [ V ( x ) = − Zα/x ] the equation (15) can be solved analytically [19]in terms of the Whittaker functions. The result is com-monly referred to as the Dirac-Coulomb Green function.The radial Dirac-Coulomb Green function is representedby the form (16), with the functions φ and φ ∞ given by[8] φ C ( x ) = φ C, + ( x ) φ C, − ( x ) ! , φ ∞ C ( x ) = φ ∞ C, + ( x ) φ ∞ C, − ( x ) ! , (18) φ C, ± ( x ) = √ ± εx / (cid:20) ( λ − ν ) M ν − ( / ) , λ (2 cx ) ∓ (cid:18) κ − αZc (cid:19) M ν +( / ) , λ (2 cx ) (cid:21) , (19) φ ∞ C, ± ( x ) = √ ± εx / (cid:20) (cid:18) κ + αZc (cid:19) W ν − ( / ) , λ (2 cx ) ± W ν +( / ) , λ (2 cx ) (cid:21) , (20)and ∆ C,κ ( E ) = − c Γ(1 + 2 λ )Γ( λ − ν ) , (21)where ε = E/m , c = √ − ε , λ = p κ − ( αZ ) , ν = Zα ε/c , and M α,β and W α,β are the Whittaker functionsof the first and the second kind [21], respectively. Wemention the opposite sign of the present definition of theGreen function as compared to the definition of Refs. [1,8].Zeros of the Wronskian (21) correspond to the bound-state energy levels, λ − ν = − n r ( n r = 0 , , . . . is the radial quantum number), which yields the well-knownformula for the Dirac bound energies, E κ,n r = m " (cid:18) αZλ + n r (cid:19) − / . (22)The cut structure of the Dirac-Coulomb Green func-tion is defined by that of the square root √ m − E .The square root function is defined to be positive in thegap − m < E < m on the real E -axis. Outside of thegap, the sign of the square root is fixed by the condi-tion Re( √ m − E ) >
0. Special care should be takenevaluating the Green function for real energies | E | > m .Behaviour of the Green function on the real E axis isdefined by the sign of the infinitesimal addition in theenergy denominator of Eqs. (5) and (8). In our case theaddition is negative and, therefore, the cut E > m shouldbe approached from the upper half of the E plane, andthe cut E < − m from the lower half. So, e.g., startingfrom the gap − m < E < m and approaching the branchcut E > m from the upper half-plane, we have the fol-lowing prescription [22] for the analytical continuation ofthe square root: √ m − E → − i √ E − m .In the limit of Z →
0, the Dirac-Coulomb Green func-tion is reduced to the free Dirac Green function. Thecorresponding radial solutions are given by [8] φ F, ± ( x ) = (cid:18) i κ | κ | (cid:19) √ ± ε j l ± κ ( icx ) , (23) φ ∞ F, ± ( x ) = (cid:18) i κ | κ | (cid:19) √ ± ε h (1) l ± κ ( icx ) , (24)where l κ = | κ + 1 / | − / j ( z ) and h (1) ( z ) are sphericalBessel functions, and in ( · · · ) the upper value is cho-sen for the “+” component and the lower, for the “-”component. The Wronskian of the above solutions is∆ F,κ ( E ) = 1 /c .The numerical computation of the Whittaker func-tions required in calculations of QED corrections withthe Dirac-Coulomb Green function was first tackled byMohr in Refs. [8, 9] (see also the review [1]). His ap-proach enabled an accurate computation of the Whit-taker functions in a wide range of arguments, includinghigh values of the relativistic angular parameter κ . Adisadvantage of this numerical approach was that it re-quired the extended-precision arithmetic to be used ina certain range of the arguments. A more economicalvariation of this approach was reported in Ref. [12]. Itallowed a computation of Whittaker functions within thestandard double-precision arithmetics, for not very highpartial waves ( | κ | < ∼ (cid:1) FIG. 1: The one-loop self-energy correction. The doubleline represents the electron propagating in the bindingfield of the nucleus. The wavy line denotes the virtualphoton.FIG. 2: The magnetic-vertex self-energy correction.The wavy line terminated by a cross denotes theinteraction with an external magnetic field. ρ ( r ) ∝ δ ( r − R ). The computation of the Dirac Greenfunction for the homogeneously charged nuclear model ρ ( r ) ∝ θ ( r − R ) was reported in ref. [1].In practical calculations, more realistic models of thenuclear-charge distribution are often required, first of all,the two-parameter Fermi distribution. A numerical ap-proach for the computation of the Dirac Green functionwith the spherically-symmetric Fermi nuclear model wasdescribed in ref. [23]. This approach can be easily gen-erated for the case of an arbitrary central potential ap-proaching the Coulomb potential in the limit of r → ∞ ,in particular, for a wide class of screened nuclear poten-tials. B. Finite basis set representations
Using the spectral representation of the Green function(8), we can represent the radial Dirac Green function as G κ ( E, x, x ′ ) = X n φ κ,n ( x ) φ Tκ,n ( x ′ ) E − ε κ,n , (25)where φ κ,n are the two-component radial Dirac functionswith energies ε κ,n satisfying the radial Dirac equation h D,κ φ κ,n ( x ) = ε κ,n φ κ,n ( x ) . (26)The sum over n in Eq. (25) should be understood as asummation over the discrete part of the spectrum and FIG. 3: The double-vertex self-energy correction. Thewavy line terminated by a triangle denotes thehyperfine interaction.the integration over the positive and negative continuumparts of the spectrum.A very useful approach to the numerical evaluation ofthe Dirac Green function is provided by the finite basisset method. In this method, the radial Dirac solutionsare approximately represented by linear combinations of(a finite set of) two-component basis functions u i ( x ) , φ n ( x ) = N X i =1 c i u i ( x ) . (27)Within this representation, the solution of the radialDirac equation (26) is reduced to a generalized eigenvalueproblem for the coefficients c i ,12 (cid:2) h u i | h D | u k i + h u k | h D | u i i (cid:3) c k = E h u i | u k i c k , (28)where the summation over repeated indices is implied and i, k = 1 . . . N . This equation can be solved numericallyby the standard methods of linear algebra, which yieldsthe set of N eigenvectors and eigenvalues of the radialDirac equation. After that, by using Eq. (25) one obtainsa finite basis-set representation of the radial Dirac Greenfunction.The choice of the basis function u k can vary. Oneof the most successful implementations is delivered bythe dual kinetic balance (DKB) basis [24] constructedwith B -splines [25]. Within this method, the radial Diracsolutions are represented as φ κ,n ( x ) = N / X i =1 c i B i ( x )12 m (cid:18) ddx + κx (cid:19) B i ( x ) + N / X i =1 c i + N / m (cid:18) ddx − κx (cid:19) B i ( x ) B i ( x ) , (29)where { B i ( x ) } N / i =1 is the set of B -splines [25] on the inter-val (0 , R ), where R is the cavity radius, chosen to be suf-ficiently large in order to have no influence on the calcu-lated properties of the atom. The B -splines are chosen tovanish at x = 0 and x = R , thus yielding the zero bound-ary conditions for the wave functions, φ (0) = φ ( R ) = 0.It needs to be stressed that the DKB anzatz (29) as-sumes that the potential in the Dirac equation is regularat r →
0. This means that it can be used for solving theDirac equation for an extended-nucleus potential, but not for the point-nucleus Coulomb potential. The advantagesof the DKB basis is the absence of the so-called spuriousstates, the correct behaviour of the upper and lower ra-dial components at r → φ κ,n ( x ) = N / X i =1 c i (cid:18) B i ( x )0 (cid:19) + N / X i =1 c i + N / (cid:18) B i ( x ) (cid:19) . (30)It should be noted that the anzatz (30) leads to appear-ance of spurious eigenstates (highly oscillating eigenvec-tors with unphysical energies), as analytically proved inRef. [24]. In practical calculations, these spurious statesdo not cause significant problems (since their contribu-tions to integrals is very small due to rapid oscillations),but their presence is manifested in a slower convergenceof the calculated results as N → ∞ .We also mention the space-discretization method forthe solution of the Dirac equation [27, 28], which canbe regarded as a variant of the finite basis-set methodwith the basis constructed with δ -functions. For practicalcalculations, the B -spline basis has the advantages of be-ing more compact and consisting of continuous functions,while eigenvectors in the space-discretization method aredefined on a grid only. However, the space-discretizationmethod was successfully used in many calculations ofQED effects by G¨oteburg group, notably, in Refs. [29–31]. Moreover, this method apparently yields a betterconvergence than the B -spline approach in calculationsof the Wichmann-Kroll vacuum-polarization corrections[32]. C. Discussion
In Sec. II A and II B we described the two main rep-resentations of the bound electron propagator used inmodern calculations of QED effects in atomic spectra.The first one is the representation in terms of the reg-ular and irregular Dirac solutions (in what follows, theGreen’s function approach) and the second is the finitebasis set method. The other known representations ofthe Dirac-Coulomb Green function are not discussed inthe present work, since they have not been proved usefulin the calculations we consider here. In particular, theSturmian expansion of the Dirac-Green function, widely used in the literature for the description of multiphotonprocesses (see, e.g. , Refs. [33–36]), does not seem to beuseful for the calculations considered here. The main rea-sons are the numerical character of calculations and thelack of convergence of the Sturmian expansion when theenergy argument is in the complex plane.We now give a comparative discussion of the two mainapproaches. The basis-set method has several attrac-tive features. The corresponding numerical routine isrelatively simple, flexible and can easily incorporate anyspherical-symmetric potential. Moreover, this methodallows one to perform summations over a part of theDirac spectrum (e.g., over the positive or the negativepart only) and evaluate sums over spectrum with energydenominators different from the one in Eq. (25).Another attractive feature of the basis-set method isthat it provides an approximation to the Green functionwhich is a continuous function of the radial arguments at x ≈ x ′ . This is not the case for the exact Green function(16), whose components contain the discontinuous stepfunction θ ( x − x ′ ) (which yields a δ -function in Eq. (15)after differentiation). This feature is often referred to asthe radial ordering , since the exact Green function de-pends on x < and x > , rather than just on x and x ′ . Thisfeature complicates the numerical evaluation of matrixelements, especially for higher-order diagrams with mul-tiple radial integrations.The basis-set method has also some important draw-backs as compared to the Green-function approach. Ithas an additional parameter, the number of basis func-tions N , and the final result should be investigated forstability when N is increased. In practice, the depen-dence on the basis size often sets a limitation on theaccuracy of calculations. In addition, the number of par-tial waves (i.e., the maximal value of | κ | ) included inthe numerical evaluation is rather limited in the basis-setmethod. The typical number of partial waves employedin actual calculations with the basis-set method is ∼ and more.We conclude that the basis-set method has computa-tional advantages for a restricted (but sufficiently broad)class of problems, where the partial-wave expansion iswell converging and (or) the required numerical accuracyis not very high. The Green-function approach is prefer-able for problems where (i) high numerical accuracy isneeded, (ii) large numerical cancellations occur, (iii) thepartial-wave expansion does not converge rapidly, (iv) thecontribution of high-energy intermediate electron statesis enhanced, leading to slow convergence of the basis-setcalculations with respect to N .We now mention some of the calculations of QED cor-rections which used the above-mentioned methods forcomputing the bound electron propagators. Historically,the first was the Green-function approach elaborated,most notably, by Mohr in refs. [8, 9]. This method wasdeveloped further in calculations of the one-loop self-energy [12, 37–40], the self-energy correction to the hy-perfine splitting and g factor [13, 41, 42], the screenedQED corrections [43, 44], and the QED corrections to themagnetic shielding [45]. The B -spline basis-set methodwas used in calculations of the two-photon exchange dia-grams [46–49], the one-loop self-energy [50], the nuclearrecoil [51–53], the screened QED corrections [54, 55].The space-discretization method was extensively appliedby the G¨oteburg group in calculations of the first-orderself-energy and vacuum polarization [32], two-photon ex-change [29, 31], the self-energy corrections to the bound-electron g factor [30], to the hyperfine structure [56],to the electron-electron interaction [57]. III. GENERAL FORMULAS
In the present work we will consider actual calculationsof three contributions originating from the electron self-energy, specifically, the matrix elements of the self-energyoperator, the magnetic vertex operator, and the doublevertex operator, graphically represented in Figures 1–3.The corresponding diagrams involve one, two, and threebound electron propagators in the radiative photon loop,respectively. Calculations of self-energy diagrams withone electron propagator started already in 1970th [7–9].First calculations of the vertex diagrams with two elec-tron propagators were performed in 1990th [43, 58–61], whereas the double vertex diagrams have been tackledonly relatively recently [45, 62]. There have been no cal-culations of diagrams with more than three bound elec-tron propagators in the radiative loop performed so far.The matrix element of the one-loop self-energy opera-tor depicted on Fig. 1 yields the dominant contributionto the Lamb shift of the energy levels. It is given by h a | Σ( ε a ) | a i = i π Z ∞−∞ dω X n h an | I ( ω ) | na i ε a − ω − uε n , (31)where the summation over n is extended over the com-plete spectrum of the Dirac equation and u ≡ − i I ( ω ) isthe operator of the electron-electron interaction, definedas I ( ω, r , r ) = e α µ α ν D µν ( ω, r ) , (32)where α µ = (1 , α ) are the Dirac matrices, r = r − r ,and D µν ( ω, r ) is the photon propagator. The pho-ton propagator takes the simplest form in the Feynmangauge, where it is given by D µν ( ω, r ) = g µν e i √ ω + i r πr , (33)with r = | r | .The matrix element of the magnetic vertex operator depicted in Fig. 2 is the most problematic part of the self-energycorrection to the g factor [14]. The magnetic vertex operator, accompanied by the corresponding reducible part, isdefined by its matrix elements as h a | Λ vr ( ε a ) | a i = i π Z ∞−∞ dω X n n h an | I ( ω ) | n a i (cid:2) h n | V g | n i − h n | n i h a | V g | a i (cid:3) ( ε a − ω − u ε n )( ε a − ω − u ε n ) , (34)where V g is the effective magnetic operator responsible for the g factor [14], V g = (1 /µ a ) [ r × α ] z , with µ a being theangular momentum projection of the reference state a . We note that the scalar product h n | n i in Eq. (34) can betrivially performed due to the orthogonality of the wave functions, h n | n i = δ n n , but we find it convenient to keepit in the integral form.The double-vertex operator matrix element shown in Fig. 3 is the most problematic part of the self-energy correctionto the nuclear shielding [45, 63]. It is defined, together with the corresponding reducible parts, as (cid:10) a (cid:12)(cid:12) Λ dvr ( ε a ) (cid:12)(cid:12) a (cid:11) = 2 i π Z ∞−∞ dω ( X n n n h an | I ( ω ) | n a i ( ε a − ω − u ε n )( ε a − ω − u ε n )( ε a − ω − u ε n ) × h h n | V g | n i h n | V hfs | n i − h n | n i h n | V hfs | n i h a | V g | a i− h n | V g | n i h n | n i h a | V hfs | a i + h n | n i h n | n i h a | V g | a i h a | V hfs | a i i − X µ a ′ n h aa ′ | I ( ω ) | a ′ a i ( − ω + i h a | V g | n i ε a − ε n h n | V hfs | a i ) , (35)where V hfs is the effective magnetic operator responsible for the hyperfine interaction [15], V hfs = (1 /µ a ) [ r × α ] z /r , a ′ denotes the reference state a with a different angular momentum projection ( µ a ′ ), and the factor of 2 in the frontaccounts 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 handlethem [64] is to separate out one or two first terms of theexpansion of the electron propagators in terms of the in-teraction with the binding nuclear field. In order to getUV-finite results, the self-energy operator needs a sub-traction of the two first terms of the potential expansion,Σ( ε a ) → Σ (2+) ( ε a ) = Σ( ε a ) − Σ (0) ( ε a ) − Σ (1) ( ε a ) , (36)whereas the vertex operator needs the subtraction of thefirst term only,Λ vr ( ε a ) → Λ (1+)vr ( ε a ) = Λ vr ( ε a ) − Λ (0)vr ( ε a ) , (37)where the superscript indicates the number of interac-tions with the binding field in the electron propagator(s)inside the radiative photon loop. The double-vertex op-erator Λ dvr contains three electron propagators inside theloop and thus is UV finite.The separated terms containing zero and one interac-tion with the binding field (Σ (0) , Σ (1) , Λ (0)vr ) are regu-larized by using the dimensional regularization and cal-culated in momentum space. Their calculation does notinvolve bound electron propagators and thus is beyondthe scope of the present paper; we refer the reader for theoriginal investigations, Ref. [12] for the self-energy ma-trix element, and Ref. [14] for the magnetic vertex matrixelement. IV. ANGULAR INTEGRATION
The integration over the angular variables in the aboveformulas is conveniently carried out with help of thefollowing representation of the matrix elements of theelectron-electron interaction operator, h ab | I ( ω ) | cd i = α L max X L = L min J L ( abcd ) R L ( ω, abcd ) , (38)where J L contains all the dependence on the angular mo-menta projections, R L are the radial integrals definedin 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 theelectron state n . The function J L is given by J L ( abcd ) = X m L ( − L − m L + j c − µ c + j d − µ d L + 1 × C Lm L j a µ a ,j c − µ c C Lm L j d µ d ,j b − µ b , (39)where C jµj µ ,j µ denotes the Clebsch-Gordan coefficientand µ n is the angular momentum projection of the elec-tron state n .Substituting Eq. (38) into Eq. (31) and performing thesum of two Clebsch-Gordan coefficients, we immediately obtain the result for the matrix element of the self-energyoperator, h a | Σ( ε a ) | a i = iα π Z ∞−∞ dω × X n,L ( − j a − j n + L j a + 1 R L ( ω, anna ) ε a − ω − uε n . (40)In order to perform the angular integrations in themagnetic vertex operator, we first apply the Wigner-Eckart theorem to the matrix element of the magneticinteraction V g (which is the rank-1 spherical tensor), h n | V g | n i = ( − j − µ √ C j − µ ,j µ ( n || V g || n ) , (41)where ( . || . || . ) denotes the reduced matrix element. Nowwe can perform the angular integration in the magneticvertex matrix element as X µ µ h an | I ( ω ) | n a i h n | V g | n i = X L X L × R L ( ω, an n a ) ( n || V g || n ) , (42)where µ and µ are the angular momentum projectionsof the electron states n and n , respectively, and X L arethe angular coefficients defined by X L = X µ µ ( − j − µ √ C j − µ ,j µ J L ( an n a ) . (43)Performing the summation of three Clebsch-Gordan co-efficients with help of formulas from Ref. [65], we obtain X L = ( − j − j µ a p j a ( j a + 1)(2 j a + 1) (cid:26) j j j a j a L (cid:27) , (44)where { . . . } denotes the 6 j -symbol.Analogously, performing summations of four Clebsch-Gordan coefficients with help of formulas from Ref. [65],we perform the angular integration for the double vertexmatrix elements, X µ µ µ h an | I ( ω ) | n a i h n | V g | n i h n | V hfs | n i = X L Z L R L ( ω, an n a ) ( n || V g || n ) ( n || V hfs || n ) , (45)where the angular coefficients Z L are Z L = X j k ( − j k + j a h C j k µ a j a µ a , i j j j a L j j k j a , (46)where { . . . } denotes the 9 j -symbol. In practical calcula-tions, summations over the angular momentum projec-tions can be just as well carried out numerically.The formulas above are written in terms of explicitsummations over the Dirac spectrum, assuming the spec-tral representation of the radial Green function. In orderto use the analytical representation of the Green functionin terms of regular and irregular solutions, we would needto rewrite these formulas, identifying the components ofthe radial Green function, as X n g κ,n ( x ) g κ,n ( x ′ ) E − ε n → G κ ( E, x, x ′ ) , etc. (47)This is possible but often leads to long and unnecessarycumbersome expressions, especially for complicated dia-grams with multiple radial integrations. One can avoidthis tedious work by introducing [66] the following formalrepresentation for the radial Green function G κ ( E, x, x ′ ) = ψ κ ( E, x ) τ Tκ ( E, x ′ ) , (48)where the two-component functions ψ κ and τ κ depend on one radial argument only. The price to pay is that ψ κ and τ κ have different forms depending on the ordering ofthe radial arguments x and x ′ , ψ κ ( E, x ) = 1∆ / κ × (cid:26) φ κ ( E, x ) , when x < x ′ ,φ ∞ κ ( E, x ) , when x > x ′ , (49)and τ κ ( E, x ′ ) = 1∆ / κ × (cid:26) φ ∞ κ ( E, x ′ ) , when x < x ′ ,φ κ ( E, x ′ ) , when x > x ′ . (50)Here, φ κ and φ ∞ κ are the two-component solutions of theradial Dirac equation bounded at the origin and infinity,respectively, and ∆ κ is their Wronskian, see Eqs. (16)and (17). Employing the representation (48), we can im-mediately use formulas written via summations over theDirac spectrum for calculations with the analytical rep-resentation of the Green function. The only complicationis that the Green function is discontinuous when the tworadial arguments are equal, x = x ′ . This implies that ra-dial integrations in different matrix elements cannot beperformed independently. Their computation requires aspecial procedure, described in Sec. VII. V. CHOICE OF THE INTEGRATIONCONTOUR
The formulas presented so far contained the integrationover the virtual-photon energy ω performed along the realaxis. This choice of the integration contour, however, isnot favorable for numerical calculations, since the DiracGreen function is a highly oscillating function for largeand real energy arguments and x, x ′ → ∞ . It is advanta-geous to deform the integration contour to the region oflarge imaginary ω since the Dirac Green function acquiresan exponentially damping factor in this case. Deformingthe contour of integration, one should take care aboutpoles and branch cuts of the integrand, however. The analytical structure of the Dirac Green function isoutlined in Sec. II. The branch cuts of the photon propa-gator (15) are defined by the square root function, whichshould be understood as a limit of the regularized expres-sion with a photon mass µ > p ω + i → lim µ → p ω − µ + i
0= lim µ → p ω − µ + i p ω + µ − i . (51)The photon propagator thus has two branch cuts startingfrom ω = µ − i ω = − µ + i
0. The analytical struc-ture of the integrand of the self-energy matrix element isshown in Fig. 4.Fig. 4 also presents the deformed contour of the ω inte-gration, which we found to be optimal for most practicalcalculations. Specifically, the contour C LH consists ofthe low-energy part C L and the high-energy part C H .The low-energy part of the integration contour C L con-sists of two parts, C L + and C L − , the first of which runson the upper bank of the cut of the photon propagatorand the second, on the lower bank and in the oppositedirection. On the upper bank √ ω = ω , whereas onthe lower bank √ ω = − ω . The integrands for C L + and C L − differ only by the sign of ω in the photon propagator(and the overall sign due to the opposite directions of theintegration), thus allowing the following simplification, e iω x → e iω x − e − iω x = 2 i sin( ω x ) . The high-energy part C H is parallel to the imaginary axisand consists of two parts C H − = (∆ − i ∞ , ∆ − iǫ ) andfrom C H + = (∆ + iǫ, ∆ + i ∞ ). The integrands for C H + and C H − are typically complex conjugated, so that onecan perform the integration over C H + only, take the realpart of the result and multiply by two.In the general case of an excited reference state, thelow-energy part C L is bent in the complex plane, in orderto avoid singularities coming from virtual bound stateswith energies ε n < ε a in the electron propagator. Specif-ically, the contours C L + and C L − consist of 3 sections:(0 , δ x, − iδ y ), ( δ x, − iδ y , δ x, ), and ( δ x, , ∆), as shownon Fig. 4. The parameters of the contour δ x, , δ x, , δ y ,and ∆ may be chosen differently. In our calculations, weused the following choice (assuming the reference state a to be an excited state): ∆ = Zα ε a ; δ x, = ε a − ε s ; δ x, = 2 δ x, ; δ y = δ x, /
2. If the reference state a is theground state, there is no need to bend the low-energypart of the contour in the complex plane (as there are nointermediate states with energy 0 < ε n < ε a ); so we justintegrate along the real axis (setting δ y = 0).We note that the described contour C LH resembles thecontour used by P. Mohr in his calculations [8]. Thedifference is that he did not bend the low-energy partin the complex plane and used a different choice of theparameter ∆, ∆ = ε a .Another choice of the ω integration contour frequentlyencountered in the literature (e.g., in Refs. [50, 67, 68]) is C H+ C L- C L+ C H- C L- FIG. 4: The poles and the branch cuts of the integrandof the matrix element of the self-energy operator andthe integration contour C LH in the complex ω plane.The dashed lines (green) show the branch cuts of thephoton propagator. The poles and the branch cuts ofthe electron propagator are shown by dots and thedashed-dot line (blue). The solid line (red) shows theintegration contour C LH .the standard Wick rotation from the real into the imag-inary axis, ω → iω . In this case the intermediate stateswith energy 0 < ε n ≤ ε a lead to appearance of the poleterms, which need a special treatment. Apart for the polecontributions, small energy differences ε a − ε n appear inthe denominators of the electron propagators, leading to a rapidly varying structure of the integrand for small ω inthis choice of the contour, which may lead to numericaldifficulties. VI. INFRARED DIVERGENCIES
In this section we address the infrared (IR) reference-state divergencies which appear in the QED correctionsinvolving bound-electron propagators.We start with pointing out that the bound-state QEDcorrections do not possess the standard free-QED IR di-vergences which arise when the four-momentum p of theintermediate electron states approaches the mass shell, ρ = ( m − p ) /m = ( m − p + p ) /m →
0. For thebound-state QED corrections, the intermediate electronstates are always off mass shell ( p = ε a < m ⇒ ρ > ω →
0. Specifically,IR divergences arise in the magnetic vertex operator (34)when n = n = a and in the double vertex operator (35)when n = n = a .The general approach to the treatment of the IR di-vergencies is to separate out the divergent contributions,regularize them by introducing a finite photon mass µ in the photon propagator, evaluate the integral over ω analytically, and separate out the µ -dependent divergentterms. The divergent terms should of course cancel outwhen all relevant contributions are summed together.The evaluation of the IR-divergent integrals with the fi-nite 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 isexplicitly free from any IR divergences, h a | Λ vr ( ε a ) | a i = i π Z ∞−∞ dω " X n n h an | I ( ω ) | n a i (cid:2) h n | V g | n i − δ n n h a | V g | a i (cid:3) ( ε a − ω − u ε n )( ε a − ω − u ε n ) − X µ a ′ µ a ′′ h aa ′′ | I ( ω ) | a ′ a i (cid:2) h a ′ | V g | a ′′ i − δ a ′ a ′′ h a | V g | a i (cid:3) ( − ω + i + απ X µ a ′ µ a ′′ h aa ′′ | α µ α µ ln x | a ′ a i (cid:2) h a ′ | V g | a ′′ i − δ a ′ a ′′ h a | V g | a i (cid:3) , (52)where a ′ and a ′′ denote the reference state a with a different momentum angular projection ( µ a ′ and µ a ′′ , respectively).For the double-vertex matrix element (35) the situa-tion is somewhat more complicated because there are twotypes of divergences, the one ∝ /µ coming from threevanishing denominators ( n = n = n = a ) and theother ∝ ln µ coming from two vanishing denominators( n = n = a = n ). Still, the divergences in the double- vertex matrix can be handled with help of formulas fromAppendix B analogously to that for the magnetic vertexcase.There exists also a more economic method of handlingIR divergences in actual calculations. It relies on thefact that the matrix elements (34) and (35) are defined0so that they are overall IR finite, i.e., have a well-definedlimit at µ →
0. This means that they can be numericallyevaluated with the zero photon mass. As long as the ω integration is performed after all parts of the integrandare combined together, the integrand will have a smoothsmall- ω behaviour when integrated along the low-energypart of the integration contour C LH . The would-be IRdivergences will be cancelled numerically at a given ω be-tween different parts of the integrand. It is easy to checkthat for the magnetic vertex matrix element, both meth-ods give the identical numerical results. For the doublevertex matrix element, the numerical treatment of IR di-vergences was used in the calculation of the diamagneticshielding in Refs. [45, 63].It should be mentioned that vanishing denominatorsin the electron propagators could arise not only from theintermediate states n = a , but also from the intermedi-ate states having the same energy but the opposite par-ity as the reference state (e.g., n = 2 p / for a = 2 s for the point nuclear model). Such intermediate statesdo not cause IR divergences, since the radial matrix ele-ment in the numerator vanishes due to the orthogonalityof the wave functions, as can be seen from formulas inAppendix B.A different approach for handling the IR divergencieswas used by the Notre-Dame group [67, 69, 70]. In theirworks, numerical calculations were performed with an ex-plicit regularization parameter δ shifting the position ofthe reference-state energy, ε a → ε a (1 − δ ); the numeri-cal limit of δ → VII. RADIAL INTEGRATION
In actual calculations it is very important to find an ef-ficient way to perform multiple radial integrations. Thenumber of radial integrations is two for the self-energymatrix element, three for the magnetic vertex, and fourfor the double-vertex matrix element. In what followswe will assume that the analytical representation of theDirac Green function is used, since in the basis-set repre-sentation, the radial integrations do not cause particulardifficulties.We now formulate a general numerical approach suit-able for carrying out multiple radial integrations, firstintroduced in the context of the two-loop self-energy inRef. [66]. We start with the simplest case of the self-energy matrix element, in which the radial integral istwo-dimensional. The two-dimensional radial integralscan be schematically represented to be a linear combina-tion of terms of the following structure Z ∞ dx Z ∞ dx H ( x ) I ( x < ) L ( x > ) M ( x ) , (53)where x > = max( x , x ) and x < = min( x , x ) 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 argu-ment only, some of those being x < and x > . In particular, I ( x < ) involves φ ( x < ) from the Dirac Green function and j l ( ωx < ) from the photon propagator, and L ( x > ) involves φ ∞ ( x > ) and h (1) l ( ωx > ). It is clear that if we store allfunctions on a suitably chosen radial grid, it should bepossible to compute the integral (53) just by summingup the pre-stored data.In order to do this, we introduce a 3-dimensional radialgrid { r i,j,k } on the interval (0 , r max ) as follows. First, wefill the elements of the first layer, r i, , with i = 0 , . . . , N i ,which coarsely span the whole interval, e.g. , as r i, , = r − t t , (54)where t is uniformly distributed on the interval ( t min , t min ≈ r max . Next,we introduce a finer grid of the second layer. Specifi-cally, on each interval ( r i, , , r i +1 , , ) we introduce theset of the Gauss-Legendre abscissae { r i,j, } N j j =1 . We seethat in order to perform the outer radial integration, itis sufficient to know the integrand on the grid { r i,j, } .In order to perform the inner radial integral, we have tosplit the integration interval at the point x = x , sinceit is the discontinuity point of the integrand. We achievethis by introducing a yet finer grid of the third layer, { r i,j,k } . Specifically, for fixed values of i and j , the set { r i,j,k } N k k =1 represents the Gauss-Legendre abscissae onthe interval ( r i,j, , r i,j +1 , ) if j < N j and on the interval( r i,j, , r i +1 , , ) if j = N j . Now, for each point r i,j, ofthe outer radial integration, we can perform the innerintegral splitted into parts, (0 , r i,j, ) and ( r i,j, , r max ).We conclude that when all functions in the integrandof Eq. (53) are stored on the radial grid { r i,j,k } , the two-dimensional numerical integration can be carried out byjust summing up the pre-stored numerical values. Thedescribed procedure can be easily generalized for inte-grals of higher dimensions. So, for a computation of afour-dimensional radial integral, we need a 5-dimensionalgrid { r i,j,k,l,m } introduced in the same way as { r i,j,k } .An additional complication arises from the fact thatthe regular Dirac solution φ κ ( ε, r ) in the Dirac Greenfunction has typically the exponentially-growing be-haviour for large values of the radial argument and com-plex values of ε , whereas the irregular solutions φ ∞ κ ( ε, r )is exponentially decreasing in this region. In order toavoid numerical overflow and underflow, we store the“normalized” solutions e φ and e φ ∞ , with the approximatelarge- r and small- r behaviour pulled out, φ κ ( ε, r ) = r | κ | e cr e φ κ ( ε, r ) , (55) φ ∞ κ ( ε, r ) = r −| κ | e − cr e φ ∞ κ ( ε, r ) , (56)where c = p − ( ε/m ) . When φ κ ( ε, r ) and φ ∞ κ ( ε, r )multiply together in the Dirac Green function, the resultis usually in the range accessible in the standard double-precision (8-byte) arithmetics. A similar normalization is1required also for the regular j l and irregular h (1) l Besselsolutions, originating from the photon propagator. Withthese precautions, we are able to perform calculationscompletely within the standard double-precision arith-metics typically for κ ≤
50. For κ ≤ φ κ , φ ∞ κ ) and Bessel ( j l , h (1) l ) so-lutions but the double-precision arithmetics for the ra-dial integrations. For even higher values of κ , use of theextended-precision arithmetics becomes unavoidable [71]. VIII. MAGNETICALLY-PERTURBED GREENFUNCTION
The computation of radial integrations in diagramswith various kind of potentials can be significantly ac-celerated by introducing the first-order perturbations ofthe Green function by this potentials. Such an approachwas used long ago by Gyulassy in his evaluation of thevacuum polarization [72]. More recently, similar algo-rithms were used in calculations of various self-energycorrections (in particular, in Refs. [23, 73, 74]).In this section we describe the computation of theDirac Green function perturbed by a magnetic poten-tial V g , which will be referred to as the magnetically-perturbed Green function. Specifically, we are interestedin the radial part of the magnetically-perturbed Greenfunction, defined as G κ κ ( x , x ) = Z ∞ dx x G κ ( x , x ) V g ( x ) G κ ( x , x ) , (57)where V g ( x ) = x σ x is the radial part of the magneticpotential V g ( x ). Using the representation of the radialGreen functions in terms of the regular and irregularDirac solutions, see Eq. (16), we obtain the following ex-pressions for the magnetically-perturbed Green function.For x < x , we get G κ κ ( x , x ) = φ ∞ κ ( x ) Φ κ κ ( x ) φ ∞ T κ ( x )+ φ κ ( x ) Φ ∞ ∞ κ κ ( x ) φ T κ ( x ) , + φ κ ( x ) h Φ ∞ κ κ ( x ) − Φ ∞ κ κ ( x ) i φ ∞ T κ ( x ) , (58)whereas for x > x , G κ κ ( x , x ) = φ ∞ κ ( x ) Φ κ κ ( x ) φ ∞ T κ ( x )+ φ κ ( x ) Φ ∞ ∞ κ κ ( x ) φ T κ ( x ) , + φ ∞ κ ( x ) h Φ ∞ κ κ ( x ) − Φ , ∞ κ κ ( x ) i φ T κ ( x ) , (59)where for simplicity we assumed that φ κ and φ ∞ κ are nor-malized so that their Wronskian is unity and the func-tions Φ κ κ are defined by the integralsΦ κ κ ( x ) = Z x dx x φ T κ ( x ) V g ( x ) φ κ ( x ) , (60) Φ ∞ κ κ ( x ) = Z x dx x φ T κ ( x ) V g ( x ) φ ∞ κ ( x ) , (61)Φ ∞ κ κ ( x ) = Z x dx x φ ∞ T κ ( x ) V g ( x ) φ κ ( x ) , (62)Φ ∞ ∞ κ κ ( x ) = Z ∞ x dx x φ ∞ T κ ( x ) V g ( x ) φ ∞ κ ( x ) . (63)We observe that after storing the functions φ κ ( x ) andΦ κ ,κ ( x ) on a radial grid, we are able to construct themagnetically-perturbed Green function G κ κ ( x , x ) forany radial arguments needed in our computation. Theintegral functions Φ κ κ ( x ) are evaluated by numericalintegration with help of Gauss-Legendre quadratures. Itis important that only one integral over (0 , ∞ ) needs tobe evaluated (for a given value of the energy argument)in order to store Φ κ κ ( x ) on the whole radial grid. Anal-ogously to the case of the plain Green function, all ma-nipulations with the regular and irregular solutions needto be carried out after normalizing them according toEqs. (55) and (56) in order to prevent numerical over-flow. IX. NUMERICAL CALCULATIONS
In this section we demonstrate the technique describedin previous sections with three examples of actual calcu-lations. The first one is the calculation of the one-loopself-energy correction to the Lamb shift of a hydrogen-like ion. In Table I we present numerical results for theone-loop self-energy correction to the Lamb shift of the2 s state of hydrogen-like calcium ( Z = 20), for the pointnuclear model.The many-potential part h Σ (2+) i defined by Eq. (36)is calculated in coordinate space by the method de-scribed in the present work. Specifically, the one-potential Green function was calculated by the methoddescribed in Sec. VIII. (Alternatively, it can also be cal-culated as a derivative over the nuclear charge Z , asdescribed in Ref. [12].) For the radial integration, weused a four-dimensional grid (cid:8) r i,j,k,l (cid:9) constructed as dis-cussed in Sec. VII with the number of integration points( N i , N j , N k , N l ) = (15 , , , ω integration is car-ried out along the contour C LH introduced in Sec. V us-ing the Gauss-Legendre quadratures, after mapping ofthe integration intervals to the range (0 , | κ | = 60,with the remaining tail of the expansion estimated by aleast-square fitting in polynomials in 1 / | κ | . The remain-ing zero- and one-potential part h Σ (0+1) i is calculated inmomentum space. Their computation is relatively simpleand can be performed up to essentially arbitrary accu-racy. This part is not discussed here; we refer the readerto the original work [12].As follows from Table I, the uncertainty of the final nu-merical result for the self-energy correction comes exclu-sively from the truncation of the partial-wave expansion.It can be seen that despite the inclusion of 60 partial2waves, the resulting accuracy is significantly lower thanthat of the best literature values. There are two waysdescribed in the literature that allow to achieve a betternumerical precision. One method was developed origi-nally by Mohr [8, 9, 75] and extended by Jentschura andMohr [38, 39]. This method involves a summation ofmany thousands of partial waves and usage of extended-precision arithmetics in order to obtain very accuratenumerical results. Another method was developed inRef. [76]. It involves an additional subtraction in Σ (2+) which greatly accelerates the convergence of the partial-wave expansion and allows one to obtain accurate nu-merical results with just 20-30 partial waves. Both thesemethods are difficult to extend for computations of morecomplicated diagrams, unfortunately.Table II presents our numerical results for the self-energy correction to the g factor of the 2 s state ofhydrogen-like calcium ( Z = 20), for the point nuclearmodel. The many-potential part h Λ (2+)vr i is calculated incoordinate space by the method described in the presentwork. It is important that we calculate the magneticvertex after subtracting two first terms of its potentialexpansion, not just one as in Eq. (37). This is done inorder to accelerate the convergence of the partial-waveexpansion of the matrix element, following Refs. [14, 30].The subtracted part h Λ (0+1)vr i is calculated in momen-tum space as described in Ref. [14]. The irreducible part h Λ ir i is expressed as a non-diagonal matrix element ofthe self-energy operator; its numerical values were takenfrom Ref. [14]. The total result presented in Table IIis in good agreement with the previous value obtainedin Ref. [14]. Its numerical uncertainty comes exclusivelyfrom the truncation of the partial-wave expansion. Evenmore accurate results can be achieved if one extends thepartial-wave expansion further, as was done for the 1 s state in Ref. [71], but it requires significant efforts andintensive usage of extended-precision arithmetics.In Table III we present numerical results for the self-energy correction to the diamagnetic shielding constantof the 1 s state of hydrogen-like calcium ( Z = 20), for thepoint nuclear model. The many-potential part h Λ dvr i iscalculated in coordinate space by the method describedin the present work. We observe a slow convergence ofthe partial-wave expansion of the results presented in thetable. It can probably be accelerated by separating outthe leading term of the potential expansion (i.e., the con- tribution of the free propagators) and calculating it in themomentum space, but this has not been accomplished sofar. The other contributions to the shielding constant aredefined in Ref. [63]; the corresponding numerical resultsare taken from that work. Again, we observe that thedominant uncertainty of the final result comes from thetruncation of the partial-wave expansion. X. SUMMARY
In this paper we described the technique used in mod-ern calculations of QED corrections with the bound-electron propagators, including the notoriously problem-atic diagrams with several propagators inside the radia-tive photon loop. The bound-electron propagators aredescribed by the Green function of the Dirac equationwith the binding nuclear potential. We considered twomost widely used ways to represent the Dirac Greenfunction, the representation via the regular and irregu-lar Dirac solutions and the finite basis set representation.These representations are applicable for a wide range ofbinding potentials, including the case of the nuclear fieldmodified by a spherically-symmetric screening potentialcaused by the presence of other electrons in the atom.We demonstrated that the dominant uncertainty of theobtained results usually comes from the truncation of thepartial-wave expansion. Further extension of the partial-wave expansion is possible but often associated with largetechnical difficulties. In view of this, it is important tolook for ways to accelerate convergence of the partial-wave expansion. This was accomplished for the one-loopself-energy in Ref. [76], for the self-energy correction tothe g factor in Refs. [14, 41], and for the self-energy cor-rection to the hyperfine splitting in Refs. [15, 42]. Un-fortunately, all these methods turned out to be problem-specific, i.e. , they do not allow straightforward exten-sions to more complicated corrections. It would be thusof great importance to find a more universal approach toimprove the convergence of the partial-wave expansion insuch calculations. ACKNOWLEDGMENTS
The work presented in this paper was supported by theRussian Science Foundation (Grant No. 20-62-46006). [1] Mohr, P.J.; Plunien, G.; Soff, G. QED corrections inheavy atoms.
Phys. Rep. , , 227.[2] Beyer, H.F.; Shevelko, V.P. Introduction to the Physicsof Highly Charged Ions ; Institute of Physics Publishing:Bristol, PA, USA, 2003.[3] Indelicato, P. QED tests with highly charged ions.
J.Phys. B , , 232001.[4] Gell-Mann, M.; Low, F. Bound States in Quantum FieldTheory. Phys. Rev. , , 350. [5] Sucher, J. S -Matrix Formalism for Level-Shift Calcula-tions. Phys. Rev. , , 1448.[6] Shabaev, V.M. Two-time Green’s function method inquantum electrodynamics of high- Z few-electron atoms. Phys. Rep. , , 119.[7] Desiderio, A.M.; Johnson, W.R. Lamb Shift and BindingEnergies of K Electrons in Heavy Atoms.
Phys. Rev. A , , 1267.[8] Mohr, P.J. Self-Energy Radiative Corrections in s state of hydrogen-like calcium( Z = 20), for the point nucleus, in terms of the standardscaled function F ( Zα ) = δE/ [( α/π ) ( Zα ) /n ], where δE is a contribution to the energy in relativistic units. S l denotes the sum of partial-wave expansion P l | κ | =1 ; δS l is the increment with respect to the previous line. l S l δS l h Σ (2+) i .
268 192 85 .
541 56 3 .
273 373 86 .
515 41 0 .
973 854 86 .
967 68 0 .
452 265 87 .
223 70 0 .
256 0210 87 .
675 13 0 .
451 4415 87 .
790 89 0 .
115 7520 87 .
836 31 0 .
045 4230 87 .
869 99 0 .
033 6840 87 .
881 51 0 .
011 5250 87 .
886 57 0 .
005 0660 87 .
889 17 0 .
002 61 ∞ .
894 34 (26) 0 .
005 17 (26) h Σ (0+1) i − .
387 704Total 3 .
506 64 (26)P. J. Mohr [75] 3 .
506 648 (2)Refs. [76, 77] 3 .
506 647 (5)Hydrogen-Like Systems.
Ann. Phys. (NY) , , 26.[9] Mohr, P.J. Numerical Evaluation of the 1S / -State Ra-diative Level Shift. Ann. Phys. (NY) , , 52.[10] Soff, G.; Mohr, P. Vacuum polarization in a strong ex-ternal field. Phys. Rev. A , , 5066.[11] Manakov, N.L.; Nekipelov, A.A.; Fainshtein, A.G. Vac-uum polarization by a strong coulomb field and its contri-bution to the spectra of multiply-charged ions. Zh. Eksp.Teor. Fiz. , , 1167 [Sov. Phys. JETP , 68,673].[12] Yerokhin, V.A.; Shabaev, V.M. First order self-energycorrection in hydrogen-like systems. Phys. Rev. A , , 800.[13] Yerokhin, V.A.; Shabaev, V.M. One-loop self-energy cor-rection to the 1 s and 2 s hyperfine splitting in H-like sys-tems. Phys. Rev. A , , 012506.[14] Yerokhin, V.A.; Indelicato, P.; Shabaev, V.M. Evaluationof the self-energy correction to the g factor of S states inH-like ions. Phys. Rev. A , , 052503.[15] Yerokhin, V.A.; Jentschura, U.D. Self-energy correctionto the hyperfine splitting and the electron g factor inhydrogenlike ions. Phys. Rev. A , , 012502.[16] Yerokhin, V.A. Two-loop self-energy in the Lamb shift ofthe ground and excited states of hydrogenlike ions. Phys.Rev. A , , 052509.[17] Berestetskii, V.B.; Lifshitz, E.M.; Pitaevskii, L.P. Quan-tum Electrodynamics ; Butterworth-Heinemann: Oxford,UK, 1982.[18] Rose, M.E.
Relativistic Electron Theory ; John Wiley &Sons: New York, NY, USA, 1961.[19] Wichmann, E.H.; Kroll, N.M. Vacuum Polarization in a
TABLE II: Numerical results for the self-energycorrection to the g factor of the 2 s state ofhydrogen-like calcium ( Z = 20), for the point nucleus(in units of 10 − ). l S l δS l h Λ (2+)vr i .
130 522 17 .
563 35 − .
567 173 14 .
605 25 − .
958 104 13 .
586 86 − .
018 395 13 .
115 22 − .
471 6410 12 .
489 16 − .
626 0615 12 .
379 38 − .
109 7820 12 .
343 92 − .
035 4625 12 .
328 78 − .
015 1530 12 .
321 11 − .
007 6735 12 .
316 76 − .
004 35 ∞ .
306 43 (50) − .
010 33 (50) h Λ (0+1)vr i .
914 11 h Λ ir i .
453 02Total 2325 .
673 56 (50)Ref. [14] 2325 .
674 (5)
TABLE III: Numerical results for the self-energycorrection to the magnetic shielding constant σ of the1 s state of hydrogen-like calcium ( Z = 20), for the pointnucleus, in units of the scaled function D ( Zα ) = δσ/ [ α ( Zα ) ] where δσ is a contribution tothe shielding constant. l S l δS l h Λ dvr i − .
409 22 − .
550 9 − .
141 73 − .
559 8 − .
008 94 − .
111 6 − .
551 75 − .
438 5 − .
327 010 − .
941 8 − .
503 215 − .
986 1 − .
044 320 − .
968 3 0 .
017 725 − .
945 0 0 .
023 330 − .
925 5 0 .
019 535 − .
910 4 0 .
015 1 ∞ − .
846 7 (32) 0 .
063 7 (32) h Λ der i .
782 4 h Λ vr , Zee i .
760 7 h Λ vr , hfs i − .
404 9 h Λ po i − .
217 0Total − .
925 5 (32)Strong Coulomb Field.
Phys. Rev. A , , 843.[20] Brown, G.E.; Schaefer, G.W. Expansion in angular mo-menta in bound-state perturbation theory. Proc. R. Soc.Lond. Ser. A , , 527.[21] Gradshteyn, I.; Ryzhik, I. Table of Integrals, Series andProducts ; Academic Press: New York, NY, USA, 1994.[22] Shabaev, V.M.; Yerokhin, V.A.; Beier, T.; Eichler, J.QED corrections to the radiative recombination of anelectron with a bare nucleus.
Phys. Rev. A , , Phys. Rev. A , , 012507.[24] Shabaev, V.M.; Tupitsyn, I.I.; Yerokhin, V.A.; Plunien,G.; Soff, G. Dual Kinetic Balance Approach to Basis-Set Expansions for the Dirac Equation. Phys. Rev. Lett. , , 130405.[25] De Boor, C. A Practical Guide to Splines ; Springer: NewYork, NY, USA, 1978.[26] Johnson, W.R.; Blundell, S.A.; Sapirstein, J. Finite basissets for the Dirac equation constructed from B splines. Phys. Rev. A , , 307.[27] Salomonson, S.; ¨Oster, P. Relativistic all-order pair func-tions from a discretized single-particle Dirac Hamilto-nian. Phys. Rev. A , , 5548.[28] Salomonson, S.; ¨Oster, P. Solution of the pair equationusing a finite discrete spectrum. Phys. Rev. A , ,5559.[29] Lindgren, I.; Persson, H.; Salomonson, S.; Labzowsky,L. Full QED calculations of two-photon exchange forheliumlike-systems: Analysis in the Coulomb and Feyn-man gauges. Phys. Rev. A , , 1167.[30] Persson, H.; Salomonson, S.; Sunnergren, P.; Lindgren,I. Radiative corrections to the electron g -factor in H-likeions. Phys. Rev. A , , R2499.[31] ˚Asen, B.; Salomonson, S.; Lindgren, I. Two-photon-exchange QED effects in the 1 s s S and S states ofheliumlike ions. Phys. Rev. A , , 032516.[32] Persson, H.; Lindgren, I.; Salomonson, S.; Sunnergren, P.Accurate vacuum polarization contributions. Phys. Rev.A , , 2772.[33] Zon, B.A.; Manakov, N.L.; Rapoport, L.P. CoulombGreen’s functions in the x -representation and relativis-tic polarizability of a hydrogen-like atom. Sov. J. Nucl.Phys. , , 282.[34] Manakov, N.L.; Rapoport, L.P.; Zaprjagaev, S.A. Stur-mian expansions of the relativistic Coulomb Green func-tion. Phys. Lett. A , , 139.[35] Szymanowski, C.; V´eniard, V.; Ta¨ıeb, R.; Maquet, A.Relativistic calculation of two-photon bound-bound tran-sition amplitudes in hydrogenic atoms. Phys. Rev. A , , 700.[36] Szmytkowski, R. The Dirac - Coulomb Sturmians andthe series expansion of the Dirac - Coulomb Green func-tion: application to the relativistic polarizability of thehydrogen-like atom. J. Phys. B , , 825.[37] Indelicato, P.; Mohr, P.J. Coordinate-space approach tothe bound-electron self-energy. Phys. Rev. A , ,172.[38] Jentschura, U.D.; Mohr, P.J.; Soff, G. Calculation of theElectron Self-Energy for Low Nuclear Charge. Phys. Rev.Lett. , , 53.[39] Jentschura, U.D.; Mohr, P.J.; Soff, G. Electron self-energy for the K and L shells at low nuclear charge. Phys.Rev. A , , 042512.[40] Le Bigot, E.O.; Indelicato, P.; Mohr, P.J. QED self-energy contribution to highly excited atomic states. Phys. Rev. A , , 052508.[41] Yerokhin, V.A.; Indelicato, P.; Shabaev, V.M. Self energycorrection to the bound-electron g factor in H-like ions. Phys. Rev. Lett. , , 143001.[42] Yerokhin, V.A.; Jentschura, U.D. Electron Self-Energyin the Presence of a Magnetic Field: Hyperfine Splittingand g Factor.
Phys. Rev. Lett. , , 163001. [43] Yerokhin, V.A.; Shabaev, V.M. Accurate calculation ofself-energy screening diagrams from high- Z helium-likeatoms. Phys. Lett. A , , 274, [(E) ibid., ,210, 437].[44] Artemyev, A.N.; Shabaev, V.M.; Yerokhin, V.A. Vacuumpolarization screening corrections to the ground-state en-ergy of two-electron ions. Phys. Rev. A , , 3529.[45] Yerokhin, V.A.; Pachucki, K.; Harman, Z.; Keitel, C.H.QED Theory of the Nuclear Magnetic Shielding in Hy-drogenlike Ions. Phys. Rev. Lett. , , 043004.[46] Blundell, S.A.; Mohr, P.J.; Johnson, W.R.; Sapirstein,J. Evaluation of two-photon exchange graphs for highlycharged heliumlike ions. Phys. Rev. A , , 2615.[47] Yerokhin, V.A.; Artemyev, A.N.; Shabaev, V.M.; Sysak,M.M.; Zherebtsov, O.M.; Soff, G. Two-Photon ExchangeCorrections to the 2 p / − s Transition Energy in Li-LikeHigh- Z Ions.
Phys. Rev. Lett. , , 4699.[48] Mohr, P.J.; Sapirstein, J. Evaluation of two-photon ex-change graphs for excited states of highly charged heli-umlike ions. Phys. Rev. A , , 052501.[49] Volotka, A.V.; Glazov, D.A.; Shabaev, V.M.; Tupitsyn,I.I.; Plunien, G. Many-Electron QED Corrections to the g Factor of Lithiumlike Ions.
Phys. Rev. Lett. , ,253004.[50] Blundell, S.A.; Snyderman, N.J. Basis-set approach tocalculating the radiative self-energy in highly ionizedatoms. Phys. Rev. A , , R1427.[51] Artemyev, A.N.; Shabaev, V.M.; Yerokhin, V.A. Rela-tivistic nuclear recoil corrections to the energy levels ofhydrogenlike and high- Z lithiumlike atoms in all ordersin αZ . Phys. Rev. A , , 1884.[52] Yerokhin, V.A.; Shabaev, V.M. Nuclear Recoil Effect inthe Lamb Shift of Light Hydrogenlike Atoms. Phys. Rev.Lett. , , 233002.[53] Shabaev, V.M.; Glazov, D.A.; Malyshev, A.V.; Tupitsyn,I.I. Recoil Effect on the g Factor of Li-Like Ions.
Phys.Rev. Lett. , , 263001.[54] Volotka, A.V.; Glazov, D.A.; Shabaev, V.M.; Tupitsyn,I.I.; Plunien, G. Screened QED Corrections in Lithium-like Heavy Ions in the Presence of Magnetic Fields. Phys.Rev. Lett. , , 033005.[55] Glazov, D.A.; Volotka, A.V.; Shabaev, V.M.; Tupitsyn,I.I.; Plunien, G. Evaluation of the screened QED cor-rections to the g factor and the hyperfine splitting oflithiumlike ions. Phys. Rev. A , , 062112.[56] Sunnergren, P.; Persson, H.; Salomonson, S.; Schnei-der, S.M.; Lindgren, I.; Soff, G. Radiative corrections tothe hyperfine-structure splitting of hydrogenlike systems. Phys. Rev. A , , 1055.[57] Persson, H.; Salomonson, S.; Sunnergren, P.; Lindgren,I. Two-Electron Lamb-Shift Calculations on HeliumlikeIons. Phys. Rev. Lett. , , 204.[58] Indelicato, P.; Mohr, P.J. Quantum electrodynamic ef-fects in atomic structure. Theor. Chim. Acta , ,207.[59] Persson, H.; Schneider, S.M.; Greiner, W.; Soff, G.; Lind-gren, I. Self-Energy Correction to the Hyperfine Struc-ture Splitting of Hydrogenlike Atoms. Phys. Rev. Lett. , , 1433.[60] Blundell, S.A.; Cheng, K.T.; Sapirstein, J. All-OrderBinding Corrections to Muonium Hyperfine Splitting. Phys. Rev. Lett. , , 4914.[61] Yerokhin, V.A.; Shabaev, V.M. Self-energy correction tothe ground state hyperfine splitting of Bi . JETP Lett. , 18.[62] Shabaev, V.M.; Pachucki, K.; Tupitsyn, I.I.; Yerokhin,V.A. QED Corrections to the Parity-Nonconserving 6 s -7 s Amplitude in
Cs.
Phys. Rev. Lett. , , 213002.[63] Yerokhin, V.A.; Pachucki, K.; Harman, Z.; Keitel, C.H.QED calculation of the nuclear magnetic shielding forhydrogenlike ions. Phys. Rev. A , , 022512.[64] Snyderman, N.J. Electron Radiative Self-Energy ofHighly Stripped Heavy Atoms. Ann. Phys. (NY) , , 43.[65] Varshalovich, D.A.; Moskalev, A.N.; Khersonski˘i, V.K. Quantum Theory of Angular Momentum ; World Scien-tific: Singapure, 1988.[66] Yerokhin, V.A.; Indelicato, P.; Shabaev, V.M. Evaluationof the two-loop self-energy correction to the ground stateenergy of H-like ions to all orders in Zα . Eur. Phys. J.D , , 203.[67] Sapirstein, J.; Cheng, K.T. Determination of the two-loop Lamb shift in lithiumlike bismuth. Phys. Rev. A , , 022502.[68] Sapirstein, J.; Cheng, K.T. Calculation of radiative cor-rections to hyperfine splittings in the neutral alkali met-als. Phys. Rev. A , , 022512.[69] Blundell, S.A.; Cheng, K.T.; Sapirstein, J. Radiative cor-rections in atomic physics in the presence of perturbingpotentials. Phys. Rev. A , , 1857.[70] Sapirstein, J.; Cheng, K.T. Hyperfine splitting in lithi- umlike bismuth. Phys. Rev. A , , 032506.[71] Yerokhin, V.A.; Harman, Z. One-loop electron self-energy for the bound-electron g factor. Phys. Rev. A , , 060501.[72] Gyulassy, M. Higher Order Vacuum Polarizatin for FiniteRadius Nuclei. Nucl. Phys. A , , 497.[73] Artemyev, A.N.; Shabaev, V.M.; Tupitsyn, I.I.; Plunien,G.; Yerokhin, V.A. QED Calculation of the 2 p / -2 p / Transition Energy in Boronlike Argon.
Phys. Rev. Lett. , , 173004.[74] Artemyev, A.N.; Shabaev, V.M.; Tupitsyn, I.I.; Plunien,G.; Surzhykov, A. and Fritzsche, S. Ab initio calcula-tions of the 2 p / p / fine-structure splitting in boron-like ions. Phys. Rev. A , , 032518.[75] Mohr, P.J. Self-energy correction to one-electron energylevels in a strong Coulomb field. Phys. Rev. A , ,4421.[76] Yerokhin, V.A.; Pachucki, K.; Shabaev, V.M. One-loopself-energy correction in a strong binding field. Phys.Rev. A , , 042502.[77] Shabaev, V.M.; Tupitsyn, I.I.; Yerokhin, V.A. Model op-erator approach to the Lamb shift calculations in rela-tivistic many-electron atoms. Phys. Rev. A , ,012513.[78] Johnson, W.R.; Blundell, S.A.; Sapirstein, J. Many-bodyperturbation-theory calculations of eneryg levels alongthe lithium isoelectronic sequence. Phys. Rev. A , , 2764. Appendix A: Relativistic Slater radial integral
The matrix element of the electron-electron interaction operator (32) is represented in the form (38), where R J ( ω, abcd ) is the relativistic generalization of the Slater radial integral. The explicit expression for R J can beobtained, e.g. , by reformulating formulas presented in Appendix of Ref. [78]. The result for the radial integral R J inthe Feynman gauge is written as [12] R J ( ω, abcd ) = (2 J + 1) Z ∞ dx dx ( x x ) h ( − J C J ( κ a , κ c ) C J ( κ b , κ d ) g J ( ω, x < , x > ) W ac ( x ) W bd ( x ) − X L ( − L g L ( ω, x < , x > ) X ac,JL ( x ) X bd,JL ( x ) i , (A1)where x > = max( x , x ), x < = min( x , x ), the functions W ab and X ab,JL are defined by W ab ( x ) = g a ( x ) g b ( x ) + f a ( x ) f b ( x ) , (A2) X ab,JL ( x ) = g a ( x ) f b ( x ) S JL ( − κ b , κ a ) − f a ( x ) g b ( x ) S JL ( κ b , − κ a ) . (A3)Here, g n , f n are the upper and the lower radial components of the Dirac wave function, respectively. The function g l ( ω, x < , x > ) is the radial part of the partial-wave expansion of the photon propagator, e iωx x = X l (2 l + 1) g l ( ω, x < , x > ) P l ( ξ ) , (A4)where P l ( ξ ) is the Legendre polynomial, ξ = ˆ x · ˆ x , g l (0 , x < , x > ) = 12 l + 1 x l< x l +1 > , (A5) g l ( ω, x < , x > ) = iω j l ( ωx < ) h (1) l ( ωx > ) , (A6)6and j l ( z ), h (1) l ( z ) are the spherical Bessel functions. The angular coefficients S JL ( κ a , κ b ) differ from the zero only for L = J − J , J + 1 and can be written for J = 0 as follows: S J J +1 ( κ a , κ b ) = r J + 12 J + 1 (cid:18) κ a + κ b J + 1 (cid:19) C J ( − κ b , κ a ) , (A7) S J J ( κ a , κ b ) = κ a − κ b p J ( J + 1) C J ( κ b , κ a ) , (A8) S J J − ( κ a , κ b ) = r J J + 1 (cid:18) − κ a + κ b J (cid:19) C J ( − κ b , κ a ) . (A9)In the case J = 0 there is only one nonvanishing coefficient S ( κ a , κ b ) = C ( − κ b , κ a ) . The coefficients C J ( κ b , κ a ) aregiven by C J ( κ b , κ a ) = ( − j b +1 / p (2 j a + 1)(2 j b + 1) (cid:18) j a J j b / − / (cid:19) Π( l a , l b , J ) , (A10)where the symbol Π( l a , l b , J ) is unity if l a + l b + J is even, and zero otherwise. Appendix B: Infrared divergent integrals
In this section we evaluate the infrared divergent integrals J α with α = 2 and 3, defined as J α ( abcd ) = i π Z ∞−∞ dω h ab | I µ ( ω ) | cd i ( − ω + i α . (B1) I µ ( ω ) is the electron-electron interaction operator with a finite photon mass µ , written in the Feynman gauge as I ( ω, x ) = α (cid:0) − α · α (cid:1) e i √ ω − µ + i x x . (B2)The integral over ω with α = 2 is evaluated as i π Z ∞−∞ dω − ω + i e i √ ω − µ + i x x = − πx Z ∞ µ ω sin (cid:0)p ω − µ + i x (cid:1) = − πx Z ∞ dt t sin tx ( t + µ ) / = − π Z ∞ dt cos tx ( t + µ ) / = − π Z ∞ dt cos tx − cos tt − π Z ∞ dt cos t ( t + µ ) / = 1 π (cid:16) ln µ γ + ln x (cid:17) + O ( µ ) . (B3)Therefore, J ( abcd ) = απ (cid:16) ln µ γ (cid:17)(cid:2) h a | c i h b | d i − h a | α | c ih b | α | d i (cid:3) + απ (cid:10) ab (cid:12)(cid:12)(cid:0) − α · α (cid:1) ln x (cid:12)(cid:12) cd (cid:11) , (B4)where we dropped terms vanishing in the limit µ →
0. Analogously, we obtain J ( abcd ) = α µ (cid:2) h a | c i h b | d i − h a | α | c ih b | α | d i (cid:3) − α (cid:10) ab (cid:12)(cid:12)(cid:0) − α · α (cid:1) x (cid:12)(cid:12) cd (cid:11) . (B5)One can see that infrared divergences arise from terms of the type J α ( abab ), since in this case h a | a i h b | b i − h a | α | a ih b | α | b i = 1 . We need also consider the case of c = e a and d = e b , where the state e n has the same energy as n but the opposite parity(e.g., e n = 2 p / and n = 2 s for the point nucleus). Such states do not cause any infrared divergences since h a | e a i = 0due to orthogonality and the matrix element with α vanishes because of degeneracy in energy, h a | α | e a i = h a | i [ H D , r ] | e a i = i ( ε a − ε e a ) h a | r | e a i = 0 ..