The deuteron-radius puzzle is alive: A new analysis of nuclear structure uncertainties

To shed light on the deuteron radius puzzle we analyze the theoretical uncertainties of the nuclear structure corrections to the Lamb shift in muonic deuterium. We find that the discrepancy between the calculated two-photon exchange correction and the corresponding experimentally inferred value by Pohl et al. [1] remain. The present result is consistent with our previous estimate, although the discrepancy is reduced from 2.6 $\sigma$ to 2 $\sigma$. The error analysis includes statistic as well as systematic uncertainties stemming from the use of nucleon-nucleon interactions derived from chiral effective field theory at various orders. We therefore conclude that nuclear theory uncertainty is more likely not the source of the discrepancy.


Introduction
The charge radius of the deuteron (d), the simplest nucleus consisting of one proton and one neutron, was recently determined to be r d = 2.12562(78) fm [1] using several Lamb shift (LS) transitions in muonic deuterium (µ − d). This result provides three times the precision compared with previous measurements. Furthermore, the µ − d value is 7.5 σ or 5.6 σ smaller than the world averaged CODATA-2010 [2] or CODATA-2014 [3] values, respectively, and 3.5 σ smaller than the result from ordinary deuterium spectroscopy [4]. One can also combine the measured radius squared difference r 2 d − r 2 p obtained from isotope shift experiments on ordinary hydrogen and deuterium [5] with the absolute determination of the proton radius from muonic hydrogen experiments [6,7] (dubbed as "µp+iso") to obtain r d = 2.12771 (22) fm, which is much closer to the µ − d result, but still differs from it by 2.6 σ (see Ref. [1] for details). Altogether, these significant discrepancies have been coined "the deuteron radius puzzle". Unlike with the proton-radius puzzle [6], r d from µ − d Lamb shift measurements is consistent with the electron-deuteron (e − d) scattering data due to the large uncertainty in the scattering experiments. Ongoing efforts to improve the precision in electron scattering will provide further information [8]. However, these discrepancies, compounded with the 7 σ (5.6 σ) discrepancy between the CODATA-2010 (CODATA-2014) and Email addresses: javierh@phas.ubc.ca (O. J. Hernandez), andreas.ekstrom@chalmers.se (A. Ekström), nnevodinur@triumf.ca (N. Nevo Dinur), jichen@mail.ccnu.edu.cn (C. Ji), s.bacca@uni-mainz.de (S. Bacca), nir@phys.huji.ac.il (N. Barnea) the muonic hydrogen proton radius [6,7], highlight the need to pinpoint the source of the differences. While the very recent 2S − 4P spectroscopy on ordinary hydrogen supports the small proton radius [9], the conundrum of the proton and deuteron radius puzzles is not yet fully solved and further experimental and theoretical investigations are clearly required. The deuteron charge radius r d is extracted from the LS measurement through which is valid in an α expansion up to 5 th order, where α is the fine structure constant. The term m r in Eq. (1) is the reduced mass of the µ−d system. The LS energy difference, ∆E LS , is directly measured through pulsed laser spectroscopy experiments described in detail in [1,6,7,10]. The quantum electrodynamic (QED) corrections δ QED are obtained from highly accurate theoretical calculations [11,12]. In the extraction of r d from LS measurements the main source of uncertainty is due to nuclear structure corrections coming from a two-photon exchange (TPE) diagram, δ TPE depicted in Fig. 1. Because the latter is obtained from theoretical computations it is of paramount importance that all theoretical uncertainty contributions are thoroughly investigated. Several groups have calculated δ TPE with different methods [13][14][15][16][17][18]. The two most recent and most precise computations, Ref. [17] and Ref. [18], are consistent within 0.6%. All theoretical calculations have been summarized by Krauth et al. [19] which resulted in a recommended value of δ TPE =-1.7096(200) meV. This value was also used by Pohl et al. [1] to extract r d using Eq. (1).
On the other hand, measuring ∆E LS and knowing δ QED , Eq. (1) enables the extraction of δ TPE from an experimentally determined radius. Using r d from "µp+iso" leads to an experimental value δ TPE =-1.7638(68) meV [1], which differs from the theoretical one by 2.6 σ. This disagreement motivates a reassessment of the theoretical calculation, and in particular of its assigned uncertainties.
Finally, from the radii of light nuclei, such as hydrogen and deuterium, it is possible to determine the Rydberg constant R ∞ and consequently the radius puzzle can be turned into a "Rydberg constant puzzle". In Ref. [4] two values of R ∞ were calculated using the muonic hydrogen and muonic deuterium charge radii separately and the results were found to disagree by 2.2 σ. This difference was attributed to the δ TPE contribution used to extract r d from the Lamb shift.
The purpose of this letter is to revisit our calculations of the nuclear structure corrections in µ − d and exploit chiral effective field theory and statistical regression analysis to systematically improve the theoretical uncertainty estimation in δ TPE and shed light on the deuteron radius puzzles.
State-of-the-art calculations of δ TPE in Refs. [16,18] as well as in this work, employ nucleon-nucleon (NN) potentials derived from a low-energy expansion of quantum chromodynamics called chiral effective field theory (chiral EFT). Within this approach, which also constitutes the modern paradigm of analyzing the nuclear interaction, the nuclear potential is built from a sum of pion-exchange contributions and nucleon contact terms, see, e.g., Refs. [20,21]. Power counting enables to determine the importance of individual terms in the low-energy expansion and thereby also facilitates a meaningful truncation of higher-order diagrams that build the potential. All potentials in this work employ Weinberg's dimensional power counting schemes [22,23], whereby the order ν ≥ 0 to which a diagram belongs is proportional to Q ν , where and p is a small external momentum, Λ b is the chiral symmetry breaking scale of about the rho meson mass, and m π is the pion mass. Given a power counting, contributions with a low power of ν are more important than terms at higher powers. Starting from the leading order (LO), i.e., ν = 0, higher orders will be denoted as next-to-leading order (NLO), i.e., ν = 2, next-tonext-to-leading order N 2 LO, i.e., ν = 3, and so on. It is worth noticing that, in chiral EFT, the contributions with ν = 1 vanish due to time-reversal and parity. At each order ν of the chiral EFT potential, there will be a finite set of parameters, known as low energy constants (LECs), that determine the strength of various pion-nucleon and multi-nucleon operators. The LECs are not provided by the theory itself but can be obtained from fitting to selected experimental data, such as NN and πN scattering cross sections, and other few-body ground state observables, such as radii and binding energies. Different fitting procedures exist, and we will explore a variety of them as a way to probe both statistical and systematic uncertainties.
To avoid infinities upon iteration in the Lippmann-Schwinger equation all chiral potentials are regulated by exponentially suppressing contributions with momenta p greater than a chosen cutoff value Λ, see e.g. Refs [20,21]. Non-perturbative ab initio calculations using momentum-space chiral EFT often employ NN interactions with Λ ≈ 400 − 600 MeV. Chiral EFT and effective field theories in general, unlike phenomenological models, furnish a systematic, i.e, order-byorder, description of low-energy processes at a chosen level of resolution. In this work, it provides us with an opportunity to estimate the uncertainty of δ TPE truncated up to different chiral orders. In our previous work [16,18], we probed the theoretical uncertainty stemming from the nuclear physics models by cutoff variation, i.e., varying Λ. Strictly speaking, this prescription to estimate the chiral EFT uncertainty also requires the excluded νth-order chiral contributions to be proportional to 1/Λ ν+1 when Λ is approaching the breakdown scale Λ b : a property that hinges on order-by-order renormalizability of the canonical chiral EFT formulation, which is not yet established. Also, cutoff-variation tend to either underestimate or overestimate the chiral EFT systematic uncertainty with respect to the variation range [24,25]. To this end, and to be as conservative as possible, we will augment the procedure of cutoff-variation by implementing the chiral EFT truncation-error to obtain solid systematic uncertainty estimates.
Any rigorous estimate of the theoretical uncertainty must also consider the effects of the statistical uncertainties of the LECs due to experimental uncertainties in the pool of fitted data. For example, in Ref. [26] it was found that a rigorous statistical analysis lead to a four-fold increase in the uncertainty estimates of the proton-proton fusion S -factor as compared to previous work which only probed the systematic uncertainty of the nuclear model by limited cutoff variations. Motivated by the possibility that the uncertainties were underestimated, we rigorously probe the statistical and systematic uncertainties in the nuclear structure corrections in the Lamb shift of µ − d, by propagating the uncertainties of the LECs appearing in the NN potentials up to N 2 LO [27,28].
Details on the observables associated with the LS in µ − d are explained in Section 2 and results of the statistical analysis will be shown in Section 3. In addition, we improve our estimates of the systematic uncertainty associated with the chiral EFT expansion by carrying out our calculations up to fifth-order in chiral EFT, namely N 4 LO. We then use the method detailed in Refs. [24,25,29] to estimate the systematic uncertainty associated with the chiral truncation at each order. Results will be shown in Section 4. Finally, we will examine and combine all the relevant sources of uncertainty in Section 5, before drawing conclusions in Section 6.

