Self-energy screening effects in the $g$ factor of Li-like ions

We report an investigation of the self-energy screening effects for the $g$ factor of the ground state of Li-like ions. The leading screening contribution of the relative order $1/Z$ is calculated to all orders in the binding nuclear strength parameter $Z\alpha$ (where $Z$ is the nuclear charge number and $\alpha$ is the fine-structure constant). We also extend the known results for the $Z\alpha$ expansion of the QED screening correction by deriving the leading logarithmic contribution of order $\alpha^5\ln\alpha$ and obtaining approximate results for the $\alpha^5$ and $\alpha^6$ contributions. The comparison of the two approaches yields a stringent check of consistency of the two calculations and allows us to obtain improved estimations of the higher-order screening effects.

We report an investigation of the self-energy screening effects for the g factor of the ground state of Li-like ions. The leading screening contribution of the relative order 1/Z is calculated to all orders in the binding nuclear strength parameter Zα (where Z is the nuclear charge number and α is the fine-structure constant). We also extend the known results for the Zα expansion of the QED screening correction by deriving the leading logarithmic contribution of order α 5 ln α and obtaining approximate results for the α 5 and α 6 contributions. The comparison of the two approaches yields a stringent check of consistency of the two calculations and allows us to obtain improved estimations of the higher-order screening effects.

I. INTRODUCTION
Measurements of the bound-electron g factor in light H-like ions have recently reached the fractional accuracy of few parts in 10 −11 [1][2][3]. In a combination with advanced theoretical calculations, these measurements provided the most accurate determination of the electron mass as well as one of the best tests of the bound-state quantum electrodynamic (QED) theory. Extensions of these tests towards heavier H-like ions are anticipated in the future. The main obstacle for such extensions is presently on the theory side, caused by the insufficiently known two-loop QED effects [4][5][6].
Accurate experiments were performed also on the g factors of Li-like ions [7][8][9]. They provided sensitive tests of the QED theory of the electron-correlation and relativistic nuclear recoil effects, probing QED beyond the external-field approximation. Recently, the experiments were extended further to B-like ions [10], providing the first g-factor measurement for the non-zero orbital angular momentum states. In future, a combination of the g-factor measurements in different charge states of the same element has a potential to provide an independent determination of the fine-structure constant α [11,12].
In order to match the experimental precision, theoretical investigations of atomic g factors should be performed to all orders in the nuclear binding strength constant Zα (where Z is the nuclear charge number and α is the fine-structure constant). Such calculations are often very demanding and require taking into consideration numerous effects [13,14]. A number of highly sophisticated calculations were performed during the past decade, most notably, the calculations of the self-energy and vacuum-polarization screening corrections [15,16], the two-photon exchange correction [17], and the nuclear recoil effect [18]. Despite the achieved progress, further investigations are needed in order to match the experimental precision for light Li-like ions.
In calculations performed to all orders in Zα, the electron-electron interaction is accounted for by pertur-bation theory, with the expansion parameter 1/Z. The leading term of this expansion ∝ 1/Z 0 corresponds to the hydrogenic approximation, i.e., the approximation of non-interacting electrons. The higher-order terms ∝ 1/Z 1 , 1/Z 2 , etc. are induced by the electron-electron interaction. The modification of the hydrogenic corrections by the electron-electron interaction is often referred to as the screening effect. In the present work we investigate the effect of the screening of the QED corrections, which presently induces one of the largest uncertainties in the theoretical predictions of g factors of light Li-like ions [9,19].
First calculations of the QED screening effect [20,21] included only the leading term of the Zα expansion and were applicable just for the lightest ions. The forthcoming investigations [14,22,23] approximately included contributions of higher orders in Zα, but the accuracy of these approximation was rather low, leading to errors of the screening effects ∼10% for medium-Z ions.
The first full-scale QED calculation of the self-energy and vacuum-polarization screening effects was accomplished in Refs. [15,16]. These calculations accounted for the leading screening corrections of the relative order ∝ 1/Z rigorously and the higher-order effects ∝ 1/Z 2+ approximately. Still, the numerical uncertainty of these calculations ∼1-2% was not sufficient for matching the experimental precision in the low-Z region. Moreover, the results were reported only for four ions, thus not allowing to perform a consistency check between the allorder and the Zα-expansion calculations.
The main goal of the present work is to perform an independent calculation of the self-energy screening correction for the g factor of the ground state of Li-like ions. We aim to cross-check the previously publishes results, to improve the numerical accuracy, and to perform a detailed analysis of consistency of the all-order numerical approach against the Zα-expansion calculations. To achieve this, we extend the existing Zα-expansion results by deriving the leading logarithmic contribution of order α 5 ln α and obtaining approximate results for the α 5 and α 6 contributions. Combining the two methods, we ob-tain improved estimations for the higher-order screening effects ∝ 1/Z 2+ and increase the accuracy of the theoretical description of the QED screening effects in light Li-like ions.
The relativistic units (h = c = m = 1) and the Heaviside charge units (α = e 2 /4π, e < 0) will be used throughout this paper.

