Nucleon electromagnetic form factors in the continuum limit from ( 2+1+1 )-flavor lattice QCD

Results are presented for the nucleon isovector electromagnetic form factors using 11 ensembles generated by the MILC collaboration using the 2+1+1-flavors HISQ action. They span 4 lattice spacings $a \sim$ 0.06, 0.09, 0.12 and 0.15~fm and 3 values of $M_\pi \sim 135, 225$ and 315 MeV. High-statistics estimates are used to perform a simultaneous extrapolation in the lattice spacing, lattice volume and light-quark masses. The $Q^2$ dependence over the range 0.05-1.4 ${\rm GeV}^2$ is investigated using both the $z$-expansion and the dipole form. Final $z$-expansion estimates for the isovector r.m.s. radius are $r_E = 0.769(27)(30)$ fm $r_M = 0.671(48)(76)$ fm and $\mu^{p-n} = 3.939(86)(138)$ Bohr magneton. The first error is the combined uncertainty from the leading-order analysis, and the second is an estimate of the additional uncertainty due to using the leading order chiral-continuum-finite-volume fits. The dipole estimates, $r_E = 0.765(11)(8)$ fm, $r_M = 0.704(21)(29)$ fm and $\mu^{p-n} = 3.975(84)(125)$, are consistent with those from the $z$-expansion but with smaller errors. Our analysis highlights three points. First, all data from the eleven ensembles and existing lattice data on, or close to, physical mass ensembles from other collaborations collapses more clearly onto a single curve when plotted versus $Q^2/M_N^2$ as compared to $Q^2$ with the scale set by quantities other than $M_N$. The difference between these two analyses is indicative of discretization errors, some of which presumably cancel when the data are plotted versus $Q^2/M_N^2$. Second, the size of the remaining deviation of this common curve from the Kelly curve is small and can be accounted for by statistical and possible systematic uncertainties. Third, to improve lattice estimates, high statistics data for $Q^2<0.1$ ${\rm GeV}^2$ are needed.


I. INTRODUCTION
Experiments studying electron scattering off protons and neutrons have a long history of providing an understanding of the structure of nucleons [1,2]. Quantitative understanding of the distribution of charge is described by the electric and magnetic form factors G E ðQ 2 Þ and G M ðQ 2 Þ, respectively [3]. Quantities of phenomenological interest obtained from the slope of the form factors at fourmomentum transfer squared Q 2 ¼ 0 are the electric and magnetic charge radii of the nucleons. At present there is a 6σ discrepancy between the electric charge radius of the proton obtained from electronic energy levels combined with electron scattering data [4] versus that from the Lamb shift in muonic hydrogen E μp ð2S − 2PÞ [5,6]. A second issue that needs resolution is the behavior of the ratio G E =G M at Q 2 > 1 GeV 2 [7] and whether this ratio crosses zero at about 8 GeV 2 as indicated by experiments at JLab [8,9]. In this work, we focus on determining the electromagnetic form factors in the range 0.05 ≲ Q 2 ≲ 1 GeV 2 and extracting the charge radii from them.
The electric and magnetic form factors, G E and G M , of the nucleon can be calculated directly from large-scale simulations of lattice QCD. In recent years, advances in algorithms and computing power have allowed the community to push the calculations toward physical masses for the light u and d quarks and on lattice spacings that are small enough that discretization effects are expected to be at the few percent level [10][11][12][13]. In this paper we present results from 13 calculations on 11 ensembles that cover a range of lattice spacings (0.06 ≲ a ≲ 0.15 fm), pion masses (135 ≲ M π ≲ 320 MeV) and lattice volumes (3.3 ≲ M π L ≲ 5.5). These ensembles were generated using 2 þ 1 þ 1 flavors of highly improved staggered quarks (HISQ) [14] by the MILC Collaboration [15]. This suite of calculations allows us to understand and assess various sources of systematic errors. The analysis is carried out using both the dipole Ansatz and the z expansion, which give consistent estimates for the isovector mean-square charge radii hr 2 E i and hr 2 M i and the magnetic moment μ p−n . Our final results for the isovector mean-square charge radii hr 2 E i and hr 2 M i (also for Dirac, hr 2 1 i, and Pauli, hr 2 2 i, radii) defined in Eqs. (10)- (12) and for the magnetic moment μ are given in Table IX. We also present a comparison with lattice data obtained close to the physical pion mass by other collaborations and with the Kelly parameterization of the experimental data [16] in Fig. 22. Our estimates for hr 2 E i, hr 2 M i and μ p−n are about 17%, 19%, and 16% smaller, respectively, than the phenomenological values given in Eq. (D1) and the precise experimental value in Eq. (9). Throughout this paper, we have paid attention to the size of possible statistical and systematic errors and find that a linear combination of these is large enough to explain the deviations.
We analyze the world data for G E and G M in Sec. VII and find that data from all 13 of our calculations and those from other collaborations done at or near the physical pion mass fall roughly onto a single curve when plotted versus Q 2 or Q 2 =M 2 N . However, there is a noticeable shift between the two curves when compared to the Kelly fit. The difference between the two ways of analyzing the data is a discretization artifact: Specifically, it is a consequence of the difference in values of the lattice scale obtained from different observables. The size of the difference again indicates that the present underestimate of hr 2 E i, hr 2 M i and μ should not be considered significant. Our overall conclusion is that to significantly reduce the systematics and improve the precision with which these observables can be extracted will require high-statistics data at smaller values of the lattice spacing and with Q 2 < 0.1 GeV 2 . We stress that the long-term goal of lattice QCD is to directly predict the form factors and not to reproduce the Kelly curve, a parameterization of the experimental data. Throughout this paper, we use the Kelly curve to provide a reference point for comparison and for discussing systematics and trends in the lattice data. We do not show an error band on the Kelly curve as it is negligible on the scale of the errors in the lattice data.
This paper is organized as follows. In Sec. II, we review the theory, computational approach and the status of the experimental and phenomenological results. In Sec. III, we describe the salient features of the calculation. The fits used to isolate excited-state contamination (ESC) and extract the form factors are described in Sec. IV. Fits to quantify the Q 2 behavior of the (p − n) form factors are discussed in Sec. V, and the extraction of our final results for the isovector mean-square charge radii hr 2 E i and hr 2 M i and the anomalous magnetic moment μ p−n are presented in Sec. VI. Comparisons with form factors extracted from experiments and with previous lattice QCD calculations are made in Sec. VII. We end with conclusions in Sec. VIII. Some further details of the calculations are given in four Appendixes: lattice parameters in Appendix A, analysis of nucleon mass in Appendix B, ESC in Appendix C, and a review of the experimental data for the form factors in Appendix D.

