Parameterization and applications of the low-$Q^2$ nucleon vector form factors

We present the proton and neutron vector form factors in a convenient parametric form that is optimized for momentum transfers $\lesssim$ few GeV$^2$. The form factors are determined from a global fit to electron scattering data and precise charge radius measurements. A new treatment of radiative corrections is applied. This parametric representation of the form factors, uncertainties and correlations provides an efficient means to evaluate many derived observables. We consider two classes of illustrative examples: first, neutrino-nucleon scattering cross sections at GeV energies for neutrino oscillation experiments; second, nucleon structure corrections for atomic spectroscopy. The neutrino-nucleon charged current quasielastic (CCQE) cross section differs by 3-5% compared to commonly-used form factor models when the vector form factors are constrained by recent high-statistics electron-proton scattering data from the A1 collaboration. Nucleon structure parameter determinations include: the magnetic and Zemach radii of the proton and neutron, $[r_M^p, r_M^n] = [0.739(41)(23), 0.776(53)(28)]$ fm and $[r_Z^p, r_Z^n] = [ 1.0227(94)(51), -0.0445(14)(3)]$ fm, with the dominant uncertainty propagating from the experimental data and radiative corrections, and the second error due to the fitting procedure; the Friar radius of nucleons, $[(r^p_F)^3, (r^n_F)^3] = [ 2.246(58)(2), 0.0093(6)(1)]\,{\rm fm}^3$; the electric curvatures, $[\langle r^4 \rangle^p_E, \langle r^4 \rangle^n_E ] = [1.08(28)(5), -0.33(24)(3)]~{\rm fm}^4$; and bounds on the magnetic curvatures, $[ \langle r^4 \rangle^p_M, \langle r^4 \rangle^n_M ] = [ -2.0(1.7)(0.8), -2.3(2.1)(1.1)]~{\rm fm}^4$.


Introduction
A new generation of precision measurements, including accelerator-based neutrino experiments and muonic atom spectroscopy, demands a rigorous assessment of nucleon structure parameters and their uncertainties. The electromagnetic form factors of the proton and neutron are critical inputs to searches for new physical phenomena and to new precise measurements of the elementary particles. As one example, precise neutrino-nucleus interaction cross sections are required in order to access fundamental neutrino properties at long-baseline oscillation experiments [1][2][3]; the electroweak vector form factors of the nucleons are an important input to these cross sections, and are determined by an isospin rotation of the electromagnetic form factors. As another example, muonic atom spectroscopy [4,5] has opened a new window on the determination of fundamental constants, and has revealed surprising discrepancies in comparisons of different approaches to nucleon structure [6]; it is critical to quantify uncertainties of nucleon structure inputs for the muonic atom program and also to incorporate constraints from these new measurements into other processes, such as the above-mentioned neutrino cross sections.
Recently, with Ye and Arrington [7], two of us provided a new global fit of the proton and neutron electromagnetic form factors, encompassing the entire momentum transfer (Q 2 ) range of available elastic electron scattering data. That analysis provides a comprehensive discussion of the available data, and the fit provides a general purpose tool for studying the form factors over a broad range of Q 2 . However, the fit of Ref. [7] is not optimized for relatively low-Q 2 applications, such as neutrino scattering in the GeV energy regime. First, the inclusion of very high-Q 2 data necessitates the introduction of a large number of parameters, many of which are irrelevant to low-Q 2 applications. Second, the presentation of errors in Ref. [7] (an envelope around the curve as a function of Q 2 ) does not allow a systematic propagation of errors into derived observables. Finally, since the focus of Ref. [7] was in summarizing the implications of electron scattering data in isolation, it did not incorporate the low-Q 2 constraint on the proton electric form factor that emerges from muonic atom spectroscopy. While there is not a complete consensus regarding the outcome of the so-called proton radius puzzle [8][9][10], we believe it is important to be able to consistently incorporate this data and study its impact for applications such as neutrino scattering.
In this paper, we utilize the raw data selections and uncertainty evaluations for electron scattering cross sections from Ref. [7] to present a complete and compact parametric representation and covariance matrix for the form factors suitable for GeV-and sub-GeV scale applications. Section 2 begins by describing the salient features of the data analysis and presenting the fit results. Section 3 considers several illustrative applications, beginning with a discussion of form factor radii and curvatures in Section 3.1. Section 3.2 evaluates neutrino-nucleon scattering cross sections. Section 3.3 presents central values and uncertainties for several nucleon structure parameters that are important for muonic atom spectroscopy. Section 4 provides a summary discussion. Appendix A discusses tensions between data sets. Appendix B provides details on the dispersive evaluation of two-photon exchange radiative corrections. Appendix C compares our numerical results for nucleon structure parameters to previous estimates.

