Two-photon exchange in (muonic) deuterium at N3LO in pionless effective field theory

We present a study of the two-photon-exchange (2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upgamma $$\end{document}γ-exchange) corrections to the S-levels in muonic (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}μD) and ordinary (D) deuterium within the pionless effective field theory (EFT). Our calculation proceeds up to next-to-next-to-next-to-leading order (N3LO) in the EFT expansion. The only unknown low-energy constant entering the calculation at this order corresponds to the coupling of a longitudinal photon to the nucleon–nucleon system. To minimise its correlation with the deuteron charge radius, it is extracted using the information about the hydrogen–deuterium isotope shift. We find the elastic 2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upgamma $$\end{document}γ-exchange contribution in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}μD larger by several standard deviations than obtained in other recent calculations. This discrepancy ameliorates the mismatch between theory and experiment on the size of 2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upgamma $$\end{document}γ-exchange effects, and is attributed to the properties of the deuteron elastic charge form factor parametrisation used to evaluate the elastic contribution. We identify a correlation between the deuteron charge and Friar radii, which can help one to judge how well a form factor parametrisation describes the low-virtuality properties of the deuteron. We also evaluate the higher-order 2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upgamma $$\end{document}γ-exchange contributions in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}μD, generated by the single-nucleon structure and expected to be the most important terms beyond N3LO. The uncertainty of the theoretical result is dominated by the truncation of the EFT series and is quantified using a Bayesian approach. The resulting extractions of the deuteron charge radius from the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}μD Lamb shift, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2S{-}1S$$\end{document}2S-1S transition in D, and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2S{-}1S$$\end{document}2S-1S hydrogen–deuterium isotope shift, with the respective 2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\upgamma $$\end{document}γ-exchange effects evaluated in a unified EFT approach, are in perfect agreement.

while the most accurate extraction of the deuteron charge radius, is an indirect achievement combining measurements from the spectroscopy of ordinary and muonic atoms [2]: the 2S-1S hydrogen-deuterium (H-D) isotope shift and the Lamb shift in μH. This result is driving the presently recommended value of the deuteron charge radius from the CODATA 2018 report [4]: As one can see from Eq. (1), the charge radius extractions are limited by the theory uncertainty, which for muonic atoms is almost solely due to subleading nuclear-structure effects, and in particular, the O(α 5 ) two-photon exchange (2γ exchange) discussed in this work. The initial tension between the r d (μD) and r d (μH & iso) extractions, shown above, was resolved in 2019 by amending the μD theory [5] to include the subleading O(α 6 ) electronic vacuum polarization (VP) effects [6] and the inelastic threephoton exchange (3γ exchange) [7]: r d (μD) = 2.12710(13) exp (81) theory fm = 2.12710(82) fm [6].
This distinct discrepancy for the deuteron radius -the "deuteron radius puzzle" -is strongly affected by the 2γ exchange. It is thus timely to re-evaluate the 2γ-exchange effects in a model-independent manner and try to improve their precision. While the latest developments [6,7] are certainly important, they do not provide a path to a more systematic improvement of the theory error on the side of nuclear structure.
In this work, we consider the forward 2γ-exchange contributions to D and μD, including the accompanying electronic VP contributions, within the pionless effective field theory ( / πEFT) of nuclear forces [11][12][13][14][15][16][17]. This framework allows one to represent the nuclear observables in a well-defined perturbation theory, expanding in powers of the small parameter P/m π , where P is the typical momentum scale (e.g., the size of the relative momentum between two nucleons, or that of the momentum of an external probe) and m π 139 MeV is the pion mass. The typical momentum scale in the deuteron is characterized by the binding momentum γ = √ M N B 45 MeV, where M N is the nucleon mass and B is the deuteron binding energy. The momentum scale probed by the electromagnetic interaction in μD is ∼ αm μ , which is less than an MeV. This is also well below the limiting scale of the theory set by m π . The atomic systems should thus be well-suited for the application of / π EFT. In addition, it has been shown that / π EFT provides a good description of low-energy experimental data on real deuteron Compton scattering [18,19], and can be used to investigate the deuteron electric polarizability and electromagnetic form factors (FFs) [15]. Finally, the effective-field-theory (EFT) expansion allows one to quantify the theoretical uncertainty using methods such as Bayesian inference [20]. The basis for our 2γ-exchange calculation is provided in Ref. [21], where closed analytic expressions for the unpolarized amplitudes of forward doubly-virtual Compton scattering (VVCS) off the deuteron are derived.
The paper is organized as follows. In Sect. 2, we briefly introduce the 2γ-exchange and / π EFT frameworks. In Sect. 3, we calculate the elastic finite-size, inelastic deuteronpolarizability and single-nucleon contributions to the μD Lamb shift, and compare our results to other recent predictions. In Sect. 4, we repeat the same calculation for D and use the H-D isotope shift to fix the unknown low-energy constant (LEC) l C0 S 1 that enters the VVCS amplitude. In Sect. 5, we utilize the unique possibility to cross-check the theoretical predictions for the 2γ exchange in μD with an empirical determination. The latter is extracted from the measured μD Lamb shift by fixing the deuteron charge radius to the independent value r d (μH & iso). We also compile an update for the theory prediction of the μD Lamb shift that will be used to extract r d (μD) from the measurement of the CREMA Collaboration. In Sect. 6, we discuss deuteron and proton charge radii extractions from μD, D and the H-D isotope shift. A discussion of the neutron charge radius is postponed to Appendix G. In Sect. 7, we finish with conclusion and outlook. The appendices cover: (A) the Bayesian error analysis, (B) the inclusion of nucleon FFs beyond the / π EFT framework, (C) electronic VP corrections, and updated theory com-  The leading order (LO) in α 2γ-exchange correction corresponds to the forward kinematics, shown in Fig. 1. It gives a δ(r)-function correction to the Coulomb potential, thus, only shifts the S-levels, which have a non-vanishing atomic wave function at the origin. The spin-independent forward 2γ exchange is related to the VVCS amplitude off an unpolarized deuteron: where f L (ν, Q 2 ) and f T (ν, Q 2 ) are the longitudinal and transverse scalar amplitudes with Q 2 = −q 2 and ν = p · q/M d the photon virtuality and lab frame energy, and M d the deuteron mass. The modified photon polarization vector components are defined as with q andq = q/|q| being the photon three-momentum in the lab frame and its unit vector, and ( 0 , ) the time and space components of the photon polarization vector. This description in terms of f L (ν, Q 2 ) and f T (ν, Q 2 ) is natural for the / π EFT framework, but not unique. The explicitly covariant tensor decomposition with two other scalar amplitudes T i (ν, Q 2 ) related via is widely used in, e.g., the dispersive 2γ-exchange evaluations [23,24]. We start from the covariant expression for the forward O(α 5 ) 2γ-exchange correction to the energy of a nS state in (muonic) deuterium, given in these references, and rewrite them in terms of the longitudinal and transverse amplitudes: where m is the electron or muon mass, [φ n (0)] 2 = 1/(π n 3 a 3 ) is the (Coulomb) wave function of the nS atomic state at the origin, a = 1/(Zαm r ) is the Bohr radius, Z is the nuclear charge (Z = 1 for the deuteron), and m r = m M d /(m + M d ) is the atomic reduced mass. Separating the scalar amplitudes into the deuteron-pole and non-pole parts, one splits the 2γ-exchange correction into the elastic and inelastic part [24]. The inelastic part, after doing the Wick rotation ν = iq 0 and introducing the hyperspherical coordinates, takes the form: with τ l = Q 2 /(4m 2 ). Here we assume that the pole-part is subtracted from the scalar VVCS amplitudes. The elastic part of the 2γ exchange is readily obtained via the deuteron electromagnetic FFs -charge, magnetic, and quadrupole - , and G Q (Q 2 ), resulting in [24]: where τ d = Q 2 /(4M 2 d ), and the weighting functions are defined by: Note that the contributions of point-like charge and charge radius of the deuteron are removed from the elastic part to avoid double counting [24]. This is done by subtracting the unity and the term proportional to G C (0) in Eq. (11).

