The Lamb shift of the $1s$ state in hydrogen: two-loop and three-loop contributions

We consider the $1s$ Lamb shift in hydrogen and helium ions, a quantity, required for an accurate determination of the Rydberg constant and the proton charge radius by means of hydrogen spectroscopy, as well as for precision tests of the bound-state QED. The dominant QED contribution to the uncertainty originates from $\alpha^8m$ external-field contributions (i.e., the contributions at the non-recoil limit). We discuss the two- and three-loop cases and in particular, we revisit calculations of the coefficients $B_{61}, B_{60}, C_{50}$ in standard notation. We have found a missing logarithmic contribution of order $\alpha^2(Z\alpha)^6m$. We have also obtained leading pure self-energy logarithmic contributions of order $\alpha^2(Z\alpha)^8m$ and $\alpha^2(Z\alpha)^9m$ and estimated the subleading terms of order $\alpha^2(Z\alpha)^7m$, $\alpha^2(Z\alpha)^8m$, and $\alpha^2(Z\alpha)^9m$. The determination of those higher-order contributions enabled us to improve the overall accuracy of the evaluation of the two-loop self-energy of the electron. We investigated the asymptotic behavior of the integrand related to the next-to-leading three-loop term (order $\alpha^3(Z\alpha)^5m$, coefficient $C_{50}$ in standard notation) and applied it to approximate integration over the loop momentum. Our result for contributions to the $1s$ Lamb shift for the total three loop next-to-leading term is $(-3.3\pm10.5)(\alpha^3/\pi^3)(Z\alpha)^5m$. Altogether, we have completed the evaluation of the logarithmic contributions to the $1s$ Lamb shift of order $\alpha^8m$ and reduced the overall $\alpha^8m$ uncertainty by approximately a factor of three for H, D, and He$^+$ as compared with the most recent CODATA compilation.


