# Isovector electromagnetic form factors of the nucleon from lattice QCD and the proton radius puzzle

D. Djukanovic, T. Harris, G. von Hippel, P.M. Junnarkar, H.B. Meyer, D. Mohler, K. Ottnad, T. Schulz, J. Wilhelm, H. Wittig

IIsovector electromagnetic form factors of the nucleon from latticeQCD and the proton radius puzzle

D. Djukanovic,

1, 2

T. Harris, G. von Hippel, P.M. Junnarkar, H. B. Meyer,

1, 2, 4

D. Mohler,

1, 2

K. Ottnad, T. Schulz, J. Wilhelm, and H. Wittig

1, 2, 41

Helmholtz Institute Mainz, Staudingerweg 18, D-55128 Mainz, Germany GSI Helmholtzzentrum f¨ur Schwerionenforschung, Darmstadt (Germany) School of Physics and Astronomy,University of Edinburgh, Edinburgh EH9 3JZ, UK PRISMA + Cluster of Excellence and Institute for Nuclear Physics,Johannes Gutenberg University of Mainz,Johann-Joachim-Becher-Weg 45, D-55128 Mainz, Germany Institut f¨ur Kernphysik, Technische Universit¨at Darmstadt,Schlossgartenstraße 2, 64289 Darmstadt (Dated: February 16, 2021)

Abstract

We present results for the isovector electromagnetic form factors of the nucleon computed on theCLS ensembles with N f = 2 + 1 ﬂavors of O ( a )-improved Wilson fermions and an O ( a )-improvedvector current. The analysis includes ensembles with four lattice spacings and pion masses rangingfrom 130 MeV up to 350 MeV and mainly targets the low- Q region. In order to remove anybias from unsuppressed excited-state contributions, we investigate several source-sink separationsbetween 1.0 fm and 1.5 fm and apply the summation method as well as explicit two-state ﬁts. Thechiral interpolation is performed by applying covariant chiral perturbation theory including vec-tor mesons directly to our form factor data, thus avoiding an auxiliary parametrization of the Q dependence. At the physical point, we obtain µ = 4 . stat (13) sys for the nucleon isovector mag-netic moment, in good agreement with the experimental value and (cid:104) r (cid:105) = 0 . stat (11) sys fm for the corresponding square-radius, again in good agreement with the value inferred from the ep -scattering determination [Bernauer et al., Phys. Rev. Lett., , 242001 (2010)] of the protonradius. Our estimate for the isovector electric charge radius, (cid:104) r (cid:105) = 0 . stat (22) sys fm , how-ever, is in slight tension with the larger value inferred from the aforementioned ep -scattering data,while being in agreement with the value derived from the 2018 CODATA average for the protoncharge radius. PACS numbers: 11.15.Ha, 12.38.Gc, 12.38.-t, 13.40 Gp, 14.20 DhKeywords: Lattice QCD, Electromagnetic Form Factors a r X i v : . [ h e p - l a t ] F e b . INTRODUCTION The internal structure of the nucleon still poses many open questions. Not only is thecomposition of its spin and momentum not completely understood [1–4], but even its sizeis subject to signiﬁcant uncertainty arising from discrepancies between diﬀerent determi-nations: there is a decade-old inconsistency [5, 6] between the electric charge radius of theproton as obtained from ep -scattering ( (cid:104) r (cid:105) / = 0 . (cid:104) r (cid:105) / = 0 . (cid:104) r (cid:105) / = 0 . (cid:104) r (cid:105) / using (electronic) hydrogen spectroscopy [11–13] mostly tendto give somewhat smaller values, while the newest determinations from ep -scattering givediﬀerent results, albeit with still rather large uncertainties: the A1 collaboration at MAMIuses an Initial-State Radiation setup [14] to achieve very small momentum transfers andﬁnds a large value (cid:104) r (cid:105) / = 0 . (cid:104) r (cid:105) / = 0 . ep -scattering experiment, MAGIX, is being preparedat Mainz [19]. To complement the result from muonic hydrogen spectroscopy with a resultfrom µp scattering, the MUSE collaboration aims to measure the µp scattering cross-sectionto sub-percent accuracy [20]; a similar experiment is in preparation at CERN by the COM-PASS++/AMBER collaboration [21]. At present, a signiﬁcant tension remains betweenthe determinations based on muonic hydrogen and ep -scattering, respectively. A variety ofpossible theoretical explanations have been proposed, but so far there has not been decisiveevidence in favour of any particular explanation [10, 22]. It remains notable that dispersiveanalyses favour a smaller value of the proton charge radius [23, 24], and in fact have doneso since before the muonic measurements [25].In the context of scattering experiments, the proton charge radius is determined fromthe derivative of the electric form factor G E ( Q ) at Q = 0. To determine this derivativeaccurately, high-quality data at very low momentum transfer Q and/or a parameterizationthat remains trustworthy over a large range of Q are required. In order to gain a properunderstanding of the origin of the proton radius puzzle and associated questions, theoreticaldeterminations of nucleon structure observables, and in particular of the electromagneticform factors of the nucleon, from ﬁrst principles are required. Lattice QCD calculations aretherefore instrumental in predicting the nucleon charge radii from QCD. This has generatedlively activity on this topic within the lattice QCD community [26–39], although at present,the precision of lattice results is not yet suﬃcient to rule out either the electronic or themuonic result for the proton radius.At large momentum transfers Q , the electromagnetic form factors are the subject ofanother puzzle: while polarization-transfer experiments [40] ﬁnd that the ratio of electricand magnetic form factor of the proton, µ p G E , p ( Q ) /G M , p ( Q ), decreases roughly linearlyfor large Q [41–45], experiments based on the Rosenbluth separation formula [46] ﬁndthat it is roughly constant and of order one (albeit with rapidly increasing errors at large Q , where G E , p contributes little to the total cross-section) [47–49]. This discrepancy hasbeen explained theoretically as the result of two-photon exchange contributions to the cross-section measurements [50–52], but the situation is not yet completely clariﬁed [53].In this paper, we present our lattice QCD-based determination of the isovector electro-2agnetic form factors of the nucleon from the N f = 2 + 1 CLS ensembles [54]. We usetwo state-of-the-art methods, known as the summation method and two-state ﬁts to ex-tract G E and G M from Euclidean correlation functions for a range of momentum transfers Q (cid:46) . Given the limited Q range of our data, we focus on extracting the electricand magnetic charge radii and the magnetic moment of the nucleon using a variety of meth-ods to check for consistency. Extrapolated to the physical point, our results favour a smallvalue of the electric charge radius, although the present accuracy is not suﬃcient to make adecisive statement in this regard.The paper is organized as follows: in section II we present the ensembles and operatorsused in our simulations; section III describes the methods we use to account for the pres-ence of excited-state contaminations in our data, and section IV the methods we employ toparameterize the form factor data on each lattice ensemble, while section V gives the resultsfor the extrapolation of our results to the continuum and inﬁnite-volume limits at the phys-ical pion mass. Our conclusions and a discussion of the results are contained in section VI.For completeness and ease of reference, we provide tables with the values of all measuredform factors in appendix A. Appendix B lists the priors we use to stabilize our two-stateﬁts, while the results of both dipole and z -expansion ﬁts on each ensemble are given inappendix C. Appendices D and E give the results of the extrapolations to the physical pointusing two variants of chiral perturbation theory. We consider the ratio G M ( Q ) /G E ( Q ) inappendix F. II. LATTICE SETUP

We use the CLS N f = 2 + 1 ensembles [54] that have been generated with non-perturbatively O ( a )-improved Wilson fermions [55, 56] and the tree-level improved L¨uscher-Weisz gauge action [57]. In order to prevent topological freezing [58], the ﬁelds obey openboundary conditions in time [59], with the exception of ensemble E250 which uses periodicboundary condition in the time direction. The reweighting factors needed to correct for thetreatment of the strange quark determinant during the gauge ﬁeld generation are obtainedusing the method of Ref. [60]. See Tab. I for a list of ensembles used in this work, whichcover the range of lattice spacings from 0 .

050 fm to 0 .

086 fm. Our setup in this work isidentical to that used in our paper on the isovector charges and momentum fractions of thenucleon [61], to which we refer the reader for further details.We obtain the matrix element of the vector current through the ratio [63] R J µ ( t, t s ; q ) = C J µ ( t, t s ; q ) C ( t s ; ) (cid:115) C ( t s − t ; − q ) C ( t, ) C ( t s ; ) C ( t s − t ; ) C ( t ; − q ) C ( t s ; − q ) , (1)where the nucleon two- and three-point-functions are given by C ( t ; p ) = Γ αβ (cid:88) x e − i px (cid:68) Ψ β ( x , t )Ψ α (0) (cid:69) , (2) C J µ ( t, t s ; q ) = Γ αβ (cid:88) x , y e i qy (cid:68) Ψ β ( x , t s ) J µ ( y , t )Ψ α (0) (cid:69) (3)in our setup, where the nucleon at the sink is at rest, i.e. for a momentum transfer q theinitial and ﬁnal nucleon states have momenta p (cid:48) = 0 , p = − q . (4)3 D β T /a L/a M π [MeV] M π L M N [GeV] N HP N LP N CFG t s [fm]H105 3.40 96 32 278(4) 3.90 1.020(18) 4076 48912 1019 1.0, 1.2, 1.4C101 3.40 96 48 223(3) 4.68 0.984(12) 2000 64000 2000 1.0, 1.2, 1.4S400 3.46 128 32 350(4) 4.34 1.123(15) 1725 27600 1725 1.1, 1.2, 1.4, 1.5N203 3.55 128 48 347(4) 5.42 1.105(13) 1540 24640 1540 1.0, 1.2, 1.3, 1.4S201 3.55 128 32 293(4) 3.05 1.097(21) 2092 66944 2092 1.0, 1.2, 1.3, 1.4N200 3.55 128 48 283(3) 4.42 1.053(14) 1697 20364 1697 1.0, 1.2, 1.3, 1.4D200 3.55 128 64 203(3) 4.23 0.960(13) 1019 32608 1019 1.0, 1.2, 1.3, 1.4E250 3.55 192 96 130(1) 4.04 0.928(11) 976 31232 244 1.0, 1.2, 1.3, 1.4N302 3.70 128 48 353(4) 4.28 1.117(15) 1177 18832 1177 1.0, 1.1, 1.2, 1.3J303 3.70 192 64 262(3) 4.24 1.052(17) 531 8496 531 1.0, 1.1, 1.2, 1.3TABLE I. Overview of ensembles used in this study. The quoted errors on the pion and nucleonmasses include the error from the scale setting [62]. N HP and N LP denote the number of high-precision (HP) and low-precision (LP) measurements on all of the N CFG conﬁgurations for eachvalue of the source-sink t s , respectively. For E250 N HP and N LP refer to the number of sourcesused for the two largest values of t s = 1 . , . t s = 1 . t s = 1 . Our interpolating operatorΨ α ( x ) = (cid:15) abc (cid:16) ˜ u Ta ( x ) Cγ ˜ d b ( x ) (cid:17) ˜ u c,α ( x ) , (5)for the proton is built using Gaussian-smeared [64] quark ﬁelds˜ q = (1 + κ G ∆) N G q , q = u, d, (6)using spatially APE-smeared [65] gauge links in the covariant Laplacian ∆ and tuning theparameters κ G and N G so that a smearing radius r G ∼ . γ )(1 + iγ γ ) . (7)Furthermore we use the improved conserved vector current J µ ( z ) c = 12 (cid:16) ¯ q ( z + ˆ µa )(1 + γ µ ) U µ ( z ) † q ( z ) − ¯ q ( z )(1 − γ µ ) U µ ( z ) q ( z + ˆ µa ) (cid:17) ,J µ ( z ) Imp . = 12 (cid:16) J µ ( z ) c + J µ ( z − ˆ µa ) c (cid:17) − a c cs V ∂ ν (cid:16) ¯ q ( z )[ γ µ , γ ν ] q ( z ) (cid:17) , (8)with the improvement coeﬃcient c cs V from [67]. The improvement term is implemented usingthe symmetric lattice derivative, i.e. ∂ ν φ ( x ) := φ ( x + a ˆ ν ) − φ ( x − a ˆ ν )2 a . (9)4o compute the three-point functions, we employ extended propagators in the “ﬁxed-sink”method, requiring additional inversions for each value of t s studied while allowing the mo-mentum transfer to be varied via a phase factor at the point of the current insertion [68]. Asin [61], we apply the truncated solver method with bias correction [69–71] to reduce the costof the inversions. The number of high-precision and low-precision measurements carried outon each gauge conﬁguration is indicated in Table I.For the polarization given in Eq. (7) the asymptotic value for the spectral decompositionof Eq. (1) reads R J ( t, t s ; q ) ≡ (cid:115) m N + E q E q G eﬀE ( Q , t, t s ) , (10a) R J ( t, t s ; q ) t, ( t s − t ) (cid:29) −−−−−−→ (cid:115) m N + E q E q G E ( Q ) , (10b)Re R J i ( t, t s ; q ) t, ( t s − t ) (cid:29) −−−−−−→ (cid:115) E q ( E q + M N ) G M ( Q ) (cid:15) ij q j , (10c)where G E 1 and G M are the isovector electric and magnetic Sachs form factors, respectively,with G E (0) = 1. The extraction of the magnetic form factor via Eq. (10c) amounts tosolving a system of linear equations, since, in general, several diﬀerent choices of q producethe same value of Q . Consequently there are several possibilities for obtaining an estimateof G M ( Q ). We use the following estimator, G eﬀM ( Q , t, t s ) = (cid:112) E q ( M N + E q ) q + q (cid:104) q Re R J − q Re R J (cid:105) , q (cid:54) = 0 ∨ q (cid:54) = 0 , (11)averaging over all momenta q contributing to the same Q . The resulting eﬀective formfactors for every source-sink separation for the ﬁrst non-zero momentum and a momentumclose to 0 . on the ensembles D200 and E250 are shown in Fig. 1.Unless stated otherwise, errors are computed using the jackknife method on binned datawith a bin size of two for all ensembles except E250, where the spacing between two analyzedconﬁgurations in terms of molecular dynamics time is twice as large compared to e.g. D200or C101 to begin with. For the conversion to physical units we use the lattice spacingdetermination of [62]. III. EXCITED-STATE SYSTEMATICS