II. g FACTOR
The linear Zeeman shift of the energy of an atomic state v can be written as where µ 0 = |e|/(2m) is the Bohr magneton, B = |B| is the external magnetic field, g is the g factor of the atomic state, and µ v is the angular-momentum projection on the direction of the magnetic field. In the present work we assume that the nucleus has zero spin, so that all interaction with the magnetic field comes from the electrons. The relativistic interaction of an electron with the magnetic field is represented by an operator where A(r) = (B × r)/2 is the vector potential and we choose the z axis to be directed along B. Expressing the energy shift caused by V magn in terms of the g factor and fixing the angular-momentum projection of the atomic state as µ v = 1/2, we introduce the effective operator responsible for the g factor as The matrix element of the operator V g between two Dirac wave functions is evaluated as where j and µ are the total angular momentum and its projection, respectively, C jm j1µ1,j2µ2 is the Clebsch-Gordan coefficient, and the radial integral P is given by Here, κ is the relativistic angular momentum quantum number, C L (κ a , κ b ) is the reduced matrix element of the normalized spherical harmonics (see, e.g., Eq. (C10) of Ref. [24]), and g(r) and f (r) are the upper and the lower radial components of the Dirac wave function defined as in Ref. [24].
For the point-like nucleus, the diagonal matrix element of V g with hydrogenic Dirac wave functions can be evaluated analytically as where ε v is the Dirac energy. In particular, for the case relevant for this work of v being the 2s state, where γ = 1 − (Zα) 2 .

III. GENERAL FORMULAS
We now turn to the general formulas describing the self-energy screening correction to the g factor of a Lilike ion. We will assume that the electronic configuration has the form of one valence electron state (denoted by v) over a closed shell of core electron states (denoted by c). The derivation of the formulas was first presented in Ref. [16] within the formalism of the two-time Green function method [25]. In the present work, we will reformulate this problem in order to suit our calculational approach.
We start with introducing two operators which will be building blocks in the following formulas. The first one is the electron-electron interaction operator I(ω), defined as where α µ = (1, α) are the Dirac matrices, r 12 = r 1 − r 2 , and D µν (ω, r 12 ) is the photon propagator. In the present work we use the Feynman gauge, in which the photon propagator takes the simplest form, where r 12 = |r 12 | and ǫ is a positive infinitesimal addition.
The one-loop self-energy (SE) operator Σ(ε) is defined by its matrix elements with the one-electron wave functions |a and |b , where the sum over n is carried out over the complete spectrum of the Dirac equation (implying the summation over the discrete part of the spectrum and the integration over the continuum part of the spectrum) and u = 1 − iǫ. We will split the total self-energy screening correction into four parts as ∆g sescr = ∆g po + ∆g vr,Zee + ∆g vr,scr + ∆g dvr , (11) with the individual contributions defined in the remaining of this Section.

A. Perturbed-orbital SE contribution
The perturbed-orbital SE contribution incorporates all terms that can be expressed as matrix elements of the one-loop SE operator Σ(ε). It can be represented as a sum of two parts, where the first part contains matrix elements of the SE operator with a perturbed wave functions on one side, whereas the second term has perturbed wave functions on both sides. The first term can be expressed as where and where a, b, c, and d are the one-electron states of the core or the valence electron.
Here and in what follows, ∆ ab = ε a − ε b , the prime on the summation symbol means that terms with vanishing denominator should be omitted from the summation, and each prime in I ′ (ω) and I ′′ (ω) denotes the derivative over the energy argument. The summation over µ core runs over the angular-momentum projections of the core electron states, µ core = ± 1 /2 for the (1s) 2 shell. The second term in Eq. (12) is represented by where µ c denotes the angular-momentum projection of the core electron state c, the perturbed wave functions are defined by and (ab) = (vc) or (cv).
In the notations of Ref. [16], ∆g po1 corresponds to the sum of the A, E, and G terms, and ∆g po2 corresponds to the B term. Formulas (12)- (17) were derived in Ref. [16] by the two-time Green's function method [25]. They can be also obtained by the standard Rayleigh-Schrödinger perturbation theory as demonstrated in Appendix A.

B. Perturbed Zeeman-vertex contribution
The perturbed Zeeman-vertex contribution incorporates terms that can be expressed as non-diagonal matrix elements of the Zeeman vertex operator plus the corresponding reducible part. It is given by where V g,aa ≡ a|V g |a , Σ ′ (ε) denotes the derivative of the self-energy operator over the energy argument ε, and the perturbed wave function is defined as with (ab) = (vc) or (cv). The matrix element of the Zeeman vertex operator (with the corresponding reducible part) is given by In the notations of Ref. [16], ∆g vr,Zee corresponds to the sum of the C1+H1 terms and a part of the H3 term.

