Dimensional Regularization and Two–Loop Vacuum Polarization Operator: Master Integrals, Analytic Results and Energy Shifts

We present a complete reevaluation of the irreducible two-loop vacuum-polarization correction to the photon propagator in quantum electrodynamics, i.e. with an electron-positron pair in the fermion propagators. The integration is carried out by reducing the integrations to a limited set of master integrals, which are calculated using integration-by-parts identities. Dimensional regularization is used in D = 4 − 2 ǫ dimensions, and on-mass shell renormalization is employed. The one-loop eﬀect is given to order ǫ , to be combined with the 1 /ǫ divergence of the two-loop amplitude. Master integrals are given. Final evaluations of two-loop energy shifts for 1 S , 2 S , and 2 P states are done analytically, and results are presented, with an emphasis on muonic hydrogen. For relativistic Dirac–Coulomb reference states, higher-order coeﬃcients are obtained for the Zα -expansion. We compare the results obtained to the existing literature.


I. INTRODUCTION
The two-loop vacuum-polarization correction to bound-state energy levels is an important contribution to the Lamb shift in muonic bound systems.Specifically, the irreducible two-loop effect lowers the 2S level energy in muonic hydrogen by roughly 1.25 meV in comparison to the 2P level [1,2].This energy shift needs to be compared to the proton-size discrepancy commonly referred to as the proton radius puzzle (0.88 fm versus 0.84 fm, for a recent brief summary see Ref. [3]).The proton radius puzzle corresponds to a 0.3 meV shift of the 2S states, when expressed in energy units [1,2].The two-loop effect is thus phenomenologically extremely relevant; any conceivable inaccuracy in its evaluation could lead to at least a partial theoretical explanation of the proton size discrepancy.Vacuum-polarization effects are drastically enhanced in muonic as compared to electronic bound systems [4][5][6][7].Here, we present a reevaluation of the irreducible two-loop diagrams, on the basis of modern integration-by-parts techniques, employing dimensional regularization.We note that dimensional regularization was not yet sufficiently developed in 1973 (despite Refs.[8,9] which appeared in the literature at the time) to make an evaluation of the two-loop effect using dimensional regularization techniques feasible [10].For vacuum-polarization effects in particular, the modern techniques lead to drastic simplifications of the calculations.We find both the imaginary as well as the real part of the two-loop vacuum tensor in closed analytic form.
We also derive analytic expressions for the expectation value of the two-loop potential, evaluated on (nonrelativistic) Schrödinger-Coulomb eigenstates, generalizing the treatment originally outlined in Ref. [4] to twoloop order.In general, for predominantly nonrelativistic bound systems, energy shifts can be represented in terms of a semi-analytic expansion (in powers of Zα and ln[(Zα) −2 ], where α is the fine-structure constant and Z is the nuclear charge number).As a byproduct, we de-Our aim is to evaluate the one-loop vacuumpolarization effect to order ǫ, i.e., in a form suitable for the later two-loop calculation.For simplicity we set the lepton mass (electron mass) in the vacuumpolarization loops to be m = 1.We use for a metric tensor g µν = diag(1, −1, −1, −1) (see, for example, Refs.[11,12]), where p µ = (p 0 , p), p µ = (p 0 , − p), Here, m = 1 is the electron mass.Due to current conservation, the vacuum polarization function Π µν has this tensorial structure Π µν (q 2 ) = (q µ q ν − q 2 g µν ) Π(q 2 ) . ( We follow the convention of Ref. [13] so that dΠ(q 2 ) dq 2 > 0 (5) at the one-loop and two-loop level.For completeness, we note that the opposite conventions are used elsewhere in the literature, e.g., in Refs.[14,15].The conventions employed here are consistent with those employed in Ref. [13].Summation over the Lorentz indices in D dimensions leads to the result Π µ µ (q 2 ) = (−q 2 g µ µ + q µ q µ ) Π(q 2 ) = q 2 (1 − D) Π(q 2 ) .
The threshold for pair production of Π(q 2 = s) is q 2 = s = 4.The argument of the polarization tensor is which is the four-momentum square of the photon entering the vacuum-polarization loop.

B. Calculation
For the renormalization of the two-loop diagram with a self-mass insertion we need the one-loop vacuum polarization function Π (1) , and the function Π (1a) from the same diagram with iterated electron propagators on one side of the loop.The unrenormalized tensors read as follows, Π (1)  µν (q 2 ) = −ie 2 The needed quantities are expressed in terms of two oneloop master integrals, The prefactor in these results simplifies to unity provided we set One observes that the identification ( 12) is not completely canonical in dimensional regularization.Namely [see, e.g., Eq. ( 10.173) of Ref. [15]], in the MS scheme, one normally has the relation e 2 = (4π) 1−ǫ α µ 2ǫ e γE ǫ , where γ E = 0.57721 56649 . . . is the Euler-Mascheroni constant, and µ is the renormalization scale induced by dimensional regularization.This is equivalent to Eq. ( 12) up to terms linear in ǫ.The prefactor 1/C(ǫ) helps to simplify results after integration.The one-loop integrals are as follows, where The above definition of [d D k] contains the factor µ −2ǫ .It simplifies with the factor 1/µ −2ǫ of N (ǫ) so that the master integrals are independent of µ.The normalization factor is .
(16) The latter form exhibits the dimension D = 4 − 2ǫ.With the normalization factor N (ǫ), all master integrals M occurring in our calculations are dimensionless and not only uncluttered from logarithms of µ, but also, from logarithms of 4π.
One can show that the one-loop master integrals satisfy the following differential equations, where the differentiation is with respect to q 2 .One finds that M 2 is a momentum-independent vacuum integral, For M 12 (q 2 ), one writes the result as a function of the variable The result is where we include the terms of order ǫ and ǫ 2 which will be useful in the context of the two-loop calculation.The explicit form of N 12 (q 2 ) is given in Appendix A.

C. Renormalization
Substituting Eqs. ( 18) and ( 20) into (13), we can find the expansions of P (1) (q 2 ) and P (1a) (q 2 ).But for the renormalization we need to subtract the values at q 2 = 0, In our formalism, the scaled functions P (1) and P (2) do not have powers of µ, but the scalar vacuum polarization functions Π (1) , Π (2) depend on µ.Yet, after renormalization the dependence disappears in the limit ǫ → 0. The renormalized expression for the one-loop function P R (q 2 ) is as follows, P R (q 2 ) = P (1) (q 2 ) − P (1) (0) The coefficient of the term linear in ǫ is For the explicit form of the term of order ǫ 2 , we refer to Appendix A. The renormalized value for the function P (1a) R (q 2 ) with iterated electron propagators is as follows, The term linear in ǫ is given as follows, The quadratic term in ǫ, denoted as (v), is given in Appendix A. Actually, only the terms proportional to ǫ are necessary for the renormalization of the two-loop vacuum polarization, but higher-order terms are useful for higher-loop calculations, and are thus relegated to Appendix A.
In the limit ǫ → 0 + , one recovers the known result for the one-loop vacuum polarization, The dispersion relation is The imaginary part is positive infinitesimally above the cut, and negative infinitesimally below the cut.One has the simple representation Im Π (1) On the cut, the v variable runs from v = 0 to v = 1.

B. Reduction to Master Integrals
After the calculation of traces, performed with computer algebra (using FORM, see Refs.[16,17]), and a Wick rotation, the scalar contribution Π (1) (q 2 ) and Π (2) (q 2 ) of the two diagrams are reduced to a combination of 24 and 41 Feynman integrals, respectively.In order to further reduce the expressions, we use the integration by parts (IBP) identities [18,19].A system of IBP identities is built and solved with the program SYS, which is based on an algorithm described in Refs.[20,21].These identities allow to reduce Π (1) (q 2 ) and Π (2) (q 2 ) to a linear combination of five irreducible master integrals.

C. Differential Equations and Master Integrals
A system of differential equations in q 2 is being created; we need to compute the 5 master integrals.Three of them are reducible into products of one-loop integrals, as outlined in Eq.(34.We obtain a system of differential equations in q 2 satisfied by the irreducible master integrals [22][23][24].The first of these equations is of second order, while the remaining one is a first-order differential equation, where the differentiation is with respect to q 2 .The initial conditions at q 2 = 0 for these diagrams are simple vacuum integrals.

D. Results for the Master Integrals
Due to the appearance of the factor D − 4 in the denominators of Eqs. ( 36)-( 38) (for the coefficients, see Appendix B), the master integrals have to be evaluated up to order O(ǫ).It is not sufficient to keep only the finite part.The results are as follows.For M 235 , one finds the following relation, where we note the prefactor The order-ǫ coefficient N 235 (q 2 ) is found in Appendix C.
For M 235k , one finds the relation Again, N 235k (q 2 ) can be found in Appendix C. Similar master integrals have recently been considered in Refs.[25,26].

E. Renormalization
We now substitute the master integrals (34a), (34c), (34b), ( 40) and ( 41) into the reduction formulas given in Eqs.(36) and (38).We can thus find the series expansions of the vacuum-polarization functions Π (2;1) (q 2 ) and Π (2;2) (q 2 ).For the renormalization, we need the values at q 2 = 0.The behavior, for small q 2 , is the following: Terms containing 1/q 2 appear, due to the fact that the single diagrams are not gauge invariant.Taking the gauge invariant sum P (2;1) (q 2 ) + 2P (2;2) (q 2 ), the divergent terms cancel out, so that, taking the limit q 2 → 0 one gets In order to carry out the renormalization, we need the renormalization constants where J stands for the vertex renormalization (J = 1), the wave function renormalization (J = 2), and the mass renormalization (J = m), and L enumerates the loop order.We define the one-loop quantities 2 , (45) and obtain the results We note that Z (1) 2 because of the Ward-Takahashi identity.The identity m is a simple coincidence, which occurs only at one-loop order in dimensional regularization, but to all orders in ǫ.The constants F , S 0 and S 1 in Eq. ( 45) are defined with the prefactor (−1) so that the expansion of the the self-energy insertion into the fermion line about the mass shell acquires positive terms, for ✁ q ≈ m = 1: (47) Finally, the renormalized contributions of the diagrams are found as follows, P (2;1) R (q 2 ) = P (2;1) (q 2 ) − P (2;1) (0) − 2F P R (q 2 ) , (48) P (2;2) R (q 2 ) = P (2;2) (q 2 ) − P (2;2) (0) and the renormalized two-loop function P (2) R (q 2 ) is obtained as We recall the definition of v in Eq. ( 19).In the limit ǫ → 0 + , one obtains the two-loop vacuum-polarization function We find a compact representation (see also Ref. [27]) where (53) For small positive q 2 , one obtains the expression The representation for the imaginary part just above the cut is even more compact, Im Π G. Comparison with the Literature An essential ingredient of our calculations is the master integrals, which should be compared to results communicated in Refs.[28][29][30].In Refs.[28,29], the authors define a master integral J 011 which is equivalent to our M 12 up a factor In Eq.(1) of Ref. [28], the general term of the ǫ-expansion of J 011 is given in terms of log-sine integrals which can be written in terms of Nielsen polylogarithms; the analytic continuation is shown in Eqs.(2.9)-(2.23) of Ref. [29].
Changing the variable from t = q 2 to v and using the transformations for Nielsen polylogarithm S n,p (−1/z) → S n,p (−z), we found that the results in our Eqs.(20) and (A1) agree perfectly with the corresponding terms of Eq. (1) of Ref. [28] and the expansion given in Eqs.(2.10)-(2.14) of Ref. [29].
In Eqs.(4.9), (4.10), and (E.4) of Ref. [30], analytical results for master integrals equivalent to our M 235 and M 235k were presented.These master integrals are related to ours by the following relations: one finds for the ultraviolet-convergent integral (4.9) of Ref. [30], while, for the ultraviolet-convergent linear combination of integrals given in Eq. (4.10) of Ref. [30], one has and for the ultraviolet divergent integral (E.4) of Ref. [30], The ǫ−expansion obtained from those of M 235 and M 235k agree perfectly with those presented in Ref. [30].
Our results for the complete two-loop function given in Eq. ( 52) and for its imaginary part, as given in Eq. ( 55), agree with the calculation of Ref. [31].One notes that Eq. (49) of Ref. [31] contains the imaginary part, while Eq.(57) of Ref. [31] also contains the real part of the twoloop vacuum polarization.One should note, though, that some integrals are left unevaluated in Ref. [31].Furthermore, it needs to be pointed out that the expressions in Ref. [31] also contain the reducible diagram (see diagram 3 in Fig. 1) with two iterated one-loop vacuum polarizations on the same line.Our result for the imaginary part also agrees with the result given in Eq. (5-4.200) of Ref. [32].Our results are also in agreement with the paper [33], specifically, with the results given in Eqs. ( 71), (78), ( 79), (80), and (81) of Ref. [33].Furthermore, we can refer to related investigations on the vector part of the quantum chromodynamic vacuum-polarization tensor [34][35][36].

IV. CONTRIBUTIONS TO ENERGY SHIFTS A. Muonic Hydrogen
Armed with the compact expressions given in Eqs. ( 52) and (55) for the imaginary part of the two-loop vacuumpolarization function, we can evaluate energy shifts in hydrogenlike ions.In Chap. 10 of Ref. [15], the calculation of the vacuum-polarization energy shift is outlined in detail; specifically, we refer to Eqs. (10.244) and Eq. ( 10.245) of Ref. [15].(We note the different sign conventions for the sign of the imaginary part of the renormalized scalar function Π R , as compared to Ref. [15].)From the one-loop and two-loop scalar functions, one infers the one-loop and two-loop irreducible vacuum-polarization potentials R (q 2 + iǫ) . (62) R (q 2 + iǫ) . (63) Of particular phenomenological importance for the proton radius puzzle [3] is the contribution of the irreducible two-loop vacuum-polarization effect on the 2P -2S energy interval in muonic hydrogen.We use the following results, Here, the β parameter is where Z is the nuclear charge number, α is the finestructure constant, and µ is the reduced mass of the system.For muonic hydrogen, one has Z = 1.In units where the electron mass is m = 1 and the muon mass is m µ , one has µ = m p /(1 + m p ) for ordinary hydrogen and µ = m µ m p /(m µ + m p ) for muonic hydrogen.The contribution of the irreducible two-loop vacuum polarization diagrams to the 2P -2S splitting in muonic hydrogen is obtained as Quite surprisingly, the two-loop energy shift can be expressed analytically.We define and find the result where the expressions I j (a) with j = 1, 2 read as follows, Some of the terms could also be expressed in terms of the Legendre χ function, For muonic hydrogen (Z = 1), one obtains for the nonrelativistic energy shift due to the irreducible twoloop diagrams, using CODATA values for the relevant physical constants [40,41] (see also Ref. [27]), This is consistent with the literature (e.g., Ref. [2]).The result (71) confirms the numerical value of a contribution to the proton radius puzzle in muonic hydrogen whose numerical magnitude could have potentially contributed to an explanation the discrepancy.The confirmation is obtained based on modern quantum-field theoretical methods.The reducible diagram, according to Eqs. (D8) and (D9), leads to an energy shift of For selected individual low-lying atomic reference states, an analytic integration of the one-loop and two-loop energy shifts, for arbitrary β, still in the nonrelativistic approximation, is presented in Appendix D.

B. Relativistic Ordinary Hydrogen
Given our analytic approach, we may also also study the energy shifts to low-lying S and P states for ordinary (electronic) hydrogen, with a relativistic reference state.We focus on higher-order terms in the semi-analytic expansion in powers and logarithms of Zα, for relativistic Dirac-Coulomb reference states (see Chap. 8 of Ref. [15]).
The potential R due to the reducible diagram is R (q2 + iǫ)
(73) We study the energy shifts where |ψ is the relativistic reference-state wave function, and ψ stands for the Hermitian adjoint (not for the Dirac adjoint ψ = ψ † γ 0 ).The relevant semi-analytic expansion is as follows, For the nP 1/2 levels, the coefficients are as follows, B 60 For the nP 3/2 levels, the coefficients are as follows, B The fine-structure difference is consistent with the result communicated in Eq. (7.7) of Ref. [42].

V. CONCLUSIONS
We have completed a complete reevaluation of the twoloop vacuum-polarization tensor in dimensional regularization, on the basis on integration-by-parts identities.The two-loop vacuum-polarization tensor constitutes a numerically significant contribution to the Lamb shift of muonic hydrogen which influences the determination of the proton radius from muonic hydrogen spectroscopy [2].
In Sec.II, we have discussed the evaluation of the oneloop vacuum polarization insertion into the photon propagator, evaluated to order ǫ, and ǫ 2 , and thus, in a form suitable for inclusion into higher-order loop calculation where knowledge of the terms of higher orders in ǫ is indispensable.We note that we use somewhat nonstandard conventions for the MS-charge, as given in Eq. ( 12).
In Sec.III, the irreducible two-loop vacuum polarization insertion has been evaluated, by expressing it in terms of master integrals.The renormalization has been carried out, and final results have been presented for the real and imaginary parts, in Eqs. ( 52) and (55).A comparison to the existing literature is being performed as well (Sec.III G).
In Sec.IV and Appendix D, we demonstrate that, for arbitrary reduced mass of a two-body bound system, the two-loop vacuum-polarization corrections to the energy can be evaluated analytically (for nonrelativistic reference states) and expressed in terms of dilogarithmic, and trilogarithmic, functions.This applies both to the 2P difference (see Sec. IV A) as well as to individual hydrogenic levels (see Appendix D).Higher-order coefficients for the semi-analytic expansion of the two-loop vacuumpolarization energy shift could be evaluated with the help of the fully relativistic hydrogen wave function.Results have been presented in Sec.IV B.

ACKNOWLEDGMENTS
The authors acknowledge insightful conversations with Professor Ettore Remiddi.Support from the National Science Foundation (grant PHY-2110294) is gratefully acknowledged.S.L. also acknowledges support from Italian Ministry of University and Research (MUR) via the PRIN 2022 project n. 20225X52RA -MUS4GM2 funded by the European Union via the Next Generation EU package.We first give the terms quadratic in ǫ for Eq. ( 20), for the master integral M 12 , We now give the expressions for the terms of order ǫ 2 in the one-loop effect, as discussed in Eq. (22).For the quadratic term Q 2 (v) in Eq. ( 22), one obtains the result For the term Q (v) from the integral with two fermion propagators, one finds for the term quadratic in ǫ the expression The coefficients in Eq. ( 37) are as follows,

G
(2;1) The coefficients in Eq. ( 38) are as follows, Appendix C: Higher-Order Terms in the Master Integrals We have the following results for the higher-order terms from the master integrals, with reference to Eqs. ( 40) and (41).For N 235 (q 2 ), we obtain For N 235k (q 2 ), we obtain We choose the 2P -2S interval in order to demonstrate that energy shifts due to the one-loop and two-loop vacuum-polarization effects can be evaluated analytically.For muonic hydrogen, one has the result µ = m µ m p /(m µ + m p ), where m µ and m p are the muon and proton masses, respectively.
It turns out to be advantageous to express the result in terms of the variable w defined in Eq. ( 67), with the result We can get the expansions for small and large β, For muonic hydrogen, one has β = 0.737384 . . ., while, for ordinary hydrogen, one has β = 137.110 . . . .Hence, the second of the above equations is relevant for ordinary hydrogen, while none of the expansions can be used with good accuracy for muonic hydrogen.It should be possible, though, to generalize the results reported in this Appendix to reference states other than 2S and 2P if desired.Relativistic corrections to the vacuumpolarization shift have been analyzed in Ref. [43], with an emphasis on muonic hydrogen.

Two-Loop Reducible Diagram
Next, we consider the contribution of the reducible diagram 3 in Fig. 1, Im Π R (q 2 + iǫ) We use the conventions outlined Eq. ( 67).An analytic evaluation leads to the formula (D9) The expansion for small β reads as follows, For large β, one obtains the result (D11)

Two-Loop Irreducible Diagram
We continue with the contribution of the irreducible two-loop vacuum polarization diagrams, The result is a bit more complex and given in Eq. ( 68).The expansion for small β is In Sec.IV A, we had concentrated on the 2P -2S energy difference, for muonic hydrogen.It is instructive to separately give the nonrelativistic contributions to selected energy levels with principal quantum number n and orbital angular momentum quantum number ℓ (in the nonrelativistic approximation) from the one-loop diagram, the two-loop irreducible and the two-loop reducible diagrams, E (1) (nℓ) = α π (Zα) 2 µf (1) We consider the cases |nℓ = |1S , |2S , |2P .The variable w has been defined in Eq. (67).We define the variable u as follows, For the ground state, we obtain f

FIG. 1 .
FIG.1.The diagrams concern the electronic two-loop vacuum-polarization diagrams in muonic hydrogen (the negatively charged muon line is denoted by µ − ).These are are naturally divided into diagram (1), which is a proper twoloop diagram, diagram(2), which is a one-loop diagram with a self-energy insertion, and the loop-after-loop (reducible) diagram (3).
are given in Appendix B. The gauge-dependent part of diagram 1 is given as indicated in Appendix B. For diagram 2, one has the following reduction, Appendix A: Terms of Order ǫ 2

TABLE I .
Coefficients from the reducible diagram are indicated for S and P states.

TABLE II .
Coefficients from the diagrams are indicated for S and P states.