Two-photon exchange contributions
For the calculation of δ TPE we separate terms that depend on the few-nucleon dynamics, denoted with A, from terms that exclusively depend on properties of the single-nucleon, denoted with N, as The first contribution, which is the focus of this paper, can be written as Here δ (0) is the leading, δ (1) the subleading and δ (2) the subsubleading term in the context of an η-expansion. The small parameter η is √ m r /m d with m d being the deuteron mass. The term δ A Zem is the third Zemach moment. These first four terms are calculated in the point-nucleon limit, while nucleon size (NS ) corrections are added with the additional terms δ (1) NS and δ (2) NS . The formulas of these corrections are detailed in Refs. [16,18,30,31] and not repeated here.
The single-nucleon terms, which we take from the literature, can be decomposed into The nucleon Zemach moment and polarizability terms are estimated from data and taken to be δ N Zem = −0.030(2) meV 1 and δ N pol = −0.028(2) meV [15]. The subtraction term δ N sub instead is not well constrained by experimental data. Its contribution in µH was calculated by Birse and McGovern using chiral perturbation theory to be 0.0042(10) meV [32]. With the operator product expansion approach, Hill and Paz recalculated the TPE in µH [33], obtaining a central value that agrees with Ref. [32] but has a much larger uncertainty (see also Refs. [34,35]). For the nucleonic subtraction in µD, we adopt the strategy of Ref. [19] assuming that proton and neutron subtraction terms are of the same size and assigning a 100% uncertainty. This yields δ N sub = 0.0098(98) meV [19]. If we enlarge the uncertainty to 200% to be comparable with the finding in Ref. [33], we have δ N sub = 0.0098(196) meV. The calculations for the deuteron are based on a harmonic oscillator expansion of the wave function and a discretize approach to the sum rules to compute the terms in Eq. (4), see Refs. [16,18]. We have previously shown that this approach is reliable and agrees very well with other calculations, such as from Pachucki [13] and Arenhövel [36].

