Two-loop corrections to Lamb shift and hyperfine splitting in hydrogen via multi-loop methods

We revisit the contributions of order $\alpha^2(Z\alpha)^5m$ and $\alpha^2(Z\alpha)E_F$, respectively, to the Lamb shift and to the hyperfine splitting from mixed self-energy-vacuum-polarization diagrams, involving fermionic loop. We use modern multi-loop calculation techniques based on IBP reduction and differential equations. We construct the $\epsilon$-regular basis [LeeOnishchenko2019] and explicitly demonstrate that it is compatible with the renormalization. We obtain analytic results in terms of one-fold integral involving elliptic function and dilogarithm. As a by-product, we obtain the analogous contribution for the limiting cases of heavy and light fermionic loop.


Introduction
Hydrogen atom is the simplest atomic system.Traditionally it served as a touchstone for testing the bound-state quantum electrodynamics (QED).At present, the precision of both theory and experiment for the electronic hydrogen has been increased to such an extent that comparison of calculated and measured transition energies can be used for the most accurate determination of the Rydberg constant once the contribution of proton (which is the hydrogen nucleus) structure is established.The latter gets much more pronounced in the (electronic) hydrogen cousin -the muonic hydrogen.In particular, the Lamb shift in the muonic hydrogen provides the most precise value of the proton charge radius. 1he calculations of various contributions to the Lamb shift and the hypefine splitting have a long history starting from Refs.[1][2][3][4], see also review [5] and references therein.For higher order corrections these calculations often produced only numerical results, which might have insufficient accuracy for future comparison of theory and experiment.Also, despite these efforts, the present theoretical calculations are well behind the experimental measurements for some physical observables.E.g., the 1S − 2S transition frequency measurements have reached accuracy of a few tens Hz, while the corresponding uncertainty of the available theoretical predictions is only a few tens of kHz.Therefore, new ways of Lamb shift calculations are very much welcome.Recent progress in multi-loop calculations provides an opportunity to apply the developed methods to such calculations.
In the present paper we apply the multi-loop methods to the calculation of the contribution to the Lamb shift (LS) and to the hyperfine splitting (HFS) of the diagram, depicted in Fig. 1.The double line denotes the electron propagator in the electromagnetic field of the nucleus A µ = Z|e|/r, µ × r/r 3 .Here Z|e| and µ is the nucleus charge and magnetic moment, respectively.For LS calculations we neglect the nucleus magnetic field, while for the HFS calculations we consider linear in µ contributions.Note that the magnetic moment of the nucleus is the vector operator µ = µS/S.However, since we consider only linear in µ contributions, we can formally treat it as a numeric vector up to the point where we average the result over the spin wave functions.
The leading in Zα contribution of the diagram in Fig. 1 to LS and HFS is of the order O α 2 (Zα) 4 m e and O α 2 E F , respectively.This contribution is completely determined by the ∝ n f terms 2 in the slope of the Dirac form factor and the value of the Pauli form factor of the electron at zero momentum transfer at two loops calculated long ago in Refs.[6,7].Therefore, we do not consider these contributions.Here e and m e is the electron charge and mass, respectively, is the Fermi energy, 3 S is the nucleus spin, α = e 2 /(4π) ≈ 1/137.036 is the fine structure constant.
The next order contribution, O α 2 (Zα) 5 m e and O α 2 ZαE F , respectively, to LS and HFS comes from the diagrams, depicted in Figs. 2 and 3, which we will refer to as light-by-light (LbL) contribution and "free loop" (FL) contribution, respectively.
The light-by-light contributions to LS and HFS were calculated in Refs.[8][9][10][11] and [12,13], respectively.The result for the Lamb shift was presented in terms of a fourfold integral, while that for the hyperfine splitting was presented in terms of a three-fold integral.These integrals were then treated numerically.The free loop contribution to Lamb shift was obtained in terms of two-fold integrals in Refs.[8,14], while the corresponding contribution to the hyperfine splitting was obtained in Refs.[13,15] in terms of one-fold integral involving elliptic function and logarithm (or rather, arctangent).
In a sense, the result of the present paper is the representation of all four contributions (LbL and FL to LS and HFS) in the form similar to that of Ref. [15].We use modern multi-loop methods, namely, the IBP reduction [16] and the differential equations method [17,18].In order to apply the latter, we cosider the diagrams where the mass of the fermion in the loop is different from that of the electron line.As a by-product, we also obtain the contribution of the muon loop for the ordinary hydrogen and that of electron loop for the muonic hydrogen. 2Here we follow the standard convention by introducing a formal parameter n f counting the number of electron loops.