Baryonic correlation functions suﬀer from a strong exponential growth of the relative sta-tistical noise when the distance in Euclidean time between operators is increased [72]. There-fore, for the typical source-sink separations in current lattice calculations of baryon structureobservables, it cannot be guaranteed that contributions from excited states are suﬃcientlysuppressed. Evidently, special care is required to avoid any bias from unwanted excited-statecontributions [73, 74]. Predominantly, two approaches have been widely adopted to addressthis problem: the summation method [39, 75–79] and multi-state ﬁts [39, 80–83]. In addition we have extracted the electric form factor from the spatial components of the current matrixelement, however the extracted values are less accurate compared to Eq. (10b). . . . . . . . G e ﬀ E ( t , t s ) . . . . . . G e ﬀ M ( t , t s ) t s =16 t s =18 t s =20 t s =220 . . . . . . G e ﬀ E ( t , t s ) . . . . . G e ﬀ M ( t , t s ) . . . . . . G e ﬀ E ( t , t s ) . . . . G e ﬀ M ( t , t s ) − − (cid:16) t − t s (cid:17) /a . . . . . . . G e ﬀ E ( t , t s ) − − (cid:16) t − t s (cid:17) /a . . . . . . . G e ﬀ M ( t , t s ) FIG. 1. Eﬀective form factors for ensemble D200 (upper panel) and E250 (lower panel). In eachpanel, the ﬁrst row corresponds to the smallest non-vanishing momentum in the given ensemble,i.e. Q = 0 . , . for D200 and E250, respectively, and the second row corresponds to Q ∼ . . For the four available source-sink separations t s , the eﬀective form factors aredisplayed as a function of the current insertion time t , oﬀset to the midpoint between nucleonsource and sink. The curves represent the two-state ﬁts in their respective ﬁt intervals. The grayband and black data point correspond to the estimate for the ground-state matrix element for thesummation and two-state method, respectively. The data points are displaced for better visibility. t/a . . . . . . . E { , } t/a . . . . . . c { , } × − FIG. 2. Energy levels (left) and overlap factors (right) extracted from the zero-momentum nucleontwo-point function on ensemble D200, for the ground state (blue) and the ﬁrst excited state (red).All quantities are given in lattice units.

While the former is (in its simplest incarnation) a straightforward method to apply, thelatter is more involved as one is forced to make speciﬁc assumptions and/or parameterchoices, regarding e.g. ﬁt windows, at various steps of the analysis. In this section we givedetails on our implementation of the two respective methods and discuss how errors relatedto methodology are incorporated in the ﬁnal results. The form factor values obtained withboth methods are collected in Appendix A for all ensembles.For the two-state ﬁts of the eﬀective form factors, we use priors obtained from an analysisof the two point functions. We ﬁt the two-point function with the ansatz C ( t, p ) = c ( p ) e − E ( p ) t + c ( p ) e − E ( p ) t , (12)and extract the energy gap between ground ( E ) and ﬁrst excited state ( E ), as well as theratio of the respective overlaps, i.e. ∆( p ) = E ( p ) − E ( p ) and ρ ( p ) = c ( p ) /c ( p ).In practice, the gaps and ratios depend on the choice of ﬁt ranges in time, especially thestarting timeslice. We therefore repeat the ﬁts for diﬀerent starting timeslices and obtainour best estimate as a weighted average over the region where the results have stabilized(see Fig. 2).For the nonlinear exponential ﬁts we use the VarPro method [84], which only needsinitial guesses for the energy levels. Monitoring the ground state energy, we ﬁnd that theextraction works well for all ensembles for momenta up to 1 GeV , see Fig. 3. The results,which are used in the subsequent analysis, are given in Appendix B for all ensembles.For the asymptotic limit of the ratio in Eq. (1) we obtain R as ( t, t s , Q ) = r (cid:110) ρ ( q )2 (cid:2) e − ∆( q )( t s − t ) − e − ∆( q ) t s (cid:3) + ρ ( )2 (cid:2) e − ∆( ) t − e − ∆( ) t s (cid:3)(cid:111) + r e − ∆( q ) t + r e − ∆( ) ( t s − t ) + r e − ∆( q ) t e − ∆( ) ( t s − t ) + . . . , (13)where r is proportional to the ground state matrix element G E/M ( Q ), and the last twoterms in the ﬁrst line come from the expansion of the two-point functions in Eq.(1). For each Note that r , r and r , even though the indexing might suggest otherwise, are not directly proportionalto the matrix element of the current. The ellipsis denotes terms with at least one exponential from the 2-and 3-point functions and further terms from the excitation spectrum. n . . . E n . . E FIG. 3. Ground state energy extracted in lattice units from two-state ﬁts to the nucleon two-pointfunction on ensembles D200 (left) and J303 (right), where the blue line describes the relativisticdispersion relation ( p = πL n ). The dashed, dashed-dotted lines indicate Q ≤ . , Q ≤ . , respectively.

16 18 20 22 t s /a S E

16 18 20 22 t s /a S M FIG. 4. Dependence of the summed ratios for the electric (left) and magnetic (right) eﬀective formfactors on the source-sink separation Eq. (15). Data is shown for the ﬁrst non-zero momentum onD200, i.e. Q = 0 .

089 GeV , together with a linear ﬁt using Eq. (16). value of the momentum q , the gap ∆( q ) and the terms proportional to ρ are universal,and we therefore proceed by ﬁtting the electric and magnetic eﬀective form factors simul-taneously. The ﬁts are performed up to a maximum transfer momentum of about 1 GeV .The ﬁts to the eﬀective form factors are stabilized using Gaussian priors for ∆ and ρ (seeAppendix B), whose central values are set to the results of the ﬁts to the two-point function.We monitor the impact of our particular choice for the priors on the extracted form factorvalues by varying the width of the priors in all ﬁts. To this end we multiply the errors of∆ and ρ by a factor between 1 and 5. The associated ﬁt results are labeled 1x, 2x, . . . ,5x to reﬂect their dependence on the prior width. In order for the prior to be eﬀective,we constrain the width to maximally half the mean value. The idea is to strike a balancebetween the statistical accuracy of our extracted values and any potential bias introducedby the priors. Therefore we choose the ﬁnal values for the two-state method to come fromthe analysis with the maximum prior width that (a) gives values compatible within errorswith all determinations based on a smaller width, and (b) maintains an acceptable errorincrease. We made a rather conservative choice for the latter, allowing the error to increase8y factors between 2 and 4 for prior widths 2x to 5x. For the ensembles in this work, itturns out that, for Q ≤ . the ﬁnal values come from an analysis with prior width5x.In addition to the above analysis, we perform ﬁts to the summed correlators. The sum-mation method takes advantage of the fact that in the ratios of Eq. (13), when summedover timeslices in between source and sink, the contributions from excited states are para-metrically suppressed. Summing Eq. (13) over t , omitting t skip timeslices at both ends, weobtain S ( t s ) = t s − t skip (cid:88) t = t skip R as ( t, t s , Q )= r (cid:104) − ρ ( q )2 e − ∆( q ) t s − ρ ( )2 e − ∆( ) t s (cid:105) ( t s + a − t skip ) a + (cid:104) r + r ρ ( q )2 (cid:105) ( e ∆( q )( a − t skip ) − e − ∆( q )( t s − t skip ) ) e a ∆( q ) − (cid:104) r + r ρ ( )2 (cid:105) ( e ∆( )( a − t skip ) − e − ∆( )( t s − t skip ) ) e a ∆( ) − . . . . (14)For the eﬀective form factors we may write the summed ratios as S E / M ( Q , t s ) = t s − a (cid:88) t =2 a G eﬀE / M ( Q , t, t s ) , (15)where in our analysis we use t skip = 2 a . The summation method crucially depends on thecomputation of observables for multiple source-sink separations. However, since the signalof the correlators deteriorates with larger time separations, we are limited in the number ofavailable t s . Our current data does not deviate signiﬁcantly from linear behavior (see Fig.4) for any of the ensembles and we therefore ﬁt Eq. (14) in the asymptotic limit with onlyground state contributions, i.e. S E / M ( Q , t s ) t s →∞ −−−→ C E / M ( Q ) + t s a G E / M ( Q ) + . . . (16)where C E / M is an irrelevant oﬀset.The results of the two methods are shown together in Fig. 5 for the near-physical ensembleE250. In order to assess systematics associated with excited-state eﬀects, we apply allsubsequent analyses to the data obtained by both methods, hence we distinguish betweensummation and explicit two-state ﬁts in the following. IV. PARAMETERIZATION OF THE Q DEPENDENCE