Unpolarized deuteron VVCS in pionless EFT
In our analysis, we use results from the / π EFT calculation of the unpolarized deuteron VVCS amplitudes f L (ν, Q 2 ) and f T (ν, Q 2 ) presented in Ref. [21]. This section gives a brief recap of the / π EFT framework applied to the deuteron VVCS, as well as a description of the technicalities relevant to the 2γ-exchange calculation.
/ π EFT is an EFT for nucleon interactions at low energies, where the high-energy scale is set by the pion mass m π . If the momentum transfer between two nucleons is P m π , one can treat a pion-exchange interaction as a contact one. In / π EFT nucleons are thus interacting through contact interactions [11][12][13]15]. The Lagrangian is constructed performing a non-relativistic expansion in the one-nucleon sector and writing out the relevant two-nucleon interactions [11][12][13][14][15][16][17]. To assign a particular order to a Feynman graph, one counts powers of momenta [Q = O(P)] and energies [ν = O(P 2 )] coming from the interaction vertices, nucleon propagators [O(P −2 )], and loops [O(P 5 )]. The small expansion parameter is the ratio P/m π . For the deuteron, where the typical momentum scale is the binding momentum γ , this corresponds to P/m π 1/3. Note that different momentum scales can count as different powers of the typical momentum P, depending on the problem setting. For instance, the counting we use has the photon three-momentum |q| = O(P), whereas its energy is ν = O(P 2 ), and hence also its virtuality Q = O(P). This reflects our expectation that the virtual photons in the 2γ -exchange, as viewed in the lab frame, mostly transfer three-momentum, and very little energy, to the intermediate deuteron state, and is in contrast to, e.g., a typical real Compton scattering setting where ν = |q|, implying they have to be of the same size in the counting.
Regarding the description of the deuteron state, one can use different prescriptions to perform the expansion around the deuteron pole of the nucleon-nucleon (N N) scattering amplitude. The z-parametrisation [17], chosen in Ref. [21, Sec. II B], is particularly well-suited for quantities such as the deuteron electric dipole polarizability α E1 that receive mostly long-range contributions and are thus sensitive to the correct description of the long-range tail of the deuteron wave function. This parametrisation reproduces the residue Z of the N N scattering amplitude at the deuteron pole at next-toleading order (NLO). The residue is related to the effective range ρ d in the N N triplet channel via Z = (1 − γρ d ) −1 , and is also connected to the asymptotic normalisation of the deuteron S-wave via It is therefore straightforward to see that this procedure also reproduces the correct large-distance asymptotics at NLO. Note also that it introduces a new formal expansion parameter Analysing the counting for the VVCS shows [21, Sec. II B] that the longitudinal amplitude, driven by the deuteron electric polarizability α E1 , is dominant, starting at O(P −2 ), whereas the transverse amplitude starts two orders higher at O(P). In the context of the 2γ-exchange correction, Eq. (9b) shows that the f T contribution is additionally suppressed, compared to the contribution of f L , by the factor ν 2 /Q 2 = O(P 2 ). The transverse contribution to the 2γ-exchange correction therefore starts only at O(P 2 ), or N4LO compared to the leading longitudinal contribution. It is also at N4LO that, as explained in Ref. [21], higher powers of momenta entering the / π EFT expansion render the 2γ-exchange correction naively divergent. This divergence ought to be absorbed by a four-nucleon two-lepton contact term entering at this order, and, since there is no data that would allow one to pinpoint the corresponding coupling other than the 2γ-exchange correction itself, this is where the predictive power of / π EFT is exhausted. This motivated us to calculate the longitudinal amplitude up to N3LO in Ref. [21]. We also calculated the transverse amplitude up to O(P), or its respective NLO; this allows us to quantify here the corresponding 2γ-exchange contribution.
Further details of the / πEFT framework used to calculate f L and f T at their respective N3LO and NLO can be found in Ref. [21,Sec. II]. The results for the VVCS amplitudes are given in a closed analytic form in Ref. [21,Sec. III], in terms of the longitudinal and transverse four-point functions M L ,T (ν, Q 2 ) and the inverse of the derivative of the deuteron self-energy (SE) at the deuteron pole We use them here to calculate the 2γ-exchange correction. At N3LO in the / π EFT expansion of f L , one encounters a previously undetermined coupling l C0 S 1 of a longitudinal photon to the two-nucleon system, which contributes, in particular, to G C (Q 2 ) and r d . The latter quantity was used in Ref. [21] to extract l C0 S 1 from a fit to r d (μD) in Eq. (1b). This procedure is potentially problematic due to the fact that l C0 S 1 also enters the 2γ-exchange correction, both in μD and the isotope shift.
We investigate the resulting correlations below and demonstrate explicitly that they are negligible at the current level of theoretical and experimental precision.

2γ exchange in muonic deuterium
In the following section, we will present in detail our calculation of the elastic and inelastic 2γ-exchange contributions to the Lamb shift in μD. A summary of our results can be found in Sect. 3.4.

Elastic contribution
We start by considering the elastic contribution to the 2γexchange correction based on the / π EFT deuteron FFs. Taking the N3LO result for G C (Q 2 ) in Ref. [21,Eq. (75)] and expanding in Eq. (11) also to N3LO results in This neglects the magnetic and quadrupole FFs, whose contributions are subleading in the / π EFT counting and are indeed numerically very small, see below in Table 2. The electric contact term coupling l C0 S 1 can be fixed through the deuteron charge radius: with being the isoscalar nucleon charge radius, with the proton Darwin-Foldy term 3 /8 M −2 p added to it. Previously, l C0 S 1 was chosen to reproduce the deuteron charge radius from μD spectroscopy, resulting in [21] where the uncertainty stems from the error of the deuteron radius, Eq. (1b), 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). This correlation can be practically eliminated if the deuteron radius extracted from the combination of the proton radius and the H-D 2S−1S isotope shift, as given in Eq. (2), is used as the reference data point. One has to note that the isotope shift also has a 2γ-exchange contribution, but its relative importance as well as its correlation with r 2 d is much smaller. To investigate this issue quantitatively, we perform a re-analysis of the isotope shift using the / π EFT formalism to predict the 2γ-exchange correction in ordinary D, see Sect. 4 and Appendix D. Our calculation confirms that the contribution of l C0 S 1 to the isotope shift can indeed be safely neglected. The corresponding result for the electric contact term coupling, which will be used throughout this work, is This agrees with the result that we deduced from Eq. (2) [25], but differs from Eq. (16) by about 1 σ , since the extraction via the isotope shift gives a value of r d (μH & iso) slightly different from r d (μD) in Eq. (1b). The related effect on E elastic 2S is small. The final numerical result for the elastic contribution is: where the numbers here stand for the order-by-order contributions. The uncertainty of E elastic 2S is due to higher orders in the / π EFT expansion; we quantify it as explained in Appendix A.
To study the elastic (and inelastic) contribution in detail, we split them as shown in Table 1, keeping track of different terms appearing both due to the (Z − 1) factors coming from the NLO piece of (E d ) −1 and due to new sources at each order in the longitudinal four-point function M L . This representation will also be useful below in the investigation of the theoretical uncertainty. In order to split the elastic term this way, it is convenient to rewrite the last term in Eq. (11) replacing , where the normalization of G C (0) = 1 is used. This reflects the fact that the elastic part of the VVCS amplitude is proportional to the deuteron FFs squared, and allows one to separate the contributions in the integrand without generating spurious singularities at Q = 0. One can see from the table that the most important contributions to E elastic 2S come, as expected, from the LO part of M L , with the nucleon charge radius contributions providing the most important correction at N2LO. One can also see that the only new contributions beyond NLO come either from the nucleon structure or from the N3LO contact term proportional to l C0 S 1 . The nucleon charge radius contributions may seem somewhat larger than expected at N2LO and N3LO; to judge whether this is an indication of potentially sizeable corrections to E elastic 2S beyond N3LO, it is instructive to look at the details of the deuteron charge FF at small photon virtualities. Indeed, it is evident from Eq. (11) that the 2γ-exchange integrand is strongly weighted towards low Q 2 . Therefore, it is the slopes and the curvatures of the deuteron FFs at Q 2 = 0 that will have significant influence on the elastic contribution. The slope of the charge FF, proportional to r 2 d , is reproduced at N3LO; based on that alone, a sizeable modification of the shape of G C (Q 2 ) at small Q 2 could come from higher-order coefficients in its expansion  [21]. Quantities without labels are the contributions at the respective order excluding the labelled terms listed separately. Labels indicate specific terms within M L : r 2 N , w 2 , P, and l C0 S 1 stand, in order, for the nucleon charge radii correction, the contribution proportional to the N N triplet S-wave shape parameter w 2 , the N N P-wave contribution, and the contribution proportional to Starting with the recent higher-order, N4LO in the respective counting, chiral effective theory (χ ET) calculation of Refs. [26,27], a good agreement between the N3LO / π EFT and the N4LO χ ET results for G C (Q 2 ) at low Q was pointed out in Ref. [21]. As expected, our result for E elastic