II. ELECTROMAGNETIC FORM FACTORS OF THE NUCLEON
The Dirac, F 1 , and Pauli, F 2 , form factors are extracted from the matrix elements of the electromagnetic current within the nucleon state N through the relation hNð ⃗p f ÞjV em μ ð⃗ qÞjNð ⃗p i Þi where ⃗ q ¼ ⃗p f − ⃗p i is the momentum transfer. The discrete lattice momenta are given by 2πn=La with the entries of the vector n ≡ ðn 1 ; n 2 ; n 3 Þ taking on integer values, n i ∈ f0; Lg. The spacing between the momenta is controlled by the spatial lattice size La. The normalization used for the nucleon spinors in Euclidean space is X s u N ð ⃗p; sÞū N ð ⃗p; sÞ ¼ Eð ⃗pÞγ 4 − i⃗ γ · ⃗p þ M 2Eð ⃗pÞ ; ð2Þ In the isospin symmetric limit, the difference of its matrix elements between a proton and a neutron state are related to the isovector form factors of the proton by the relation The quantity we calculate on the lattice is the left-hand side of Eq. (4), i.e., the isovector form factors of the proton. Throughout this paper, the term isovector form factors of the proton and the (p − n) form factors refer to the same quantities as defined in Eq. (4). These will henceforth be analyzed in terms of the spacelike four-momentum squared, Another common set of definitions of the electromagnetic form factors, widely used in the analysis of experimental data, are the Sachs electric, G E , and magnetic, G M , form factors that are related to the Dirac and Pauli form factors as From these, the vector charge is given by and the difference between the magnetic moment of the proton and the neutron by The anomalous magnetic moments of the proton and the neutron, in units of the Bohr magneton, are known very precisely [17]: The electric and magnetic size of the nucleon are defined as the slope of the form factors with respect to Q 2 at Q 2 ¼ 0 [18]: The form factors G E;M are normalized by their values at Q 2 ¼ 0: G E ðQ 2 ¼ 0Þ ≡ g V and G M ðQ 2 ¼ 0Þ=g V ≡ μ. This definition makes them independent of the renormalization constant Z V of the lattice vector current and improves the signal because some of the systematics cancel in the ratios. Therefore, in this work, we will use Eq. (10) when calculating hr 2 E i and hr 2 M i. Note that Z V g V ¼ 1 as the electric charge is conserved. A second independent estimate of Z V , obtained using nonperturbative lattice calculations in the regularization-independent symmetric momentum-subtraction (RI-sMOM) scheme, is given in Ref. [19], where the difference between the two estimates was shown to be ≲3%.
One similarly defines the isovector Dirac and Pauli mean-square radii as These are related to hr 2 E i, hr 2 M i and μ ≡ 1 þ κ as Our analysis of the lattice data is carried out in terms of G E and G M . After extracting hr 2 E i and hr 2 M i in Sec. VI, we also give results for hr 2 1 i and hr 2 2 i in Table IX using these relations.
The electric root-mean-square charge radius r E ≡ ffiffiffiffiffiffiffiffi ffi hr 2 E i p of the proton has been measured in three ways: (i) laser spectroscopy of the Lamb shift in muonic hydrogen [5,6,20], (ii) continuous-wave laser spectroscopy of hydrogen [21], and (iii) elastic scattering of electrons off protons [22,23]. Results using electrons, i.e., the latter two ways, are included in the committee on data (CODATA) 2014 world average [4,24]: and the third result is from muonic hydrogen. The large difference between the CODATA 2014 and muonichydrogen values was termed the "proton radius puzzle." The new CODATA 2018 value [25] resolves the puzzle in favor of the muonic-hydrogen result. The magnetic radius of the proton extracted from experiments using electrons is [4,24] r p M ¼ 0.776ð38Þ fm: Values for the isovector charge radii, extracted from the experimental data and used to compare lattice data against, are given in Eq. (D1) in Appendix D.
To reduce the uncertainty in results from electron scattering experiments, which have been done down to Q 2 ≈ 0.004 GeV 2 , new experiments to constrain the low Q 2 behavior have been initiated [26,27]. Similarly, for lattice QCD calculations to help resolve the puzzle, we need to calculate the form factors to Q 2 ≈ 0.004 GeV 2 to extract r E with better than 1% accuracy.
A challenge to the direct extraction of hr 2 i i from the lattice data is that the value of the smallest momenta, 2π=La, is large in typical lattice simulations. In our calculations, it is ≳220 MeV, and the range of Q 2 values, given in Table I, are in the range 2-10M 2 π . It is, therefore, traditional to fit the data for the G i to an Ansatz and then use the fit to evaluate the derivative given in Eq. (10). Both using an Ansatz and estimating its parameters from fits to data with ffiffiffiffiffiffi Q 2 p ≳ 200 MeV introduce systematic uncertainties when evaluating the derivative at Q 2 ¼ 0. We estimate the dependence of hr 2 i i on the choice of the Ansatz by comparing results for each ensemble obtained using two different fits: the dipole model and the z expansion.
Two alternate approaches are, first, to calculate the form factors at fixed Q 2 and extrapolate these to the continuum limit and then fit the Q 2 behavior. Unfortunately, the values of Q 2 are different on each ensemble. Second, combine the dipole or the z-expansion parameterization of the Q 2 behavior with the chiral-continuum-finite-volume (CCFV) Ansatz for one overall fit. This combined fit is discussed in Sec. VI D. The central analysis presented here consists of first fitting the data versus Q 2 using the dipole model and the z expansion to extract hr 2 i i and μ on each ensemble and then getting the physical results from a CCFV fit in a, M π and M π L that addresses the associated systematics.
It is important to note that both the electron scattering experiments and lattice QCD calculations suffer from paucity of data close to Q 2 ¼ 0 that impacts the extraction of the charge radii. However, there is a large range, 0.004 ≲ Q 2 ≲ 1 GeV 2 , over which accurate experimental data exist. Thus, more than just extracting the charge radii, our goal is to directly compare the lattice and the experimental data over this range of Q 2 as discussed in Sec. V.
An Ansatz that is commonly used to fit the experimental data is the dipole. It arises if one assumes an exponentially falling charge distribution. The resulting form factor is characterized by a single parameter, the mass M, and normalized to F 1 ¼ G E ¼ g V at Q 2 ¼ 0. It goes as Q −4 in the Q 2 → ∞ limit in accord with perturbation theory [28]. The second Ansatz is a model-independent parameterization called the z expansion [29,30]: where the a k are fit parameters and z is defined as with t cut ¼ 4M 2 π denoting the nearest singularity in G E;M ðQ 2 Þ. In terms of z, the domain of analyticity of G E ðQ 2 Þ is mapped into the unit circle with the branch cut at Q 2 ¼ −4M 2 π [30]. We analyzed the data with two choices of the parametert 0 :t 0 ¼0 andt 0 mid ¼ f0.12;0.20;0.40gGeV 2 for the M π ≈f135;220;315gMeV ensembles. By choosing the value oft 0 to lie in the middle of the range of Q 2 at which we have data, we reduce the maximum value of z and hope to improve the stability of the estimates. In practice, for our data, we find that the quality of the fits and the results from different truncations of the series are insensitive to the choice oft 0 . The final results for the charge radii and magnetic moment are obtained from fits usingt 0 mid . The values of Q 2 ¼ ⃗p 2 − ðE − mÞ 2 for the 13 calculations are given in Table I. Note that in four cases there are only five nonzero values. The data for G E ðQ 2 Þ and G M ðQ 2 Þ versus z witht 0 mid are shown in Fig. 8. As discussed in Sec. V B, we restrict our fits to Q 2 ≤ 1 GeV 2 because the reliability of some of the higher Q 2 data is questionable.
To implement the perturbative behavior G i ðQ 2 Þ → Q −4 as Q 2 → ∞ [28] in the z expansion requires Q n G i ðQ 2 Þ → 0 for n ¼ 0, 1, 2, 3. These constraints can be incorporated into the z expansion as four sum rules [31]: For n ¼ 0 it reduces to P k max k¼0 a k ¼ 0. Using these sum rules ensures that the a k are not only bounded but must also decrease at large k [31].
A key issue in the z-expansion analysis is the value of k max required to obtain results with a certain precision. The analysis of the experimental data carried out in Appendix D shows that results stabilize for k max ≈ 4 with and without sum rules. For the lattice data, the choice has to take into account the number of values of Q 2 at which data have been generated to not overparameterize the fit. For our data and fits without priors, the a k fluctuate and the higher-order coefficients (k ≥ 4) are ill determined due to the overparameterization of the fits. To avoid the resulting large fluctuations in a k , we put a bound on them as suggested in [31]. For G E and G M =5, we constrain ja k j ≲ 5.0 for all k by using Gaussian priors with central value zero and width five. With this constraint, results for hr 2 E i, hr 2 M i and μ do not change significantly for k max ≤ 3 and stabilize for k max ≥ 4 as shown in Fig. 9. The convergence of estimates from fits with sum rules is slower and occurs for k max ≥ 7 as also shown in Fig. 9. We, therefore, use the fits without sum rules and with k max ¼ 4 for our final results. Since hr 2 E i, hr 2 M i and μ are best extracted from data at small Q 2 , the sum rule constraints imposed to guarantee the large Q 2 behavior are not essential for their determination. In short, results with sum rules are used only as consistency checks.
Overall, the fits to G E ðQ 2 Þ are more stable than those to G M ðQ 2 Þ. The main reason is the extra data point at G E ðQ 2 ¼ 0Þ which pins down the sign of the slope of G E ðQ 2 Þ at small Q 2 . Using a value for G M ð0Þ, derived from the ratio G M ðQ 2 Þ=G E ðQ 2 Þ as discussed in Sec. IV B, greatly improved the stability of fits to G M ðQ 2 Þ.

III. LATTICE METHODOLOGY
The parameters of the 13 calculations done on 11 HISQ ensembles are the same as used in Ref. [19] for the calculation of isovector charges. To keep this work selfcontained, the lattice parameters of the calculations and the number of measurements made are summarized in Table XII in Appendix A. The parameters used to generate the Wilson-clover quark propagators using the multigrid algorithm [32] are also given in Table XIII. We remind the reader that two ensembles, a06m310 and a06m220, have been analyzed twice with different smearing parameters giving a total of 13 calculations. Also, compared to Refs. [33,34], six ensembles (a12m220S, a12m220, a12m220L, a09m310, a09m220 and a09m130W) have been simulated afresh with randomly chosen source points on each configuration to increase their statistical independence, and data at a larger number of momenta have been accumulated.
To increase the statistics cost-effectively, we used a truncated solver with bias correction method [35,36]. We also used the coherent source method [37,38] to construct sequential propagators from the sink time slice, at which a zero-momentum nucleon state is inserted.
The details of our strategy for the calculations and the analysis have been published in earlier works [19,33,34]. Here we provide a brief summary of the points relevant to the calculation of the electric and magnetic form factors: (i) In our calculation, the nucleon operator used is with color indices fa; b; cg, charge conjugation matrix C ¼ γ 0 γ 2 , and q 1 and q 2 denoting the two different flavors of light Dirac quarks. The quark propagator is smeared both at the source and the sink using a gauge-invariant Gaussian smearing procedure [39] described in Appendix A. To improve the signal, the amplitude A 0 0 with which the nucleon interpolating operator at the source time slice couples to the ground state j0 0 i with energy E 0 and momentum p 0 should be large while the coupling to excited states should be small. We find that for the smearing parameters given in Table XIII, the signal in all ten momentum channels analyzed is good in most cases. Note that the nonrelativistic projection ð1 AE γ 4 Þ=2, inserted to improve the signal [33,40,41], as well as the smearing at the source and the sink time slices, breaks Lorentz covariance. (ii) All errors are determined using a single-elimination jackknife method over configurations; i.e., we first construct the bias corrected average for each configuration and then carry out the fits to the two-and three-point functions within the same jackknife procedure over these configuration averages.
(iii) To control excited-state contamination, we use the same toolkit as in Ref. [19]. The two-point functions are fit keeping four states in the spectral decomposition. The amplitudes and the masses obtained from these fits are input into the analysis of threepoint functions. The results for the masses are given in Table XIV in Appendix B. (iv) On each ensemble, we calculate the three-point functions at multiple values of source-sink separation τ. These values of τ, given in Table XII, are the same as in Ref. [19]. (v) The insertion of the vector current at definite momenta p is carried out on each time slice t between the source and the sink and for each value of τ. These data for the three-point functions, C ð3ptÞ Γ ðt; τ; p 0 ; pÞ, at a large number of values of t and τ are fit using three states in the spectral decomposition: where the source point is translated to t ¼ 0, the operator is inserted at time t, and the nucleon state is annihilated at the sink time slice τ, which numerically is also the source-sink separation. In this relation, the numbers refer to the state jni, and a state with superscript prime denotes that it could have nonzero momentum p 0 . The momentum p at the sink is fixed to zero in all three-point functions. (vi) With our data, the term h2 0 jO Γ j2i could not be resolved. So, in all the fits we set the contribution of the term with h2 0 jO Γ j2i equal to zero and call these 3 Ã -state fits. There is a caveat to the a12m220S ensemble: The p 0 ¼ 0 data are analyzed using 3 Ãstate fits, while the p 0 ≠ 0 data are fit using two states because the 3 Ã -state fits for Q 2 ≠ 0 are unstable.  Table I. These are obtained using the nucleon ground-state energy Eð ⃗pÞ extracted using 4-state fits to the twopoint functions. (viii) To extract the desired matrix element h0 0 jO Γ j0i using Eq. (20), the masses M i , energies E i , and the amplitudes jA i j and jA 0 i j are taken from the fit to the two-point function within one overall jackknife procedure. This procedure assumes that the full tower of excited states and their parameter are determined accurately by the fits to the two-point functions.
(ix) Off-diagonal terms with nonzero momentum transfer such as jA 0 i jjA j jhi 0 jO Γ jji are related to jA 0 j jjA i jhj 0 jO Γ jii by a combination of Lorentz boost, parity and Hermitian transformation provided the tower of states and the coupling to them are the same on either side of the operator. Since the latter is not guaranteed, we treat all such pairs of matrix elements as independent free parameters in the fits. (x) The data for three-point functions at nonzero momentum transfer are not symmetric about the midpoint, τ=2, between the source and the sink. Nevertheless, in the simultaneous 3-state fit to the data with multiple source-sink separations τ and intermediate times t, we skip the same t skip points adjacent to the source and the sink for every τ to remove points with the largest ESC. Two considerations motivated this choice: (i) The time slice of the onset of the plateau in the nucleon effective mass plot is roughly independent of the momentum as shown in Refs. [19,34], and (ii) because we choose the values of t skip to be as small as possible based on the stability of the covariance matrix used in the fits. The values of t skip used here are the same as in Ref. [19]. (xi) The vector current in the continuum theory is conserved; however, the local vector current used in our lattice calculations is not. The renormalization constant Z V for this current has been determined in two ways: (i) nonperturbatively in the RI-sMOM scheme and then converted to MS using perturbation theory and (ii) measured directly from the matrix element of V 4 at Q 2 ¼ 0, i.e., 1=g V . The two sets of values are compared in Ref. [19] and differ by up to 3%. This size of difference is not unreasonable in our clover-on-HISQ formulation which has discretization effects starting at Oðα s aÞ. Our final results are obtained using the ratios G i ðQ 2 Þ=G E ð0Þ, method (ii), in which some of the systematics cancel. The discretization errors in hr 2 E i, hr 2 M i and μ are addressed by the continuum extrapolation, a part of the CCFV fit. The key input, other than statistical precision of the three-point data, that impacts the stability of the nstate fits to control ESC and obtain the ground-state matrix elements is the energy of the first excited state since the terms with h1 0 jO Γ j0i and h0 0 jO Γ j1i give the dominant contribution. Once the ground-state matrix elements have been determined, the procedure for obtaining the form factors from them is described in the next section.

