Precision calculation of the electromagnetic radii of the proton and neutron from lattice QCD

We present lattice-QCD results for the electromagnetic form factors of the proton and neutron including both quark-connected and -disconnected contributions. The parametrization of the $Q^2$-dependence of the form factors is combined with the extrapolation to the physical point. In this way, we determine the electric and magnetic radii and the magnetic moments of the proton and neutron. For the proton, we obtain at the physical pion mass and in the continuum and infinite-volume limit $\sqrt{\langle r_E^2 \rangle^p} = 0.820(14)$ fm, $\sqrt{\langle r_M^2 \rangle^p} = 0.8111(89)$ fm, and $\mu_M^p = 2.739(66)$, where the errors include all systematics.

Introduction.The so-called "proton radius puzzle", i.e., the observation of a large tension in the proton's electric (charge) radius extracted either from atomic spectroscopy data of muonic hydrogen [1,2] or, alternatively, from corresponding measurements on electronic hydrogen [3] as well as ep-scattering data [4,5], has gripped the scientific community for more than 10 years and triggered a vigorous research effort designed to explain the discrepancy.
While the situation regarding the electric radius is awaiting its final resolution, one also finds discrepant results for the proton's magnetic radius.Specifically, there is a tension of 2.7 σ between the value extracted from the A1 ep-scattering data alone and the estimate from the corresponding analysis applied to the remaining world data [24].Clearly, a firm theoretical prediction for basic properties of the proton and the neutron, such as their radii and magnetic moments, would be highly desirable in order to assess our understanding of the particles that make up the largest fraction of the visible mass in the universe.
In this letter we present our results for the radii and magnetic moment of the proton computed in lattice QCD.Compared with previous lattice studies [25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41][42], our calculation is the first to include the contributions from quark-disconnected diagrams while controlling all sources of systematic uncertainties arising from excited-state con-tributions, finite-volume effects and the continuum extrapolation.We determine the proton's magnetic radius ⟨r 2 M ⟩ p with a total precision of 1.1 %, which is competitive with recent analyses of ep-scattering data [4,15,16,24].Moreover, our lattice QCD estimate for the proton's magnetic moment is in good agreement with experiment.Our result for the electric radius, which has a total precision of 1.7 %, is consistent with the value determined in muonic hydrogen within 1.5 standard deviations.
Lattice setup.Our aim is to compute the electric and magnetic Sachs form factors G E (Q 2 ) and G M (Q 2 ) of the proton and neutron.The electric form factor at zero momentum transfer yields the nucleon's electric charge, i.e., G p E (0) = 1 and G n E (0) = 0, whereas the magnetic form factor at Q 2 = 0 is identified with the magnetic moment, G M (0) = µ M .The radii can in turn be extracted from the slope of the form factors at zero momentum transfer, .
The only exception to this definition is the electric radius of the neutron, where the normalization factor 1/G(0) is dropped.
For our lattice determination of these quantities, we use the ensembles generated by the Coordinated Lattice Simulations (CLS) [43] effort with 2 + 1 flavors of nonperturbatively O(a)-improved Wilson fermions [44,45] and a tree-level improved Lüscher-Weisz gauge action [46], correcting for the treatment of the strange quark determinant using the procedure outlined in Ref. [47].Table I shows the set of ensembles entering the analysis: they cover four lattice spacings in the range from 0.050 fm to 0.086 fm, and several pion masses, including one slightly below the physical value (E250).Further details on our setup of the simulations and the measurements of the two-and three-point functions of the nucleon can be found in the accompanying paper [48].To extract the effective form factors from two-and threepoint correlation functions, we employ the ratio method [49] and the same estimators for the effective electric and magnetic Sachs form factors as in Ref. [41].For further technical details, we again refer to the companion paper [48].The effective form factors are constructed for the isovector (u−d) and the connected isoscalar (u+d) combinations, as well as for the light and strange disconnected contributions.In the isovector case, the disconnected contributions cancel.The full isoscalar (octet) combination u + d − 2s, on the other hand, is obtained from the connected and disconnected pieces as Note that the disconnected part only requires the difference l − s between the light and strange contributions, in which correlated noise cancels and which can be computed efficiently by the one-end trick [50][51][52].
We express all dimensionful quantities in units of the gradient flow time t 0 [53] using the determination of t sym 0 /a 2 from Ref. [54].Only in the final step, i.e., after the extrapolation to the physical point, are the radii converted to physical units by means of the FLAG estimate of t 0,phys = 0.14464(87) fm for N f = 2 + 1 from Ref. [55].
Excited-state analysis.Due to the strong exponential decay of the signal-to-noise ratio for baryonic correlation functions with increasing source-sink separation [56,57], an explicit treatment of the excited-state systematics is required in order to extract the ground-state form factors from the effective ones [58].In this work, we employ the summation method [59][60][61].It exploits the fact that the contributions from excited states to effective form factors are parametrically more strongly suppressed when the insertion of the electromagnetic current is summed over timeslices between the source and sink.In the asymptotic limit, the slope of the summed correlator ratio with respect to the source-sink separation t sep yields the ground-state form factor [34,41].
In our analysis, we monitor the stability of fit results for different starting values t min sep of the source-sink separation.Rather than selecting one particular value of t min sep on each ensemble, we perform a weighted average over t min sep , where the weights are given by a smooth window function [62,63] [cf.eq. ( 18) in the accompanying paper [48]].
This averaging strategy is illustrated in fig. 1 for the isoscalar combination at the first non-vanishing momentum on ensemble E300.One finds that the window averages agree within their errors with what one can identify as plateaux in the blue points.This is observed for practically all other ensembles and momenta employed in the analysis, and hence we conclude that the window method reliably isolates the asymptotic value.Moreover, it reduces the human bias arising from manually picking one particular value for t min sep on each ensemble, because we use the same window parameters in units of t 0 on all ensembles.Since our window average does not yield a significantly smaller error in comparison with the individual points entering the average (cf.fig.1), we are confident that our error estimates are sufficiently conservative to exclude any systematic bias in estimating the ground-state form factors.
Further crosschecks on our excited-state analysis, including a detailed comparison with an alternative approach based on two-state fits to the effective form factors, can be found in appendix B of the accompanying paper [48].We did not find any evidence that our preferred strategy presented above introduces a systematic bias or underestimates the errors, but is on the contrary rather conservative in this regard.
Direct BχPT fits.To extract the radii from the form factors, we need to describe the Q 2 -dependence of the latter.As in Refs.[34,41], we employ two different methods: our preferred procedure is to combine the parametrization of the Q 2 -dependence with the extrapolation to the physical point (M π = M π,phys , a = 0, L = ∞) by fitting our form factor data directly to the expressions resulting from covariant baryon chiral perturbation theory (BχPT) [64].This is presented in the following.Alternatively, we have implemented the more traditional strategy of first performing a generic parametrization of the Q 2 -dependence on each ensemble, followed by extrapolating the resulting radii to the physical point.A crosscheck of our main analysis with this two-step approach can be found in the accompanying paper [48].
For our main analysis, we fit our form factor data to the full expressions of Ref. [64] without explicit ∆ degrees of freedom.The fits are performed for the isovector  and isoscalar channels separately, but for G E and G M simultaneously.This allows us to properly treat the correlations between different Q 2 and also between G E and G M .Different gauge ensembles, on the other hand, are treated as statistically independent.In the isovector channel, we include the contributions from the ρ meson in the expressions for the form factors, while in the isoscalar channel, we include the leading-order terms from the ω and ϕ resonances.The physical pion mass M π,phys is fixed in units of √ t 0 using its value in the isospin limit [65], i.e., we employ √ t 0,phys M π,phys = 0.09881 (59).Here, we neglect the uncertainty of M π,iso in MeV since it is completely subdominant compared to that of √ t 0,phys which enters in the conversion of units.
We perform several such fits with various cuts in the pion mass (M π ≤ 0.23 GeV and M π ≤ 0.27 GeV) and the momentum transfer (Q 2 ≤ 0.3, . . ., 0.6 GeV 2 ), as well as with different models for the lattice-spacing and/or finite-volume dependence, in order to estimate the corresponding systematic uncertainties.The variations of the results due to the cuts are in most cases much smaller than their statistical errors and will be included in our quoted systematic errors.We reconstruct the proton and neutron form factors as linear combinations of the BχPT formulae for the isovector and isoscalar channels, evaluating the low-energy constants as determined from the separate fits in these channels.For further technical details of the fits, we refer to the accompanying paper [48].
The major benefits of the direct approach compared to the two-step procedure have been previously observed in our publication on the isovector electromagnetic form factors [41]: Firstly, including multiple ensembles as well as G E and G M in one fit significantly reduces the resulting errors on the radii.Secondly, it greatly increases the number of degrees of freedom in the fit, which has a stabilizing effect with regard to lowering the momentum cut.
Model average and final results.Since we do not have a strong a priori preference for one specific setup of the direct fits, we determine our final results and total errors from averages over different fit models and kinematic cuts.For this purpose, we use weights derived from the Akaike Information Criterion [66][67][68][69][70].In order to estimate the statistical and systematic uncertainties of our model averages, we adopt a bootstrapped variant of the method from Ref. [71].Our procedure is explained in more detail in the accompanying paper [48].As our final results, we obtain ⟨r 2 E ⟩ p = (0.672 ± 0.014 (stat) ± 0.018 (syst)) fm 2 , (5) ⟨r 2 M ⟩ p = (0.658 ± 0.012 (stat) ± 0.008 (syst)) fm 2 , (6) We note that the precision of the magnetic radius of the proton, ⟨r 2 M ⟩ p = (0.8111 ± 0.0074 (stat) ± 0.0050 (syst)) fm, is commensurate with that of its electric counterpart, ⟨r 2 E ⟩ p = (0.820 ± 0.009 (stat) ± 0.011 (syst)) fm.
To further compare our results to experiment we perform model averages of the form factors themselves.The results are plotted in fig. 2 for the proton.One observes that the slope of the electric form factor as obtained from our calculation is closer to the PRad measurement [6] than to that of the A1 Collaboration [4].The magnetic form factor, on the other hand, agrees well with the A1 data.Moreover, our estimates reproduce within their errors the experimental results for the magnetic moments  [4] obtained using Rosenbluth separation, and the green diamonds the corresponding data by PRad [6].The experimental value of the magnetic moment [72] is depicted by a red cross.
both of the proton and of the neutron [72].The plots for the neutron corresponding to fig. 2 in this letter are contained in fig.7 of the accompanying paper [48].
In fig.3, our results for the electromagnetic radii and magnetic moment of the proton are compared to recent lattice determinations and to the experimental values.We note that the only other lattice result including disconnected contributions is ETMC19 [39], which, however, has not been extrapolated to the continuum and infinitevolume limits.Our estimate for the electric radius is larger than the results of Refs.[38][39][40], while Ref. [32] quotes an even larger central value.
We stress that any difference between our estimate and previous lattice calculations is not related to our preference for direct fits to the form factors over the conventional approach via the z-expansion, as the latter yields consistent values for the radii (cf. the accompanying paper [48]).For the magnetic radius, our result agrees with that of Refs.[38,39] within 1.2 combined standard deviations, while that of Ref. [31] is much smaller.Our statistical and systematic error estimates for the electric radius and magnetic moment are similar or smaller compared to other lattice studies, while being substantially smaller for the magnetic radius.As a final remark we note that the lack of a data point at Q 2 = 0 complicates the extraction of the magnetic low-Q 2 observables in most recent lattice determinations, especially for z-expansion fits on individual ensembles.By contrast, the direct approach -in addition to combining information from several ensembles and from G E and G M -is more constraining at low Q 2 , allowing for considerably less variation in the form factors in that regime.We believe this to be responsible, to a large extent, for the small errors we achieve in the magnetic radii.
Conclusions.We have performed the first lattice QCD calculation of the radii and magnetic moment of the proton to include the contributions from quark-connected 0 .7 5 0 .8 0 0 .8 5 .Comparison of our best estimates for the electromagnetic radii and the magnetic moment of the proton with other lattice calculations, i.e., Mainz21 [41], ETMC20 [40], ETMC19 [39], PACS19 [38], and CSSM/QCDSF/UKQCD14 [31,32].Only ETMC19 and this work include disconnected contributions.The Mainz21 values have been obtained by combining their isovector results with the PDG values for the neutron [72].We also show this estimate using our updated isovector results (cf. the accompanying paper [48]).The experimental value for µ p M is taken from PDG [72].The two data points for ⟨r 2 E ⟩ p depict the values from PDG [72] (cross) and Mainz/A1 [4] (square), respectively.The two data points for ⟨r 2 M ⟩ p , on the other hand, show the reanalysis of Ref. [24] either using the world data excluding that of Ref. [4] (diamond) or using only that of Ref. [4] (square).
and -disconnected diagrams and present a full error budget.The overall precision of our calculation is sufficient to make a meaningful contribution to the debate surrounding the proton radii.Our final estimates are listed in eqs.( 5) to (10).
As an important benchmark, we reproduce the experi-mentally very precisely known magnetic moments of the proton and neutron [72] within our quoted uncertainties.A detailed discussion of our results for the neutron radii can be found in the accompanying paper [48].Our result for the electric (charge) radius of the proton is much closer to the value inferred from muonic hydrogen spectroscopy [2] and the recent ep-scattering experiment by PRad [6] than to the A1 ep-scattering result [4].For the magnetic radius, on the other hand, our estimate is well compatible with the analyses [4,24] of the A1 data and exhibits a 2.8 σ tension with the other collected world data [24].The analyses of combined A1+PRad data [16] and A1 data alone [15], based respectively on dispersive and dispersively improved fit ansätze, arrive at values of ⟨r 2 M ⟩ p that are significantly larger than the A1-data analysis [24] and in tension with our result.This could partly be due to unaccounted-for isospin-breaking effects.The shape of the magnetic form factor determined in our lattice calculation, however, agrees very well with the measurement by A1 [4] over the whole range of Q 2 under study.
Returning to the (electric) proton radius puzzle, our lattice results lend further support to the emerging consensus that the issue has essentially been settled [73][74][75].Meanwhile, the situation regarding the magnetic radius remains to be clarified.
Acknowledgments.This research is partly supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through project HI 2048/1-2 (project No. 399400745) and through the Cluster of Excellence "Precision Physics, Fundamental Interactions and Structure of Matter" (PRISMA + EXC 2118/1) funded within the German Excellence Strategy (project ID 39083149).Calculations for this project were partly performed on the HPC clusters "Clover" and "HIM-ster2" at the Helmholtz Institute Mainz.Other parts were conducted using the supercomputer "Mogon 2" offered by Johannes Gutenberg University Mainz (https: //hpc.uni-mainz.de),which is a member of the AHRP (Alliance for High Performance Computing in Rhineland Palatinate, https://www.ahrp.info)and the Gauss Alliance e.V.The authors also gratefully acknowledge the John von Neumann Institute for Computing (NIC) and the Gauss Centre for Supercomputing e.V. (https: //www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer JUWELS at Jülich Supercomputing Centre (JSC) through projects CHMZ21, CHMZ36, NUCSTRUCLFL, and GC-SNUCL2PT.
Our programs use the QDP++ library [76] and deflated SAP+GCR solver from the openQCD package [77], while the contractions have been explicitly checked using the Quark Contraction Tool [78].We thank Simon Kuberski for providing the improved reweighting factors [79] for the gauge ensembles used in our calculation.Moreover, we are grateful to our colleagues in the CLS initiative for sharing the gauge field configurations on which this work is based.

Figure 1 .
Figure 1.Isoscalar electromagnetic form factors at the lowest non-vanishing momentum (Q 2 ≈ 0.067 GeV 2 ) on ensemble E300 as a function of the minimal source-sink separation entering the summation fit.Each blue point corresponds to a single fit starting at the value given on the horizontal axis.The associated weights in the average are represented by the red diamonds, with the gray curves and bands depicting the averaged results.

Figure 2 .
Figure 2. Electromagnetic form factors of the proton as a function of Q 2 .The orange curves and bands correspond to our final results at the physical point with their full uncertainties obtained as model averages over the different direct fits.The light orange bands indicate the statistical errors.The black diamonds represent the experimental ep-scattering data by the A1 Collaboration[4] obtained using Rosenbluth separation, and the green diamonds the corresponding data by PRad[6].The experimental value of the magnetic moment[72] is depicted by a red cross.

Figure 3
Figure3.Comparison of our best estimates for the electromagnetic radii and the magnetic moment of the proton with other lattice calculations, i.e., Mainz21[41], ETMC20[40], ETMC19[39], PACS19[38], and CSSM/QCDSF/UKQCD14[31,32]. Only ETMC19 and this work include disconnected contributions.The Mainz21 values have been obtained by combining their isovector results with the PDG values for the neutron[72].We also show this estimate using our updated isovector results (cf. the accompanying paper[48]).The experimental value for µ p M is taken from PDG[72].The two data points for ⟨r 2 E ⟩ p depict the values from PDG[72] (cross) and Mainz/A1[4] (square), respectively.The two data points for ⟨r 2 M ⟩ p , on the other hand, show the reanalysis of Ref.[24] either using the world data excluding that of Ref.[4] (diamond) or using only that of Ref.[4] (square).

Table I .
[48]view of the ensembles used in this study.Further details are contained in table I of the accompanying paper[48].