A reassessment of nuclear effects in muonic deuterium using pionless effective field theory at N3LO

We provide a systematic assessment of the order-$\alpha^5$ nuclear contributions to the Lamb shift of muonic deuterium, including the accompanying radiative corrections due to vacuum polarization, up to next-to-next-to-next-to-leading order (N3LO) within the pionless effective field theory (EFT). We also evaluate higher-order corrections due to the single-nucleon structure, which are expected to be the most important corrections beyond N3LO. We find a correlation between the deuteron charge and Friar radii, which can be useful to judge the quality of charge form factor parametrisations. We refine the theoretical description of the $2\gamma$-exchange contribution, especially in the elastic contribution and the radiative corrections, ameliorating the original discrepancy between theory and experiment in the size of $2\gamma$-exchange effects. Based on the experimental Lamb shift of muonic deuterium, we obtain the deuteron charge radius, $r_d(\mu\text{D})=2.12763(13)_\text{exp}(77)_\text{theory}$~fm, which is consistent with (but less precise than) the value obtained by combining the H-D isotope shift with the muonic hydrogen Lamb shift. The theory uncertainty is evaluated using a Bayesian procedure and is dominated by the truncation of the pionless EFT series.


I. INTRODUCTION
Recent advances in the spectroscopy of light muonic atoms, by the CREMA Collaboration at the Paul Scherrer Institute, offer a new leap in precision for studies of nucleon and nuclear structure [1]. Their determinations of the charge radii of the proton, r p , and the deuteron, r d , from the Lamb shift of muonic hydrogen (µH) [2,3] and muonic deuterium (µD) [4], have been particularly intriguing. The former disagreed, quite spectacularly, with the previous determinations of r p using the conventional methods of hydrogen (H) spectroscopy and electron-proton (ep) scattering, viz., the "proton-radius puzzle" (see, e.g., [5,6]). The µD determination of r d disagrees with CODATA-2014 [7], but also with a more recent determination from deuterium (D) [8]. Initially, it even disagreed with µH, given the precise determination of the radius difference, r 2 d − r 2 p , from the H-D isotope shift. The latter discrepancy has recently been resolved on the theory side by adding several missing contributions to the µD Lamb shift [9,10]: most notably, electron vacuum polarization (eVP) corrections to the two-photon exchange (2γ exchange), as well as a threephoton-exchange (3γ-exchange) contribution. These are certainly interesting improvements, motivating a more systematic account of nuclear effects in µD and light muonic atoms in general.
In this work, we give a reassessment of the finite-size and polarizability effects in µD using the pionless effective field theory (/ πEFT) of nuclear forces [11][12][13][14][15][16][17] at next-to-next-tonext-to-leading order (N3LO), including effects generated by the structure of individual nucleons. The / πEFT framework gives a better control of theoretical uncertainties in a welldefined perturbation theory, with the small parameter given by P/m π ∼ γ/m π 1/3, where the typical momentum scale P in the deuteron is characterized by the binding momentum γ = √ M N B 45 MeV, with M N the nucleon mass, B the deuteron binding energy, and m π 139 MeV the pion mass. To quantify the uncertainty due to omitted higher-order terms in the / πEFT expansion, we use the methods of Bayesian statistics along the lines of Refs. [18,19]. The momentum scale probed in µD is of the order of α m 1 MeV (with α the fine-structure constant and m the muon mass), well below the limiting scale of / πEFT, set by m π . The expansion in the atomic momentum appears naturally as part of the usual expansion in α.
Our present work is based on the recent N3LO evaluation of the deuteron charge form factor (FF) and the forward longitudinal amplitude of doubly-virtual Compton scattering (VVCS) [20], from which we derive all of the nuclear effects at O(α 5 ), as well as their interference with the eVP. At N3LO, all of the nuclear contributions are governed by only one free parameter -the low-energy constant l C0 S 1 , which describes the coupling of the two-nucleon system to a Coulomb photon, and will be determined here very precisely from the isotope shift. The nuclear corrections to the µD Lamb shift are then a free-parameter-free prediction of N3LO / πEFT. Already at O(α 5 ), we find an appreciable difference with previous evaluations of the "elastic" 2γ contribution to the µD Lamb shift. We then combine our results with the known O(α 6 log α) and O(α 6 ) corrections, coming from the Coulomb distortion [21] and the 3γ-exchange [9], to obtain the full nuclear-structure contribution to the µD Lamb shift. It results in a slightly different extraction of the deuteron radius from the experimental µD Lamb shift, see Fig. 1. This extraction is now indeed consistent with the extraction from arXiv:2206.14066v2 [nucl-th] 28 Oct 2022 the isotope shift, but the latter is, of course, more precise. In what follows, we give more details of this analysis, with emphasis on novel aspects of the 2γ contribution in / πEFT. The "elastic" and "inelastic" 2γ contributions are evaluated in Sections II and III, respectively. Hadronic contributions beyond N3LO are discussed in Section IV. The total results are compiled in Section V. A summary is given in Section VI.