IV. EXTRACTING FORM FACTORS FROM MATRIX ELEMENTS
The following ratios R μ of the three-point to the twopoint correlation functions: give the desired ground-state matrix elements (ME) h0 0 jO Γ j0i, introduced in Eq. (20), in the limits t → ∞ and ðτ − tÞ → ∞. In the calculation of the nucleon threepoint functions, we use the spin projection operator P 3 ¼ ð1 þ γ 4 Þð1 þ iγ 5 γ 3 Þ=2. With this P 3 , and the vector current defined in Eqs. (3) and (4) with Euclidean γ μ , the following quantities have a signal and give either the electric or the magnetic form factors: Note that, in practice, these ratios are used only to plot the data. Our results are obtained by making n-state fits to the correlation functions. Exploiting the cubic symmetry under spatial rotations, we construct two averages over equivalent three-point correlators before doing fits to get the ground-state matrix elements: over ReðC 1 Þ and ReðC 2 Þ for G M ðQ 2 Þ and over ImðC 1 Þ, ImðC 2 Þ and ImðC 3 Þ for G E ðQ 2 Þ. We label these form factors as G V i M and G V i E . Together with G V 4 E extracted from Eq. (24), they constitute the three form factors analyzed. Their extraction is straightforward as each of the three is given by a distinct three-point function. At the same time, it is important to note that the discretization artifacts and the excited-state contaminations in each can be very different.
The data for the ratio defined in Eq. (21) and the results of 3 Ã fits to the three three-point correlators are illustrated in Figs. 25-31 and 34-35, respectively. The ideal expected behavior of all three-point functions with large t and τ − t is a flat region near τ=2 that becomes independent of τ. Our data show that this is not manifest even at τ ≈ 1.4 fm. We, therefore, use 3 Ã -state fits to data at the various values of t and τ to obtain estimates of the ground-state matrix elements. Results for the three sets of form factors, G V 4 E , G V i E and G V i M , extracted from these matrix elements using Eqs. (22), (23) and (24) are given in Tables II-IV for the 13 calculations.
A. Extraction of G E ðQ 2 Þ The pattern of the ESC in the extraction of G V i E versus G V 4 E can be, and is found to be, very different as shown in Figs. 25 and 26 for the two physical mass ensembles. The data for G V 4 E show a clear monotonic but slow convergence from above and a flattish region near the middle. The estimates of the τ → ∞ values given by the 3 Ã fits are found to be stable under variations in t skip and the values of τ included in the fits.
The data for G V i E show much larger ESC and the ME h0 0 jO Γ j1i and h1 0 jO Γ j0i are an order of magnitude larger for n 2 ¼ 1 as compared to those from G V 4 E . The resulting pattern versus t is essentially linear for each τ. As τ is increased, this "line" rotates toward becoming flat, but the rotation is slow. The pivot point is approximately the point of intersection of the various τ lines and converges to the ground-state estimate as t and ðτ − tÞ → ∞.
The difference in the shape of the ESC between G V i E and G V 4 E can be explained by the behavior of the transition matrix elements under parity transformation and Hermitian conjugation. The imaginary parts of the matrix elements of V i at nonzero momentum pick up a negative sign under the combined transformations. As a result, for example, the term jA 0 0 jjA 1 jh0 0 jO Γ j1ie −E 0 t−M 1 ðτ−tÞ has opposite sign to that of its partner jA 0 Thus, each such pair of terms give a "sinh"-like correction, that makes the data look like a straight line at an angle to the extracted ground-state result. On the other hand, the matrix elements in the related pairs of terms from the real parts of V i and V 4 have the same sign and therefore exhibit a "cosh"-like correction. Even in this case, the magnitudes of the two ME in such pairs of terms are not the same. Therefore, in fits to the three-point data using Eq. (20), we leave all the matrix elements as free parameters. In fact, in practice, it is the product of the amplitudes and the ME, such as jA 0 0 jjA 1 jh0 0 jO Γ j1i, that are free parameters in the fits. Thus, in effect, only the energies are free parameters and these are taken from the two-point functions.   Table XII. The results are obtained using 4-state fits to the two-point functions and 3 Ã -state fits to the three-point functions (2-state fits for the a12m220S ensemble) as described in the text. The value G E ð0Þ ¼ 1=Z V given in the first row provides one estimate of the renormalization constant for the vector current. The momentum transfer Q 2 , in units of GeV 2 , associated with each ⃗ n is given in Table I (5) 0.290(7) ⃗ n a09m130W a 06m310 a06m310W a 06m220 a06m220W a 06m135 (0,0,0) 1.052 (6) 1.043 (6) 1.035 (11) 1.050 (7) 1.039 (9) 1.042 (10) (12) 0.380 (20) It is also evident from Figs. 25 and 26 that the ESC in G V i E is the largest at the smallest nonzero momentum; i.e., the "angle" the data make with the horizontal line is the largest. On the other hand, the ESC in G V 4 E increases with momentum. By comparing the data in the two figures, we also conclude that the ESC increases with decreasing a for both G V i E and G V 4 E . A consequence of this difference in the ESC behavior is that the errors in G V i E are 3-10 times larger than in G V 4 E (see data in Tables II and III). Also, since one cannot extract a value for G E ðQ 2 ¼ 0Þ using the operators V i due to kinematic constraint, the fits to G V i E versus Q 2 , discussed in Sec. V B, are less stable because they are not anchored at Q 2 ¼ 0. As a result, the extraction of the electric charge radius from the G V i E data has much larger errors. Because of these two reasons, it has been common to analyze only G V 4 E ðQ 2 Þ. With our high-statistics data, we are able to compare the ESC, the efficacy of the 3 Ã fits, and the discretization errors between G V 4 E and G V i E . A comparison of results for G V i E and G V 4 E is presented in Fig. 1 for the 13 calculations. As stated above, the errors in G V i E are much larger than those in G V 4 E ; however, there are two additional noteworthy patterns. First, the data for G V i E for Q 2 ≲ 0.2 GeV 2 on the a12m220, a09m220, a06m220 and the two physical mass ensembles a09m130W and a06m135 have the largest errors and mostly lie below those from G V 4 E . On the other hand, the data for Q 2 ≳ 0.2 GeV 2 overlap in most cases. Our conclusion, based on these data, is that for Q 2 ≳ 0.2 GeV 2 the two measurements can be considered to have the same mean but with different variance.
The pattern of data at Q 2 ≲ 0.2 GeV 2 is puzzling and we do not have an explanation for the larger errors or the systematic differences. In particular, we cannot discern whether they are due to residual ESC, statistical fluctuations and/or different discretization errors. In summary, while our high-statistics data have allowed us to quantify the larger errors and fluctuations in G V i E , we do not have a resolution for the difference. Operationally, using a weighted average of the nonzero Q 2 data from G V i E and G V 4 E , i.e., assuming that the differences are statistical fluctuations, gives results that are essentially identical to those from G V 4 E . We, therefore, analyze only the data from G V 4 E in the rest of this paper. To establish full control over all systematics, future calculations should demonstrate consistency between G V i E and G V 4 E .  momentum transfer, the convergence is monotonic from below as shown in Fig. 29 for n 2 ¼ 1. The ESC is observed to grow with decreasing a and M π . The pattern of convergence changes with Q 2 : For small n 2 it is from below but by about n 2 ¼ 6, it has changed to from above in most cases as illustrated in Fig. 30. As a result, removing ESC increases the value of G M ðQ 2 Þ at small momentum transfers and decreases it at larger momenta. Consequently, if ESC is not removed, both the magnetic charge radius and the magnetic moment extracted are underestimated.

B. Extraction of
The results of the 3 Ã fits to the data for the bare form factor G M ðQ 2 Þ are summarized in Table IV. A key shortcoming of the analysis of the lattice G M ðQ 2 Þ is the lack of data at Q 2 ¼ 0. To overcome this, we note that the ratio G M =G E , shown in Fig. 12, is, within errors, linear in Q 2 for Q 2 ≲ 0.6 GeV 2 . We, therefore make a linear fit to the ratio of the form factor data, G M ðQ 2 Þ=G E ðQ 2 Þ, with momenta up to ⃗ n ¼ ð2; 1; 1Þ to obtain an estimate for the renormalized G M ðQ 2 ¼ 0Þ. The corresponding unrenormalized values, which we call derived G M ðQ 2 ¼ 0Þ, are also given in Table IV. These values are indistinguishable from those obtained from taking a ratio of the two correlators and then making a linear fit versus Q 2 to these data. Including these values of G M ðQ 2 ¼ 0Þ improved the stability of the z-expansion fits. Note that the extrapolation of G M =G E , inclusion of the extrapolated value of G M ð0Þ, and the fit to G M are done within a single jackknife loop; therefore, the statistical errors are accounted for correctly.
To estimate the importance of using the derived point G M ðQ 2 ¼ 0Þ, which anchors the fits to data, especially on ensembles with largish values of the minimum Q 2 , we performed the following test. We fit the nonzero Q 2 data for G V 4 E to extract the value and the slope at Q 2 ¼ 0 for each ensemble. Comparing the value for g V from this fit with the data given in Table II, we find the magnitude of the The first row gives data for the M π ≈ 310 MeV ensembles; the second row for the M π ≈ 220 MeV ensembles; the third for the two physical mass ensembles a09m130W and a06m135; and the data for the remaining three calculations are shown in the fourth row. The solid black line shows the Kelly fit to the experimental isovector, G p−n E , data.
difference for the dipole and z 4 fit is between 0.01 and 0.04 for the 13 calculations. The difference in the slope, hr 2 E i, compared to the data in Table V is up to 9% for the dipole fit and up to 20% for the z 4 fit. Based on this test, it is not unreasonable that an uncertainty of similar size can be present in the extraction of μ and hr 2 M i. Thus, to get highprecision results without resorting to a derived value for G M ð0Þ or without using priors requires having data at smaller values of Q 2 .
C. Dependence of G E ðQ 2 Þ and G M ðQ 2 Þ on the lattice parameters In Figs. 2-6, we explore the dependence of the renormalized form factors G V 4 E ðQ 2 Þ=g V and G V i M ðQ 2 Þ=g V , which we henceforth label G E ðQ 2 Þ=g V and G M ðQ 2 Þ=g V for brevity, as a function of the pion mass, lattice spacing, lattice volume and the smearing size. The significant features are as follows.
(i) The dependence of G E ðQ 2 Þ=g V on the pion mass, keeping the lattice spacing roughly constant, is shown in Fig. 2. The data show a steeper fall off as the quark mass is lowered. The behavior of G M ðQ 2 Þ=g V is similar as shown in Fig. 4. (ii) The data for G E ðQ 2 Þ=g V do not show any significant dependence on the lattice spacing a for fixed pion mass as shown in Fig. 3. A similar insensitivity to change in a is exhibited by G M ðQ 2 Þ=g V as shown in Fig. 5. Estimates for hr 2 M i from z-expansion fits without including our derived value for G M ð0Þ are, in many cases, unstable even for the z 3 or z 3þ4 fits; i.e., estimates for r M become negative. We conclude that the fits in these cases are overparameterized. Including the derived value of G M ð0Þ and imposing the constraint on a k discussed in Sec. II greatly improved the z-expansion fits. On the other hand, the dipole fits give consistent estimates with or without using a value for G M ð0Þ. Our final results for both types of Q 2 fits are obtained including the G M ð0Þ points.
Lastly, the comparison of the lattice data with the Kelly fit to the experimental data is shown in Figs. 3 and 5. Both G E ðQ 2 Þ=g V and G M ðQ 2 Þ=g V move toward the Kelly curve as M π and a are reduced. However, G E ðQ 2 Þ=g V from the two physical mass ensembles still shows significant deviations from the Kelly fit. The data for G M ðQ 2 Þ=g V show a different curvature from the Kelly curve and points with Q 2 ≲ 0.2 GeV from the physical mass ensembles move below the Kelly curve. This change in behavior in G M ðQ 2 Þ=g V results in an underestimate of both hr 2 M i and the magnetic moment μ p−n as discussed in Sec. VI.

