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

We present results for the isovector (p - n) electromagnetic form factors of the nucleon using 11 ensembles of gauge configurations generated by the MILC Collaboration using the highly improved staggered quark action with 2 + 1 + 1 dynamical flavors. These ensembles span four lattice spacings a approximate to 0.06, 0.09, 0.12 and 0.15 fm and three values of the light-quark masses corresponding to the pion masses M-pi approximate to 135, 225 and 315 MeV. High-statistics estimates using the truncated-solver method allow us to quantify various systematic uncertainties and perform a simultaneous extrapolation in the lattice spacing, lattice volume and light-quark masses. We analyze the Q(2) dependence of the form factors calculated over the range 0.05 less than or similar to Q(2) less than or similar to 1.4 GeV2 using both the model-independent z expansion and the dipole Ansatz. Our final estimates, using the z-expansion fit, for the isovector root-mean-square radius of nucleon 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 estimates from the dipole Ansatz, r(E) = 0.765(11)(8) fm, r(M) = 0.704(21)(29) fm and mu(p-n) = 3.975(84)(125) Bohr magneton, are consistent with those from the z expansion but with smaller errors. Our analysis highlights three points. First, all our data for form factors from the 11 ensembles and existing lattice data on, or close to, physical mass ensembles from other collaborations collapse 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 ways of analyzing the data 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 for < r(E)(2)>, < r(M)(2)> and mu, high-statistics data for Q(2) < 0.1 GeV2 are needed.