Statistical uncertainty estimates
To probe the uncertainties in δ TPE due to statistical uncertainties in the LECs -originating from the uncertainties in the pool of fitted data -we use the N k LO sim potentials from Ref. [28], with k from 0 to 2 2 . The presently available sim potentials employ seven different cutoff values Λ = 450, 475, . . . , 575, 600 MeV, and for each cutoff a potential was simultaneously fit to six increasingly larger energy-ranges of the SM99 world database of NN scattering cross sections, along with πN scattering cross sections, and ground state properties 3 of the 2,3 H and 3 He systems. The various subsets of the NN scattering data are delimited by the maximum kinetic energy in the laboratory frame of reference; T max Lab = 125, 158, 191, 224, 257, 290 MeV. The statistical covariance matrix of the LECs for each sim interaction was determined numerically to machine precision using forward-mode automatic differentiation. The covariance matrices make it possible to propagate the statistical uncertainties of the LECs to, e.g., δ TPE , and the various Λ and T max Lab cuts gauge the systematic uncertainties in the fitting procedure and the cutoff choice.
We compute the covariance matrix of the nuclear structure corrections relevant to the µ − d system using the linear approximation. For two quantities A and B, matrix elements of the covariance matrix are then obtained from Here, J A is the row vector of partial derivatives of the quantity A, with respect to the set of LECs α, and analogously for J B . The covariance matrix Cov(α) of LECs is provided from [28]. The vector components J A,i = ∂A ∂α i are obtained from a univariate spline fit to ten numerical function evaluations of A in a small neighborhood of the optimal value of the LEC α i . To benchmark our procedure for calculating the derivatives in this fashion we compared our results for the deuteron ground state properties, such as the ground state energy E 0 , the root meansquare point-nucleon radius r, and the quadrupole moment Q d , with existing computations based on automatic differentiation algorithms and obtained excellent agreement for all quantities, i.e. better than 0.005% relative uncertainty. The linear correlation coefficient between A and B is given by where σ A,stat ≡ √ Cov(A, A), and similarly for σ B,stat , are the statistical uncertainties of A and B, respectively. A value A regression analysis provides an exact quantification of the statistical correlations that are present in the chiral EFT description of the nuclear polarization effects. This allows us to determine constraints between different observables predicted within the N k LO sim models and serves as a valuable check of the uncertainty propagation formalism. In this work, we focus on δ A TPE and its components, such as the leading dipole term δ (0) D1 and the small magnetic term δ (0) M . Since they are related to the electric dipole and magnetic dipole response, we also study the electric dipole polarizability α E and the magnetic susceptibility β M . Detailed expressions can be found in Refs. [16,18]. We also compute other observables of interest, such as the ground-state energy E 0 and the mean square point-proton distribution radius r, related to r d by where r p/n are the proton and neutron radius, respectively, while m is the proton mass and r 2 2BC is a contribution due to twobody currents [37,38] and other relativistic corrections beyond the Darwin-Foldy term 3 4m 2 . Furthermore, we also study other ground-state observables, such as the electric quadrupole moment Q d and the magnetic dipole moment µ d , along with the D-wave probability P D . The linear correlation coefficients between these observables are presented in Fig. 2.  We observe that the quadrupole moment Q d of the deuteron is strongly correlated with the P D [39]. The magnetic moment of the deuteron ground state µ d can be calculated analytically and depends only on the P D probability [39]. As expected, from our numerical analysis, the correlation ρ(µ d , P D ) is strongly negative. Furthermore, from this correlation we also expect to see a correlation between P D and δ (0) M which we indeed observe. Based on this analysis, we see that the magnetic properties of the deuteron are strongly related to P D , indicating that they are largely determined by the D-wave component of the deuteron. By contrast, the electric properties are found to be related to r. For example, a correlation is observed between r and α E . This correlation is predicted on the basis of the zero-range model of the deuteron [40] and has been observed to hold in heavier systems [41,42]. In fact, we find that δ A TPE , which is strongly dominated by the dipole term δ (0) D1 , is correlated to r and α E . At this point the production of expected correlations within the formalism of the statistical uncertainty propagation served as a way to inspect the validity of our statistical analysis, but it could also be of guidance if one decided to use alternative fitting procedures in the future.
The statistical uncertainties of r, the electric dipole polarizability α E , and δ A TPE are 0.02%, 0.05%, and 0.05% respectively, i.e., negligible compared to the size of the systematic uncertainty, as we will show in the next Section. Importantly, if one uses the separately (or sequentially) optimized N 2 LO sep potentials of Ref. [28], which fit the LECs in the πN, NN and 3N sectors separately, the deduced statistical errors would be larger. This originates from neglecting the statistical covariances between the πN and NN (and 3N) sectors of the chiral EFT potentials while also employing LECs with rather large statistical uncertainties. Most microscopic interactions are constrained to data in such a sequential manner, in particular the chiral EFT potentials up to 5 th order that we employ in the next Section. However, the sub-leading πN LECs employed in those potentials were separately optimized using a novel Roy-Steiner extrapolation of the πN scattering data [43]. The resulting πN LECs exhibit very small statistical uncertainties [44], and a forward error propagation to δ A TPE would most likely lead to similarly reduced statistical uncertainties.