Dependence on lattice size
Simulations on large lattices are not only important for reducing finite-volume effects but also provide the simplest solution to obtaining data at smaller Q 2 for fixed a and M π . To demonstrate the improvement possible, we compare data from the a12m220S, a12m220 and a12m220L ensembles in Fig. 31 in Appendix C. As the data move to smaller Q 2 with increasing L, the statistical quality of the signal also improves for a fixed number of measurements.
In Fig. 32, we show G E and G M versus Q 2 for these three ensembles. The data on the two larger volumes, a12m220 In Fig. 33, we compare the results of three fits to G E ðQ 2 Þ and G M ðQ 2 Þ given in Tables II and IV versus Q 2 for these three ensembles. For the z-expansion fits, the results for hr 2 E i and hr 2 M i from the two larger volumes are consistent within 1σ, while those on a12m220S differ. We find no significant difference in the dipole fits. These comparisons indicate that finite-volume corrections are smaller than the statistical errors on the two larger volumes corresponding to M π L ≳ 4.4. For this reason, we carry out CCFV fits including (11-point fit) and discarding the a12m220S point (10 Ã -point fit). Operationally, the fits are insensitive to the 12m220S point due to the larger errors in it. Nevertheless, our final results, presented in Sec. VI, are from the 11-point fit.
The bottom line is that increasing L for fixed M π and a improves the analysis in a number of ways because the values of Q 2 for a given ⃗ n decrease. First, the statistical errors for a fixed number of measurements decrease. The reduction in errors roughly compensates for the increase in cost of each measurement due to a larger volume. Second, with the decrease in Q 2 , the ESC in G V 4 E becomes smaller, while that in G M becomes easier to control using n-state fits. Lastly, the extraction of hr 2 E i, hr 2 M i and μ improves, since the fit parameters are determined from data with values of Q 2 closer to zero.

Dependence on smearing size
In Figs. 34 and 35 in Appendix C, we compare the ESC in G V 4 E and G V i M for two different smearing sizes using data from the a06m310 and a06m220 ensembles. The data show that the ESC is smaller with the larger smearing size.
The results of the dipole, z 4 and z 5þ4 fits to G E ðQ 2 Þ and G M ðQ 2 Þ versus Q 2 for these two ensembles are shown in Fig. 36. Results for hr 2 E i and hr 2 M i are consistent within 1σ for the two smearings. The data in Fig. 6, however, show that estimates of μ can differ by about 5% between the two calculations with different smearing size. This level of difference can be explained by a combination of statistical and systematic uncertainties in either calculation.

Dependence on lattice scale setting
The two places the lattice scale enters our calculation is in converting Q 2 a 2 to physical units and in the CCFV fits. In Table XII, we give the values of a for the HISQ ensembles obtained by the MILC Collaboration using the Sommer scale r 1 [15,42]. In Table XIV in Appendix B, we give the value of M N obtained on each ensemble using these values of a and fit them using the leading-order CCFV fit defined in Eq. (B1). The result in the continuum limit is M N ¼ 976ð20Þ MeV. The deviation of about 4% from the experimental value indicates a systematic uncertainty of 2%-6% in the scale obtained from r 1 versus M N , the latter analyzed using the leadingorder CCFV fit. The question then is, how does this difference impact the analysis of the form factors and the extraction of hr 2 E i, hr 2 M i and μ? The lattice data plotted in Figs. [2][3][4][5] show that the dependence of the form factors on M 2 π and a is small. To explore the dependence further, we remove the use of a taken from the analysis of the Sommer scale r 1 on the HISQ ensembles by plotting the data versus Q 2 =M 2 N in Fig. 7  (bottom) where the lattice values of M N are used to construct the dimensionless ratio Q 2 =M 2 N for the lattice data and M N ¼ 939 MeV for the Kelly curve. The relative movement between the data and the Kelly curve, when plotted versus Q 2 =M 2 N as compared to Q 2 , brings the data closer together onto a single curve as can be seen by comparing the top and bottom set of panels. For the physical mass ensembles, the size of the relative movement of data depends only on the discretization errors, i.e., the value of M N at that value of a, assuming finite-volume corrections are negligible. Presuming a cancellation of some of the systematics when the data are plotted versus Q 2 =M 2 N , this comparison indicates that the observed larger deviation from the Kelly curve, when the data are plotted versus Q 2 , can be explained partly as a systematic effect due to discretization errors, i.e., variations in the lattice scale set using different observables. This systematic is avoided if data at a given Q 2 are first extrapolated to M π ¼ 135 MeV and then to a ¼ 0 before comparing with the Kelly curve. An attempt at doing this is described in Sec. VI D.
It is important to note that hr 2 E i, hr 2 M i and μ extracted for each ensemble are unchanged whether one calculates them using Q 2 or Q 2 =M 2 N as the independent variable in Eq. (10). The result would be different if the product M 2 N hr 2 E i is calculated on each ensemble and extrapolated to the continuum limit first and the result divided by the experimental value for M N . We discuss this analysis in Sec. VI C.
Having made clear that part of the noticeable spread in the behavior of the form factors shown in Fig. 7 can be accounted for, in a large part, as due to discretization errors, the question is-what is the more robust way of analyzing the data? Should we use the scale set using r 1 or work with dimensionless variables in units of M N ? While our analysis has exposed this systematic, our conclusion is that a larger dataset, or the use of a lattice action with much smaller discretization errors or a better-determined extrapolation Ansatz, are needed to significantly reduce such systematics. Having highlighted the size of this systematic uncertainty, most of the analysis presented below is carried out versus Q 2 . We provide comparison with results plotted versus Q 2 =M 2 N at appropriate places and analyze data for M 2 N hr 2 E i and M 2 N hr 2 M i in Sec. VI C.

V. CHARACTERIZING THE Q 2 BEHAVIOR OF THE FORM FACTORS
In order to extract the charge radii defined in Eq. (10) and the magnetic moment in Eq. (8), we need to parameterize the form factors versus Q 2 . The two fits we explore are the dipole and the z expansion truncated at some power k as discussed in Sec. II. Since the dipole Ansatz is the solution to an exponentially falling charge distribution (thus a model) and the z expansion involves a truncation plus a constraint on the size of the coefficients a k , it behooves us to first test these Ansätze on the high-precision experimental data as discussed next.

A. Experimental data for the form factors and their Q 2 behavior
Electromagnetic form factors of nucleons are extracted from differential cross sections measured in the scattering of electrons off nuclei. The process of going from measurements of the differential cross sections to nucleon form factors is nontrivial and involves modeling [1,31,43]. As already stated, we have two reasons to analyze the experimental data: to compare them against the lattice data over the range 0 < Q 2 ≲ 0.8 GeV 2 and to test the efficacy of the dipole and z-expansion fit Ansätze. For these purposes, we have collected together compiled experimental data for the proton and the neutron in Appendix D (see Figs. 37 and 38). From these, we have determined the Kelly parameterization for the isovector combinations G p E − G n E and G p M − G n M . Henceforth, for brevity, we will continue to use G E ðQ 2 Þ=g V and G M ðQ 2 Þ= g V to represent the (p − n) combinations when comparing the lattice and the experimental data.
Next, we test the fit Ansätze on the experimental data. The results for hr 2 E i and hr 2 M i for the proton are shown in Fig. 37. Based on the χ 2 =DOF, the dipole fit works surprisingly well for G E ðQ 2 Þ, and the deviation from the data is less than a percent over the range 0 < Q 2 ≤ 1 GeV 2 . This difference is much smaller than the precision of our lattice data. For G M ðQ 2 Þ, the deviation is larger (up to 6%) and the χ 2 =DOF of the fit is poor. In the z-expansion fits with constraints, results for hr 2 E i, hr 2 M i and μ stabilize for k ≥ 5 as shown in Fig. 39.
Based on this analysis, and as noted in Appendix D, one should not expect a match between our lattice and the experimental data to better than about 5% or be able to resolve differences between the dipole and the z-expansion fits at or below this level. These comparisons provide a framework for our lattice analyses using the z expansion: extract hr 2 E i and hr 2 M i from k ¼ 4 to avoid overparameterization for some of the ensembles.
In the final estimates, we have assigned an additional systematic uncertainty to account for the fact that the CCFV fits have been made using just the leading-order corrections. This is discussed further in Sec. VI.