II. CHARGE RADIUS AND ELASTIC 2γ EXCHANGE
We start with the finite-size contributions to the µD Lamb shift at N3LO in / πEFT. Note that at this order / πEFT gives very simple results, many of which can be computed analytically.
The main nuclear-structure effect in the Lamb shift is the finite-size contribution of O(α 4 ): with n the principal quantum number and m r = mM d /(m + M d ) the atomic reduced mass. The squared charge radius, in the case of the deuteron, is given by the slope of the charge FF: . The N3LO / πEFT result for r 2 d reads, order-by-order: where Z = 1.67893 (30) is the residue of the NN scattering amplitude at the deuteron pole [17], and r 2 0 = 1 /2 r 2 p + 3 /4 M −2 p + r 2 n is the isoscalar nucleon charge radius, with the proton Darwin-Foldy term 3 /8 M −2 p added to it. In Ref. [20], l C0 S 1 was chosen to reproduce the deuteron charge radius from µD spectroscopy [4], r d = 2.12562(78) fm, resulting in l C0 S 1 = −2.32(41) × 10 −3 , where the uncertainty in the brackets stems from the error of the deuteron radius and the uncertainty of Z. However, the extraction of r 2 d from µD spectroscopy depends on the theory result for the 2γ-exchange correction (even though the contribution of l C0 S 1 to the 2γexchange correction is small). To avoid this interdependence, it is reasonable to use the deuteron charge radius from the H-D 1S − 2S isotope shift instead. The isotope shift also has a 2γ-exchange contribution, but it has a much smaller correlation with r 2 d , and the contribution of l C0 S 1 to the isotope shift can safely be neglected at the current level of precision. The value of r d obtained from the H-D isotope shift, using the µH value of r p , is [22]: and the corresponding result for the electric contact term coupling is: This value is used in the calculation of the 2γ-exchange correction presented below.
The nuclear effects of O(α 5 ) can be described by the forward 2γ-exchange, see The forward O(α 5 ) 2γ-exchange correction to the energy of a nS state in µD is expressed through the forward VVCS off an unpolarised deuteron, which is parametrised by two scalar amplitudes f L (ν, Q 2 ) and f T (ν, Q 2 ) -the longitudinal and transverse amplitudes, functions of the photon virtuality Q 2 and the lab-frame energy ν [23]. This 2γ-exchange correction reads [24]: where is the (Coulomb) wave function of the nS atomic state at the origin, a = 1/(Zαm r ) is the Bohr radius, and Z is the nuclear charge (equal to 1 for the deuteron). The / πEFT counting assigns Q = O(P), ν = O(P 2 ); the leading terms of f L and f T are, respectively, O(P −2 ) and O(P 0 ) [20]. The factor ν 2 /Q 2 = O(P 2 ) in Eq. (5) further suppresses the transverse contribution, which is N4LO relative to the leading longitudinal one. The N3LO result for f L (ν, Q 2 ) obtained in Ref. [20] thus allows us to calculate the forward 2γ-exchange correction up to N3LO in the / πEFT counting. The elastic contribution to the 2γ exchange is given in terms of the elastic electromagnetic FFs [24], and, neglecting the magnetic and quadrupole contributions (they both arise at higher orders in the / πEFT counting, and their smallness, < 1 µeV, is confirmed using empirical deuteron FF parametrisations [25,26]), reads: with the auxiliary functionsγ 2 yields, order-by-order: with the uncertainty dominated by the higher-order terms in the / πEFT expansion, and estimated following the Bayesian approach in Refs. [18,19]. This significantly deviates from the result obtained in Ref. [24], E elastic 2S = −0.417(2) meV, which used the t 20 deuteron FF parametrisation of Abbott et al. [25], also adopted in the recent Ref. [27]. At the same time, the / πEFT result is confirmed using the recent chiral effective theory (χET) fit for the deuteron charge FF [28,29], with which we obtained E elastic 2S = −0.4456(18) meV, with the uncertainty corresponding to the χET uncertainty of the charge FF. These results, together with the results for the inelastic contribution, are summarised in Table I.
The fact that the two low-energy effective theories give coinciding results vindicates their choice as tools to investigate the low-momenta properties of the deuteron. It also appears that the discrepancy in E elastic 2S between the EFTs and the Abbott parametrisation [25] is due to the poor low-Q properties of the latter. To elucidate this issue, we consider the Friar radius r Fd , given in terms of the charge FF as: This integral coincides, up to a constant factor, with the leading term in the expansion of E elastic 2S in powers of m [30]. Up to N3LO in / πEFT, it can be calculated analytically, the N3LO result being: The dependence of both r 2 d , Eq. (2), and r 3 Fd , Eq. (9), on l C0 S 1 can be represented as a correlation line, which is shown in Fig. 4, with the band corresponding to the estimated N3LO relative uncertainty of 1% (γ/m π ) 4 . One can see that the point corresponding to the parametrisation of Abbott et al. [25] deviates from the correlation line, while the values calculated with the charge FF resulting from / πEFT, the χET fit of [28,29], and the empirical parametrisation of [26] cluster very closely together. Note that it is not only the value of the deuteron charge radius in the parametrisation of [25], r d = 2.094 (9), that tends to be lower than the other results, it is also the higher derivatives of the charge FF being smaller, resulting in a smaller r 3 Fd . The deuteron Friar radius as an integral low-Q feature, and, in particular, its correlation with the deuteron charge radius, can be an important criterion for the quality of a FF parametrisation. πEFT results, with the band showing the estimated 1% N3LO uncertainty. The red disc, purple cross, green diamond, and blue square show the values obtained, respectively, from / πEFT, the χET form factor [29], the parametrisation of Ref. [26], and the parametrisation of Ref. [25]. The dash-dotted lines indicate the isotope shift value of r 2 d and the corresponding value of r 3 Fd obtained at N3LO in / πEFT.