Calculation
We consider the diagrams depicted in Figs. 2 and 3.The contribution of each set is gauge invariant and can be calculated independently.Note that for the HFS calculations one should replace one of the two Coulomb exchanges in those diagrams with the magnetic exchange corresponding to the contribution of the nuclear magnetic moment.
The external electron legs on the diagrams denote the bound electron wave function, however, with the precision that we pursue here, the bound-state effect is properly taken into account by the factor |ψ n,l (0)| 2 = δ l,0 (meZα) 3 πn 3 , where δ is a Kronecker symbol, n and l are a principal and azimuthal quantum number, correspondingly.Indeed, the characteristic loop momenta in the discussed diagrams are k i ∼ m e , while that of the wave function |p| ∼ m e Zα m e .Therefore, we calculate the diagrams in Figs. 2 and 3 for the external electron legs corresponding to free electron with momentum p = m e τ , where τ = (1, 0) is the time ort.
The LbL contribution is both UV and IR finite, while the FL contribution is UV and IR divergent.In order to obtain the finite result we add one-loop counter terms, depicted in Fig. 3 (second line).The two-loop counter terms, as well as the IR subtraction are expressed via scaleless integrals which are equal to zero in dimensional regularization.In particular, denoting the momentum of the Coulomb line as k 1 = (0, k 1 ), we see that the two-loop counter terms contain only one scaleless denominator (p − k 1 ) 2 − m 2 e = k 2 1 .The Lamb shift energy correction can be represented as where u Ôu corresponds to the diagrams in Figs. 2 and 3 with both incoming and outgoing momenta equal to p = m e τ .
In the case of hyperfine splitting we should replace one Coulomb exchange with the magnetic exchange.Since we use the dimensional regularization, we should avoid using γ 5 and ε ijk because those objects are poorly generalized to the generic spacetime dimension d = 4 − 2 .Therefore we rewrite all formulas involving vector product (a × b) k = ε ijk a i b j in terms of antisymmetric tensors a i b j − a j b i .We obtain where u Ôi,j HFS u corresponds to the sum of the diagrams in Figs. 2 and 3 in which one of the two Coulomb exchanges Z|e|γ 0 /k 2 1 is replaced by the "magnetic exchange" Differential system and boundary conditions.
In order to apply the differential equations method, we decouple the mass of the bound electron and the mass of the particle in the loop.We put the latter to 1 while keeping the mass of the bound electron as a free parameter m.
For the light-by-light contribution we consider the integral family where The δ-function in Eq.( 2.3) corresponds to zero energy transfer on the nucleus.Note that n 8 is not positive.
For the free-loop contribution we consider the family where where n 7 and n 8 are not positive.
Making the IBP reduction [16,19] with LiteRed [20], we reveal 14 master integrals for light-by-light contribution and 8 master integrals for free-loop contribution.Four master integrals are common for the two bases.Thus in the merged basis we have 14 + 8 − 4 = 18 master integrals presented in Fig. 4. The first eight graphs correspond to FL basis, while the graphs 3 − 6 and 9 − 18 correspond to LbL basis 4 .
Using the IBP reduction, we construct the differential equations [17,18] for the master integrals with respect to m 2 : where j = (j 1 , ..., j 18 ) is a column of master integrals.We use Libra [21] to manipulate the differential system (2.7).We find it convenient to work with the master integrals in 2 − 2 dimensions and later express the master integrals in 4 − 2 dimensions via lowering dimensional recurrence relation [22].Note that the differential system (2.7) can not be reduced to -form due to integrals j 3−6 which appear both in FL and LbL contributions.Those four master integrals can be expressed via hypergeometric functions.In particular, Due to this reason the corresponding block of matrix M in Eq. (2.7) can not be reduced in -form.All other blocks can be reduced to -form.With some trial and error method we have succeeded to obtain the (A + B)-form of the differential equation, where M 1 (m 2 , ) = A(m 2 ) + B(m 2 ) and the matrix A has nonzero entries only in the columns 3-6.We fix the boundary conditions at the point m 2 = 0.The general solution has the form where C 0 is a column of the boundary constants and Using Libra, we relate the boundary constants C 0 to specific asymptotic coefficients of original master integrals j 1−18 at m 2 → 0. We compute these coefficients by using expansionby-regions method [23] and direct integration of Feynman parametrization.We find that the boundary constants C 0 are expressed in terms of the following 5 nonzero constants: 2 , 1, 2 + 3 2 ; 1 − , + 2; 1 , (2.12) where [j k ] m ν denotes the coefficient in front of m ν in the small-mass asymptotics of j .Nevertheless, using the PSLQ and our experience, we have been able to establish that the -expansion of j (2−2 ) 15 can be written in terms of Goncharov's polylogarithms at fourth root of unity.The two first terms of -expansion read where 2k+1) 2 = 0.915965594... is the Catalan constant.For our present goal we will need only the leading term of this expansion.
-regular basis Since the differential system (2.7) for the master integrals can not be reduced to -form, we should rely on the Frobenius method for the calculation of the evolution operator U (m2 , 0) in Eq. (2.11).The dependence of this operator on constitutes substantial technical difficulties and blocks the way to high-precision numerical calculation suitable for the application of PSLQ algorithm, [24], for the recognition of the master integrals in the point m = 1.Thus we choose to switch to the -regular basis [25].After finding this basis, we simply put = 0.
Note that the counter-term diagrams in the second line of Fig. 3 can also be expressed in terms of the three-loop master integrals in Fig. 4 although they have only two loops.To this end we multiply the corresponding integrals by 1 = . Then the contribution of counter-terms is expressed via the master integrals with unit mass tadpole loop, namely, via j 1 , j 2 .Luckily, Γ[1 − d/2] in the denominators cancels with the same Γ in δZ which stands for the cross in the couterterm diagrams.Therefore, the finite sum of all diagrams in Fig. 3 is expressed via the integrals of the family (2.5) with rational coefficients 5 .Then we are in position to state the existence of -regular basis [25].
Let us describe how we construct this -regular basis.
3. Then we use the following rule of thumb: if for, say, J i function both the boundary constant C i and the right-hand side of the equation are zero at = 0 is zero, then the function itself is also zero at = 0. Therefore, we redefine J i → J i and modify respectively the boundary constant C i → C i and the differential system.The latter modification is reduced to the multiplication/division by the i-th column/row of the matrix A + B.
4. We also use substitutions of the form J i → J i + k<i c k J k , where the coefficients c k are rational numbers chosen so as to nullify as many entries on the i-th row of the matrix M (m 2 , 0) as possible.
This approach works perfectly for all rows except for the rows 3 − 6 corresponding to the equations irreducible to -form.For those 4 master integrals we use the presumable finiteness of the sum of diagrams in Figs.2,3 as a guiding line for finding the relations between j 3−6 near d = 2.We find that Thus we choose J(m 2 , ) = Q(m 2 , )/ 2 as an element of -regular basis.Indeed, we see that after this choice both the boundary constants and the matrix M (m 2 , ) have finite limit → 0.
The result of our approach is the differential system in which we can simply put = 0.The higher orders in which we miss with putting to zero are guaranteed not to appear in our final results, which is the rationale behind the notion of -regular basis.The resulting system has the form where the matrix M is strictly lower triangular except for the diagonal 4 × 4 block with indices 3-6.It is essential that M does not depend on .The singular points of the differential system are m 2 = 0, 1, ∞.Again, we write the solution in the form The almost lower triangular structure of the matrix M allows us to write the general solution in terms of polylogarithms and/or one-fold integrals of j 3−6 multiplied by polylogarithms.However we choose here to construct the Frobenius expansions near each of the three singular points of the differential system (2.16), m 2 = 0, 1, ∞.We match the obtained expansions pairwise in the points which belong to the intersection of convergence regions of the corresponding two expansions.For example the expansions near m 2 = 0 and m 2 = 1 are connected via the relation (2.17) Then the boundary constants at m 2 = 1 are expressed as We calculate 1000 terms of series of U (1/2, 1) and U (1/2, 0) and compute more than 300 digits for C 1 .In order to use PSLQ recognition, we need to have a basis of appropriate transcendental numbers and we extract all but one required nontrivial constants from Ref. [15].Therein the result for the FL contribution to HFS was expressed in terms of the weighted sum of the integrals (cf.Eq. ( 11) of Ref. [15]) where K(x) and E(x) are complete elliptic integrals.Moreover, using PSLQ it is easy to establish a relation 6c 2 − 5c 3 + 2c 4 = 0 overlooked in Ref. [15].Using some guess work, we find the following basis sufficient for our purpose B = {1, π, ln 2, π 2 , ζ 3 , e 0 , 1/e 0 , e 1 , e 2 } , e 0 = 3 4π where Li 2 (x) is a dilogarithm.The constants c 1−4 in Eq. (2.18) are expressed as linear combinations with rational coefficients of {1, e 0 , 1/e 0 , e 1 }: (2.20) The benefit of using e 0 , e −1 0 , and especially e 1 instead of c 1−4 is that their form allowed us to guess the form of the last nonstandard constant e 2 , which was deduced from e 1 by noting that arctan (x) = Im Li 1 (ix) and then replacing Li 1 → Li 2 .To obtain the contributions to the Lamb shift and the hyperfine splitting we need a few coefficients of the expansion of J reg (m 2 ) near the point m 2 = 1.This expansion has the form where U (M |m 2 , 1) has the form of generalized power series in 1 − m 2 .In particular, U (m 2 , 1) contains ln(1 − m 2 ).However, we have checked that these logarithms disappear when U (m 2 , 1) is multiplied by C 1 so that the specific solution has no branching at the point m 2 = 1 as it should be.The power series (2.21) converge when m 2 ∈ (0, 2).To obtain the results for J reg (m 2 ) at m 2 > 2 we pass to the variable z = m 2 −1 m 2 and again match the power series near z = 0 (m 2 = 1) and z = 1 (m 2 = ∞) at z = 1/2.In this way we obtain the high-precision numerical result for the column of boundary constants C ∞ which define the coefficients in the asymptotic expansion of J reg (m 2 ) at m 2 → ∞.In order to define the analytic form of C ∞ we again use PSLQ recognition with the following basis: The nontrivial constants i 0 and i −1 0 were conjectured by examining the large-mass asymptotics of j 3 from Eq. (2.8).