2S
perfectly agrees with what one obtains using the χ ET charge FF from Ref. [27]: where we neglected the magnetic and quadrupole contributions, and evaluated the uncertainty using the uncertainty of the χ ET result for G C (Q 2 ). On the other hand, using the recent empirical deuteron FFs from Ref. [28], Carlson et al. obtained a considerably smaller value [24]: with the uncertainty estimated using the different FF parametrisations derived in [28]. The same result has been adopted in Ref. [29]. We repeat this calculation, separating the contributions from the charge, magnetic and quadrupole FFs. The results are presented in Table 2 (using parametrisation II of   Abbott et al.), along with values obtained by us based on the Sick and Trautmann parametrisation [9], as well as the χ ET and / π EFT FFs. One can see that the contributions of both the magnetic and quadrupole FFs to the elastic part of the 2γexchange correction can be safely neglected at the current level of precision. While the values of E elastic 2S obtained in / π EFT, χ ET, and with the Sick and Trautmann parametrisation of the deuteron FF agree, the parametrisation of Abbott et al. gives a significantly smaller value for the elastic contribution. The left panel of Fig. 2 shows that this discrepancy is due to the behaviour of the parametrisation of Ref. [28] being very different from the other three calculations (which would all overlap) at low Q.
Expanding the integrand in Eq.
Therefore, the bulk of the difference can be further traced down to the deuteron charge radius and the 4th moment of the deuteron charge density: G C (0) = −r 2 d /6 and G C (0) = r 4 d /60. Also interesting are two further quantities related to  (77) the elastic 2γ-exchange contribution, namely, the cubic and the Friar radii, defined respectively as [30]: In / πEFT at N3LO, the considered moments have the following analytic expressions, obtained using G C (Q 2 ) from Ref. [21,Eq. (75)] (note again that the integrand in r 3 Fd has to be expanded up to N3LO): (23c) Table 3 shows the values of these quantities for the considered FFs. It is evident that parametrisation II of Abbott et al. [28] gives smaller values for all radii. Smaller r d and r 4 d lead to a significantly smaller value of the integrand at low Q, as seen in the left panel of Fig. 2, and consequently a smaller E elastic 2S as well as smaller Friar and cubic radii. Note that, neglecting recoil corrections, the elastic contribution can be approximated through the Friar radius as [30] This approximation, however, results in a noticeable underestimation of E elastic

2S
. The / π EFT value, for instance, turns out to be E elastic, F 2S = −0.4323 meV, which has to be compared to Eq. (18). We therefore conclude that at the present level of theoretical precision it is important to retain the full weighting functionγ 2 (τ d , τ l ) in Eq. (11) instead of only taking the leading Friar radius term.
The dependence of both r 2 d and r 3 Fd on l C0 S 1 can be represented as a linear correlation between these quantities. We show the correlation line in the right panel of Fig. 2, where we also plot a ±1% ∼ (γ /m π ) 4 band as a simple estimate of terms beyond N3LO in the / πEFT expansion. One can see that the N4LO χ ET result lies almost on the correlation line, very close to the / π EFT results fixed by the H-D 2S−1S isotope shift, see Sect. 4 and Appendix D. The parametrisation of Ref. [9] lies some distance from the line, albeit reasonably close to it, whereas that of Ref. [28] is much further away. It would be interesting to see if this correlation line can be reproduced in a χ ET calculation.
The above considerations indicate that the FF parametrisation of Ref. [28], used in Refs. [24,29], might not adequately describe the behaviour of the deuteron charge FF at low virtualities. The agreement between the N3LO / πEFT and N4LO χ ET calculations, see Ref. [21,Sec. IV] for a detailed comparison of the FFs, is not entirely surprising as both these EFTs are expected to well reproduce low-momenta/longrange properties of the deuteron, and both calculations are of sufficiently high orders in the respective expansions. This vindicates our choice of the / π EFT as the calculational framework. One can also conclude that the correlation shown in Fig. 2 can serve as a diagnostic criterion for a realistic parametrisation of the deuteron charge FF. Furthermore, one can note that the / π EFT expression for the deuteron charge FF at N3LO, as given in Ref. [21,Sec. IV], can serve as an analytic one-parameter fit to the electron-deuteron scattering data in the low-Q 2 range that is to be covered in the planned DRad experiment [31].

