The two-photon exchange contribution to muonic hydrogen from chiral perturbation theory

We compute the spin-dependent and spin-independent structure functions of the forward virtual-photon Compton tensor of the proton at one loop using heavy baryon effective theory including the Delta particle. We compare with previous results when existing. Using these results we obtain the leading hadronic contributions, associated to the pion and Delta particles, to the Wilson coefficients of the lepton-proton four fermion operators in NRQED. The spin-independent coefficient yields a pure prediction for the two-photon exchange contribution to the muonic hydrogen Lamb shift, $\Delta E_{\rm TPE}(\pi\&\Delta)=34(13)$ $\mu$eV. We also compute the charge, $\langle r^n \rangle$, and Zemach, $\langle r^n \rangle_{(2)}$, moments for $n \geq 3$. Finally, we discuss the spin-dependent case, for which we compute the difference between the four-fermion Wilson coefficients relevant for hydrogen and muonic hydrogen.


Introduction
The spin-dependent and spin-independent structure functions of T µν , the forward virtualphoton Compton tensor of the proton, carry important information about the QCD dynamics. They test the Euclidean region of the theory since Q 2 ≡ −q 2 > 0. For Q 2 ∼ m 2 π = 0, the behavior of T µν is determined by the chiral theory, and can be obtained within a chiral expansion using Heavy Baryon Effective Theory (HBET) [1]. If one works within a large N c ideology (where N c is the number of colours) the Delta particle should be incorporated in the HBET Lagrangian [2] (see also [3,4]). This produces a double expansion in ∼ m π /m ρ and ∼ ∆/m ρ , where ∆ = M ∆ − M N . Note that this creates a new expansion parameter m π /∆ ∼ 1/2; the associated corrections will be incorporated in our computation together with the pure chiral result.
Within this framework we compute the spin-dependent and spin-independent structure functions of the forward virtual-photon Compton tensor of the proton at one loop in Heavy Baryon Chiral Perturbation Theory (HBχPT) including the Delta particle. T µν cannot be directly related to cross sections obtained at fixed energies, as it tests the Euclidean regime. Nevertheless, it is possible to obtain it (up to eventual subtractions) from experiment through dispersion relations, i.e., through specifically weighted averages of measured cross sections over all energies. Possible constructions are the so-called generalized sum rules, which, for large energies, can be related with the deep inelastic sum rules. These have been studied in Ref. [5] for the spin-dependent case. The spin-independent case has been briefly discussed in Ref. [6]. We will not enter into this interesting line of research in this paper.
Instead, our main motivation for obtaining the chiral structure of T µν is that T µν appears in the matching computation between HBET and non-relativistic QED (NRQED) that determines c pl i 3 and c pl i 4 (l i = e or µ), the Wilson coefficients of the lepton-proton four-fermion operators in the NRQED [7] Lagrangian. As soon as hadronic effects start to become important in atomic physics, these Wilson coefficients play a major role. They appear in the hyperfine splitting (spin-dependent) and Lamb shift (spin-independent) in hydrogen and muonic hydrogen (see Refs. [8,9,6]). Therefore, their determination allows us to relate the energy shifts obtained in hydrogen and muonic hydrogen. Even more important, these Wilson coefficients usually carry most of the theoretical uncertainty in these splittings. This is particularly so in the case of the muonic hydrogen Lamb shift. At present, it is the limiting factor for improving the precision of the determination of the electromagnetic proton radius from the measurements taking place at PSI [10,11] of the muonic hydrogen spectra. This necessity to improve our knowledge (of the spin-independent) lepton-proton four-fermion Wilson coefficient has led us to compute this quantity in HBχPT including the Delta particle. Fortunately enough, this object is chiral enhanced. Therefore, the oneloop chiral computation yields a pure prediction, without the need of new counterterms, of ∆E TPE , the (hadronic) two-photon exchange contribution to the muonic hydrogen Lamb shift: ∆E L = E(2P 3/2 ) − E(2S 1/2 ). These results have been used in the recent determination of the muonic hydrogen Lamb shift and the proton radius performed in Ref. [12]. One of the main motivations of this paper is to give the details of the hadronic-related part of that analysis.
We profit this analysis to revisit the distinction between Born and non-Born terms of T µν and ∆E TPE . Such distinction produces the so-called Zemach (or Born) and polarizability

Effective Field Theories
In this section, we will present the main building blocks of the HBET and NRQED Lagrangians needed for our analysis (see also Ref. [8]).

