Two-photon exchange corrections to the $g$ factor of Li-like ions

We report calculations of QED corrections to the $g$ factor of Li-like ions induced by the exchange of two virtual photons between the electrons. The calculations are performed within QED theory to all orders in the nuclear binding strength parameter $Z\alpha$, where $Z$ is the nuclear charge number and $\alpha$ is the fine-structure constant. In the region of low nuclear charges we compare results from three different methods: QED, relativistic many-body perturbation theory, and nonrelativistic QED. All three methods are shown to yield consistent results. With our calculations we improve the accuracy of the theoretical predictions of the $g$ factor of the ground state of Li-like carbon and oxygen by about an order of magnitude. Our theoretical results agree with those from previous calculations but differ by 3-4 standard deviations from the experimental results available for silicon and calcium.


I. INTRODUCTION
Modern Penning-trap experiments based on the continuous Stern-Gerlach effect provide very precise measurements of the Zeeman splitting of energy levels in one-and few-electron ions [1][2][3]. The linear Zeeman splitting is usually parameterized in terms of the g factor of the atomic system. The fractional accuracy of the recent measurements of the g factors of H-like and Li-like ions has reached few parts in 10 −11 [4,5]. Combined with dedicated theoretical calculations, these measurements provided the determination of the electron mass [6] and one of the best tests of the bound-state quantum electrodynamics (QED) [7]. Extension of these tests towards heavier ions are anticipated in the future [8], which might open new ways for determination of the fine-structure constant α [9,10] and searches for physics beyond the Standard Model [11].
In view of the very high accuracy of the measurements, theoretical investigations of atomic g factors often need to be carried out without any expansion in the nuclear binding strength parameter Zα (where Z is the nuclear charge number). In such calculations, the electron-electron interaction has to be treated by perturbation theory. The starting point of the perturbation expansion is the hydrogenic approximation, i.e., the approximation of non-interacting electrons. The electron-correlation corrections come from the exchange of virtual photons between the electrons. An exchange by each photon leads to the suppression of the corresponding correction by a parameter of 1/Z. The first-order perturbation correction ∼ 1/Z 1 is due to the one-photon exchange. This correction is relatively simple and was calculated for Li-like ions first in Ref. [12] and later reproduced in Refs. [13,14].
The QED calculation of the two-photon exchange correction ∼ 1/Z 2 is a difficult task. First calculations of this correction were accomplished in Refs. [15,16]. In these studies, results were reported for just four ions and their numerical uncertainty was significant on the level of the current experimental precision. In the present work we will perform an independent calculation of the two-photon exchange correction for the ground state of Li-like ions. Our goals will be to cross-check the previous calculations, to improve the numerical accuracy, and to study the Z-dependence of the twophoton correction in the low-Z region, checking the consistency of the applied method with the Zα-expansion calculations performed recently in Ref. [17].