C. Perturbed screened-vertex contribution
The perturbed screened-vertex contribution is a part that can be expressed in terms of non-diagonal matrix elements of the screened (i.e., two-electron) vertex operator plus the corresponding reducible part. We represent it as a sum of the vertex and the reducible parts, ∆g vr, scr = ∆g ver, scr + ∆g red, scr .
In the notations of Ref. [16], ∆g ver, scr corresponds to the sum of the C2+H2+F terms, and ∆g red, scr corresponds to a part of the H3 term.

D. Double-vertex contribution
The double-vertex contribution is comprised of the matrix element of the double-vertex operator plus the corresponding reducible parts, all of them containing the third power of ω in the denominator. It is represented as where the operator Λ dvr consists of four parts, The first term in the sum is the double-vertex operator, which is defined by its matrix element as .
The second term is the derivative of the screened-vertex operator, whose matrix element is The third term is the derivative of the Zeeman-vertex operator, which is .
The last term in Eq. (28) is the second derivative of the SE operator, In the notations of Ref. [16], the four terms in the right-hand-side of Eq. (28) correspond to the D, I2, I1, and I3 terms, respectively.

IV. DIVERGENCIES
General formulas for the individual contributions presented in the previous Section contain divergencies, both of the ultraviolet (UV) and infrared (IR) kind. The UV divergencies appear in contributions containing the first and the second power of ω in the denominator(s) inside the radiative photon loop. According to the standard procedure [26], UV divergencies are covariantly regularized by isolating one or two first terms of the expansion of the bound-electron propagators in terms of the interaction with the binding nuclear field. These terms are calculated in momentum space within the dimensional regularization, whereas the remainder is calculated in coordinate space using the partial-wave expansion of the bound-electron propagators. The UV divergencies are identified in terms of one-loop renormalization constants and cancelled when all individual contributions are added together. The cancellation of UV divergencies was demonstrated in Ref. [16] and does not need to be repeated here. In practical calculations, it is sufficient just to replace the free SE operator and the free one-loop vertex operator by their renormalized expressions.
We now turn to the IR divergencies, which have not been discussed in detail in Ref. [16]. These divergencies occur when the denominators of the electron propagators inside the radiative photon loop vanish at ω → 0. As we will show below, the IR divergencies originate from terms of the form with β ≥ 2. In the present investigation, we will encounter IR divergent terms with β = 2 and β = 3.
It should be noted that the term with β = 1 (appearing, e.g., in the one-loop SE matrix element) is IR safe.
In order to show this, it is sufficient to rotate the left half of the ω integration contour on the right half-axis, where the small addition −i0 indicates that this part lies on the lower bank of the cut of the photon propagator. On the upper bank of the cut of the photon propagator, √ ω 2 = ω, whereas on the lower bank, √ ω 2 = −ω. Therefore, which is obviously converging at ω → 0. In order to evaluate the IR divergent integrals J 2 and J 3 , we regularize the divergencies by introducing a fi-nite photon mass µ in the photon propagator, evaluate the integral over ω analytically, and separate out the µdependent divergent terms, as described in Ref. [27]. The results for the IR divergent integrals (omitting terms vanishing in the limit µ → 0) are given by where γ is Euler's constant. We now demonstrate the cancellation of the IR divergencies in the sum (11). It is convenient to express the IR-divergent parts of individual contributions in the form where i runs over the contributions described in Sec. III. The perturbed-orbital SE contribution (12) does not contain any IR divergences. In the perturbed Zeeman-vertex contribution (18), IR divergencies intrinsically present in the vertex and reducible parts cancel each other, so that the total expression is finite and does not require a separate treatment. The other contributions in Sec. III contain IR divergencies, identified as follows: It can be easily seen that the sum of all IR contributions (38) -(43) vanishes. In actual calculations, the IR-divergent contributions were isolated by introducing point-by-point subtractions in the integrand and then evaluated analytically according to Eqs. (35) and (36). Specifically, the matrix element of the double-vertex operator (29) is represented as (without the IR part accounted for by Eq. (40)) where a ′ and c ′ (c ′′ ) denote the electron states that differ from a and c only by the angular-momentum projection, µ a ′ and µ c ′ (µ c ′′ ), correspondingly.