HBET
Our starting point is the SU(2) version of HBET coupled to leptons where the Delta particle is kept as an explicit degree of freedom. The degrees of freedom of this theory are the proton, neutron and Delta, for which the NR approximation can be taken, and pions, leptons (muons and electrons) and photons, which will be taken relativistic. This theory has a cut-off µ << M p , m ρ , which is much larger than any other scale in the problem. The Lagrangian can be split in several sectors. Nevertheless, the fact that some particles will only enter through loops, since only some specific final states are wanted, simplifies the problem. The Lagrangian can be written as an expansion in e and 1/M p and can be structured as follows L HBET = L γ + L l + L π + L lπ + L (N,∆) + L (N,∆)l + L (N,∆)π + L (N,∆)lπ , (2.1) representing the different sectors of the theory. In particular, the ∆ stands for the Delta particle: the spin 3/2 baryon multiplet (we also use ∆ = M ∆ − M p , the specific meaning in each case should be clear from the context). The photonic Lagrangian reads (the first corrections to this expression scale like α 2 /M 4 p ) where d 2,R stands for the hadronic contribution. The second term will not be considered any further in this paper, since we are mainly interested in the lepton-proton four-fermion operators.
The leptonic sector can be approximated to (iD µ = i∂ µ − eA µ ) where l i = e, µ.

2
The Lagrangian of a heavy baryon at O(1/M 2 p ) coupled to electromagnetism reads where iD 0 p = i∂ 0 + Z p eA 0 , iD p = i∇ − Z p eA. For the proton Z p = 1 (for the neutron Z p = 0 and for all indices p → n).
The Delta particle mixes with the nucleons at O(1/M p ) (O(1/M 2 p ) terms are not needed in our case). The only relevant interaction in our case is the p-∆ + -γ term, which is encoded in the second term of where T stands for the delta 3/2 isospin multiplet, N for the nucleon 1/2 isospin multiplet and the transition spin/isospin matrix elements fulfill (see [13]) The baryon-lepton Lagrangian provides new terms that are not usually considered in HBET. The relevant term in our case is the interaction between the leptons and the nucleons (actually only the proton): The above matching coefficients fulfill c pl i 3,R = c p 3,R and c pl i 4,R = c p 4,R up to terms suppressed by m l i /M p , which will be sufficient for our purposes.
Let us note that with the conventions above, N p is the field of the proton (understood as a particle) with positive charge if l i represents the leptons (understood as particles) with negative charge.
The hadronic interactions are organized according to their chiral counting. Since a single chiral loop already produces a factor 1/(4πF 0 ) 2 ∼ 1/M 2 p , we only need the free-particle pionic Lagrangian: We do not need to account for pion self-interactions, and the pion-baryon interactions are only needed at O(m π ), the leading order, which is known [14,15]: where, T µ a is the Rarita-Schwinger spin 3/2 field and S µ = i 2 γ 5 σ µν v ν is the spin operator (where we take v µ = (1, 0)).
This finishes all the needed terms for this paper, since the other sectors of the Lagrangian would give subleading contributions.

NRQED(µ)
In the muon-proton sector, by integrating out the m π and ∆ scales, an effective field theory for muons, protons and photons appears. In principle, we should also consider neutrons but they play no role at the precision we aim. The effective theory corresponds to a hard cut-off ν << m π and therefore pions and Deltas have been integrated out. The Lagrangian is equal to the previous case but with neither pions nor Deltas, and with the following modifications: N µ , where it is made explicit that the the muon has become NR. Any further difference goes into the Wilson coefficients, in particular, into the Wilson coefficients of the baryon-lepton operators. In summary, the Lagrangian reads with the following definitions: iD 0 µ = i∂ 0 − Z µ eA 0 , iD µ = i∇ + Z µ eA and Z µ = 1. L e stands for the relativistic leptonic Lagrangian in Eq. (2.3) and L N e for Eq. (2.7), both for the electron case only.
Our main interest is the determination of c

NRQED(e)
If we focus in the electron-proton sector, things go quite as in the previous section. After integrating out scales of O(m π , ∆), an effective field theory for electrons coupled to protons (and photons) appears. This effective theory has a cut-off ν << m π and pions, Deltas and muons have been integrated out, but the electron is still relativistic. After integrating out scales of O(m e ) in the electron-proton sector, we still have an effective field theory for electrons coupled to protons and photons. Nevertheless, now the electrons are NR. The Lagrangian is quite similar to the one in Subsec. 2.2 but without a light fermion and with the replacement µ → e. It reads We will perform the matching to this theory directly from HBET. At O(α 2 ) this matching can be symbolically represented by the same figure as in the case of the muon, namely Fig. 1.