Since the magnetic moment µ is deﬁned as the intercept of G M and the electric andmagnetic radii are determined by the slope of the form factors at zero momentum transfer, µ = G M (0) = κ + 1 , (17a) (cid:104) r , M (cid:105) = − G E , M (0) ∂ G E , M ( Q ) ∂Q (cid:12)(cid:12)(cid:12)(cid:12) Q =0 , (17b) Note that the ellipsis denotes terms after the expansion of the ratio in Eq. (13) to more than one expo-nential in the nucleon energy. . . . . . . Q [GeV ] . . . G E ( Q ) . . . . . . Q [GeV ] G M ( Q ) FIG. 5. Comparison of the summation (blue circles) and two-state (red triangles) method forthe priors in Appendix B on ensemble E250. The black band, which corresponds to the param-eterization of [85], is displayed to enable a ﬁrst comparison to phenomenology. The continuumextrapolation of the lattice data is discussed in section V. a description of the Q dependence is necessary. In analogy to [39] we perform two analyses.In the ﬁrst analysis we parameterize the Q dependence using either a dipole or a z -expansionansatz [86] and subsequently perform chiral and continuum extrapolations. In the secondanalysis we use covariant Baryon Chiral Perturbation Theory (BChPT) [87–90] to ﬁt theavailable form factor data for all ensembles simultaneously. The latter approach combines thechiral extrapolation and the ﬁt to the Q dependence, i.e. without intermediary extractionof the radii using a separate ansatz of the Q behavior for each ensemble. We use theexpressions in [90] as they include vector meson degrees of freedom, e.g. ρ -mesons, in orderto extend the description of the form factors to larger values of Q .For the dipole ﬁts we use the ansatz G dipoleE / M ( Q ) = a E / M (cid:16) Q M (cid:17) , (18)where a E = 1. The dipole mass and a M are the ﬁt parameters. The dipole ﬁt is performedseparately for the electric and the magnetic form factor, thus the dipole mass is allowed tobe diﬀerent for G E and G M . The ﬁts are performed with cuts in Q of 0.6 GeV and 0.9GeV , and the corresponding results are collected in Appendix C.A model-independent description of the Q dependence of G E , M can be obtained byemploying the z -expansion [86]. The form factors may be decomposed as G E ( Q ) = ∞ (cid:88) k =0 a k z ( Q ) k , (19a) G M ( Q ) = ∞ (cid:88) k =0 b k z ( Q ) k , (19b)with z ( Q ) = (cid:112) t cut + Q − √ t cut − t (cid:112) t cut + Q + √ t cut − t . (19c)10 .

15 0 .

20 0 .

25 0 .

30 0 . M π [GeV] . . . . . h r E i [ f m ] .

15 0 .

20 0 .

25 0 .

30 0 . M π [GeV] . . . . . . h r M i [ f m ] .

15 0 .

20 0 .

25 0 .