II. ELECTRONIC STRUCTURE CORRECTIONS TO THE g FACTOR
In the present work we assume the nucleus to be spinless and the electron configuration to be a valence electron v beyond a closed shell of core electrons denoted by c. Contributions to the g factor can be formally obtained as corrections induced by the effective magnetic interaction [18] V g (r) = 1 where α is the vector of Dirac matrices and µ v is the angular momentum projection of the valence electron.
To the zeroth order in the electron-electron interaction, the g factor of the ground state of a Li-like ion is given by the expectation value of the magnetic potential V g on the hydrogenic Dirac wave function of the valence 2s state. For the point nucleus, the result is known in the closed form g Dirac = 2s|V g |2s = 2 3 2 1 − (Zα) 2 + 2 Corrections to the g factor of a Li-like ion due to the presence of core electrons are evaluated by perturbation theory in the electron-electron interaction, with the expansion parameter 1/Z. The leading correction of order 1/Z 1 is induced by the one-photon exchange between the valence and core electrons. The corresponding correction was calculated in Ref. [12] (see also Ref. [13]) and can be written as where the summation runs over the angular-momentum projection of the core electron µ c and Λ 1ph (abcd) = n =a a|V g |n nb|I(∆ db )|cd ε a − ε n + 1 4 ab|I ′ (∆ db )|cd d|V g |d − b|V g |b .
The electron-electron interaction operator I(ω) is defined as where α µ a = (1, α a ) is the four-vector of Dirac matrices acting on r a , D µν is the photon propagator, and ω is the photon energy. In the present work we use the photon propagator in the Feynman and Coulomb gauges. In the Feynman gauge, the electron-electron interaction takes the simplest form, where r 12 = |r 12 | = |r 1 − r 2 |, and |ω| should be understood as |ω| = √ ω 2 + iǫ, where ǫ is a positive infinitesimal addition. The electron-electron interaction operator in the Coulomb gauge reads It can be easily seen that the one-photon exchange correction ∆g 1ph given by Eq. (3) can be obtained from the corresponding correction to the Lamb shift, (8) by perturbing this expression with the magnetic potential V g . Specifically, one perturbs the one-electron wave functions, and energies, ε a → ε a + δε a , δε a = a|V g |a . In this work, we use this approach in order to obtain formulas for the two-photon exchange corrections to the g factor. We start with the two-photon exchange correction for the Lamb shift, graphically represented in Fig. 1. The Feynman diagrams for the g factor in Fig. 2 are obtained from the Lambshift diagrams by inserting the magnetic interaction V g in all possible ways. The corresponding formulas for the g factor are obtained by using formulas for the two-photon exchange correction for the Lamb shift derived in Refs. [19,20] and perturbing them with the magnetic potential V g . The two-photon exchange correction to the g factor is conveniently represented as a sum of the direct ("dir"), exchange ("ex"), and the three-electron ("3el") contributions, obtained as perturbations of the corresponding Lamb-shift corrections. Furthermore, each of the three contributions is sub-divided into the irreducible ("ir") and reducible ("red") parts. The reducible parts are induced by the intermediate states degenerate in energy with the energy of the reference state of the ion. We thus represent the total two-photon exchange correction to the g factor as the sum of three irreducible and three reducible contributions, ∆g 2ph = ∆g ir,dir + ∆g ir,ex + ∆g 3el,ir + ∆g red,dir + ∆g red,ex + ∆g 3el,red .
We now examine each of these terms one by one.

A. Direct irreducible part
The direct irreducible contribution comes from the ladder ("lad") and crossed ("cr") diagrams. For the Lamb shift, this contribution is given by Eq. (32) of Ref. [20]. Changing the variable ω → −ω in the ladder part and using the property I(ω) = I(−ω), we write the expression as where∆ an ≡ ε a − ε n (1 − i0). Furthermore, where µ 1 and µ 2 denote the angular-momentum projections of the states n 1 and n 2 , respectively. The summations over n's run over the complete spectrum of the Dirac equation, implying the sum over the corresponding relativistic angular quantum numbers κ n and the principal quantum numbers of the discrete spectrum and the integration over the continuum part of the spectrum. The terms excluded from the summation over n 1 and n 2 in Eq. (12) will be accounted for by the reducible part.

C. Direct reducible part
The reducible part of the two-electron diagrams for the Lamb shift is given by Eq. (41) of Ref. [20]. Separating the direct contribution, we write it as We note that the crossed (cv) term in the above expression ex-actly coincides with the one excluded from the summation in Eq. (12). The ladder (cv) and (vc) terms in the above expression are very similar to those excluded from the summation in Eq. (12) but differ by signs of i0. Specifically, the terms excluded from the summation over n 1,2 in Eq. (12) contained poles both at ω = i0 and ω = −i0 (or at ω = −∆ + i0 and ω = −∆ − i0), thus "squeezing" the integration contour between the two poles, causing singularities. By contrast, the ladder terms in Eq. (32) have double poles, from one side of the integration contour. Therefore, the integration contour can be "moved away" from the pole (assuming a finite photon mass in the case of ω = 0), so there is no real singularities in Eq. (32). Taking into account that the ladder (cv) term and the crossed (cv) term cancel each other (as proven in Ref. [19]), the expression in simplified further to yield where P denotes the principal value of the integral and F ′ (ω) = ∂F (ω)/(∂ω).
Corrections to the g factor arise through perturbations of the Lamb-shift formulas with the magnetic potential V g . We divide them into two parts, where the first term is induced by perturbations of the wave functions and the second, by perturbations of the energies. The perturbed-wave-function part is immediately obtained from Eq. (33) as where δF lad,dir is defined by Eq. (15). The derivation of the energy-perturbed reducible part is more difficult and is carried out by perturbing the formulas given by Eq. (47) of Ref. [19]. This derivation is potentially problematic, because vanishing contributions to the Lamb shift may induce nonzero magnetic perturbations. For example, the energy difference ∆ an = ε a − ε n induces a perturbation a|V g |a − n|V g |n . If ε n = ε a and µ n = µ a , the energy difference vanishes but the magnetic perturbation survives. In order to avoid potential ambiguities, we fix the reducible part by requirement of the gauge invariance of the total correction to the g factor. The result for the direct energy-perturbation reducible part is where It should be pointed out that the second integration by parts, leading to the second derivative of the photon exchange operator I ′′ (ω) in Eq. (36) is potentially troublesome. The reason is that the imaginary part of the first derivative I ′ (ω) is discontinuous at ω = 0. Specifically, Im I ′ (0 + ) = −Im I ′ (0 − ) = 0. This discontinuity leads, in principle, to appearance of additional off-integral terms in Eq. (36). We found, however, that their numerical contributions are completely negligible for the case under consideration in the present paper. The same holds for the exchange reducible part.