V. COMPUTATION OF INDIVIDUAL CONTRIBUTIONS
A. Perturbed-orbital SE contribution The calculation of the perturbed-orbital SE contributions, given by Eqs. (13) and (16), is naturally reduced to a computation of non-diagonal matrix elements of the SE operator. In the present work, we use the numerical approach developed in Ref. [28], which has an important advantage of a rapid convergence of the partial-wave expansion. The perturbed wave functions in Eqs. (13) and (16) were calculated with help of the finite basis-set for the Dirac equation constructed with B-splines [29]. We do not use the dual kinetic balance (DKB) method [30] in the present work, since our calculations are performed with the point nuclear model, for which the DKB approach is not applicable.
The calculation of the perturbed-orbital SE contributions is simplified by the fact that the matrix element of the SE operator is diagonal in the relativistic angular-momentum quantum number κ and the angularmomentum projection µ of the external wave functions, where (. . .) does not depend on the angular-momentum projections. Therefore, Eq. (13) involves the perturbed wave functions of just one angular symmetry. Eq. (16) contains a summation over several angular symmetries of the perturbed wave functions (three in the general case but just two in our case). We also note that the perturbed wave functions |δ po a and |δ scr a contain imaginary parts which contribute when combined with the imaginary part of the SE operator. For actual calculations of Eq. (16), it is convenient to move the summation over µ c into the definition of one of the perturbed wave functions. For the first matrix element in Eq. (16), this can be done immediately. In order to do this in the second matrix element, we fix the angular-momentum projection of the c state in the magnetic perturbed wave function as µ c = 1/2, and move the summation over µ c into the definition of |δ scr c , together with the appropriate factor s µc that carries the dependence of the magnetic matrix element n|V g |c on µ c , where s µc = (−1) µc−1/2 C 10 jnµc,jc−µc C 10 In this way, we reduce the number of matrix elements of the SE operator to be computed by half.
In order to check our numerical procedure for computation of the perturbed wave functions in Eqs. (13) and (16), we replace the self-energy operator by the vacuumpolarization potential, Σ(ε) → V VP , thus reproducing the corresponding vacuum-polarization screening corrections, calculated previously in Refs. [16,23].

B. Perturbed Zeeman-vertex contribution
The perturbed Zeeman-vertex contribution, defined by Eq. (18), can be considered as a non-diagonal generalization of the Zeeman-vertex correction for the hydrogenlike atom, specifically, ∆g vr given by Eq. (12) of Ref. [31]. The only difference of Eq. (18) as compared to the hydrogenic case is that the reference-state wave function on the right-hand side of the matrix element is replaced by the perturbed wave function |δa given by Eq. (19). It should be mentioned that the perturbed wave function has contributions from several angular symmetries (three in the general case but only two in our case of κ v = κ c = −1) and an imaginary part, which contributes to the final result.
Similarly to the perturbed-orbital SE contribution, it is convenient to move the summation over the angularmomentum projection of the core electron µ c in Eq. (18) into the definition of the perturbed wave function. For the first matrix element v| . . . |δv , it can be done immediately. For the second matrix element c| . . . |δc , some manipulations are needed. Specifically, we observe that the dependence of the vertex matrix element on the angular-momentum projections of the external wave functions can be factorized out as where (. . .) does not depend on the angular-momentum projections. Therefore, we can fix the momentum projection µ c = 1/2 in the matrix element and redefine the perturbed wave function as where s µc is given by Eq. (48). The numerical evaluation of the perturbed Zeemanvertex contribution is similar to the calculation of the diagonal matrix elements for the hydrogenic atoms, described in details in Refs. [31,32]. Specifically, the whole contribution is separated into three parts, ∆g vr,Zee = ∆g (0) vr,Zee + ∆g (1) vr,Zee + ∆g where the superscript indicates the number of interactions with the binding Coulomb potential in the electron propagators. The first two terms in the right-handside of the above equation are evaluated in momentum space, without any partial-wave expansion. Only the last term containing two and more interactions with the Coulomb field is calculated in coordinate space. Thanks to the separation of the one-potential term ∆g (1) vr,Zee , the partial-wave expansion of the remainder ∆g (2+) vr,Zee converges rapidly and can be calculated to high accuracy. The contributions ∆g (0) vr,Zee and ∆g (1) vr,Zee need some generalization as compared to the diagonal case because of different angular symmetries of the perturbed wave function; the corresponding formulas were already derived in our calculation of the self-energy correction to the magnetic shielding [33,34].
We note that as long as we calculate the vertex and the reducible part together and use an appropriate ωintegration contour (consisting of the low-and highenergy parts), the integrand in Eq. (20) has a smooth small-ω behaviour. The would-be IR divergences in the vertex and reducible parts are cancelled numerically at a given ω in this approach. Alternatively, the contribution with ε n1 = ε n2 = ε a can be separated out and evaluated analytically with help of formulas from Sec. IV. We checked the equivalence of both methods in order to test the consistency of our numerical procedure.
We mention here some of the further cross-checks of the numerical procedure made in order to eliminate possible errors: (i) by replacing |δa → |a in Eq. (20) we reproduced known results for the vertex and reducible diagonal matrix elements for the 1s and 2s hydrogenic states [31]; (ii) by replacing Λ Zee (ε a ) + V g,aa Σ ′ (ε a ) → V g in Eq. (20) we reproduced known results for the onephoton exchange correction to the g factor, both in coordinate and momentum space.
The free vertex part contains only free electron propagators; it is renormalized and calculated in momentum space. The many-potential part contains one and more interactions with the binding Coulomb field; it is calculated in coordinate space using the partial-wave expansion of the electron propagators. For the reducible contribution, we separate the zero-potential and one-potential contributions, ∆g red,scr = ∆g (0) red,scr + ∆g (1) red,scr + ∆g The zero-and one-potential contributions are calculated in momentum space. The separation of the one-potential contribution improves the convergence of the partialwave expansion in the many-potential reducible contribution.
The computation of the many-potential part is very similar to that for the Lamb-shift case. We introduce perturbed wave functions defined by Eqs. (17) and (26) and calculate them by using the finite basis-set method [29]. Note that several different symmetries of the perturbed wave functions contribute to the the final result (three in the general case but just two in our case).
Contrary to the many-potential part, the evaluation of the free part turned out to be different from our previous calculations for the Lamb shift [35,36]. The difference is that for the Lamb shift, analytical formulas for the basic angular integrals were derived using averaging over the angular-momentum projections of the valence state. In the g-factor calculations, the angular-momentum projection of the valence-electron state is fixed. Moreover, we need to account for the case when the angular symmetry of the perturbed wave function is different from the angular symmetry of the reference state. For this reason, we developed a generalized procedure for performing angular integrations in the momentum space. The evaluation of the free screened-vertex matrix elements is described in Appendix C.
Our numerical calculations of the many-potential vertex and reducible parts were performed using the analytical representation of the Dirac-Coulomb Green function in terms of the Whittaker functions. Integrations over the radial variables were carried out by the numerical approach described in detail in the recent review [27]. The partial-wave expansion was extended up to |κ| = 30; the remaining tail of the expansion was estimated by a polynomial fitting in 1/|κ|. We mention here several cross-checks of the computational procedure made in order to eliminate possible errors: (i) we checked that in the diagonal case, our calculations reproduce known results for the Lamb shift [35]; (ii) we also checked that by replacing the radiatively corrected vertex by the plain vertex we reproduce known results for the one-photon exchange correction to the g factor, both in coordinate and momentum space. Specifically, as a part of this test, we checked that the replacement Γ R,µ → γ µ in Eq. (C1) yields the matrix element of the electron-electron interaction operator I(∆) in the coordinate-momentum representation.