30 0 . M π [GeV] µ FIG. 6. Dipole (red squares) and z -expansion results (blue circles) for the radii and the magneticmoment, obtained with the summation method for a momentum cut of 0.9 GeV . For bettervisibility, dipole and z -expansion data points corresponding to the same ensemble have been slightlyseparated from each other horizontally. The parameter t is the value of − Q which is mapped to z = 0. On each ensemble weset t cut = 4 M π , where M π denotes the pion mass of the respective ensemble. In general westabilize the ﬁts using priors, where the coeﬃcients are constrained using Gaussian distri-butions. The mean and width of these distributions are chosen such that coeﬃcients do notchange drastically with increasing order of the expansion. To that end we ﬁrst perform un-constrained ﬁts to order k = 1. For all subsequent ﬁts with orders k > z -expansionwe add Gaussian priors for the ﬁrst two coeﬃcients, centered around the means of the un-constrained ﬁt, with widths of 5 times the corresponding error estimate. For the remainingcoeﬃcients we choose Gaussian priors centered around zero with twice the maximum of theﬁtted coeﬃcients up to ﬁrst order . Throughout we always enforce G E (0) = 1. Addition-ally we performed ﬁts with stronger constraints on the large k behavior of the coeﬃcients a k coming from the fall-oﬀ of the Sachs form factors for large space like momentum transfer[91]. These conditions may be implemented in the form of sum rules [92] ∞ (cid:88) k = n k ( k − . . . ( k − n + 1) a k = 0 , n = 0 , , , . (20)For a given range in Q the optimal value for t is (see Ref. [86]) given by t opt0 ( Q ) = t cut (cid:16) − (cid:112) Q /t cut (cid:17) . (21) We have checked that this procedure gives consistent results with ﬁts that use linearly extrapolated valuesfor G M (0) and with ﬁts ﬁrst removing a residual monopole dependence. k . . . . h r E i k . . . . . h r M i FIG. 7. Dependence of the extracted electric and magnetic radii on the maximum order k of the z -expansion for t = 0 (circle) and t = − .

190 GeV (diamond) for the two-state method with(ﬁlled symbols) and without (empty symbols) constraints of Eq. (20) on ensemble D200. We perform the ﬁts with t = 0 GeV and t = t opt0 (0 . ) for momentum cuts of Q ≤ . and Q ≤ . ; the results are given in Appendix C. The dependence of theextracted electric and magnetic radii on the maximum order of the z -expansion is shownin Fig. 7. We ﬁnd that the ﬁts stabilize around k max of 5 for the ansatz with weakerassumptions about the bounds on the z -expansion parameters. Fits using the large- Q constraints of Eq. (20) in general converge more slowly, i.e. for larger values of k max ; onceplateaued, they however give compatible results. Since we are interested in the low- Q region, we choose the ﬁts without imposing Eq. (20). In Fig. 6 we show the extracted valuesfor data with Q ≤ . for the summation method. We see that the values obtainedfrom the diﬀerent parameterizations of the Q dependence are consistent, while the dipoleansatz gives smaller errors. We demand that at least four non-vanishing values of Q enterthe dipole and z -expansion ﬁts, respectively. This implies, for instance, that ensemble H105is excluded from the analysis when Q ≤ . is applied.In the direct approach based on covariant BChPT, we ﬁt the form factor data for both G E and G M on all ensembles simultaneously, and thus the intermediate parameterization ofthe Q dependence is avoided. While the ensembles are treated as statistically independent,we do take the correlation among diﬀerent Q and between G E and G M within an ensembleinto account. Statistical errors are derived from the covariance matrix estimated in theleast-square ﬁts that yield the results for charge radii, the magnetic momentum and otherﬁt parameters. Since this method uniﬁes the description of the Q and the M π dependence,it is presented in section V. V. CHIRAL AND CONTINUUM EXTRAPOLATION

In this section we present our chiral and continuum extrapolation, in order to arrive atresults for the isovector electric and magnetic radii, as well as for the magnetic moment. Weadopt two strategies, the ﬁrst of which follows up on the intermediate results of section IVand is presented in subsection V A. The second achieves a simultaneous description of the Q and the M π dependence of the results of section III and is presented in subsection V B. Themodel-averaging procedure used to arrive at our ﬁnal results is explained in subsection V C.12 . HBChPT extrapolation of the radii and the magnetic moment As detailed in Sec. IV, we perform ﬁts to the Q dependence of the Sachs form factorsusing a dipole and z -expansion ansatz. In this way we obtain estimates for the magneticmoment and the electromagnetic radii on each ensemble. For the chiral and continuumextrapolation of this data set we perform ﬁts based on Heavy Baryon Chiral PerturbationTheory (HBChPT) [26]. We apply momentum cuts to the data between Q ≤ . and Q ≤ . as well as a pion mass cut of M π ≤ .

29 GeV. We ﬁt the low-energy constants κ , E , B and B c , while setting the remaining constants to their phenomenological values ,i.e. g A = 1 . , F π = 0 .

092 GeV , c v = − .

26 GeV − , (22a)∆ = 0 .

294 GeV , g πN ∆ = 1 . . (22b)The results for these ﬁts extrapolated to the value of the pion mass in the isospin limitof QCD (134.8 MeV [93]) are given in Appendix D. We have analyzed variations of the ﬁtfunction, amending the formulae with lattice spacing terms of O ( a ) and/or including ﬁnitevolume eﬀects [95], i.e. κ = κ HB + a κ a + κ L M π (cid:16) − M π L (cid:17) e − M π L , (cid:104) r (cid:105) = (cid:104) r (cid:105) HB + a (cid:104) r (cid:105) a + (cid:104) r (cid:105) L e − M π L , (cid:104) r (cid:105) = (cid:104) r (cid:105) HB + a (cid:104) r (cid:105) a + (cid:104) r (cid:105) L e − M π L . (23)Here, the subscripts HB, a and L are used to distinguish the estimates in HBChPT fromthe coeﬃcients describing lattice artifacts and ﬁnite-volume eﬀects, respectively. Fits to theHBChPT expression fail if they are applied to the charge radii and the magnetic moment de-termined from the dipole ansatz. Even excluding the results from ensemble E250, which arehardest to accommodate in the ﬁt, does not improve the HBChPT description signiﬁcantly.For the z -expansion extraction on the other hand, a good ﬁt can be achieved for momentumcuts between Q ≤ . and Q ≤ . using HBChPT with and without latticeartifacts as parameterized in Eq. (23). We ﬁnd simultaneous ﬁts of ﬁnite lattice spacing andﬁnite volume dependence to be unstable. Let us stress again that we only include ensembleswith at least four data points remaining after momentum cuts are applied. This eﬀectivelylimits the lower bound on the cut in Q , which still allows for a chiral and continuum extrap-olation, to roughly 0 . . We also note that, even though for each ensemble a diﬀerentnumber of data points enters in the z -expansion, the relative weight in the HBChPT ﬁtdoes not reﬂect the density of available Q points at low momentum transfers. In this sensethe two-step process, ﬁrst performing z -expansion ﬁts and subsequently extrapolating usingHBChPT, masks the relative paucity of data points at small momentum transfer for someensembles. Note that when leaving g πN ∆ as a free parameter, the ﬁt is in general not improved. For the radii we choose a simpliﬁed ansatz due to the subtleties in unambiguously deﬁning ﬁnite volumecorrections, c.f. Ref. [94]. .

02 0 .

04 0 .

06 0 .

08 0 .

10 0 . . . . . . . . . µ .

02 0 .

04 0 .

06 0 .

08 0 .

10 0 . . . . . . . . . µ .

02 0 .

04 0 .

06 0 .

08 0 .

10 0 . . . . . . . h r M i [ f m ] .

02 0 .

04 0 .

06 0 .

08 0 .

10 0 . . . . . . . h r M i [ f m ] .

02 0 .

04 0 .

06 0 .

08 0 .

10 0 . M π [GeV ] . . . . . . . h r E i [ f m ] .

02 0 .

04 0 .

06 0 .

08 0 .

10 0 . M π [GeV ] . . . . . . . h r E i [ f m ] FIG. 8. HBChPT ﬁts to the radii and the magnetic moment, extracted via the z -expansion ofthe Sachs form factors determined with the summation (left panel) and two-state method (rightpanel), with Q ≤ . and M π ≤ .

28 GeV. Red points correspond to PDG values [96] for µ and (cid:104) r (cid:105) . For (cid:104) r (cid:105) we show the result of a reanalysis of available world data from Ref. [92], eitherbased exclusively on the Mainz/A1 measurement [7] (green diamond) or excluding it from worlddata (red diamond). The gray bar depicts our ﬁnal result of the model-averaged covariant BChPTanalysis, where the width indicates statistical error, and the black bar includes systematic eﬀects. B. Direct BChPT ﬁts

As an alternative to the intermediate determination of the Q dependence via z -expansionor dipole ﬁts, we perform direct ﬁts of the covariant BChPT expressions of [90] to our formfactor data. In this way we obtain a combined description of the Q and the M π dependence.The ﬁt depends linearly on the four LECs d , ˜ c , d x and G ρ [90]. It turns out that animportant advantage of this approach to extracting the electromagnetic radii compared tothe combined z -expansion and HBChPT analysis is its stability against considerably loweringthe momentum cuts applied.For the direct ﬁts we obtain results for various momentum cuts between Q ≤ . and Q ≤ . for both the summation method and the two-state method. We perform14he ﬁts with and without terms parameterizing the lattice spacing and/or ﬁnite volumedependence, G E ( Q ) = G E ( Q ) χ + a Q G a E + Q G L E e − M π L ,G M ( Q ) = G M ( Q ) χ + a G a M + κ L M π (cid:16) − M π L (cid:17) e − M π L + Q G L M e − M π L . (24)Fits leaving κ L as a free parameter are unstable, and we therefore ﬁx κ L to the value fromHBChPT [95] , i.e. κ L = − M N g A πF π . (25)As a cross-check, we perform ﬁts where the lattice artifacts enter multiplicatively, i.e. G E ( Q ) = G E ( Q ) χ (cid:0) a Q G a E + Q G L E e − M π L (cid:1) ,G M ( Q ) = G M ( Q ) χ (cid:0) a G a M + Q G L M e − M π L (cid:1) + κ L M π (cid:16) − M π L (cid:17) e − M π L . (26)In total we have six models, i.e. without any lattice artifact and either including dis-cretization or ﬁnite volume eﬀects, for the additive parameterization of Eq. (24) or themultiplicative one of Eq. (26), respectively.Within our statistical errors, discretization and ﬁnite-volume eﬀects are hardly signiﬁcant,implying that for most cuts the corresponding coeﬃcients are compatible with zero . Weperform simultaneous correlated ﬁts of G E ( Q ) and G M ( Q ) for each model with the givendata set and cuts applied for every ensemble. The errors for the ﬁt parameters are estimatedusing derivatives of the χ function. The direct method leads to more stable results incomparison to the two-step procedure in which ﬁts to the z -expansion are performed ﬁrst,followed by the chiral and continuum extrapolations using HBChPT. While the two methodsgive consistent results, the direct method has smaller errors, especially for the magnetic formfactor (see Fig. 8). Moreover, direct ﬁts allow for more stringent cuts in the momentumtransfer, and the analysis is more driven by data in the low Q region, where the radiiand magnetic moment are deﬁned. Finally, the inﬂuence of priors for the z -expansion iseliminated in the direct ﬁts. For these reasons we restrict the following presentation of theﬁnal results to the direct method, however noting that the same procedure applied to theHBChPT extractions give consistent results, albeit with larger errors.The quality of the direct covariant BChPT description is illustrated in Fig. 9, wherewe present a typical ﬁt to the extracted form factors for summation data, correspondingto the model without lattice artifacts or ﬁnite-volume corrections. The data is describedrather well over the ﬁt range in Q for all ensembles, already suggesting that, within ourstatistical accuracy, lattice artifacts are not discernible. Additionally, Fig. 9 illustrates thevery diﬀerent density of low- Q data points for each ensemble. For the magnetic moment,most recent lattice determinations [30, 35, 82, 97] lie below the experimental value. Forour most chiral ensemble, E250, we observe (top right panel of Fig. 9) that the directﬁt lies somewhat above the data points for the magnetic form factor, while still being Note that in Ref. [95] f is used instead of F π and that we are using the expression for the isovectormagnetic moment. Similar to the HBChPT ﬁts we ﬁnd simultaneous ﬁts of ﬁnite volume and lattice spacing dependence arenot stable for all applied cuts and we do not include them in our ﬁnal estimate. . . . . . . . . . . . G E ( Q ) E250 0 . . . . . . G M ( Q ) E2500 . . . . . . . . . . . G E ( Q ) D200 0 . . . . . . G M ( Q ) D2000 . . . . . . . . . . . G E ( Q ) C101 0 . . . . . . G M ( Q ) C1010 . . . . . . . . . . . G E ( Q ) J303 0 . . . . . . G M ( Q ) J3030 . . . . . . Q [GeV ] . . . . . G E ( Q ) H105 0 . . . . . . Q [GeV ] G M ( Q ) H105

FIG. 9. The summation-method data points for the Sachs form factors, and the blue band describ-ing the corresponding direct covariant BChPT ﬁt with momentum cut Q ≤ . , pion masscut 0 .

28 GeV and without lattice artifacts. The data point for G M (0) is obtained from a linearﬁt to the ratio of G M and G E and is not used in the direct covariant ChPT ﬁt. The ﬁt dependslinearly on the four LECs d , ˜ c , d x and G ρ (c.f. [39]). .

00 0 .

25 0 . Q [GeV ] G M ( Q ) G E ( Q ) D200C101J303 E250N200 0 . . . . M π [GeV] µ ( M π ) FIG. 10. The left panel shows the ratio of the magnetic and electric form factors for summa-tion method data. The right panel shows the direct covariant BChPT result for the pion-massdependence of the magnetic moment as a blue shaded area, together with the data points obtainedfrom a linear extrapolation to Q = 0 of the summation method data on the left panel. The blackdiamond corresponds to the PDG value [96]. Note that the displayed data points were not used inobtaining the covariant BChPT result. compatible within the uncertainties. Thus the possible presence of a non-negligible sourceof systematic error in calculations of the magnetic form factor at the physical point meritsfurther investigation.From phenomenological dipole ﬁts to experimental data, the ratio of the electric andmagnetic form factor is known to show a rather constant behavior over a large range of Q . Indeed we observe similar behavior in our data (see Fig. 10), where the two-state andsummation data are rather ﬂat with the exception of J303 showing signs of a light upwardslope. This pattern motivates a linear extrapolation of the ratio for the magnetic momentup to Q ≤ . ( Q ≤ .

29 GeV for E250), where the results are given in Tab. XIIand shown in Fig. 10 (right). It is reassuring to see that our ﬁt, while not using these points,does reproduce the estimates from a linear extrapolation of the ratio rather well. C. Model average and ﬁnal result

As we have no a priori preference for the results from the summation method or two-stateﬁts, we will treat both data sets on an equal footing. We obtain our ﬁnal estimates andtotal errors from averages over ﬁt models and kinematic cuts using weights derived fromthe Akaike Information Criterion (AIC) [98, 99]. In this context the momentum and pionmass cuts applied can be reinterpreted in terms of a model selection problem [100]. Onemay introduce for each would-be-cut data point an additional ﬁt parameter that matches therespective data point exactly, thus giving no contribution to the χ . The corresponding datapoint is eﬀectively excluded from the least-square ﬁt, while the model weight is decreasedvia the penalty term for additional parameters in the AIC. The AIC readsAIC i = χ ,i + 2 n f + 2 n c , (27)where χ ,i denotes the minimum of the weighted least square, for the i-th model, n f thenumber of ﬁt parameters, and n c the number of cut data points. In this way we obtain acriterion that takes the goodness of ﬁt into account. At the same time it penalizes increasing17he number of ﬁt parameters while it favors including more actual data points [39]. For theweighting of diﬀerent models on the same input data set we use w AIC i = e − AIC i (cid:80) j e − AIC j , (28)i.e. we normalize the AIC obtained for all models for summation and two-state data sep-arately. Finally, we apply a ﬂat weight function to the estimates from summation andtwo-state ﬁts. We adopt the procedure from [101], which we brieﬂy sketch in the following,for estimating the systematic and statistical error of the model-averaged values. To that endwe treat the model-averaged estimate as a random variable with the following cumulativedistribution function (CDF) P x ( y ) = y (cid:90) −∞ n (cid:88) i w i N ( y (cid:48) ; x i , σ i ) dy (cid:48) (29)i.e. the weighted sum of Gaussian distributions where the mean x i and variance σ i is givenby the best estimate and ﬁt error of each model, and the weight w i is obtained as explainedabove. This eﬀectively smoothens the otherwise rugged distribution of model postdictionsand allows for a more robust estimate of the distribution parameters (see Fig. 11). The ﬁnalvalue and total error are easily read oﬀ from the distribution in Eq. (29) as the median,and the 1- σ percentiles, respectively. Under the assumptions that a rescaling of all errors .

00 3 .

25 3 .

50 3 .

75 4 .

00 4 .

25 4 . κ . . . . . . P κ ( y ) . . . . . . . h r i [fm ] . . . . . . P h r E i ( y ) .

50 0 .

55 0 .

60 0 .

65 0 .

70 0 .

75 0 . h r i [fm ] . . . . . . P h r M i ( y ) FIG. 11. Cumulative distribution function of all ﬁtted models, where dash-dotted and short-dashed lines indicate median and 68% percentiles, respectively. i only aﬀects the statistical error but not the systematic one, we can further separate thestatistical and systematic errors, c.f. [101].In our previous work based on N f = 2 ensembles [39], we used the spread in the centralvalues as an estimate of systematic errors. While this procedure is robust, it is also veryconservative and susceptible to overestimating the true error due to systematics. Therefore,in order to not be overly conservative and still be able to incorporate systematic errors in arobust way, we adopt the above model averaging procedure using AIC weights and obtainas our ﬁnal results κ = 3 . ± . ± . , (cid:104) r (cid:105) = 0 . ± . ± .

022 fm , (cid:104) r (cid:105) = 0 . ± . ± .

011 fm , (30)where the ﬁrst and second errors refer to statistical uncertainty and the total systematicerror, respectively.One may even proceed further and estimate the individual contributions for every varia-tion to the total systematic error. That is achieved by building the CDF in Eq. (27) not overall variations but rather ﬁrst iterating over a particular feature, e.g. a momentum cut, andperforming the analysis for every variant of that feature separately. From this we then builda secondary CDF like Eq. (27) and extract the variation-speciﬁc systematic error. Repeatingthis analysis for all variations we obtain the following systematic error budget, δκ = 0 . exc ± . artifacts ± . Q ± . m π ± . method (= 0 . ,δ (cid:104) r (cid:105) = 0 . exc ± . artifacts ± . Q ± . m π ± . method fm (= 0 .

022 fm ) ,δ (cid:104) r (cid:105) = 0 . exc ± . artifacts ± . Q ± . m π ± . method fm (= 0 .

012 fm ) . (31)We note that, due to correlations, the individual terms added in quadrature (as indicated oneach line by the number in brackets) need not exactly reproduce the total error of Eq. (30).For the magnetic moment and the electric radius, the dominant source of systematic errorremains the excited state contribution, while for the magnetic radius all systematic eﬀectsconsidered here have comparable size. Moreover, in the current analysis the magnetic radiusis least aﬀected by systematic errors.In Fig. 12 we compare our work to recent lattice determinations and to the phenomeno-logical values for the isovector magnetic moment and the isovector electromagnetic radii.While we postpone the comparison of our results with phenomenological determinations ofthe radii to section VI, we remark that our value for the magnetic moment is in good agree-ment with the experimentally precisely known diﬀerence of proton and neutron magneticmoments. As for the comparison with other lattice calculations, we note that our estimate iscompatible with the determinations from [97, 102], while there is a sizeable diﬀerence to thevalues from [30, 35, 82]. We stress that the diﬀerence is not related to the issue of preferringdirect ﬁts to the form factor data over the more conventional route via the z -expansion, asthe latter shows a trend to higher values for the radius for our data. Our error estimates forthe statistical and systematic errors are comparable in size with the other lattice determina-tions. For the isovector magnetic moment we see good agreement with phenomenology and[35, 102]. We note that the missing data point for Q = 0 complicates the extraction of thelow- Q observables in most recent lattice determinations. Especially the z -expansion ﬁts,at least for orders n ≥

2, tend to overﬁt the dependence of the form factor at low Q . Inorder to remedy this, either priors are introduced or mock data points at Q = 0, e.g. linear19 . . µ This workPhenomenologyPNDME20LHPC17PACS18ETMC18ETMC20 0 . . h r E i [fm] . . . h r M i [fm] FIG. 12. Comparison of our best estimate (downward-pointing triangle) for the isovector quantities µ , (cid:104) r E (cid:105) , (cid:104) r M (cid:105) to other lattice calculations, i.e. PNDME [82] (circle), ETMC [30, 97] (diamond),PACS [35] (upward-pointing triangle), LHPC [102] (square). The phenomenological value for µ isderived form the PDG values [96]. The two data points for (cid:104) r E (cid:105) are derived from CODATA 2018(cross) or Mainz/A1 [7] (square) values for the proton electric charge radius, respectively, whilethe values for the neutron are taken from [96] for both. The two data points for (cid:104) r M (cid:105) depict thevalues inferred from the proton results taken from the reanalysis of [92] including only data from[7] (square) or excluding the Mainz/A1 data set from the analysis (cross), while taking the valuesfor the neutron magnetic radius from PDG [96] for both. For ease of comparison, the blue bandrepresents our ﬁnal result with the full uncertainty, with the light band indicating the statisticalerror. extrapolations of the ratio of the isovector form factors, are used to stabilize the description.We note that the direct approach, in this sense, has less freedom and by itself allows forconsiderably less variation in the form factors at low Q (see Fig. 9). We believe this to beresponsible, in large part, for the small errors we ﬁnd in the isovector magnetic radius. VI. CONCLUSIONS

We have calculated the isovector electromagnetic form factors of the nucleon in latticeQCD with dynamical up, down and strange quarks. The electromagnetic radii and themagnetic moment have been extracted accounting for systematic eﬀects due to excited states,ﬁnite volume and non-zero lattice spacing. Our ﬁnal estimates are listed in Eq. (30), with adetailed systematic error budget given in Eq. (31).As an important benchmark, we reproduced the experimental value of the magneticmoment with an overall precision of 3.6%. The precision of the present calculation is signiﬁ-cantly higher than that of our earlier study in two-ﬂavor QCD [39], especially concerning themagnetic properties. For the isovector electric charge radius, our result is in good agreementwith the phenomenological estimate inferred from the 2018 CODATA recommended value ofthe proton radius . By contrast, after adding all errors in quadrature, we ﬁnd a 2.4 σ tensionwith the result from ep -scattering [7]. For the isovector magnetic radius, on the other hand, The central value for the latter is very close, in comparison to its uncertainty of 2.3 ‰ , to that extractedfrom muonic hydrogen [9], which is yet 4.9 times more precise. ep -scattering based determination [7],and exhibits a sizeable tension with the other collected world data [92]. For ease of com-parison, we translate our estimate for the isovector (cid:104) r (cid:105) into a result for the proton radiuswith the help of the experimental determination of the (squared) neutron charge radius, (cid:104) r (cid:105) = − . [96]. After combining all errors we obtain (cid:104) r (cid:105) / = 0 . Q , i.e. at large volumes and at physical pion massare necessary. We plan to extend our analysis to such ensembles as they become available.A promising strategy to further stabilize multi-state ﬁts of the three-point functions wouldbe to perform a dedicated study of the excitation spectrum. Moreover, an analysis of theexcited-state contributions in chiral eﬀective theory, as has been done for the axial formfactor [103], and the expression for the ﬁnite volume dependence, would be highly desirableto improve the assessment of the related systematic errors. Acknowledgments [1] J. Ashman et al. (European Muon Collaboration), Phys. Lett.

B206 , 364 (1988), [,340(1987)].[2] C. A. Aidala, S. D. Bass, D. Hasch, and G. K. Mallot, Rev. Mod. Phys. , 655 (2013),1209.2803.[3] X.-S. Chen, W.-M. Sun, X.-F. Lu, F. Wang, and T. Goldman, Phys. Rev. Lett. , 062001(2009), 0904.0321.[4] X. Ji, F. Yuan, and Y. Zhao, Nature Rev. Phys. , 27 (2021), 2009.01291.[5] R. Pohl et al., Nature , 213 (2010).

6] J.-P. Karr, D. Marchand, and E. Voutier, Nature Rev. Phys. , 601 (2020).[7] J. C. Bernauer et al. (A1 Collaboration), Phys. Rev. Lett. , 242001 (2010), 1007.5076.[8] P. J. Mohr, B. N. Taylor, and D. B. Newell, Rev. Mod. Phys. , 1527 (2012), 1203.5425.[9] A. Antognini et al., Science , 417 (2013).[10] C. E. Carlson, Prog. Part. Nucl. Phys. , 59 (2015), 1502.05314.[11] A. Beyer, J. Alnis, K. Khabarova, A. Matveev, C. G. Parthey, D. C. Yost, R. Pohl, T. Udem,T. W. H¨ansch, and N. Kolachevsky, Annalen Phys. , 671 (2013).[12] S. Thomas, H. Fleurbaey, S. Galtier, L. Julien, F. Biraben, and F. Nez, Annalen Phys. ,1800363 (2019), 1903.04252.[13] H. Fleurbaey, S. Galtier, S. Thomas, M. Bonnaud, L. Julien, F. Biraben, F. Nez, M. Abgrall,and J. Gu´ena, Phys. Rev. Lett. , 183001 (2018), 1801.08816.[14] M. Mihoviloviˇc et al., Phys. Lett. B771 , 194 (2017), 1612.06707.[15] M. Mihoviloviˇc et al. (2019), 1905.11182.[16] A. H. Gasparian (PRad Collaboration), JPS Conf. Proc. , 020052 (2017).[17] W. Xiong et al., Nature , 147 (2019).[18] A. Gasparian et al. (PRad Collaboration) (2020), 2009.10510.[19] S. Grieser, D. Bonaventura, P. Brand, C. Hargens, B. Hetz, L. Leßmann, C. Westph¨alinger,and A. Khoukaz, Nucl. Instrum. Meth. A , 120 (2018), 1806.05409.[20] R. Gilman et al. (MUSE Collaboration) (2017), 1709.09753.[21] C. Dreisbach et al. (COMPASS++/AMBER working group), PoS DIS2019 , 222 (2019).[22] J. Krauth et al., in (2017), pp. 95–102, 1706.00696.[23] M. Hoferichter, B. Kubis, J. Ruiz de Elvira, H. W. Hammer, and U. G. Meißner, Eur. Phys.J. A , 331 (2016), 1609.06722.[24] H.-W. Hammer and U.-G. Meißner, Sci. Bull. , 257 (2020), 1912.03881.[25] M. A. Belushkin, H. W. Hammer, and U. G. Meißner, Phys. Rev. C , 035202 (2007),hep-ph/0608337.[26] M. G¨ockeler, T. R. Hemmert, R. Horsley, D. Pleiter, P. E. L. Rakow, A. Sch¨afer, andG. Schierholz (QCDSF Collaboration), Phys. Rev. D71 , 034508 (2005), hep-lat/0303019.[27] S. N. Syritsyn et al., Phys. Rev.

D81 , 034507 (2010), 0907.4194.[28] C. Alexandrou, M. Constantinou, S. Dinter, V. Drach, K. Jansen, C. Kallidonis, and G. Kout-sou, Phys. Rev.

D88 , 014509 (2013), 1303.5979.[29] C. Alexandrou, M. Constantinou, K. Hadjiyiannakou, K. Jansen, C. Kallidonis, G. Koutsou,and A. Vaquero Aviles-Casco, Phys. Rev.

D96 , 034503 (2017), 1706.00469.[30] C. Alexandrou, S. Bacchio, M. Constantinou, J. Finkenrath, K. Hadjiyiannakou, K. Jansen,G. Koutsou, and A. V. A. Casco (2018), 1812.10311.[31] P. E. Shanahan, A. W. Thomas, R. D. Young, J. M. Zanotti, R. Horsley, Y. Nakamura,D. Pleiter, P. E. L. Rakow, G. Schierholz, and H. St¨uben, Phys. Rev.

D90 , 034502 (2014),1403.1965.[32] P. E. Shanahan, A. W. Thomas, R. D. Young, J. M. Zanotti, R. Horsley, Y. Nakamura,D. Pleiter, P. E. L. Rakow, G. Schierholz, and H. St¨uben (CSSM, QCDSF/UKQCD Collab-orations), Phys. Rev.

D89 , 074511 (2014), 1401.5862.[33] T. Bhattacharya, S. D. Cohen, R. Gupta, A. Joseph, H.-W. Lin, and B. Yoon, Phys. Rev.

D89 , 094502 (2014), 1306.5435.[34] T. Yamazaki, Y. Aoki, T. Blum, H.-W. Lin, S. Ohta, S. Sasaki, R. Tweedie, and J. Zanotti,Phys. Rev.

D79 , 114505 (2009), 0904.2039.

35] E. Shintani, K.-I. Ishikawa, Y. Kuramashi, S. Sasaki, and T. Yamazaki, Phys. Rev.