D. Exchange reducible part
The reducible exchange correction for the Lamb shift is given by Eq. (41) of Ref. [20]. We write this correction as We observe that the crossed (cc) and (vv) terms in the above expression exactly coincide with the two terms excluded from the summation in Eq. (22). The ladder (cv) and (vc) terms in the above expression are similar to those excluded from the summation in Eq. (22) but differ from them by the signs of i0. We evaluate the above expression by integrating by parts and taking the principal value of the integral, separating the pole contribution. The result is It should be mentioned that the individual terms in the brackets under the integral in the above formula contain singularities at ω = 0 and ω = ∆. When the ladder and exchange terms are combined together, however, the singularities disappear and the principal value of the resulting integral becomes well defined and can be calculated numerically. The reducible exchange contribution for the g factor is the sum of perturbations of the wave functions and perturbations of the energies, ∆g red,ex = ∆g red,ex,wf + ∆g red,ex,en , where The three-electron irreducible contribution to the Lamb shift is given by Eq. (14) of Ref. [20], Here, "1", "2", and "3" label the three electrons of the ions (in arbitrary order), the operators P and Q permute the initial-state and the final-state electrons, (−1) P and (−1) Q are the sign of permutations, and the prime on the sum means that the terms with the vanishing denominator are excluded from the summation. Furthermore, ∆ ab ≡ ε a − ε b , and I abcd (∆) ≡ ab|I(∆)|cd . The corrections to the g factor are obtained as first-order perturbations of Eq. (41) by V g . It is convenient to split the whole contribution into the perturbations of the external wave functions ("pwf"), external energies ("en"), and the propagator ("ver"), ∆g 3el,ir = ∆g 3el ir,pwf + ∆g 3el ir,en + ∆g 3el ir,ver , where the operator Ξ acts on energy denominators ∆ 1 , ∆ 2 as follows: The three-electron reducible correction to the Lamb shift is given by Eq. (19) of Ref. [20], The corresponding corrections to the g factor arise as perturbations of the wave functions and energies,

III. MBPT APPROXIMATION
In this section we obtain formulas for the two-photon exchange correction to the g factor in the approximation of the relativistic many-body perturbation theory (MBPT). The corresponding formulas can be obtained from the QED expressions by (i) using the Coulomb gauge in the photon propagators and neglecting the energy dependence, I(ω) → I Coul (0), and (ii) restricting the summations over the Dirac spectrum to the positive-energy part. Under these assumptions, all reducible contributions vanish and the integrals over ω can be performed by the Cauchy theorem. The ω integral for the crossed contribution vanishes, so the total two-electron correction comes only from the ladder irreducible part. Performing the ω integrations in Eqs. (21) and (25), we obtain the two-electron contribution in the MBPT approximation as Here, the prime on the sum means that the terms with the vanishing denominator should be omitted, and the summation over n 1 and n 2 is performed over the positive-energy part of the Dirac spectrum. We note that Eq. (49) can be also obtained directly by perturbing the two-photon MBPT correc-tion for the Lamb shift, given by Eq. (43) of Ref. [20]. The three-electron MBPT correction is immediately obtained from Eqs. (42)-(45), after the substitution I(ω) → I Coul (0) and the restriction of the summations to the positive-energy part of the spectrum.
We note that the standard formulation of MBPT assumes the restriction of all summations over the Dirac spectrum to the positive-energy part. The consistent treatment of the negative-energy spectrum is possible only within the QED theory. However, it can be easily observed that one can include some negative-energy contributions already in the MBPT formulas, namely, in those cases when it does not lead to the so-called continuum dissolution, i.e., vanishing energy denominators. Specifically, one can include the negativeenergy spectrum in the three-electron contributions, Eqs. (42)-(45), and in the magnetic perturbations of the wave functions, Eq. (9). We will refer to this variant of the MBPT as "MBPTneg". We will demonstrate that such partial inclusion of the negative-energy spectrum within MBPT is crucially important to approximately reproduce the QED results in the region of small nuclear charges, whereas the standard MBPT yields a very much different result, even in the limit of Z → 0. Previously the same conclusion was drawn by the St. Petersburg group [21,22].
The connection between the QED and MBPT formulas was extensively used in this work for checking the numerical procedure for the ω integrations. Specifically, after neglecting the energy dependence of the photon propagators in the Coulomb gauge, we checked that the numerical ω integration yields the same result as the analytical integration by the Cauchy theorem.