B. Analysis of the lattice QCD data for the form factors
A comparison of the form factors G E ðQ 2 Þ=g V and G M ðQ 2 Þ=g V from all 13 simulations with the Kelly parameterization of the experimental (p − n) data is shown in Fig. 7. The data for G E ðQ 2 Þ=g V lie above the Kelly curve with those from the two physical mass ensembles being the closest as shown in Fig. 3, whereas the data for G M ðQ 2 Þ=g V lie about the Kelly curve for Q 2 ≳ 0.2 GeV 2 and then fall below it for smaller Q 2 as highlighted in Fig. 5. In both cases, these deviations from the Kelly curve impact the slope at Q 2 ¼ 0; i.e., both hr 2 E i and hr 2 M i come out smaller than the phenomenological estimates. More importantly, the very precisely measured magnetic moment G M ð0Þ ¼ μ p − μ n is underestimated by about 16%. As remarked above in Secs. IVA and IV B, removing the ESC using the 3 Ã fits increases the value of all three; nevertheless, the final results presented in Sec. VI are smaller than the experimental values. Furthermore, deviations of the lattice form factors from the Kelly curve are apparent over a range of Q 2 .
As discussed in Sec. IV C 3, part of the difference between the Kelly curve and the data is due to the mismatch FIG. 8. The data for G E and G M =ðμ ¼ 4.7058Þ plotted versus z for the 13 calculations. The vertical red line on the left corresponds to Q 2 ¼ 0, while on the right to Q 2 ¼ 1 GeV 2 except for a15m310 (Q 2 ¼ 1.4 GeV 2 ) and a12m220 (Q 2 ¼ 0.8 GeV 2 ). In the two physical mass cases, the right vertical red line lies outside the panel. in the scale set by r 1 and M N . This is highlighted in Fig. 7, where data are plotted versus Q 2 (top panels), evaluated using the lattice scale set by r 1 , and versus the dimensionless variable Q 2 =M 2 N (bottom panels). Note that this change of variable does not impact the results for hr 2 E i, hr 2 M i and μ p−n on each individual ensemble and thus their extrapolated values using the same CCFV Ansatz.
The data for G E and G M =ðμ ¼ 4.7058Þ versus z are shown in Fig. 8. Our overall strategy for extracting hr 2 E i, hr 2 M i and μ p−n is the following: We first determine by eye the largest value of Q 2 up to which the data are smooth in z. Next, since we are interested in the value and slope of the fits at Q 2 ¼ 0, we restricted the data to Q 2 ≤ 1 GeV 2 , except for the a15m310 (Q 2 ≤ 1.4 GeV 2 ) and a12m220 (Q 2 ≤ 0.8 GeV 2 ) ensembles. The allowed range 0-Q 2 j max , where Q 2 j max is the largest value allowed by the cuts defined above, is marked by the two vertical red lines in Fig. 8. With these cuts, the points at all Q 2 are retained for most of the ensembles. Only the high Q 2 data for a15m310, a12m310, a12m220S and a12m220 ensembles are removed. These show a break in the smooth behavior in z as is clear from Fig. 8. Going back to the ESC analysis, the reliability of these points is questionable since the data have large errors and the ESC fits were poor.
The results from the z-expansion fits are stable for k ≥ 4 as shown in Fig. 9. Results from fits including the sum rules are similar, except that stability is reached only for k ≥ 7. Estimates from fits with and without sum rules are consistent; however, the errors are larger with the sum rules. The values and χ 2 =DOF of the dipole fits have been stable under increase in statistics for all 13 calculations. On the other hand, the results from the z-expansion fits required high statistics to exhibit convergence with the order of the truncation.
The results from seven fit Ansätze are collected together in Tables V-VII  M i in units of fm 2 from the seven fits to the isovector form factor G M ðQ 2 Þ. The derived value for G M ð0Þ is included in the fits as discussed in the text. The rest is the same as in  10. Results of the dipole, z 4 and z 5þ4 fits to the unrenormalized isovector G E ðQ 2 Þ versus Q 2 (GeV 2 ) for ten ensembles. The top two panels show data from the a15m310 and a12m310 ensembles; the second row from a12m220 and a12m220L ensembles; the third row from a09m310 and a09m220; the fourth row from a06m310 and a06m220; and the fifth row from the two physical mass ensembles a09m130 and a06m135. Estimates of the dipole mass M E (GeV) and the charge radius r E (fm) from the three fits are given in the labels. The numbers within the square parentheses are the χ 2 =DOF of the fit. Data points without circles around them are not included in the fits as explained in the text.
FIG . 11. Results of the dipole, z 4 and z 5þ4 fits to the unrenormalized isovector G M ðQ 2 Þ versus Q 2 (GeV 2 ) for ten ensembles. The rest is the same as in Fig. 10.
used only as consistency checks. The dipole, z 4 and z 5þ4 fits and results are shown in Figs. 10 and 11 for ten ensembles.
The data for the ratio μ p−n × G E ðQ 2 Þ=G M ðQ 2 Þ are shown in Fig. 12. Experimental data indicate that this ratio for the proton is estimated to cross zero around Q 2 ¼ 8 GeV 2 [44]. Our data for the isovector combination do show a negative slope over the region Q 2 ≲ 0.6 GeV 2 ; nevertheless, data at larger Q 2 are needed to determine if and where the ratio crosses zero. As discussed in Sec. IV B, we have used this "linear" behavior at small Q 2 to estimate G M ð0Þ from the ratio, including which helped stabilize the fits to G M ðQ 2 Þ. In the next section, we discuss the CCFV fits used to get the physical estimates.
VI. RESULTS FOR hr 2 E i, hr 2 M i AND μ To obtain results for hr 2 E i, hr 2 M i and μ in the limits a → 0, M π → 135 MeV and M π L → ∞, we make a simultaneous (CCFV) fit in these three variable to the data given in Tables V-VII. Given the spread in the lattice parameters of the 11 ensembles analyzed, we include the leading-order correction term in each of the three variables, i.e., fits with four free parameters, c E;M;μ i . The fit Ansatz for the electric mean-square charge radius used is where the mass scale λ is chosen to be M ρ ¼ 775 MeV and the form of the chiral and FV corrections are taken from Refs. [45][46][47]. For the magnetic mean-square charge radius, we use where the leading dependence on M π is taken from Refs. [45,46]. Lastly, the Ansatz used for the magnetic moment is where the forms of the chiral and finite-volume correction terms are taken from Refs. [45,48]. We express all masses in units of GeV and the lattice spacing in femtometers. In all three CCFV fit Ansätze, Eqs. (25)- (27), heavy baryon chiral perturbation theory (χPT) has been used only to determine the form of the leading-order chiral correction. For example, for μ, χPT predicts the slope c μ 3 of the linear dependence on M π as M N g 2 A =ð4πF 2 π Þ with F π ¼ 92.2 MeV [49]; however, we leave c μ 3 a free parameter. For hr 2 E i and hr 2 M i, we do not have data at enough values of M π to test the contribution of the different terms in the χPT prediction [45] as discussed later in this section. To avoid overparameterization of the fit we, therefore, include only the nonanalytical term in Eqs. (25) and (26). Our focus is on obtaining estimates at M π ¼ 135 MeV, and this is achieved by relying on the data from the two physical mass ensembles to anchor the chiral part of the fit.
In Tables V-VII, we also give the results of the CCFV fits for the following four combinations of the 13 data points: (i) 13-point fit.-All 13 calculations are considered to be independent, even though the a06m310 and a06m220 ensembles have been analyzed twice with different smearing sizes. (ii) 11-point fit.-We use the average of the two values for hr 2 E i, hr 2 M i and μ on the a06m310 and a06m220 ensembles as these have been analyzed twice. In this averaging, we assume maximum correlation between data. (iii) 10-point fit.-The coarsest ensemble, a15m310, is removed from the 11 data points defined above. (iv) 10 Ã -point fit.-We remove the smallest volume ensemble, a12m220S, from the 11 data points defined above. For each of these fits, we give results with (labeled extrap c X 4 ≠ 0) and without (labeled extrap c X 4 ¼ 0) the finitevolume correction. The values of the coefficients are given in Table VIII. In the limit a → 0 and M π L → ∞, only the terms proportional to c X 1 and c X 3 contribute. From these fits, we observe the following: (i) Of our estimate hr 2 E i ≈ 0.59 fm 2 , roughly half comes from c E 1 and the other half from the chiral correction  (iii) There is a significant dependence of μ on the lattice spacing a. As a result, we get a low value, μ ≈ 4 Bohr magneton, in the continuum limit. (iv) The coefficient of the finite-volume term is poorly determined, which is reflected in the larger error estimates with c X 4 ≠ 0. In all cases, the two types of results overlap. To be conservative, we quote all final results including the finite-volume term. In Figs. 13-15, we show the CCFV fits versus a, M π and M π L for three analyses: dipole, z 4 and z 5þ4 . In addition, we  (27) for the 11-point fit used to obtain hr 2 E i, hr 2 M i and μ in the continuum limit from the dipole and z 4 data.
0.32 (18) 0.024(6) −0.14ð19Þ  Table V. In each panel, the CCFV fit (pink band) is shown versus a single variable with the other two variables set to their values at the physical point. The extrapolated values of hr 2 E i are shown using the symbol red star. Fits in a single variable (a or M π ) are shown as gray bands and the corresponding extrapolated value by a black star. The solid red line is the prediction of χPT using the expressions given in Ref. [45].
show fits versus a single variable a or M 2 π (gray bands). When the pink and gray bands are close or overlap, it means that the dominant sensitivity of the CCFV fit is with respect to the single variable of the gray band.
For hr 2 E i, we also show the fit using the χPT expression given in Ref. [45] as a solid red line in Fig. 13. The variation with M 2 π in the data obtained with the dipole, z 4 and the z 5þ4 fits is small over the range 135 < M π < 350 MeV and consistent with the prediction of χPT as discussed in Sec. VI A and shown in Fig. 16: Over the range 135 < M π < 350 MeV, the decrease in the log part is partially compensated for by the increase in the analytical contribution. Only below M π ∼ 135 MeV does the singular term start to dominate.
The shape of the CCFV fit bands are similar for the dipole, z 4 and the z 5þ4 , except that the z-expansion data and the fits have larger errors. A visual overview of all 13 individual results for hr 2 E i and of the four CCFV fits is presented in Fig. 17. The variation with a and M π in the 13 individual calculations is small and somewhat smaller in the dipole than in the various z-expansion estimates.
For hr 2 M i, the variation with a, M π and M π L for each of the three cases, the dipole, z 4 and the z 5þ4 , is small as shown in Figs. 14 and 18. Again, as shown in Fig. 16, to resolve the expected 1=M π chiral behavior [see Eq. (26)] requires data below M π ∼ 135 MeV. Also, the CCFV fits to the data from the three Q 2 fits, shown in Fig. 14, are similar.
The largest variation in μ is versus a as shown in Fig. 15. Again, the fit bands are similar for the dipole, z 4 and the z 5þ4 data. The positive slope versus a lowers the continuum limit result with respect to the experimental value μj expt ¼ 4.7058. The size of the difference between lattice data and experimental results suggests that discretization and other systematic errors in G M ðQ 2 Þ are underestimated. From the summary of the results presented in Fig. 19, it is clear that the largest uncertainty is in the smallest volume, a12m220S, and the two physical mass, a09m130W and a06m135, points.  Table VI. The rest is the same as in Fig. 13.  FIG. 15. The 11-point CCFV fits for μ p−n to the dipole (top), z 4 (middle) and z 5þ4 (bottom) data given in Table VII. The rest is the same as in Fig. 13.   FIG. 16. The prediction of chiral perturbation theory for the isovector hr 2 E i (left) and μhr 2 M i (right) using the expressions given in Ref. [45]. The contribution of the subterms, "constant" (blue dot-dot-dash line), "log" (black dash line), "analytical" (green dot-dash line) and the "1=M π ", defined in the text, are shown separately. Their sum is shown by the solid red line. The red plus sign marks the physical point M π ¼ 135 MeV. The values of the LEC used in these fits are given in the text.