D99 ,014510 (2019), 1811.07292.[36] K.-I. Ishikawa, Y. Kuramashi, S. Sasaki, N. Tsukamoto, A. Ukawa, and T. Yamazaki (PACSCollaboration), Phys. Rev.

D98 , 074510 (2018), 1807.03974.[37] J. R. Green, J. W. Negele, A. V. Pochinsky, S. N. Syritsyn, M. Engelhardt, and S. Krieg,Phys. Rev.

D90 , 074507 (2014), 1404.4029.[38] A. J. Chambers et al. (QCDSF, UKQCD, CSSM Collaborations), Phys. Rev.

D96 , 114509(2017), 1702.01513.[39] S. Capitani, M. Della Morte, D. Djukanovic, G. von Hippel, J. Hua, B. J¨ager, B. Knippschild,H. B. Meyer, T. D. Rae, and H. Wittig, Phys. Rev.

D92 , 054511 (2015), 1504.04628.[40] B. D. Milbrath et al. (Bates FPP Collaboration), Phys. Rev. Lett. , 452 (1998), [Erratum:Phys. Rev. Lett.82,2221(1999)], nucl-ex/9712006.[41] M. K. Jones et al. (Jeﬀerson Lab Hall A Collaboration), Phys. Rev. Lett. , 1398 (2000),nucl-ex/9910005.[42] O. Gayou et al. (Jeﬀerson Lab Hall A Collaboration), Phys. Rev. Lett. , 092301 (2002),nucl-ex/0111010.[43] V. Punjabi et al., Phys. Rev. C71 , 055202 (2005), [Erratum: Phys. Rev.C71,069902(2005)],nucl-ex/0501018.[44] A. J. R. Puckett et al., Phys. Rev. Lett. , 242301 (2010), 1005.3419.[45] A. J. R. Puckett et al., Phys. Rev.

C85 , 045203 (2012), 1102.5737.[46] M. N. Rosenbluth, Phys. Rev. , 615 (1950).[47] R. C. Walker et al., Phys. Rev. D49 , 5671 (1994).[48] L. Andivahis et al., Phys. Rev.

D50 , 5491 (1994).[49] I. A. Qattan et al., Phys. Rev. Lett. , 142301 (2005), nucl-ex/0410010.[50] P. A. M. Guichon and M. Vanderhaeghen, Phys. Rev. Lett. , 142303 (2003), hep-ph/0306007.[51] J. Arrington, W. Melnitchouk, and J. A. Tjon, Phys. Rev. C76 , 035205 (2007), 0707.1861.[52] A. Afanasev, P. G. Blunden, D. Hasell, and B. A. Raue, Prog. Part. Nucl. Phys. , 245(2017), 1703.03874.[53] J. C. Bernauer et al. (2020), 2008.05349.[54] M. Bruno et al., JHEP , 043 (2015), 1411.3982.[55] B. Sheikholeslami and R. Wohlert, Nucl. Phys. B259 , 572 (1985).[56] J. Bulava and S. Schaefer, Nucl. Phys.