IV. NUMERICAL EVALUATION
We now turn to the numerical evaluation of the two-photon exchange corrections. Since the calculation of the threeelectron contributions is relatively straightforward, we concentrate mainly on the two-electron terms. The direct and exchange irreducible contributions given by Eqs. (21) and (25) represent the main computational difficulty. It is advantageous to deform the contour of the ω integration in them, in order to escape strong oscillations of the photon propagators ∝ e i|ω|r12 for large real values of ω. Deforming the contour, one needs to take into account the branch cuts of the photon propagators and the pole structure of the Dirac propagators. The analytical structure of the integrand as a function of complex ω is shown in Fig. 3 for the direct part and in Fig. 4 for the exchange part, respectively.
For the evaluation of the direct irreducible contribution, we use two different choices of the ω-integration contour. The first choice is the standard Wick rotation, ω → iω, which splits the correction into the pole contribution and the integral along the imaginary ω axis. This contour was used in the previous Lamb-shift calculations [20,23]. The advantage of this choice is that the analysis of the pole terms is the simplest. There are, however, also some difficulties. The first problem is the rapidly-varying structure of the integrand in the vicinity of ω = 0, due to poles of the electron propagators lying near the imaginary axis. The second difficulty is that the contributions with n 2 = v in Eq. (21) contain singular terms ∼ 1/ω 2 , which need to be integrated by parts before the numerical evaluation.
In order to achieve a more regular behaviour of the inte-grand for small ω, we adopted the contour C D , shown in Fig. 3. This contour is convenient for the numerical evaluation, especially for low values of Z. Its disadvantage is the presence of a pole on the low-energy part of the integration contour and thus the need to evaluate the principal value of the integral. We find, however, that the contour C D is very similar to the contour C X used in the evaluation of the exchange contribution and discussed in detail below. Because of this similarity, we were able to use essentially the same numerical procedure both for the direct and the exchange part. We checked that our numerical evaluation of the integral along the contour C D leads to the same results as the integration along the Wick-rotated contour.
For the numerical evaluation of the exchange irreducible part, we use the contour C X depicted in Fig. 4. This contour was suggested in Ref. [23] for the Lamb shift and later used for the g factor and hyperfine structure in Refs. [15,24]. As can be seen from Fig. 4, the deformation of the contour from (−∞, ∞) to C X leads to appearance of pole terms at ω = 0 and ω = ∆. In the case of the Lamb shift, the pole terms are identified as follows, for the ladder contribution, where∆ i denotes ∆ i with the infinitesimal imaginary addition according to Eqs. (21) and (25). For the crossed contribution, the corresponding equation reads For the g factor, formulas with squared energy denominators are required, which can be obtained by a formal differentiation of the above formulas over ∆ 1 and ∆ 2 . For a numerical evaluation, the integral over C X is represented as a sum of three pieces, where δ is a free parameter 0 < δ < ∆. A typical value of δ = ∆/2 was used. An advantage of the contour C X is that the integrand has a more regular behaviour at the end points of the intervals, ω = 0 and ω = ∆, because Im I(ω) ∝ ω as ω → 0. There are, however, singularities inside the intervals along the real axis, and thus the infinitesimal imaginary terms i0 should be retained for them. For the g factor, we encounter single, double, and even triple poles on the integration contour. Specifically, for the v = 2s reference state, the singularities arise from the intermediate states n 1 and/or n 2 = 2p 1/2 , whose energy ε 2p 1/2 < ε 2s is separated from the reference-state energy by the finite nuclear size effect. To deal with these singularities, we introduce subtractions obtained by expanding the integrand in the Taylor series in the vicinity of the poles. The subtractions remove singularities and make the integrand a regular and smooth function suitable for the numerical integration. The subtracted terms are then re-added, with the principal-value integrals calculated analytically. The corresponding formulas are summarized in Appendix A.
In order to check our numerical procedure, we performed calculations also by using a different integration contour, namely, the contour C irr suggested in Ref. [20] (shown in Fig. 5 of that work). A very good agreement of numerical results obtained with two different contours was used as a confirmation of the internal consistency of the numerical procedure.
For the numerical evaluation of the reducible direct and exchange contributions, we used the ω-integration contour con- sisting of three sections: The parameter δ of the contour was taken to be δ > ∆, which allowed us to evaluate the principal value of the integrals at points ω = 0 and ω = ±∆.
The summations over the Dirac intermediate states were performed by using the dual kinetic balance basis-set method [25], with the basis set constructed with the B-splines. The standard two-parameter Fermi model was used to represent the nuclear charge distribution, with the nuclear radii taken from Ref. [26]. The infinite partial-wave summation over the relativistic angular momentum quantum number κ was performed up to |κ max | = 25, with the remaining tail estimated by the polynomial fitting of the expansion terms in 1/|κ|. The largest numerical uncertainty was typically induced by convergence in the number of basis functions. Our final values were typically obtained by performing calculations with N = 85 and N = 105 B-splines and extrapolating the results to N → ∞ as δg = δg N =105 + 0.93 (δg N =105 − δg N =85 ), where the numerical coefficient was obtained by analysing the convergence pattern of our numerical results. An example of the convergence study with respect to N is presented in Table I.
In the present work we perform calculations in the Feynman and the Coulomb gauge. The expressions for the matrix elements of the electron-electron interaction in the Coulomb gauge are summarized in Appendix B. We note that in the present work (unlike, e.g., in Ref. [20]) we use the expression for the Coulomb-gauge matrix element [Eq. (B3)] that does not rely on commutator relations for the wave functions. This expression is valid for the general case when the wave functions in the matrix element are not eigenfunctions of the Dirac Hamiltonian, in particular, when they are the magnetic perturbations of the Dirac wave functions. Another advantage of this expression is that it allows a numerical evaluation of the Coulomb-gauge radial integrals and their derivatives for very small but nonvanishing photon energies ω. The region of small but nonzero ω is usually numerically unstable for the expressions based on the commutator relations, especially for the second derivative of the photon propagator I ′′ (ω).