I. INTRODUCTION
Already for a few years, there exists a discrepancy in the determination of the proton charge radius by means of the spectroscopy of ordinary and muonic hydrogen (see, e.g., [1,2]), commonly known as the proton radius puzzle. There are different contributions to the uncertainty of the determination of the proton radius by those methods. The largest uncertainty originates from the hydrogen spectroscopy and a serious experimental activity in this direction is in progress (see, e.g., [3][4][5]). The second largest uncertainty comes from the Quantum electrodynamics (QED) theory of the 1s Lamb shift in hydrogen [2]. There are a few theoretical problems which require clarification. They relate to two-loop and threeloop radiative corrections. Some higher-order contribu- * Electronic address: savely.karshenboim@mpq.mpg.de tions have not been cross-checked, and some not studied at all. In particular, the two-loop contributions of order α 2 (Zα) 5 m [6,7] are well established, while at the next order in Zα the contributions for the virtual light-by-light scattering have not been studied properly (see, e.g., a discussion on a previously missed term in [8]). Meanwhile, the results for the pure self-energy contribution of order α 2 (Zα) 6 m [9,10] are to some extent controversial (see, e.g., [2]). One more challenge is related to the next-toleading order three-loop contribution (order α 3 (Zα) 5 m); the existing estimation [2] does not have a solid ground.
Besides the proton radius puzzle, an improvement of the theoretical prediction of the 1s Lamb shift is essential for the determination of the Rydberg constant [1,2], precision tests of the bound-state QED, constraints on light neutral particles, such as a dark photon from physics of simple atoms (see, e.g., [11]), and interpretation of the currently ongoing 1s − 2s He + experiments [12,13].
The Lamb shift of the atomic energy levels is a QED effect, that can be experimentally studied in light hydrogen-like atoms with high accuracy (cf. [14]). The theoretical prediction of this phenomenon involves the values of the input parameters, such as the Rydberg constant and the proton charge radius, that limit the accuracy of the calculations. A separate input to the uncertainty originates in the computation of various highorder QED effects. The dominant contributions to the QED error budget come from the radiative corrections in the external-field approximation. We follow the standard convention and parametrize these corrections as (see, e.g., [2,15]) ns (Zα) corresponds to the i-loop radiative insertions and the relevant contributions are at the one-, two-, and three-loop level. The four-loop contributions are neglected in (1). The uncertainty due to the unknown leading four-loop term, which is expected at the level of a few units of α 4 /π 4 (Zα) 4 m, is essentially below the uncertainty of the higher-order two-loop and three-loop terms. The latter are at the level of ten units of α 2 /π 2 (Zα) 6 m and α 3 /π 3 (Zα) 5 m, respectively (see below).
Theory of the one-loop contributions is firmly established (see [2,15] for details). The largest and most important contribution, related to the electron self-energy, has been calculated directly for Z = 1, 2 [16], i.e., for H and He + . We consider below the two-loop and three-loop radiative corrections.
The functions F (i) can be expanded at low Zα and at two and three loops, the results read Here, we focus on the 1s state and the F (i) coefficients are always meant to be related to the aforementioned, 1s state. It is not clear a priori which logarithmic terms are present in (2) and (3). Sometimes a special study is required. For example, it was believed [2] until recently that C 63 = 0, while the presence of B 72 = 0 was rather disputable. Both issues have been recently resolved in [17] and we discuss it also below.
A number of the two-loop (B ... ) and three-loop (C ... ) coefficients have been known with a sufficient accuracy. These include B 40 , B 50 , B 63 , B 62 , B 61 , and C 40 . Estimations with a credible uncertainty have also been available for B 60 , C 50 , and C 63 . A concise summary concerning all these coefficients can be found in [2]. Some of the corrections have been revisited since publication of [2]. These include, e.g., B 61 [8] and B 72 , C 63 , and C 62 [17].
In the case of the two-loop corrections, rather than B 60 , we use G 60 (Zα) defined as i.e., it is equal to B 60 with all the higher-order (in Zα) corrections included. G 60 (Zα) is more appropriate if one uses results of numerical calculations.  (34) 0 −0.01 We begin with the two-loop logarithmic coefficient B 61 . First calculated in [18,19], the result was applied in [2]. After original publication, the diagrams with the light-by-light (LbL) scattering block (see Fig. 1a) have been studied, and a correction to the previous result was found [8] due to the LbL diagrams overlooked in [18,19]. The LbL contributions are the most difficult for the numerical calculations. The analytic calculations have been available only for the order α 2 (Zα) 5 m [7] and absent for a while for the next order in Zα until the publication of [8]. Those diagrams receive a contribution from soft photons responsible for the appearance of a long-distance potential. After integrating out the hard modes (i.e. with momenta comparable with m), effective local operators appear, which give rise to the two-photon vertices shown in diagrams b and c in Fig. 1. The remaining photons are soft (i.e. their momenta are much smaller than m). There are two possible soft pairs of photons, those which connect the nucleus and the electron loop (see Fig. 1b) and those which connect the electron line and the electron loop (see Fig. 1c). The former case was covered by [8], while we consider the latter here. The bottom part of the diagram in Fig. 1a, i.e., the electron loop in the Coulomb field of the nucleus, is the known virtual Delbrück scattering amplitude (see, e.g., [20,21] and references therein). The upper part, i.e., the electron line and two soft photons connecting the electron line and the electron loop, can be drastically simplified within the soft-photon kinematics, where the energy transfer (q 0 ) is comparable with the momentum transfer (q) and Zαm ≪ |q 0 | ∼ |q| ≪ m. The integral over q 0 simplifies considerably within a kind of staticelectron kinematics discussed in details in [22].
The resultant integral induces an effective potential that behaves as r −4 (cf. [23], see also [8]). In the 1s state, the expectation value of this potential diverges at the short end of the interval 1/m ≪ r ≪ 1/(Zαm). The manifestation of the divergence in the perturbation theory is a logarithmically enhanced correction to the hydrogen energy levels. On the technical side, the calculation is similar to that in [8] if we use the effective field theory approach. To confirm our result, we also considered diagrams with triple photon exchange and extracted the logarithmically divergent part; both methods gave the same correction to the B 61 coefficient. The total result for the logarithmic LbL contribution, including the one from [8], reads Result (5) is over 2.5 times larger than that of the previous [partial] computation [8]. In standard notation (cf. [2]) it corresponds to the B 61 coefficient.
While considering the non-logarithmic part of the α 2 (Zα) 6 m correction, i.e., the coefficient B 60 and higherorder terms, one has to distinguish three groups of diagrams and treat them differently. One group originates from the 'pure' self-energy (SE) diagrams, i.e., the diagrams without any closed electron loops (see Fig. 2a). The remaining groups, on the other hand, include the closed electron loops. The second group contains the loops in the so-called free-loop approximation, i.e., all appearing closed electron loops are due to the vacuum polarization (see Fig. 2b). The last group contains virtual LbL scattering subdiagrams (see, e.g., Fig. 1a). (In high-Z atomic physics, those diagrams are referred to as the vacuum polarization in the presence of the Coulomb field of a nucleus.) b a The most accurately computed results exist for the free-loop approximation diagrams studied in [24]. The result reads G free 60 (Z = 1) = −15.0(4) , G free 60 (Z = 2) = −13.9(1) .
The contributions beyond the free-loop approximation at order α 2 (Zα) 6 m belong to two groups. One is due to radiative corrections to the Wichmann-Kroll contribution. (The Wichmann-Kroll contribution by itself is of order α(Zα) 6 m.) We estimate it as B rWK 60 (ns) = 0.13 ± 0.13 .
The estimation is based on a similarity of the behavior of a radiative correction to the Wichmann-Kroll potential and the Källen-Sabry potential in the so-called t channel.
The other group arises due to Coulomb corrections to the LbL contribution of order α 2 (Zα) 5 m. We have already considered their logarithmic part above in (5). We estimate the non-logarithmic part as The B 60 term beyond the free-loop approximation was previously estimated in [24]. However, it was based on incorrect assumptions about the logarithmic contributions for the diagrams beyond the free-loop approximation, and thus we do not take into consideration those estimates. The quantitatively largest contribution in our consideration of the B 60 term beyond the free-loop approximation comes as a 'tail' of the logarithmic B 61 term. The summary of the individual contributions to G 60 (1s) is given in Table III.
The estimate above is obtained by a suggestion that a natural magnitude of the constant accompanying a logarithm is π, which is inspired by the value of the imaginary part of the logarithm of a negative real number. In a term with several logarithms, we substitute each of them by π,   [24]. The pure SE value as well as the contribution beyond the free-loop approximation is a result of this letter.
which produces a combinatoric factor. Often the terms beyond the leading logarithmic term are estimated by 50% of its value. Using π and combinatoric factors, we estimate the subleading terms in the case of the leading single-logarithm as below 50%, but in the case of leading double-or triple-logarithmic term as above 50%. We think that it is more realistic than a naive 50%-estimate for all the separate cases.