A. Lessons from chiral perturbation theory
It is instructive to compare our data to the predictions of chiral perturbation theory given in Ref. [45]. In Fig. 16 we show the chiral expansion for the isovector hr 2 E i (μhr 2 M i) as a sum of three (four) terms: those independent of M π (labeled constant), proportional to ln M 2 π =M 2 N (labeled log) and terms proportional to powers of M 2 π (labeled analytical). For μhr 2 M i, there is a fourth term proportional to 1=M π . In making these plots, the low-energy constants (LEC) used are For hr 2 E i, the constant term is negative (≈ − 0.49 fm 2 ). In the range M π ¼ 135-350 MeV, hr 2 E i is approximately constant at 0.9 fm 2 ; in this interval, the growth in the ln M 2 π =M 2 N term compensates for the decrease in the analytical term. Below M π ≈ 135, the log term drives the rise in the sum. As shown in Fig. 13, lattice data are significantly smaller in magnitude but show a similar small variation between M π ¼ 135-350 MeV. Because of this small variation, and having data at only three M π values, even including one additional analytical term, say the leading one proportional to M 2 π , in our CCFV fits would overparameterize the fit. Furthermore, over this range, a simple M 2 π term would equally well mimic the sum of the log and the analytical terms. To avoid overparameterization, we have included only one of the possible M π -dependent terms, the log, in our CCFV fits.
In the case of μhr 2 M i [see Fig. 16 (right)], the log and analytical terms are small and the log shows little variation. The largest contributions come from the 1=M π and constant terms. Since the 1=M π term is dominant and has the largest variation with M π , we only include it in the CCFV fit defined in Eq. (26).
In short, even though the χPT-based expressions used for both hr 2 E i and hr 2 M i, given in Eqs. (25) and (26), keep only the leading chiral correction term from the expressions in Ref. [45], the variation in our data at three values of M π spanning 135-315 MeV is small, and including more terms would result in overparameterization. In this situation, having lattice data at M π ≈ 135 MeV is crucial for controlling the uncertainty in the chiral fit to the lattice data.  Table VI from the 13 calculations and the four CCFV fits. In each case we show results for three Q 2 fits: the dipole, z 4 and z 5þ4 . Also, note that the errors we quote in the CCFV fit results are comparable to those in the two physical mass points.
B. Extraction of final results for hr 2 E i, hr 2 M i and μ To choose between z-expansion analysis with and without sum rules, we note that the convergence of estimates given in Tables V-VII, and shown in Fig. 9, is better without sum rules, i.e., k ≥ 4 versus k ≥ 7 versus. Also, hr 2 E i, hr 2 M i and μ should ideally be extracted from the smallest possible range of Q 2 near Q 2 ¼ 0 containing a sufficiently large number of data points. On the other hand, the sum rules are included to get the right asymptotic behavior as Q 2 → ∞. Our final results are therefore obtained as follows: We take the z 4 result without sum rules for the central value and the first error in it represents the analysis uncertainty, i.e., the composite of ESC, Q 2 and CCFV fits. We also quote a second systematic uncertainty to account for having used just the leading-order CCFV fits. This is taken to be the largest of the following: (i) The difference between the two values on the physical mass ensembles, a09m130W and a06m135. The second error estimate for hr 2 M i is given by this difference.
(ii) The difference between the value at a06m135 and the continuum value given by the CCFV fit. This gives the second error estimate for hr 2 E i and for μ p−n , which show the largest variation versus a.
(iii) For the z expansion, we also considered the difference between the z 3 (z 5 ) and z 4 values. These turn out to be smaller than the estimates from the previous two cases. The final results, obtained by applying this prescription to the 11-point CCFV fit values summarized in Tables V-VII,  are given in Table IX. For completeness, we also give, in Table IX, the results for the Dirac and Pauli radii derived using Eq. (12).
The central values for r E p−n , r M p−n and μ p−n are about 17%, 19% and 16% smaller, respectively, than the phenomenological values given in Eq. (D1) and the precise experimental value in Eq. (9). Estimates from the dipole and z-expansion fits, given in Table IX, are consistent; however, the errors in hr 2 E i and hr 2 M i from the z-expansion fits are much larger, and about half the difference from the experimental or phenomenological values. The errors in the dipole fits are small compared to the difference between the lattice and the phenomenological and experimental estimates. As discussed in a number of places above, FIG. 19. A summary of the results for μ presented in Table VII from the 13 calculations and the four CCFV fits. In each case we show results for three Q 2 fits: the dipole, z 4 and z 5þ4 . differences of this magnitude can be accounted for if a linear combination of the statistical and the various systematic errors is taken. The encouraging results from our analysis are (i) the data for both G E ðQ 2 Þ and G M ðQ 2 Þ are seen to converge toward the Kelly parameterization as a and M π are decreased; (ii) the stability of the z-expansion fits improves with statistical precision-however, constraints on the coefficients a k are still needed; (iii) while it is hard to test the nonanalytical chiral behavior predicted in Eqs. (25) and (26) with data at only three values of M 2 π , having data at the two physical mass ensembles anchors the CCFV fit and provides control over the uncertainty in values obtained from the fits.
A weakness of the lattice analysis is that G M ð0Þ cannot be calculated directly due to kinematic constraints. We have motivated the use of a derived value of G M ð0Þ to stabilize fits to G M ðQ 2 Þ. Looking ahead, the most significant improvement needed for extracting all three quantities with higher precision is generating data at smaller values of Q 2 . This, unfortunately, requires ensembles with larger spatial volumes and/or new approaches such as a lattice formulation of the Dirac action with twisted boundary conditions [50,51]. Both options are beyond the scope of this work as they require new simulations.
Two variants of the above analysis are described briefly next.
C. Alternate analysis: Fig. 20 versus M 2 π . For comparison, the phenomenological values for the isovector mean-square charge radii given in Eq.
If some of the systematics cancel in the product, then one may get smaller variation with a and M π . The data for M 2 N hr 2 E i and M 2 N hr 2 M i in Fig. 20 show that there is significant variation with M π . Lacking a wellmotivated fit Ansatz, a reasonable option is to take the average of the values from the two physical mass ensembles. These estimates are again low: M 2 N hr 2 E ij dipole ¼ 13.75ð32Þ, M 2 N hr 2 E ij z 4 ¼ 12.16ð1.22Þ, M 2 N hr 2 M ij dipole ¼ 11.34ð38Þ, and M 2 N hr 2 M ij z 4 ¼ 9.40ð1.50Þ. Similar values are obtained on multiplying the results given in Table IX by FIG. 20. The data for the dimensionless quantities M 2 N hr 2 E i and M 2 N hr 2 M i are plotted versus M 2 π . The top two panels show the data obtained using the dipole fit, and the lower two using the z 4 fit.
Given that the errors are also similar and because these estimates neglect possible a dependence, we do not find this variant of the analysis as providing an obvious improvement.

D. Combined Q 2 -CCFV fit
We also carried out a combined Q 2 -CCFV fit to the G E ðQ 2 Þ and G M ðQ 2 Þ data from the 13 calculations using a product of the z expansion for the Q 2 behavior and the functional forms for the CCFV Ansatz given in Eqs. (25) and (26): Here η represents the vector of variables in the CCFV fit, and each coefficient d k of the z expansion has a CCFV expansion of the form given in Eq. (25) or in Eq. (26). For example, for the four-term CCFV Ansatz given in Eq. (25), η ¼ ð1; a; logðM 2 π =λ 2 Þ; logðM 2 π =λ 2 Þe −M π L Þ, and the combined fit has 20 parameters for the z 4 analysis. In performing these fits, we used Gaussian priors with mean 0 and width 5, in their appropriate units, for all the parameters. The resulting central values of the parameters were within this range.
The central values of the results with and without the finite-volume term are consistent; however, the errors with the finite-volume correction term included are about a factor of 2 larger. In Fig. 21, we show, for the z 4 case, the combined fits neglecting the finite-volume correction term. The results of these combined fits are summarized in Table IX and found to be consistent with those obtained by doing the z 4 and CCFV fits separately (labeled z 4 fit).
All the data included in the comparison are shown in the four panels in Fig. 22. From the plot versus Q 2 , we draw the following overall conclusions: FIG. 22. Comparison of the data for the renormalized isovector G E ðQ 2 Þ and G M ðQ 2 Þ from collaborations that have published results at M π ≈ 135 MeV. The lattice parameters of the various calculations are given in Table X. The data are plotted as a function of Q 2 (top row) and Q 2 =M 2 N (bottom row). The solid line is the Kelly fit to the experimental isovector data. Results for r E , r M , μ and the nucleon mass from published calculations at or near the physical pion mass. The quantity used to set the lattice scale is given in the third column, with r 2 0 Fðr 0 Þ and r 1 extracted from the heavy quark potential [42]. ETMC'18 [57] results are derived from a single fit in Q 2 to the combined 2-and (2 þ 1 þ 1)-flavor data, i.e., neglecting the dependence on the number of flavors N f and the difference in the lattice spacing a. The LHPC'17 [54] results are from a single ensemble and taken from their analysis using the summation method to control ESC. The calculation of the scale used in LHPC'17 is given in Ref. [58] and that by the PACS Collaboration in Ref. [59]. (i) The G E ðQ 2 Þ data approach the Kelly curve from above, while G M ðQ 2 Þ from below for Q 2 <0.2GeV 2 . (ii) No significant dependence on the number of flavors or the lattice spacing a is manifest. (iii) The PACS'18A data at Q 2 < 0.1 GeV 2 , obtained using a large volume, show a qualitatively different behavior and lie closer to the Kelly curve. In this range of Q 2 , the data are almost linear and highly correlated. They give a larger slope in both G M ðQ 2 Þ and G E ðQ 2 Þ and thus larger hr 2 E i and hr 2 M i. (iv) We can also compare data at Q 2 ≈ 0.05 and 0.1 GeV 2 from our a09m130W and a06m135 ensembles, from ETMC'18 [57], and the low-error LHPC'17 [54] points with that from PACS'18A. Given the size of the statistical and systematic errors in individual data points, it is not clear if the observed small differences are significant at these two Q 2 values. Two points become clear on plotting the data versus Q 2 =M 2 N , as shown in the bottom two panels of Fig. 22. First, the collapse of all the data into a single curve over the whole range Q 2 =M 2 N ≲ 0.8 GeV 2 becomes even more pronounced. Second, the deviation of this common curve from the Kelly curve is smaller. Thus, not only do all our data from the 13 calculations fall on a common curve when plotted versus Q 2 =M 2 N , as shown in Fig. 7, but so do data from four other collaborations using different lattice actions and volumes. A priori, a common curve would suggest that all the systematics cancel, and the apparent differences between the various calculations when the data are plotted versus Q 2 was largely a consequence of how the lattice scale is set.
Deviations from the Kelly curve are, however, significant in G E ðQ 2 Þ for Q 2 > 0.1 GeV 2 . Data for G M ðQ 2 Þ data undershoot for Q 2 < 0.2 GeV 2 and are consistent with the Kelly curve above it. These differences are a 2 − 3σ effect and comparable to the size of the shift when the data are plotted versus Q 2 =M 2 N or Q 2 . While a full understanding of how the different systematics contribute and whether there is one that dominates requires future more detailed calculations, we remind the reader that, during the course of our analyses, we have pointed out systematics, for example due to ESC and the deteriorating signal in both the two-and three-point correlation functions at large ⃗ q 2 , that could give rise to uncertainties of this size.
We have already shown that the pattern of ESC in our data is sensitive to the value of Q 2 . In particular, as discussed in Sec. IV, the ESC in correlators from which G E is extracted increases with momentum and the convergence is from above. On the other hand, it is large at small ⃗ q 2 in correlators from which we get G M and the convergence is from below. Thus, possible residual ESC could account for the observed deviation from the Kelly curve.
For G E , there is a clear benefit to performing calculations at small Q 2 . As illustrated in Figs. 25 and 26 for the physical mass ensembles, the ESC in G E ðQ 2 Þ is still small for ⃗ n 2 ¼ 2 corresponding to Q 2 ≈ 0.1 GeV 2 . [Note that for ⃗ q 2 ¼ 0, the ESC is essentially zero as the vector charge is conserved and the local current has no OðaÞ correction in forward matrix elements.] There is, however, an increase in the ESC with decreasing a as shown in the bottom panels in Fig. 27. On the other hand, the ESC in G M ðQ 2 Þ is large at small Q 2 as shown in Fig. 29, and the resulting larger errors in G M ðQ 2 Þ reflect that uncertainty. In contrast, the PACS'18A calculation claims that the ESC is removed in both form factors by using a simple but tuned smearing of sources (exponentially falling along the axis in Coulomb gauge), for generating quark propagators as opposed to the excited-state pattern observed using a gauge-invariant Gaussian smearing in our and the other four calculations summarized in Table X. Clearly, the efficacy of the exponential source used by PACS'18A to remove essentially all ESC needs to be validated.
The collapse of the data into a single curve indicates that finite-volume corrections are already small for M π L ≥ 4, and the main advantage of the large volume used in the PACS'18A [56] study is access to data at low Q 2 . These data for Q 2 < 0.1 GeV 2 represent a qualitative change in the behavior of both G M ðQ 2 Þ and G E ðQ 2 Þ which leads to larger values for hr 2 E i and hr 2 M i. Note that since the PACS'18A estimate M N ¼ 942ð11Þ MeV is close to M phy N ¼ 939 MeV, their data do not move with respect to the Kelly curve when plotted versus Q 2 =M 2 N or Q 2 . The authors attribute the much smaller errors, compared to the much higher statistics PACS'18 [55] calculation, to the use of the all-mode-averaging method and to a better tuned smearing Ansatz (exponential in Coulomb gauge) for the quark sources used to calculate the quark propagators. Since the advantage of simulations on large volume lattices to get data at low Q 2 , and thus reliable estimates for hr 2 E i and hr 2 M i, is obvious, it is important to validate the relatively low-statistics PACS'18A calculation.