III. INELASTIC 2γ EXCHANGE
The inelastic part of the 2γ-exchange contribution is given in terms of the non-pole parts of the VVCS amplitudes as: which arises from Eq. (5) after a Wick rotation and integration in hyperspherical coordinates. The respective pole parts, as well as the Thomson term entering f T , are assumed to be subtracted from the VVCS amplitudes. The evaluation with the / πEFT VVCS amplitudes is straightforward, resulting in, order-by-order: for the LO+NLO transverse contribution. We include the latter in the final result despite its smallness, both in the / πEFT counting and numeric, for the sake of comparing with other calculations that typically include it. It is also in a very good  agreement with the recent calculations of Refs. [27,31]. We therefore get: The uncertainty here is calculated in the same Bayesian approach as for E elastic

2S
, with the uncertainty of the transverse part neglected. The comparison of the results for E inel 2S is shown in Table I. Our result agrees both with the recent covariant dispersive calculation [27], as well as with the calculation of Ref. [31], within the uncertainties. The latter has a slightly larger in magnitude central value. The value obtained by us is appreciably smaller than the N2LO / πEFT result of Emmons et al. [32], which is not unexpected, given that the nucleonsize effects, included in our calculation but not in Ref. [32], suppress the magnitude of the 2γ-exchange correction. The difference, however, is accommodated by the uncertainty of the N2LO calculation. The data-driven evaluation of Carlson et al. [24] obtains an even larger E inel 2S , albeit with a large uncertainty making it compatible with any of the other results.

