High-accuracy calculation of the deuteron charge and quadrupole form factors in chiral effective field theory
A. A. Filin, D. Möller, V. Baru, E. Epelbaum, H. Krebs, P. Reinert
HHigh-accuracy calculation of the deuteron charge and quadrupole form factors inchiral effective field theory
A. A. Filin, ∗ D. M¨oller, † V. Baru,
2, 3, 4, ‡ E. Epelbaum, § H. Krebs, ¶ and P. Reinert ∗∗ Ruhr-Universit¨at Bochum, Fakult¨at f¨ur Physik und Astronomie,Institut f¨ur Theoretische Physik II, D-44780 Bochum, Germany Helmholtz-Institut f¨ur Strahlen- und Kernphysik and Bethe Centerfor Theoretical Physics, Universit¨at Bonn, D-53115 Bonn, Germany Institute for Theoretical and Experimental Physics NRC “Kurchatov Institute”, Moscow 117218, Russia P.N. Lebedev Physical Institute of the Russian Academy of Sciences, 119991, Leninskiy Prospect 53, Moscow, Russia (Dated: September 21, 2020)We present a comprehensive analysis of the deuteron charge and quadrupole form factors basedon the latest two-nucleon potentials and charge density operators derived in chiral effective fieldtheory. The single- and two-nucleon contributions to the charge density are expressed in termsof the proton and neutron form factors, for which the most up-to-date empirical parametrizationsare employed. By adjusting the fifth-order short-range terms in the two-nucleon charge densityoperator to reproduce the world data on the momentum-transfer dependence of the deuteron chargeand quadrupole form factors, we predict the values of the structure radius and the quadrupolemoment of the deuteron: r str = 1 . +0 . − . fm , Q d = 0 . +0 . − . fm . A comprehensive andsystematic analysis of various sources of uncertainty in our predictions is performed. Following thestrategy advocated in our recent publication Phys. Rev. Lett. , 082501 (2020), we employ theextracted structure radius together with the accurate atomic data for the deuteron-proton mean-square charge radii difference to update the determination of the neutron charge radius, for which wefind: r n = − . +0 . − . fm . Given the observed rapid convergence of the deuteron form factors inthe momentum-transfer range of Q (cid:39) − . − , we argue that this intermediate-energy domainis particularly sensitive to the details of the nucleon form factors and can be used to test differentparametrizations. PACS numbers: 13.75.Cs, 12.39.Fe, 13.40.Ks, 13.40.Gp, 14.20.Dh ∗ arseniy.fi[email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] ∗∗ [email protected] a r X i v : . [ nu c l - t h ] S e p I. INTRODUCTION
Chiral effective field theory (EFT) is becoming a precision tool for analyzing low-energy few-nucleon reactions andnuclear structure [1–4]. The chiral expansion of the nucleon-nucleon (NN) force has been recently pushed to fifthorder (N LO) [5] and even beyond [6]. The last-generation chiral EFT NN potentials of Ref. [7] provide an excellentdescription of the neutron-proton and proton-proton scattering data, which, at the highest considered order, is evenbetter than the one achieved using so-called high-precision phenomenological potentials such as the CD Bonn [8], NijmI, II and Reid93 [9] and AV18 [10] models. The essential feature of these chiral NN forces is the usage of a semi-localregulator [11, 12], see also Refs. [13, 14], which allows one to significantly reduce the amount of finite-cutoff artifacts inthe long-range part of the interaction. For an alternative regularization approach using a non-local cutoff see Ref. [15].The chiral NN potentials of Ref. [7] also provide a clear evidence of the two-pion exchange, which is determined in aparameter-free way by the chiral symmetry of QCD along with the empirical information on pion-nucleon scatteringfrom the recent analysis in the framework of the Roy-Steiner equations [16, 17]. In the most recent work of Ref. [18],the potential of Ref. [7] was updated to include also the charge-independence-breaking and charge-symmetry-breakingNN interactions up through N LO.In parallel with these studies, a simple and universal algorithm for quantifying truncation errors in chiral EFT withoutreliance on cutoff variation was formulated in Ref. [11] and validated in Ref. [12]. This approach has been successfullyapplied to a variety of low-energy hadronic observables, see e.g. Refs. [19–28]. In Refs. [29–32], it was re-interpretedand further scrutinized within a Bayesian approach.These developments provide a solid basis for applications beyond the two-nucleon system and offer highly nontrivialpossibilities to test chiral EFT by pushing the expansion to high orders. In this paper, we focus on the charge andquadrupole elastic form factors (FFs) of the deuteron.The electromagnetic FFs of the deuteron certainly belong to the most extensively studied observables in nuclearphysics, see Refs. [33–35] for review articles. A large variety of theoretical approaches ranging from non-relativisticquantum mechanics to fully covariant models have been applied to this problem since the 1960s, see Ref. [36] for anoverview. The electromagnetic structure of the deuteron has also been investigated in the framework of pionless [37]and chiral [38–44] EFT.In spite of the extensive existing theoretical work, there is a strong motivation to take a fresh look at the deuteronFFs in the framework of chiral EFT. First of all, the calculation of the deuteron charge FF with unprecedentedaccuracy, by employing consistent NN interactions and charge density operators up to the fifth order in the chiralexpansion, provides direct access to the structure radius of the deuteron and through that to the neutron chargeradius, as elaborated in Ref. [45]. Similarly, the quantitative description of the quadrupole FF, supplemented withthe comprehensive error analysis, opens the possibility to extract the quadrupole moment of the deuteron that is knownvery accurately and thus probes our understanding of the nuclear forces and currents. In this context, it is worthmentioning the tendency of modern nuclear interactions derived in chiral EFT to significantly underpredict the radii ofmedium-mass and heavy nuclei, see e.g. [46]. The existing calculations for A ≥
16 systems do, however, not take intoaccount contributions to the three-nucleon force beyond third order of the chiral expansion (N LO), exchange currentsand relativistic corrections and also suffer from uncertainties intrinsic to truncations of the many-body Hilbert space.It is, therefore, of great importance to test the role of these effects in consistent calculations of electromagnetic few-nucleon processes at high orders in chiral EFT along with a careful error analysis. Focusing on the few-nucleon sectorhas an advantage of avoiding potential uncertainties associated with many-body methods. In particular, no additionalsoftening of the interactions by using e.g. Similarity Renormalization Group transformation [47] is necessary for thelight nuclei like H, H, He and He. It is also interesting and important to test the performance and applicabilityrange of the newest high-precision chiral NN potentials of Refs. [7, 18] and the charge density operators by studying themomentum-transfer ( Q ) dependence of the deuteron FFs and their convergence with respect to the chiral expansion.This provides a rather non-trivial test of the applicability range of chiral EFT since the deuteron FFs decrease byseveral orders of magnitude with increasing values of Q . Therefore, a correction to the charge operator that is smallat Q = 0 may, potentially, have a large impact at higher- Q values.In this paper, we perform a detailed analysis of the deuteron charge and quadrupole FFs in chiral EFT. We includeall contributions to the charge-density operator at fourth order (N LO) relative to the leading single-nucleon operatorand take into account the short-range operators at N LO. The strength of the N LO short-range operators is adjustedto obtain the best fits to the experimental data for the deuteron charge and quadrupole FFs. We demonstrate thatboth the single- and two-nucleon charge density operators can be expressed in terms of the nucleon FFs and exploitthis fact in the calculation of the deuteron FFs. This allows us to avoid reliance on the strict chiral expansion for thenucleon FFs by employing the corresponding empirical parametrizations. Since the errors related to the truncationof the chiral expansion are still very small in the momentum range of Q (cid:39) − . − , this intermediate energydomain appears to be particularly sensitive to the nucleon FFs and thus can be used to test the consistency of theemployed up-to-date nucleon FFs with the deuteron FFs.Once the two NN contact terms in the charge density operator are determined from a fit to the world data on thedeuteron FFs, we arrive at a parameter-free prediction for the quantities at Q = 0, namely the structure radiusand the quadrupole moment of the deuteron. It is worth mentioning that the nucleon FFs do not contribute to theextracted deuteron observables at Q = 0. We perform various consistency checks of our theoretical approach anddemonstrate that (i) our results show only a mild residual cutoff dependence; (ii) the results for the deuteron FFs, thestructure radius and the quadrupole moment are basically insensitive to the choice of off-shell parameters enteringthe NN potentials and the charge density operator. However, this is only true as long as the NN potentials and thecharge density are calculated consistently, which implies that the nucleon FFs must be included both in the one andtwo-body charge density operators, as advocated below. Finally, we perform a detailed error analysis of the obtainedresults by addressing various sources of uncertainties.In Ref. [45], we already employed this approach to extract the structure radius from the charge deuteron FF. Here,we provide additional details of the calculation and update the analysis of Ref. [45] in the following aspects: (a) weemploy the latest version of the NN potential from Ref. [18] that includes the relevant isospin breaking corrections,(ii) we carry out a combined analysis of both the charge and quadrupole deuteron FFs, (iii) relying on our Bayesianestimate of the truncation error from the chiral expansion, in the fits to the FF data we extend the momentum rangeto Q ∼ − as compared to Q ∼ − used in Ref. [45].Our paper is organized as follows. In Section II, we discuss a general formalism to calculate the form factors of thedeuteron. Sections III and IV are devoted to the chiral expansion and regularization of the charge density operator. InSection III we also give a short overview of the nucleon FFs used as input in our calculations. Section V deals with thetreatment of the relativistic corrections. Next, the notation for various contributions to the form factors, their chiralorder and relations to the structure radius and the quadrupole moment are specified in Section VI. Our results forthe momentum-transfer dependence of the charge and quadrupole FFs are presented in Section VII. After fixing theshort-range charge density operator from the best fit to the experimental data we extract the values of the deuteronstructure radius, the neutron charge radius and the deuteron quadrupole moment and analyze various sources ofuncertainties. Also, we discuss the convergence of the chiral expansion for both the deuteron FFs and the extractedquantities at Q = 0. The main results of our study are summarized in Section VIII, where we also discuss their impacton the determination of the neutron charge radius using high-accuracy atomic data on the deuteron-proton chargeradius difference. II. FORMALISMA. Elastic electron-deuteron scattering
The kinematics of elastic electron-deuteron scattering is visualized in Fig. 1 (a) and can be defined as d ( P, λ d ) + e − ( p e , ν ) → d ( P (cid:48) , λ (cid:48) d ) + e − ( p (cid:48) e , ν (cid:48) ) , (1)where variables in brackets denote the momentum and spin projection of the corresponding particle. Throughout thiswork, we focus on the one-photon-exchange mechanism, see Fig. 1 (b), which provides a direct relation between theelectron-deuteron scattering observables and the deuteron form factors. Each additional photon exchange is suppressedby one power of the fine-structure constant. Thus, in line with the conclusions of Ref. [48], these corrections willbe neglected below — see Sec. II E for a more detailed discussion. The one-photon-exchange amplitude of elasticelectron-deuteron scattering, see Fig. 1 (b), can be factorized into leptonic and hadronic parts [49]: M = e ¯ u ( p (cid:48) e , ν (cid:48) ) γ µ u ( p e , ν ) 1 k (cid:104) P (cid:48) , λ (cid:48) d | J µ | P, λ d (cid:105) , (2)where e is the magnitude of the electron charge, u and ¯ u are the spinors of the initial and final electrons normalized as¯ u ( p, ν ) u ( p, ν ) = 2 m e with m e being the electron mass, γ µ are the Dirac matrices and k = P (cid:48) − P is the four-momentum ( a ) ( b )= + · · · P, λ d P ′ , λ ′ d p ′ e , ν ′ p e , ν k FIG. 1. Diagrams representing elastic electron-deuteron scattering. Diagram (a) shows a general contribution to the elas-tic electron-deuteron scattering process and the corresponding kinematics. Diagram (b) visualizes the one-photon-exchangecontribution, while the ellipses refer to multi-photon-exchange processes suppressed by powers of the fine structure constant.Single, double and wiggly lines correspond to electrons, deuterons, and photons respectively. of the exchanged photon. For convenience, we define a quantity Q , which is positive in the space-like region, and thecorresponding dimensionless variable η via Q := − k µ k µ = − k = − ( P (cid:48) − P ) ≥ , η := Q m d , (3)where m d = 1 . (cid:104) P (cid:48) , λ (cid:48) d | J µ | P, λ d (cid:105) can be expressed as [33, 51] (cid:104) P (cid:48) , λ (cid:48) d | J µ | P, λ d (cid:105) = − e G ( Q ) ( ξ ∗ ( P (cid:48) , λ (cid:48) d ) · ξ ( P, λ d )) ( P (cid:48) + P ) µ − e G ( Q ) ( ξ µ ( P, λ d ) ( ξ ∗ ( P (cid:48) , λ (cid:48) d ) · k ) − ξ ∗ µ ( P (cid:48) , λ (cid:48) d ) ( ξ ( P, λ d ) · k ))+ e G ( Q ) ( ξ ( P, λ d ) · k ) ( ξ ∗ ( P (cid:48) , λ (cid:48) d ) · k ) ( P (cid:48) + P ) µ m d , (4)where dimensionless, real, Lorentz-scalar functions G ( Q ), G ( Q ), and G ( Q ) parametrize the photon-deuteroninteraction, and the deuteron polarization four-vectors, ξ ( P, λ d ) and ξ ( P (cid:48) , λ (cid:48) d ), satisfy the following constraints ξ ( P, λ d ) · P = 0 , ξ ( P (cid:48) , λ (cid:48) d ) · P (cid:48) = 0 . (5) B. The electromagnetic form factors of the deuteron
In practice, instead of the scalar functions G i ( Q ) from Eq. (4), one usually introduces the deuteron charge, magneticand quadrupole form factors G C ( Q ), G M ( Q ) and G Q ( Q ), respectively, which are related to G i ( Q ) via the followingequations: G C ( Q ) = G ( Q ) + 23 η G Q ( Q ) ,G M ( Q ) = G ( Q ) ,G Q ( Q ) = G ( Q ) − G ( Q ) + (1 + η ) G ( Q ) . (6)At Q = 0, these form factors are normalized according to [33] G C (0) = 1 , G M (0) = m d m p µ d (cid:39) . , G Q (0) = m d Q d (cid:39) . , (7)where G C (0) = 1 corresponds to the electric charge conservation, Q d = (0 . ± . [52, 53] is the deuteronquadrupole moment, µ d = 0 . m p stands for the proton mass. The derivative of G C ( Q ) with respect to Q taken at Q = 0 is related to thedeuteron charge radius, as discussed in Section VI. C. From observables to form factors
Using the one-photon exchange approximation, the unpolarized elastic electron-deuteron differential cross section inthe laboratory frame reads dσd
Ω ( Q , θ ) = dσd Ω (cid:12)(cid:12)(cid:12)(cid:12) NS (cid:2) A ( Q ) + B ( Q ) tan ( θ/ (cid:3) , (8)where a no-structure pointlike cross section, dσd Ω (cid:12)(cid:12) NS , is defined as the product of the Mott differential cross section, σ Mott , multiplied with the recoil factor dσd Ω (cid:12)(cid:12)(cid:12)(cid:12) NS = σ Mott (cid:16) Em d sin ( θ/ (cid:17) , σ Mott = (cid:16) α E (cid:17) cos ( θ/ ( θ/ . Here E is the energy of the incoming electron, θ is the scattering angle of the electron in the laboratory frame and α is the fine-structure constant. The elastic structure functions A and B are related to the deuteron form factors givenin Eq. (6) via A ( Q ) = G ( Q ) + 23 ηG ( Q ) + 89 η G ( Q ) ,B ( Q ) = 43 η (1 + η ) G ( Q ) . (9)While the unpolarized electron-deuteron scattering cross section in Eq. (8) provides access to the magnetic FF via itsrelation to the structure function B ( Q ), it does not allow one to extract the charge and quadrupole FFs individuallyas they contribute to A ( Q ) in a linear combination. A complementary information on these form factors can beextracted from polarization data. In particular, the experimentally measurable tensor analyzing power T ( Q , θ )gives additional relation: −√ (cid:2) A ( Q ) + B ( Q ) tan ( θ/ (cid:3) T ( Q , θ ) =83 ηG C ( Q ) G Q ( Q ) + 89 η G ( Q ) + 13 η (cid:0) η ) tan ( θ/ (cid:1) G ( Q ) . (10)Therefore, all three deuteron FFs can be extracted individually from a combined analysis of the structure functions A ( Q ) and B ( Q ) together with the polarization observable T . D. Experimental data base
In Ref. [55], a rigorous extraction of the charge, quadrupole and magnetic deuteron form factors from the availableworld data for elastic electron-deuteron scattering was performed in the 4-momentum transfer range of Q = 0 − − .This analysis also includes polarization data of Ref. [56] from JLab. In addition, there is one more recent measurementof tensor polarization observables in elastic electron-deuteron scattering from Novosibirsk [57]. Therefore, in whatfollows, we employ the world data for the deuteron form factors extracted in Refs. [55, 57] as experimental inputexcept for the data point for G Q at Q = 2 .
788 fm − given in Table 1 of Ref. [55], for which we believe the uncertaintyhave been misprinted. Indeed, unlike the data point at Q = 2 .
788 fm − shown in Fig. 1 in Ref. [55] (see the squarewith the strongly asymmetric uncertainty), the uncertainty quoted in Table 1 is symmetric and an order of magnitudesmaller than the one shown in the plot. The error for G Q at this energy is also significantly smaller than those forthe other energies within the same experiment.In a recent review article [35], a parametrization of the world data on the deuteron form factors was provided thathas much smaller uncertainties than in the previous extractions. While we do not use this parametrization in our fits,we will use it for the sake of comparison. E. A comment on the two-photon exchange corrections
Unlike the extensive investigations of the two-photon exchange (TPE) contributions to electron-proton scattering,there are very few works focusing on the study of the TPE corrections for the deuteron electromagnetic FFs. Specif-ically, in Ref. [48] a gauge invariant set of diagrams for the TPE corrections to electron-deuteron scattering wasidentified and estimated under certain assumptions for the photon momentum in the loops. As a result, the effectof the TPE on the charge and quadrupole form factors was found to be very small (less than 1%). Meanwhile, intheir previous investigation [58], the authors found an order of magnitude larger effect from TPE on the deuteronFFs when only one subset of diagrams was included. A significant suppression of the TPE corrections in Ref. [48]is therefore presumably related to the restoration of gauge invariance once the complete set of diagrams is included.The enhanced role of TPE effects was also claimed in Ref. [59], which might again be related to the incomplete set ofdiagrams considered in that work. In the current study we, therefore, rely on the conclusions of Ref. [48] and neglectthe TPE contributions. It would be interesting to have a fresh look at this in future studies.
F. Deuteron form factors in the Breit frame
Deuteron form factors are Lorentz-scalars and can be calculated in any frame, but for practical calculations it isconvenient to choose the Breit frame. In the Breit frame, the kinematic variables take the simple form k = (0 , k ) , P = (cid:18) P , − k (cid:19) , P (cid:48) = (cid:18) P , + k (cid:19) , P = (cid:114) m d + k m d (cid:112) η, k = Q , (11)where the direction of the photon momentum k is chosen along the positive z axis. The polarization vectors ofthe incoming and outgoing deuterons in the Breit frame can be derived by boosting the corresponding rest-framepolarization vectors. For the incoming deuteron, one obtains ξ µ ( P, ±
1) = (cid:18) , ∓ √ , − i √ , (cid:19) , ξ µ ( P,
0) = (cid:16) −√ η, , , (cid:112) η (cid:17) , (12)where the second argument of ξ µ denotes the spin projection of the deuteron onto the z -axis. Similarly, the polarizationvector of the outgoing deuteron in the Breit frame reads ξ ∗ µ ( P (cid:48) , ±
1) = (cid:18) , ∓ √ , + i √ , (cid:19) , ξ ∗ µ ( P (cid:48) ,
0) = (cid:16) √ η, , , (cid:112) η (cid:17) , (13)where the sign of the zeroth component of the polarization vector is opposite from that of the incoming deuteron. Asexpected, these definitions of ξ explicitly satisfy the constraints in Eq. (5).To calculate the deuteron FFs, we express them in terms of the matrix elements (cid:104) P (cid:48) , λ (cid:48) d | J µ | P, λ d (cid:105) defined in Eq. (4).First, we simplify Eq. (4) using the relations ξ ∗ ( P (cid:48) , λ (cid:48) d ) · ξ ( P, λ d ) = ( − δ λ (cid:48) d ,λ d + 2 η δ λ (cid:48) d , δ λ d , ) ,ξ ( P, λ d ) · k = ( − m d ) √ η (cid:112) η δ λ d , ,ξ ∗ ( P (cid:48) , λ (cid:48) d ) · k = ( − m d ) √ η (cid:112) η δ λ (cid:48) d , , (14)which can be derived using the explicit form of the deuteron polarization vectors in the Breit frame given in Eqs. (12)and (13). Simplifying the zeroth and three-vector components in Eq. (4) one obtains (cid:10) P (cid:48) , λ (cid:48) d (cid:12)(cid:12) J (cid:12)(cid:12) P, λ d (cid:11) = 2 P (cid:110) G ( Q ) δ λ d ,λ (cid:48) d + 2 ηδ λ d , δ λ (cid:48) d , (cid:0) G ( Q ) − G ( Q ) + (1 + η ) G ( Q ) (cid:1)(cid:111) , (cid:10) P (cid:48) , λ (cid:48) d (cid:12)(cid:12) J i (cid:12)(cid:12) P, λ d (cid:11) = 2 P √ η G ( Q ) (cid:16) ξ i ( P, λ d ) δ λ (cid:48) d , − ξ ∗ i ( P (cid:48) , λ (cid:48) d ) δ λ d , (cid:17) . (15) ( a ) ( b )+ += ( c ) P, λ d P ′ , λ ′ d k FIG. 2. The matrix element (cid:104) P (cid:48) , λ (cid:48) d | J µ | P, λ d (cid:105) written as a sum of single-nucleon contributions (a) and (b) and the two-nucleon contribution (c). Single, double and wiggly lines refer to nucleons, deuteron particles and photons, respectively. Blackdots and the gray rectangle denote the full photon-nucleon interaction vertex and the two-nucleon current operator. Using Eqs. (6), (12) and (13), we finally obtain G C ( Q ) = 13 e P (cid:0) (cid:104) P (cid:48) , | J B | P, (cid:105) + (cid:104) P (cid:48) , | J B | P, (cid:105) + (cid:104) P (cid:48) , − | J B | P, − (cid:105) (cid:1) , (16) G Q ( Q ) = 12 eη P (cid:0) (cid:104) P (cid:48) , | J B | P, (cid:105) − (cid:104) P (cid:48) , | J B | P, (cid:105) (cid:1) , (17) G M ( Q ) = 1 √ ηe P (cid:28) P (cid:48) , (cid:12)(cid:12)(cid:12)(cid:12) J xB + iJ yB √ (cid:12)(cid:12)(cid:12)(cid:12) P, (cid:29) , (18)where J µB = ( J B , J xB , J yB , J zB ) are contravariant components of the four-vector current in the Breit frame. G. Matrix elements of the electromagnetic current
In the Breit frame, the deuteron form factors are expressed in terms of the matrix elements of the electromagneticcurrent convolved with the deuteron wave functions, (cid:104) P (cid:48) , λ (cid:48) d | J µ | P, λ d (cid:105) , according to Eqs. (16)-(18). The matrixelements read 12 P (cid:104) P (cid:48) , λ (cid:48) d | J µB | P, λ d (cid:105) = (cid:90) d l (2 π ) d l (2 π ) ψ † λ (cid:48) d (cid:18) l + k , v B (cid:19) J µB ψ λ d (cid:18) l − k , − v B (cid:19) , (19)where J µB is the four-vector current calculated in the Breit frame, ψ λ is the deuteron wave function with the polarization λ and the deuteron in the final (initial) state moves with the velocity v B ( − v B ) with v B = k / (2 (cid:112) k / m d ) =ˆ k (cid:112) η/ (1 + η ) and the momenta are defined in Eq. (11). This matrix element is visualized in Fig. 2, where diagrams(a) and (b) involve the single-nucleon electromagnetic current while diagram (c) corresponds to the matrix elementof the two-nucleon current.In this paper, we calculate the deuteron FFs in the framework of chiral EFT utilizing an expansion around thenon-relativistic limit and taking into account relativistic corrections as required by power counting. Specifically, westart with the expressions for the single- and two-nucleon charge density operators, whose chiral expansion will besummarized in the next section. Using the deuteron wave functions at the corresponding order in the chiral expansionand employing consistently regularized expressions for the charge density operators in the partial wave basis, wecalculate numerically the corresponding convolution integrals. III. CHIRAL EXPANSION OF THE CHARGE DENSITY OPERATOR
The nuclear electromagnetic charge and current operators have been recently worked out to N LO in chiral EFT byour group using the method of unitary transformation [61–63] and by the JLab-Pisa group employing time-ordered See Refs. [35, 49, 60] for related studies using manifestly covariant approaches. perturbation theory [64–66], see also Ref. [67] for a pioneering study along this line. Following our works on thederivation of the electromagnetic currents [61–63] and nuclear forces [7, 11, 12, 68–72], in this study we employ theWeinberg power counting for the operators constructed in chiral EFT. The hierarchy of the operators is based on theexpansion parameter q ∈ { p/ Λ b , M π / Λ b } with p being a typical soft scale and Λ b ∼ m N M π (with M π for the pionmass) referring to the breakdown scale of the chiral expansion. This implies that the contributions to the charge andcurrent operators appear at orders q − (LO), q − (NLO), q (N LO), q (N LO) and q (N LO). Notice that theJLab-Pisa group employed the counting scheme with m N ∼ Λ b used in the single-nucleon sector, so that their NLOcorrections appear already at order q − . We further emphasize that the expressions for the two-nucleon charge andcurrent densities in Refs. [61–63] and [64–66] do not completely agree with each other. The differences are, however,irrelevant for the calculation of the deuteron charge and quadrupole form factors. For a comprehensive review of theelectroweak currents and a detailed comparison between the two sets of calculations see Ref. [73]. A. Single-nucleon contributions to the charge density operator
At the chiral order we are working, the single-nucleon contributions to the charge density operator in the kinematics N ( p ) + γ ( k ) → N ( p (cid:48) ) take a well-known form (see Refs. [63, 74] and references therein) ρ = e (cid:18) − k m N (cid:19) G E ( k ) + ie G M ( k ) − G E ( k )4 m N ( σ · k × p ) . (20)Here, G E ( k ) and G M ( k ) are the electric and magnetic form factors of the nucleon respectively, and e is the absolutevalue of the electron charge. The single-nucleon form factors can be written in terms of the isospin projectors andthe corresponding form factors of the proton and neutron G E ( k ) = G p E ( k ) 1 + τ G n E ( k ) 1 − τ ,G M ( k ) = G p M ( k ) 1 + τ G n M ( k ) 1 − τ . (21)For convenience, we also introduce the isoscalar nucleon form factors which are relevant for electron-deuteron scattering G S E ( k ) := G p E ( k ) + G n E ( k ) , G S M ( k ) := G p M ( k ) + G n M ( k ) . (22)In order to facilitate the comparison with phenomenological studies, it is also convenient to decompose the singlenucleon charge density from Eq. (20) into ρ = ρ Main1N + ρ DF1N + ρ SO1N , (23)with ρ Main1N = eG E ( k ) , ρ DF1N = e (cid:18) − k m N (cid:19) G E ( k ) , ρ SO1N = ie G M ( k ) − G E ( k )4 m N σ · k × p , (24)where, apart from the main contribution, DF and SO stand for the Darwin-Foldy and spin-orbit contributions,respectively. Terms involving order- O ( m − N ) corrections to the charge density are beyond the accuracy of our study.The chiral expansion of the electromagnetic FFs of the nucleon is well known to converge slowly as they turn outto be dominated by contributions of vector mesons [75, 76], which are not included as explicit degrees of freedom inchiral EFT. To minimize the impact of the slow convergence of the EFT expansion of the nucleon FFs on two-nucleonobservables, the following two approaches can be employed: • Instead of looking at the individual FFs of the deuteron G C and G Q , one calculates the ratio G C /G Q as donee.g. in Refs. [40, 41]. This is advantageous if one can neglect the contribution of the magnetic form factor G M ( k ) in Eq. (20). However, in addition to this, one also needs to assume either that the contributions fromtwo-nucleon charge densities can be neglected altogether or that they scale with G E ( k ) in the same way asthe one-body densities. Then, the quantity G E ( k ) drops out in the ratio G C /G Q . In this study, we showthat two-nucleon charge density operators should indeed be proportional to G E ( k ), see Sec. III C for details.We also note that due to the numerical smallness of the SO contribution, which is the only term proportionalto G M ( k ), considering this ratio may, in practice, indeed provide quite accurate results. On the other hand,formally, this approximation is not valid at the accuracy level of our analysis. • Instead of relying on the strict chiral expansion of the nucleon FFs one can employ empirical parametrizationsextracted from experimental data, as done e.g. in Ref. [42].In this work, we utilize the second approach and use up-to-date parametrizations extracted from experimental dataas will be described in the next section. The uncertainty of our results associated with the single-nucleon FFs will beaddressed in Section VII E 2.
B. Input for nucleon form factors
The electromagnetic form factors of the proton and neutron probe the charge and magnetization distributions of thenucleons via the interaction of electromagnetic currents and have been investigated experimentally for more than 70years using electron scattering — see e.g. Refs. [81, 84–87] for selected review articles.The most recent extraction of the proton form factors was carried out in Refs. [78, 79], where a global analysis of allexisting data was done including the corrections for different normalization of various data, and effects from TPE. Theresults of Ref. [78] are shown in Fig. 3 (left panel) by red bands confined by solid lines. These fits were constrainedat low Q by the latest CODATA-2018 values for the proton charge radii [91] and by the magnetic radii from theParticle Data Group (PDG) [50] while at high Q a power-law falloff was enforced. Another global analysis of theproton data was carried out by the A1 collaboration in Refs. [77, 92], where specific functional form for the formfactors was assumed to fit the world data and no constraints on the proton radii were imposed. Apart from somedifferences in G p M at low Q and very large differences in the estimated uncertainties, the extracted electric andmagnetic form factors of the proton in Refs. [78, 79] and [77] are essentially consistent with each other, cf. red andgray bands in Fig. 3.A determination of the neutron form factors is much more complicated than for the proton, since there are nofree-neutron targets and it is, therefore, necessary to analyze experimental data on nuclear targets like H or He.A reliable extraction of the neutron form factors from such data requires a detailed understanding of the nuclearcorrections (involving nuclear wave functions, final state interaction, meson exchange currents etc.). The results ofthe most up-to-date parametrization of the neutron FFs carried out in Ref. [78] are presented in Fig. 3 (see red bandsbetween solid lines in the right panel).Already in Refs. [93, 94], it was pointed out that analyticity and unitarity put strong constraints on the nucleon FFs.Using the spectral-function-based dispersive approach, the nucleon FFs were obtained in Ref. [80] from a simultaneousfit to the data for all four FFs in both space-like and time-like regions including the constraints from meson-nucleonscattering data, unitarity, and perturbative QCD. The results of this analysis for the so-called “superconvergenceapproach” (SC) are shown as blue bands confined by the dashed lines in Fig. 3. An update of the analysis of Ref. [80]based on the fit to the most recent MAMI data for electron-proton scattering and simultaneously to the world data forthe neutron form factors was made in Ref. [82] and shown in Fig. 3 by black long-dashed lines. Another strategy wasused in the latest dispersive analysis of Ref. [83]. First, the world experimental data on electron proton scattering werecorrected in Ref. [83] for the TPE contributions, which were calculated including the nucleon and ∆(1232)-resonanceintermediate states. Then, the corrected data were fitted using the proton FFs evaluated in the dispersive approach.No updates of the neutron FFs were made. The comparison of the results of the dispersive approach with those fromthe analysis in Refs. [78, 79] reveals that the electric and magnetic proton FFs from Ref. [83] are compatible with The difference between the nucleon form factor parametrizations presented in the original work of Ref. [79] and its update Ref. [78]lies in the value for the proton charge radius used as input. Ref. [78] employs the more recent (CODATA-2018) value consistentwith the measurements from muonic hydrogen Lamb shift [88] as well as with the latest atomic hydrogen measurements of the Rydbergconstant [89] and the Lamb shift [90], while Ref. [79] relies on the larger value for the proton charge radii taken from CODATA-2014 [54].The effect of the proton charge radius on the shape of the proton FFs is relevant only at very low Q (lower than 1 fm − ). At larger Q ,the shape of the proton form factor is strongly constrained by other experimental data. The difference in G p M might be at least partly related to the fact that the world average value for the magnetic radii of the proton [50]used as input in Ref. [79] has some tension with the value extracted by the A1 collaboration in Ref. [92]. � � � � � � � ������������� � [ �� - � ] � � � / � � � � ��������� ●●●● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ● ●● ●● ● ●● ● ● ● ●● � � � � � � � �������������� � [ �� - � ] � � � / � � � � � � � � � ����������������� � [ �� - � ] � � � / ( μ � � � ) ●● ● ● ●● ● ● ●● ● ● ● ● ● ● � � � � � � � ����������������� � [ �� - � ] � � � / ( μ � � � ) FIG. 3. (Color online) The proton (left panel) and neutron (right panel) form factors normalized to the dipole form factor G D ( Q ) = (1 + Q / Λ D ) − with Λ D = 0 .
71 GeV . The gray bands between the dotted lines correspond to the proton form fac-tors extracted in Ref. [77] from a combined fit to all data including polarized ratio measurements. The uncertainty correspondsto the combined statistical and systematic uncertainties taken in quadrature, among which the sensitivity to the functionalform of the spline used in the fits is the largest. The red bands between the solid lines represent the results of the global analysisof the world proton observables and the neutron FFs from Ref. [78], see also Ref. [79] for the published version and text forthe details. The blue bands between the dashed lines show the results of the SC approach of Ref. [80] from a simultaneousdispersive analysis of all four FFs (data are also shown as dots — see Refs. [79–81] for more details) in both the space-like andtime-like regions. The dashed lines show an update of the analysis of Ref. [80] which is based on the fit to the MAMI datafor electron-proton scattering and simultaneously to the world data for the neutron form factors [82]. The dot-dashed linesrepresent the results for the proton FFs extracted using the dispersive approach of Ref. [83] from a global analysis of the worlddata for electron-proton scattering. No errors for the nucleon FFs were given in Refs. [82] and [83]. the band from Ref. [78] at small and intermediate Q , although they visibly deviate from each other at Q larger than4 fm − (cf. dot-dashed curve with the red band). A closer look at the small momentum range, which is particularlysensitive to the proton charge radius, shows a very good agreement between the results of these analyses, see thezoomed plot for G p E in Fig. 3. This is not surprising given that the value for the proton charge radius predicted inRef. [83] is consistent with the latest CODATA-2018 update employed in Ref. [78] as input, see also Ref. [95] for amini-review on the status of the proton radius puzzle.Last but not least, the lattice QCD simulations for the nucleon FFs are already approaching the accuracy compatiblewith the experimental precision. For example, in Ref. [96], the electromagnetic FFs of the nucleon are computedincluding both the connected and disconnected contributions for the pion masses basically at the physical point. Theresulting isoscalar and isovector nucleon FFs were found to overshoot the experimental data by about one standarddeviation which, as proposed in Ref. [96], could be due to small residual excited state contamination. Furthersimulations should help in resolving this issue.In this work, we will employ a set of different parametrizations for the proton and neutron FFs as input to calculatethe deuteron FFs and, in this way, to make a complementary nontrivial test of our understanding of the nucleonFFs. Indeed, since our current study is aimed at a high-accuracy systematic investigation of the nuclear effects up toN LO in chiral EFT, the comparison of the calculated deuteron form factors with data should provide useful insights1 � � � � � � � ������������� � [ �� - � ] ( � � � + � � � ) / � � � ����������� � � � � � � � ������������������� � [ �� - � ] ( � � � + � � � ) / � � FIG. 4. (Color online) Isoscalar nucleon electric (left panel) and magnetic (right panel) form factors normalized to the dipoleform factor G D ( Q ). For remaining notation see Fig. 3. into the consistency of the single nucleon input with the elastic scattering data on the deuteron. Since the up-to-datedispersive results from Refs. [82, 83] are given without uncertainties, we will use the results of Ref. [78, 79] as ourcentral input, while the FFs from Ref. [80] will be employed as a consistency check.Finally, since the deuteron FFs involve only the isoscalar combinations of the nucleon FFs, see Eq. (22), we plot thesecombinations in Fig. 4. Notice further that the charge and quadrupole FFs of the deuteron are sensitive predominantlyto the isoscalar electric FF of the nucleon, while the isoscalar magnetic FF contributes only through the numericallysmall spin-orbit correction. For the isoscalar electric form factor, the dispersive results of Refs. [80, 82] are essentiallyconsistent with each other as well as with those from the analysis [78] at least for Q < ∼ . − . C. Two-nucleon contributions to the charge density operator
The charge density operator is dominated by the LO single-nucleon contribution, while the first two-nucleon (2N)terms appear only at N LO [61, 62]. The dominant contributions to the 2N charge density operator stem from one-loop diagrams of the one-pion exchange (OPE), two-pion exchange and contact types, whose explicit expressions areparameter-free and can be found in Refs. [61, 62]. All these terms are of isovector type and, therefore, do not contributeto the deuteron form factors. In addition to the already-mentioned static (i.e. order-(1 /m N ) ) contributions, one alsohas to consider tree-level one-pion exchange diagrams with a single insertion of the kinetic energy or 1 /m N -correctionsto the leading pion-nucleon vertex. In the two-nucleon kinematics N ( p ) + N ( p ) + γ ( k ) → N ( p (cid:48) ) + N ( p (cid:48) ) (25)with auxiliary three-momenta defined as q = p (cid:48) − p and q = p (cid:48) − p , the isoscalar one-pion exchange chargedensity can be written as [62] ρ π = (1 − β ) eg A F π m N ( τ · τ ) ( σ · k )( σ · q ) q + M π +(2 ¯ β − eg A F π m N ( τ · τ ) ( σ · q )( σ · q )( q · k )( q + M π ) + (1 ↔ , (26)where the dimensionless quantities ¯ β and ¯ β parametrize the unitary ambiguity of the long-range contributions tothe nuclear forces and currents at N LO. The explicit form of the corresponding unitary transformations is givenin Eq. (4.23) of Ref. [62]. Further, g A is the axial-vector coupling constant of the nucleon, F π is the pion decayconstant and (1 ↔
2) stands for a contribution resulting from interchanging the nucleon labels. Notice that the OPEcontribution has also been taken into account in phenomenological studies, where it represents a part of the so-calledmeson-exchange currents, see e.g. [97].2It is important to emphasize that all terms of the OPE charge density in Eq. (26) are proportional to unobserv-able unitary-transformation-parameters ¯ β and ¯ β . The same parameters also appear in the 1 /m N - and 1 /m N -contributions to the two- [98] and three-nucleon forces at N LO [69], see also Ref. [99] for a related discussion. Thisunitary ambiguity reflects the fact that nuclear forces and currents are not directly measurable and, in general, scheme-dependent. In contrast, observable quantities such as e.g. the form factors must, of course, be independent of thechoice of ¯ β , ¯ β and other off-shell parameters. This can only be achieved by using off-shell consistent expressions forthe nuclear forces and currents. In particular, to be consistent with the new semilocal momentum-space regularizedNN potentials of Refs. [7, 18] which we employ to calculate the deuteron wave function (DWF) for our analysis, theso-called minimal nonlocality choice with ¯ β = 1 / , ¯ β = − / β , ¯ β . The residual dependence on these parameters should be of a higher order, whichprovides a useful tool to check consistency of the calculations. In Section VII E 5, we will demonstrate that thedeuteron FFs calculated with different sets of ¯ β , ¯ β yield consistent results.An important consequence of the unitary ambiguity associated with ¯ β and ¯ β is that one can use unitary transfor-mations to reshuffle the contributions to observables between the charge density and DWF. One can even completelyeliminate the isoscalar 2N charge density operator at N LO. As will be shown below, this also holds true for theshort-range corrections at N LO. The expression in Eq. (26) is thus to be understood as the contribution induced byacting with the unitary operator specified in Eq. (4.23) of Ref. [62] on the isoscalar part of the leading single-nucleoncharge density operator ρ Main, LO1N = e , where the electric nucleon FF at leading order (labeled by the superscript LO)was set to unity. Since we do not rely on the chiral expansion of the nucleon FFs in our analysis, it is more consistentand appropriate to define the isoscalar OPE contribution as the one induced by ρ Main1N rather than ρ Main, LO1N , whichgeneralizes the expression in Eq. (26) to ρ π = (1 − β ) G S E ( Q ) eg A F π m N ( τ · τ ) ( σ · k )( σ · q ) q + M π (28)+(2 ¯ β − G S E ( Q ) eg A F π m N ( τ · τ ) ( σ · q )( σ · q )( q · k )( q + M π ) + (1 ↔ . While this expression is equivalent to Eq. (26) up to terms of a higher order, using Eq. (28) ensures that our results forthe deuteron FFs are independent of the parameters ¯ β and ¯ β to a very high degree, as will be explicitly demonstratedin Section VII E 5.Although the pionic contributions to the isoscalar charge density at N LO have not been worked out yet, the com-plete expression for the contact operators at N LO is derived and given in Appendix B. The expression for theantisymmetrized isoscalar contact operators at N LO reads: ρ Cont = 2 e (cid:18) A + B + C (cid:19) σ · σ + 34 1 − τ · τ k + 2 e C − τ · τ (cid:18) ( k · σ )( k · σ ) − k ( σ · σ ) (cid:19) +2 e ( A − B − C ) 1 − σ · σ (cid:18) τ · τ + 34 (cid:19) k , (29)where the first (second) line in Eq. (29) contributes to the isospin-0-to-isospin-0 (isospin-1-to-isospin-1) channel and A , B , and C denote the corresponding LECs. These LECs contribute to the deuteron FFs in two linear combinations A + B + C/ C . The expression in Eq. (29) agrees with the isoscalar part of the result published in Ref. [101], whilethe corresponding isovector terms are different, see Appendix B. Notice further that the contact operator relevant forthe quadrupole moment of the deuteron (the term ∼ C in Eq. (29)) was first derived in Ref. [37]. Notice, however, that the leading isovector contributions to the 2N charge density at N LO cannot be eliminated by means of unitarytransformations [61, 100]. U = e AT + BT + CT , (30)where the anti-Hermitean generators read T = (cid:0) p (cid:48) + p (cid:48) − p − p (cid:1) ,T = (cid:0) p (cid:48) + p (cid:48) − p − p (cid:1) ( σ · σ ) , (31) T = σ · ( p − p + p (cid:48) − p (cid:48) ) σ · ( p (cid:48) − p (cid:48) − p + p ) + σ · ( p (cid:48) − p (cid:48) − p + p ) σ · ( p − p + p (cid:48) − p (cid:48) ) . Here, p i ( p (cid:48) i ) denote the initial and final momenta of the nucleons. However, in Refs. [7, 18], the freedom to performsuch unitary transformations has already been exploited to eliminate the redundant contact interactions contributingto the S and S partial waves and the mixing angle (cid:15) at N LO. Therefore, to be consistent with the choice ofthe off-shell behavior adopted in the NN potentials of Ref. [18], the short-range contributions to the charge densityin Eq. (29) have to be taken into account explicitly. Here, we follow the same procedure as in the case of the OPEcharge density and employ the short-range charge density operator induced by applying the unitary transformationin Eq. (30) to the charge density operator ρ Main1N from Eq. (24): δ ˆ ρ = ˆ U † ˆ ρ Main1N ˆ U − ˆ ρ Main1N (cid:39) (cid:104) ˆ ρ Main1N , A ˆ T + B ˆ T + C ˆ T (cid:105) , (32)where square brackets denote a commutator and ˆ X indicates that the quantity X is to be regarded as an operatorrather than a matrix element with respect to momenta of the nucleons. Evaluating the commutator in the givenkinematics yields the generalization of Eq. (29) for the contact isoscalar charge density ρ Cont = 2 eG S E ( k ) (cid:20) (cid:18) A + B + C (cid:19) σ · σ + 34 1 − τ · τ k + C − τ · τ (cid:18) ( k · σ )( k · σ ) − k ( σ · σ ) (cid:19) +( A − B − C ) 1 − σ · σ (cid:18) τ · τ + 34 (cid:19) k (cid:21) , (33)where the nucleon FF G S E ( k ) coming from ρ Main1N accounts for a non-pointlike nature of the NN γ vertex. The linearcombinations of the LECs ( A + B + C/
3) and C will be determined from the deuteron data as discussed in Section VI C.The combination ( A − B − C ) corresponds to the isospin-1-to-isospin-1 transition and should be determined fromother processes. For the complete expression including isovector terms the reader is referred to Appendix B.Finally, we emphasize that the above expressions do not provide the complete contribution to the 2N charge densityoperator at N LO. It is, however, conceivable that most (if not all) of the corrections of the one- and two-pion exchangerange, which still have to be worked out, are of isovector type and, therefore, do not contribute to the deuteron FFs.We expect that isoscalar long-range contributions at N LO not considered in our study are, to some extent, effectivelytaken into account by the short-range operators for not too high values of the momentum transfer. For the sake ofbrevity, we will refer to all results based on the short-range part of the 2N charge density operator in Eq. (33) asbeing N LO.
IV. REGULARIZATION OF THE CHARGE DENSITY OPERATOR
We now discuss regularization of the charge-density operators introduced in the previous section. The single-nucleoncharge density operator requires no regularization. However, two-nucleon contributions (both OPE and contact) have To eliminate the redundant terms in NN potential, the parameters A , B , C have to be taken formally of the order O ( m N / Λ b ) ratherthan O (( m N / Λ b ) ). This is the reason for the apparent mismatch in the chiral order of the off-shell contact terms in the NN potential(N LO) and the corresponding short-range charge density operators (N LO). LO contribute to the deuteron FFs. Still, it is important for our analysis to employ a properregulator chosen in a way compatible with the NN potentials of Refs. [7, 18]. In particular, since the contributionof the single-nucleon charge density to the deuteron FFs drops off rapidly with increasing values of the momentumtransfer, the calculated FFs at larger Q -values become sensitive to the two-nucleon charge density operator whichdepends on the regulator.We start with the OPE operators given by Eq. (28). These operators contain single and squared pion propagators.The regularization of the contributions with the single pion propagator is defined in Ref. [7] and can be effectivelywritten as a substitution: 1 p + M π → p + M π exp (cid:18) − p + M π Λ (cid:19) , (34)where Λ is a fixed cutoff chosen consistently with the employed NN potential in the range of 400–550 MeV. Apart from the single pion propagator, the OPE charge density, Eq. (28), also contains the pion propagator squared.The prescription for regularizing the squared pion propagator can be obtained from Eq. (34) by taking a derivativewith respect to M π , as done in Ref. [7], which yields1( p + M π ) → (cid:32) p + M π ) + 1Λ ( p + M π ) (cid:33) exp (cid:18) − p + M π Λ (cid:19) . (35)Using the regularization procedure specified above, the regularized expression for the isoscalar part of the OPE chargedensity takes the form ρ π, reg2N = (1 − β ) G S E ( Q ) eg A F π m N ( τ · τ ) ( σ · k )( σ · q ) q + M π exp (cid:18) − q + M π Λ (cid:19) +(2 ¯ β − G S E ( Q ) eg A F π m N ( τ · τ )( σ · q )( σ · q )( q · k ) × (cid:32) q + M π ) + 1Λ ( q + M π ) (cid:33) exp (cid:18) − q + M π Λ (cid:19) + (1 ↔ . (36)As a next step, we consider the regularization of the contact charge density given by Eq. (33). To ensure consistencybetween regularizations of potential and charge density and avoid ambiguity due to the dependence of the chargedensity operator on the photon momentum, we exploit the fact that both the off-shell contact NN potential andthe short-range charge density operators can be generated by the same unitary transformation acting on the kineticenergy term and the single-nucleon charge density, respectively. The contact part of the NN potential is regularizedin Ref. [7] via a non-local Gaussian cutoff V regCont = V Cont exp (cid:32) − ( p (cid:48) − p (cid:48) ) + ( p − p ) (cid:33) . (37) In Ref. [7], also the results for Λ = 350 MeV are given. However, for such a soft cutoff one already observes a substantial amount offinite-regulator artifacts, and the description of NN data deteriorates noticeably. For this reason we do not use this cutoff value in ouranalysis. regularized generators T i T reg i = T i exp (cid:32) − ( p (cid:48) − p (cid:48) ) + ( p − p ) (cid:33) , i = 1 , , . (38)Then, by acting with this unitary transformation on the single-nucleon charge density ρ Main1N from Eq. (24), we obtainthe consistently regularized 2N short-range charge density operator: ρ regCont = 2 eG S E ( k ) (cid:18) ( A + B ( σ · σ )) F (cid:18) p − p , p (cid:48) − p (cid:48) , k (cid:19) + CF (cid:18) p − p , p (cid:48) − p (cid:48) , k (cid:19)(cid:19) , (39)where the functions F and F are defined as F i ( p , p (cid:48) , k ) = E i (cid:18) p − k , p (cid:48) (cid:19) + E i (cid:18) p + k , p (cid:48) (cid:19) + E i (cid:18) p (cid:48) − k , p (cid:19) + E i (cid:18) p (cid:48) + k , p (cid:19) , (40)with E ( p , p (cid:48) ) = (cid:0) p − p (cid:48) (cid:1) exp (cid:18) − p + p (cid:48) Λ (cid:19) ,E ( p , p (cid:48) ) = [( σ · p )( σ · p ) − ( σ · p (cid:48) )( σ · p (cid:48) )] exp (cid:18) − p + p (cid:48) Λ (cid:19) . (41)Note that, similarly to the procedure used for obtaining the contact interactions, the regularized OPE contribution tothe 2N charge density (as given in Eq. (36)) can also be derived by regularizing the long-range unitary transformation(as given in Eq. (4.23) of Ref. [62]) and acting with it on the single-nucleon charge density ρ Main1N from Eq. (24).Equations (36) and (39) provide the final expressions for the 2N charge density operator at N LO and N LO used inthe calculation of the deuteron FFs.
V. RELATIVISTIC CORRECTIONS
Although the deuteron FFs are Lorentz-invariant, the individual ingredients (charge density operators and deuteronwave functions) do depend on the reference frame. At N LO and below, all frame-dependent corrections are irrelevant,but starting from N LO, the relativistic corrections to each ingredient have to be systematically taken into account.Frame dependence of the charge density operator is automatically accounted for by the kinematics, because theoperator is calculated explicitly including all relevant 1 /m N corrections. In this section, we will focus on the relativisticcorrections to the deuteron wave functions stemming from the motion of the initial and final deuterons.The DWF is typically calculated for the deuteron at rest. However, a calculation of the deuteron FFs always involvesat least one moving deuteron. Our calculation is carried out in the Breit frame, where the initial and final deuteronsare moving in opposite directions. To account for this motion, the rest-frame DWF needs to be boosted. To thechiral order we are working (N LO), the DWF boost corrections have to be considered only when calculating theconvolution integrals of the DWF with the leading single-nucleon charge density ρ Main1N from Eq. (24).Since subleading corrections to the single-nucleon charge density as well as the first contributions to the two-nucleoncharge density appear only at N LO, see Sections III and VI C, the corresponding DWF-boost corrections are beyondthe scope of our study. It is reassuring that the relativistic corrections to the OPE charge density operator considered in Ref. [104] were found to have a tinyeffect on the deuteron charge radius and quadrupole moment. /m N , see also Ref. [34] for a review. Alternatively, boosted DWFwas calculated in Refs. [105–109] utilizing the 1 /m N expansion of the generators of the Poincar´e group. This is theapproach we follow in our analysis. For a deuteron moving with the velocity v , the boosted DWF operator has theform [109] ψ ( p , v ) (cid:39) (cid:18) − v (cid:19) (cid:20) −
12 ( v · p )( v · ∇ p ) − i m N v · ( σ − σ ) × p (cid:21) ψ ( p , , (42)where p is the relative momentum of two nucleons, and ψ ( p ,
0) is the rest-frame DWF which is normalized as (cid:90) d p (2 π ) | ψ ( p , | = (cid:90) d p (2 π ) | ψ ( p , v ) | = 1 . (43)Then, to the order we are working, the boost-corrected matrix element (19) evaluated with the leading density ρ Main1N reads 12 P (cid:10) P (cid:48) , λ (cid:48) d (cid:12)(cid:12) ρ Main1N (cid:12)(cid:12)
P, λ d (cid:11) = e G S E ( k ) (cid:90) d p (2 π ) ψ † λ (cid:48) d (cid:18) p + k boosted , (cid:19) ψ λ d (cid:18) p − k boosted , (cid:19) , (44)where k boosted = k (cid:113) − v B = k √ η , (45)and we used the fact that the spin dependent term in Eq. (42) vanishes for spin-1-to-spin-1 transitions relevant forthe deuteron FFs.The term on the rhs of Eq. (44) is related to the length contraction of that part of the relative nucleon momentum inthe deuteron which is parallel to k . As a consequence of this contraction, the matrix element must be evaluated withthe deuteron wave function taken in its rest frame but with the Breit momentum k replaced by k boosted .Finally, we remind the reader on the ambiguity of the relativistic corrections to the NN potential associated with theemployed form of the Schr¨odinger equation. The corrections to the kinetic energy of relative motion of the nucleonsare most easily taken into account by replacing the nonrelativistic expression p /m N with 2 (cid:112) p + m N − m N insteadof using the Taylor expansion, since otherwise the spectrum of the 2N Hamiltonian is unbounded from below. Insteadof solving the corresponding relativistic Schr¨odinger equation, it is more convenient to rewrite it in the equivalentnonrelativistic form as explained in Ref. [99]. This choice was adopted in the Nijmegen partial wave analysis [110]and is made in the chiral NN potentials of Refs. [7, 11, 12, 18]. While rewriting the Schr¨odinger equation does affectthe m − N and m − N contributions to the NN potential, the deuteron wave function remains unchanged so that we candirectly employ the DWF from Refs. [7, 18]. VI. ANATOMY OF THE CALCULATION
In this section, we summarize the analytic expressions for the deuteron charge and quadrupole form factors, as wellas for the charge radius and the quadrupole moment. We discuss the individual contributions to these quantitiesfrom different types of the charge density introduced in the previous sections. We define the structure radius of thedeuteron and argue, following Ref. [45], that a high-accuracy calculation of this quantity along with high-precisionatomic data for the 1S-2S hydrogen-deuterium isotope shift provide access to the neutron charge radius.7
A. The charge form factor and structure radius of the deuteron
The deuteron charge form factor G C can, up to N LO, be written as G C ( Q ) = G MainC ( Q ) + G DFC ( Q ) + G SOC ( Q ) + G BoostC ( Q ) + G π C ( Q ) + G ContC ( Q ) , (46)where G MainC ( Q ), G DFC ( Q ), and G SOC ( Q ) arise from charge densities defined in Eq. (24), G BoostC ( Q ) is a relativisticcorrection due to initial and final deuteron motion, G π C ( Q ) stems from the one-pion-exchange charge density inEq. (36), and G ContC ( Q ) is generated by the contact charge density in Eq. (39). The main contribution G MainC ( Q )can be factorized as G MainC ( Q ) = ( G p E ( Q ) + G n E ( Q )) G matterC ( Q ) , (47)where G p E ( Q ) and G n E ( Q ) are the electric FFs of the proton and neutron, while G matterC ( Q ) is a functional of thedeuteron wave function.Charge conservation restricts the behavior of the charge form factor at Q = 0. In particular, G C (0) = G MainC (0) = G matterC (0) = G p E (0) = 1, while all other contributions to G C vanish at Q = 0.The deuteron charge radius can be expressed as a derivative of the charge form factor with respect to Q at Q = 0 r d = − ∂G C ( Q ) ∂Q (cid:12)(cid:12)(cid:12)(cid:12) Q =0 . (48)Taking derivatives of all terms in Eq. (46), we get the complete set of contributions to the deuteron charge radius upto N LO r d = r m + r p + r n + r + r + r + r π + r , (49)where the deuteron matter radius r m , the proton charge radius r p and the neutron charge radius r n are defined as: r m = − ∂G matterC ( Q ) ∂Q (cid:12)(cid:12)(cid:12)(cid:12) Q =0 , r p = − ∂G p E ( Q ) ∂Q (cid:12)(cid:12)(cid:12)(cid:12) Q =0 , r n = − ∂G n E ( Q ) ∂Q (cid:12)(cid:12)(cid:12)(cid:12) Q =0 , (50)and the remaining corrections to the deuteron charge radius are calculated as r i = − ∂G i C ( Q ) ∂Q (cid:12)(cid:12)(cid:12)(cid:12) Q =0 , i = { DF , SO , Boost , π, Cont } . (51)Since the r DF -term and the charge radii of the individual nucleons are not related to the two-body dynamics of thedeuteron, they can be conveniently subtracted from the deuteron charge radius. The resulting quantity is referred toas the deuteron structure radius and is defined as (see, e.g. Ref [111]) r = r d − ( r p + r n + r ) . (52)The deuteron-proton mean-square charge radii difference r d − r p in Eq. (52) can be extracted experimentally with anextremely high precision from spectroscopic measurements of the 1S-2S hydrogen-deuterium isotope shift [111]. Inparticular, a series of very precise measurements of the 1S-2S isotope shift, accompanied with an accurate theoreticalQED analysis (see Ref. [112] for the latest update up through O ( α )), resulted in the extraction of the deuteron-protonmean-square charge radii difference [111] r d − r p = 3 . . (53)8Due to its high accuracy, this difference provides a tight link between r d and r p and thus is important in connectionwith the light nuclear charge radius puzzle. For many years, the values for r p extracted from electron and muonexperiments showed more than a 5 σ discrepancy [113]. The very recent atomic hydrogen measurements [89, 90],however, claim consistency with the analogous muonic hydrogen experiments. The recommended value for the protonroot-mean-square charge radius has been changed to r p = 0 . σ larger than the spectroscopic measurement on the muonic deuterium [114] butstill 2.9 σ smaller than the r d value from electronic deuterium spectroscopy [115].As follows from Eq. (52), the deuteron-proton charge radii difference from Eq. (53) allows one to extract the difference r − r n to a very high accuracy. The neutron charge radius can be deduced from measurements of the coherentneutron-electron scattering length extracted from data on neutron scattering off Pb,
Bi and other heavy atoms.The value for the neutron charge radius quoted by the PDG is r n = − . , where the estimated error wasincreased by a scaling factor of 1 . r n = − . , which is consistent with the PDG result,was employed based on the measurement on Pb from Ref. [116]. Using this neutron radius and Eq. (53) for thedeuteron-proton charge radii difference, the value of r str = 1 . Pb and
Bi quoted in the mostrecent investigation of Ref. [116] differ from each other by 0.0090 fm , which is much larger than even the increaseduncertainty given by the PDG. Therefore, a different logical chain was adopted in Ref. [45], namely, (a) by employingthe nuclear forces and currents derived up through fifth order in chiral EFT, a very accurate determination of r str isbecoming possible based on the analysis of the deuteron charge form factor; (b) by using the predicted value for thedeuteron structure radius together with the atomic data for the deuteron-proton charge radii difference, the chargeradius of the neutron was for the first time extracted from light nuclei. In this investigation, we follow the same logic toupdate the analysis of Ref. [45]. In particular, we employ the updated NN potentials which include isospin breakingcorrections up through N LO and provide a statistically perfect description of neutron-proton and proton-protonscattering data up to the pion production threshold [18] to extract the structure radius from a combined analysis ofthe charge and quadrupole deuteron FFs in the range of momentum transfer up to Q = 6 fm − . Then, we updatethe value for the neutron charge radius, see Sec.VII C for the results. B. The quadrupole form factor and quadrupole moment of the deuteron
Deuteron quadrupole form factor can be decomposed in the same way as the charge form factor, namely: G Q ( Q ) = G MainQ ( Q ) + G DFQ ( Q ) + G SOQ ( Q ) + G BoostQ ( Q ) + G π Q ( Q ) + G ContQ ( Q ) , (54)where the individual terms originate from different charge-density contributions in full analogy with Eq. (46). Thedeuteron quadrupole moment is defined as the value of the quadrupole form factor at Q = 0, namely Q d = 1 m d G Q (0) . (55)Taking the Q = 0 limit in Eq. (54) yields the individual contributions to the deuteron quadrupole moment, whichread Q d = Q + Q SO + Q Boost + Q π + Q Cont , (56)where we used the fact that G DFQ (0) = 0 and defined the individual terms as Q = 1 m d G MainQ (0) , Q i = 1 m d G i Q (0) , i = { SO , Boost , π, Cont } . (57)The analytic expressions for various contributions to the deuteron charge and quadrupole form factors as well as tothe structure radius and the quadrupole moment are collected in Appendix A.9 C. Calculational setup
The deuteron FFs at different chiral orders are calculated as follows: • LO:The main contribution to the single-nucleon charge density ρ Main1N in Eq. (24) is convoluted with the LO deuteronwave function. • NLO:Same as LO but using the NLO deuteron wave function. • N LO:Same as LO but using the N LO deuteron wave function. • N LO:The contributions ρ Main1N , ρ DF1N and ρ SO1N from Eq. (24) to the single-nucleon charge density and the OPE contri-bution in Eq. (36) are convoluted with the N LO deuteron wave functions; the relativistic boost corrections tothe single-nucleon contributions are calculated as explained in Section V. • N LO:Same as N LO but using the N LO + deuteron wave function and including the 2N short-range charge densityoperators from Eq. (39).Unless specified otherwise, all results presented below are based on the semilocal momentum-space NN potentials ofRef. [7], updated to incorporate a more complete treatment of isospin-breaking corrections [18]. In particular, theupdated potentials take into account the charge dependence of the pion-nucleon coupling constants. The determinationof the pion-nucleon coupling constants from NN data in Ref. [18] leads to the average value of g πN , which is about 1%larger than the one employed in Ref. [7], and the resulting change in the deuteron wave function leads to a visible effecton the quadrupole FF of the deuteron at higher Q -values. Clearly, in all cases, the same cutoff value chosen from theset Λ = { , , , } MeV is used in the regularized NN potential and in the 2N charge density. For single-nucleon FFs, we employ the most up-to-date parametrization by Ye et al. [78] for our central results. We propagatethe uncertainty in the determination of these FFs to estimate its impact on the deuteron FFs in Section VII E 2. Inthe same section we also consider the impact of using different single-nucleon FFs parametrizations.It remains to specify the values of the various parameters used in the expressions for the 2N charge density operatorin Eqs. (36) and (39). Following Refs. [7, 18], we employ the value of g A = 1 .
29 for the effective axial-vector couplingconstant, which accounts for the Goldberger-Treiman discrepancy, F π = 92 . m N = 2 m p m n / ( m p + m n ) = 938 .
918 MeV for the nucleon mass and M π = (2 M π ± + M π ) / .
03 MeV for thepion mass. Notice that the expressions for the 2N charge density are taken in the isospin limit as the correspondingisospin-breaking corrections start to contribute at N LO, which is beyond the accuracy of our analysis. Finally,the two linear combinations of LECs entering the short-range part of the 2N charge density operator at N LO aredetermined from the best combined fit to the experimental data on the momentum-transfer dependence of the chargeand quadrupole deuteron FFs as described in Section VII. This then allows us to make a parameter-free predictionfor the structure radius and the quadrupole moment of the deuteron.
VII. RESULTS FOR CHARGE AND QUADRUPOLE DEUTERON FORM FACTORS
In this section, we present our results for the deuteron charge and quadrupole form factors. We fix two LECs appearingin the N LO contact charge density by fitting the calculated FFs, G thC ( Q ) and G thQ ( Q ), to the corresponding worldexperimental data for Q < − . Using the LECs extracted from the best fit, we predict the structure radius andthe quadrupole moment of the deuteron. Following Ref. [45], we use the predicted structure radius to extract theneutron charge radius from the precisely measured deuteron-proton charge-radii difference. We provide a detailedanalysis of various uncertainties, discuss several important consistency checks, and discuss the role of the individualcontributions to the charge and quadrupole deuteron form factors.0 A. Fitting procedure
The values of the LECs appearing in the N LO contact charge density of Eq. (39) are determined from a χ -fitof our theoretical predictions for G thC ( Q ) and G thQ ( Q ) to the experimental data. The analytic expressions for theindividual contributions to G thC ( Q ) and G thQ ( Q ) are given in Appendix A, and the experimental data set used in thefit is described in Section II D. In the infinite cutoff limit, G thC ( Q ) depends only on one combination of the LECs,namely A + B + C/
3, while G thQ ( Q ) depends only on the LEC C . Once the regularization is applied, both G thC ( Q )and G thQ ( Q ) in general depend on the two mentioned linear combinations of the LECs, see Eqs. (A28) and (A32) inAppendix A. The function χ ( A + B + C/ C ) to be minimized is defined as follows χ = (cid:88) i (cid:0) G thC ( Q i ; A + B + C/ C ) − G expC ( Q i ) (cid:1) ∆ G C ( Q i ) + (cid:88) i (cid:0) G thQ ( Q i ; A + B + C/ C ) − G expQ ( Q i ) (cid:1) ∆ G Q ( Q i ) , (58)where { Q i } are the set of momenta, for which experimental data are available, and the summations are performed for Q i below Q max = 6 fm − . The intrinsic systematic uncertainty related to the choice of Q max will be discussed below.Following Refs. [31, 118], the uncertainties ∆ G C ( Q i ) and ∆ G Q ( Q i ) in χ include, apart from the experimental errors,also theoretical uncertainties added in quadrature∆ G X ( Q i ) = ∆ G expX ( Q i ) + ∆ G th,truncX ( Q i ) + ∆ G th,nuclFFX ( Q i ) , ( X = C and Q) . (59)In this way, we take into account uncertainties from the truncation of the chiral expansion and from the parametrizationof the nucleon form factors. As the expansion parameter in chiral EFT increases with the momentum transfer, thetruncation errors also grow with Q , as discussed in Section VII E 1. Thus the inclusion of the truncation errors directlyin the objective function allows us to use the deuteron data in a larger range of Q , namely up to Q max = 6 fm − andeven higher. The uncertainty related to the parametrization of the nucleon FFs is yet another source of the theoreticaluncertainty which we include directly in the fit, see Section VII E 2 for details. Other kinds of uncertainties such asthe ones associated with the choice of Q max and with the π N and 2N LECs used in the NN potential are estimatedseparately and discussed below.Our central fit is performed for the cutoff Λ = 500 MeV and Q max = 6 fm − . Assuming that the experimental datapoints are independent , the resulting χ and χ / d.o.f. values are χ = 15 . , χ / d.o.f. = 0 . . (60)The low value of χ / d.o.f. may signal an overestimation of the truncation errors, but it can also be caused byneglecting correlations when estimating truncation errors at similar values of the momentum transfer. The value of χ / d.o.f., therefore, does not allow for a straightforward statistical interpretation. The obtained values of the tworelevant linear combinations of the LECs read A + B + C − ±
64) GeV − (cid:39) ( − . ± .
15) 1 F π Λ b ,C = ( − ±
35) GeV − (cid:39) ( − . ± .
08) 1 F π Λ b , (61)where the error corresponds to the 1 σ deviation of the χ and Λ b = 650 MeV refers to the breakdown scale of thechiral expansion, see Sec.VII E 1 for a discussion. Notice that the both linear combinations of the LECs come out ofa natural size, see the second equalities in Eq. (61). This is an important consistency check of our calculations, whichis also fulfilled for the contact interactions entering the employed NN potentials, see Fig. 7 of Ref. [3]. Finally, thecorrelation matrix for A + B + C/ C reads ρ = (cid:18) − . − . (cid:19) . (62) Note that for the number of degrees of freedom we take just the number of data points minus the number of free parameters. Correlationsbetween data points are neglected. ���������������� | � � ( � ) | ������������ � � ( � ) �� � � �� � � ���������� � � � � � ���������� � [ �� - � ] � � ( � ) � � � � ���������� � [ �� - � ] � � ( � ) �� � � �� FIG. 5. (Color online) Left panel: the deuteron charge and quadrupole FFs calculated at N LO for the cutoff choice ofΛ = 500 MeV (solid red lines) along with the estimated truncation error (68% degree-of-belief) shown by the light-shadedband. Bands between dashed (red) lines correspond to a 1 σ error in the determination of the two short-range contributionsat N LO. Right panel: the same form factors divided by the scaling functions as defined in Eq. (63) and (64). Open violetcircles and green triangles are experimental data from Refs. [57] and [55], respectively. Black solid circles correspond to theparametrization of the deuteron FFs from Refs. [35, 119].
B. Results for the deuteron form factors
The results for the deuteron charge and quadrupole FFs from the best fit to data up to Q = 6 fm − , evaluated forthe cutoff Λ = 500 MeV, are visualized in Fig. 5 together with the N LO truncation errors and statistical uncertaintyof the LEC’s in ρ regCont from Eq. (39). The plot contains two theoretical uncertainty bands: the light-shaded bandstands for the estimated truncation error corresponding to the 68% degree-of-belief interval, while the band betweenlong-dashed (red) lines corresponds to a 1 σ error in the determination of the two short-range contributions at N LO.In principle, these two uncertainty bands are not fully independent since the truncation error is also included in theestimate of the 1 σ error for the LECs in the charge density operator as discussed in previous Section. In this way,however, the truncation error is estimated more conservatively.Since the variation of the FFs at small Q -values is difficult to see on the logarithmic scale, we also plot the rescaledFFs using a linear scale in the right panels of Fig. 5. Specifically, following Ref. [35], we define the rescaled chargeand quadrupole FFs via G scaledC ( Q ) = G C ( Q ) (cid:32) (cid:88) i =0 a i exp( − b i Q ) (cid:33) − , (63)2 TABLE I. Deuteron structure radius squared and deuteron quadrupole moment predicted at N LO in χ EFT (2nd column)and the individual contributions to the corresponding uncertainties from the truncation of the chiral expansion (3rd column),the statistical error in the short-range charge density operator extracted from G C ( Q ) (4th column), the statistical uncertaintyin π N LECs from the Roy-Steiner analysis (RSA) of Ref. [16, 17] propagated through the variation of the deuteron wavefunctions (5th column), the statistical uncertainty in 2N LECs and π N coupling constants f i from the analysis of the 2N dataof Ref. [7, 18] (6th column) and the choice of the maximal energy in the fit (7th column). The total uncertainties evaluated asa sum of presented uncertainties in quadrature are quoted in the 8th column.central truncation ρ regCont π N LECs RSA 2N LECs and f i Q -range total r [fm ] 3 . ± . ± . ± . ± . +0 . − . . − . Q d [fm ] 0 . ± . ± . ± . ± . +0 . − . . − . with a = 0 . a = 0 . a = 0 . b = 3 .
149 fm , b = 1 .
183 fm , b = 0 .
346 fm , b = 0 .
036 fm and a = 1 − a − a − a and G scaledQ ( Q ) = G Q ( Q ) m d Q d (cid:32) (cid:88) i =0 a i exp( − b i Q ) (cid:33) − , (64)with Q d = 0 . , a = 0 . a = 0 . a = 0 . b = 1 .
483 fm , b = 0 .
475 fm , b = 0 .
222 fm , b = 0 .
085 fm and a = 1 − a − a − a . In these plots, along with the comparison of our theoretical results with theexperimental data, we also show the results of the parametrization of the deuteron FFs provided in Refs. [35, 119].While the results for G thC ( Q ) and G thQ ( Q ) are generally quite consistent with this parametrization within errors, amore close look in G scaledC ( Q ) reveals a discrepancy in the range of intermediate Q ’s from 1 fm − to 2 fm − , where theuncertainty from the chiral expansion is still very small. Meanwhile, as will be discussed in Sec. VII E 2, this range ofthe transferred momentum is especially sensitive to the choice of a parametrization of the nucleon FFs. In particular,the inclusion of the uncertainty for the parametrization from Ref. [78] results in the reduction of the discrepancy withRefs. [35, 119]. Nevertheless, the shape of G thC ( Q ) in the range of Q ’s from 1 fm − to 3.5 fm − appears to changemore rapidly as compared to the parametrization by Sick et al. C. Prediction for structure radius and quadrupole moment. Extraction of neutron charge radius.
Using the fitted values of the LECs from Eq. (61) and the theoretical expressions for r str and Q d collected in Ap-pendix A, we make a parameter-free prediction for the deuteron structure radius and the quadrupole moment, whichread r str = 1 . +0 . − . fm , Q d = 0 . +0 . − . fm , (65)where the uncertainties are obtained as a sum of all individual uncertainties given in Table I taken in quadrature, seeSec. VII E for discussion. As advocated in Ref. [45], the knowledge of the deuteron structure radius provides access tothe neutron charge radius, which measures the charge distribution inside the neutron. Using Eqs. (52), (53) and (65),we find r n = − . +0 . − . fm , (66)which is consistent with our previous determination in Ref. [45]. In Section VII D we discuss some differences betweenthe current result and the result of Ref. [45]. D. Comparison to PRL 124, 082501 (2020) (Ref. [45])
While this study is performed along the lines with Ref. [45], there are several updates incorporated in the currentanalysis. These updates can be summarized as follows: (i) the updated SMS potentials of Ref. [18] that include isospin3breaking corrections are employed to calculate the deuteron wave functions; (ii) we now simultaneously fit two linearcombinations of the LECs and use data for both the charge and quadrupole FFs; (iii) our central result is based onthe fit to data up to Q max = 6 fm − ; (iv) statistical uncertainty of the 2N LECs in the NN potential is propagated ina more reliable way.The small difference in the predicted value for the deuteron structure radius and, consequently, also for the neutroncharge radius as compared to Ref. [45] is largely caused by increasing the fitting range up to Q max = 6 fm − . Forsuch value of Q max , both r and Q d are basically saturated with Q max , that is they do not show any significantdeviations in their magnitudes when Q max is increased further. To estimate the error related with the Q max dependenceconservatively, we vary Q max from 3 fm − to 7 fm − . The resulting uncertainties are shown in Table I. The “saturation”of r and Q d above Q max = 6 fm − also explains the asymmetry of the Q max related uncertainties.In addition, we want to make a remark about a finite-cutoff effect, which was neglected in Ref. [45]. In the infinitecutoff limit, G thC ( Q ) at N LO depends only on one linear combination of the LECs, namely A + B + C/
3. On theother hand, for a finite cutoff, both combinations of the LECs A + B + C/ C contribute to both G thC and G thQ ,which, therefore, can be written schematically as G thC ( Q ) = G thC,1 ( Q ) + ( A + B + C/ G thC,2 ( Q ) + C G thC,3 ( Q ) , (67) G thQ ( Q ) = G thQ,1 ( Q ) + ( A + B + C/ G thQ,2 ( Q ) + C G thQ,3 ( Q ) . (68)While the expressions for G thX,2 ( Q ) and G thX,3 ( Q ) with X = C, Q are very different a priori, as can be seen fromAppendix A, in the actual calculations it occurs numerically that the momentum-transfer dependence of G thC,2 ( Q )and G thC,3 ( Q ) (and similarly of G thQ,2 ( Q ) and G thQ,3 ( Q )) is basically identical. In practice, this means that even for afinite cutoff, the FFs in Eqs. (67) and (68), that depend on both linear combinations of the LECs, largely decouple,so that one can study G C independently from G Q . For this reason, in Ref. [45], only the charge FF was considered, inwhich the very last term in Eq. (67) was not included as being redundant. However, because this decoupling is onlyapproximate, in this study we make a combined analysis of both G thC ( Q ) and G thQ ( Q ). By comparing the structureradius extracted in this study with that of Ref. [45], we conclude that they are completely consistent and that theeffect of considering both G C and G Q simultaneously is negligible. On the other hand, since the LECs A + B + C/ C contribute also to other reactions, it is important to extract them individually. This goal can only be achievedif a combined analysis of G thC ( Q ) and G thQ ( Q ) is performed, which allows one to fix A + B + C/ C separately. E. Error analysis
1. Truncation error
We start from the discussion of the chiral expansion for the deuteron form factors which is important for the truncationerror estimate. The convergence pattern of the chiral expansion for the charge and quadrupole deuteron form factorsis shown in Fig. 6 for the cutoff Λ = 500 MeV. Up to and including N LO, the calculation does not involve any freeparameters, while at N LO, two linear combinations of the LEC’s are adjusted to achieve an overall best descriptionof the deuteron FFs in the range of Q -values up to 6 fm − . As a general pattern, the chiral expansion of both formfactors converges quite well.For a given value of the cutoff Λ, truncation errors can be estimated from the convergence pattern of the chiralexpansion using the algorithm formulated in Ref. [11]. This simple approach has, however, a disadvantage of notdirectly providing a statistical interpretation of the estimated errors. We therefore follow here the Bayesian approachdeveloped in Refs. [29–31, 120], which allows one to estimate truncation errors for a given degree-of-belief (DoB)interval. Throughout this analysis, we employ the Bayesian model ¯ C . − specified in Ref. [32] and assume thecharacteristic momentum scale p that determines the expansion parameter q = max (cid:18) p Λ b , M eff π Λ b (cid:19) (69)to be given by | k | / LO, it is easy to see that thedeuteron wave function is being probed at the momentum | k | / | k | , see Ref. [41] for a discussion. The4 ���������������� | � � ( � ) | ������������ � � ( � ) �� � � �� � � � � � � ���������� � [ �� - � ] � � ( � ) � � � � � � ���������������� � [ �� - � ] � � �� � � �� ( � ) FIG. 6. (Color online) Convergence of the chiral expansion for the charge (upper panel) and quadrupole (lower panel) deuteronFFs for the cutoff Λ = 500 MeV. The curves correspond to different chiral orders, namely, black dotted (LO), yellow dashed(NLO), green dot-dashed (N LO), blue long-dashed (N LO) and red solid (N LO). For remaining notation see Fig. 5. quantity M eff π in Eq. (69) serves to model the expansion of few-nucleon observables around the chiral limit, while Λ b denotes the breakdown scale of chiral EFT.In Fig. 5, we show the charge and quadrupole FFs calculated at N LO for the cutoff Λ = 500 MeV along with thetruncation error corresponding to the 68% degree-of-belief (DOB) interval estimated using Eq. (21) from [32] with h = 10, c < = 0 . c > = 10 and assuming Λ b = 650 MeV and M eff π = 200 MeV [121].The truncation errors for the structure radius and the quadrupole moment given in Table I are estimated in exactlythe same way. To make this uncertainty estimate conservatively the truncation error in these quantities is, like inthe deuteron FFs, included twice: (i) by performing the Bayesian analysis for r and Q d explicitly and (ii) throughthe statistical uncertainty in the short-range charge density extracted from the fit to G expC ( Q ) and G expQ ( Q ) usingEqs. (58) and (59). We also provide in Table II the results for the deuteron structure radius and the quadrupolemoment at different orders of the chiral expansion along with the corresponding truncation errors, which show arather natural pattern of convergence for the considered cutoff value of Λ = 500 MeV.
2. Uncertainty from parametrizations of the nucleon form factors
In Fig. 7, we demonstrate the effect of the uncertainties from the nucleon FFs on the deuteron charge and quadrupoleFFs. Our central results, as given by red bands (between solid lines) in Fig. 7, rely on the nucleon FFs extracted5
TABLE II. Convergence pattern of the chiral expansion and the truncation errors for the deuteron structure radius and thequadrupole moment. All results are obtained for the cutoff Λ = 500 MeV and Q max = 6 fm − . Truncation errors for r str arerecalculated from errors estimated for r using the Bayesian approach as described in this Section.LO NLO N LO N LO N LO r [fm ] 3 . ± . . ± .
13 3 . ± .
029 3 . ± .
008 3 . ± . r str [fm] 1 . ± . . ± .
03 1 . ± .
007 1 . ± . . ± . Q d [fm ] 0 . ± .
10 0 . ± .
01 0 . ± .
006 0 . ± . . ± . ���������������� | � � ( � ) | ������������ � � ( � ) �� � � �� � � ���������� � � � � � ���������� � [ �� - � ] � � ( � ) � � � � ���������� � [ �� - � ] � � ( � ) �� � � �� FIG. 7. (Color online) Effect of the uncertainty from various parametrizations of the nucleon form factors (see Fig. 4) onthe deuteron charge (upper panel) and quadrupole (lower panel) FFs. Red bands (between two solid lines) correspond to thenucleon FFs extracted from the analysis of Ref. [78]; blue bands (between the dashed line) based on the nucleon FFs fromRef. [80]. For remaining notation see Fig. 5. from a recent global analysis of electron scattering data on H, H and He targets carried out in Refs. [78, 79] usingthe proton charge radius from CODATA-2018 as input, see Sec. III B for details. The uncertainty from the nucleonFFs, as given in Ref. [78], is included in the statistical uncertainty of our calculation, see Eq. (59).To investigate the sensitivity of the results to parametrizations of the nucleon FFs, we refitted G C ( Q ) and G Q ( Q )using the nucleon FFs from the dispersive analysis of Ref. [80] (the SC approach), where constraints from unitarityand analyticity were included. The results are shown as blue bands between dashed lines in Fig. 7. On the one hand,the results obtained using the parametrizations of Ref. [78] and Ref. [80] are generally consistent with each other asone may already expect from the comparison of the isoscalar nucleon FFs in Fig. 4. On the other hand, the rangewhere the calculated deuteron FFs appear to be especially sensitive to the details of the nucleon FFs corresponds to6the intermediate momentum transfers of Q (cid:39) − . − . In this range, the errors related to the truncation of thechiral expansion are still very small, which can be used to test the consistency of the employed up-to-date nucleonFFs with the deuteron FFs. In the regime of intermediate momenta, G C ( Q ) based on the one-nucleon input fromRef. [80] is systematically lower than that for the nucleon FFs from Ref. [78]. This can be seen from Fig. 7 especiallyif one compares the theoretical results with the parametrization from Refs. [35, 119]: red bands based on the nucleonFFs from Ref. [78] are essentially consistent with this parametrization while the blue bands between dashed lineslie systematically lower. This might be related to the fact that the analysis of Ref. [80] was done before the newhigh-precision data from Mainz [77, 92] have become available. Meanwhile, the updated versions of the dispersiveapproach [82, 83] including the MAMI data produce larger values for the proton electric and magnetic FFs at smalland intermediate momenta and, as shown in Fig. 3, are in a good agreement with the analysis of Ref. [78]. Since theresults of Refs. [82, 83] are given without errors and no updates for a combined dispersive analysis of the proton andneutron FFs was provided in Ref. [83], we refrain from using these results in the current investigation.It is important to emphasize that our results for the structure radius and, therefore, also for the neutron charge radiusare only very weakly sensitive to the details of the nucleon FFs used in the fits. This can be understood as follows.The quality of the fits to the world data for the deuteron charge form factor (at least for Q max ∼ − and higher)increases significantly if the momentum-transfer range around Q ∼ − , where G C becomes small and changes itssign, is well reproduced. Therefore, the contact interaction in the charge density at N LO is adjusted predominantlyto reproduce this area. Meanwhile, the comparison of Figs. 5 and 7 reveals that by far the largest source of theuncertainty at Q ∼ − stems from the truncation of the chiral expansion while the nucleon FFs in this Q -rangehave only a minor impact on the statistical uncertainty. Therefore, the structure radius is insensitive to the choice ofthe parametrization of the nucleon FFs.
3. Statistical uncertainty of the LECs determined from π N and NN data
The chiral SMS NN potential involves two groups of LECs: (i) the π N LECs from the Roy-Steiner analysis ofRef. [16, 17], and (ii) the 2N LECs and π N coupling constants, which are adjusted to achieve the best fit of theneutron-proton and proton-proton scattering data in Ref. [18]. We consider uncertainties coming from each group.To account for the statistical uncertainty of the π N LECs from the Roy-Steiner analysis, we generated a sample of 50N LO + NN potentials with normally distributed π N LECs. Then, the propagation of this uncertainty is performedthrough the variation in the deuteron wave functions. By re-fitting the deuteron FF data we, therefore, extractedthe impact of this uncertainty on r and Q d , as shown in Table I. The resulting uncertainty from these π N LECsappears to be very small.The errors from the statistical uncertainty in the 2N LECs and π N coupling constants extracted in Ref. [18] werealso propagated to r and Q d and the corresponding results are given in Table I. These errors correspond to themaximum deviations from the central values of r m and Q , which are compatible with the variation of the χ inthe range [ χ , χ + 1] for the description of the neutron-proton and proton-proton data as done in Ref. [18].This approach is similar to what was used to estimate the uncertainties of the asymptotic deuteron wave functionnormalization A S and the S NN scattering length in Ref. [7]. Note also that in the present work, the method oferror propagation from 2N LECs is different from what was done in Ref. [45], where a covariance matrix was used. Wefound that the covariance-matrix approach overestimates the corresponding uncertainties for r . For the deuteronquadrupole moment, however, both approaches give very similar error estimates. Q -range dependence As long as the truncation error is included in the uncertainty employed in the fitting procedure, as done in Eq. (59),all data available can, in principle, be included in the fits. This procedure allows us to utilize the deuteron data inthe range of Q up to Q max = 6 fm − and even higher. The effective weight of the data points at higher transferredmomenta is reduced as compared to data points with similar experimental errors at lower Q because the truncationerror increases with growing values of Q . To estimate (conservatively) the error for the extracted deuteron quantitiesrelated with the truncation of the Q -range in the fits, we consider the variation of Q max from 3 fm − to 7 fm − andinclude this error in the uncertainty budget, as shown in Table I. The results for both r and Q d appear to be quitestable to this variation.7 ���������������� | � � ( � ) | ������������ � � ( � ) �� � � �� � � ���������� � � � � � ���������� � [ �� - � ] � � ( � ) � � � � ���������� � [ �� - � ] � � ( � ) �� � � �� FIG. 8. (Color online) Residual cutoff dependence versus the truncation error for the deuteron charge and quadrupole FFs atN LO. Light-shaded blue bands between two solid lines correspond to the cutoff variation in the range of Λ = 400 . . .
550 MeV.For remaining notation see Fig. 5.
5. Consistency checks
We are now in the position to perform several consistency checks of our calculations.As already pointed out in Sec. III C, the two-body charge density from OPE is proportional to unobservable unitary-transformation parameters ¯ β and ¯ β . The observables must be independent of these parameters at least approx-imately, i.e. up to higher order effects. The results presented above are based on the minimal nonlocality choice,Eq. (27), which is consistent with the employed chiral NN potentials of Ref. [18] and also with their predecessors fromRef. [7]. To check the sensitivity of the deuteron FFs to ¯ β and ¯ β , we developed an approximately phase-equivalentversion of the 2N potential using a different choice of the unobservable phases, namely ¯ β = ¯ β = 1 /
2, by re-doingthe fit of NN data using exactly the same protocol as in Ref. [18]. For this particular choice of ¯ β , ¯ β , the OPEcontribution to the charge density vanishes exactly: ρ π = 0. Repeating the fits of the calculated deuteron FFs tothe world data, we find for the central values r = 3 . , Q d = 0 . , (70)which should be compared with the values in Table I. As expected, the dependence on ¯ β and ¯ β for r turns out tobe very small, that is much smaller than the truncation error of the chiral expansion at N LO. For the quadrupolemoment, the dependence on these parameters is also consistent with the truncation error. Note that to achieve theindependence of ¯ β , ¯ β to such a high degree, we found it to be crucial for the nucleon FFs to be included not only8 �� ��� � � �� � � �� � � �������������� ������ ����� � � � � � [ � � � ] �� ��� � � �� � � �� � � ���������������������������������� ������ ����� � � [ � � � ] FIG. 9. (Color online) Convergence of the chiral expansion for the deuteron structure radius squared (left panel) and thequadrupole moment (right panel). Error bars correspond to the truncation errors in a given chiral order. Green circles standfor the results for Λ = 400 MeV, orange squares for Λ = 450 MeV, red diamonds for Λ = 500 MeV, and purple trianglesfor Λ = 550 MeV. The band bounded by solid horizontal gray lines on the right panel correspond to the quadrupole momentextracted in Refs. [52, 53]. The bands bounded by the dashed horizontal red lines on both panels correspond to the total uncertainties of our central results as given in Table I. in the one-body but also in the two-nucleon OPE charge density. This can be understood as follows: as discussedin Sec. III C, the derivation of the two-nucleon charge density operators relies on taking the commutator of theleading one-body charge density with the generators of the unitary transformation. Because the one-body density isproportional to the nucleon FF, the same should also hold for the two-body densities. If one neglects the nucleonFFs in the OPE charge density, a sizable violation of the ¯ β , ¯ β independence would immediately reveal itself in thedeuteron quantities. Specifically, in this case one gets r = 3 . and Q d = 0 . , and one sees that thedifference with the values given in Table I exceeds the truncation error significantly. The effect on these quantitiesof neglecting the nucleon FFs in the short-range two-body charge-density operator at N LO is of basically the samesize. Also, we would like to emphasize that to observe ¯ β , ¯ β independence in a large range of momentum transfersit is important to follow the procedure, as described above: first, construct the phase-equivalent NN potentials forsome choice of ¯ β , ¯ β by fitting NN data and then calculate corresponding deuteron wave functions. If one applies theunitary transformation to the existing wave function, then the unitary equivalence will hold only at small momentumtransfers, as shown in Ref. [122].Since the chiral expansion for the deuteron FFs is expected to converge more rapidly for not too soft values of thecutoffs, our central results are obtained for the cutoff Λ = 500 MeV, for which we also carried out a detailed erroranalysis as described in the previous sections. In Fig. 8, as a consistency check, we confront the cutoff dependence of G C ( Q ) and G Q ( Q ) from the variation of the cutoff from 400 to 550 MeV with the truncation error. We conclude thatfor G C ( Q ), the cutoff dependence lies well within the truncation error, while for the quadrupole FF they are essentiallycompatible with each other except for the region of small Q where the cutoff dependence is a little larger. We remindthe reader, that the truncation error corresponds to the 68% DOB interval. In Fig. 9, we show the convergence patternof the chiral expansion for r and Q d along with the truncation errors and the cutoff dependence. While the resultsfor r converge quite rapidly for any cutoff value, the quadrupole moment, in line with the discussion above, showsa lower rate of convergence. We further emphasize that the uncertainty of our result for the quadrupole moment atthe highest considered order is dominated by the statistical errors in the NN LECs and by the uncertainty associatedwith the choice of the Q -range in the fit. Both of these error sources are considerably larger than the truncationuncertainty.9 TABLE III. Impact of the individual contributions to the charge density operator on the charge form factor of the deuteron G C ( Q ) at N LO for the cutoff Λ = 500 MeV.Q [fm − ] Main SO Darwin Boost 1 π CT Full0 . . . . . . . . . . . − . . − . − . . . . . − . . − . − . . . . . − . . − . − . . . . . − . . − . − . . . . . − . . − . − . . . . . − . . − . − . − . . − . . . . − . − . − . G Q ( Q ) /m d ) at N LO for the cutoff Λ = 500 MeV. The values are given in fm .Q [fm − ] Main SO Darwin Boost 1 π CT Full0 . . − . . . . . . . . − . − . − . . . . . . − . − . − . . . . . . − . − . − . . . . . . − . − . . . . . . . − . − . . . . . . . . − . . . . . . . . − . . . . . F. Role of the individual contributions at N LO The role of the individual charge density contributions to the deuteron charge and quadrupole form factors calculatedat N LO can be seen in Tables III and IV. Unlike the results in Fig. 6, which illustrate the convergence of the chiralexpansion, all contributions in Tables III and IV were evaluated using the deuteron wave function at N LO + . Theresults for the one-pion-exchange charge density were obtained using the minimal nonlocality choice for the parameters¯ β and ¯ β from Eq. (27). For the charge FF, the most important correction beyond the main one stems from theCT contribution which dominates in the whole domain of momenta Q considered apart from the region of small Q ( < ∼ − ), where the Darwin term is equally important. Next in importance are the 1 π and boost correctionswhich, however, cancel each other to a large extent. The contribution from the SO is basically negligible. For thequadrupole FF at low Q ( < ∼ − ), the dominant corrections beyond the main term originate from the 1 π , CTand SO contributions, in the order of their importance, where the first two interfere constructively while the SOis destructive. While the boost correction is negligible for all Q -values, the Darwin term, being negligible at smallmomentum transfers, provides a sizable contribution for Q >
VIII. SUMMARY AND CONCLUSIONS
In spite of the extensive progress in the understanding of the deuteron structure that has been achieved since morethan 50 years, there is still strong motivation to reanalyze the deuteron form factors in the framework of chiral EFT.Being largely governed by the leading-order single-nucleon charge density, the charge and quadrupole deuteron formfactors (FFs) are qualitatively described in most of the calculations reported in the literature (at least in some rangeof the momentum transfer). However, as long as higher-order corrections are concerned, the existing calculations showlack of systematics, consistency and controlled error estimate.In this paper, the deuteron charge and quadrupole form factors are calculated using consistently regularized two-0nucleon potentials and the charge density in chiral effective field theory. This allowed to extract the important staticproperties of the deuteron, namely the structure radius and the quadrupole moment, with unprecedented accuracyand to reliably estimate various sources of uncertainty. Our analysis provides a first step towards the understandingof radii of medium-mass and heavy nuclei, which are currently known to be significantly underpredicted.The novel aspects of our study include: • For the first time, the calculation of the deuteron FFs is pushed beyond N LO, which allows one to reduce theuncertainty from the truncation of chiral expansion and thus to extend the range of momenta considered. Toachieve this goal we (i) employed the most recent two-nucleon potentials up through N LO + [18], which utilizea complete treatment of isospin-breaking effects and provide a statistically perfect description of NN data belowpion production threshold and (ii) implemented the charge density operator at N LO, supplemented with the2N short-range operators at N LO. • Regularization of the charge density operators is carried out consistently with the two-nucleon potential usingthe same unitarity transformations for the charge density operators and the nuclear forces. Specifically, the two-nucleon charge density operators are generated by taking the commutator of the leading one-body charge densitywith the generators of the unitary transformation that incorporate the regulator as discussed in Sec. III C. As aconsistency check, we have demonstrated that the residual cutoff dependence of the deuteron charge FF and theextracted structure radius is much weaker than the error estimated from the truncation of the chiral expansionat N LO. The cutoff dependence of the quadrupole moment (and in general of the quadrupole FF at low Q ) atthe highest considered order is somewhat larger than the estimated truncation error, but still of the same sizeas the total uncertainty of our result. Furthermore, the short-range charge density operators contributing tothe charge and quadrupole FFs of the deuteron come out of a natural size. • Instead of relying on the strict chiral expansion of the nucleon FFs known to converge slowly, we employed themost up-to-date phenomenological parametrizations of experimental data from the global analysis of Refs. [78,79]. The nucleon form factors from the dispersive approach of Ref. [80] have also been used as a consistencycheck. We emphasize that making a reliable calculation of the deuteron FFs requires the inclusion of the nucleonFFs both in the one- and two-nucleon charge density operators, the feature that becomes obvious in the way wegenerate the two-body charge density by means of the unitary transformation. We have verified this conclusionby explicitly checking the insensitivity of our results for the FFs to the choice of unobservable unitary phases¯ β and ¯ β , which holds true to a very high degree of accuracy when keeping the nucleon FFs in the OPE chargedensity. The same conclusion applies when the nucleon FFs are neglected in the contact two-nucleon chargedensity at N LO. • A comprehensive and systematic analysis of various sources of uncertainties in the calculated deuteron FFsis performed. Specifically, we estimated the uncertainty from (i) propagating the statistical errors of the π Nand NN low-energy constants (LECs) entering the two-nucleon potentials (ii) truncation of the chiral expansionevaluated using Bayesian methods (iii) statistical uncertainties in the N LO short-range charge density operators,(iv) employed parametrizations of the nucleon FFs and (v) fixing the maximum value of the momentum transfers Q max in the fits of the short-range charge operators.Pushing the calculation to N LO and using the consistently regularized charge density operators together with thephenomenological nucleon form factors is found to result in a very good description of the deuteron form factors atleast up to Q (cid:39) − . Having adjusted the two short-range operators to achieve the best fit of the world data onthe charge and quadrupole FFs of the deuteron, we predict the deuteron structure radius and quadrupole moment tohave the values of r str = 1 . +0 . − . fm , Q d = 0 . +0 . − . fm . (71)Equipped with this prediction for the structure radius, we employ the high-accuracy data for the hydrogen-deuteriumisotope shift in Eq. (53) to extract the mean-square neutron charge radius, for which we obtain r n = − . +0 . − . fm . (72)This result is consistent with our previous determination in Ref. [45] but deviates by about 1 . σ from the currentvalue r n = − . given by the Particle Data Group [50] and deviates by about 1 . σ from the very recentdetermination r n = − . ± . (stat . ) ± . (syst . ) fm from the collective analysis of the nucleon form factors ofRef. [123].1 ACKNOWLEDGMENTS
We are grateful to U.-G. Meißner for a careful reading of the manuscript and valuable comments and to Z. Ye forproviding us with the unpublished results for the nucleon form factors from Ref. [78]. We also thank H.-W. Hammerfor providing us with the parametrization of the nucleon form factors from Ref. [80] and I. Sick for the parametrizationof the deuteron form factors from Ref. [35]. We are grateful to M. Hoferichter and J. Ruiz de Elvira for the informationon the central values and covariance matrix of the N LO π N LECs from the Roy-Steiner analysis. This work wassupported in part by DFG and NSFC through funds provided to the Sino-German CRC 110 “Symmetries and theEmergence of Structure in QCD” (NSFC Grant No. 11621131001, Grant No. TRR110), the BMBF (Grant No.05P18PCFP1) and the Russian Science Foundation (Grant No. 18-12-00226).
Appendix A: Analytic expressions for the contributions to G C , G Q , r str and Q d . In this appendix we list the analytic expressions for individual contributions to the deuteron charge form factor G C (Eq. (46)), quadrupole form factor G Q (Eq. (54)), structure radius squared r (Eqs. (52), (49)) and quadrupolemoment Q d (Eq. (56)). Results are given in momentum space or in coordinate space depending on which form issimpler for practical calculations. All contributions are grouped according to the charge density operator which theyare obtained from. Specifically, we distinguish the following types of contributions: main (LO) contributions, Darwin-Foldy-type contributions, spin-orbit corrections, deuteron boost corrections, pion-exchange current contributions, andcontact contributions. Results are expressed in terms of deuteron wave functions and single-nucleon form factors.The deuteron WFs are normalized according to: ∞ (cid:90) p (cid:16) u ( p ) + w ( p ) (cid:17) dp = ∞ (cid:90) (cid:16) u ( r ) + w ( r ) (cid:17) dr = 1 , (A1)and the deuteron D -state probability is: P D = ∞ (cid:90) p w ( p ) dp = ∞ (cid:90) w ( r ) dr. (A2)We also introduce the following common combinations of the deuteron wave functions: C ( r ) ≡ u ( r ) + w ( r ) , Q ( r ) ≡ u ( r ) w ( r ) − w ( r ) √ . (A3)Note that all momenta in this section are three-dimensional. For a vector x , we use x to denote x ≡ | x | .
1. Main (LO) contributions
Main contributions stem from the LO charge density operator in Eq. (24); G MainC ( k ) = G S E ( k ) G matterC ( k ) , G MainQ ( k ) = G S E ( k ) G (0)Q ( k ) , (A4)where we introduced the following auxiliary functions: G matterC ( k ) ≡ (cid:90) ∞ C ( r ) j (cid:18) kr (cid:19) dr, G (0)Q ( k ) ≡ √ m d k (cid:90) ∞ Q ( r ) j (cid:18) kr (cid:19) dr. (A5)2The LO contribution to the deuteron structure radius (so called deuteron matter radius) reads: r m = r = 14 ∞ (cid:90) (cid:16) p (cid:16) u (cid:48) ( p ) + w (cid:48) ( p ) (cid:17) + 6 w ( p ) (cid:17) dp = 14 (cid:90) ∞ (cid:16) u ( r ) + w ( r ) (cid:17) r dr. (A6)The LO contribution to the deuteron quadrupole moment reads: Q = ∞ (cid:90) (cid:32) p u (cid:48) ( p ) w (cid:48) ( p )5 √ − p w (cid:48) ( p ) + 3 pw ( p ) u (cid:48) ( p )5 √ − w ( p ) (cid:33) dp = √ (cid:90) ∞ Q ( r ) r dr. (A7)
2. Darwin-Foldy-type of contributions
The Darwin-Foldy-type (DF) of contributions stem from the charge density operator ρ DF1N in Eq. (24). Since the DFcharge density operator differs from the LO operator only by a pre-factor, the resulting DF contributions to the formfactors are trivially related to the LO ones. Specifically, the DF contributions to the deuteron charge and quadrupoleform factors read: G DFC ( k ) = G S E ( k ) (cid:18) − k m N (cid:19) G matterC ( k ) , G DFQ ( k ) = G S E ( k ) (cid:18) − k m N (cid:19) G (0)Q ( k ) . (A8)Deuteron structure radius does, by definition, exclude the Darwin-Foldy contribution, but the deuteron charge radiusreceives a constant correction r DF = 3 / (4 m p ) = 0 . , where m p is a proton mass. Finally, the Darwin-Foldyterm does not contribute to the deuteron quadrupole moment since it is proportional to the photon momentum k ,while the quadrupole moment is defined at k = 0.
3. Spin-orbit contributions
The spin-orbit contributions to the deuteron form factors stemming from the charge density operator ρ SO1N of Eq. (24)read: G SOC ( k ) = (cid:0) G S E ( k ) − G S M ( k ) (cid:1) G angC ( k ) , G SOQ ( k ) = (cid:0) G S E ( k ) − G S M ( k ) (cid:1) G angQ ( k ) , (A9)where G angC ( k ) ≡ m N (cid:90) ∞ w ( r ) r ∂∂r (cid:18) j (cid:18) kr (cid:19)(cid:19) dr = − k m N (cid:90) ∞ w ( r ) r j (cid:18) kr (cid:19) dr, (A10) G angQ ( k ) ≡ ( −
1) 6 √ k m d m N (cid:90) ∞ w ( r ) (cid:18) ∂∂r (cid:18) u ( r ) r (cid:19) − √ r ∂w ( r ) ∂r (cid:19) j (cid:18) kr (cid:19) dr. (A11)The corresponding contributions to the deuteron structure radius and quadrupole moment read: r = − m N (2 µ n + 2 µ p − P D , Q SO = (1 − µ n − µ p ) Q angular , (A12)where µ p and µ n are the magnetic moments of the proton and the neutron, respectively, in units of nuclear magnetons,and Q angular ≡ ( −
1) 6 √ m N (cid:90) ∞ w ( r ) (cid:18) ∂∂r (cid:18) u ( r ) r (cid:19) − √ r ∂w ( r ) ∂r (cid:19) r dr. (A13)3
4. Boost corrections
Corrections to the deuteron form factors which appear due to the motion of initial and final deuterons are discussedin the Section V. The final expressions for the boost corrections to the charge and quadrupole form factors have theform: G BoostC ( k ) = G S E ( k ) G angC ( k ) + G S E ( k ) (cid:0) G matterC ( k ) − G matterC ( k ) (cid:1) , (A14) G BoostQ ( k ) = G S E ( k ) G angQ ( k ) + G S E ( k ) (cid:16) G (0)BoostedQ ( k ) − G (0)Q ( k ) (cid:17) , (A15)where the boosted momentum k boosted is defined by Eq. (45) and the boosted version of G (0)Q is G (0)BoostedQ ( k ) ≡ √ m d k (cid:90) ∞ w ( r ) (cid:18) u ( r ) − w ( r )2 √ (cid:19) j (cid:18) k boosted r (cid:19) dr. (A16)Boost corrections do not contribute to the deuteron structure radius and quadrupole moment.
5. One-pion-exchange contributions
One-pion-exchange (OPE) contributions to the deuteron form factors originate from the charge density operator givenby Eq. (36). In momentum space, the expressions for the OPE contributions involve six-dimensional integration andare somewhat cumbersome. The Fourier transform to coordinate space makes these expressions much shorter andthe number of integrations reduces to one. Below we give the OPE contributions in coordinate space. For the sake ofcompactness, we introduce the functions ¯ h ( x ) and ¯ h ( x ) that correspond to the Fourier transforms of the regularizedsingle and squared pion propagators, respectively:¯ h ( r ) ≡ (cid:90) d l (2 π ) F ( l , Λ) e i l · r l + M π , ¯ h ( r ) ≡ (cid:90) d l (2 π ) F ( l , Λ) e i l · r ( l + M π ) , (A17)where F ( l , Λ) and F ( l , Λ) are the corresponding momentum-space regulators. Without regularization (i.e. when F ( l , Λ) = F ( l , Λ) = 1), the functions ¯ h ( r ) and ¯ h ( r ) take a simple form:¯ h unreg1 ( r ) = e − M π r πr , ¯ h unreg2 ( r ) = e − M π r πM π . (A18)For the regulator employed in the SMS NN potentials of Ref. [7] with F SMS1 ( l , Λ) ≡ exp (cid:18) l + M π Λ (cid:19) , F SMS2 ( l , Λ) ≡ exp (cid:18) l + M π Λ (cid:19) (cid:20) l + M π Λ (cid:21) (A19)we get the following closed form of the function ¯ h ( r )¯ h SMS1 ( r ) = exp( − M π r ) erfc (cid:0) M π Λ − Λ r (cid:1) − exp( M π r ) erfc (cid:0) M π Λ + Λ r (cid:1) πr . (A20)The function ¯ h enters the final result only under a derivative operator. To simplify the expressions even further werewrite ¯ h (cid:48) ( r ) in terms of ¯ h ( r ). We employ the relation l ( l + M π ) F ( l , Λ) = − ∇ l (cid:18) l + M π F ( l , Λ) (cid:19) , (A21)4which is fulfilled by both the unregularized and SMS-regularized pion propagators. Substituting the relation inEq. (A21) in the definition of ¯ h , taking the derivative and integrating by parts leads to the following relation incoordinate space: ¯ h (cid:48) ( r ) = (cid:16) − r (cid:17) ¯ h ( r ) . (A22)Using the simplifications above, the OPE contribution to the deuteron charge form factor can be written as: G π C ( k ) = G S E ( k ) g A F π m N ∞ (cid:90) dr (cid:32) (2 ¯ β − kj (cid:18) kr (cid:19) (cid:16) C ( r ) (cid:0) r ¯ h (cid:48)(cid:48) ( r ) + 4¯ h (cid:48) ( r ) (cid:1) + 4 √ Q ( r ) (cid:0) r ¯ h (cid:48)(cid:48) ( r ) + ¯ h (cid:48) ( r ) (cid:1)(cid:17) +(1 − β ) kj (cid:18) kr (cid:19) (cid:16) C ( r ) + 4 √ Q ( r ) (cid:17) ¯ h (cid:48) ( r ) (cid:33) , (A23)where j n ( x ) are the spherical Bessel functions. The OPE contribution to the deuteron quadrupole form factor reads: G π Q ( k ) = G S E ( k ) g A m d F π m N ∞ (cid:90) dr (cid:40) (2 ¯ β − × (cid:32) k r j (cid:18) kr (cid:19) (cid:16) − C ( r ) (cid:0) ¯ h (cid:48) ( r ) − r ¯ h (cid:48)(cid:48) ( r ) (cid:1) + √ Q ( r ) (cid:0) h (cid:48) ( r ) − r ¯ h (cid:48)(cid:48) ( r ) (cid:1) + 9 w ( r ) ¯ h (cid:48) ( r ) (cid:17) − k j (cid:18) kr (cid:19) (cid:16) C ( r ) (cid:0) r ¯ h (cid:48)(cid:48) ( r ) + ¯ h (cid:48) ( r ) (cid:1) + √ Q ( r ) (cid:0) h (cid:48) ( r ) − r ¯ h (cid:48)(cid:48) ( r ) (cid:1)(cid:17) (cid:33) +(1 − β ) (cid:18) k r j (cid:18) kr (cid:19) w ( r ) ¯ h (cid:48) ( r ) − k j (cid:18) kr (cid:19) (cid:18) C ( r ) − Q ( r ) √ (cid:19) ¯ h (cid:48) ( r ) (cid:19) (cid:41) . (A24)Finally, the OPE contributions to the deuteron structure radius and quadrupole moment have the form: r π = − g A F π m N ∞ (cid:90) dr r (cid:32) (2 ¯ β − (cid:16) C ( r ) (cid:0) r ¯ h (cid:48)(cid:48) ( r ) + 4¯ h (cid:48) ( r ) (cid:1) + 4 √ Q ( r ) (cid:0) r ¯ h (cid:48)(cid:48) ( r ) + ¯ h (cid:48) ( r ) (cid:1)(cid:17) +2(1 − β ) (cid:16) C ( r ) + 4 √ Q ( r ) (cid:17) ¯ h (cid:48) ( r ) (cid:33) , (A25) Q π = g A F π m N ∞ (cid:90) dr r (cid:32) (2 ¯ β − (cid:16) − C ( r ) (cid:0) r ¯ h (cid:48)(cid:48) ( r ) + 4¯ h (cid:48) ( r ) (cid:1) + 2 √ Q ( r ) (cid:0) r ¯ h (cid:48)(cid:48) ( r ) + ¯ h (cid:48) ( r ) (cid:1) + 27 w ( r ) ¯ h (cid:48) ( r ) (cid:17) − (1 − β )¯ h (cid:48) ( r ) (cid:16) C ( r ) − √ Q ( r ) − w ( r ) (cid:17) (cid:33) . (A26)Our analytic expressions for OPE contributions agree with the ones of Ref. [97] after the following notational changesare performed: ¯ h → M π π h, g A M π πF π → f , ¯ β → µ −
14 ¯ β → ν . (A27)5
6. Contact charge density contributions
Contact N LO contributions to the deuteron form factors stem from the corresponding short-range charge densityoperators in Eq. (39). The contact contribution to the deuteron charge form factor is given by G ContC ( k ) = 1 π G S E ( k ) ∞ (cid:90) p dp ∞ (cid:90) p (cid:48) dp (cid:48) F Λ (cid:18) p − k , p (cid:48) (cid:19) × (cid:2) F uuG C ( p, p (cid:48) , k ) u ( p ) u ( p (cid:48) ) + F uwG C ( p, p (cid:48) , k ) w ( p ) u ( p (cid:48) ) (cid:3) + ( k → − k ) , (A28)where F Λ ( p, p (cid:48) ) = exp (cid:18) − p + p (cid:48) Λ (cid:19) , (A29) F uuG C ( p, p (cid:48) , k ) = (cid:18) A + B + C (cid:19) kp (cid:32) Λ + Λ (cid:32)(cid:18) p − k (cid:19) − p (cid:48) (cid:33)(cid:33) , (A30) F uwG C ( p, p (cid:48) , k ) = √ C (cid:18) Λ kp + Λ (4 p − k )3 kp + Λ ( k − p )3 p (cid:19) , (A31)and ( k → − k ) means that the same contribution, but with opposite sign of k should be added. The contact contri-bution to the deuteron quadrupole form factor reads G ContQ ( k ) = m d π G S E ( k ) ∞ (cid:90) p dp ∞ (cid:90) p (cid:48) dp (cid:48) F Λ (cid:18) p − k , p (cid:48) (cid:19) (cid:2) F uuG Q ( p, p (cid:48) , k ) u ( p ) u ( p (cid:48) )+ F uwG Q ( p, p (cid:48) , k ) w ( p ) u ( p (cid:48) ) + F wwG Q ( p, p (cid:48) , k ) w ( p ) w ( p (cid:48) ) (cid:3) + ( k → − k ) , (A32)where F uuG Q ( p, p (cid:48) , k ) = ( − C Λ k p (cid:32) k (cid:18) p − k (cid:19) + k ( k − p )Λ + 3Λ (cid:33) ,F uwG Q ( p, p (cid:48) , k ) = ( A + B ) ( − √ k p (cid:32) k p (cid:0) ( k − p ) − p (cid:48) (cid:1) Λ − kp (cid:0) k − kp + 12( p − p (cid:48) ) (cid:1) Λ +3 (cid:0) k − kp + 4( p − p (cid:48) ) (cid:1) Λ + 36Λ (cid:33) + C √ k p (cid:16) k p (cid:16) ( k − p ) + 4 p (cid:48) (cid:17) Λ − kp (cid:0) k − kp + 12 (cid:0) p + p (cid:48) (cid:1)(cid:1) Λ + 3 (cid:0) k + 4 (cid:0) p + p (cid:48) (cid:1)(cid:1) Λ (cid:17) ,F wwG Q ( p, p (cid:48) , k ) = C p (cid:48) (cid:0) k p Λ − kp Λ + 3Λ (cid:1) k p . (A33)Next, the contact charge density contribution to the deuteron structure radius has the form r = 1 π ∞ (cid:90) p dp ∞ (cid:90) p (cid:48) dp (cid:48) F Λ ( p, p (cid:48) ) (cid:2) F uur ( p, p (cid:48) ) u ( p ) u ( p (cid:48) ) + F uwr ( p, p (cid:48) ) w ( p ) u ( p (cid:48) ) (cid:3) , (A34)6where F uur ( p, p (cid:48) ) ≡ − (cid:18) A + B + C (cid:19) (cid:32) − p + p (cid:48) )Λ + ( p − p (cid:48) ) Λ (cid:33) , (A35) F uwr ( p, p (cid:48) ) ≡ √ C (cid:32) p Λ + p (cid:0) p (cid:48) − p (cid:1) Λ (cid:33) . (A36)Finally, the contact contribution to the quadrupole moment reads: Q Cont = 1 π ∞ (cid:90) p dp ∞ (cid:90) p (cid:48) dp (cid:48) F Λ ( p, p (cid:48) ) (cid:2) F uuQ ( p, p (cid:48) ) u ( p ) u ( p (cid:48) ) + F uwQ ( p, p (cid:48) ) w ( p ) u ( p (cid:48) ) + F wwQ ( p, p (cid:48) ) w ( p ) w ( p (cid:48) ) (cid:3) , (A37)where F uuQ ( p, p (cid:48) ) ≡ ( − C ) (cid:18) − p + p (cid:48) )3Λ + 2( p + p (cid:48) )15Λ (cid:19) , F wwQ ( p, p (cid:48) ) ≡ C p p (cid:48) Λ , (A38) F uwQ ( p, p (cid:48) ) ≡ √ p (cid:18) ( A + B ) (cid:18) + 3( p (cid:48) − p )Λ (cid:19) + C (cid:18) − + p + p (cid:48) Λ (cid:19)(cid:19) . (A39) Appendix B: Complete expressions for the contact charge density at N LO including isovector terms
In this appendix we present the N LO contact charge density operators including isovector contributions. The isovectorcomponents do not contribute to the deuteron observables in the single-photon approximation, but have to be takeninto account when calculating the FFs and charge radii of heavier nuclei. Charge-density operators presented here arederived using the same procedure as used for derivation of Eq. (33), but keeping the isovector terms. After calculatingand antisymmetrizing the commutators of the LO charge density with the generators of the unitary transformationEq. (31) we obtain the following result for the N LO contact charge density: ρ ( A + B + C/ = 2 e (cid:18) A + B + C (cid:19) σ · σ + 34 (cid:20) G S E ( k ) 1 − τ · τ k + G V E ( k ) (cid:18) ( τ − τ ) k · ( p − p (cid:48) ) − i ( τ × τ ) k · ( p + p (cid:48) ) (cid:19) (cid:21) , (B1) ρ ( A − B − C )Cont,AS = 2 e ( A − B − C ) 1 − σ · σ (cid:20) (cid:18) G S E ( k ) τ · τ + 34 + G V E ( k ) ( τ + τ ) (cid:19) k + G V E ( k ) (cid:18) ( τ − τ ) k · ( p − p (cid:48) ) + i ( τ × τ ) k · ( p + p (cid:48) ) (cid:19) (cid:21) , (B2) ρ ( C )Cont,AS = 2 e C (cid:20) G S E ( k ) 1 − τ · τ (cid:18) ( k · σ )( k · σ ) − k ( σ · σ ) (cid:19) + G V E ( k ) ( τ − τ ) (cid:18) ( k · σ ) σ · ( p − p (cid:48) ) + ( k · σ ) σ · ( p − p (cid:48) ) − k · ( p − p (cid:48) )( σ · σ ) (cid:19) (B3) − G V E ( k ) i ( τ × τ ) (cid:18) ( k · σ ) σ · ( p + p (cid:48) ) + ( k · σ ) σ · ( p + p (cid:48) ) − k · ( p + p (cid:48) )( σ · σ ) (cid:19) (cid:21) . Notice that all isoscalar operators are proportional to G S E ( k ), while all isovector ones are proportional to G V E ( k ).Finally we would like to make a remark about the S → S contact operator in the first line of Eq. (B2), whichinvolves the isospin operator ( τ + τ ) . This structure is remarkable in several ways. First, from all presentedisovector terms, this is the only one which is allowed by the Pauli principle in S-to-S-wave transitions. Second, this7structure ensures that correct nucleon form factors appear in all isospin-1-to-isospin-1 channels, namely: G S E ( k ) τ · τ + 34 + G V E ( k ) ( τ + τ ) G p E for pp → ppG p E + G n E for pn → pn G n E for nn → nn (B4)Our derivation of the contact charge density operator demonstrates that the isovector structure in Eq. (B2) should beproportional to the same linear combinations of LECs as corresponding isoscalar part. This is in contrast to Ref. [101],where an extra LEC associated with the isovector terms was introduced. [1] E. Epelbaum, H. W. Hammer and U.-G. Meißner, Rev. Mod. Phys. , 1773 (2009).[2] E. Epelbaum and U.-G. Meißner, Ann. Rev. Nucl. Part. Sci. , 159 (2012).[3] E. Epelbaum, H. Krebs and P. Reinert, Front. in Phys. , 98 (2020).[4] R. Machleidt and D. R. Entem, Phys. Rept. , 1 (2011).[5] D. R. Entem, N. Kaiser, R. Machleidt and Y. Nosyk, Phys. Rev. C , no. 1, 014002 (2015).[6] D. R. Entem, N. Kaiser, R. Machleidt and Y. Nosyk, Phys. Rev. C , no. 6, 064001 (2015).[7] P. Reinert, H. Krebs and E. Epelbaum, Eur. Phys. J. A , no. 5, 86 (2018).[8] R. Machleidt, Phys. Rev. C , 024001 (2001).[9] V. G. J. Stoks, R. A. M. Klomp, C. P. F. Terheggen and J. J. de Swart, Phys. Rev. C , 2950 (1994).[10] R. B. Wiringa, V. G. J. Stoks and R. Schiavilla, Phys. Rev. C , 38 (1995).[11] E. Epelbaum, H. Krebs and U.-G. Meißner, Eur. Phys. J. A , no. 5, 53 (2015).[12] E. Epelbaum, H. Krebs and U.-G. Meißner, Phys. Rev. Lett. , no. 12, 122301 (2015).[13] A. Gezerlis, I. Tews, E. Epelbaum, S. Gandolfi, K. Hebeler, A. Nogga and A. Schwenk, Phys. Rev. Lett. , no. 3,032501 (2013).[14] M. Piarulli, L. Girlanda, R. Schiavilla, R. Navarro P´erez, J. E. Amaro and E. Ruiz Arriola, Phys. Rev. C , no. 2,024003 (2015).[15] D. R. Entem, R. Machleidt and Y. Nosyk, Phys. Rev. C , no. 2, 024004 (2017).[16] M. Hoferichter, J. Ruiz de Elvira, B. Kubis and U.-G. Meißner, Phys. Rev. Lett. , no. 19, 192301 (2015).[17] M. Hoferichter, J. Ruiz de Elvira, B. Kubis and U.-G. Meißner, Phys. Rept. , 1 (2016).[18] P. Reinert, H. Krebs and E. Epelbaum, [arXiv:2006.15360 [nucl-th]].[19] S. Binder et al. [LENPIC Collaboration], Phys. Rev. C , no. 4, 044002 (2016).[20] S. Binder et al. [LENPIC Collaboration], Phys. Rev. C , no. 1, 014002 (2018).[21] E. Epelbaum et al. [LENPIC Collaboration], Phys. Rev. C , no. 2, 024313 (2019).[22] R. Skibi´nski et al. , Phys. Rev. C , no. 6, 064002 (2016).[23] D. L. Yao, D. Siemens, V. Bernard, E. Epelbaum, A. M. Gasparyan, J. Gegelia, H. Krebs and U.-G. Meißner, JHEP , 038 (2016).[24] D. Siemens, V. Bernard, E. Epelbaum, A. M. Gasparyan, H. Krebs and U.-G. Meißner, Phys. Rev. C , no. 5, 055205(2017).[25] J. E. Lynn, D. Lonardoni, J. Carlson, J. W. Chen, W. Detmold, S. Gandolfi and A. Schwenk, J. Phys. G , no. 4, 045109(2020).[26] N. Nevo Dinur, O. J. Hernandez, S. Bacca, N. Barnea, C. Ji, S. Pastore, M. Piarulli and R. B. Wiringa, Phys. Rev. C , no. 3, 034004 (2019).[27] A. N. Hiller Blin, Z. F. Sun and M. J. Vicente Vacas, Phys. Rev. D , no. 5, 054025 (2018).[28] D. Lonardoni, S. Gandolfi, J. E. Lynn, C. Petrie, J. Carlson, K. E. Schmidt and A. Schwenk, Phys. Rev. C , no. 4,044318 (2018).[29] R. J. Furnstahl, N. Klco, D. R. Phillips and S. Wesolowski, Phys. Rev. C , no. 2, 024005 (2015).[30] J. A. Melendez, S. Wesolowski and R. J. Furnstahl, Phys. Rev. C , no. 2, 024003 (2017).[31] S. Wesolowski, R. J. Furnstahl, J. A. Melendez and D. R. Phillips, J. Phys. G , no. 4, 045102 (2019).[32] E. Epelbaum, J. Golak, K. Hebeler, H. Kamada, H. Krebs, U.-G. Meißner, A. Nogga, P. Reinert, R. Skibinski, K. Topol-nicki, Y. Volkotrub and H. Witala, Eur. Phys. J. A , no. 3, 92 (2020).[33] M. Garcon and J. W. Van Orden, Adv. Nucl. Phys.
293 (2001).[34] R. A. Gilman and F. Gross, J. Phys. G , R37 (2002).[35] L. E. Marcucci et al. , J. Phys. G , 023002 (2016).[36] D. R. Phillips, Nucl. Phys. A , 52 (2004).[37] J. W. Chen, G. Rupak and M. J. Savage, Nucl. Phys. A , 386 (1999).[38] D. R. Phillips and T. D. Cohen, Nucl. Phys. A , 45 (2000).[39] M. Walzl and U.-G. Meißner, Phys. Lett. B , 37 (2001).[40] D. R. Phillips, Phys. Lett. B , 12 (2003). [41] D. R. Phillips, J. Phys. G , 365 (2007).[42] M. P. Valderrama, A. Nogga, E. Ruiz Arriola and D. R. Phillips, Eur. Phys. J. A , 315 (2008).[43] M. Piarulli, L. Girlanda, L. E. Marcucci, S. Pastore, R. Schiavilla and M. Viviani, Phys. Rev. C , no. 1, 014006 (2013).[44] E. Epelbaum, A. M. Gasparyan, J. Gegelia and M. R. Schindler, Eur. Phys. J. A , 51 (2014).[45] A. A. Filin, V. Baru, E. Epelbaum, H. Krebs, D. M¨oller and P. Reinert, Phys. Rev. Lett. , no. 8, 082501 (2020).[46] A. Cipollone, C. Barbieri and P. Navr´atil, Phys. Rev. C , no. 1, 014306 (2015).[47] S. K. Bogner, R. J. Furnstahl, P. Maris, R. J. Perry, A. Schwenk and J. P. Vary, Nucl. Phys. A , 21 (2008).[48] Y. B. Dong, Phys. Rev. C , 025208 (2009).[49] R. G. Arnold, C. E. Carlson and F. Gross, Phys. Rev. C , 1426 (1980).[50] M. Tanabashi et al. [Particle Data Group], Phys. Rev. D , no. 3, 030001 (2018).[51] R. G. Arnold, C. E. Carlson and F. Gross, Phys. Rev. C , 363 (1981).[52] T. E. O. Ericson and M. Rosa-Clot, Nucl. Phys. A , 497 (1983).[53] D. M. Bishop and L. M. Cheung, Phys. Rev. A , 381 (1979).[54] P. J. Mohr, D. B. Newell and B. N. Taylor, Rev. Mod. Phys. , no. 3, 035009 (2016).[55] D. Abbott et al. [JLAB t20 Collaboration], Eur. Phys. J. A , 421 (2000).[56] D. Abbott et al. [JLAB t(20) Collaboration], Phys. Rev. Lett. , 5053 (2000).[57] D. M. Nikolenko et al. , Phys. Rev. Lett. , 072501 (2003).[58] Y. B. Dong and D. Y. Chen, Phys. Lett. B , 426 (2009).[59] A. P. Kobushkin, Y. D. Krivenko-Emetov and S. Dubnicka, Phys. Rev. C , 054001 (2010).[60] F. Gross, Phys. Rev. C , no. 2, 024001 (2020).[61] S. K¨olling, E. Epelbaum, H. Krebs and U.-G. Meißner, Phys. Rev. C , 045502 (2009).[62] S. K¨olling, E. Epelbaum, H. Krebs and U.-G. Meißner, Phys. Rev. C , 054008 (2011).[63] H. Krebs, E. Epelbaum and U.-G. Meißner, Few Body Syst. , no. 2, 31 (2019).[64] S. Pastore, R. Schiavilla and J. L. Goity, Phys. Rev. C , 064002 (2008).[65] S. Pastore, L. Girlanda, R. Schiavilla, M. Viviani and R. B. Wiringa, Phys. Rev. C , 034004 (2009).[66] S. Pastore, L. Girlanda, R. Schiavilla and M. Viviani, Phys. Rev. C , 024001 (2011).[67] T. S. Park, D. P. Min and M. Rho, Nucl. Phys. A , 515 (1996).[68] V. Bernard, E. Epelbaum, H. Krebs and U.-G. Meißner, Phys. Rev. C , 064004 (2008).[69] V. Bernard, E. Epelbaum, H. Krebs and U.-G. Meißner, Phys. Rev. C , 054001 (2011).[70] H. Krebs, A. Gasparyan and E. Epelbaum, Phys. Rev. C , 054006 (2012).[71] H. Krebs, A. Gasparyan and E. Epelbaum, Phys. Rev. C , no. 5, 054007 (2013).[72] E. Epelbaum, A. M. Gasparyan, H. Krebs and C. Schat, Eur. Phys. J. A , no. 3, 26 (2015).[73] H. Krebs, [arXiv:2008.00974 [nucl-th]].[74] J. L. Friar, J. Martorell and D. W. L. Sprung, Phys. Rev. A , 4579 (1997).[75] B. Kubis and U.-G. Meißner, Nucl. Phys. A , 698 (2001).[76] M. R. Schindler, J. Gegelia and S. Scherer, Eur. Phys. J. A , 1 (2005).[77] J. C. Bernauer et al. [A1 Collaboration], Phys. Rev. C , no. 1, 015206 (2014).[78] Z. Ye, private communication.[79] Z. Ye, J. Arrington, R. J. Hill and G. Lee, Phys. Lett. B , 8 (2018).[80] M. A. Belushkin, H.-W. Hammer and U.-G. Meißner, Phys. Rev. C , 035202 (2007).[81] V. Punjabi, C. F. Perdrisat, M. K. Jones, E. J. Brash and C. E. Carlson, Eur. Phys. J. A , 79 (2015).[82] I. T. Lorenz, H.-W. Hammer and U.-G. Meißner, Eur. Phys. J. A , 151 (2012).[83] I. T. Lorenz, U.-G. Meißner, H.-W. Hammer and Y.-B. Dong, Phys. Rev. D , no. 1, 014023 (2015).[84] S. Pacetti, R. Baldini Ferroli and E. Tomasi-Gustafsson, Phys. Rept. , 1 (2015).[85] D. Drechsel and T. Walcher, Rev. Mod. Phys. , 731 (2008).[86] C. F. Perdrisat, V. Punjabi and M. Vanderhaeghen, Prog. Part. Nucl. Phys. , 694 (2007).[87] J. Arrington, C. D. Roberts and J. M. Zanotti, J. Phys. G , S23 (2007).[88] R. Pohl et al. , Nature , 213 (2010).[89] A. Beyer et al. , Science , 79 (2017).[90] N. Bezginov, T. Valdez, M. Horbatsch, A. Marsman, A. C. Vutha and E. A. Hessels, Science , no. 6457, 1007 (2019).[91] E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor (2019), “The 2018 CODATA Recommended Values of theFundamental Physical Constants” (Web Version 8.0). Database developed by J. Baker, M. Douma, and S. Kotochigova.Available at http://physics.nist.gov/constants, National Institute of Standards and Technology, Gaithersburg, MD 20899.[92] J. C. Bernauer et al. [A1 Collaboration], Phys. Rev. Lett. , 242001 (2010).[93] G. H¨ohler, E. Pietarinen, I. Sabba Stefanescu, F. Borkowski, G. G. Simon, V. H. Walther and R. D. Wendling, Nucl.Phys. B , 505 (1976).[94] P. Mergell, U.-G. Meißner and D. Drechsel, Nucl. Phys. A , 367 (1996).[95] H. W. Hammer and U.-G. Meißner, Sci. Bull. , 257 (2020).[96] C. Alexandrou, M. Constantinou, K. Hadjiyiannakou, K. Jansen, C. Kallidonis, G. Koutsou and A. Vaquero Aviles-Casco,Phys. Rev. D , no. 3, 034503 (2017).[97] J. L. Friar, Phys. Rev. C , 796 (1980).[98] E. Epelbaum, W. Glockle and U.-G. Meißner, Nucl. Phys. A , 362 (2005).[99] J. L. Friar, Phys. Rev. C , 034002 (1999).[100] H. Hyuga and H. Ohtsubo, Nucl. Phys. A , 348 (1978). [101] D. R. Phillips, Ann. Rev. Nucl. Part. Sci. , 421 (2016).[102] H. Krebs, PoS CD2018 , 098 (2019).[103] E. Epelbaum, [arXiv:1908.09349 [nucl-th]].[104] H. Arenhovel, F. Ritz and T. Wilbois, Phys. Rev. C , 034002 (2000).[105] R. A. Krajcik and L. L. Foldy, Phys. Rev. D , 1777 (1974).[106] J. L. Friar, Annals Phys. , 380 (1977).[107] F. Ritz, H. Goller, T. Wilbois and H. Arenhovel, Phys. Rev. C , 2214 (1997).[108] S. J. Wallace, Phys. Rev. Lett. , 180401 (2001).[109] R. Schiavilla and V. R. Pandharipande, Phys. Rev. C , 064009 (2002).[110] V. G. J. Stoks, R. A. M. Klomp, M. C. M. Rentmeester and J. J. de Swart, Phys. Rev. C , 792 (1993).[111] U. D. Jentschura et al. , Phys. Rev. A , 042505 (2011).[112] K. Pachucki, V. Patk´oˇs and V. A. Yerokhin, Phys. Rev. A , no. 6, 062511 (2018).[113] R. Pohl, R. Gilman, G. A. Miller and K. Pachucki, Ann. Rev. Nucl. Part. Sci. , 175 (2013).[114] R. Pohl et al. [CREMA Collaboration], Science , no. 6300, 669 (2016).[115] R. Pohl et al. , Metrologia , no. 2, L1 (2017).[116] S. Kopecky, M. Krenn, P. Riehs, S. Steiner, J. A. Harvey, N. W. Hill and M. Pernicka, Phys. Rev. C , 2229 (1997).[117] L. V. Mitsyna, V. G. Nikolenko, S. S. Parzhitski, A. B. Popov and G. S. Samosvat, Nucl. Phys. A , 1 (2009).[118] B. D. Carlsson, A. Ekstr¨om, C. Forss´en, D. F. Str¨omberg, G. R. Jansen, O. Lilja, M. Lindby, B. A. Mattsson andK. A. Wendt, Phys. Rev. X , no.1, 011019 (2016).[119] I. Sick, private communication.[120] J. A. Melendez, R. J. Furnstahl, D. R. Phillips, M. T. Pratola and S. Wesolowski, Phys. Rev. C , no. 4, 044001 (2019).[121] E. Epelbaum, PoS CD2018 , 006 (2019).[122] J. Adam, H. Goller and H. Arenhovel, Phys. Rev. C48