IV. PURE SELF-ENERGY TWO-LOOP
The situation concerning the pure SE part of B 60 is more complicated than that of the closed-electron-loop contributions. A partial calculation exists, and it is accompanied by a plausible estimate of the unknown part of the contribution [9], B pure SE 60 = −61.6(9.2) .
The large magnitude of the B pure SE 60 coefficient is due to an enhancement of the low-momentum contribution, while the uncertainty comes from the unknown highmomentum one. Suggesting that the missing highmomentum contribution is not enhanced, we arrive at the result of ±π 3 B 63 for the missing contribution that coincides with the uncertainty in (9). Consequently, our estimation of the magnitude of the unknown terms in (8) is consistent with that in [9].
There exist essentially three approaches to calculation of the higher-order two-loop contributions. One suggests the use of an Zα expansion, in which case the accuracy is limited by (9) [9]. It is also necessary to know the higher-order logarithmic terms, such as [17] The size of the logarithmic contribution is smaller than of the uncertainty above, but not negligible. The second approach uses exact in Zα numerical calculations at Z = 1, 2. In the case of two-loop contributions, that approach has been successfully applied for the contributions with closed electron loops in the free-loop approximation [24] (see above), but its application to pure SE diagrams has proved challenging. Only the results at medium Z such as Z = 10, 12, 15, 17, 20, 25, 30 [10] are available. Those still can be extrapolated to Z = 0, 1, 2.
(The third approach includes a fit using the result of (10) as a data point at Z = 0 for G SE 60 (Z) from (4).) Due to the low accuracy such an extrapolation is possible for F (2) (Zα), only because a number of the leading coefficients, such as (B 40 , B 50 , B 63 , B 62 , and B 61 ) is known (see [2,15] and references therein). In the meantime the 'data area' (Z = 10,12,15,17,20,25,30) is relatively far from the 'target area' (Z = 0, 1, 2) and it contains relatively few data points. The logarithmic terms go through a bigger change on their way from the data area to the target area than within the data area. Accordingly, from the point of view of a phenomenological fit, we have to consider nearly coinciding fits with different extrapolation expectations. Since we need not only to fit the data but eventually also to extrapolate, we have to maintain the correct shape of the fit function (see (2)).
To deal with logarithmic terms at orders α 2 (Zα) 7 m and α 2 (Zα) 8 m, a calculation of some and an estimation of others is necessary. We present the summary of the leading logarithmic terms at each order in Zα in Table IV. The coefficients B 84 and B 93 are calculated in this letter using techniques developed in [17,26,27]. To estimate the subleading terms we use several approaches. As for an 'estimation' we understand a constraint with a relatively large uncertainty such as in (9), which allows us to obtain more than one estimate for each subleading coefficient. Ultimately, we choose the most conservative constraint for the coefficients. The summary of our estimations is given in Table V. The coefficients can be used for both for a low-Zα expansion and fits. We have to explain how we used the constraints for the fits. We consider the constraints as additional data and include their deviation from the related central values, measured in the units of their uncertainties, into the final χ 2 we have to minimize. That is, e.g., similar to a treatment applied in [2] for the various not very precise theoretical corrections. Such a fitting procedure allows one to easily combine theoretical constraints with the existing 'true' data.  Often, for the higher-order terms, it is possible to estimate the magnitude plausibly, but not the sign of a coefficient and therefore frequently the central values of estimations are zero. Once B 71,70 are estimated, we can find the result of the low-Zα expansion of G SE 60 (Z = 1, 2) (see Table VI). TABLE VI: Higher-order two-loop pure self-energy contribution to the 1s Lamb shift in hydrogen and the helium ion.
The combined fit includes the numerical data from [10] and the value of B60 from [9], and less accurate numerical results from [25].
To compare those low-Zα results with the numerical data, we need to fit them. We found that by setting B 72 = 0, the fit results for B 60 are shifted by 8 − 10 from the value B 60 of (9). That would put the fit and that value of B 60 in disagreement and would not allow a combined fit. All the previously used fits have ignored the double-logarithmic B 72 term. Consequently, the fits found in the literature use an unrealistic shape with no estimation of systematic effects (see, e.g., [10]). Consequently, a comparison with the previously performed fits is meaningless. We have performed a fitting of the numerical data [10] ourselves, using realistic approximation functions. We present the result in Table VI including the results of the combined fit, i.e., a fit which includes the low-Zα constraints and numerical data from [10]. We consider a difference between the low-Z value and the fit over the numerical data, which is somewhat below 2 σ as a fair agreement which validates the use of a combined fit.
The summary for the calculation of the two-loop contributions in the external field approximation is given in Table I.