V. RESULTS AND DISCUSSION
Numerical results of our calculations of the two-photon exchange corrections for the ground state of Li-like ions are pre- sented in Tables II and III. Table II contains a breakdown of our calculations in two gauges for Z = 14 and Z = 83 and demonstrates the gauge invariance of our numerical results. Table III presents the final results of our calculations for Z = 6 -92. It also compares results of the QED calculation with those obtained within the standard MBPT and the MBPT with the partial inclusion of the negative-energy spectrum ("MBPT-neg"). We observe that the standard MBPT yields the two-photon exchange correction by about three times larger than the complete QED results. The disagreement is evidently present even in the limit of Z → 0. On the contrary, the MBPT-neg approach closely reproduces the QED treatment in the region of low values of Z. Previously the same conclusion was reached by the St. Petersburg group [21,22].
(53) Fig. 5 shows the comparison of the QED, MBPT-neg, and NRQED results. We observe that all three methods yield results converging to each other in the limit Z → 0. The difference between the results by different methods scales as ∝ (Zα) 2 as expected. It is interesting that the MBPT-neg approach does not yield any significant improvement over the NRQED treatment for low-and medium-Z ions. The QED calculation of the two-photon exchange correction for the g factor of Li-like ions was previously carried  Numerical results for the two-photon exchange correction to the g factor of the ground state of Li-like ions, in units of 10 −6 . QED results are obtained in the Feynman gauge. "MBPT" labels results obtained within the standard relativistic many-body perturbation theory. "MBPT-neg" labels results obtained within MBPT supplemented by the correction from the negative-continuum part of the Dirac spectrum. out in Refs. [15,16]. Unfortunately, the numerical results were presented only for four ions and mostly in the form of the total electron-electron interaction correction. The only ion for which the calculations are directly comparable is silicon, Z = 14, for which we find some tension. Our calculation yields −6.8787 (1) × 10 −6 , whereas Ref. [16] reported −6.876 × 10 −6 . As an additional cross-check, we performed calculations for the two-photon exchange correction to the ground-state hyperfine splitting of Li-like bismuth (which is another example of a magnetic perturbation potential) and found agreement with results listed in Table I of Ref. [15].
Having obtained results for the two-photon exchange correction, we are now in a position to update the theoretical predictions for the ground-state g factor of Li-like ions. Table IV presents a compilation of all known binding corrections to the g factor of the ground state of four Li-like ions, C 3+ , O 5+ , Si 11+ , and Ca 17+ . As compared to the analogous compilation in our previous investigation [17], we introduced several improvements. The two-photon exchange correction (i.e., the electronic-structure contribution of relative order 1/Z 2 ) is computed in the present work. Beside this, we included the one-loop 1/Z and 1/Z 2+ QED effects from our recent work [28] and the nuclear recoil corrections of rel-ative orders 1/Z 0 , 1/Z 1 , and 1/Z 2+ calculated by Shabaev et al. [29]. Furthermore, we added the two-loop (Zα) 5 effects calculated recently by Czarnecki and co-workers [30,31] and by us [32]. The two-loop results in those studies were reported for the 1s hydrogenic state. We here convert them to the 2s state by assuming the 1/n 3 scaling. Having in mind that the result for the nonlogarithmic (Zα) 5 term is not complete, we ascribe the uncertainty of 20% to the two-loop (Zα) 5 correction. The uncertainty due to higher-order two-loop effects was evaluated on the basis of available one-loop results, with the extension factor of 2.
One of the largest uncertainties of the theoretical predictions comes from the higher-order electronic-structure correction ∼ 1/Z 3+ . The values in Table IV for this correction are obtained within NRQED in Ref. [17]. Theoretical estimates for their uncertainties are obtained by taking the relative deviation of the NRQED and full QED results for the 1/Z 2 correction and multiplying it by the extension factor of 2.
The comparison presented in Table IV shows agreement of the present theoretical g-factor values with previous theoretical predictions. In particular, for C 3+ and O 5+ , our results are in excellent agreement with, but 5-10 times more precise than our earlier results in Ref. [17]. For Si 11+ and Ca 17+ , our IV: Binding corrections to the g-factor of the ground state of Li-like ions, in 10 −6 . The sum of all binding contributions is the difference of the atomic g factor and the free-electron g factor, ge = 2.002 319 304 361 5 (6) [27]. total g-factor values are in agreement with the previous theoretical results of the St. Petersburg group [5,16]. We note, however, that some tension exists between the calculations on the level of individual contributions. Specifically, the total electron-electron interaction correction for silicon in Ref. [5] is reported as 314.812 (3) × 10 −6 , whereas our calculation yields 314.806 (2)×10 −6 . This deviation disappears when the electron-structure correction is combined with the 1/Z QED contribution.
Table IV also compares the obtained theoretical predictions with experimental results available for two Li-like ions, Si 11+ and Ca 17+ . In both cases theoretical values deviate from the experimentally observed g factors, by 3.1 σ for silicon and 4.2 σ for calcium. The discrepancy grows with the increase of Z. Such effect could be caused by some unknown contribution missing in theoretical calculations.
The largest uncertainty in the theoretical predictions for Li-like silicon and calcium presently comes from the higherorder electron correction ∼ 1/Z 3 . This contribution cannot be calculated rigorously to all orders in Zα but needs to be treated by approximate methods or within the Zα expansion. It should be pointed out that our calculations and those by the St. Petersburg group [5,16] use different approaches for handling this correction. Our result is based on the Zα expansion, whereas the St. Petersburg group used screening potentials in the two-photon-exchange calculations and explicitly computed the three-photon-exchange contribution within the Breit approximation [5]. The Zα-expansion approach is most suitable in the low-Z region, whereas the screening-potential method is advantageous for high-Z ions. For light ions, the Zα-expansion results can be improved further, by performing the NRQED calculation of the next-order term in the α expansion.
Summarizing, we performed calculations of the two-photon exchange corrections to the g factor of the ground state of Lilike ions without an expansion in the nuclear binding strength parameter Zα. The calculations were carried out in two gauges, the Feynman and the Coulomb ones, thus allowing an explicit test of the gauge invariance. In the low-Z region, the obtained results were checked against those delivered by two different and independent methods, namely, the relativistic many-body perturbation theory with a partial inclusion of the negative-energy continuum and the nonrelativistic quantum electrodynamics. It was demonstrated that all three methods yield consistent results in the limit of small nuclear charges.
Our calculation improves the overall accuracy of theoretical predictions of the g factor of Li-like ions, especially in the low-Z region. An agreement with previous theoretical calculations is found. However, the theoretical predictions are shown to systematically deviate from the experimental results for Li-like silicon and calcium, by approx. 3 and 4 standard deviations, respectively. The reason for these discrepancies is not known at present, but is likely to be on the theoretical side. We conclude that further work is needed in order to find the reasons behind the observed discrepancies.
In this section we present explicit formulas used in the present work to numerically evaluate integrals with poles separated by the infinitesimal small additions from the integration contour. The evaluation procedure is as follows. First, use the Sokhotsky-Plemelj formula to convert integrals with poles near the integration contour to the principal-value integrals. Next, we expand the integrand in a Taylor series in the vicinity of the poles and determine the subtractions that remove singularities from the integrand. Finally, we re-add the subtractions and the perform the principal-value integrals analytically. In the Lamb-shift calculations, we encounter the integrals of three types, evaluated as follows Here, δ ∆∈(ab) is 1 if ∆ ∈ (ab) and zero otherwise. Note that the integrals in the right-hand-side of the above identities are regular, without any need to assume the principal value. In order to determine the subtractions in the integrand, we use, e.g., for Eq. (A2), that and then expanded the terms in the right-hand side in the vicinity of the corresponding poles.
Calculations for the g factor require formulas with higher powers of denominators. Such formulas can be obtained by formal differentiation of the above identities over ∆'s.