Presentation of form factors
In this section, we begin by recalling definitions and conventions, discuss our data selection and fit procedure, and present the fit results.

Definitions
The Dirac and Pauli form factors, F N 1 and F N 2 , respectively, are defined as matrix elements of the electromagnetic current: where q µ = p µ − p µ , Q 2 = −q 2 = −(p − p) 2 and N stands for p (proton) or n (neutron). In the presence of radiative corrections, the on-shell form factorsF i are IR divergent and we define IR finite form factors [11]. 1 We will present results in terms of the Sachs electric and magnetic form factors, which are related to the Dirac-Pauli basis by For some applications, it is convenient to work with the isoscalar and isovector linear combinations,  Table 1: Number of data points from each data set included in fits below the momentum transfer Q 2 max . The total χ 2 and number of degrees of freedom of each fit are also displayed. 4 The proton electric charge radius can be inferred from the measurement of the 2S-2P Lamb shift in muonic hydrogen [4]; we employ the updated value [5] r p E µH = 0.84087(39) fm .
The neutron electric charge radius is determined from neutron scattering length measurements on heavy targets [82,83], which yield We do not include external constraints on r p M or r n M . We will consider two types of fits: first, a fit of separate proton and neutron data to their respective form factors; second, a fit of combined proton and neutron data to isospin-decomposed form factors. For our default proton fit (line 1 of Table 1), we employ the "Mainz" ep-scattering data set in combination with the proton electric charge radius. For our default neutron fit (lines 2 and 3 of Table 1), we consider en-scattering data with Q 2 ≤ 1 GeV 2 in combination with the neutron electric charge radius. For our default isospin-decomposed fit (line 4 of Table 1), we consider all of the above proton and neutron data. Finally, we also consider an isospin-decomposed fit (line 5 of Table 1) that includes all of the above data, as well as neutron data with 1 GeV 2 < Q 2 ≤ 3 GeV 2 , and "World" and "Pol" data with 0 < Q 2 ≤ 3 GeV 2 . The total number of data points for each of these fits is summarized in Table 1. We also show the total χ 2 and number of degrees of freedom for the fit. 5