Systematic uncertainty estimates
Here we present the systematic uncertainties in our calculations due to various truncations introduced in chiral EFT. First, we address the truncation via T max Lab in the energy range of the NN scattering data used in the fit of the LECs. In Fig. 3, the calculated values of δ A TPE are plotted as a function of the maximum lab energy, T max Lab , in the fit of the N 2 LO sim potentials for various choices of the cutoff Λ. The error bars indicate the statistical uncertainties, computed as detailed in the previous Section, which were on average found to be 0.001 meV, or 0.06%. In Fig. 3 it is clear that the statistical uncertainties are small in comparison to the systematic uncertainties due the variation of the cutoff Λ and T max Lab . Furthermore, the range of the calculated values of δ A TPE for different Λ decreases at the largest T max Lab energies, indicating that the nuclear dynamics as described by the LECs become better constrained with more data.
Next, we address the uncertainties coming from truncating chiral EFT at the order ν. The common approach to gauge this uncertainty is by varying the cutoff Λ over a range of values. However, this approach to uncertainty estimation suffers  from several deficiencies, such as the arbitrariness in the chosen Λ−range. Furthermore, often the residual Λ dependence underestimates uncertainties as discussed in Refs. [24,25]. To address these deficiencies and give a conservative estimate of the systematic uncertainties, in addition to cutoff variation, we follow [25,29] and include an uncertainty estimate based on the expected size of the next higher-order contribution in the chiral EFT expansion. This approach is in semi-quantitative agreement with a Bayesian uncertainty analysis. Assuming that an observable A(p), associated with an external momentum scale p, and computed non-perturbatively from chiral EFT, follows the same order-by-order pattern as chiral EFT itself, then it can be expressed as where A 0 is the leading order value, Q is the small expansion parameter, given in Eq. (2), and c ν (p) is an observable-and interaction-specific expansion coefficient determined a posteriori. The uncertainty in A due to truncation at some finite order ν, i.e., LO, NLO, N 2 LO,. . ., can be estimated by This expression rests on a prior assumption of independent expansion coefficients c ν with a boundless and uniform distribution.
To estimate the typical momentum scale p of the nuclear structure corrections, we compute the average energy value of the largest term in δ A TPE , namely the leading order dipole correction [16] where S D1 (ω) is the dipole response function. The average value ω D1 is calculated as Given that we obtain ω D1 ≈ 7 MeV, which corresponds to a momentum scale p smaller than m π , the chiral convergence parameter Q for our uncertainty estimates is always taken to be m π /Λ b . For a solid estimation of truncation errors, it is also crucial to adopt a suitable Λ b . Here we follow [25,29] and set Λ b to 600 MeV, as a choice shown by Bayesian analyses [29,45]   In Fig. 4, the uncertainty estimates for δ A TPE are displayed for the N 2 LO sim potentials as a function of the cutoff Λ. The blue band shows a conservative spread of δ A TPE due to variations of T max Lab from 125 to 290 MeV, which, with respect to the central value, amounts to about 0.004 meV or 0.2%. The green band also includes the uncertainty stemming from the chiral truncation. Due to the fact that the two systematic uncertainties are not independent from each other, the chiral truncation error is calculated from Eq. (10) and added to each point in T max Lab and Λ. The green band encompasses the maximum and minimum values of δ A TPE so obtained. The largest contribution of the truncation uncertainty alone is found to be 0.007 meV for the N 2 LO sim potentials. The overall systematic uncertainty including cutoff variation, T max Lab variations and chiral truncation error amounts to 0.011 meV or 0.65% for the N 2 LO sim potentials and thus dominates with respect to the 0.06% statistical uncertainty.
To study the convergence of δ A TPE with respect to the chiral orders greater than N 2 LO, in addition to the sim potentials, we also carry out calculations using the chiral potentials available up to N 4 LO. Two groups of chiral interactions have been constructed using different fitting procedures and slightly dif-ferent operatorial form in the potentials, but identical powercounting. We will use all orders available from Ref. [24] 4 and denote them as N k LO EK M , and those from Ref. [46], which will be denoted as N k LO EMN . For the N k LO EK M family of potentials we will explore the cutoffs (R 0 , Λ) = (0.8, 600), (1.0, 600) and (1.2, 400) [fm, MeV], where R 0 is a coordinate-space regulator, and for the N k LO EMN family of potentials we will use Λ = 450, 500 and 550 MeV. The use of a higher order in chiral EFT will allow us not only to update our results with respect to our previous work [16,18], but also to get a more reliable estimate the chiral convergence uncertainty using Eq. (10). Our goal is to provide an updated value of δ TPE with its overall uncertainty. This will be discussed in the next Section.

