Improved estimates of the nuclear structure corrections in $\mu$D

We calculate the nuclear structure corrections to the Lamb shift in muonic deuterium by using state-of-the-art nucleon-nucleon potentials derived from chiral effective field theory. Our calculations complement previous theoretical work obtained from phenomenological potentials and the zero range approximation. The study of the chiral convergence order-by-order and the dependence on cutoff variations allows us to improve the estimates on the nuclear structure corrections and the theoretical uncertainty coming from nuclear potentials. This will enter the determination of the nuclear radius from ongoing muonic deuterium experiments at PSI.


Introduction
The root-mean-square charge radius of the proton was recently determined by spectroscopic measurements of the 2S-2P atomic shift, i.e., the Lamb shift (LS) [1], in muonic hydrogen [2,3], where the proton is orbited by a muon instead of an electron as in ordinary hydrogen. With respect to the CODATA-2010 compilation [4], which is based on the combined electron proton scattering data and the spectroscopic measurements in the ordinary hydrogen atom, the accuracy was improved tenfold and a proton radius value smaller by 7 standard deviations was observed. This large deviation between the muonic and the electronic measurements constitutes the so-called "proton radius puzzle". It has attracted a lot of attention since 2010, from both theoretical and experimental viewpoints. Several beyondthe-standard-model theories, including lepton universality violations, have attempted to solve this puzzle (see e.g. [5] for a review). For example, the authors of Refs. [6][7][8] investigated the possibility of the existence of new interaction mediators that can explain not only the proton radius puzzle, but also the (g − 2) muon anomaly. As yet, none of these theories have been either verified or ruled out by experiments. Alternative explanations are being sought after either through novel aspects of hadronic structure [9][10][11], or from renewed analyses of the electron scattering data, e.g., Refs. [12,13]. To date, no commonly accepted explanation of the puzzle has been found. Various new dedicated experiments have been planned to measure electron [14,15] and muon [16] scattering on the proton. In addition, experimental reexamination of ordinary hydrogen spectroscopy is under way, e.g., Ref. [17]. A complementary experimental program based on high-precision spectroscopic measurements on various muonic atoms aims to study the systematics of the discrepancy with ordinary atoms as a function of the atomic mass A and charge number Z. In particular, the CREMA collaboration [18] plans to measure the Lamb shift and isotope shifts in several light muonic atoms. The deuteron is the lightest compound nucleus, made up of one proton and one neutron, and it plays an important role in few-body nuclear physics. The Lamb shift of its muonic atom, µD, is currently being measured at PSI. This measurement will provide a solid and independent test of the systematic uncertainties in the µH experiment. Furthermore, assuming a new interaction mediator that violates lepton universality, the µD experiment may help to constrain the possible couplings of this new interaction to the proton and the neutron. Therefore, it is important to compare the deuteron charge radius r 2 ch 1/2 d extracted from the µD Lamb shift with the values determined from previous and ongoing experiments on eD scattering [19,20], as well as from the precision measurements on the H/D isotope shift [21,22].
The extraction of the nuclear charge radius from LS measurements relies heavily on theoretical input. For the deuteron, the 2S-2P energy transition is related to r 2 ch 1/2 d by [23] where α is the fine-structure constant and m r is the reduced muon mass. Quantum electrodynamic (QED) corrections δ QED , as well as nuclear structure corrections δ pol and δ Zem , are obtained from theoretical calculations. δ QED , which originates from vacuum polarization, lepton self-energy, and relativistic recoil effects, is known with very good accuracy [24]. The uncertainty of the extracted radius is by far limited by the uncertainty of the nuclear structure corrections. These corrections arise from the two-photon exchange (TPE) (see e.g. Fig.1 in Ref. [25]), in which the virtual photons transfer energy and momentum to the nucleus. The TPE contribution is traditionally separated into the sum of the elastic Zemach term δ Zem , proportional to the third Zemach moment [23], and the inelastic polarization term δ pol , i.e., δ TPE = δ Zem + δ pol . The inelastic part is further separated into δ pol = δ A pol + δ N pol , the sum of the nuclear polarization δ A pol and the intrinsic nucleon polarizability δ N pol . For µD, these nuclear structure corrections have been most recently estimated by Carlson et al. [26] using forward dispersion relations to analyze experimental elastic deuteron form factors and inelastic electromagnetic deuteron scattering data. Due to lack of data in the most relevant low-energy quasielastic regime, this analysis suffers from a 35% uncertainty. Theoretical estimates of nuclear corrections were recently made by Friar [27], utilizing the zero-range model. The uncertainty was roughly estimated to be 1-2% due to missing higher-order corrections, mainly from two-body operators and S-D mixing in the deuteron ground state [28]. However, this uncertainty is obtained by dimensional analysis, and is thus only an approximation. A pioneering calculation of δ TPE in µD was made by Pachucki [29], who used a high precision nucleon-nucleon (NN) force of phenomenological nature, namely the AV18 potential [30]. Pachucki quoted ∼ 1% uncertainty on δ TPE , mostly due to higher-order atomic-physics terms in the α expansion, but did not include the uncertainty from nuclear potentials.
In order to obtain an improved estimate, one may rely on microscopic calculations based on nuclear Hamiltonians constructed with state-of-the-art NN potentials [30][31][32]. The NN force is an effective potential emerging from the low-energy quantum chromodynamics (QCD). As such it depends on the adopted resolution scale, a set of fitting parameters, and truncation orders. Realistic NN potentials are typically fitted to reproduce NN scattering data (mostly np and pp) with high accuracy χ 2 ≈ 1 up to pion production threshold. In order to provide a rigorous estimate of the theoretical uncertainty one should consider the effect of varying all the fitting parameters within a certain confidence interval. Such a task is at the moment out of hand. Alternatively, one could explore various realistic NN potentials as a sample of the possible model space.
An early attempt to estimate the theoretical error of the nuclear polarization in µD using several of the NN forces available at that time was made by Leidemann and Rosenfelder [33], who found a potential dependence less than 2%. However contributions from Coulomb distortion were missing in that calculation.
The purpose of this Letter is to provide new calculations and improved estimates of the theoretical uncertainty in the nuclear structure corrections using state-of-the-art nuclear potentials from chiral effective field theory (χEFT) [31,32]. The χEFT nuclear potentials result from a systematic expansion of the interaction in powers of a soft momentum scale over a hard momentum scale Q/Λ χ [31]. The soft scale Q is associated with the typical momenta of the system or dynamics under consideration, and the hard scale Λ χ determines at which scale the effective theory breaks down. The short range part of the χEFT potential stems from contact terms in the Lagrangian, which are controlled by free parameters in the χEFT Lagrangian, commonly known as low energy constants (LECs). Their values are set by fitting the interaction to the NN scattering data at min-imum χ 2 /datum for given resolution scales, determined by the regularization cutoff Λ ∼ O(Λ χ ). For the chiral forces better χ 2 /datum and smaller Λ dependence in predictions is obtained with the increase of the expansion order. One can improve the description of the np scattering data, going from a χ 2 /datum of 36.2 at next-to-leading order (NLO) calculation up to 1.1 at next-to-next-to-next-to-leading order (N 3 LO) in the energy range of 0-290 MeV [34]. By performing a systematic study in χEFT order-by-order and by varying the resolution scale Λ, we can better assess the theoretical error for δ TPE stemming from nuclear forces.