D. Double-vertex contribution
The computation of the double-vertex contribution is the most complicated part of the calculation. We start our discussion with the last three terms in the right-hand side of Eq. (28). These terms are induced by derivatives of the Zeeman-vertex, screened-vertex, and self-energy operators. Each of these operators were already examined, so we need only to evaluate the derivative. In actual calculations, we find it convenient to convert the derivative ∂/(∂ε a ) to ∂/(∂ω) and to apply integration by parts, moving the derivative to the photon propagator. Specifically, we use the following identities, and ∞ −∞ dω ∂ 2 ∂ 2 ε a n =a The advantage of using the above formulas is that the derivative over the photon propagator can be easily evaluated analytically, contrary to the derivative over the electron propagator. An additional bonus from this transformation is that the behaviour of the transformed integrand is smoother in the low-ω region. All double-vertex contributions (29)-(32) contain the third power of ω in the denominator and thus are convergent in the ultraviolet region. Therefore, in principle, one does not need to separate out the zero-potential contributions in them. However, we find it advantageous to do so, since this subtraction improves the convergence of the partial-wave expansion drastically in the low-Z region.
Specifically, we separate out the contributions of the free-electron propagators from the last three terms on the right-hand side of Eq. (28) and evaluate them in momentum space, without any partial-wave expansion. Formulas for these zero-potential contributions are easily obtained by differentiation of the corresponding expressions for the zero-potential Zeeman-vertex, screenedvertex, and self-energy operators. A comparison of the numerical results obtained in the high-Z region with and without the subtraction was employed as a useful check of our numerical procedure. In the low-Z region, the convergence of the partial-wave expansion becomes increasingly slower. Even with the subtraction, we had to extend the partial-wave expansion up to |κ max | = 60, in order to reach the desired numerical accuracy. The computation was carried out with help of the analytical representation of the Dirac-Coulomb Green function [27].
The evaluation of the matrix elements of the doublevertex operator (29) is the most computationally intensive part. The straightforward approach is to compute it as it stands, after separating IR divergencies according to Eq. (44). However, this turns out to be applicable only in the high-Z region; for lower Z, the convergence of the partial-wave expansion becomes excruciatingly slow. The standard way to accelerate the convergence would be to separate out the zero-potential double-vertex contribution and to calculate it in momentum space, which is a very cumbersome task. Fortunately, this is not really needed. It turns out that the convergence of the partial-wave expansion can be greatly accelerated if one separates out only a relatively simple part of the full zeropotential contribution; this is the approach used in the present work.
Specifically, we introduce the subtraction term ∆g (s) dver which we obtain from ∆g dver by applying the following prescriptions: (i) all bound-electron propagators are replaced by the free electron propagators, (ii) only the direct contribution is taken (P vP c = vc and QvQc = vc), (iii) only the Coulomb part of the electron-electron interaction is taken (I(∆) → α/r 12 ). We note that under these restrictions the contribution in which the radiative loop is attached to the core-electron line vanishes after the summation over the angular-momentum projections of the core electrons, so that only the contribution with the radiative loop attached to the valence electron line survives.
The subtraction term defined in this way is represented in momentum space as (see Sec. III of Ref. [31]) ∆g (s) dver = 4im dp 1 dp 2 dl (2π) 6 V core (l) where q = p 1 − l − p 2 , V core is the Fourier-transformed potential of the charge density of the core electrons, and Λ (0) dver is the free double-vertex operator defined as Here, p 1 , l, and p 2 are 4-vectors with a fixed time component, p 1 = (ε v , p 1 ), l = (ε v , l), p 2 = (ε v , p 2 ), p / = p µ γ µ , and γ µ = (γ 0 , γ) are the Dirac matrices. We now observe that the subtraction term (56) can be obtained from the one-potential vertex contribution given by Eq. (43) of Ref. [31], by replacing the nuclear Coulomb potential V C (q) = −4πZα/q 2 with the potential of the core charge density V core (q). We, therefore, just use the formulas derived in Sec. IIIB of Ref. [31] in order to compute the double-vertex subtraction contribution (56) in momentum space.
The remainder represented by the difference ∆g dver − ∆g (s) dver was calculated in coordinate space using the partial-wave expansion of the electron propagators. The low-energy part of the remainder was computed using the finite basis-set method. The high-energy part of the remainder was evaluated with the analytical representation of the Dirac Green function. The radial integrations were computed by the numerical approach described in detail in the review [27].
The number of partial waves included into the computation varied from |κ max | = 50 in the high-and medium-Z region to |κ max | = 75 for Z ≤ 10. In order to crosscheck our numerical procedure, we calculated the lowenergy part of the remainder in two ways, using the analytical representation of the Green function and the finite-basis set B-spline representation.