B874 , 188 (2013), 1304.7093.[57] M. L¨uscher and P. Weisz, Commun. Math. Phys. , 59 (1985), [Erratum: Commun. Math.Phys.98,433(1985)].[58] S. Schaefer, R. Sommer, and F. Virotta (ALPHA Collaboration), Nucl. Phys. B845 , 93(2011), 1009.5228.[59] M. L¨uscher and S. Schaefer, JHEP , 036 (2011), 1105.4749.[60] D. Mohler and S. Schaefer, Phys. Rev. D , 074506 (2020), 2003.13359.[61] T. Harris, G. von Hippel, P. Junnarkar, H. B. Meyer, K. Ottnad, J. Wilhelm, H. Wittig, andL. Wrang (2019), 1905.01291.[62] M. Bruno, T. Korzec, and S. Schaefer, Phys. Rev. D , 074504 (2017), 1608.08900.[63] C. Alexandrou, T. Korzec, G. Koutsou, M. Brinet, J. Carbonell, V. Drach, P.-A. Harraud,and R. Baron (European Twisted Mass Collaboration), PoS LATTICE2008 , 139 (2008),0811.0724.[64] S. G¨usken, U. L¨ow, K. H. Mutter, R. Sommer, A. Patel, and K. Schilling, Phys. Lett.

B227 ,

66 (1989).[65] M. Albanese et al. (APE Collaboration), Phys. Lett.

B192 , 163 (1987).[66] G. M. von Hippel, B. J¨ager, T. D. Rae, and H. Wittig, JHEP , 014 (2013), 1306.1440.[67] A. G´erardin, T. Harris, and H. B. Meyer, Phys. Rev. D99 , 014519 (2019), 1811.08209.[68] G. Martinelli and C. T. Sachrajda, Nucl. Phys.

B316 , 355 (1989).[69] G. S. Bali, S. Collins, and A. Sch¨afer, Comput. Phys. Commun. , 1570 (2010), 0910.3970.[70] T. Blum, T. Izubuchi, and E. Shintani, Phys. Rev.

D88 , 094503 (2013), 1208.4349.[71] E. Shintani, R. Arthur, T. Blum, T. Izubuchi, C. Jung, and C. Lehner, Phys. Rev.

D91 ,114511 (2015), 1402.0244.[72] G. Lepage, in

Theoretical Advanced Study Institute in Elementary Particle Physics (1989),pp. 97–120.[73] J. Green, PoS

LATTICE2018 , 016 (2018), 1812.10574.[74] K. Ottnad, in (2020), 2011.12471.[75] S. Gusken, Nucl. Phys. Proc. Suppl. , 361 (1990).[76] L. Maiani, G. Martinelli, M. Paciello, and B. Taglienti, Nucl. Phys. B , 420 (1987).[77] T. Doi, M. Deka, S.-J. Dong, T. Draper, K.-F. Liu, D. Mankame, N. Mathur, and T. Streuer,Phys. Rev. D , 094503 (2009), 0903.3232.[78] J. Bulava, M. Donnellan, and R. Sommer, JHEP , 140 (2012), 1108.3774.[79] S. Capitani, M. Della Morte, G. von Hippel, B. J¨ager, A. J¨uttner, B. Knippschild, H. Meyer,and H. Wittig, Phys. Rev. D , 074502 (2012), 1205.0180.[80] B. Yoon et al., Phys. Rev. D , 074508 (2017), 1611.07452.[81] C. Chang et al., Nature , 91 (2018), 1805.12130.[82] Y.-C. Jang, R. Gupta, H.-W. Lin, B. Yoon, and T. Bhattacharya, Phys. Rev. D , 014507(2020), 1906.07217.[83] R. Gupta, Y.-C. Jang, B. Yoon, H.-W. Lin, V. Cirigliano, and T. Bhattacharya, Phys. Rev.D , 034503 (2018), 1806.09006.[84] G. H. Golub and V. Pereyra, SIAM Journal on Numerical Analysis , 413 (1973).[85] Z. Ye, J. Arrington, R. J. Hill, and G. Lee, Phys. Lett. B , 8 (2018), 1707.09063.[86] R. J. Hill and G. Paz, Phys. Rev. D82 , 113005 (2010), 1008.4619.[87] B. Kubis and U.-G. Meißner, Nucl. Phys. A , 698 (2001), hep-ph/0007056.[88] T. Fuchs, J. Gegelia, and S. Scherer, J. Phys. G , 1407 (2004), nucl-th/0305070.[89] M. R. Schindler, J. Gegelia, and S. Scherer, Eur. Phys. J. A , 1 (2005), nucl-th/0509005.[90] T. Bauer, J. C. Bernauer, and S. Scherer, Phys. Rev. C86 , 065206 (2012), 1209.3872.[91] G. P. Lepage and S. J. Brodsky, Phys. Rev.

D22 , 2157 (1980).[92] G. Lee, J. R. Arrington, and R. J. Hill, Phys. Rev.

D92 , 013013 (2015), 1505.01489.[93] S. Aoki et al., Eur. Phys. J. C , 112 (2017), 1607.00299.[94] B. C. Tiburzi, Phys. Rev. D , 014510 (2008), 0710.3577.[95] S. R. Beane, Phys. Rev. D , 034507 (2004), hep-lat/0403015.[96] P. Zyla et al. (Particle Data Group), PTEP , 083C01 (2020).[97] C. Alexandrou, K. Hadjiyiannakou, G. Koutsou, K. Ottnad, and M. Petschlies, Phys. Rev.D , 114504 (2020), 2002.06984.[98] H. Akaike, B. N. Petrov, and F. Csaki, Second international symposium on informationtheory (1973).[99] H. Akaike, IEEE Transactions on Automatic Control , 716 (1974).[100] W. I. Jay and E. T. Neil, Bayesian model averaging for analysis of lattice ﬁeld theory results (2020), 2008.01069. , 034504 (2018), 1711.11385.[103] O. B¨ar, Phys. Rev. D , 034515 (2020), 1912.05873.[104] R. G. Edwards and B. Joo (SciDAC, LHPC, UKQCD Collaborations), Nucl. Phys. Proc.Suppl. , 832 (2005), [,832(2004)], hep-lat/0409003.[105] M. L¨uscher and S. Schaefer, Comput. Phys. Commun. , 519 (2013), 1206.2809.[106] D. Djukanovic, Comput. Phys. Commun. , 106950 (2020), 1603.01576. Appendix A: Form factor data