V. NEXT-TO-LEADING THREE-LOOP CONTRIBUTIONS
The three-loop theory is more complicated and less advanced than the two-loop one. Only its leading term to the Lamb shift (order α 3 (Zα) 4 m, coefficient C 40 ) is known [28][29][30]. The next-to-leading one (order α 3 (Zα) 5 m, C 50 ) has been calculated only partially [31] and boldly estimated in [2]. After improvement of the accuracy of B 60 above, C 50 [2] becomes the largest source of the QED uncertainty for the 1s Lamb shift in hydrogen.
The α 3 (Zα) 5 m contribution can be represented as a set of two-photon exchange diagrams (see Fig. 3). The related expression may be written in terms of the integral over the loop momentum q (cf. [26]) where T (q 2 ) is a radiative correction to the skeletondiagram integral, which is related to a virtual forward Compton scattering amplitude. The calculation of the radiative factor T (q 2 ) is very complicated. Here, we calculate its asymptotics at high and low q and estimate the total integral, by integrating those asymptotic expressions. As mentioned previously, a part of the contributions, i.e., the diagrams with closed electron loops in free-loop approximation except for graphs with the two-loop pure self-energy with one electron vacuum-polarization insertion, have already been considered in [31]. Here we estimate the unknown diagrams by integrating the asymptotics of the related integrand. The complete three-loop result is C total 50 (ns) = −3.3(10.5) .
To verify our method, we have also found the contributions of order α(Zα) 5 m and α 2 (Zα) 5 m. Our estimation is in a perfect agreement with known results [6,32]. In each case of interest (one loop, two loops, three loops) the asymptotics of T (q 2 ) are of the same sign for high and low q. That is an important requirement for a reliable estimation of the integral through the asymptotics of the integrand.
The present situation with the three-loop contributions is summarized in Table II. The C 50 uncertainty is reduced by a factor of 3. This makes the C 50 uncertainty comparable with the CODATA's C 63 one in [2]. Fortunately, the latter was eliminated in [17], where it was found that