VI. NRQED EXPANSION
Within the nonrelativistic quantum electrodynamics (NRQED), the bound-electron g factor of light atoms is represented by an expansion in powers of the finestructure constant α. The leading QED contribution to the bound-electron g factor is of order α 3 and comes from the anomalous magnetic moment of the electron. The corresponding contribution was derived by Hegstrom [38] and calculated in Refs. [19][20][21]. In the present work we calculate the leading logarithm of the next-order contribution, α 5 ln α.
The α 5 ln α QED contribution can be easily obtained as a straightforward generalization of the corresponding result for the hydrogenic ions derived in Ref. [39]. So, the QED contribution to the bound-electron g factor of light atoms is were a numerates the electrons in the atom, the operator Q a is given by Eq. (7) of Ref. [19], and . F denotes the so-called Fermi (spin-factorized) radial matrix element, defined for an arbitrary one-electron operator H as (see Ref. [19] for details) where ψ is the full wave function of the reference state (i.e., the antisymmetrized product of the spacial and the spin functions), φ is the spacial part of the wave function, and S = a σ a is the spin operator. The matrix elements of Q a were calculated for the ground state of Li-like ions in Ref. [19] and fitted to the 1/Z-expansion form as The Fermi matrix element of the δ function are calculated in the present work and listed in Table I. We represent the results in a form of the 1/Z expansion as a δ 3 (r a ) The result (59) can be extended further by observing that all terms in the α expansion for hydrogenic atoms proportional to the expectation value of the δ function can be immediately generalized to the many-electron case by the substitution δ 3 (r) → a δ 3 (r a ) F . We, therefore, obtain the following approximate representation for the radiative QED contribution to the bound-electron g factor of light atoms, where a 40 and a 50 are the hydrogenic coefficients for the 2s state, a 40,se (2s) = −10.707 716 [39], a 40,vp (ns) = −16/15 [40], a 50,se (ns) = 23.282 005 [41], and a 50,vp (ns) = 127π/216 [40,42], with subscripts "se" and "vp" labelling contributions originating from the self-energy and vacuum-polarization, correspondingly. The hydrogenic coefficient a 40 (ns) is weakly ndependent (a 5% difference between the 1s and 2s states), so we expect the unknown screening contribution to it to be within 10-20%. The next-order coefficient a 50 is nindependent, so the corresponding result for few-electron atoms is exact.