Forward virtual Compton tensor T µν
The electromagnetic current reads J µ = i Q iqi γ µ q i , where i = u, d (we will not consider the strange quark in this paper) and Q i is the quark charge. The form factors (which we will understand as pure hadronic quantities, i.e. without electromagnetic corrections) are then defined by the following equation: where q = p ′ − p and F 1 , F 2 are the Dirac and Pauli form factors, respectively. The states are normalized in the following (standard relativistic) way: where s is an arbitrary spin four-vector obeying s 2 = −1 and p · s = 0. More suitable for a NR analysis are the Sachs form factors: Nevertheless, the main object of interest of this paper is the forward virtual-photon Compton tensor, which has the following structure (ρ = q · p/M p ≡ v · q, although we will usually work in the rest frame where ρ = q 0 ): It depends on four scalar functions, which we call structure functions. We split the tensor into the symmetric (spin-independent), T µν S = T νµ S (the first two terms of Eq. (3.6)), and antisymmetric (spin-dependent) pieces, T µν A = −T νµ A (the last two terms of Eq. (3.6)). We have computed this tensor at one loop in HBχPT. The diagrams that contribute are listed in Figs. 2, 3 and 4 (without closing the loop with the muon). The first figure refers to diagrams without Delta contributions (pure chiral), the second to the tree-level Delta contribution, and the last to one-loop chiral diagrams involving the Delta particle. Expressions in D = 4 − ǫ and four dimensions for each diagram can be found in Appendix C. Summing them up we can reconstruct the tensor structure of T µν (in other words, check gauge invariance). In principle, more diagrams, besides those drawn should be considered but they do not contribute to the structure functions at the order we aim in this work.
It is also common to split T µν into two components, which we label "Born" and "pol": The Born term is defined as the contribution coming from the intermediate state being the proton (somewhat the elastic contribution). The associated structure functions can be For the spin-dependent case, the only contribution is the term proportional to G M , which comes from the A Born 1 term (this is the only term that contributes to the Born (Zemach) piece of the hyperfine splitting). For the spin-independent case we only need G Following common practice we define the electromagnetic charge density as The inverse of its Fourier transform allows us to obtain the even powers of the moments of the charge distribution of the proton, The odd powers are obtained (defined) through the relation: . (3.19) An analytic expression of this quantity is relegated to Eq. (4.12). Note that, by using dimensional regularization, we can eliminate all the terms proportional to integer even powers of q 2 in this expression. For k > 1, this integral is dominated by the chiral result and can be approximated by (3.20) Finally, let us note that, by construction, both T µν Born and T µν pol comply with current conservation. The separation (definition) of the Born and polarizability terms is in general ambiguous, see, for instance, the discussion in Refs. [19,20]. In our case, as far as we give an explicit definition for T µν Born , this ambiguity disappears. In what follows we consider the computation of T µν pol .

Computation of T µν pol
We split each S pol i /A pol i in the following way: S pol i,π and A pol i,π encode the contributions only due to pions. They are produced by the diagrams listed in Fig. 2. Summing them up we can reconstruct the tensor structure of T µν . In D dimensions the structure functions read

24)
A pol 2,π (q 2 , q 0 ) = where the loop functions J i have been defined in D-dimensions in Eq. (B.1). These structure functions reduce to the following expressions in D = 4: and For D = 4 we can compare with previous results in the literature. S pol 1,π and S pol 2,π were originally computed in [6]. We agree with those results, which were obtained with different methods, either by dispersion relations or through a diagrammatic computation assuming gauge invariance. In the case of real photons (for q 2 = 0 in the Coulomb gauge) we recover the results of [16]. S pol 1,π has also been checked in the limit q 0 → 0 in Ref. [21], and S pol 1/2,π for all q 0 and q 2 in Ref. [22].