Radiative corrections to ep scattering
One goal of the current work is a more robust accounting for radiative corrections to unpolarized ep cross sections in the fits. 6 Besides the standard QED corrections on the electron line, there are three types of radiative corrections that must be applied to scattering data in order to extract the IR finite "Born" form factors defined after Eq. (1). They may be classified as "hadronic vertex", "hadronic vacuum polarization" and "two-photon exchange" (TPE). The first, hadronic vertex, type of correction involves soft radiation and the shape of the event distribution as a function of the inelasticity ∆E = E elastic e − E e (where E e is the scattered electron energy and E elastic e is the elastic limit). This correction is calculable from QED in the soft limit ∆E m π , but is numerically enhanced by large logarithms in that limit. In the "Mainz" data, the soft-photon tail was analyzed in detail but neglected higher-order corrections that are larger than stated systematic uncertainties [11]; however, the bulk of these corrections is absorbed by floating 4 The Mn/Mp factor results from a conventional expression of µn in units of the nuclear magneton, e /2Mp. This difference is insignificant compared to other uncertainties in electron scattering fits, cf. footnote 4 of Ref. [7]. 5 The number of degrees of freedom, n dof , is equal to sum of the number of data points from the respective row in the table, minus the number of form factor parameters (for definiteness, we count kmax − 4 parameters for each form factor that are not fixed by sum rules, cf. Table 2). Note that nuisance parameters in the data set (floating normalization and systematic parameters [7]), are subject to χ 2 constraints (i.e., each nuisance "parameter" is accompanied by a corresponding "data point"), and we do not include them in counting n dof . 6 We follow the analysis of [7,84,85] by omitting radiative corrections to the form factor ratios from polarization data, which are expected to be small compared to other uncertainties. normalization parameters. In the "World" data, uncorrelated uncertainties were included in the data set [7] to account for possible model dependence in the treatment of the radiative tail. In all cases, the error budgets from Ref. [7] are assumed to contain any residual error from approximate treatment of this correction; we include here the small discrepancy between the MS convention defined after Eq. (1), and the commonly-used Maximon-Tjon convention [86] for the soft subtraction (see Appendix B of Ref. [11] for related discussion). The second, hadronic vacuum polarization, type of correction was omitted from the Mainz data set [17], and treated nonuniformly in the World data set. As with the hadronic vertex correction, the error budgets from Ref. [7] are assumed to contain any residual error from approximate treatment of this correction (see Section 1 of Ref. [14] for related discussion).
The third, TPE type of correction, remains a significant contributor to the error budget for epscattering. As with the hadronic vertex correction, the soft-photon part of the TPE correction is computable without uncertainty in QED, while the remaining hard-photon part is removed according to Here dσ expt is the experimental cross section after extracting leptonic QED corrections, and the abovementioned hadronic vertex, vacuum polarization, and soft-photon TPE effects. The resultant dσ Born is identified with the tree-level (Mott) cross section computed using Born form factors. At arbitrary Q 2 , we account for differences between the true TPE correction and the previous default model employed in Ref. [7], called there "SIFF Blunden" [87], by writing (for notational simplicity, we henceforth suppress the subscript "TPE" on δ) For data below Q 2 = 1 GeV 2 , we consider the dispersive analysis from Refs. [88][89][90][91][92][93], which determines δ = δ dispersive ; details are provided in Appendix B. We take the discrepancy between the default model and the dispersive analysis as an uncertainty and allow x = 1 ± 1. Above Q 2 = 1 GeV 2 , we consider the phenomenological correction from Ref. [85], δ = δ default + yδ AMT , which is designed to improve agreement between polarization measurements and TPE-corrected unpolarized Rosenbluth measurements at high Q 2 ; the explicit form of δ AMT is provided in Appendix B. We take the discrepancy between the default model and this phenomenological ansatz as an uncertainty and conservatively allow y = 1 ± 2. Since the ansatz involving δ AMT is purely phenomenological, we perform fits with y = 1 ± 2 enforced as an uncorrelated error, as in Ref. [7], whereas x = 1 ± 1 is enforced as a correlated error. We have verified that taking x uncorrelated or y correlated does not significantly alter the results. 7 As a practical summary, our treatment of radiative corrections follows Ref. [7] above Q 2 = 1 GeV 2 , with the additional parameter x to describe TPE corrections below Q 2 = 1 GeV 2 .

Fit parameters and procedure
Having defined our data sets and treatment of radiative corrections, let us determine the relevant parameters for the z-expansion analysis. We use data with Q 2 ≤ 1 GeV 2 for our default fits, and choose t 0 as in Eq. (5) to minimize the maximum size of |z| in this Q 2 range. We enforce sum rules on the coefficients, and choose the number of free parameters in the z expansion, n max = k max − 4 = 4, sufficiently large so that terms of order |z| nmax+1 are small compared to experimental precision. 8 Our results do not 7 The results for the default proton and iso (1 GeV 2 ) fits are x = 1.41(52) and x = 1.43 (52), respectively. For the iso (3 GeV 2 ) fit treating y as a correlated error, we obtain x = 2.17 (41) and y = 1.74 (60). Loosening the Gaussian bounds on x by a factor of 5, the results for the default proton and iso (1 GeV 2 ) yield small changes to x = 1.55(61) and x = 1.58 (60), respectively. If we do the same for x and y for the iso (3 GeV 2 ) fit, we obtain x = 2.40 (45) and y = 1.81 (64). 8 For the isovector threshold tcut = 4m 2 π and choice of t0 = −0.21 GeV 2 , we have |z| 5 < 0.0033 when 0 < Q 2 < 1 GeV 2 , and |z| 5 < 0.042 when 0 < Q 2 < 3 GeV 2 . For the isoscalar threshold tcut = 9m 2 π and choice of t0 = −0.28 GeV 2 , we have |z| 5 < 0.0007 when 0 < Q 2 < 1 GeV 2 , and |z| 5 < 0.019 when 0 < Q 2 < 3 GeV 2 . Table 2: Parameter choices for the z expansion of the form factors in this paper. Throughout the paper, we use the charged pion mass m π = 0.13957 GeV for the evaluation of t cut . The values for t 0 are obtained by rounding t opt 0 (1 GeV 2 ; 4m 2 π ) ≈ −0.21, and t opt 0 (1 GeV 2 ; 9m 2 π ) ≈ −0.28.
change significantly when k max is increased; we illustrate this in Section 3 by recomputing observables using k max → k max + 1. For all form factors, we have enforced Gaussian bounds, |a k | ≤ 5, |b k | ≤ 5 (k = 1, . . . , k max ) on the coefficients (i.e., a term a 2 k /5 2 is included in the χ 2 function). Our results do not change significantly when this bound is increased by a factor of two. In Table 2, we summarize the choices of z-expansion parameters used in our fits. For each choice of data set in Table 1, the fit returns form factors expressed as central values, errors, and correlations for the indicated number of free parameters.