VII. RESULTS AND DISCUSSION
Numerical results of our all-order calculation of the self-energy screening correction to the g factor of the ground (1s) 2 2s state of Li-like ions are presented in Table II. Our calculation is performed for the point nucleus. We note significant numerical cancellations between individual contributions, present throughout the whole Z region. For Z = 82, our result is in perfect agreement with that of Ref. [16] but significantly more accurate.
In order to analyze our all-order numerical results for the self-energy screening correction, it is convenient to separate out its leading α and Z dependence, by introducing the function G(Zα) as follows ∆g sescr = α 2 (Zα) G(Zα) .
The Zα expansion of the function G follows from the results obtained in Sec. VI, where the leading coefficient g 30 = −274/(2187π) = −0.039 880 . . ., g 51 = −2.6557/(8π) × (32/9) = −0.3757, g 50 ≈ −2.6557/(8π) × (−10.708), and g 60 = −2.6557/(8π) × 23.282. Fig. 1 shows the comparison of the all-order numerical results for the scaled function G(Zα) with the predictions based on the Zα-expansion (65). We conclude that our all-order results converge to the prediction of the Zα-expansion as Z → 0. We also observe that the inclusion of the approximate higherorder contributions g 50 and g 60 significantly improves the agreement between the all-order and Zα-expansion results. The deviation of the all-order results from the Zα-expansion in the low-Z region is consistent with an additional contribution to g 50 , δg 50 ≈ 0.2. So far we addressed the self-energy screening correction of the relative order 1/Z as compared to the leading, one-electron contribution. The complementary vacuumpolarization screening correction was calculated previously in Ref. [16] and recently reproduced in Ref. [23].
In order to complete the calculation of the QED screening effects, we need also to estimate the contribution of the higher-order screening ∝ 1/Z 2+ . For low-Z ions, this can be done immediately with help of the NRQED formula (63). Such approach, however, would result in large uncertainties for high-Z ions. In order to avoid this, we devised an estimate which is equivalent to the one obtained from Eq. (63) in the low-Z region but applicable also for high-Z ions. Specifically, we estimate the QED screening correction of order 1/Z 2 as We also devise an alternative estimation of the 1/Z 2 QED screening contribution, based on the one-electron QED results. Specifically, we introduce the 1/Z 0 higherorder remainder function H (0) , by representing the allorder results for the one-electron self-energy and vacuumpolarization of the 2s state as After that, our second approximation for the 1/Z 2 QED screening correction is obtained as where c 0 = 1 is the first coefficient of the 1/Z expansion (62).
The QED screening effects of order 1/Z 3 and higher are relevant only for the lightest ions and can be accounted for in the leading order of the α expansion,  where the matrix element contains all terms of its 1/Z expansion starting with 1/Z 3 , The column "SE(fns)" presents our results for the shifts of the point-nucleus correction due to the finite nuclear size. It was calculated by taking the corresponding hydrogenic correction for the 2s state, obtained as in Ref. [43], and scaling it by δ 3 (r) → a δ 3 (r a ) F . We assumed an uncertainty of 25% of this approximation. This correction is relevant only for high-Z ions; its uncertainty is negligible on the level of the total theoretical uncertainty.
The column "VP" lists results for the vacuumpolarization screening correction of the relative order 1/Z, as calculated in Ref. [23]. The uncertainty of this part comes from uncalculated higher-order contributions in the so-called magnetic-loop vacuum-polarization. It was estimated by multiplying the results for the magnetic-loop correction obtained in the free-loop approximation (see Ref. [23] for details) by the factor of 2 (Zα) 2 .
The column "h.o." presents results for the QED screening correction of the relative order 1/Z 2+ . They were obtained as a half-sum of the two approximations given by Eqs. (66) and (69) plus the higher-order screening contribution given by Eq. (70). The quoted uncertainty was obtained as twice the difference between the two approximations.
In Table III we compare our final theoretical values for the QED screening effect with previous results obtained by Glazov, Volotka and co-authors [9,17]. We mention that our present approach is equivalent to theirs for the 1/Z part of the QED screening effects but differs in estimating the higher-order 1/Z 2+ screening contributions. In our approach, we base our estimate on the NRQED results, whereas Refs. [9,17] approximately accounted for the higher-order screening by using various screening potentials. We observe the general consistency of the two calculations. This consistency is a rather strict test, because of delicate numerical cancellations between the individual terms required to get the final result. Still, the estimated error bars of the two calculations do not overlap in most cases, the deviation being on the level of 2-3 σ. This indicates that further work is needed in order to fully cross-check the QED screening calculations.
In Table IV we present a compilation of all known binding corrections to the g-factor of the ground state of Lilike silicon, 28 Si 11+ . As compared to the analogous table in our previous investigation [19], the one-loop QED effects of the relative order 1/Z and 1/Z 2+ and the nuclear recoil contributions were updated. The screened QED effects are calculated in the present work, whereas the nuclear recoil corrections of relative orders 1/Z 0 , 1/Z 1 , and 1/Z 2+ were calculated by Shabaev et al. [18]. It should be mentioned that Shabaev et al. found a mistake in the previous calculations of the recoil effect [20,21], which resulted in a small shift of theoretical values for this correction.
Our final theoretical value presented in Table IV is in agreement with the previous theoretical results of Glazov et al. [9] and Volotka et al. [17]. However, our result disagrees with the recent experimental value [9] by about five standard deviations. It is interesting that the theoretical predictions, with time, are moving away from the experimental result, while steadily increasing in estimated accuracy.
Considering the comparison of the present theoretical prediction for silicon with the latest theoretical result of Glazov et al. [9], we note two small deviations. One is the difference of +0.005 (2) × 10 −6 in the QED screening effect, whereas the other is the difference of −0.003 (3) × 10 −6 in the 1/Z 3+ electron-correlation effect. These two differences partially cancel each other, resulting in the total difference of the two total theoretical values of +0.002 (4) × 10 −6 , well within the quoted error bars.
Commenting on the disagreement with the experimental result for silicon, we note that not all effects in the theoretical prediction has been confirmed by independent calculations so far. Apart from the already mentioned small inconsistencies between the two calculations of the QED screening correction and the 1/Z 3+ electron-correlation effect, the two-photon exchange QED effect has so far been calculated by one group only [17]. The other contributions in the theoretical prediction seem to be under a better control. In particular, the nuclear recoil effect was recently calculated rigorously within QED by Shabaev and co-workers [18]. They found a mistake in the earlier calculation by Yan [20,21]. After correcting this mistake, the extrapolation of Yan's results (as in Ref. [19]) yields a result for silicon that agrees with that of Shabaev and co-workers up to few parts in 10 −10 . So, this effect cannot be responsible for the deviation from the experimental value.