Inelastic contribution
The calculation of the longitudinal contribution to inelastic part of the 2γ-exchange correction with the known f L (ν, Q 2 ) is straightforward. The only technical complication is that the longitudinal term in the integral of Eq. (10) goes as f L (0, Q 2 )/Q 3 ∝ 1/Q for Q → 0 when one sets x = 0. This singularity, however, is spurious, and can be avoided by the red disc, purple cross, green diamond, and blue square show the values obtained, respectively, from / πEFT at N3LO, the N4LO χET form factor [27], the parametrisation of Ref. [9], and the parametrisation of Ref. [28] subtracting from f L (ν, Q 2 ) its static part: The integration over x in the integral of f L (0, Q 2 ) can be done analytically, resulting in a f L (0, Q 2 )/Q 2 ∝ Q 0 behaviour for Q → 0. At the same time, the difference in the square brackets is O(x 2 ) at small x and therefore cancels the singularity in the weighting function. The longitudinal contribution then results in: One can see that the coefficients in front of l C0 S 1 here and in Eq. (14) partially cancel each other. The resulting contribution of the N3LO contact term to the 2γ exchange in μD is rather small. The numerical order-by-order result for E inel,L 2S , using l C0 S 1 as obtained from the H-D isotope shift, is: The uncertainty here is due to higher-order terms in the / π EFT expansion, calculated as explained in Appendix A.
The individual terms of the inelastic contribution are shown in Table 1, in an analogy to what is shown for the elastic part. While the bulk of E inel,L 2S is given by the LO part of M L , the most important correction comes from the nucleon charge radii, with the second-biggest correction driven by the NLO term of M L . The remaining mechanisms all give much smaller contributions.
The above results, Eqs. (26), (27), and Table 1, are obtained with the substitution |q| → Q in the expressions for M L ; using |q| = Q 2 + ν 2 brings the total value to −1.507 meV. This gives an estimate of the relativistic corrections at N4LO. The smallness of the effect corroborates the choice of counting scheme in our calculation, namely, that the energy transfer is suppressed and ν/Q = O(P). This statement can be made more quantitative by observing that shrinking the x integration interval in Eq. (10) to , is small in accordance with the prediction of the / πEFT counting: It is also in a very good agreement with the existing dispersive χ ET-based evaluations [29,32]. Despite the smallness of the transverse contribution, we add it to the total inelastic contribution, since it is included in most of the alternative calculations, thus having The uncertainty of the transverse contribution is neglected. Based on the observations above, we conclude that the / π EFT counting used by us works well for the present calculation. We also do not expect any higher-order corrections that would change the pattern that one sees at N3LO; a quantification of this statement follows through the Bayesian procedure in Appendix A.
In Table 4, we compare our E inel 2S result with other recent evaluations. Our result agrees with the recent covariant dispersive calculation [29] as well as with the value quoted in Ref. [32] within the uncertainties. The latter has a slightly larger in magnitude central value. These two results obtain the deuteron response functions at N3LO in the χ ET expansion to calculate E inel 2S from a dispersive integral. The datadriven evaluation of Carlson et al. [24] also uses a dispersive approach, but extracts information on the deuteron response functions from experimental data. It calculates an even larger E inel 2S with a large uncertainty that makes it compatible with all other results. In addition, we compare the results in the point-nucleon limit, where the contributions from the nucleon charge radii are removed (in which case we also omit the contribution of l C0 S 1 ). Our result here is compatible with the earlier N3LO χ ET result [32], as well as with that obtained from the N2LO / π EFT deuteron longitudinal response function in the point-nucleon limit [33].