VI. SUMMARY AND CONCLUSIONS
The summary on the theoretical accuracy of the 1s Lamb shift calculation for light hydrogen-like atoms with Z = 1, 2 is presented in Table VII. The uncertainty from the external-filed contributions, considered in this paper, is due to α 8 m terms and consists of two sources, one is two-loop's G 60 and the other is three-loop's C 50 . A comparison with the existing calculations of other authors is given in the introduction, in Tables I and II in terms of the related coefficients and absolute values of the contributions for hydrogen. As one can see from there both two-loop and three-loop uncertainties are reduced approximately by factor of three.  RR16 stands for the α(Zα) 6 m 2 /M radiative-recoil contribution, which is known only in the leading logarithmic approximation (see (14)).
The uncertainty in Table VII comes from an estimation of subleading terms. For its estimation we use here the approach with π's and combinatoric coefficients, as explained above, and the uncertainty is somewhat above 50% (cf. [2]). The key uncertainty sources in [2] have also included pure recoil corrections, but their uncertainty (of about 0.7 kHz for H) has recently been eliminated [35] by a direct calculation of the recoil corrections for Z = 1, 2. Concluding, we have revisited the theory of the α 8 m contributions to the Lamb shift of the 1s state in hydrogen and deuterium atoms and helium ions. We completed the calculation of the logarithmic terms, considered a controversy in the non-logarithmic two-loop contribution and improved its accuracy by approximately a factor of three, and obtained a complete approximate result for the three-loop terms, which is more reliable and three times more accurate than a previous bold estimation.
The most accurate experimental results are available for the 1s − 2s transition in hydrogen and deuterium [14,36]. Experimental efforts to measure the 1s−2s transition in the helium ion are underway [12,13]. Since the weight of the individual contributions to the theoretical uncertainty budget varies substantially (see Table VII), combining the hydrogen and helium-ion experimental results would be beneficial not only for hydro-gen and helium-ion spectroscopy but also for various applications including precision tests of bound-state QED, determination of the Rydberg constant, and constraints on new light neutral particles.
A complete and detailed derivation, covering the technical side of the computations of our new results presented in this letter, is under preparation and will be published elsewhere.