Calculation details
To calculate the deuteron polarization we essentially follow the derivation provided in our recent paper [25]. Within the multipole expansion formalism we take into account the leading electric dipole term δ (0) D1 , the Zemach related contribution δ (1) Z3 , and other multipole corrections including a monopole δ (2) R2 , a quadrupole δ (2) Q and an interference between two rank-1 operators δ (2) D1D3 . The small parameter in the expansion is η = √ m r /m d where m r is the reduced muon mass and m d is the deuteron mass. According to our power counting, the contribution of a term with superscript (1) (or (2)) to δ A pol is suppressed by η (or η 2 ) relative to that of the leading δ (0) D1 term. We also include the Coulomb distortion corrections δ (0) C , relativistic longitudinal and transverse corrections δ (0) L and δ (0) T , magnetic dipole corrections δ (0) M , and finite-nucleon-size corrections including δ (1) Z1 , δ (1) R1np and δ (2) NS . Therefore, we calculate δ A pol as (2) Detailed formulas relating (most of) these corrections are found in [25] and are not repeated here. Except for δ (1) Z3 , δ (1) Z1 and δ (1) R1np , which are ground-state observables, each of the above contributions can be written as a sum of terms of the form where g a (ω) is an energy-dependent weight function (different for each of the corrections), ω th is the threshold excitation energy, and S a (ω) is the response function. S a (ω) is given by whereÔ a is the transition operator (that can be different in the above corrections), |0 and | f are the ground and final eigenstates of the Hamiltonian H with eigenvalues E 0 and E f respectively. The integration in Eq. (4) carries over the excited scattering states of the deuteron. The Coulomb corrections δ (0) C used here contain only the term of order α 6 ln α. Therefore we have where S D1 (ω) is the electric-dipole response function with operator D 1 ≡ R p Y 1 (R p ) [25]. The terms of higher orders in α, which were included in δ (0) C of Ref. [25], only correct δ A pol by 0.25%, and are thus neglected in this analysis. Because the deuteron is a nucleus with the total angular momentum J 0 = 1 in the ground state, the magnetic term δ (0) M , which is negligible in the µ 4 He + , has to be included here. It relates to the magnetic response function S M1 (ω) by where g p = 5.586, g n = −3.826, and m p is the proton mass. The magnetic dipole (M1) operator is defined by s p − s n [29]. The sum of the terms δ (1) Z3 and δ (1) Z1 cancels exactly the elastic Zemach term δ Zem : where β = 4.120 fm −1 and λ = 0.01935 fm 2 are used as in Refs. [25,27] to reproduce the proton charge radius r 2 ch 1/2 p = 0.8409 fm [3] and the neutron charge radius squared r 2 ch n = −0.1161 fm 2 [35]. The n-th moment r n (2) is defined by with ρ p denoting the point-proton density in deuteron. The cancellation of δ Zem was explicitly applied in Refs. [27,29], whose results for nuclear structure effects are provided as the combination of δ A pol and δ Zem . In this work, we will also provide our final result as a combination of the two.
Unlike the µ 4 He + case of Ref. [25], a pp charge correlation is not present in µD, as it contains only one proton. On the other hand, δ (1) R1np , which denotes a contribution from the np charge correlation, has to be included, and is calculated as where the operator |R p − R n | denotes the relative distance between the proton and the neutron.
To calculate |0 and | f we expand them on the harmonic oscillator (HO) basis [36] and then diagonalize the Hamiltonian matrix, whose dimension is set by N max = 2n + ℓ, where n is the HO quantum number and ℓ is the relative angular momentum between the proton and the neutron. We use the HO basis because it is very flexible and allows the use of both local potentials in coordinate space, like the AV18, as well as non-local forces, like those derived in χEFT, within the same framework. From the diagonalization of the Hamiltonian matrix we obtain a set of eigenstates |µ and eigenvalues E µ . Once the size of the basis is large enough, the lowest eigenvalue state coincides with the ground-state of the deuteron. All the other states |µ can be regarded as a discretization of the two-body continuum. In terms of this discrete basis, Eq. (3) becomes a sum over the transition probabilities from the ground state to the discretized excited states with weighted functions g a as where ω µ = E µ − E 0 is the excitation energy, and the finite sum runs over all the calculated states. As proven in Ref. [37], even if one is approximating continuum states with discrete states, Eq. (10) becomes exact when the size of the basis becomes large enough (with increasing N max ). Furthermore, for the two-body problem the basis size remains tractable so that a direct diagonalization is always possible and one does not need to introduce the Lanczos algorithm as in Ref. [37]. We calculate the deuteron ground-state properties and polarizability using the AV18 potential [30] and NN forces from χEFT. The AV18 results serve as checks with Ref. [29] and the χEFT results are genuine new. We use chiral potentials developed by Epelbaum et al. [38] at different order in χEFT, starting from NLO to N 3 LO. For these, a certain range of cutoff variation will be explored, going from Λ = 400 to 600 MeV for the cutoff to regularize the Lippmann-Schwinger equation [39] and fromΛ = 600 to 700 MeV for the spectral function cutoff (see [38] for details). Hereafter, we shall refer to these potentials either individually as N k LO(Λ,Λ) or collectively as N k LO-EGM. We will also use the chiral potential at N 3 LO developed by Entem and Machleidt at a fixed cutoff Λ = 500 MeV [34] (N 3 LO-EM). By doing so, we will provide a systematic study of the convergence and control the theoretical uncertainty.