VIII. CONCLUSIONS
We have presented calculations of the isovector electric and magnetic form factors G p−n E and G p−n M , using 13 calculations on 11 ensembles of 2þ1þ1 flavors of HISQ [14] fermions generated by the MILC Collaboration [15]. These ensembles are at four lattice spacings, a ≈ 0.06, 0.09, 0.12 and 0.15 fm, and three values of pion masses, M π ≈ 135, 220 and 310 MeV, and the spatial lattice size covers the range 3.3 ≲ M π L ≲ 5.5. Each of these ensembles have been analyzed using Oð10 5 Þ measurements using the truncated-solver method with bias correction. Using these high-statistics data, we demonstrate control over excited-state contamination and perform a simultaneous fit in lattice spacing a, pion mass M π and lattice size M π L to get results at the physical point that can be compared with experimental values.
Our work constitutes three improvements: (i) The much higher statistics allowed us to understand and control ESC better by keeping three states in the spectral decomposition of the three-point correlation functions. (ii) Calculations at multiple values of a and M π show that the variations in the data versus these two parameters is small for Q 2 ≳ 0.1 GeV 2 as illustrated in Figs. 2-5 and 7. (iii) We have presented first results with a CCFV fit to control the lattice artifacts due to discretization, chiral and finite-volume effects. The data for hr 2 E i and hr 2 M i and the CCFV fits in Figs. 13 and 14 show the variation versus M π is small and consistent with the predictions of chiral perturbation theory [45] shown in Fig. 16. In the χPT prediction, the nonanalytical term in M π , included in Eqs. (25) and (26), becomes significant only for M π < 135 MeV, whereas over the range 350 > M π > 135 MeV, its growth is compensated for by the decrease in the analytical corrections. With such competing contributions in M π , the data on the two physical pion mass ensembles at a ≈ 0.09 and 0.06 fm play a significant role in controlling the uncertainty of our final results. (iv) The data for μ exhibit significant dependence on a as shown in Fig. 15. The result in the continuum limit obtained from the CCFV fit is smaller than the precise experimental value by ∼16%. Our final results for the mean-square charge radii, hr 2 E i and hr 2 M i (or equivalently the Dirac, hr 2 1 i, and Pauli, hr 2 2 i, radii derived from them), and the magnetic moment μ are given in Table IX. Using the dipole Ansatz and the z expansion to fit the Q 2 dependence give consistent results; however, the errors from the latter approach are about 2-3 times larger. The central values for hr 2 E i, hr 2 M i and μ are about 17%, 19% and 16%, respectively, smaller than the phenomenological values given in Eq. (D1) and the precise experimental value in Eq. (9). The trend in the data for the form factors, however, is toward the experimental values as the Q 2 , lattice spacing and the light-quark mass are decreased.
With higher-precision data, the major improvement observed has been in the z-expansion estimates. Including constraints on the fit parameters, ja k j ≲ 5, the z-expansion fits for different truncations became consistent. Based on an analysis of the experimental data with the same fit Ansätze and evaluation of various systematics in the lattice calculations, the extraction of charge radii and magnetic moment could have Oð10%Þ errors due to the modeling of the Q 2 behavior. Errors of similar size could also be due to statistics and ESC fits. Keeping in mind these estimates of the magnitude of possible systematics, the total uncertainty in estimates given in Table IX, especially for the dipole fit, are likely underestimated. Consequently, we do not consider the current deviations from the experimental values significant.
The magnitude of the systematic associated with what variable is used to set the lattice scale is exposed by plotting the data versus Q 2 and Q 2 =M 2 N . As shown in Fig. 7 (bottom), the collapse of our data from the 13 calculations onto a common curve is better when plotted versus Q 2 =M 2 N . In Fig. 22, we further show that both G E and G M from all lattice calculations done close to the physical pion mass also collapse onto this curve. This collapse of the data is remarkable considering that the number of quark flavors, lattice size and lattice spacing are different in the various calculations. On the other hand, the shift in the data when G E and G M are plotted versus Q 2 =M 2 N as compared to Q 2 is a discretization effect; i.e., the scale obtained from M N is different from that obtained by the various collaborations using the quantities shown in Table XI. This suggests that there are systematic uncertainties that are at least as large as the shift. The deviation of the combined lattice data from the Kelly curve is significantly reduced when plotted versus Q 2 =M 2 N . Our analysis indicates that the improvement in estimates reported by the PACS'18A [56] Collaboration from a calculation with a large lattice size, M π L ¼ 7.41, is due to data at smaller Q 2 ; i.e., the agreement between data from different collaborations presented in Sec. VII indicates that finite-volume corrections are already small for M π L ≈ 4. The PACS'18A data show no movement with respect to the Kelly curve because their central value for the nucleon mass is consistent with the physical value. This may be because the lattice scale is set using the Omega baryon mass M Ω , which is likely correlated with M N , rather than indicating that discretization errors are already small at a ≈ 0.09 fm. Given the low statistics, we consider it important to validate, in future calculations, their data and confirm that fits to data between 0.01 < Q 2 < 0.1 GeV 2 give estimates of hr 2 E i and hr 2 M i that are consistent with the experimental values.
To conclude, our analysis highlights three points. First, all lattice data are remarkably consistent and the form factors show little dependence on the number of flavors, lattice spacing, quark mass or the lattice volume, at least for data with M π ≲ 300 MeV and M π L ≳ 4. Second, the size of the remaining deviations in G E and G M between the lattice data and the Kelly curve are consistent with the various quantifiable systematics such as excited-state contamination and deteriorating signal at large ⃗ q. Third, current results provide confidence that there are no hidden systematics that afflict the calculations of form factors on the lattice.
With the lattice methodology in place, improved estimates for form factors will be obtained in future highstatistics calculations that provide data at Q 2 < 0.1 GeV 2 and use nucleon interpolating operators that have smaller excited-state contamination. In this Appendix, we summarize, in Table XII, the parameters of the 11 ensembles used in the calculation. Two ensembles, a06m310 and a06m220, have been analyzed twice with different smearing sizes as listed in Table XIII, where we give the parameters used in the generation of the clover propagators. These two tables are essentially the same as in Ref. [19] and have been reproduced here to keep the discussion self-contained. XII. Parameters, including the Goldstone pion mass M sea π , of the 11 (2 þ 1 þ 1)-flavor HISQ ensembles generated by the MILC Collaboration and analyzed in this study are quoted from Ref. [15]. All fits are made versus M val π and finite-size effects are analyzed in terms of M val π L. Estimates of M val π , the clover-on-HISQ pion mass, are the same as given in Ref. [40] and the error is governed mainly by the uncertainty in the lattice scale. In the last four columns, we give, for each ensemble, the values of the source-sink separation t sep used in the calculation of the three-point functions, the number of configurations analyzed, and the number of measurements made using the high-precision (HP) and the low-precision (LP) truncation of the inversion of the clover operator. The smearing size used in the calculation of the quark propagator is given in Table XIII. TABLE XIII. The parameters used in the calculation of the clover propagators. The hopping parameter for the light quarks, κ l , in the clover action is given by 2κ l ¼ 1=ðm l þ 4Þ. m l is tuned to achieve M val π ≈ M sea π . The parameters used to construct Gaussian smeared sources [39], fσ; N KG g, are given in the fourth column where N KG is the number of applications of the Klein-Gordon operator and the width of the smearing is controlled by the coefficient σ, both in Chroma convention [60]. The resulting root-mean-square radius of the smearing in lattice units,