SUMMARY
We performed a calculation of the self-energy screening effects for the g factor of the ground state of Li-like ions. The contribution of the relative order 1/Z was calculated rigorously within QED, to all orders in the binding nuclear strength parameter Zα. The higher-order screening contribution ∝ 1/Z 2+ was calculated approximately, using results for the coefficients of the α expansion obtained in this work in the framework of nonrelativistic QED. In the result, we were able to improve the theoretical accuracy of the QED screening effects in light Li-like ions and, as a consequence, the accuracy of the theoretical prediction of the g factor of Li-like silicon. The total theoretical result for Li-like silicon is in agreement with previous theoretical calculations but differs by about five standard deviations from the experimental result. We conclude that further theoretical work is needed in order to investigate small deviations between different theoretical calculations and a much larger discrepancy with the experimental result.   − n ′ a|V g |n n|U |a ab|I|cd (ε a − ε n ) 2 + a|V g |n nb|I|cd a|U |a (ε a − ε n ) 2 + a|U |n nb|I|cd a|V g |a (ε a − ε n ) 2 . (A4) The second part Λ red (abcd) contains terms with derivatives of the electron-electron interaction, Λ red (abcd) = n ′ a|V g |n nb|I ′ |cd ε a − ε n d|U |d − b|U |b + n ′ a|U |n nb|I ′ |cd ε a − ε n d|V g |d − b|V g |b + n ′ a|U |n n|V g |a ε a − ε n ab|I ′ |cd + 1 4 ab|I ′′ |cd d|V g |d − b|V g |b d|U |d − b|U |b .
The integrals of three and four spherical harmonics R 3 and R 4 are evaluated in terms of Clebsch-Gordan coefficients by standard formulas [47]. For specific cases, the angular coefficients can be evaluated analytically, as it was done in Ref. [35]. In the present work we prefer to evaluate all sums of Clebsch-Gordan coefficients numerically. The specific cases of J = 0 and J = 1 were simplified and evaluated separately.
We now consider the matrix element of the zero-potential two-electron vertex operator for the Coulomb gauge of the exchanged photon. In this case the 4-vector potential A µ = (A 0 , A i ) becomes In order to perform the angular integrations in the momentum integrations, we use the following representation where the functions R i are expressed in terms of functions R i as follows: R ac 3 = 1 q 2 p 2 2 R ac 1 − p 2 1 R ac 2 + R ac 3 p 2 1 − p 12 + R ac 4 p 12 − p 2 2 , (C26) R ac 5 = 1 q 2 − p 2 1 R ac 1 + p 2 2 R ac 2 + R ac 5 p 2 1 − p 12 + R ac 6 p 12 − p 2 2 . (C27) We observe that the matrix element in the Coulomb gauge involves the same angular coefficients (C14)-(C16) as in the Feynman gauge.