Results
To estimate the possible spread in δ TPE due to the variation in nuclear Hamiltonian models, we first present few deuteron ground-state observables and follow their evolution with the nuclear forces. Our numerical results for the energy E 0 , the rootmean-square structure radius r 2 str 1/2 d , the quadrupole moment Q d , the electric polarizability α E , and the magnetic susceptibility β M are presented in Table 1. The first three are also presented graphically in Fig. 1. With the AV18 and N 3 LO-EM we reproduce the experimental binding energy, 2.224573(2) MeV [40], up to 5 decimal digit as in Refs. [30] and [34]. For the chiral N k LO(Λ,Λ) potentials we use the cutoff values provided by Epelbaum [41]. We note that the range of the explored cutoffs does not coincide exactly with the cutoff range spanned in Ref. [38]. Only for N 2 LO there is an exact one-to-one correspondence of cutoff sets used; there we agree with Epelbaum's values within the third decimal digit for energy, radius and quadrupole moment.
The electromagnetic observables in Table 1 (all but E 0 ) have been calculated assuming that nucleons are point-like and using only one-body operators with no relativistic corrections (RC) or meson-exchange currents (MEC). The nonrelativistic onebody operators correspond to the dominant contributions to the electromagnetic observables in the χEFT expansion.
The charge radius is related to the structure radius by finite nucleon size corrections as where the structure radius is separated into the leading component from the nonrelativistic one-body current and the subleading part from RC+MEC. In the χEFT language, RC and MEC (also called two-body currents) enter at higher order in the Q/Λ χ expansion with respect to the leading component, specifically at N 2 LO and N 3 LO, respectively, see e.g. [42]. Similar corrections enter in other electromagnetic observables as well.
Relativistic and two-body corrections should be evaluated consistently within a theory as to satisfy gauge invariance. They have already been derived in χEFT, see e.g. [43][44][45], where pion-exchange and contact two-nucleon currents appear, but are not applied in this work. Instead, following Refs. [30,34], we show in Table 2   Comparing our α E calculated with the AV18 potential with the result of [49] we get an agreement of about 0.15%. Our N 3 LO-EGM and N 3 LO-EM results agree within 1% among the AV18 result and the pion-less EFT results of [50]. For the magnetic susceptibility using the AV18 potential we obtain β M = 0.0679 fm 3 , which agrees within 0.1% with previous calculations using the same potential and within 1% with other 1 In Ref. [  calculations using different potentials [51]. The latter is comparable with the sensitivity to the NN potential that we observed in Table 1. In fact, our N 3 LO-EGM and N 3 LO-EM results for β M scatter at ∼ 1% and are within few percent of the results based on pion-less EFT from Ref. [50]. In Fig. 1, we first show the binding energy. One observes that as the chiral order gets higher, convergence is approached and the band due to the cutoff variation becomes order-by-order smaller. One also notes, that despite the fact that the deuteron binding energy is a prediction of N k LO-EGM potentials [38], it agrees well with the values from the AV18 and N 3 LO-EM potentials, where the deuteron binding energy has been included in the fit [30,34]. Overall, they deviate from the experimental value by only 0.4%.
For the structure radius the situation is different. While quite a nice convergence in the chiral order is achieved, the band obtained by the cutoff variation becomes larger at N 3 LO. Interpreting the cutoff band as an estimate of the theoretical error, it means the error is larger for the highest chiral order, which also results in a larger deviation from the experimental value. This was already observed in Ref. [38] and attributed to the effect of two-body currents. For r 2 str 1/2 d and Q d we find that the values obtained with the chiral potential N 3 LO-EM and AV18 fall marginally out of the cutoff band spanned by the N 3 LO-EGM potentials. This could imply that a slightly different parameterization of the two-body currents is needed for each potential. If we introduced the two-body operators, they would come with new low-energy constants (LECs). The new LECs would have to be calibrated in a way to compensate the cutoff dependence, thus reducing the theoretical errors.
A full compilation of the nuclear structure corrections is displayed in Table 3. We first discuss the calculation with the AV18 potential, comparing our results with Pachucki [29]. As already pointed out in [25] we find very good agreement with Pachucki for the leading order dipole and Coulomb correction terms. Small but non-negligible differences appear in the relativistic corrections. These differences can be traced to the leading-order truncation in the ω/m r expansion of the weight function g a made in [29], whereas we have used the full form given in [25]. Our results for δ (2) R2 , δ (2) Q and δ (2) D1D3 are smaller than those obtained by Pachucki in [29] by ∼ 8%, although the same formulas are used. As an independent check, we use the deuteron electric-quadrupole response function calculated by Arenhövel [52,53] from the Paris potential and obtain δ Q = 0.060 meV, which is only 1% smaller than our corresponding value. 2 This 1% discrepancy is also consistent with the difference of deuteron structure radius between the two calculations and compatible with the sensitivity to different potentials. Our calculated magnetic dipole contribution δ (0) M is ∼ 50% smaller than that in [29], despite the fact that the same formula, i.e., Eq. (6), is used. As noted above, our values for the related magnetic susceptibility β M stand in line with previous calculations. We have also integrated the deuteron M1 response function from Arenhövel [52,53] and obtained δ (0) M = 0.0067 meV. The evolution with EFT orders of the total nuclear elastic/inelastic contribution δ A pol + δ Zem , the dipole component δ (0) D1 , and the magnetization term δ (0) M are presented in Fig. 2, using the same values of the cutoffs (Λ,Λ) as in Table 1. We find similar convergence pattern and cutoff dependence for δ A pol + δ Zem and δ (0) D1 . This is a reflection of the dipole dominance in δ A pol . We also note that the cutoff bands increase with the chiral order, similarly to what we observe for r 2 str 1/2 d and Q d . In fact, we found 2 The difference is calculated with δ Q in more digits than is presented in Table 3. strong correlation between r 2 str 1/2 d and δ (0) D1 (see also [54]). For the magnetic correction δ (0) M , we observe that the cutoff dependence is large starting from the NLO and it is of comparable size as at the N 3 LO. This is again due to the effect of the missing two-body currents, which appear already at NLO in the case of magnetic transition [42]. However, given the size of this term, such further refinements are negligible here.  Comparing the spread of each term in Table 3 due to the potentials, we find that the variation of δ A pol with potential models mostly comes from the variation in the dipole term δ (0) D1 . In fact the contribution of all other terms to the spread in δ A pol is about 1µeV. Averaging over all potentials we can estimate δ A pol as and δ Zem as where the standard deviation represents an estimate of the uncertainty in the high precision nuclear Hamiltonians at a 1σ level and can be used to estimate confidence intervals for the polarization and the Zemach moment. From N 2 LO to N 3 LO in the χEFT expansion, the value of δ Zem +δ A pol changes by 0.3%. This can be considered the systematic error due to the χEFT truncation and needs to be included in the total error budget.
The atomic physics error from further corrections of order (Zα) 6 was estimated by Pachucki to be 1% [29].
Combining all the uncertainties in a quadrature sum, the overall accuracy of δ Zem + δ A pol is expected to be about 1.16%. 5 The proton two-photon exchange contribution to the µD Lamb shift is estimated to be −0.043(3) meV [29]. However, Ref. [26] suggests that only the inelastic part of such contribution is related to δ N pol . Including also the neutron effect, the overall δ N pol is estimated as in Ref. [26] to be −0.027(2) meV 3 . Therefore, we provide a total nuclear/hadron two-photon exchange contributions to µD Lamb shift as δ TPE = δ A pol + δ Zem + δ N pol = −1.690 ± 0.020 meV (14) The elastic contribution δ Zem depends on the input of the proton radius. However, due to cancellation of δ Zem in the elastic and inelastic parts of δ TPE , the only proton-radius dependence in δ TPE comes from the δ (2) NS term in δ A pol . Using the CODATA value 0.8775 fm [4], instead of the proton radius 0.8409 fm from the µH experiment [3], will change δ Zem by −0.007 meV, but δ TPE by only −0.0016 meV. This is only 0.1% of δ TPE , thus negligible in the quadrature sum of the error.

Conclusions
A solid understanding of the theoretical indetermination is crucial to constrain any explanation of the proton radius puzzle based on physics beyond standard model. In this work we have made an attempt to constrain the nuclear structure corrections to the µD Lamb shift utilizing state-of-the-art nuclear potentials derived from chiral EFT. Considering the contributions to the nuclear polarization up to order O(η 3 ), where η = m µ /m d , we obtain an estimate for δ TPE = −1.690 meV. Combining the spread in the results due to the nuclear potentials (0.5%) and the convergence with χEFT orders (0.3%), we evaluate the nuclear physics error to be about 0.6%. This should be added in quadrature to the ∼ 1% error coming from atomic physics.
Our predicted δ TPE differs from that obtained by Pachucki [29] by only 0.6%. However, this excellent agreement can be regarded as partly accidental, since differences of several individual contributions are larger than 0.6%.
Our improved estimates will benefit the determination of r ch 2 1/2 d from the ongoing µD Lamb shift measurement.