11
The spin-dependent structure functions, A pol 1,π and A pol 2,π , agree with the ones given in Eqs. (30) and (34) of [5], up to a normalization factor. They follow from summing up all the contributions of the diagrams in Fig. 2 that have an antisymmetric contribution, i.e. diagrams (2), (4) and (5) of Fig. 2.
We now move to contributions involving Delta particles. We first consider tree-level Delta mediated contributions. The corresponding diagram is pictured in Fig. 3, and the associated contributions read: Eq. (3.32) agrees with [6] and, in the limit q 0 → 0, with the leading order expression of [21] up to normalization. Eq. (3.32) differs from the expression obtained in Ref. [6] using dispersion relations by a local term. For the spin-dependent terms we are in agreement with [5].
The last set of diagrams that we consider are those with one internal chiral loop and virtual Delta particles. They are drawn in Fig. 4 producing the following D-dimensional expressions for the structure functions: The results for D = 4 dimensions are: where we have defined Z as The D = 4 expressions for S pol 1,π∆ and S pol 2,π∆ agree with Eqs. (51) in [15] for the case of real photons, i.e. q 2 = 0 and in the Coulomb gauge.
Summing up all the contributions of the diagrams in Fig. 4 which have an antisymmetric contribution, i.e. diagrams (2), (4) and (5), we get the spin-dependent part that agrees with Eqs. (33) and (36) in [5], up to a normalization factor.
In all expressions we use principal value prescriptions, the Dirac delta contributions associated to the propagators have gone into the Born term. Nevertheless, from the point of view of the effective theory this splitting between the polarizability and Born term is quite arbitrary.
4 Matching HBET to NRQED: The matching between HBET and NRQED can be performed in a generic expansion in 1/M p , 1/m µ and α. We have two sorts of loops: chiral and electromagnetic. The former are always associated to 1/(4πF 0 ) 2 factors, whereas the latter are always suppressed by α factors. Any scale left to get the dimensions right scales with m π or ∆. In our case we are only concerned with obtaining the matching coefficients of the lepton-baryon operators of In what follows, we will assume that we are doing the matching to NRQED(µ). Therefore, we keep the whole dependence on m l i /m π . The NRQED(e) case can then be derived by expanding m e versus m π .
At O(α 2 ), the contribution to c pl i 3 (see Fig. 1) from matching HBET to NRQED can be written in a compact way in terms of the structure functions of the forward virtual-photon Compton tensor. It reads [23] This result keeps the complete dependence on m l i and is valid both for NRQED(µ) and NRQED(e). This contribution is usually organized in the following way This goes beyond the aimed accuracy of our calculation and so we neglect c pl i 3,R . The second term in Eq. (4.2) corresponds to Eq. (4.1) assuming the proton to be point-14 like. With the precision needed it reads in the MS scheme (see [24]  The third term in Eq. (4.2) is generated by the spin-independent Born contribution to T µν in Eq. (3.12). We symbolically picture it in Fig. 5. At leading order in the NR expansion it reads 3 Note again that this result holds for both NRQED(e) and NRQED(µ). In other words, the exact dependence on m l i is kept (at leading order in the NR expansion). The linear dependence in the lepton mass makes this contribution much smaller for the case of hydrogen. G (0) E = 1. We take the expression for G E from Eq. (3.13). The use of effective field theories and dimensional regularization is a strong simplification, which we have already used when writing Eq. (4.5). This guarantees that only low energy modes contribute to the integral, and 2 In this expression we have computed the loop with the proton being relativistic to follow common practice. Nevertheless, this assumes that one can consider the proton to be point-like at the scales of the proton mass. To stick to an standard EFT approach one should consider the proton to be NR. Then one would obtain The difference between both results is of the order of c 3,R , and gets absorbed into this coefficient (which we do not know anyhow). Therefore, the value of c pli 3 , will be the same no matter the prescription used. In practice there could be some difference due to truncation, but always of the order of the error of our computation. 3 In Ref. [9] we named this object c pli 3,Zemach .
15 that we only need the non-analytic behavior of G (2) E in q 2 around m π and ∆. In other words, even though some point-like contributions are still encoded in G (2) E , they do do not contribute to the integral. The analytical behavior in q 2 produces scaleless integrals, which are zero in dimensional regularization. This is a reflection of the factorization of the different scales. Therefore, we do not need to introduce the point-like interactions to regulate the infrared divergences of the integrals at zero momentum, as it is done if trying to compute this object directly from the experimental data. We will come back to this issue when we discuss the Zemach moments.
The computation of c pl i 3,Born was made in Ref. [9]. Here we give a simplified expression: where the first line is due to scales of O(m π ) and the terms proportional to B n are due to scales of O(∆), where (this corrects Eq. (61) of Ref. [9]) H n is the n harmonic number, and Γ(n) is the Euler Γ function. Eq. (4.6) encapsulates all the non-analytic dependence in the light quark masses and in the splitting between the nucleon and the Delta mass (proportional to powers of 1/N c in the large N c limit) of c pl i 3,Born . This expression is the leading contribution to the Zemach term in the chiral counting (supplemented with a large N c counting). This is a model independent result. Other contributions to the Zemach term are suppressed in the chiral counting. c pl i 3,Born can be related with (one of) the Zemach moments: The Zemach moments can be determined in a similar way as the moments of the charge distribution of the proton. For even powers we have the relation 4 The odd powers are obtained (defined) through the relation: .
(4.10) Again, using dimensional regularization, we can eliminate all the terms proportional to integer even powers of q 2 in this expression. For k ≥ 1 this integral is dominated by the chiral result and can be approximated by (4.11) It is possible to get an analytic result for these integrals. We obtain (y ≡ mπ ∆ ) In Table 1 we give our predictions for some selected charge and Zemach moments, 5 both in the effective theory with only pions and in the effective theory with pions and Deltas. The even powers are obtained by direct Taylor expansion of Eq. (3.13). The odd powers are obtained from Eq. (4.12). We have also numerically checked the values of r 2k+1 directly using Eq. (3.19). In order to estimate the error of the charge/Zemach moments and the other quantities we compute in this paper we proceed as follows. We count m π ∼ Λ QCD m q and ∆ ∼ Λ QCD Nc . We then have the double expansion mπ Λ QCD ∼ mq Λ QCD and ∆ Λ QCD ∼ 1 Nc . We still have to determine the relative size between m π and ∆. We observe that m π /∆ ∼ N c mq Λ QCD ∼ 1/2. Therefore, we associate a 50% uncertainty to the pure chiral computation. For all Zemach moments we observe good convergence, with the contribution due to the Delta being much smaller than the pure chiral result, and well inside the 50% uncertainty. Leaving aside the Delta, the splitting with the next resonances suggest a mass gap of order Λ QCD ∼ 500-770 MeV depending on whether one considers the Roper resonance or the ρ. Therefore, we assign  Table 1: Values of r n in fermi units. The first two rows give the prediction from the effective theory: the first row for the effective theory with only pions and the second for the theory with pions and Deltas. The third row corresponds to the standard dipole fit of Ref. [25] with r 2 = 0.6581 fm 3 . The fourth and fifth rows correspond to different parameterizations of experimental data [26,27], with the latest fit being the more recent analysis based on Mainz data. For completeness, we also quote r 3 (2) = 2.71 fm 3 from Ref. [28].
mπ Λ QCD ∼ 1/3 and ∆ Λ QCD ∼ 1/2, as the uncertainties of the pure chiral and the Delta-related contribution respectively. We add these errors linearly for the final error. This gives the expected size of the uncomputed corrections but numerical factors may change the real size of the correction. In particular, huge discrepancies with these estimates may signal the failure of HBχPT for obtaining some of the observables considered in this paper.
The chiral prediction is expected to give the dominant contribution of r n for n ≥ 3. For n = 2 it could also give the leading chiral log. For smaller n the chiral corrections are subleading. Note that for all n ≥ 3, these expressions give the leading (non-analytic) dependence in the light quark mass as well as in 1/N c . This is a valuable information for eventual lattice simulations of these quantities where one can tune these parameters. In Table 1 we also compare with the standard dipole ansatz [25], and with different determinations using experimental data of the electric Sachs form factor fitted to more sophisticated functions [26,27]. 6 The latest fit claims to be the more accurate. Nevertheless, we observe large differences, bigger than the errors. This is specially worrisome for large n, since the chiral prediction is expected to give the dominant contribution of r n for n ≥ 3. In this respect, we believe that the chiral result may help to shape the appropriated fit function and, thus, to discriminate between different options, as well as to assess uncertainties. The impact of choosing different fit functions can be fully appreciated, for instance, in the different values of the electromagnetic proton radius obtained in Ref. [29] versus Refs. [30,31] from direct fits to the ep scattering data. Such values differ by around 3 standard deviations. On the other hand, even if on general grounds one may expect the charge/Zemach moments will be more and more sensitive to the chiral region for n → ∞, large fractions of the experimental numbers are determined by the subtraction terms included to render these objects finite (for odd powers of n). We stop the discussion here but the reason for such large discrepancies should be further investigated.
As we have already mentioned, c pl i 3,Born can be related with (one of) the Zemach moments: Note again that the terms proportional to "1" and r 2 vanish in dimensional regularization.
We can now obtain (m r = m µ M p /(m µ + M p )) the Born contribution to ∆E TPE , from the effective field theory. We quote our results in Table 2. The pure chiral result was already obtained in Ref. [9]. The π&∆ result corrects the evaluation made in that reference due to the error in its Eq. (61). Note that the new result is much more convergent, since the correction associated to the Delta is much smaller. On the other hand our result is now much more different with respect to standard values obtained from dispersion relations. We quote two of them in Table 2. One may wonder whether such difference is due to relativistic corrections. An estimate of the relativistic effects can be obtained from the analysis made in Ref. [32], which, however, is based on dipole form factors parameterizations. The difference between the relativistic and NR expression was found to be small (∼ 3µeV). It should be checked that this feature holds with different parameterizations. If so, the difference seems to be mainly due to the computation of the Zemach correction (see Table 1 and the discussion above). Therefore, as stated above, the reason for such large discrepancies should be investigated. In the mean time we will stick to our model independent prediction from the effective theory.  Table 2: Predictions for the Born contribution to the n = 2 Lamb shift. The first two entries correspond to dispersion relations. The last two entries are the predictions of HBET: The 3rd entry is the prediction of HBET at leading order (only pions) and the last entry is the prediction of HBET at leading and next-to-leading order (pions and Deltas).

3,pol
Finally, we consider the polarizability correction. It is obtained from Eq. where For the case of only pions we havem π = m µ /m π and For the case of Delta at tree level we havem ∆ = m µ /∆ and For the case of loops including the Delta we havem ∆π = m µ /∆ and where t = m π /∆. Note that the imaginary part of these expressions comes only from the Wick rotation of k 0 and will vanish upon integration. The pure pion contribution was already found in Ref. [6]. Our full prediction for the polarizability term including the Delta effects reads ∆E pol = c   Table 3 we compare our determination with previous results. Most of them are obtained by a combination of dispersion relations plus some modeling of the subtraction term that we discuss below. The analysis of Ref. [22] has a different status. In this reference the polarizability correction was computed using BχPT with only pions. Such computation treats the baryon relativistically. The result incorporates some subleading effects, which are sometimes used to give an estimate of higher order effects in HBχPT. Nevertheless, the computation also assumes that a theory with only baryons and pions is appropriate at the proton mass scale. This should be taken with due caution. Still, it would be desirable to have a deeper theoretical understanding of this difference, which may signal that relativistic corrections are important for the polarizability correction. In any case, the BχPT computation differs of our chiral result by around 50% (this means around 1.5 times the error we use for the chiral contribution, once the Delta is incorporated in the calculation), which we consider reasonable.
It is also worth discussing the LEX approximation used in Ref. [22]. This approximation consists in setting q 0 = 0 everywhere except in the denominator in Eq. (4.1). For the pure chiral result, this approximation works remarkable well (18.51(exact) vs 17.85 (LEX)). Nevertheless, such success does not survive the incorporation of the ∆ particle. For the Delta tree-level contribution we find (-1.58(exact) vs 0 (LEX)). The real problem appears from the one-loop pion-Delta result. For such contribution there are 1/q 0 singularities in the tensor that only cancel if the complete expression is used. Doing the LEX approximation leads to divergent expressions. Even more worrisome is the fact that, at present, there are no theoretical justification for using the LEX approximation for the integral in Eq. (4.1). It is not correct to assume that the photon energy that appears in the integral, q 0 , corresponds to the energy in the atomic system. It rather reflects virtual fluctuations of order of the pion and muon mass (as well as of the ∆ scale). Since those particles are relativistic at those scales it is theoretically incorrect, a priori, to neglect q 0 . In any case, on the light of the good agreement for the pure chiral case, it would be interesting to see whether one could find a theoretical justification for such behavior.
It is also interesting to consider the limit m l i ≪ m π , which is relevant for the hydrogen atom. In this limit Eq. (4.15) approximates, with logarithmic accuracy, to (µeV) DR + Model [33] [34] [35] [36] BχPT [22](π) HBET [6](π) [12](π&∆) ∆E pol 12 (2) Table 3: Predictions for the polarizability contribution to the n = 2 Lamb shift. The first four entries use dispersion relations for the inelastic term and different modeling functions for the subtraction term. The number of the fourth entry has been taken from [22]. The 5th entry is the prediction obtained using BχPT. The 6th entry is the prediction of HBET at leading order (only pions) and the last entry is the prediction of HBET at leading and next-to-leading order (pions and Deltas).
These logs can be obtained by computing the ultraviolet behavior of the diagram in Fig. 6. This contribution is proportional to c A 1 and c A 2 or, in other words, the polarizabilities of the proton (see [37,38]). For the pure pion cloud, the polarizabilities were computed in Ref. [16]. The contribution due to the ∆ can be found in Ref. [39]. The scale in the logarithm is compensated by the next scale of the problem, which can be m π or ∆. For contributions which are only due to the ∆ or pions, the scale is unambiguous. In the case where pions and ∆ are both present in the loop we will choose the pion mass (the difference being beyond the logarithmic accuracy). It is known that the pure chiral prediction of α M . Nevertheless, this object is comparatively small, and even more so for 5α M , the combination that appears in the logarithmic approximation. Whereas the experimental number reads 5α [40], the pure chiral result gives (5α M )(π) ≃ 60 × 10 −4 fm 3 , and after the inclusion of the Delta we obtain (5α Again the inclusion of the Delta deteriorates the agreement but the difference is of the order of one sigma according to our error analysis. We take this as an indication that effective field theory result will not be very far off from the real number for the case of muonic hydrogen and that, maybe, the pure chiral result compares better with experiment than after the inclusion of the Delta. Nevertheless, we will not make any assumption in this respect and stick to the complete prediction of the effective theory. It is also customary to split the polarizability term (note that the Born term has already been subtracted from it) in what is called the inelastic and subtraction term: It is argued that the inelastic term does not require further subtractions and can be obtained through dispersion relations. On the other hand the subtraction term can not be directly obtained from experiment. This fact has been used in Ref. [41] to emphasize that the polarizability term is affected by huge theoretical uncertainties. In this paper, we can avoid making any assumption about the dispersion relation properties of these quantities. This is possible within the framework of effective field theories. In this setup the splitting between the inelastic and subtraction terms is unmotivated, and to some extent artificial (as it was the splitting between the Born and polarizability term). Let us elaborate on this point and see what effective field theories have to say in this respect. The main problem comes, as it has already been pointed out in Ref. [22], from the diagram in Fig. 3. This diagram yields a finite (an small) contribution to c pl i 3,pol (and therefore to the energy shift, see Eq. (4.23)). Nevertheless, when splitted into c pl i 3,inel and c pl i 3,sub , each term diverges in the following way If we set the ultraviolet cutoff to the ρ mass, ν = m ρ , the energy shift of each term is one order of magnitude bigger ∼ −11.37 µeV than the exact result for the sum. Obviously such contribution is fictitious and may alter the value of the individual terms. On the other hand it is possible to perform this splitting for the case of the pion and pion-Delta loop. We obtain the following: They are of the same magnitude. Their size is barely one order of magnitude smaller than the total polarizability term. For the case of the pion loop it is possible to obtain analytic expressions in the limit m µ = m π , which is a rather good approximation: where G ≃ 0.9160 is the Catalan's constant. For these quantities the LEX approximation works quite well, both for the pion and the pion-Delta loop case. We find 7 which is again asking for a theoretical explanation of this relatively good agreement. For comparison we show different values obtained for the subtraction and inelastic term obtained in the literature in Table 4.  Table 4: Values for the subtraction and inelastic terms that one can find in the literature. (1) This number is the adjusted value of Ref. [36], given in [22].

∆E
We now combine the contribution from the Born and polarizability term and summarize our final results for ∆E TPE : ∆E TPE = ∆E Born + ∆E pol = 28.59(π) + 5.86(π&∆) = 34.4(12.5)µeV . (4.33) We would like to emphasize that this result is a pure prediction of the effective theory. It is also the most precise expression that can be obtained in a model independent way, since O(m µ α 5 m 3 µ Λ 3 QCD ) effects are not controlled by the chiral theory and would require new counterterms. Our number is only marginally bigger than ∆E TPE = 33(2)µeV [21]. This number is the one used in Ref. [11] for its determination of the proton radius. It is obtained as the sum of the elastic and inelastic terms from Ref. [35] and the subtraction term from Ref. [21]. Note that this evaluation is model dependent. Even though the low energy behavior of the forward virtual Compton tensor was computed to O(p 4 ), this does not reflect in an improved determination of the polarizability correction, since an effective dipole form factor is used, not only at the ρ mass scale, but also at the chiral scale. This problem also introduces a model dependence in its error estimate. Other existing determinations [33,34,35] yield quite similar numbers but suffer from the same systematic uncertanties. In this respect our calculation is model independent and have completely different systematics. The fact that we obtain similar numbers is comforting for the reliability of the proton radius determinations obtained in Refs. [11,12]. On the other hand one should not forget that the individual contributions are quite different, and the reasons for that should be further investigated.
Yet it is quite remarkable that the total sum gives such similar numbers. 7 For m µ = m π an analytic expression can be found for the pion-loop case [22]: We proceed in the same way as in the spin-independent case. We will assume that we are doing the matching to NRQED(µ). Therefore, we keep the whole dependence on m l i /m π . The NRQED(e) case can then be derived by expanding m e versus m π . We match HBET and NRQED order by order in a generic expansion in 1/M p , 1/m µ and α. We have two sorts of loops: chiral and electromagnetic. The former are always associated to 1/(4πF 0 ) 2 factors, whereas the latter are always suppressed by α factors. Any scale left to get the dimensions right scales with m π or ∆. At O(α 2 ), the contribution to c pl i 4 (see Fig. 1) from matching HBET to NRQED can be written in a compact way in terms of the structure functions of the forward virtual-photon Compton tensor. In Euclidean space it reads consistent with the expressions obtained long ago in Ref. [42]. This result keeps the complete dependence on m l i and is valid both for NRQED(µ) and NRQED(e), i.e. for hydrogen and muonic hydrogen. Similarly to the spin-independent case, this contribution can be organized in the following way Within the effective field theory framework the contribution from energies of O(m ρ ) or higher in Eq. (5.1) are encoded in c pl i 4,R ≃ c p 4,R (analogously to c 3 ). The other terms (associated to energies of O(m π )) were computed with O(α 2 × (ln m q , ln ∆, ln m l i )) accuracy in Ref. [8]. We quote them here for ease of reference 8 (c (p) Summing up the three terms one has Parametrically, the three contributions, Eqs. (5.3), (5.4) and (5.6), are of the same order. Nevertheless, the polarizability and the point-like term are much smaller. This is consistent with the fact that polarizability correction seems to be small [44,45,46], if determined through dispersion relations. As already discussed in Ref. [8], the effective field theory computation gives a double explanation to this fact. On the one hand, this is due to the smallness of the numerical coefficient of the polarizability term, but there also seems to be some large N c rationale behind. Since g πN ∆ = 3/(2 √ 2)g A in the large N c limit, the polarizability term vanishes (see [5]) except for the tree-level-like Delta contribution (the last term in Eq. (5.6)). Nevertheless, the latter also vanishes against the κ p -dependent point-like contribution (which effectively becomes the result of a point-like particle) in the large N c limit, since b F 1 = 3/(2 √ 2)κ V and κ p = κ V /2 [47]. Note also that the point-like term and the tree-level-like Delta contribution are suppressed by 1/π factors with respect the Born contribution.
This discussion also illustrates that splitting the total contribution into different terms may introduce spurious effects that vanish in the total sum. We have also seen a similar thing but in a different context for the case of the spin-independent computation.
Our computation allows us to relate c Note that we have already used the fact that c pl i 4,Born cancels in the difference, as it is independent of the lepton mass. The experimental and theoretical results discussed before suggest that c pl i 4,Born is the leading contribution to the Wilson coefficient. Therefore, such contribution can be obtained from c ple 4 , which can be determined from the hyperfine splitting of Hydrogen. In Ref. [8] it was estimated to be c ple 4 ≃ −48α 2 . By considering differences in Eq. (5.9) the ultraviolet behavior gets regulated and the logarithmic divergences vanish. This makes these contributions to be very small an negligible compared with the uncertainties. For the point-like contribution we obtain Overall we obtain c plµ 4 ≃ −46α 2 . The bulk of this contribution is expected to come from the Born term, which in turn is related to the Zemach magnetic radius, by the following relation The chiral log result compares well (∼ 30%) with existing predictions (∼ 1.04-1.08 fm) from hydrogen hyperfine [48,49], from dispersion relations [28,27], or from the muonic hydrogen hyperfine [11]. Note that in the case of the determinations of r Z from the hyperfine splitting (either from hydrogen or muonic hydrogen) one needs to control the relativistic hadronic affects associated to the Born term as well as the polarizability correction. On the other hand if we are only interested in the hyperfine splitting it may make more sense to consider c pl i 4 as a whole. We relegate a more detailed discussion to future work.

Conclusions
We have computed the spin-dependent and spin-independent structure functions of the forward virtual-photon Compton tensor of the proton at one loop in HBχPT including the Delta particle. We have given D-dimensional expressions too. Those are relevant for future higher order loop computations. We have compared our results with previous computations. The D = 4 expressions for the spin-dependent structure functions were computed in [5]. We agree with their results. The D = 4 expressions for the pure chiral (without Delta contributions) spin-independent structure functions were computed in [6]. We agree with their results too. The Delta-associated contributions to the spin-independent structure functions are new. We also profit to present all these results obtained throughout the years in a unified form.
We have used these results to determine the leading chiral and large N c structure of c pl i 3 and c pl i 4 , or, in other words, to determine their non-analytic dependence on m q and N c . The fact that we have full control over the quark mass dependence makes our result very useful for eventual lattice determinations of these quantities. By fine tunning the mass in simulations we can identify the results obtained in this paper and up to which mass the chiral is good approximation. One could also vary N c to check the theory.
These Wilson coefficients appear in the hyperfine splitting (spin-dependent) and Lamb shift (spin-independent) in hydrogen and muonic hydrogen. c pl i 3 , the relevant Wilson coefficient for the Lamb shift, is chiral enhanced. Therefore, the one-loop chiral result is a pure prediction of the effective theory, which we use to determine ∆E TPE = 28.6(π) + 6.1(π&∆) = 34.4(12.5)µeV , (6.1) the energy shift associated to the (hadronic) two-photon exchange of the Lamb shift. These results have been used in the recent determination of the muonic hydrogen Lamb shift and the proton radius performed in Ref. [12]. We would like to emphasize that Eq. (6.1) is the most precise expression that can be obtained in a model independent way, since O(m µ α 5 m 3 µ Λ 3 QCD ) effects are not controlled by the chiral theory and would require new counterterms. Our final number is quite similar to previous estimates existing in the literature. Nevertheless, those computations require the splitting of the two-photon contribution into different terms. Some of them are then computed using different dispersion relations, whereas one last term requires modeling its Q 2 dependence. In contrast, we have used the same method for all computations contributing to our result, yielding a parameter-free prediction. On the other hand one should not forget that the individual contributions are quite different, and the reasons for that should be further investigated. In this respect we have discussed what the effective theory has to say about the separation into Born, polarizability, inelastic and subtraction term. The Born contribution is related with the Zemach moments. In this paper we have also given the prediction of the effective theory for some charge, r n , and Zemach, r n (2) moments. Finally, we have also discussed the chiral dependence of the spin-dependent four-fermion Wilson coefficient, c pl i 4 , and obtained the relation between c pe 4 and c pµ 4 given by the effective theory.

B Master Integrals
We follow the notation of [14,15] and assume a negative infinitesimal imaginary part for all the propagators.
where in D-dimensions: and in D = 4 − ǫ: For the function J 0 we get when |ω| < m and for the case where ω < −m we get the analytically continued function All the other functions are related to Eq. (B.7)/Eq. (B.6) and Eq. (B.4) by: J 7 (q 0 , m) = m 2 J 2 (q 0 , m) + (D + 2)J 6 (q 0 , m). (B.14) We also define the derivative function

C Amplitudes for the diagrams
Throughout this work we use the normalizationū(p)u(p) = 2M p and we define ∆ = M ∆ − M N .m 2 = m 2 π − q 2 x(1 − x) and the function Z has been defined in Eq. (3.44). We work in the rest frame where v = (1, 0).

C.1 Pion loops
Here we collect the amplitudes of all the diagrams contributing to the proton polarizability through a loop of pions, represented in Fig. 2, plus the ones with a crossed photon lines or permutations, which are assumed to be implicit in the representation. For all the diagrams here we consider the overall factor A = 2M p We assume a positive infinitesimal imaginary part for the propagators of h 14 − h 19 . Diagrams with only 1 pion are zero due to the fact that we are working in the static limit.

C.2 Pion loops which include a ∆ excitation
Here we collect the amplitudes of all the diagrams contributing to the proton polarizability with a ∆ particle and through a loop of pions, represented in Fig. 4, plus the ones with a crossed photon lines or permutations, which are assumed to be implicit in the representation. For all the diagrams here we consider the overall factor A = − 8 3 M p g 2 πN∆ h ∆ 4 (q 2 , q 0 ) =