Total uncertainty estimates
First, systematic and statistical uncertainties of the various nuclear interactions are combined into σ Nucl , which is detailed in Table 1. For the N k LO EK M and N k LO EMN potentials at all orders we include systematic uncertainties from chiral convergence and cutoff variation. Systematic errors stemming from T max Lab variations cannot presently be estimated at N 3 LO and at N 4 LO, thus we include the corresponding uncertainty evaluated from N 2 LO sim . For the lowest orders in the EK M and EMN potentials, systematic errors from T max Lab variations are taken from the corresponding order of the sim interactions. The sim potentials contain all of the above and statistical uncertainties, estimated consistently at each order. Statistical uncertainties are found to be negligible in the sim potentials and, while at present a consistent evaluation is not possible, they are also expected to be small in the N k LO EK M and N k LO EMN families. Thus, we take the N 2 LO sim statistical values also for the The present calculation is performed to sub-subleading order in the η-expansion, thus uncertainties deriving from higher order η contributions need to be estimated in the total error budget. These higher order η contributions are estimated to provide a 0.3% effect based on a different approach to the computation of δ A TPE , which allows to include higher order electromagnetic multipoles [47]. So far we have concentrated on δ A TPE , which is the only term in δ TPE with explicit dependence on the nuclear dynamics. For a complete discussion on δ TPE we should consider the additional nucleonic terms in Eq. (5), namely δ N pol , δ N Zem , δ N sub and their respective uncertainties, using the values quoted in Section 2. Finally, our δ TPE formulas are valid in an α expansion up to 5 th order. Higher order terms in the α expansion were estimated first by Pachucki [13] to be of 1%. Here, we will keep this value and refer to it as the atomic physics uncertainty, as in Refs. [16,18], and add it to the other uncertainties in quadrature. Atomic, single-nucleon, and η-expansion uncertainties of δ TPE are included in σ Tot , on top of the nuclear physics uncertainties, given in Table 1. In Fig. 5 and in Table 1, we show the convergence of δ TPE and its overall uncertainties with respect to all chiral orders from LO up to N 4 LO using δ N sub from Ref. [19] (we will include the larger uncertainties of Ref. [33] later in our analysis). In particular, in Table 1, one can appreciate the difference between σ Nucl and σ Tot . One can readily see that uncertainty bands decrease as the order of the chiral expansion increases, as expected. Beside LO calculations, where the three potential families somehow differ, all results at higher orders are quite stable around the same value, independently of the potential used. Interestingly, N 4 LO results are almost identical to N 3 LO results, indicating convergence of the chiral expansion for this observable. The uncertainty estimates at N 3 LO and N 4 LO are compatible with our previous estimates in Refs. [16,18], even though slightly larger, mostly due to the inclusion of the systematic error using Eq. (10). Furthermore, Table 1 shows that, although the nuclear physics errors are dominant at lower order in the chiral expansion, at N 4 LO the leading uncertainty is not stemming from nuclear physics, but rather from the other sources.
Results are also compared to the experimentally inferred δ TPE correction [1] and theoretical compilation [19]. We find that the N 4 LO band is consistent with the theoretical compilation and encompasses also our result δ TPE = −1.709 meV from Ref. [16] based on the AV18 potential [48], which is also included in the theory summary by Krauth [19]. We also observe, though, that our N 4 LO band is not compatible with the experimental determination of δ TPE .
Finally, based on our analysis, we provide an updated value of δ TPE = −1.715 meV with its itemized uncertainty budget in Table 2   uncertainties from cutoff variation and chiral truncation are obtained from our N 4 LO studies by taking the combined range of the EMN and EK M bands. While here we studied the chiral convergence of the potential, the same should in principle be done regarding the chiral expansion of electromagnetic currents [37,38] leading to further systematic corrections. In our formalism, δ A TPE is related to electromagnetic multipoles, where the electric dipole dominates. The latter is protected by the Siegert theorem [36], so that two-body currents are implicitly included via the use of the continuity equation [36]. There exists corrections to the Siegert term. We have estimated the magnitude of those by integrating an E1-response function provided by Arenhövel [49], which included two-body currents and relativistic corrections as in Ref. [36] for the AV18 potential. Their effect on the leading dipole correction δ (0) D1 was both found to be of the order of 0.05%, thus negligible.
Despite the disputed single-nucleon TPE uncertainty [33][34][35], it is evident from Table 2 that the atomic physics error remains a major source of uncertainty. It is approximately 1% from a reasonable, but rough, estimate by Pachucki et al. [13]. The estimate is based on taking 50% of relativistic and higher order corrections, which are the smallest contributions to δ A TPE . A more thorough estimate of α 6 effects requires going to third order in perturbation theory and study three-photon exchange effects, which is beyond the scope of this work. Here, we have shown that uncertainties stemming from the chiral EFT description of the nuclear interaction alone are not capable of explaining the discrepancy between the calculated δ TPE and the corresponding experimental extraction by Pohl et al. [1].
Since the deuteron point-nucleon radius r based on CODATA was used in the fitting procedure, e.g., of N 2 LO sim potentials, one may suspect this yields a biased δ TPE in muonic atoms. However, in order to remedy the discrepancy between the µ − d and "µp+iso" values of r d , one may just vary r by ∼ 0.1%, see Eq. (8). Due to the linear correlation between r and δ TPE , this would only lead to a maximum variation δ TPE of the order of ∼ 0.1%, which is negligible with respect to the required ∼ 3% change needed to explain the discrepancy between its theory and experimental determination.