APPENDIX B: NUCLEON SPECTRUM
Our strategy for obtaining the nucleon spectrum is to (i) use the largest possible Euclidean time interval between the source and sink allowed by the stability of the covariance matrix and (ii) keep as many states in the spectral decomposition without overparameterizing the fits. Comparison of results from 2-, 3-and 4-state fits and the corresponding Euclidean time interval used in the fits is given in Ref. [19]. Since some of the results have changed slightly due to additional statistics, the masses of the nucleon ground and three excited states pertinent to this study are summarized in Table XIV. As shown in Ref. [19], the ground-state masses are found to be stable under changes in the number of states kept in the spectral decomposition of the two-point function and the Euclidean time interval used in the fits since the data exhibit a reasonable plateau in the effective mass plot for all the ensembles. Nevertheless, we note the possibility that even the ground-state energy can shift if there exist low-mass-excited states that are not captured by our 4-state fits.
On the other hand, with current statistics, the excitedstate energies are sensitive to the details of the fits. The main reason is the small number, 6-10, of points at short times that are available to determine the six excited-state parameters in a 4-state fit before the ground-state dominates the two-point function. This is particularly true of the a15m310, a09m310 and a06m310W ensembles. In general, even the first excited-state mass depends on the fit range and the number of states kept. Within our strategy of using the largest possible time interval and up to four states in the spectral decomposition, the first excited-state mass is also found to be consistent between the 3-and 4-state fits as shown in Ref. [19].
In the fits to the nucleon two-point function, we find a strong correlation between the excited-state energies and the amplitudes. Thus, we often got a number of values of the fit parameters with similar χ 2 =DOF, in fact regions of nearly the same χ 2 =DOF. We, therefore, resorted to the empirical Bayesian procedure discussed in Ref. [19]. Even using priors poses a challenge: what priors to choose in the 3-and 4-state fits, especially when there are near flat directions in the parameter space. On varying the priors, we found that their choice can impact the second and third excited-state parameters. To address such variation, we iterated the values of the priors and used wide widths. Our main goal was to establish a stable first excited-state energy under a reasonably large range of variations. To summarize, empirical Bayesian priors with wide widths were used for the amplitudes and energies of the three excited states in the 4-state fit. With the current data, we find that estimates of the parameters of the third and fourth states are not fully independent of the input priors and expect much higherprecision data are needed to get stable estimates. In each fit, the propagation of errors, taking into account possible correlations between parameters, is included using the jackknife procedure; however, an uncertainty accounting for the spread in the parameters under variations of the priors is not included in the final results.
Overall, the estimates for the excited-state masses extracted from fits to two-point function are larger than values expected based on phenomenological arguments. For example, the first excited-state mass listed in Table XIV, while stable, is larger than the Nπ and Nππ multiparticle states and the N(1440), "Roper": For our physical mass ensembles these should all be between 1.2 and 1.5 GeV. This is of concern because a reliable estimate of the first excited-state energy is the key variable in 2-or 3 Ã -state fits to control ESC in the three-point functions, where we find that the j0i ↔ j1i transition is the dominant artifact.
The two physical mass ensembles give estimates for M N that are about 13 MeV larger than the physical value M phy N ¼ 939 MeV. To investigate the dependence of M N on the lattice spacing, pion mass and lattice size, we have carried out two fits to the ground-state nucleon mass M 0 given in Table XIV: and where the second relation enforces M phys MeV. We keep the first two correction terms, proportional to fa; a 2 g and fM 2 π ; M 3 π g, to quantify possible curvature in the data while not overparameterizing.  (17) The values of a for the HISQ ensembles used to convert the lattice data to GeV are taken from Ref. [15] and given in Table XII. The two fits are shown in Fig. 23 versus M 2 π , and the values of the fit parameters, with and without the a 2 term, are given in Table XV. The constrained CCFV fit, with and without the a 2 term, is plotted versus a in Fig. 24.
With data at just three values of the pion mass and four values of a, most of the fit parameters, c i , are poorly determined as is evident from Table XV. Qualitatively, we note that the χPT predicted value for the coefficient c 4 ¼ 3g 2 A =ð32πF 2 π Þ ¼ −5.716 using g A ¼ 1.276 and F π ¼ 92.2 MeV, whereas both fits give smaller values. Both curvature coefficients c 2 and c 4 have large uncertainty. The coefficient of the finite-volume correction, c 5 , has been estimated in Ref. [49] and its sign is negative, consistent with our fits.
The unconstrained fit gives a large value, M N ¼ 976ð20Þ MeV, at the physical point. This could be due to unresolved systematics in extracting M 0 and/or the need for including higher-order chiral corrections in the CCFV Ansatz. In short, these fits highlight the need for data at more values of a, M π and M π L to constrain the fit parameters and allow us to carry out a defensible CCFV analysis of lattice calculations of the nucleon mass.
In addition to impacting the analysis of the ESC, the nucleon spectrum also has bearing on the shape of the form factors versus Q 2 if the scale a used to convert Q 2 to physical units is different when extracted from different quantities. In this work, we examine the consequence of the mismatch between the lattice scale a set on each HISQ ensemble using r 1 [15] versus that from M N calculated using Wilson-clover fermions as discussed above. The variation when the form factors are plotted versus Q 2 (with a from r 1 ) and versus Q 2 =M 2 N (with a from M N ) is shown in Fig. 7 and discussed in Sec. V B.

APPENDIX C: ESC IN THE EXTRACTION OF THE FORM FACTORS
In this Appendix, we show the data and the 3 Ã -state fits used to control the ESC in the extraction of the electric and magnetic form factors. There are three sets of figures:  (ii) The improvement in the quality of the signal with increase in the lattice size L is shown in Fig. 31 using data from the three ensembles a12m220L, a12m220 and a12m220S. A study of finite-size effects in the form factors using these three ensembles is made in Fig. 32 and in the extraction of hr 2 E i and hr 2 M i using the dipole, z 4 and z 5þ4 fits in Fig. 33. (iii) A comparison of the ESC with two different smearing sizes is shown in Figs. 34 and 35 for the a06m310 and a06m220 ensembles, respectively.

APPENDIX D: SUMMARY OF EXPERIMENTAL FORM FACTORS
In this Appendix, we collect in one place the experimental data for the form factors for the proton and the neutron. In Fig. 37, we show the data for G p E ðQ 2 Þ and G p M ðQ 2 Þ compiled by Higinbotham [43,61,62] from the cross sections provided in the Lee-Arlington-Hill supplemental material [31], who rebinned the original data obtained by the A1 Collaboration using the MAMI beam at Mainz [3,63]. The neutron data, G n E ðQ 2 Þ, are collected from Refs. [64][65][66] and G n M ðQ 2 Þ from Refs. . These are shown in Fig. 38. From these data, we evaluate the isovector form factors G p−n E ðQ 2 Þ and G p−n M ðQ 2 Þ, against which our lattice data are compared.
The construction of the isovector form factors is done as follows: We first fit the four sets of experimental data for Q 2 ≲ 1 GeV 2 using the Kelly parameterization as shown in Figs. 37 and 38. Using the resulting Kelly fits, we then construct the isovector combinations G p E − G n E and G p M − G n M . This parameterization is used throughout the paper to compare the lattice data against. From this procedure we get r p−n E j exp ¼ 0.929ð27Þ; r p−n M j exp ¼ 0.849ð11Þ: whereas, using the parameter values given in the original Kelly fit [16], we get It is clear that our simple analysis gives larger error estimates. Lattice results for the isovector combination of the radii should be compared to the values given in Eq. (D1) since these are also obtained using the same lesssophisticated analysis. The results of the dipole fits to the proton data shown in Fig. 37 give r p E ∼ 0.833 and r p M ∼ 0.795, which are roughly consistent with the careful analysis of electron-experiment results [43] given in Eqs. (13) and (14) and the Kelly fits shown in the bottom row of Fig. 37. Overall, the dipole Ansatz does a good job of fitting the experimental G E ðQ 2 Þ data, and the deviation is less than 1% for Q 2 < 1 GeV 2 . The dipole fit to G M ðQ 2 Þ is less good as shown by the large χ 2 =DOF.
The convergence of the z-expansion fits, including constraints on the a k , versus k is shown in Fig. 39. Estimates with k ≥ 5 are stable for all three quantities. The results from z-expansion fits, also shown in Fig. 37, are marginally larger than those from the dipole and differ by a few percent from those in Eqs. (13) and (14). The errors in the z-expansion estimates are larger, especially with the inclusion of the sum rules. The overall lesson from this exercise of analyzing the experimental data using the methodology used to analyze the lattice data is that an uncertainty of Oð5%Þ could be present in our analysis using either the dipole or the z-expansion fits.  The horizontal gray band is the τ → ∞ value, and the colored lines are the fit result for τ ¼ 12, 14, 16. The range of the y axis is chosen to be the same for the left panels whereas the total interval Δy ¼ 0.3 is kept the same for the right panels.
FIG. 30. Data and the 3 Ã -state fits to the unrenormalized isovector G M extracted from the Re V i channels using Eq. (22). The data for p 2 ¼ 5ð2π=LaÞ 2 and for eight ensembles are plotted versus t − τ=2. The y-axis total interval Δy ¼ 1.4 is the same for all the plots. The rest is the same as in Fig. 25. The pattern of the ESC changes with momentum as can be seen by comparing with Fig. 29.
FIG. 31. Data and the 3 Ã -state fits to the unrenormalized isovector form factors G V 4 E and G V i M for the three ensembles a12m220L, a12m220 and a12m220S. The first two rows show G E for n 2 ¼ 1 and n 2 ¼ 5, while the last two show G M . The plots for a12m220L are the same as in Figs. 27-30. Note that the data for fixed n but different L cannot be compared since Q 2 , and thus the value of the form factor, changes with the lattice size L. For fixed a and M π , the data shift to smaller values of Q 2 for a given n 2 as listed in Table I. Consequently, the quality of the signal improves with L.  Table XIII. Plots in the first and third rows show data with n 2 ¼ 1 and in the second and fourth rows with n 2 ¼ 5.
FIG. 37. The experimental data for the electric (left) and magnetic (right) form factors G p E ðQ 2 Þ and G p M ðQ 2 Þ for the proton are plotted versus Q 2 . These data [43] are a rebinned version of the data from the A1 Collaboration at Mainz [3] provided by Higinbotham [43]. The top row shows the results of the seven fits used by us to analyze the lattice data. The bottom row shows the same data fit with the Kelly parameterization where "Kelly 2004" refers to using the parameters given in Ref. [16].
FIG. 38. The data for the electric (left) and magnetic (right) form factors of the neutron, G n E ðQ 2 Þ and G n M ðQ 2 Þ, plotted versus Q 2 (GeV 2 ). The G n E ðQ 2 Þ data are compiled from Refs. [64][65][66], and the G n M ðQ 2 Þ data from Refs. . Also shown are the fits with the Kelly parameterization.