Fit results
For the proton fit (line 1 of Table 1), the form factor coefficients are For the neutron fits (lines 2 and 3 of Table 1), the form factor coefficients are [a n 1 , a n 2 , a n 3 , a For the isospin-decomposed fit (line 4 of Table 1), the isovector form factor coefficients are and the isoscalar form factors are given by Whereas in Ref. [7] we considered the most inclusive data set, here we have chosen the default proton data set to contain the most recent precise measurements and to minimize internal data tensions. For definiteness we have included neutron data up to the same Q 2 max = 1 GeV 2 . We remark that the Mainz data set predicts r p E = 0.879 (18) fm when the µH charge radius constraint is removed [14]; this value is in only mild tension, 2.2 σ, with Eq. (9). 9 The absence of more severe internal data tensions does not guarantee the absence of potentially underestimated systematics; for a fuller discussion we refer to Ref. [14].
Plots in Appendix A compare our G p E , G p M , G n E , G n M , and G V E , G V M form factors against those of our previous global fit in Ref. [7] and to the BBBA2005 parameterization of Ref. [94]. In Supplementary  Table 3: Electric and magnetic radii of proton and neutron using form factor parameters and bounds of  (83) (38) Material, we provide values for the coefficients and covariance matrices suitable for precise evaluation of charge radii and other physical quantities 10 , as well as values for the coefficients from a fit with k max = 9 which we use in applications to estimate the error from z-expansion truncation.

Illustrative applications
Having determined the form factor coefficients, errors and correlations, let us illustrate with some relevant physical examples. We begin in Section 3.1 by evaluating form factor radii and curvatures. Section 3.2 discusses neutrino scattering applications and Section 3.3 considers nucleon structure parameters for atomic spectroscopy.  The nucleon radii, defined in Eqs. (7) and (8), are presented in Table 3, where each line represents the result of the fit using the corresponding data set in Table 1. For each entry in the table, the first number represents the fit with default k max = 8 (as in Table 2), and the second number represents the fit with k max = 9.

Form factor radii and curvatures
As expected, the output proton and neutron electric radii are driven by the precise external constraints on these quantities and the k max dependence is insignificant and not displayed. The determination of r p E using low-Q 2 data released in 2019 by the PRad experiment is discussed in Appendix A.
The proton and neutron magnetic radii are consistent between fits, and represent best values for these quantities obtained from electron scattering data plus external charge radius constraints. The proton magnetic radius from the default fit, r p M = 0.739(41)(23) fm, should be compared to (and does not essentially alter) our previous extraction from 2010 Mainz data, r p M = 0.776(34)(17) fm. 11 The difference results from omitting the muonic hydrogen constraint in Eq. (9), omitting the floating TPE correction in 10 The linear combination of coefficients that defines the charge radius is more precisely determined with the form factor parameters and the covariance matrix from the Supplementary Material than by evaluation using Eq. (13) and neglecting correlations. 11 A naive average with the analogous fit to world data without Mainz data, r p M = 0.914(35) fm, is used to arrive at the Eq. (12), omitting sum rule constraints on the coefficients, and restricting to Q 2 max = 0.5 GeV 2 . 12 The neutron magnetic radius from our default fit, r n M = 0.776(53)(28) fm, represents a new extraction. Our value may be compared to the result 0.89(3) fm from Ref. [96] that performed a z-expansion fit to ep and en scattering data and ππ → NN data, utilizing a data set for G p M from Ref. [85] that did not include 2010 Mainz data. The PDG recommended value r n M = 0.864(9) fm [95] results from a naive average of this result and the result r n M = 0.862(9) fm from Ref. [97] that performed a global fit of spacelike and timelike data to model spectral functions.
The form factor curvatures, r 4 from Eqs. (7) and (8), are presented in Table 4. Only the proton electric curvature is determined to be nonzero with statistical significance. As for the radii, different fit variations are consistent within uncertainties. We provide previous estimates of curvature in Table 9 of Appendix C.
For both radii and curvatures, the iso (1 GeV 2 ) fit yields a modest reduction in uncertainty compared to the separate proton and neutron fits, which can be traced to more data and a higher threshold t cut (hence smaller range of |z| and smaller number of relevant form factor coefficients) in the isoscalar channel. There is further reduction in the errors using the iso (3 GeV 2 ) fit due to the inclusion of more data. As discussed in Appendix A, this additional data introduces a significant and unresolved tension with the Mainz data set; we focus on the p, n and iso (1 GeV 2 ) fits as our default results.

Neutrino-nucleon scattering
The elementary signal process for neutrino oscillation experiments is charged current quasielastic scattering, 13 ν + n → − + p .
Neglecting power corrections to four-fermion theory of order Q 2 /M 2 W (M W is the W ± boson mass), the cross section in the laboratory frame is: 14 where G F is the Fermi constant, V ud is a CKM matrix element, M = (M p + M n )/2 is the average nucleon mass, m is the final-state lepton mass, E ν is the incoming neutrino energy, and the difference in Mandelstam variables can be written as s − u = 4E ν M − Q 2 − m 2 . The three structure-dependent functions A, B, and C are given by PDG recommended value r p M = 0.851(26) fm [95]. Our current fit corresponds to the "Alternate approach" described in Section VI.C.3 of Ref. [14], which yielded r p M = 0.792(49) fm (line 7 of Table XIV). 12 Removing the µH constraint from our default fit shifts the central value by ∼ 0.7 σ: r p M = 0.739(41) fm → 0.768 (42) fm; further removing the floating TPE parameter, the result would be 0.774(41) fm. 13 For a classic review see Ref. [98]. For recent reviews see Refs. [1][2][3]. 14 Our sign convention assumes negative axial charge FA(0) ≡ gA < 0, hence the negative sign before B(Q 2 ). For antineutrino-proton scattering, this sign is positive. Our expression corresponds to Ref. [98] and differs from Ref. [99] in the axial form factor contribution to the function A.  (18) represents the "Born" cross section for the quasielastic process, analogous to Eq. (11) for the case of ep scattering. Soft radiation effects and two-boson exchange contributions have been subtracted and are to be treated separately. It is important to include such radiative corrections and to account for collinear and hard-photon emission in a practical experiment; however, our focus here is to determine the Born cross section for the quasielastic process.
The axial form factor F A is taken from Ref. [100]. In order to illustrate the utility of the new vector form factors, we will use the standard PCAC ansatz for the pseudoscalar form factor F P (whose effects are suppressed by powers of the lepton mass), The isovector vector form factors F V 1 , F V 2 are determined either by taking the difference of proton and neutron form factors, or by directly implementing the isospin-decomposed fit. We have ignored secondclass form factors in Eq. (20), and isospin-violating corrections to the relation of F p,n i to F V i . These effects are suppressed by the fine structure constant α or by (m d − m u )/Λ QCD , and are expected to be small compared to other uncertainties.
To illustrate the relevant range of Q 2 for neutrino beams in the GeV energy regime, we display the ν µ n CCQE cross section as a function of Q 2 , fixing E ν = 5 GeV, in the left-hand side of Fig. 1. The cross section is dominated by Q 2 1 GeV 2 , and is relatively insensitive to the detailed form factor behavior at larger momentum transfers. For comparison, the right hand side of Fig. 1 shows the ν τ n CCQE cross section; this rare process accesses somewhat larger Q 2 . In both cases, there can be residual sensitivity to higher-order coefficients in the z expansion that are poorly constrained by the chosen electron scattering data set. This sensitivity can be determined in practice by recomputing observables using different values of k max .
Our CCQE cross sections for muon (anti)neutrino are displayed as a function of neutrino energy in Fig. 2, using our default isospin-decomposed (iso 1 GeV 2 ) fit. The current large uncertainty of the axial form factor dominates the error budget. We remark that the deviation of central values between our fit and the commonly-used BBBA2005 model [94] is comparable to the axial form factor uncertainty at E ν 1 GeV. The cross section depends strongly on lepton flavor at energies near or below the muonproduction threshold, as shown in Fig. 3.   9.6(9) 6.47 (47) Tables 5 and 6 show the CCQE total cross sections at two benchmark points E ν = 1 GeV and E ν = 3 GeV. In addition to axial and vector form factor uncertainty, we include a z expansion truncation uncertainty estimated from the shift in central value when the default fit with k max = 8 is replaced by the fit with k max = 9. We also compare our evaluation with Ref. [100], where the BBBA2005 parameterization was used for vector form factors.

Spectroscopy of electronic and muonic atoms
Modern spectroscopy experiments with ordinary and muonic hydrogen [4,5,[101][102][103][104][105][106] are sensitive to the internal structure of the proton. In particular, the small size of muonic atoms enhances sensitivity to structure-dependent effects and makes measurements with muons attractive in searches for new physics and precise studies of proton and nuclear dynamics. The leading structure-dependent effect, which is proportional to r 2 E , shifts energy levels at order m α 4 and enters via the exchange of one virtual photon between the lepton ( = e or = µ) and the proton. This effect does not depend on the spin state of the energy level. The leading spin-dependent contribution of order m α 5 arises from the two-photon exchange. It contributes to the hyperfine splitting of energy levels [5]. Modern measurements of the Lamb shift in muonic hydrogen [4,5], of the hydrogen-deuterium isotope shift [107] and the 1S-2S transition in hydrogen [108] are sensitive to spin-independent two-photon exchange contributions as well. For both ordinary and muonic hydrogen, the bulk of the two-photon exchange contributions is determined by certain structure parameters, "moments", expressed as Q 2 integrals over products of elastic form factors. In this Section, we compute the Friar and Zemach radii governing spin-independent and spin-dependent two-photon exchange, respectively. Some previous results are compiled in Appendix C.

Lamb shift
The leading structure-dependent contribution to the Lamb shift is proportional to the (cube of the) Friar radius r F : where G E (0) = dG E /dQ 2 | Q 2 =0 . We evaluate r 3 F exploiting the fit of proton and neutron data as well as isospin-decomposed fits and present our results in Table 7. The first error is from the extracted form factor covariance matrix, and the second is the shift in central value when the default fit with k max = 8 is replaced by the fit with k max = 9. We note that removing the µH constraint from our default proton fit shifts (r p F ) 3 = 2.246(58) fm 3 → 2.97(35) fm 3 .

Hyperfine splitting
The first measurements of the 1S hyperfine splitting in muonic hydrogen with ppm precision are being planned by the CREMA [104], FAMU [105], and J-PARC [106] collaborations. The leading nucleonstructure contribution to the hyperfine splitting of S energy levels is given by the two-photon exchange diagram. The bulk of the correction is proportional to the Zemach radius r Z [109], which can be expressed as a convolution of nucleon electric and magnetic form factors, Similar to the Friar radii, we present Zemach radii evaluated using the fits from Sec. 2.5 in Table 8. These results provide a first rigorous error estimate. We note that removing the µH constraint from our default proton fit shifts r p Z = 1.0227(94) fm → 1.0426(132) fm.

Summary
We have presented a compact representation of the proton and neutron vector form factors in terms of z-expansion coefficients, including central values, errors and correlations. The results can be used to evaluate both central values and error bars for many derived quantities that are sensitive to GeV and sub-GeV momentum transfers. In our default fits we employed the following data: (i) the high-statistics Mainz data set for ep cross sections; (ii) en elastic scattering data at momentum transfers Q 2 ≤ 1 GeV 2 ; and (iii) precise external constraints on the proton and neutron electric charge radii. We considered two types of fits to this data. First, we performed separate proton and neutron fits, i.e., proton data fit to proton form factors and neutron data fit to neutron form factors. Second, we performed a fit of both proton and neutron data to isospin-decomposed form factors. For proton structure observables, there is only a slight reduction in uncertainty when the proton fit is replaced by the isospin-decomposed fit; for simplicity we use the proton fit as our final result: r p M = 0.739(41)(23) fm; r 4 p E = 1.08(28)(5) fm 4 ; r 4 p M = −2.0(1.7)(8) fm 4 ; (r p F ) 3 = 2.246(58)(2) fm 3 ; r p Z = 1.0227(94)(51) fm. For neutron structure observables, the abundance and precision of proton data relative to neutron data leads to a significant reduction in uncertainty when using the isospin-decomposed fit; we thus use the isospin-decomposed fit as our final result:  (77). We present the uncertainty coming from vector form factors and the truncation uncertainty as the last two errors respectively in results of this paragraph.
Significant tensions exist between the default data set and other world ep data. In Ref. [7], we quantified this tension as a function of Q 2 . Without knowing the source of discrepancy it remains unclear how to rigorously address this tension in the fit, and how to propagate it to derived observables. In this paper, we bypass this issue by focusing on the internally consistent default data set, but present results for comparison also from the global fit that includes "World" and "Pol" data as in Table 1. This combined iso (3 GeV 2 ) fit is similar to our global fit from Ref. [7]; the detailed comparison is discussed in Appendix A. The iso (3 GeV 2 ) fit includes more data than our default fit and it is thus not surprising that this fit predicts smaller uncertainties in derived observables. However, the fit does not address internal data set tensions and the iso (3 GeV 2 ) uncertainties are likely underestimates.
A primary goal of this work is to provide a consistent framework for applications such as neutrino event generators to propagate form factor constraints and uncertainties into cross section predictions. The framework is readily adapted to new data. Our new precise vector form factors have small uncertainty but deviate significantly from commonly-used parameterizations; such deviations become comparable to the dominant axial form factor uncertainty for larger neutrino energies. It is important to address these discrepancies with future experimental and/or lattice QCD data. We remark that the axial form factor was extracted under a specific (BBBA2005) assumption for the vector form factors. This ansatz can be justified given the current large uncertainty of elementary target neutrino data. However, correlations between vector and axial form factors should be accounted for when future more precise data becomes available.  Table 1 [top]; the n fits of lines 2 and 3 in Table 1 [middle]; and the iso (1 GeV 2 ) fit of line 4 in Table 1 [bottom]. The purple bands are the results of the iso (3 GeV 2 ) fit of line 5 in Table 1. The red dotted curves correspond to the global fit of Ref. [7], and the blue dash-dotted curves are the BBBA2005 result of Ref. [85]. Figure 4 compares the form factors from our default p, n, and iso (1 GeV 2 ) fits to our iso (3 GeV 2 ) fit. We also compare to our previous global fit from Ref. [7]. That global fit corresponds with the iso (3 GeV 2 ) fit after the following modifications: (i) inclusion of the µH constraint (9); (ii) improved treatment of TPE correction (12); (iii) omission of data above Q 2 = 3 GeV 2 ; (iv) choice of form factor expansion parameters t 0 and k max optimized for 0 < Q 2 < 1 GeV 2 . Note that the error band from Ref. [7] Figure 5: Comparison of G p E from fits with and without PRad data. In both plots, the black, long dashdotted curve is our default (proton) fit. On the left-hand side, blue points are the tabulated PRad form factors with statistical errors; the blue, dash-dotted curve is the PRad extraction; and the red, dotted curve is our extraction from PRad data. On the right-hand side, we compare our default fit to the fit when the µH constraint is replaced by PRad data (purple, dashed curve).

A.1 Mainz and World+Pol data sets
includes an ad hoc "data tension" error to account for the tension between Mainz and other World data. Since we have in mind applications to neutrino cross sections, we compare also to the commonly-used BBBA2005 parameterization. The BBBA2005 parameterization resulted from a fit to data preceding the A1 experiment, and is in severe tension with our default fit for G p M .

A.2 PRad and Mainz data sets
The PRad Collaboration recently presented new measurements of elastic electron-proton scattering at JLAB [110]. At two beam energies E = 1.1 GeV and 2.2 GeV, 33 and 38 measurements were taken in the range of Q 2 up to 0.016 GeV 2 and 0.058 GeV 2 , respectively. PRad announced a result r p E = 0.831 ± 0.007 stat ± 0.012 syst fm fitting to a rational functional form for G p E , Notably, this is within 1σ of r p E µH in Eq. (9) from muonic hydrogen spectroscopy. The PRad Collaboration employed particular assumptions in fitting the cross sections, which are detailed in the Supplementary Material of Ref. [110]. To extract form factors, the measured scattering cross sections were fit to the following reduced cross section where τ = Q 2 /(4M 2 p ), = [1 + 2(1 + τ ) tan 2 (θ/2)] −1 , n is a normalization parameter for a given beam energy, and G K M is the Kelly parameterization for the proton magnetic form factor [111]. The PRad Collaboration showed that the cross sections vary by less than 0.2% when different models for the magnetic form factor are used. 15 In Fig. 5, we compare G p E from four fits: 15 Note that after factoring out the normalization parameter from the reduced cross section, the ansatz in Eq. (25) does not strictly reproduce the correct anomalous magnetic moment. Since the parameter τ is small in the range of Q 2 covered by the PRad experiment, the fits are insensitive to the replacement G K M → G K M /n.
1. (blue, dash-dotted in left-hand plot) Fitting with the above rational functional form to the provided form factor tabulations with statistical-only errors, we have reproduced the PRad results for the proton charge radius and reduced χ 2 .
2. (red, dotted in left-hand plot) Following a modified version of the PRad procedure, we also fit directly to the tabulated PRad cross sections with statistical and systematic uncertainties added in quadrature. Fixing the magnetic form factor to the dipole form G p M (Q 2 ) = µ p G D (Q 2 ), we employed the z expansion (without sum rules) with k max = 3 for G p E . For each beam energy, a separate normalization parameter multiplying the entire reduced cross section is used. We did not apply additional TPE corrections. The extracted radius value is r p E = 0.836 (19) fm, with χ 2 = 23.88 and n dof = 68.
3. (black, long dash-dotted in both left-and right-hand plots) G p E obtained from our default protononly fit to Mainz data with r p E µH constraint.
4. (purple, solid/dashed in right-hand plot) G p E obtained from the proton-only fit using the z expansion with k max = 8 to combined Mainz and PRad data (statistical and systematic uncertainties added in quadrature) without r p E µH constraint. TPE corrections are applied to both data sets. The extracted radius value for this fit is r p E = 0.843(11) fm, with χ 2 = 503.24 and n dof = 720.
In the right-hand plot of Fig. 5, we show that the combined fit without the µH constraint and our default fit with the µH constraint lie within the 1σ uncertainty bands of each other for the entire Q 2 range of the PRad data. It would be desirable to include the PRad data directly into our fits, alongside the data in Table 1; we refrain from doing so since the PRad uncertainties are systematics dominated and uncertainty correlations are not yet available [110]. We have shown, however, that this data will not significantly alter our fits once the precise external µH constraint is imposed. Taking the PRad errors at face value (i.e., neglecting correlations), we remark that the z-expansion fit to PRad data results in a significantly larger uncertainty for r p E than is obtained using Eq. (24), comparable to the ∼ 0.020 fm uncertainty of our default proton fit, while the combined fit returns a radius uncertainty that is a factor of two smaller than either data set in isolation.

B Two-photon exchange
In this Appendix, we provide pertinent details of two-photon exchange corrections that were discussed in Section 2.3.
For momentum transfers Q 2 1 GeV 2 , dispersion relations have been used to constrain TPE corrections using available experimental data for inelastic cross sections. At relatively small momentum transfer and small scattering angles, the contribution from all inelastic intermediate states was evaluated [112] on top of the proton state [113] accounting for unpolarized proton inelastic structure functions in the resonance region. To calculate the TPE correction at large scattering angles, the data-driven dispersion relation framework was recently developed [88][89][90][91][92]. The imaginary part of TPE amplitudes is evaluated from on-shell information in the physical region of electron-proton scattering. The real part of TPE amplitudes requires information from the unphysical region as input. Novel methods of analytical continuation [89,92] allow us to overcome this complication. Contributions from proton and πN intermediate states are evaluated for Q 2 1 GeV 2 . At low momentum transfer and backward scattering angles, the relative contribution of inelastic intermediate states is found to be much smaller than the elastic contribution to TPE. At larger electron beam energies and momentum transfer, the intermediate states with higher invariant mass, e.g., ππN , become kinematically enhanced and prevent making a rigorous prediction in the absence of exclusive experimental data. At small momentum transfer Q 2 0.25 GeV 2 and scattering angles, we account for all inelastic intermediate states [112]. At large angles and momentum transfer, proton and πN states are included [89,90,92]. The intermediate region is described by interpolation between these two calculations as in Ref. [93]. We denote this dispersive result as δ dispersive and provide the corresponding correction for each point in the Mainz data set in the Supplementary Material.
At larger momentum transfers, Q 2 1 GeV 2 , the explicit form of the phenomenological TPE modification is as follows [85]: which is negative (since 0 ≤ ε ≤ 1) and increases the inferred Born cross section. As discussed in the main text, this correction serves to improve agreement between polarization measurements and TPE-corrected unpolarized Rosenbluth measurements at high-Q 2 .

C Comparison to literature
In this Section, we provide some existent results for form factor curvatures, Friar radii and Zemach radii.