Results
The discussed corrections to the Lamb shift and hyperfine splitting have the following form where i = F L, LbL for different contributions and we have recovered the recoil factor accounting the main effect of the finite nuclear mass M .For the case of electron loop in electron atom we have the following result: One is tempted to obtain also the results for the contribution of the electron loop in the muonic hydrogen by making a substitution m e → m µ in Eq. (3.1) and using the large-m asymptotics of J reg .Then the results for functions B i and C i have the following form However, for muonic atom the characteristic momenta of bound muon (or the inverse size of its wave funcetion), m µ Zα is of the same order, or as the characteristic loop momenta m e even for Z = 1.Therefore, the applicability condition for the approximation used in the present paper is violated and Eq.(3.4) may be considered at most as an order of magnitude estimate.

Conclusion
In the present paper we revisit the contributions of order α 2 (Zα) 5 m to the Lamb shift and of the order α 2 (Zα)E F to the hyperfine splitting from mixed self-energy-vacuumpolarization diagrams depicted in Figs. 2 and 3. We construct the -regular basis [25] and explicitly demonstrate that its elements taken at = 0 are sufficient to express the renormalized results.The results have the form of one-fold integral involving elliptic function and dilogarithm and agree with previously known numerical results.We also obtain the contribution of the same set of diagrams with different mass of the fermion in the loop and in the fermion line, which allows us to determined the corresponding contribution of the muonic loop in the conventional hydrogen as well as the estimate of the electron loop contribution in the muonic hydrogen.