IV. HADRONIC CONTRIBUTIONS BEYOND N3LO
In addition to the N3LO / πEFT result, we evaluate the corrections that arise due to the nucleon structure (such as the higher-order terms in the expansion of the nucleon FFs and the nucleon polarisabilities) at orders beyond N3LO. Those corrections are problematic in the strict / πEFT expansion, since they are leading to divergent contributions to E fwd nS and, as such, demand the inclusion of a two-nucleon-two-lepton contact term to regularise this divergence. We circumvent this problem, plugging in the full nucleon FFs in the nucleon-photon vertices at LO, similarly to what is done in χET [27,29]. This allows us to calculate the correction to the inelastic part of the 2γ-exchange contribution due to higherorder terms in the nucleon FFs: The respective effect on the elastic part of the 2γ-exchange contribution is very small, and we neglect it. The (inelastic and subtraction) contributions of the nucleon polarisabilities are calculated in a similar fashion. Namely, the inelastic contribution can be inferred from a dispersive integral with input from empirical deuteron structure functions at energies starting from the pion production threshold [24]; on the other hand, both this term and the nucleon subtraction term can be calculated by re-scaling the sum of the respective proton and neutron contributions with the appropriate wave function factor [21]. For the inelastic contribution, we adopt the dispersive result of [24], Note that this result coincides with that of the re-scaling, using the single-nucleon values from [34], which gives −0.030(2) meV. For the subtraction contribution, we adopt the covariant chiral perturbation theory (χPT) result for the proton subtraction contribution [35] and its neutron counterpart, and use the re-scaling procedure, getting E hadr, subt 2S = 0.009(6) meV, which agrees well with the value adopted in Ref. [ = −0.032(6) meV. (16) Since the hadronic contribution is larger than our N3LO uncertainty estimate for E inel 2S , Eq. (12), we add it to the total. We expect these hadronic corrections to be the only sizeable contributions beyond N3LO, whereas purely nuclear effects at N4LO or higher in the / πEFT expansion should be fully contained by the given uncertainty estimate.

V. TOTAL 2γ EFFECT AND THE DEUTERON RADIUS DETERMINATION
To arrive at our final number for the nuclear-structure correction, we add two further contributions. The first of them is the eVP correction to the forward 2γ exchange, also calculated by us using the / πEFT VVCS amplitudes, with a negligibly small uncertainty. This result agrees with the first evaluation of this contribution by Kalinowski [10]. The second contribution is the off-forward 2γ exchange (or Coulomb distortion). Despite being a subleading effect in the QED expansion, [O(α 6 log α)], it is needed for a meaningful extraction of r d from the µD Lamb shift and for a comparison with empirical results for the nuclear-structure effect. We adopt the value from the theory compilation [21], It is derived using modern nucleon-nucleon potentials; since the deuteron electric dipole polarizability obtained in / πEFT [20] is in agreement with the results obtained with the applied deuteron potentials [36], we expect it to be consistent with / πEFT.
We now compare the theory prediction to the empirical extraction of the 2γ-exchange contribution from the combined µD, µH, and H-D measurements. We start from the theory of the µD Lamb shift compiled in Ref. [22], together with an updated NLO hadronic vacuum polarization contribution [37], see Ref. [38] for details 1 : The first term here contains QED and other structureindependent effects. The second one is the finite-size contribution. The last two are the nuclear-structure effects from, respectively, 2γ and 3γ-exchange. For the latter, we are taking the sum of the elastic [22, ##r3 and r3'] and inelastic contributions [9]. This updated µD theory is then compared with the experimental result from the CREMA collaboration [4], and using the value of r d (µH, iso) from Eq. (3), we obtain the following empirical extraction of the 2γ-exchange contribution: The prediction obtained in the / πEFT framework, Eq. (19), is fully consistent with this updated empirical result, but has an about 3.5 times larger uncertainty.

VI. SUMMARY
An accurate interpretation of the 2γ-exchange correction to the µD Lamb shift is now possible, thanks to the unprecedented experimental precision achieved by the CREMA collaboration [4]. Here we employ a nuclear EFT which, among other good features, allows one to systematically improve the theoretical description of nuclear effects and to quantify the associated uncertainty.
We report on a calculation of the 2γ-exchange contribution to the Lamb shift in µD, performed at N3LO in / πEFT. Based on our earlier calculation of the deuteron VVCS amplitudes [20], we calculate the elastic and inelastic 2γexchange contributions at O(α 5 ). The resulting elastic contribution shows a significant discrepancy with the data-driven result [24], which uses the empirical deuteron form factors of Abbott et al. [25]. We argue that the latter parametrisation is not suitable for the low-Q regime; the correlation between the deuteron charge and Friar radii can be used to judge the quality of the low-Q description. The resulting inelastic contribution, on the other hand, is consistent with the recent evaluations [27,31]. In addition, we evaluate the eVP corrections to the 2γ-exchange contribution, for the first time in the EFT framework.
The systematic nature of the / πEFT expansion allows us to apply a Bayesian procedure to quantify the theoretical uncer-tainty along the lines of Refs. [18,19] (in this context see also Ref. [39], which applies Bayesian reasoning to study the uncertainties arising from a different expansion in muonic atoms and ions).
We also evaluate the higher-order contributions coming from the single-nucleon structure, which are expected to be the dominant terms beyond N3LO. This is achieved by departing from a strict / πEFT formalism; we believe that purely nuclear higher-order contributions play only a very minor rôle compared to these hadronic corrections.
We have updated the original empirical extraction of the 2γexchange contribution [4], resulting in a somewhat improved uncertainty. Our N3LO / πEFT prediction is fully consistent with it, but for now has a larger uncertainty. The same, in turn, is valid for the two alternative extractions of the deuteron charge radius, i.e., the value obtained from combining the µH Lamb shift with the H-D isotope shift [22], and the one obtained from the µD Lamb shift using our / πEFT calculation, Eq. (23). Thus, the µH and µD extractions are, via the H-D isotope shift, consistent with each other; the original inconsistency was mainly caused by the missing eVP contribution, as already noted in [10].
As for an outlook, we note that, at the current level of precision, the single-nucleon effects are becoming appreciable, despite being a few orders of magnitude smaller than the nuclear effects. The nucleon structure is even more important in more tightly bound nuclei, such as helium. This emphasises the importance of investigating the nucleon structure, such as the nucleon polarisabilities. Furthermore, treating the nuclear and nucleon structure consistently, within the same EFT framework, would be an important step to an improved theoretical description.