In this appendix, we present the results of extracting the isovector electromagnetic formfactors of the nucleon with either the summation method or the two-state method, bothdescribed in section III, for every gauge ensemble listed in Table I.The summation-method results quoted below are obtained from a ﬁt using Eq. (16) as ﬁtansatz. For the two-state-method, we perform ﬁts for diﬀerent prior widths δ , where we useinteger factors between one and ﬁve multiplying the initial error estimate. However, for theprior to have an eﬀect, we limit the width to at most 50% of the central value. We arrive atthe numbers listed in the tables below using the data for which (a) all values up to a givenwidth are within 2 σ , and (b) the error does not increase by more than a factor of (1 + δ ),i.e. factors of [2 , , ,

4] for prior widths [2 x, x, x, x ]. Q [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.089 0.803(16) 0.803(11) 3.49(22) 3.36(11)0.174 0.671(20) 0.663(12) 2.71(16) 2.688(77)0.255 0.556(25) 0.558(21) 2.40(15) 2.381(78)0.334 0.473(28) 0.465(29) 2.09(15) 2.04(10)0.410 0.420(27) 0.389(42) 1.70(13) 1.69(14)0.484 0.363(31) 0.359(34) 1.62(14) 1.61(12)0.624 0.294(43) 0.272(39) 1.52(17) 1.14(16)0.692 0.301(43) 0.256(27) 1.16(16) 1.11(11)0.757 0.303(58) 0.207(20) 1.13(22) 0.949(99)0.821 0.319(60) 0.197(20) 1.21(23) 0.93(10)0.884 0.158(81) 0.139(29) 0.44(34) 0.74(13) G E and G M for ensemble D200. [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.087 0.7921(91) 0.755(22) 3.53(13) 3.35(21)0.171 0.668(13) 0.638(27) 3.00(10) 2.81(10)0.252 0.559(15) 0.501(53) 2.47(11) 2.35(15)0.329 0.459(25) 0.420(47) 2.26(11) 2.08(15)0.404 0.420(18) 0.381(35) 1.99(11) 1.79(14)0.476 0.363(23) 0.324(42) 1.686(90) 1.53(16)0.615 0.323(26) 0.276(13) 1.58(14) 0.72(31)0.682 0.295(30) 0.246(16) 1.58(14) 1.35(21)0.747 0.266(37) 0.245(34) 1.42(15) 1.11(18)0.810 0.271(37) 0.203(35) 1.27(14) 1.00(17)0.872 0.320(77) 0.155(48) 1.12(32) 0.80(20) G E and G M for ensemble C101. Q [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.193 0.632(32) 0.590(58) 2.75(26) 2.58(31)0.370 0.453(33) 0.370(88) 2.32(21) 2.32(40)0.536 0.275(49) 0.233(69) 1.73(22) 1.58(27)0.692 0.142(72) 0.206(44) 1.41(39) 0.89(34)0.840 0.209(62) 0.166(72) 0.82(24) 0.82(27)0.980 0.168(55) 0.140(39) 0.63(19) 0.72(19)1.244 0.20(12) 0.111(45) 0.10(47) 0.54(21) G E and G M for ensemble H105. Q [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.156 0.723(19) 0.7293(75) 3.17(18) 2.999(63)0.303 0.561(22) 0.5656(77) 2.66(16) 2.462(49)0.441 0.470(31) 0.4685(92) 2.15(14) 1.957(51)0.573 0.315(44) 0.381(12) 1.57(18) 1.695(52)0.699 0.295(38) 0.330(11) 1.42(16) 1.464(47)0.820 0.331(46) 0.315(12) 1.56(20) 1.329(55)1.048 0.36(11) 0.206(18) 0.77(35) 0.925(72)1.156 0.22(10) 0.172(21) 0.60(41) 0.825(85) G E and G M for ensemble N200. [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.156 0.6946(80) 0.694(17) 3.033(80) 3.061(57)0.304 0.512(11) 0.515(23) 2.428(69) 2.408(59)0.444 0.394(14) 0.419(20) 2.040(77) 1.980(62)0.577 0.322(19) 0.329(21) 1.679(89) 1.54(13)0.705 0.253(17) 0.266(15) 1.485(80) 1.11(20)0.827 0.190(22) 0.208(18) 1.272(94) 0.91(22)1.059 0.172(40) 0.162(24) 1.01(16) 0.897(89)1.170 0.074(41) 0.059(54) 0.66(16) 0.50(17)1.277 0.165(63) 0.147(30) 1.05(29) 0.75(12)1.381 0.144(69) 0.144(15) 0.80(30) 0.723(88) G E and G M for ensemble N203. Q [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.256 0.575(21) 0.613(12) 2.62(15) 2.435(80)0.491 0.352(24) 0.417(16) 1.78(12) 1.851(71)0.708 0.149(36) 0.276(29) 1.32(16) 1.39(13)0.912 0.135(57) 0.212(27) 1.12(25) 1.16(11)1.105 0.123(55) 0.165(36) 1.30(25) 0.95(12)1.287 0.092(77) 0.103(34) 0.56(33) 0.56(14) G E and G M for ensemble N302. Q [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.146 0.722(24) 0.715(18) 3.32(23) 2.92(12)0.284 0.532(25) 0.544(25) 2.65(17) 2.36(10)0.415 0.445(35) 0.453(18) 2.54(19) 2.12(12)0.539 0.348(49) 0.373(18) 2.17(23) 1.753(84)0.658 0.347(41) 0.326(18) 1.95(19) 1.586(71)0.772 0.360(52) 0.288(14) 1.86(23) 1.381(64)0.988 0.310(94) 0.241(18) 1.64(40) 1.196(85)1.090 0.40(12) 0.181(39) 2.09(50) 0.93(14)1.190 0.08(14) 0.117(50) 0.80(59) 0.64(20)1.287 0.10(16) 0.112(32) 0.53(63) 0.46(16) G E and G M for ensemble J303. Q [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.340 0.448(54) 0.522(17) 2.06(32) 2.27(10)0.643 0.405(60) 0.370(12) 1.38(23) 1.588(73)0.919 0.59(23) 0.238(29) 2.31(69) 1.11(12)1.175 0.44(29) 0.114(33) 2.3(1.2) 0.47(15) G E and G M for ensemble S201. [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.246 0.604(21) 0.607(14) 2.84(16) 2.746(79)0.471 0.422(23) 0.412(18) 2.26(13) 2.090(58)0.680 0.273(38) 0.312(17) 1.78(17) 1.616(55)0.877 0.138(63) 0.204(29) 0.94(28) 1.20(11)1.062 0.113(53) 0.181(18) 1.04(26) 1.128(54)1.239 0.106(85) 0.112(35) 0.99(37) 0.772(99) G E and G M for ensemble S400. Q [GeV ] G E (sum) G E (two-state) G M (sum) G M (two-state)0.040 0.885(22) 0.889(18) 3.68(47) 3.66(31)0.079 0.812(31) 0.817(17) 3.43(39) 3.16(18)0.117 0.732(36) 0.718(42) 2.88(33) 2.80(27)0.155 0.658(41) 0.675(35) 3.34(36) 3.02(22)0.191 0.633(36) 0.641(26) 2.86(30) 2.56(16)0.227 0.554(39) 0.542(76) 2.19(24) 2.28(26)0.297 0.457(44) 0.511(34) 2.27(24) 2.26(20)0.331 0.398(44) 0.441(39) 2.01(24) 2.08(17)0.365 0.395(54) 0.438(29) 1.84(34) 1.88(14)0.398 0.368(49) 0.399(35) 1.65(27) 1.75(15)0.431 0.280(68) 0.364(40) 1.47(31) 1.71(19)0.463 0.288(56) 0.361(37) 1.71(27) 1.71(13)0.494 0.266(56) 0.301(49) 1.42(25) 1.47(19)0.556 0.33(11) 0.285(46) 1.13(42) 1.10(26)0.586 0.256(72) 0.301(40) 0.75(31) 1.27(18)0.616 0.239(75) 0.273(47) 0.92(30) 1.12(29)0.646 0.207(84) 0.240(62) 1.45(33) 1.31(22)0.675 0.178(82) 0.239(54) 0.86(34) 0.96(29)0.704 0.180(83) 0.079(73) 0.90(31) 0.77(27)0.732 0.15(11) 0.174(52) 0.32(36) 0.91(21)0.788 0.17(12) 0.188(67) 0.68(40) 0.64(46)0.816 0.09(11) 0.121(64) 0.52(38) 0.16(51)0.843 0.10(11) 0.109(50) -0.05(37) 0.67(26)0.870 -0.00(13) 0.081(56) -0.24(50) 0.17(73)0.923 -0.02(12) 0.02(13) 0.22(44) 0.53(41) G E and G M for ensemble E250. ppendix B: Priors The two-state method, as we apply it, requires a certain amount of prior information inorder to stabilize the ﬁts. In this appendix, we summarize the priors, which are extractedfrom the nucleon two-point functions, for the ratio of overlaps ρ and for the energy diﬀerencebetween ground and ﬁrst excited state; for details see Sec. III. The energy gap is given inlattice units. n D200 C101 H105 N200 N203 N302 J303 S201 S4000 1.35(8) 1.17(7) 1.16(13) 1.37(3) 1.08(5) 1.47(5) 1.74(9) 1.23(4) 1.08(4)1 1.33(7) 1.14(5) 1.24(13) 1.42(3) 1.07(5) 1.55(5) 1.89(11) 1.30(4) 1.17(4)2 1.40(8) 1.13(5) 1.31(11) 1.48(3) 1.04(6) 1.52(5) 2.14(18) 1.40(4) 1.22(4)3 1.39(9) 1.16(6) 1.43(14) 1.54(4) 1.03(7) 1.51(5) 2.28(24) 1.41(5) 1.37(6)4 1.68(13) 1.19(7) 1.26(14) 1.54(4) 1.48(14) 1.83(20) 2.09(20) 1.57(7) 1.40(7)5 1.59(12) 1.20(9) 1.39(18) 1.62(5) 1.37(15) 1.50(9) 2.12(22) 1.69(8) 1.44(8)6 1.66(15) 1.19(7) 1.83(50) 1.66(6) 1.33(17) 1.42(9) 2.52(42) 1.51(9) 1.66(14)8 1.76(23) 1.32(18) 1.69(51) 1.76(8) 1.62(4) 1.40(10) 1.99(17) 1.74(7) 1.70(13)9 2.41(45) 1.29(13) 1.55(35) 1.82(7) 1.64(5) 1.69(7) 2.23(24) 1.89(10) 2.30(21)10 2.25(34) 1.31(10) 1.65(36) 1.77(7) 1.65(5) 1.78(9) 2.50(33) 1.93(11) 2.20(29)11 2.18(51) 1.35(11) 1.48(27) 1.94(9) 1.74(6) 1.86(10) 2.40(33) 2.14(17) 3.04(44)12 2.07(34) 1.42(26) 1.57(21) 2.00(14) 1.72(8) 1.81(13) 2.04(19) 3.28(1.28) 3.50(91)TABLE II. Overlap factors ρ ( p ) for all momenta on ensembles D200, C101, H105, N200, N203,N302, J303, S201 and S400. n E250 n E250 n E250 n E250 n E250 n E2500 1.24(9) 6 1.47(12) 13 1.57(6) 20 1.65(8) 27 1.75(12) 35 1.74(17)1 1.22(10) 8 1.49(14) 14 1.54(6) 21 1.59(8) 29 1.65(12) 36 1.73(15)2 1.28(9) 9 1.50(18) 16 1.70(8) 22 1.63(12) 30 1.76(15)3 1.48(11) 10 1.61(5) 17 1.58(7) 24 1.68(11) 32 1.74(14)4 1.28(11) 11 1.58(5) 18 1.64(7) 25 1.74(11) 33 1.73(14)5 1.35(11) 12 1.66(6) 19 1.57(8) 26 1.71(10) 34 1.72(16)TABLE III. Overlap factors ρ ( p ) for all momenta on ensemble E250. D200 C101 H105 N200 N203 N302 J303 S201 S4000 0.233(18) 0.172(15) 0.253(44) 0.396(12) 0.178(12) 0.238(11) 0.241(19) 0.431(22) 0.445(23)1 0.230(17) 0.135(8) 0.214(35) 0.409(13) 0.199(14) 0.204(9) 0.216(17) 0.415(19) 0.404(19)2 0.202(17) 0.212(19) 0.213(20) 0.423(14) 0.232(18) 0.221(11) 0.200(16) 0.401(17) 0.403(21)3 0.194(19) 0.160(13) 0.270(27) 0.445(16) 0.237(21) 0.264(14) 0.195(16) 0.416(21) 0.354(22)4 0.151(13) 0.241(24) 0.295(31) 0.435(19) 0.181(16) 0.189(12) 0.215(23) 0.408(22) 0.375(34)5 0.163(17) 0.158(14) 0.244(32) 0.434(20) 0.185(20) 0.251(21) 0.213(25) 0.401(21) 0.372(34)6 0.188(21) 0.161(13) 0.243(12) 0.456(23) 0.167(17) 0.326(32) 0.190(21) 0.550(57) 0.334(34)8 0.177(16) 0.476(55) 0.260(20) 0.476(31) 0.228(8) 0.454(54) 0.273(26) 0.597(48) 0.375(37)9 0.167(9) 0.299(47) 0.267(25) 0.493(20) 0.269(10) 0.618(44) 0.255(25) 0.694(27) 0.364(31)10 0.189(13) 0.370(38) 0.299(30) 0.519(23) 0.312(14) 0.634(47) 0.260(24) 0.709(30) 0.548(58)11 0.173(11) 0.442(26) 0.497(80) 0.494(24) 0.481(22) 0.629(54) 0.268(30) 0.665(31) 0.388(32)12 0.193(11) 0.298(76) 0.596(79) 0.495(34) 0.445(15) 0.688(39) 0.357(40) 0.622(59) 0.557(34)TABLE IV. Energy gap ∆( p ) for ensembles D200, C101, H105, N200, N203, N302, J303, S201and S400 in lattice units. n E250 n E250 n E250 n E250 n E250 n E2500 0.223(28) 6 0.219(26) 13 0.377(24) 20 0.369(28) 27 0.372(38) 35 0.311(51)1 0.246(29) 8 0.222(32) 14 0.290(22) 21 0.378(32) 29 0.315(40) 36 0.414(59)2 0.233(27) 9 0.197(32) 16 0.377(27) 22 0.203(14) 30 0.225(36)3 0.204(22) 10 0.380(20) 17 0.364(28) 24 0.208(23) 32 0.355(46)4 0.259(32) 11 0.368(20) 18 0.383(26) 25 0.422(36) 33 0.422(52)5 0.251(31) 12 0.202(14) 19 0.376(31) 26 0.417(38) 34 0.416(57)TABLE V. Energy gap ∆( p ) for ensemble E250 in lattice units. ppendix C: Dipole and z -expansion results This appendix presents the results for the magnetic moment and the electromagnetic radiifor the dipole and z -expansion ﬁts for all ensembles, both for the summation and two-statedata given in appendix A, for a momentum cut of Q ≤ . . The pion mass indicatedin the ﬁrst column of the tables uniquely identiﬁes the corresponding gauge ensemble viaTable I. M π [GeV] µ (sum) µ (two-state) (cid:104) r (cid:105) (sum) (cid:104) r (cid:105) (two-state) (cid:104) r (cid:105) (sum) (cid:104) r (cid:105) (two-state)0.130 3.61(40) 3.87(25) 0.680(67) 0.630(42) 0.528(99) 0.556(64)0.203 4.22(31) 4.12(15) 0.614(43) 0.644(24) 0.630(78) 0.624(47)0.223 4.21(17) 4.19(25) 0.628(29) 0.691(30) 0.561(43) 0.600(84)0.262 3.80(35) 3.63(19) 0.597(52) 0.541(24) 0.299(63) 0.374(31)0.278 4.25(52) 4.09(81) 0.629(59) 0.800(105) 0.511(103) 0.584(175)0.283 4.21(32) 3.93(11) 0.529(39) 0.501(13) 0.446(65) 0.423(23)0.293 - - 0.571(115) 0.482(20) - -0.347 3.96(13) 4.12(12) 0.599(20) 0.606(25) 0.425(30) 0.464(34)0.350 4.19(38) 4.10(20) 0.549(37) 0.544(22) 0.385(68) 0.406(34)0.353 4.47(56) 3.52(23) 0.614(39) 0.515(20) 0.557(112) 0.364(53)TABLE VI. Dipole ﬁts for µ , (cid:104) r (cid:105) and (cid:104) r (cid:105) (in fm ) on every ensemble. M π [GeV] µ (sum) µ (two-state) (cid:104) r (cid:105) (sum) (cid:104) r (cid:105) (two-state) (cid:104) r (cid:105) (sum) (cid:104) r (cid:105) (two-state)0.130 3.72(58) 4.10(34) 0.782(176) 0.814(106) 0.822(629) 1.01(42)0.203 4.35(32) 4.12(17) 0.610(71) 0.652(42) 0.759(151) 0.692(104)0.223 4.29(18) 4.21(24) 0.645(39) 0.662(43) 0.573(92) 0.626(116)0.262 4.06(41) 3.69(22) 0.599(59) 0.521(34) 0.444(119) 0.432(89)0.278 3.70(56) 3.63(71) 0.575(66) 0.622(82) 0.313(195) 0.395(232)0.283 4.20(32) 3.89(12) 0.518(46) 0.479(19) 0.457(101) 0.397(50)0.293 3.37(1.13) 3.58(49) 0.639(104) 0.472(33) 0.390(259) 0.367(155)0.347 3.85(13) 3.94(12) 0.558(20) 0.469(26) 0.392(40) 0.384(49)0.350 3.80(39) 3.76(21) 0.474(37) 0.469(25) 0.263(111) 0.313(58)0.353 3.93(42) 3.34(25) 0.487(39) 0.463(26) 0.403(106) 0.300(90)TABLE VII. z -expansion ﬁts for µ , (cid:104) r (cid:105) and (cid:104) r (cid:105) (in fm ) on every ensemble. ppendix D: HBChPT ﬁts Here we summarize the results for the physical values of the magnetic moment µ = κ + 1and the electromagnetic radii of the HBChPT ﬁts to the z -expansion data, applying apion mass cut M π ≤ .

28 GeV, as explained in Section V. Note that the value for AIC isnot corrected for the cut data points. The last column indicates which of the correctionsappearing in Eq. (23) is included in the ﬁt χ / DOF κ (cid:104) r (cid:105) [fm ] (cid:104) r (cid:105) [fm ] p-value AIC Q -cut [GeV ] correction0.45 3.08(40) 0.856(42) 0.88(12) 0.89 11.62 0.6 -0.52 2.62(62) 0.82(10) 0.94(28) 0.76 16.60 0.6 O ( a )0.53 3.56(75) 0.92(15) 0.86(38) 0.76 16.63 0.6 O ( e − m π L )0.61 3.27(35) 0.871(37) 0.777(91) 0.82 14.71 0.7 -0.80 3.03(55) 0.858(94) 0.82(22) 0.61 20.36 0.7 O ( a )0.45 3.84(50) 0.89(10) 0.91(26) 0.89 17.63 0.7 O ( e − m π L )0.68 3.23(35) 0.869(37) 0.772(91) 0.76 15.44 0.8 -0.88 2.97(54) 0.856(93) 0.80(21) 0.53 21.05 0.8 O ( a )0.54 3.78(49) 0.89(10) 0.91(26) 0.83 18.32 0.8 O ( e − m π L )0.65 3.21(34) 0.872(36) 0.774(90) 0.78 15.18 0.9 -0.85 2.97(54) 0.851(92) 0.82(21) 0.56 20.76 0.9 O ( a )0.53 3.75(48) 0.87(10) 0.93(26) 0.83 18.28 0.9 O ( e − m π L )TABLE VIII. HBChPT ﬁts for z -expansion extractions from summation data. χ / DOF κ (cid:104) r (cid:105) [fm ] (cid:104) r (cid:105) [fm ] p-value AIC Q -cut [GeV ] correction1.00 3.25(26) 0.875(36) 0.772(88) 0.43 16.03 0.6 -0.99 2.77(49) 0.781(84) 0.74(21) 0.42 18.97 0.6 O ( a )0.66 3.88(74) 1.24(21) 0.83(45) 0.66 17.29 0.6 O ( e − m π L )0.79 3.26(24) 0.877(30) 0.748(74) 0.65 16.65 0.7 -0.70 2.88(42) 0.804(68) 0.72(16) 0.69 19.59 0.7 O ( a )0.74 4.03(58) 0.92(12) 0.94(35) 0.66 19.90 0.7 O ( e − m π L )0.95 3.34(23) 0.875(28) 0.732(68) 0.49 18.41 0.8 -0.63 2.86(40) 0.792(62) 0.68(15) 0.75 19.08 0.8 O ( a )0.90 4.13(55) 0.92(12) 0.96(32) 0.52 21.16 0.8 O ( e − m π L )0.89 3.31(22) 0.875(27) 0.735(66) 0.55 17.76 0.9 -0.59 2.87(38) 0.794(61) 0.69(15) 0.79 18.68 0.9 O ( a )0.90 4.01(53) 0.90(11) 0.91(30) 0.52 21.18 0.9 O ( e − m π L )TABLE IX. HBChPT ﬁts for z -expansion extractions from two-state data. ppendix E: Covariant B χ PT ﬁts

Here we summarize the results for the physical values of the magnetic moment µ = κ + 1and the electromagnetic radii of the direct covariant ChPT ﬁts as discussed in Section V,applying a pion mass cut of M π ≤ . χ / DOF κ (cid:104) r (cid:105) [fm ] (cid:104) r (cid:105) [fm ] (cid:104) r E (cid:105) [fm] (cid:104) r M (cid:105) [fm] p-value AIC Q -cut correction1.63 3.75(11) 0.818(22) 0.663(27) 0.905(12) 0.814(17) 0.00 112.11 0.6 -1.68 3.83(19) 0.816(33) 0.652(35) 0.903(18) 0.808(21) 0.00 115.87 0.6 O ( a )1.68 3.82(12) 0.790(31) 0.666(42) 0.889(18) 0.816(26) 0.00 116.45 0.6 O ( e − M π L )1.67 3.77(32) 0.800(53) 0.662(63) 0.894(29) 0.814(39) 0.00 115.84 0.6 ∗ O ( a )1.69 3.79(12) 0.770(45) 0.671(36) 0.878(26) 0.819(22) 0.00 116.65 0.6 ∗ O ( e − M π L )1.49 3.81(11) 0.818(22) 0.669(28) 0.905(12) 0.818(17) 0.01 91.44 0.5 -1.55 3.80(20) 0.819(34) 0.670(36) 0.905(19) 0.818(22) 0.01 95.44 0.5 O ( a )1.56 3.87(12) 0.797(33) 0.667(45) 0.893(19) 0.817(28) 0.01 96.36 0.5 O ( e − M π L )1.54 3.82(32) 0.798(53) 0.668(64) 0.893(30) 0.817(39) 0.01 95.12 0.5 ∗ O ( a )1.56 3.85(12) 0.783(47) 0.670(37) 0.885(27) 0.819(23) 0.01 96.43 0.5 ∗ O ( e − M π L )1.42 3.74(13) 0.800(26) 0.651(35) 0.895(15) 0.807(22) 0.04 64.64 0.4 -1.48 3.62(22) 0.803(40) 0.668(45) 0.896(22) 0.817(28) 0.03 68.20 0.4 O ( a )1.52 3.80(13) 0.793(36) 0.664(58) 0.891(20) 0.815(35) 0.02 69.60 0.4 O ( e − M π L )1.48 3.68(33) 0.777(59) 0.653(69) 0.881(34) 0.808(43) 0.03 68.16 0.4 ∗ O ( a )1.52 3.79(13) 0.776(51) 0.639(55) 0.881(29) 0.799(34) 0.02 69.61 0.4 ∗ O ( e − M π L )1.75 3.76(13) 0.797(30) 0.660(42) 0.893(17) 0.812(26) 0.01 56.95 0.3 -1.88 3.76(25) 0.805(48) 0.661(53) 0.897(27) 0.813(32) 0.00 60.87 0.3 O ( a )1.91 3.81(14) 0.781(45) 0.655(79) 0.883(25) 0.810(49) 0.00 61.55 0.3 O ( e − M π L )1.87 3.80(35) 0.776(67) 0.658(77) 0.881(38) 0.811(48) 0.00 60.71 0.3 ∗ O ( a )1.90 3.82(14) 0.769(60) 0.631(72) 0.877(34) 0.794(46) 0.00 61.41 0.3 ∗ O ( e − M π L )TABLE X. Covariant BChPT ﬁts for summation data. / DOF κ (cid:104) r (cid:105) [fm ] (cid:104) r (cid:105) [fm ] (cid:104) r E (cid:105) [fm] (cid:104) r M (cid:105) [fm] p-value AIC Q -cut correction1.35 3.499(77) 0.791(18) 0.653(21) 0.889(10) 0.808(13) 0.03 94.66 0.6 -1.35 3.68(13) 0.779(26) 0.633(25) 0.883(15) 0.795(16) 0.04 95.44 0.6 O ( a )1.42 3.509(90) 0.806(34) 0.614(47) 0.898(19) 0.783(30) 0.02 100.06 0.6 O ( e − M π L )1.34 3.74(18) 0.755(36) 0.647(38) 0.869(21) 0.804(24) 0.04 95.32 0.6 ∗ O ( a )1.37 3.530(86) 0.726(55) 0.554(45) 0.852(33) 0.744(30) 0.03 96.70 0.6 ∗ O ( e − M π L )1.19 3.554(81) 0.787(20) 0.662(24) 0.887(11) 0.814(15) 0.16 74.45 0.5 -1.21 3.66(13) 0.776(27) 0.649(28) 0.881(16) 0.806(18) 0.14 77.29 0.5 O ( a )1.23 3.599(95) 0.826(38) 0.647(52) 0.909(21) 0.804(33) 0.12 78.55 0.5 O ( e − M π L )1.18 3.76(19) 0.750(38) 0.659(42) 0.866(22) 0.812(26) 0.17 75.67 0.5 ∗ O ( a )1.24 3.617(84) 0.821(78) 0.615(47) 0.906(43) 0.784(30) 0.11 79.02 0.5 ∗ O ( e − M π L )1.46 3.529(91) 0.779(24) 0.661(29) 0.882(14) 0.813(18) 0.03 66.33 0.4 -1.51 3.61(15) 0.753(39) 0.650(34) 0.868(22) 0.806(21) 0.02 69.20 0.4 O ( a )1.34 3.58(10) 0.908(53) 0.645(66) 0.953(28) 0.803(41) 0.08 63.06 0.4 O ( e − M π L )1.42 3.73(20) 0.705(50) 0.656(47) 0.839(30) 0.810(29) 0.04 66.10 0.4 ∗ O ( a )1.48 3.566(98) 0.844(94) 0.569(63) 0.919(51) 0.754(42) 0.03 68.08 0.4 ∗ O ( e − M π L )1.77 3.45(10) 0.772(29) 0.642(39) 0.879(16) 0.801(24) 0.01 57.65 0.3 -1.87 3.54(17) 0.740(50) 0.633(45) 0.860(29) 0.796(28) 0.00 60.71 0.3 O ( a )1.65 3.58(12) 0.959(69) 0.727(96) 0.979(35) 0.852(56) 0.02 54.87 0.3 O ( e − M π L )1.77 3.64(21) 0.690(60) 0.639(56) 0.830(36) 0.800(35) 0.01 58.05 0.3 ∗ O ( a )1.88 3.50(11) 0.86(11) 0.584(85) 0.926(58) 0.764(55) 0.00 60.97 0.3 ∗ O ( e − M π L )TABLE XI. Covariant BChPT ﬁts for two-state data. ppendix F: Ratio G M /G E A common observation among lattice determinations of the isovector form factors is thatthe ratio of the magnetic and electric form factor exhibits a rather ﬂat behavior. One maytherefore hope to extract the magnetic moment as the intercept of a linear ﬁt to that ratioover a restricted range in Q . We perform linear ﬁts with an upper limit of Q ≤ .

29 GeV for E250 and Q ≤ . for the remaining ensembles. Note that S201, S400 and N302do not allow for a linear extrapolation with a cut oﬀ Q ≤ . . We stress that thesepoints do not enter in our ﬁnal ﬁts, but rather serve as a consistency check (see Fig. 9 andFig. 10). Ensemble µ (summation) µ (two-state)D200 4.16(28) 3.73(18)C101 4.36(17) 4.55(70)H105 3.60(61) 3.61(73)N200 4.11(28) 3.91(7)N203 4.04(12) 4.02(28)J303 4.11(38) 3.78(16)E250 4.05(52) 4.06(39)TABLE XII. Extrapolated values of the magnetic moment using the ratio G M ( Q ) G E ( Q ) ..