Figure 1 .
Figure 1.Two-loop self-energy diagram in Furry representation.Double lines denote electron propagators in the external electromagnetic field A µ = Z|e|/r, µ × r/r 3 .

Figure 2 .
Figure 2. Gauge invariant set of diagrams which corresponds to the light-by-light contribution to the Lamb shift.Numbers correspond to the enumeration of denominators in Eq. (2.4).

Figure 3 .
Figure 3. Gauge invariant set of diagrams which corresponds to the free-loop contribution to the Lamb shift.The second line of diagrams is a one-loop counter-terms.Numbers correspond to the enumeration of denominators in Eq. (2.6).

Figure 4 .
Figure 4. Basis of master integrals.Thick and thin black lines correspond to the denominators with mass 1 and m, respectively.Thick gray line corresponds to δ(−D 9 ), dashed lines correspond to massless denominators.
last constant are trivially expressed in terms of alternating multiple zeta values.It is not obvious which transcendental numbers might enter the expansion of the last constant j (2−2 ) 15 [8][9][10][11][13][14][15]f Eq.(2.19)for the constants e 0 , e 1 , e 2 , the above expressions provide a one-fold integral representation for the corresponding contributions to Lamb shift and hyperfine splitting.The numerical values for B F L , B LbL , C F L , C LbL agree with the results of Refs.[8][9][10][11][13][14][15].For the case of muon loop in electron atom we have the following result: