Nuclear structure corrections to the Lamb shift in $\mu^3$He$^+$ and $\mu^3$H

Measuring the 2S-2P Lamb shift in a hydrogen-like muonic atom allows one to extract its nuclear charge radius with a high precision that is limited by the uncertainty in the nuclear structure corrections. The charge radius of the proton thus extracted was found to be 7-sigma away from the CODATA value, in what has become the yet unsolved"proton radius puzzle". Further experiments currently aim at the isotopes of hydrogen and helium: the precise extraction of their radii may provide a hint at the solution of the puzzle. We present the first ab initio calculation of nuclear structure corrections, including the nuclear polarization correction, to the 2S-2P transition in $\mu^3$He$^+$ and $\mu^3$H, and assess solid theoretical error bars. Our predictions reduce the uncertainty in the nuclear structure corrections to the level of a few percents and will be instrumental to the on-going $\mu^3$He$^+$ experiment. We also support the mirror $\mu\,^3$H system as a candidate for further probing of the nucleon polarizabilities and shedding more light on the puzzle.


I. INTRODUCTION
The root-mean-square (RMS) charge radius of the proton r p ≡ r 2 p was recently determined with unprecedented precision from laser spectroscopy measurements of 2S-2P transitions in muonic hydrogen µH, where the electron is replaced by a muon [1,2]. The extracted r p differs by 7σ from the CODATA value [3], which is based in turn on many measurements involving electron-proton interactions. This discrepancy between the 'muonic' and 'electronic' proton radii (r p (µ − ) and r p (e − ), respectively) is known as the "proton radius puzzle," and has attracted much attention (see, e.g., Ref. [4] for an extensive review and Ref. [5] for a brief summary of current results and ongoing experimental effort). In an attempt to solve the puzzle, extractions of r p (e − ) from the ample electronproton (ep) scattering data have been reanalyzed by, e.g., Refs. [6][7][8][9], while several planned experiments aim to remeasure ep scattering in new kinematic regions relevant for the puzzle [10,11]. r p extracted from electronic hydrogen is also being reexamined, both theoretically [12] and experimentaly [13][14][15], as well as the Rydberg constant [15,16], which is relevant for several radius extraction methods. A few of the theoretical attempts to account for the discrepancy between r p (e − ) and r p (µ − ) include new interactions that violate lepton universality [17][18][19] and novel proton structures [20][21][22][23][24]. Yet the puzzle has not been solved. Answers may be provided (see, e.g. Refs. [25,26]) by a planned experiment at PSI [27] to scatter electrons, muons, and their antiparticles off the proton using the same experimental setup.
Alternatively, it will be insightful to study whether the puzzle also exists in other light nuclei, and whether it depends on the atomic mass A, charge number Z, or the number of neutrons N . In particular, the CREMA collaboration plans to extract high-precision charge radii from Lamb shift measurements that were recently performed in several hydrogen-like muonic systems [5,28], namely, µD, µ 3 He + , and µ 4 He + . These measurements may unveil a dependence of the discrepancy on the isospin of the measured nucleus and, in particular, probe whether the neutron exhibits a similar effect as the puzzling proton. To obtain some control over these issues, it is advisable that nuclei with different N/Z ratios will be mapped out. It is the purpose of this Letter to perform an ab initio calculation of nuclear structure corrections (including nuclear polarization), with solid error estimates, for the µ 3 He + system and for its nuclear mirror, µ 3 H. The Lamb shift [29] is the 2S-2P energy difference consisting of where, in decreasing order of magnitude, the three terms include: quantum electro-dynamics (QED) contributions from vacuum polarization, lepton self-energy, and relativistic recoil in δ QED ; finite-nucleus-size contributions in δ FS (R c ), where R c ≡ R 2 c is the nuclear RMS charge radius; and contributions from two-photon exchange (TPE) between the lepton and the nucleus in δ TPE . The last term can be divided into the elastic Zemach term and the inelastic polarization term, i.e., δ TPE = δ Zem + δ pol . Additionally, each of these terms is separated into contributions from nuclear (δ A ) and nucleonic (δ N ) degrees of freedom, δ TPE = δ A Zem + δ N Zem + δ A pol + δ N pol . In light muonic atoms, δ QED ≈ 10 2 -10 3 meV and is estimated from theory with a precision better than 10 −3 meV [30][31][32][33].
At leading order δ FS (R c ) = m 3 r (Zα) 4 12 R 2 c , with m r the reduced mass of the muonnucleus system, while higher-order contributions are at the sub-percentage level [30]. The limiting factor for the attainable precision of R c extracted from Eq. (1) is by far the uncertainty in δ TPE . This was confirmed in two recent papers that reviewed the theory in µD [32], and in µ 4 He and µ 3 He [33]. Ref. [32] covers all the theoretical contributions to the Lamb shift in µD, including a summary of recent efforts by several groups [34][35][36][37] to accurately obtain δ TPE in µD and reliably estimate its uncertainty, which comes out an order of magnitude larger than the uncertainties in the other terms. Ref. [33] details all the contributions for the two helium isotopes. Many terms are recalculated, not including the polarization correction δ pol . For µ 4 He + , ab initio nuclear calculations were recently applied in Ref. [38], improving on decades-old estimates of δ pol . For three-body nuclei, the only calculations of δ pol are outdated; based on old and simplistic nuclear models, their results are either inaccurate [39] or imprecise [40], reinforcing the need for modern, accurate, ab initio calculations for the three-body nuclei.

II. NUCLEAR STRUCTURE CONTRIBUTIONS
The nuclear Zemach term δ A Zem enters Eq. (1) as the elastic nuclear-structure contribution to δ A TPE 1 . This term is of order (Zα) 5 and is defined as where r 3 (2) is the 3rd nuclear Zemach moment 2 . Friar & Payne showed [43] that the first-order corrections in δ A pol contain a part that cancels δ A Zem exactly. Calculation of this part can thus be avoided, as was done in Ref. [34], providing only the sum δ A TPE = δ A pol + δ A Zem . However, following Refs. [35,38,44], we calculate explicitly all the parts of δ A pol , including the Zemach term, as detailed below. This is done in order to: (a) allow comparison with other values in the literature, and (b) provide theoretical support for the alternative way of extracting R c from Eq. (1) where the Zemach term is phenomenologically parameterized as [30] 1 δ A Zem was derived by Friar [41] as the first-order Zα correction to δ FS (Rc) and is called the 'Friar' term in Ref. [32]. 2 We refer only to charge-charge Zemach moments; for more details see, e.g., Ref. [42].
As in Refs. [35,38], the energy correction due to nuclear polarization is obtained as a sum of contributions Detailed formulas pertaining to most of the terms in Eq. (4) are found in [38] and are not repeated here. The largest contribution comes from the leading term, δ D1 , related to the electric dipole. To this we add relativistic longitudinal and transverse corrections δ C . Here we follow Ref. [35] and include in δ (0) C only the logarithmically enhanced term from the next order in Zα. We generalize the treatment in Ref. [35] of the magnetic term δ M by using the impulse approximation operator that includes the orbital angular momentum [45]. First-order corrections δ (1) R3 and δ (1) Z3 are related to a proton-proton correlation term and to the 3rd nuclear Zemach moment, respectively. Finally, at the next order we have the monopole δ N S , which we elaborate on below.

III. NUCLEON-SIZE CORRECTIONS
The TPE in the point-nucleon limit is expressed as the interaction of photons with the structureless charged protons, while the neutrons are ignored. In this limit, the point-proton density operator iŝ where τ 3 a is the isospin projection operator. When the finite nucleon sizes are considered,ρ p (R) must be convoluted with the proton's internal charge distribution, and a similar convolution is applied to the point-neutron density operator Following Refs. [38,44], we apply a low-momentum expansion for the nucleon form factors, parameterized here by their mean square charge radii, r 2 n/p ≡ r 2 n/p . We adopt r 2 n = −0.1161(22) fm 2 [46]. For the proton, we may use either r p (e − ) = 0.8775(51) fm [3] or r p (µ − ) = 0.84087(39) fm [2]. In fact, until the "proton radius puzzle" is resolved (or when R c and other properties of the nuclei under consideration are measured using muons), we should use r p (e − ) for comparison with the literature, which is based on data obtained with electrons, and r p (µ − ) for predictions in muonic systems.
The leading NS correction δ (1) N S is the sum of nucleonnucleon correlations in δ (1) R1 and Zemach-like terms in δ The former is expressed as which includes proton-proton (pp) and neutron-proton (np) correlations. It is an NS correction to the pointnucleon contribution δ R3pp in Ref. [38]). For the calculation of Zemach-like terms using point-nucleons we define with i, j denoting either p or n. The 3rd nuclear Zemach moment is thus calculated as where the first term is the point-nucleon limit and the second is the (approximated) NS correction. Accordingly, the point-nucleon Zemach term δ (10) Lastly, the nucleonic TPE correction δ N TPE also enters Eq. (1). We defer the treatment of this hadronic contribution to a dedicated section below.

IV. METHODS
Most of the above contributions can be written as sum rules of several nuclear responses with various energydependent weight functions [35,38]. They were evaluated using the newly developed Lanczos sum rule method [47]. Ground-state observables of 3 He and 3 H, as well as Lanczos coefficients, were obtained using the effective interaction hyperspherical harmonics (EIHH) method [48,49].
As only ingredients we employed in the nuclear Hamiltonian either one of the following state-of-the-art nuclear potentials: (i) the phenomenological AV18/UIX potential with two-nucleon [50] plus three-nucleon [51] forces; and (ii) the chiral effective field theory χEFT potential with two-nucleon [52] plus three-nucleon [53] forces.
It is of utmost importance to have realistic uncertainty estimates for our nuclear TPE predictions. These terms are the least well known in Eq. (1), and their uncertainties determine the attainable precision of R c extracted from Lamb shift measurements. We considered many sources of uncertainty, namely: numerical; nuclear model; isospin symmetry breaking; higherorder nucleon-size corrections; missing relativistic and Coulomb-distortion corrections; higher multipoles, terms of higher-order in Zα; and the effect of meson-exchange currents on the magnetic contribution. Their individual and cumulative effect on δ A pol , δ A Zem , and δ A TPE have been estimated and applied to the results given below. More details about these uncertainty estimates are given in the Supplementary Materials [54].

A. Results
We first compare a few observables we have calculated for the 3 He and 3 H nuclei with corresponding theoretical and experimental values available in the literature. In Table I we present the ground-state binding energy BE, charge radius R c , and magnetic momentμ gs , as well as the electric dipole polarizability α E . In general, good agreement is found with other calculations.
Our results do not include isospin-symmetry breaking (ISB), except for the Coulomb interaction between protons in 3 He. Calculations by other groups shown in Table I usually do not include ISB effects; notable exceptions are Ref. [56], which includes the T = 3/2 isospin channel in the ground-state wave function, and Ref. [55] that provides results either including or excluding it. One observes that including ISB alters BE by a few keV. In addition, the 3 He BE, not used in the calibration of the Hamiltonians, is overestimated at a sub-percentage level, and this is slightly worsened when ISB is included. As discussed in Ref. [62], changes in BE shift the threshold of sum rules, affecting mostly sum rules with inverse energy dependence, such as α E discussed below. For the other observables in Table I, the estimated uncertainty stemming from ISB is 1%.
Charge radii R c shown in Table I where we omit contributions from the spin-orbit radius (negligible for s-shell nuclei) and meson-exchange currents. The last term in Eq. (11) is the Darwin-Foldy term, where m p is the proton mass, taken from Refs. [3,46]. In Table I, we show only R c values obtained using r p (e − ) and experimental values obtained only with electrons. As a direct result of Eq. (11), using r p (µ − ) would decrease R c by 0.016 (0.018) fm for 3 He ( 3 H). We note that the uncertainty, currently governed by nuclearmodel dependence, is slightly larger than the effect of varying r p . It should also be noted that our R p values agree with the hyperspherical harmonics calculations of the Pisa group [56] for both nuclear potentials, suggesting a small ISB effect, while the Monte-Carlo calculations of Ref. [59] show less agreement and hint at a larger ISB effect. Considering that radii were not included in the calibration of the nuclear Hamiltonians, it would be interesting to further investigate their sensitivity to the theoretical apparatus. In particular, work is in progress to include meson-exchange currents [68]. Currently, for 3 He the AV18/UIX charge radius is in better agreement with the experimental value, while for 3 H the experimental error bar is larger than the nuclear-model dependence, and calls for a more precise measurement. Concerning the magnetic moments, our results are comparable with the other impulse approximation calculations presented in Table I, which deviate from experiment due to the absence of meson-exchange currents. However, we do not include meson-exchange currents in δ A TPE , since the contribution of the magnetic term δ (0) M is small enough to make these corrections negligible.
The electric dipole polarizability α E is an inverseenergy-weighted sum rule of the dipole response and is therefore closely related to δ A pol . Our results are in agreement with previous calculations, especially the recent Ref. [62]. As in [38], α E is found to be nuclear-model dependent. We provide first results for the unmeasured α E of 3 H with the AV18/UIX potential, which lies within the uncertainty estimates of [62].
We now turn to the Zemach terms, first listing available values in the literature. In Refs. [30,31] Borie calculated δ A Zem , following Friar [41], using a Gaussian distribution that fits the nuclear-chargeradius obtained from electron experiments. The result 4,5 was δ A As explained above, all of these results should be compared with our calculation that uses r p (e − ) as input and yields δ A  (19)(16) meV, where the first uncertainty results from nuclear-model dependence and the second includes all other sources. Our result is in agreement with these references (based on comments made in Refs. [33,69], we assume that the error-bars in Ref. [33] are not exhaustive). However, for the muonic systems considered here we use r p (µ − ), which gives We note that with the given error-bars this result is also in agreement with Refs. [31,33,69]. The use of Eq. (3) is adopted from Refs. [30,31], where 6 C 3 He = −1.35(4) meV fm −3 . The results of Ref. [69] can also be used to extract C 3 He = δ A Zem /R 3 c = −1.42(4) meV fm −3 from the e− 3 He scattering data. Our calculations of δ A Zem and R c with either value of r p give C 3 He [r p (e − )] = −1.383(05)(20) meV fm −3 and C 3 He [r p (µ − )] = −1.388(05)(21) meV fm −3 , which both agree with Refs. [30,31,69]. Evidently, the nuclear-model dependence is diminished for this value, since it is proportional to the geometrical ratio r 3 (2) /R 3 c . Similarly to R c discussed above, the difference between δ A Zem results obtained with the two nuclear potentials stems from the different point-proton distributions, and this largely cancels out in C, reducing its total relative uncertainty compared to δ A Zem . Repeating the above procedures we obtain predictions for 4 Ref. [31] is the arXiv version of Ref. [30], which has been updated since publication; in particular, δ A Zem 3 He was increased by ∼20% with respect to the published version. 5 The result is given using our sign convention. 6 See footnote 5.  Next, the nuclear polarization correction to the Lamb shift -δ A pol -is obtained by summing up the terms in Eq. (4). Their values for µ 3 He + and µ 3 H, calculated with the two nuclear potentials, are shown 7 in Fig. 1. Here, the NS corrections are obtained using only r p (µ − ). Taking the mean value of the two nuclear potentials we obtain where we retain the use of first and second brackets for uncertainties from nuclear-model dependence and from all other sources, respectively. Our result for µ 3 He + agrees with Rinker's −4.9 ± 1.0 meV obtained fourty years ago [40]. The µ 3 H case was rarely studied. We note, however, that a comparison with the simplistic calculation of Ref. [39] reveals a similar ratio of ∼9 between δ A pol of µ 3 He + and of µ 3 H, both in Ref. [39] and in our work. 7 The numerical values are detailed in the Supplementary Materials [54].

V. HADRONIC TPE
The last ingredient in δ TPE is the contribution from two-photon exchange with the internal degrees of freedom of the nucleons that make up the nucleus, i.e., δ N TPE = δ N Zem + δ N pol . Since it is dictated by the hadronic scale, about 10 times higher than the nuclear interaction, this contribution can be approximated as the sum of TPE effects with each of the individual nucleons. The various terms that contribute to δ N TPE are estimated based on previous studies of µH, as recently done for µD in Ref. [32]. Specifically, as suggested by Birse and McGovern [70], we adopt values of δ Zem and δ pol in µH that are combinations of results from Refs. [21,71], as detailed below.
As Friar showed in Ref. [72], the intrinsic Zemach term of each proton contributes to δ TPE of the nucleus as an additional NS correction, not accounted for in the NS corrections detailed above 8 . We denote this term δ N Zem and find its contribution to be proportional to the analogous term in µH by We take δ Zem (µH) = 0.0247(13) meV 9 and obtain In Ref. [36], δ N pol of µD was extracted from electron scattering data. Here, we resort to estimating δ N pol by relating it to the proton polarization correction in µH via [34,37,75] assuming that the neutron polarization contribution is the same as that of the proton. Here we use δ pol (µH) = 9.3(1.1) µeV 10 . Based on current knowledge of the nucleon polarizabilities [76], we assign an additional 8 In our notations this term appears as an NS correction to δ (1) R3 . 9 We use the same value as in [32]. Here, δ Zem (µH) stands for the elastic + non-pole parts of δ TPE (µH), and not for the nonrelativistic limit that is related to the proton's 3rd Zemach moment (see Refs. [73,74]). 10 δ pol (µH) = δ p inelastic + δ p subtraction . For the former we follow Ref. [36] and adopt 13.5 µeV, which is an average of three values from Ref. [71], and for the latter we use −4.2(1.0) µeV from Ref [21]. 20% uncertainty to the neutron polarization contribution. Another possible error in δ N pol arises from neglecting medium effects and nucleon-nucleon interferences in Eq. (19). These effects can be estimated by comparing the calculated δ N pol (µD) with the result evaluated in Ref. [36] from scattering data. This yields a ∼29% correction. Until this correction is calculated rigorously in other light muonic atoms, we estimate it to be of a similar size, multiplied by A/2, making it the dominant source of uncertainty in our δ N TPE . Eventually, we obtain Summing up the results in Eqs. (18) and (20) we obtain the total contribution from the nucleon degrees of freedom In µ 3 He + , δ N TPE is ∼5% of δ A TPE , i.e., about twice the overall uncertainty in δ A TPE . For µ 3 H we obtained that δ N TPE is ∼9% of δ A TPE , which is more than three times the uncertainty in the latter. Therefore, our precision in predicting δ A TPE can be important not only for the determination of R c from muonic Lamb shift measurements, but also for probing δ N TPE , if these measurements reveal discrepancies with electronic experiments that may indicate exotic contributions to δ N TPE . A study of the Lamb shift in µ 3 H will be especially sensitive to the nucleon polarizabilities, since their relative contribution is much larger in this case.

VI. SUMMARY
We have performed the first ab initio calculation of δ A Zem and δ A pol for both µ 3 He + and µ 3 H, using state-of-the-art nuclear potentials. Many possible sources of uncertainty have been considered, yet the resulting uncertainties of a few percents are much lower than in previous estimates of δ A pol and δ A TPE , which relied on imprecise data and simplistic models. In addition, our δ A Zem calculations agree with previous estimates and with recent analysis of e − 3 He scattering, and provide predictions towards 3 H measurements. They were also adapted for muonic systems by incorporating r p (µ − ) -the proton radius measured with muons.
Ultimately, our results will allow two alternative ways of extracting a much more precise R c from a recent measurement [5,28,77] of the Lamb shift in µ 3 He + , and from an analogous measurement we encourage to conduct in µ 3 H. The precision of the charge radii of 3 He and 3 H could be thus improved by factors of ∼5 and ∼50, respectively, which could have interesting implications for nuclear physics.
Finally, we estimate the hadronic contribution δ N TPE in these systems, and find it to be larger than our uncertainty estimates in δ A TPE . Therefore, this combined theoretical and experimental effort may not only shed some light on the "proton radius puzzle," but could also probe the elusive nucleon polarizabilities tightly connected to it.
Zα expansion: Except for the logarithmically enhanced Coulomb distortion contribution, we include in our calculation of δ A pol all terms of order (Zα) 5 . Since (Zα) is small for these systems, the missing contribution from all the higher-order terms can be approximated by the first unaccounted-for term in the series, (Zα) 6 , which is naturally estimated to be (Zα) 1.46% (0.73%) of the size of δ A pol in 3 He ( 3 H).
Magnetic MEC contribution: The magnetic dipole term δ (0) M is calculated using the impulse approximation (IA) operator, and is therefore missing significant corrections, mainly due to meson-exchange currents (MECs) 13 . The same IA operator was also used to calculate the magnetic momentμ gs in the nuclear ground state. The results, presented in Table I, show a deviation of the IAμ gs calculated with the AV18/UIX (χEFT) potential from the very precise experimental values by 23% (21%) for 3 He and 15% (13%) for 3 H. The same relative errors were therefore assumed also for the small IA value obtained for δ (0) M . The total uncertainties were obtained as a quadrature sum of all the above, where the last five sources do not affect the Zemach terms. The values we thus obtained for the individual and total relative uncertainties estimated for δ A pol , δ A Zem , and δ A TPE in µ 3 He + and µ 3 H are given in Table S1 below. We remind the reader that δ A TPE ≡ δ A Zem + δ A pol can be obtained directly from our results, by summing all terms in δ A pol except for the Zemach terms, due to their cancellation. Table S1. Estimated relative uncertainties, in percents, assigned to the calculated nuclear TPE corrections to the 2S-2P Lamb shift in µ 3 He + and µ 3 H. The presented values are rounded. The total uncertainties are obtained from a quadrature sum.