Single-nucleon effects beyond N3LO
The results of Sects. 3.1 and 3.2 show that single-nucleon contributions generated by the hadron structure, such as the nucleon FFs, are the most important corrections beyond the LO and NLO nuclear-structure effects. They are also potentially the most problematic, since they tend to be enhanced by factors of q 2 compared with the corresponding amplitude with point-like nucleons. For instance, an N4LO correction with two insertions of the nucleon charge radius operator in the LO M L diagrams, shown in Ref. [21, Fig. 7], would be enhanced by a factor of q 4 and would lead to a contribution to E 2S that is divergent at large Q. Another potentially sizeable single-nucleon effect, first appearing also at N4LO, is that of the nucleon polarizabilities. Their inclusion into the deuteron VVCS amplitude also leads to a similar divergent contribution. A / π EFT consideration would therefore introduce fournucleon and two-lepton contact terms at N4LO to regularise the divergence generated by the single-nucleon terms. Those contact terms, as pointed out in Sect. 2.2, limit the predictive powers of / π EFT in the study of the 2γ-exchange corrections to N3LO. In this section, we quantify these hadron struc-ture effects, expected to be the most important ones beyond N3LO, using alternative methods that go beyond the / π EFT expansion.
Starting from the nucleon FF, one alternative that can improve the bad behaviour of the nucleon FF correction would be to insert the full nucleon FFs in the nucleon charge operator vertex, replacing its LO term according to: This procedure obviously represents a departure from the strict / π EFT treatment. It provides, however, a viable workaround and allows one to estimate the effects generated by the higher-order terms in the expansion of the nucleon FFs. It also is routinely used in χ ET calculations of electromagnetic processes in nuclei, since the nucleon FFs do not converge well in a chiral expansion, either, see Refs. [27,29] for recent examples. The specific substitution of Eq. (30), strictly speaking, breaks the electromagnetic gauge invariance. The violating terms are, however, of higher orders than we consider. The modified VVCS amplitudes can be found in Appendix B.
The N3LO / π EFT prediction for the deuteron charge radius, given in Eq. (15), does not change with Eq. (30), as long as we make sure that the parametrisation of the isoscalar nucleon FF agrees with our choice of r 2 0 . We chose the nucleon FF parametrisations from Borah et al. [35], since their slopes are constrained by the nucleon radii used by us: the proton charge radius from μH spectroscopy given in Eq. (1a), and the neutron charge radius [36,37]: The elastic 2γ-exchange correction resulting from inserting the full nucleon FFs can be calculated using Eq. (11) with the re-summed deuteron FFs given in Ref. [21,Eqs. (77) and (78)  . The latter is the inelastic contribution with point-like nucleons (calculated up to N3LO, with the contribution of l C0 S 1 omitted). Values are in meV. To compare with Ref. [32], we sub-tract the subleading O(α 6 log α) Coulomb correction from their "ηless" result. The uncertainty given here for their prediction is obtained using the relative uncertainties of individual error sources from Ref. [34, This amounts to a correction of E hadr, FF 2S = −0.0129 meV with the nucleon FFs from Ref. [35]; the parametrisation of Ref. [38] gives E hadr, FF 2S = −0.0121 meV. In the following, we will adopt This effect is within our N3LO uncertainty estimate; one can also notice that it is significantly larger than a similar difference obtained in a χ ET calculation replacing linearised (expanded in Q 2 ) nucleon FFs by a realistic parametrisation [29]. At the same time, the difference due to the different nucleon FF parametrisations is negligibly small. The replacement of the charge operator by the nucleon FFs in the contributions to E inel,L 2S beyond NLO would also give a negligible effect on the total result.
Coming to the other effect we consider here, that of the nucleon polarizabilities, it consists of two parts, the inelastic and the subtraction hadronic corrections. The first one of the two can be calculated from a dispersive relation, using the empirical deuteron structure functions, at energies starting from the pion production threshold, as done in Ref. [24]: See also Ref. [42] for a similar evaluation. Another method to calculate it, similar to inserting the nucleon FFs in the consideration above, is to apply a rescaling procedure to the inelastic 2γ-exchange effect in μH and the analogous inelastic 2γ-exchange effect between a muon and a neutron (μn), adding them together and correcting for the different atomic wave functions, rescaling the sum by the fac- 2 , as done in Ref. [5] and references therein. Using the single-nucleon values from Ref. [41], we obtain a value of −0.030(2) meV, where we added uncertainties linearly to be conservative. This perfectly agrees with the dispersive evaluation given by Eq. (35), which is an indication that the rescaling procedure works well in this setting.
The second part of the single-nucleon polarizability effect, the subtraction contribution, cannot be directly related to empirical data. It has to be either modelled or predicted from baryon chiral perturbation theory (χ PT). With the rescaling procedure described above, and the covariant χ PT results for the proton VVCS subtraction function [40] and its neutron counterpart, we obtain for the subtraction-function contribution to μD 1 : which agrees well with the value adopted in Ref. [5]: 0.0098(98) meV. As one can see from Table 5, our predictions agree with the dispersive estimates from Ref. [41]. It is also instructive to compare our result for the proton subtraction contribution, 0.0035(26) meV, to predictions in the framework of heavy-baryon χ PT: 0.0042(10) meV [43] and 0.0029(12) meV [44]. The above considerations take into account the most significant higher-order nucleon structure corrections that start to appear at N4LO in the / π EFT expansion. One can notice that each one of the corrections, E hadr, FF 2S = −0.013(1) meV from Eq. (34), and the nucleon polarizability corrections, E hadr, subt 2S + E hadr, inel 2S = −0.019(6) meV from Eqs. (35) and (36), is separately smaller or of the size of the estimated N3LO uncertainty of the inelastic contribution, 0.016 meV, Eq. (26). Their total, however, is about twice as large as that uncertainty. Nevertheless, we expect the higher-order nuclear effects, as well as the relativistic corrections, to be much smaller, and we expect the remaining higher-order effects to be within our N3LO uncertainty estimate. Erring on the side of caution, we refrain from going as far as performing an N4LO adjustment of the uncertainty.

Summary of results
We conclude this section by summarizing our / πEFT predictions of the nuclear-structure effects on the 2S level in μD from the forward 2γ exchange, and including the accompanying electronic VP contributions. At N3LO, we derived the dominant 2γ-exchange effects coming from the elastic deuteron charge FF G C and the non-pole part of the deuteron VVCS amplitude: see Sects. 3.1 and 3.2 for details. The uncertainties have been quantified through the Bayesian error estimate described in Appendix A. As mentioned above, the value of E inel 2S contains the transverse contribution.
In Fig. 3, our / π EFT predictions are compared to datadriven and χ ET results. The disagreement with Carlson et al. [24] for E elastic 2S is due to the deuteron charge FF parametrisation from Ref. [28]. As one can see from Table 2, our prediction is in good agreement with the data-driven approach if the Sick & Trautmann parametrisation [9] is used instead.
Beyond N3LO, we also take into account the singlenucleon effects discussed in Sect. 3.3. They can be split into the nucleon-polarizability contribution, the single-nucleon subtraction-function contribution, and the insertion of the nucleon FFs in the nucleon charge operator vertex of / π EFT. In total, they amount to: On top of the above forward 2γ-exchange effects, there are the electronic VP corrections to the 2γ-exchange, described in Appendix C: (their uncertainty also being negligibly small). In total this adds up to: In Sect. 5, we will discuss all the relevant deuteron-structure effects, including also the Coulomb distortion from the offforward 2γ exchange [5] and the 3γ -exchange effect [7].

Hydrogen-deuterium isotope shift
In this section, we will use the isotope shift between 1S and 2S states in H and D: where h is the Planck constant, to get a prediction for the deuteron charge radius, cf. Eq. (2), and, in turn, determine the LEC l C0 S 1 as given by Eq. (17). The empirically measured value of the isotope shift is very precise [45], To extract from it r d and l C0 S 1 , we will update the theoretical prediction for the isotope shift. Our notation generally follows the work of Jentschura et al. [45]. It is, along with most of the features of the consideration in this section, such as a list of all contributions relevant to the isotope shift, presented in Appendix D. Here, we focus on our / π EFT result for the 2γ-exchange correction to the S-levels in D. The pertinent calculation proceeds analogously to Sect. 3, where the 2γ exchange in μD is evaluated, hence its details are largely omitted.

2γ exchange in deuterium
The longitudinal part of the inelastic contribution to the 2S−1S shift in D is: In the second line, we show our numerical order-by-order result with the LEC l C0 S 1 determined in the following Sect. 4.2. Note that all forward 2γ-exchange contributions scale through the atomic wave function at the origin as 1 /n 3 . Thus, to deduce the shift of the n th S-level in D, one simply has to multiply the isotope shift value by − 8 /7n 3 . The uncertainty of our result is obtained in a simplified way by multiplying the total by (γ /m π ) 4 . This is justified by the smallness of the N2LO and N3LO contributions (with the NLO contribution given by (Z − 1) times the LO result plus a small correction, cf. Table 1 for the case of μD).  Tables 2 and 4 The transverse 2γ-exchange contribution appears to be relatively more important in D than in μD: The uncertainty is obtained here by multiplying the total with (γ /m π ) 3 , where the usual NLO factor of (γ /m π ) 2 is multiplied with another γ /m π to take into account that the transverse amplitude is well reproduced already at NLO [21, Sec. V]. The full N3LO / π EFT prediction for the inelastic contribution to the forward 2γ exchange is then given by The inelastic part, calculated in the same way as done for μD [24], gives [46] ν D 9, hadr. inel = 0.148(11) kHz.
This is in perfect agreement with the result from rescaling the single-nucleon values obtained in Ref. [41]: ν D 9, hadr. inel = 0.145(12) kHz. The subtraction part is calculated by us in the same way as done for μD by rescaling the single-nucleon values from χ PT: The subtraction function contributions found in Ref. [41] tend to be smaller, cf. Table 6.
The off-forward 2γ-exchange correction, known as the Coulomb distortion, can be estimated by rescaling the results for μD presented in Ref. [47].  This has to be compared to ν D 9 = 18.70 (7) kHz used in Ref. [48] and based on Ref. [49].
The N3LO / π EFT prediction for the elastic contribution to the 2S−1S shift in D is where the uncertainty is estimated as above for ν D 9,L . This is slightly bigger than the pure Friar-radius contribution appearing in Ref. [48], which gives ν D (b) = 0.507 kHz. Adding all 2γ-exchange corrections to the 2S−1S transition in D together, we find One can see that the elastic and inelastic contributions proportional to l C0 S 1 partially cancel each other, making the total slightly less sensitive to the value of the N3LO contact term, similarly to what happens in μD. In any case, the effect of it (assuming the maximal magnitude of l C0 S 1 10 −2 ) is at most 0.01 kHz, which is far smaller than the total uncertainty of the isotope shift. Therefore, the contribution of l C0 S 1 to the 2γ exchange in the isotope shift can be safely neglected (at the current level of precision), and the deuteron charge radius extracted from the isotope shift is a good quantity to determine l C0 S 1 .
4.2 2γ exchange in isotope shift and determination of low-energy constant l C0 S 1 For the isotope shift, we also need the 2γ-exchange correction to the 2S−1S transition in H. For the elastic contribution, we use the results from Ref. [41]: which is in perfect agreement with the Friar-radius contribution ν H (b) = 0.035 kHz appearing in Ref. [48]. For the inelastic contribution, it is important that ν H 9 is consistent with the single-proton contributions entering D through ν 9, hadr. inel and ν 9, hadr. subt . Therefore, we will use the subtraction-function contribution predicted by χ PT, see Table 6, and the inelastic contributions from Ref. [41]: This compares to ν H 9 = 0.061(11) kHz used in Ref. [48] and based on Ref.
[50]. Using instead the subtractionfunction contribution from Ref. [41], we would find: In total, the 2γ-exchange correction to the 2S−1S transition in H amounts to The combined results for the isotope shift are given in Eqs. (D19) and (D23) of Appendix D.
In the appendix, we give a full updated list of all contributions entering the isotope shift, together with a comparison to the values used in Ref. [45]. Besides theoretical updates, e.g., of the VP and recoil contributions, we discuss the impact of refined values for the electron, proton and deuteron masses, and the role of the Rydberg constant. Our final result for the theoretical prediction of the 2S−1S deuterium-hydrogen isotope shift reads: Note that, in the calculation of the 2γ-exchange corrections, we used the value of the proton charge radius r p (μH) published by the CREMA Collaboration, Eq. (1a). This value is consistent with the nucleon FF parametrisations from Ref. [35], used in Sect. 3.3 to estimate the single-nucleon effects beyond N3LO in / π EFT. The proton finite-size corrections to the isotope shift use instead a refined value [51], extracted from the Lamb shift measurement of the CREMA Collaboration [1,2] accounting for the recent updates of the μH theory [52][53][54][55]: The effect of the updated r p value on the 2γ-exchange corrections would be negligibly small compared to the estimated theoretical uncertainties of the latter. The LEC l C0 S 1 is small (again, a reasonable estimate of its maximal magnitude being 10 −2 ). It is therefore justified to use the N3LO / π EFT prediction for the deuteron radius, given in Eq. (15), as an exact relation to express l C0 S 1 in Eq. (58) through r d . We can then extract r d by comparing our theory prediction and the experimental value for the isotope shift (44): where the error is completely dominated by the theory. Our result for r d is in perfect agreement with the previous extraction in Eq. (2). Setting l C0 S 1 = 0 in Eq. (58) leads to the same result, which proves that the error generated by applying Eq. (15) as an exact relation can indeed be safely neglected. A comparison and consistency check of state-ofthe-art deuteron charge radius extractions from μD, D and the H-D isotope shift can be found in Sect. 6.1. From Eq. (15), we then find: where the uncertainties in the brackets stem from our extracted value of the deuteron radius, the uncertainty of Z = 1.6893(30) [56], and the isoscalar nucleon charge radius r 0 = 0.5586(10) fm, respectively.

Muonic deuterium Lamb shift
In this section, we will extract an empirical value for the 2γ-exchange effects in the μD Lamb shift from the highprecision Lamb shift measurement by the CREMA Collaboration [3] and the deuteron radius determined from the H-D isotope shift. The empirical value will be compared to our / π EFT prediction. A theory compilation for the μD spectrum, including a review of recent theoretical predictions for the 2γ-exchange effects, can be found in Ref. [5]. At the end of this section, we will present an updated theory prediction of the μD Lamb shift, based on our / π EFT prediction for the 2γ exchange, taking into account all recent theory improvements since the publication of Ref. [5].

Empirical 2γ exchange
The theory prediction for the μD Lamb shift reads [5,Eq. (18)]: −6.11025 (28) r d fm Here, the first term is deuteron-radius independent, the next two terms are deuteron-radius dependent, and the last term contains deuteron-structure effects from 2γ exchange. Note that the prefactor in front of the radius-dependent finite-size term also contains radiative corrections, such as the electronic VP corrections partially discussed in Appendix C, see Ref. [5] for details. 2 The empirical value measured by the CREMA collaboration is: With the theory prediction for the Lamb shift in Eq. (62), the empirical value in Eq. (63), and r d (μH & iso) from Eq. (2), one obtains an empirical value for the 2γ-exchange effects in the μD Lamb shift [3]: In the following, we update this value based on the improved hadronic VP [55] and electronic light-by-light scattering contributions [57], as well as r d (μH & iso) from Eq. (60). For the effect of LO and NLO hadronic VP [55], combined with the mixed electronic and muonic VP, as well as the electronic VP loop in the SE correction [52], we use (2P − 2S): 11.64 (32) μeV. This reduces the uncertainty of the old value 11.12(71) μeV (sum of # #12, 13, 14, 30, 31 in Ref. [5, Table  1]), thereby improving the uncertainty of the deuteron-radius independent term by a factor 2. In addition, we include the inelastic 3γ -exchange, calculated for the first time in Ref. [7]. Compared to Eq. (62), the elastic 3γ -exchange contribution (# #r3, r3 in [5, Table 2]) has been removed from the radius-dependent term, so that the sum of elastic and inelastic 3γ -exchange (2P − 2S): 2.19(88)(27) μeV [7], is now listed as an individual term. The updated theory prediction for the Lamb shift in μD then reads [51]:  Table  2] the entry #r8 has been included with a wrong sign. This work −1.7585 (56) In Sect. 3.4, we summarized our / π EFT results for the deuteron-structure effects in the μD Lamb shift originating from the forward 2γ exchange, including the accompanied electronic VP contributions, and compared to other theory predictions. Our final result is given in Eq. (42). For a meaningful comparison to the empirical value for the 2γ-exchange effect, Eq. (66), we need to add effects from off-forward 2γ exchange (the Coulomb distortions). Formally of a subleading O(α 6 ln α), they are, however, numerically important. We use the recommended value from the theory compilation in Ref. [5]: derived from modern deuteron potentials (χ ET potential and AV18 model [58]). This value should be consistent with the / π EFT framework, since the deuteron electric dipole polarizability from / π EFT [21] is in agreement with predictions from the applied deuteron potentials [59]. Combining Eqs. (42) and (67), our final result for the 2γ-exchange structure effects on the 2S-level in μD reads: which is larger than the value accounted for in Ref. [5,Eq. (17)], but agrees with Ref. [6] within errors, cf. Table  7. It is also in agreement with the empirical value, Eq. (66), but more than a factor 3 less precise. Our new theory compilation will be used in Sect. 6.1 to extract r d (μD) from the experimental value for E 2P−2S .
where the uncertainty budget remained the same as in the original extraction from Ref. [3], see Eq. (1b). In addition, we consider the extraction from the measured 2S−1S transition in D [60]: Note that the entering Rydberg constant, R ∞ in Eq. (E4), is strongly driven by r p (μH). The third extraction from the H-D isotope shift and r p (μH) has been presented in Sect. 4.2: All results are shown in Fig. 4, together will older extractions, results from electron-deuteron scattering and the CODATA recommended values. We can see that the spectroscopy of ordinary and muonic hydrogen isotopes, after the recent theory updates, cf. Ref. [6], gives consistent results for the deuteron charge radius.

Proton charge radius
Analogously to the calculation of r d (μH & iso), we can use the isotope shift and r d (μD) to extract the proton charge radius: While a previous extraction along these lines, r p (μD & iso) = 0.8356 (20) fm [3], had been in tension with r p (μH), the result presented here based on the state-of-the-art theory predictions agrees. This, again, nicely shows the consistency between the spectroscopic analyses of ordinary and muonic hydrogen isotopes.

Proton-deuteron squared charge radii difference
Assuming m e M p ∼ M d , we can find an approximation for the nuclear-size correction to the H-D isotope shift, Eq. (D22), which is related to the often quoted difference of squared proton and deuteron charge radii. The best such approximation turns out to be We give here the charge radius difference exactly, based on Eqs. (59) and (60), and use the relation in Eq. (73) only to estimate the uncertainty, which is dominated by the theory of the isotope shift: Using Eq. (73) instead, the central value would decrease by about 1.5 σ : 3.820 13 +78 −31 fm 2 . These results are in good agreement with the difference between charge radii extracted from the Lamb shift in muonic hydrogen isotopes: From Eq. (74), we can see that the larger CODATA '14 recommended value for the proton charge radius, r p = 0.8751(61) fm [10], would impose a larger value for the deuteron radius inconsistent with the μD Lamb shift.

Conclusion and outlook
In this work, we calculated the 2γ-exchange corrections to the S-levels in ordinary and muonic deuterium in the / π EFT framework. The calculation was performed at N3LO, with the only unknown LEC l C0 S 1 appearing at this order extracted using the H-D isotope shift, where the correlation between that LEC and the 2γ-exchange correction is negligible. In addition, we evaluated the contribution of the nucleon structure, i.e., the effect of the nucleon polarizability and of the shape of the nucleon FFs, which are the most important single-nucleon effects beyond N3LO. We also included the accompanying electronic vacuum polarization contributions.
Our predictions for the elastic contribution to the 2γ exchange in μD from / π EFT at N3LO and χ ET at N4LO appear to be several standard deviations larger than the evaluations [24,29] based on the deuteron charge FF parametrisation of Ref. [28], cf. Fig. 3 and Table 2. This suggests that the latter parametrisation does not adequately describe the behaviour of the deuteron charge FF at low virtualities. The correlation between the Friar radius r Fd and the deuteron charge radius r d in / π EFT, cf. Fig. 2, through the LEC l C0 S 1 could serve as a diagnostic criterion for a realistic parametrisation of the deuteron charge FF. We also point out that the / π EFT expression for the deuteron charge FF at N3LO [21, Sec. IV] can be used for an analytic one-parameter fit to the electron-deuteron scattering data in the low-Q 2 range relevant to the planned DRad experiment [31].
Supplementing the μD theory [5] with a few missing electronic VP effects [6] and the inelastic 3γ exchange [7], together with the shift of the elastic contribution found in this work, the past discrepancy between theory and experiment on the size of 2γ-exchange effects, see Table 7, is now completely resolved.
The uncertainty of the theoretical result for the 2γexchange correction was quantified using Bayesian inference. While our N3LO / π EFT prediction was not yet able to improve the theoretical precision, the improved understanding of the elastic contribution is of utmost importance. In addition, by calculating the 2γ-exchange correction to μD and D, we were able to perform a few consistency checks. In particular, we showed that extractions of the deuteron charge radius from the μD Lamb shift, the 2S−1S transition in D and the 2S−1S H-D isotope shift, cf. Fig. 4, are now in excellent agreement. nications and the collaboration on the calculation of the deuteron VVCS amplitude, as well as for reading our manuscript and providing extremely valuable comments. We thank V. Baru and A. Filin for discussing the details of their work and for sharing with us the results of their χET calculation of the deuteron charge form factor. We thank C. Carlson

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: No data is associated with this manuscript. The associated program codes will be available from the authors upon request.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A: Quantification of uncertainty
The uncertainty of an EFT calculation is in many cases, including the present work, dominated by the unknown higher-order terms rather than input parameters. To quantify the uncertainty, we follow the Bayesian approach developed in Refs. [20,63] and references therein. The results and the details specific to our evaluation are presented in this section.
We start with the EFT expansion of a generic observable A in powers of the expansion parameter ξ (which is ξ = γ /m π in / π EFT): where the parameter A 0 sets the scale of A, and c n are the expansion coefficients. The uncertainty of A caused by a truncation at n = k is given by the unknown remainder: The typical prior assumption used in EFT calculations is the naturalness of the expansion coefficients, i.e., c n should be at most O(1). One may refine this assumption using the calculated expansion coefficients as done in, e.g., Ref. [64] that assigns (written out here for an N3LO calculation, k = 3): where A LO , A NLO , etc. are the contributions to A at the respective order. This also implicitly assumes that A k is dominated by its first term, A 0 c k+1 ξ k+1 ; we will employ this assumption in the following (note that it can be relaxed, in particular, in the Bayesian approach [20,63]). Looking at the expansion of E 2S , we note that the expansion of the LSZ factor, (E d ) −1 ∝ 1 + (Z − 1) + 0 + 0 + · · · , see Ref.
[21, Eq. (49)], results in every term appearing at a given order in the four-point function acquiring a factor of Z at the next order, as is also explicitly shown in Table 1. This induces correlations between the coefficients in the expansion. It is therefore natural to slightly modify the expansion for the purpose of quantifying the uncertainty, taking out the known Z factor along with the scale factor A 0 . This amounts to working with the expansion of the (integrals of) the four-point functions, and the simple estimate of Eq. (A3) becomes: where A 0 and c n are now the normalization and the expansion coefficients of the integrals of the four-point function, which can be deduced from Moving to the Bayesian quantification, we use the prior probability density functions (PDFs) introduced in Ref. [20] as Sets A, B, and C. Each one of the sets of priors consists of two PDFs, pr(c) and pr(c n |c), where the first is associated with the scale c typical of the coefficients c n , while the second describes the probability distribution of c n given the scale.
In the leading-term approximation for the uncertainty, the Bayes' theorem with the given prior PDFs gives the PDF of the uncertainty at order k, given the known c 0 , . . . c k : where we specialize to the case of E 2S and the quantities A 0 and c n again pertain to the expansion of the (integrals of) the four-point function (with the Z factored out , are given in Table 8. One can see that the expansion coefficients of the elastic contribution are somewhat larger than one, whereas the inelastic part, as well as the total, have smaller expansion coefficients. There are partial cancellations between the higher-order elastic and inelastic contributions, which, together with the LO total term being about five times larger than its elastic counterpart, suppresses the expansion coefficients of the total 2γ-exchange correction. With these coefficients, we start with the estimate of Eq. (A4), which gives: The uncertainties deduced this way are a priori not Gaussiandistributed quantities (as one sees explicitly from, e.g., the Bayesian PDFs below, cf. Fig. 5), therefore the usual addition of uncertainties in quadrature may not be an adequate way to, e.g., calculate the uncertainty of a sum given the uncertainties of its constituents; we therefore estimate the uncertainty of E sum 2S independently.
Proceeding to the Bayesian estimates, we apply Eq. (A5) with the priors of Ref. [20] to the results of our calculation. The specific parameters that we use for the priors are: The resulting PDFs for A 3 are shown in Fig. 5 One can see that the simple estimate of Eq. (A7) gives a more conservative uncertainty than is obtained from the Bayesian framework, which was also the tendency observed in Ref. [20]. It is also evident that the priors Set A and B result in somewhat bigger uncertainties than Set C; the difference is rather small for the elastic contribution, but is of the order of 25% for the inelastic term and the total 2γ-exchange correction. Seeing that the current uncertainty estimate uses the leading omitted term approximation, one can expect an about 30% change in the uncertainty once higher-order terms are taken into account. The above difference between Sets A and B, on the one hand, and Set C, on the other hand, is thus not unexpectedly big. One can also note that the uncertainty of the elastic term is probably overestimated, since higher-order terms are unlikely to change that term as much as suggested by the projected uncertainty, given that the value of r 2 d is fixed at N3LO. This is also supported by the agreement between the / π EFT and χ ET results for the elastic contribution. We take the more conservative results of Sets A and B as our uncertainty estimate for the calculated 2γ-exchange contributions. Finally, one may notice that, in particular, the relative smallness of the higher-order coefficients in Table 8, might mean that the assignment ξ = 1/3 overestimates the size of the expansion parameter (and thus also the truncation uncertainty). As demonstrated in, e.g., Ref. [20] with χ ET calculations, Bayesian tools could also be applied to quantify the (assigned) value of the expansion parameter. Such a study would ideally involve additional quantities obtained in the same / π EFT framework, and will be presented elsewhere; in this context, see Ref. [65] that studies the uncertainties arising from a different expansion employed in calculations of 2γ-exchange corrections in muonic atoms and ions.

Appendix B: Deuteron VVCS amplitudes with insertion of nucleon form factors
In this appendix, we provide the expressions for the longitudinal deuteron VVCS amplitude at LO and NLO resulting from the procedure outlined in Sect. 3.3, where one inserts the full nucleon FFs. The expressions for the four-point function (see Ref. [21,Sec. II D] for the definition) with the inserted FFs read: Here, the kinematic functions are [21]: and the barred nucleon isoscalar and isovector electric FFs are:

Appendix C: Electronic vacuum polarization corrections to finite-size and polarizability contributions
In this appendix, we consider the one-loop electronic VP given by: with τ e = Q 2 /4m 2 e and m e the electron mass, and its well known corrections to the deuteron structure effects.
We start with the corrections to the O(α 4 ) deuteron radius term to illustrate the approach, before calculating the corrections to the O(α 5 ) 2γ-exchange effect, relevant for this paper. The one-loop electronic VP correction to the deuteron charge radius term, see Fig. 6, is described by the following potential: Note that this is a contribution to the Breit potential [66, Ch. IX, §83], where the retardation effects can be neglected at this order in α, and hence Q 2 is replaced by q 2 . We shall make use of the dispersion relations (DRs) for the VP and the FF: where ffl denotes the principal-value integration. The oneloop expression for the absorptive part of electronic VP reads: The DR for VP is once-subtracted to ensure the correct normalisation of the electromagnetic field. Similarly to ensure the correct normalisation of the deuteron charge, G C (0) = 1, we can use the once-subtracted relation for the charge FF: In first-order perturbation theory to the unperturbed Coulomb wave functions, one finds for the Lamb shift: It is clear that the dominant effect comes from the small-t region in the first integral, which starts from the threshold of e + e − production. Unfortunately, we cannot simply expand (a) (b) (c) Fig. 7 a One-photon exchange with vacuum polarization; b onephoton exchange with finite-size correction; and c elastic and inelastic two-photon exchange G C around 0 before integration, since the integral will eventually diverge. Instead we use again the DR for G C given in Eq. (C3a). We then change the variable t → 4m 2 e u 2 and perform the integration over u. Afterwards, only integrals over t remain, which start from the threshold of hadron (e.g., π + π − ) production t 0 . Assuming that 2m e t 0 ≤ t , we can expand up to O(4m 2 e /t ). Neglecting terms which are suppressed by additional factors of m 2 e , we obtain: with the auxiliary function: at κ = αm r /2m e . Our formula agrees numerically with Ref. [67,Eq. (28)]. In Eq. (C7b), we used the deuteron radius determined through the isotope shift to illustrate the quantitative size of the effect, where the uncertainity is just propagated from the error of the radius in Eq. (60). A similar subleading correction stems from the interference of one-photon exchange potentials with electronic VP, and finite-size corrections, see Fig. 7a, b, respectively. The latter can be approximated with a delta-function potential proportional to the deuteron radius. To calculate this effect at second order in perturbation theory, we need to know the matrix elements of the delta-function and Yukawa-type potentials between the μD Coulomb wave functions: and the energy levels of the Coulomb potential: with n the principal quantum number. For the discrete spectrum, we obtain: with the deuteron radius in fm units. For the continuous spectrum, we apply: and to get: In total, the interference of the one-photon-exchange potentials in Fig. 7a, b amounts to: In addition, there is a correction to the μD atomic wave function that can be calculated at second order in perturbation theory from the interference of the one-photon exchange potential with VP insertion, Eq. (C9), and the forward 2γ-exchange potential: see Fig. 7a, c, respectively. Since the latter is a delta-function potential just like our approximated one-photon exchange potential with finite-size correction, Eq. (C10), the calculation of E (2) VP 2γ 2S proceeds analogously to the calculation of E (2) VP FF 2S above. We therefore present here only the results: The formula agrees numerically with the wave function correction in Refs. [6,Eqs. (17) and (18)] and [54, Table II].
The sum of electronic VP corrections to the 2γ exchange, amounts to: which is about a factor one-and-a-half larger than our error estimate for E fwd 2S . Our result is comparable to the results of Ref. [6,Eq. (19)]: E eVP 2S = −0.0265(3) meV.
nation of the fine structure constant α from a rubidium recoil measurement [76] and the best caesium recoil measurement [77] has no effect on the isotope shift. The next by size Set (ii) of Ref. [45] includes ten contributions with frequency shifts ν j . We checked that we reproduce the numbers for this set from Ref. [45] individually, using the constants and the formalism given in the CODATA 2006 report [48], up to 0.01 kHz. The new CODATA 2018 constants do not have a significant effect here. However, there have been considerable improvements of the theory, as we show in our updated evaluation below (old values in square brackets): 1. One-loop SE and electronic VP: the result changes marginally, to −0.520 kHz. This is within the uncertainty estimate in Eq. (D5). For a discussion of the pure SE contribution to G 60 see Ref. [81]. 3. Three-loop SE, electronic VP, and combined effects: Including C 50 and C 62 from Ref. [82,83], see Ref. [80]: in addition to the leading C 40 term, has no effect at the precision given above. Note that the four-loop QED contribution has been calculated in Ref. [84], but can be neglected at the present level of precision. 4. Salpeter recoil correction: In Ref. [45], pure recoil corrections at first order in the electron-nucleus mass ratio and expanded up to (Zα) 7 log 2 (Zα) −2 in the (Zα) expansion were included. The uncertainty was estimated assuming that the first neglected higher-order term proportional to (Zα) 7 log(Zα) −2 is of natural size. In Ref. [85], it was shown that the previously neglected coefficient D 71 multiplying the (Zα) 7 log(Zα) −2 term is about a factor 16 larger than the coefficient D 72 multiplying the supposedly more important (Zα) 7 log 2 (Zα) −2 term. This resolved a discrepancy between numerical all-order and analytical (Zα)-expansion results. Here, we use the allorder result reported in Ref. [86]. The two values in Eq. (D15) are the recoil corrections for point-like and Gaussian distributed nuclear charges, respectively. The uncertainty is due to the latter contribution of the nuclearfinite size. As suggested in Ref. [85], we estimate that the uncertainty of the dimensionless parameter δ fns P due to the values of the nuclear radii is given by 2δ R /R δ fns P, where for δ R we take the difference between the proton and deuteron radii in the Gauss nuclear model and the experimental values from Eqs. (59) and (2). Furthermore, we include an uncertainty due to approximations made in the calculation of δ fns P (19 × 10 −7 for H, 65 × 10 −7 for D) and the difference between the nuclear models in Ref. [85,86]. 6. Radiative recoil corrections: 4 ν 6 = −5.38 (32) kHz [−5.38 (11) kHz].
As suggested in Ref. [87], the more conservative uncertainty estimate from CODATA is used [10]. 7. Nuclear SE: While the numerical value of the sum of muonic and hadronic VP contributions remains unchanged, it has in fact been improved, including in addition the effect of NLO hadronic VP and the muonic VP correction to the electromagnetic electron vertex [55]. We checked that the scatter between theoretical predictions and experimental extraction of the hadronic VP contribution to the muon anomalous magnetic moment a μ , used as input in Ref. [55], agrees with an updated value based on the new experimental result from Fermilab [88] and the presently recommended Standard Model prediction collected in the Theory Initiative White Paper [89]. It also covers the recent lattice QCD prediction from the BMW Collaboration [90]. The further improvement in precision of the scatter is not relevant at the level of precision required for the isotope shift. 9. Nuclear-polarizability correction: from combining our results in Eqs. (51) and (55). The uncertainty is dominated by ν D 9 . The decrease of the central value is mainly due to our new / π EFT prediction of the inelastic contribution to the 2γ exchange in D. The error estimate in Ref. [45] seems to be too low by about a factor of 4; based on Eqs. (17) and (18) from Ref. [48], one should get 0.08 kHz instead of 0.02 kHz. In addition, we think that the previously used work [49] might have underestimated the uncertainty. The error estimate presented here is a conservative choice, given the smallness of higher-order terms (in particular, the effect of singlenucleon contributions). We also do not take into account the correlation between the proton contributions in H and D, which would lead to a reduction of the uncertainty estimate. 10. Compensation of the Darwin-Foldy term for the deuteron: where we added the above uncertainties quadratically. Finally, the smallest correction comes from Set (iii) in Ref. [45]. In the following, we will split the set and separate the contributions that we deduce directly by scaling the leading non-relativistic nuclear-size correction: Here, we use the recent update from Ref. [87], including the 3γ -exchange effects due to finite nuclear size.    Leading non-relativistic nuclear-size correction f D iv = −1 369 346 r d fm 2 uncertainty budget in Eq. (58), we can safely ignore this correlation in the above extraction of r n . Note that even though the value of l C0 S 1 in Eq. (G3) agrees well with Eq. (16), the corresponding values for r d disagree, since it also depends on r n through r 0 , cf. Eq. (15).