Appendix B: Matrix elements of the electron-electron interaction in Coulomb gauge
The electron-electron interaction operator in the Coulomb gauge is given by Eq. (7). The matrix element of this operator is conveniently expressed in the standard two-body-operator form that separates the angular and radial parts, with C lm l1m1,l2m2 being the Clebsch-Gordan coefficients, and R Coul L is the radial integral. Formulas for the radial integral in the Coulomb gauge were derived in Ref. [33]. For our purposes it is convenient to express them in a form similar to the Feynman-gauge radial integrals [34], where g l (ω, x 1 , x 2 ) = iω j l (ωx < ) h Furthermore, X ac,ll ′ (r) = g a (r) f c (r) S ll ′ (−κ c , κ a ) + f a (r) g c (r) S ll ′ (κ c , −κ a ) , (B6) W ac (r) = g a (r) g c (r) + f a (r) f c (r) , and the standard angular coefficients S ll ′ and C l are defined by Eqs. (A7)-(A10) of Ref. [34]. We note that the function g ret l has a finite limit at ω → 0, as the 1/ω 2 term cancels with the first term of the small-argument expansion of the spherical Bessel functions. The limiting form is The presence of the spurious singularity at ω → 0 leads to numerical instabilities in the computation g ret l at small ω, especially when evaluating the derivatives I ′ (ω) and I ′′ (ω). In order to facilitate computations for small ω, we introduce regularized functions j l and h l , separating the first term of the small-argument expansion, as follows j l (z) ≡ z l (2l + 1)!! + j l (z) , We thus obtain a regular representation for g ret l , which is suitable for a numerical evaluation, , when x 1 > x 2 . (B11) In the computation of the second derivative of g ret l over ω, we had to separate out two first terms of the the small-argument expansion of the spherical Bessel functions, in order to achieve an explicit cancelation of singular terms.
We note that when the matrix element ab|I Coul (ω)|cd is calculated with eigenfunctions of the one-particle Dirac Hamiltonian h D , it can be simplified by using the commutator relation −i α · ∇ e iωx12 = h D , e iωx12 , where [ , ] denotes commutator. In this case, we immediately have ab|I Coul (ω)|cd = α ab| 1 and thus the Coulomb-gauge matrix element is expressed in terms of the Feynman-gauge matrix elements. The above expression is convenient by its simplicity but it has a spurious singularity at ω → 0 that might lead to numerical instabilities in practical calculations. This form of the Coulomb matrix element proved to be very useful for demonstrating the gauge invariance of photon-exchange corrections [35].