Conclusions
In this work, we have explored the uncertainties in δ TPE corrections using state-of-the-art chiral potentials from LO up to N 4 LO. We have calculated the statistical uncertainties up to N 2 LO using a set of simultaneously optimized chiral potentials. From this we conclude that the uncertainty due to variances of the LECs are negligible compared to the systematic uncertainties due to cutoff variation and chiral truncation. We have also found that going beyond N 3 LO in chiral EFT does not change the overall results of δ TPE , which also indicates a high theoretical accuracy of our final result. In conclusion, the rigorous uncertainty quantification presented here weakens the disagreement between the calculated two-photon exchange correction and the corresponding experimentally inferred value by Pohl et al. [1] from 2.6 σ to within 2 σ (or 1.7 σ if using the larger single nucleon uncertainties of Ref. [33]). Breaking down the total uncertainty budget in the calculation of δ TPE shows that atomic physics and single-nucleon physics need to be addressed to further reduce the theoretical uncertainty. It is important to remark that the deuteron-radius puzzle is still alive, in that the large discrepancy between the spectroscopic measurements on muonic deuterium and on ordinary deuterium still exists and it does not seem to be simply explained from nuclear physics uncertainties in the few-body dynamics.

Acknowledgments
We would like to thank Angelo Calci for providing us with the chiral potentials at N 4 LO. We are grateful to Randolf Pohl and Hartmuth Arenhövel for useful discussions.