Nucleon form factors in dispersively improved Chiral Effective Field Theory II: Electromagnetic form factors

We study the nucleon electromagnetic form factors (EM FFs) using a recently developed method combining Chiral Effective Field Theory ($\chi$EFT) and dispersion analysis. The spectral functions on the two-pion cut at $t>4 M_\pi^2$ are constructed using the elastic unitarity relation and an $N/D$ representation. $\chi$EFT is used to calculate the real functions $J_\pm^1 (t) = f_\pm^1(t)/F_\pi(t)$ (ratios of the complex $\pi\pi \rightarrow N \bar N$ partial-wave amplitudes and the timelike pion FF), which are free of $\pi\pi$ rescattering. Rescattering effects are included through the empirical timelike pion FF $|F_\pi(t)|^2$. The method allows us to compute the isovector EM spectral functions up to $t \sim 1$ GeV$^2$ with controlled accuracy (LO, NLO, and partial N2LO). With the spectral functions we calculate the isovector nucleon EM FFs and their derivatives at $t = 0$ (EM radii, moments) using subtracted dispersion relations. We predict the values of higher FF derivatives with minimal uncertainties and explain their collective behavior. We estimate the individual proton and neutron FFs by adding an empirical parametrization of the isoscalar sector. Excellent agreement with the present low-$Q^2$ FF data is achieved up to $\sim$0.5 GeV$^2$ for $G_E$, and up to $\sim$0.2 GeV$^2$ for $G_M$. Our results can be used to guide the analysis of low-$Q^2$ elastic scattering data and the extraction of the proton charge radius.

The electromagnetic form factors (EM FFs) describe the nucleon's elastic response to external EM fields and reveal the spatial distribution of charge and magnetization inside the strongly interacting system. They are among the most basic characteristics of nucleon structure and have been the object of extensive theoretical and experimental study. The electric and magnetic FFs at invariant momentum transfers t < 0 are measured in elastic electron-nucleon scattering. Experiments at |t| ≤ 1 GeV 2 have been performed at many facilities, most recently at the Mainz Microtron (MAMI) [1][2][3] and at Jefferson Lab (JLab) [4][5][6]; see Refs. [7,8] for a review of the other data. The derivative of the proton electric FF at t = 0, or charge radius, governs the nucleon structure corrections to the energy levels of hydrogen atoms (electronic or muonic) and is measured with high precision in atomic physics experiments. Recent experimental results have raised interesting questions regarding the precise value of the proton charge radius and the extrapolation of the elastic scattering data to t → 0; see Refs. [9][10][11] for a review. Theoretical calculations of the nucleon FFs at |t| 1 GeV 2 are needed to guide the analysis of the experimental data and help answer these questions. Dedicated measurements of the proton electric FF at extremely low momentum transfers |t| ≥ 10 −4 GeV 2 are planned at JLab [12]. Knowledge of the low-t EM FFs is also required for constructing the peripheral transverse densities and generalized parton distributions (GPDs) of the nucleon [13,14].
In a previous article we have proposed a method for calculating the nucleon FFs of G-parity-even operators by combining Chiral Effective Field Theory (χEFT) and dispersion analysis (dispersively improved χEFT, or DIχEFT) [15]. It starts from the dispersive representation of the FFs as analytic functions of t and constructs the spectral functions on the two-pion cut at t > 4M 2 π using the elastic unitarity relation [16,17]. An N/D representation is employed to separate the coupling of the ππ system to nucleon from the effects of ππ rescattering in t-channel. χEFT is used to calculate the real functions J(t) = f (t)/F π (t) -the ratios of the complex ππ → NN partial-wave amplitudes (PWAs) and the timelike pion FF, which are free of ππ rescattering and show good convergence. Rescattering effects are included through the timelike pion FF |F π (t)| 2 , which is taken from sources outside of χEFT [experimental data, Lattice QCD (LQCD)]. The formulation permits firstprinciples calculations of the two-pion spectral functions of the FFs with controlled accuracy. The spectral functions can then be used to evaluate the FFs and related quantities of interest (nucleon radii, transverse densities) with subtracted dispersion relations. The method results in a dramatic improvement compared to conventional χEFT calculations of nucleon FFs and their spectral functions [18][19][20][21][22]. In Ref. [15] the method was applied to the nucleon isoscalar-scalar FF, where the ππ system in the t-channel is in the I = J = 0 state. The resulting spectral functions and FFs were found to be in good agreement with those of empirical dispersion theory and Roy-Steiner equations [23,24].
In the present article we use the DIχEFT method to calculate the nucleon EM FFs at low t and study their properties. We construct the isovector EM spectral functions on the two-pion cut by combining the elastic unitarity relation in the I = J = 1 channel, the N/D representation; χEFT calculations of the J functions at LO, NLO and partial N2LO accuracy; and the timelike pion FF |F π | 2 measured in e + e − annihilation experiments. Realistic spectral functions are obtained up to t ∼ 1 GeV 2 , which includes the ρ meson region essential for EM structure. With these spectral functions we evaluate the isovector FFs and their derivatives (radii) using subtracted dispersion relations. We obtain the individual proton and neutron FFs by supplementing the calculated isovector spectral functions with an empirical parametrization of the isoscalar ones [14]. Excellent agreement with the low-t EM FF data is achieved.
In particular, the method allows us to predict the higher derivatives of the EM FFs at t = 0 and explain their collective behavior. They are given by wellconvergent dispersion integrals, which can be evaluated reliably with the DIχEFT spectral functions, with minimal model dependence. The higher derivatives are governed by two disparate dynamical scales -the vector meson mass, M 2 V (V = ρ, ω), and the two-pion threshold, 4M 2 π -and exhibit a rich structure. The values of the higher derivatives therefore differ qualitatively from what one would estimate based on a single dynamical scale (naturalness). Recent attempts to fit the low-t proton electric FF data and extract the charge radius have engendered a debate regarding the values of higher FF derivatives and their role in the t → 0 extrapolation [1,2,[25][26][27][28][29][30]. Our predictions for the higher derivatives can be compared with those obtained in form factor fits (regarding order-of-magnitude, collective behavior) and used to discriminate between different fits. Our method incorporates the exact analytic structure of the FF in t and the multiple dynamical scales governing its behavior, which are essential in the analysis of low-t FF data and the extraction of charge radius.
The plan of the article is as follows. In Sec. II we describe the steps of the DIχEFT calculation, expanding on the general description of the method in Ref. [15] and emphasizing the aspects that are new or specific to the EM FFs. This includes the unitarity relations and N/D representation in the I = J = 1 channel, the LO χEFT calculation of the J-functions, the estimate of higherorder corrections, and the parametrization of the timelike pion FF. In Sec. III we present the results and their interpretation. This covers the isovector nucleon EM spectral functions, the nucleon EM radii, higher derivatives (moments) of the EM FFs and their structure, and the spacelike nucleon FFs at low |t|. In Sec. IV we summarize the results and comment on possible further applications of the method.
A similar method for calculating nucleon FFs, combining χEFT and dispersion theory, was described in Ref. [31] and applied to the EM FFs in Ref. [32]. The differences from our approach are mainly in the technical implementation of the ππ rescattering effects in the unitarity relations (N/D method vs. Omnes function) and the estimates of higher-order chiral corrections. Nucleon FFs were also studied in an extension of χEFT with explicit vector meson degrees of freedom in Ref. [33]. The low-t nucleon FFs and their derivatives were also calculated in heavy-baryon χEFT with ∆'s, in the context of a study of two-photon exchange corrections to muonic hydrogen in Ref. [34].

A. Nucleon EM form factors in DIχEFT
The transition matrix element of the EM current between nucleon (proton, neutron) states is parametrized by the invariant FFs F 1 (t) and F 2 (t) (Dirac and Pauli FFs; we use the conventions of Ref. [14]). The electric and magnetic FFs are defined as where τ ≡ −t/(4m 2 N ). At zero momentum transfer the electric FF gives the electric charge of the nucleon, G p,n E (0) = (1, 0), and the magnetic FF gives the total (pointlike plus anomalous) magnetic moment in units of nuclear magnetons, G p,n M (0) = µ p,n = (2.79, −1.91). The isovector and isoscalar components are defined as The FFs are analytic functions of t and satisfy dispersion relations. They express the FFs at complex t as an integral over their singularities at real t > 0, corresponding to (unphysical) timelike processes in which the current couples to the nucleon through exchange of a hadronic system in the t-channel. In the case of the isovector EM FFs the lowest-mass hadronic state is the ππ state, and the dispersion integrals start at t = 4M 2 π (two-pion threshold) The real functions Im G V i (t ) (t > 4M 2 π ) are referred to as the spectral functions. At 4M 2 π < t < 16M 2 π the ππ state is the only possible state contributing to Im G V i (t ). It is known from e + e − annihilation experiments that higher states (4π etc.) do not couple strongly to the EM current, and it is expected that the ππ channel in the nucleon spectral functions remains dominant up to t ∼ 1 GeV 2 . In this region the ππ channel includes the ρ resonance at t = M 2 ρ ∼ 0.6 GeV 2 , which has a decisive influence on the EM FFs. In the isoscalar case the dispersion integral starts with the 3π state at t = 9M 2 π , with the dominant strength at t < 1 GeV 2 coming from the ω resonance at t = M 2 ω ∼ 0.6 GeV 2 . The contribution of higher-mass states to the isovector and isoscalar FFs is constrained by the total charge and magnetic moment and has been determined empirically through fits of spacelike FF data [35,36]; the exact composition of these multi-hadron states is poorly known and will not be needed in the following applications.
In the region 4M 2 π < t < 16M 2 π the isovector spectral functions on the two-pion cut can be obtained from the elastic unitarity relations [16,17] where k cm ≡ t/4 − M 2 π is the center-of-mass momentum of the ππ system in the t-channel, f 1 ± (t) are the ππ → NN PWAs in the normalization of Ref. [35], and F * π (t) is the complex-conjugate timelike pion EM FF. Equations (5) and (6) are valid strictly in the region up to the 4π threshold, 4 M 2 π < t < 16 M 2 π ; if contributions from 4π and higher states are neglected they can effectively be used up to t ∼ 1GeV 2 = 50 M 2 π . The expressions on the right-hand side of Eqs. (5) and (6) are real because the complex functions f 1 ± (t) and F π (t) have the same phase on the two-pion cut (Watson theorem) [37]. The unitarity relations can therefore be written in a manifestly real form as [16,17,38] where The functions J 1 ± (t) are real for t > 4M 2 π and thus have no right-hand cut; their only singularities are left-hand cuts at t < 4M 2 π − M 4 π /m 2 N , the threshold resulting from the singularity of the nucleon Born term in the ππ → NN PWAs. Equations (7)-(9) are equivalent to a particular N/D representation of the PWAs [39], in which the numerator functions J 1 ± (t) contain the lefthand cut and the denominator function 1/F π (t) contains the right-hand cut.
To evaluate the spectral functions in the representation of Eqs. (7)-(9), following Ref. [15], we calculate the real functions J 1 ± (t) in χEFT and multiply them with the empirical pion FF modulus |F π | 2 . Advantages of this approach are: (a) The χEFT calculation of J 1 ± (t) is free of ππ rescattering and shows good convergence. Rescattering effects are entirely contained in |F π (t)| 2 , which is taken from other sources. In traditional "direct" χEFT calculations of the spectral functions the ππ rescattering effects result in large higher-order corrections and render the perturbative expansion impractical. (b) The functions J 1 ± (t) are dominated by the scales M π and m ∆ − m N associated with the Born graph singularities, while |F π (t)| 2 dominated by the chiral-symmetrybreaking scale Λ χ . The representation Eqs. (7)-(9) is therefore consistent with the idea of separation of scales. (c) The squared modulus |F π (t)| 2 can be imported directly from the e + e − → π + π − data or from LQCD calculations without determination of the phase [15,40]. For further discussion of the method we refer to Ref. [15]. contributions from the ∆ Born term. The setup of the χEFT calculation is described in Ref. [41] and summarized in Ref. [14] (fields, chiral Lagrangian, power counting, values of couplings). The interactions of the spin-3/2 ∆ field are formulated with consistent vertices [42][43][44][45], and the extended-on-mass-shell (EOMS) scheme is used to maintain the standard power counting [46] (diagrams with pion loops do not enter in the present calculation). The LO diagrams contributing to the ππ → NN PWAs in the I = J = 1 channel are shown in Fig. 1. They include the N and ∆ Born terms, Fig. 1a and c, and the Weinberg-Tomozawa contact term, Fig. 1b, which appears as the result of chiral invariance of the dynamics with relativistic baryons. We take the results for the LO πN → πN amplitude of Ref. [41] (the first relativistic χEFT calculation of πN scattering with explicit ∆), and project onto the I = J = 1 channel to get the PWAs f 1 ± (t). The pion EM FF at this order is just F π = 1 (pointlike), see Fig. 1d. The χEFT results for J 1 ± (t) therefore coincide with f 1 ± (t) at this order. Analytic expressions for J 1 ± (t) are given in Appendix A. Numerical results of the LO approximation will be shown below.

C. Estimates of higher-order corrections
At NLO accuracy, corrections to the I = J = 1 ππ → NN PWAs arise from the NLO contact term in the chiral Lagrangian with the low-energy constant (LEC) c 4 . The value of this LEC has to be adjusted consistently with our unitarity-based approach [14]. In standard χEFT calculations c 4 receives large contribution from ρ meson exchange. Since in our formulation the effect of the ρ is included explicitly through |F π (t)| 2 , we have to remove it from the value of c 4 to avoid double-counting (see Fig. 2). Using the estimate for the ρ contribution of Ref. [47], c ρ 4 ≈ 1.63 GeV −1 , and subtracting it from the empirical c 4 reported in Refs. [41,48], we obtain the range Note that these values are much smaller than the original c 4 and consistent with zero, which means that the NLO corrections to the isovector spectral functions in our formulation are very small. The analytic expressions for the NLO corrections to J 1 ± (t) are given in Appendix A.
At N2LO accuracy pion loop corrections appear, and the structure of the χEFT expressions becomes considerably richer. The ππ → NN PWAs and the pion FF now involve ππ rescattering in the t-channel and become complex at t > 4M 2 π , in such a way that their phases cancel and the functions J 1 ± (t) of Eq. (9) remain real. Also, πN and π∆ s-channel intermediate states appear in the πN amplitude. Following Ref. [14] we perform a simple estimate of the N2LO corrections to the spectral functions of the electric FF, by taking the N2LO tree-level amplitudes and fixing the LECs through the charge sum rule. We require that the unsubtracted dispersion relation for the isovector electric FF reproduce the isovector charge when the integration is restricted to the region This condition gives N2LO contact term contributions with sign opposite to that of the LO and NLO results, which provides a crucial curvature in the electric spectral function and allows us to extend the calculations up to t ∼ 1 GeV 2 . In the language of traditional dispersion analysis these contact terms represent the negative contributions from the ρ , which compensate the excess charge that would be produced by the ρ alone. In the magnetic FF no new tree-level amplitudes with LECs enter at N2LO level, so that the described method of estimating of the corrections cannot be applied. Our calculations of G V M are therefore limited to NLO accuracy, and are expected to describe the empirical spectral functions only at t 1 GeV 2 . Numerical results of the NLO and N2LO approximations will be shown below.

D. Timelike pion EM form factor
For evaluating the timelike pion FF entering in our calculation we use a Gounaris-Sakurai parametrization of the e + e − → π + π − exclusive annihilation data [49], including effects of ρ-ω mixing [50], with the parameters determined in Ref. [25]. The squared modulus |F π (t)| 2 is shown in Fig. 3. One clearly sees the ρ resonance at t ∼ 0.6 GeV 2 and the rapid variation resulting from ρω mixing. The fact that |F π (t)| 2 reaches a value of ∼2 at t ∼ 0.2 GeV 2 , and ∼ 10 at t ∼ 0.4 GeV 2 , shows that ππ rescattering is very substantial already at moderate t and justifies our approach of incorporating these effects empirically.
Since |F π (t)| 2 is determined very accurately from the annihilation data we neglect the effect of its uncertainty on the spectral functions. In the following we quote only the uncertainties of the spectral function resulting from the χEFT calculation of J 1 ± (t).

A. Isovector EM spectral functions
The spectral functions of the isovector EM nucleon FFs are evaluated using the DIχEFT method and parameters described in Sec. II. Figure 4 shows the results of the χEFT calculation of the real functions J 1 ± (t) at t > 4M 2 π , Eq. (9). For a better view the plot shows the functions multiplied by the kinematic factors of Eqs. (7) and (8) LO and NLO results are close because the adjusted LEC c 4 is small, Eq. (11). In J 1 + (t) the N2LO corrections, estimated as described in Sec. II C, are negative and cause the function to decrease and turn negative at larger t. (c) The χEFT results for J 1 ± (t) show reasonable agreement with the functions extracted from an analysis of πN scattering data using Roy-Steiner equations [51]. In both J 1 + and J 1 − the LO and NLO approximation agree with the Roy-Steiner result up to t ∼ 0.2 GeV 2 . In J 1 + the negative N2LO corrections extend the region of agreement up to t ∼ 1 GeV 2 . Figure 5 shows the isovector spectral functions, obtained by multiplying the χEFT results for J 1 ± (t) with the empirical |F π (t)| 2 cf. Eqs. (7) and (8). The spectral functions clearly show the effects of ππ rescattering, which are not suitable for perturbative χEFT treatment and included through the empirical pion FF in our approach. Note that the enhancement through |F π (t)| 2 is large even near the two-pion threshold t ∼ 4 M 2 π , cf. Sec. II D and Fig. 3. The convergence pattern of the spectral functions follows from that of the χEFT calculation of J 1 ± (t). In
Dashed lines: LO approximation. Blue bands: NLO approximation. Red band: NLO+N2LO, estimated as described in Sec. II C. Solid orange lines: Roy-Steiner analysis results [51]. . both Im G E and Im G M , the LO and NLO approximations are in good agreement with the Roy-Steiner results up to t ∼ 0.2 GeV 2 . In Im G E the negative N2LO correction (estimated) is sufficient to reproduce the Roy-Steiner result up to t ∼ 1 GeV 2 .

B. Nucleon EM radii
The nucleon's isovector electric and magnetic radii are given by the dispersion integrals The factor 1/t 2 ensures convergence of the integral over the range t 1 GeV 2 ; see Fig. 6. The integrals can therefore be evaluated with the DIχEFT spectral functions. Table I  LO and NLO spectral functions are larger than the Roy-Steiner result around the ρ peak. For the magnetic radius our result is larger by a factor ∼2 than the phenomenological and LQCD results, because our magnetic spectral function is likewise too large around the ρ peak. DIχEFT allows us to calculate the isovector nucleon EM FFs and radii, which are matrix elements of Gparity even operators. Experiments measure the individual proton and neutron FFs and radii. In view of the questions concerning the proton charge radius measurements it is interesting to compare our results directly with the experimental results for the proton and neutron charge radii. To do so, we supplement the DIχEFT results for the isovector spectral functions with an empirical parametrization of the isoscalar spectral functions. We use a two-pole parametrization with the ω pole at M 2 ω = 0.61 GeV 2 and a second pole at M 2 2 ≈ M 2 φ = 1 GeV 2 [14] Im where the coefficients a ω i (including their uncertainties) are taken from the dispersive FF fit of Ref. [36], and the coefficients a (2) i are adjusted to reproduce the total charge and magnetic moment. The second pole is an effective pole representing the overall strength of the spectral function at t ∼ 1 GeV 2 ; the details of the strength distribution at these values of t are not important for estimating the nucleon radii. The proton and neutron radii obtained in this way are summarized in Table II. In the proton and neutron electric radii, the estimated uncertainty is dominated by the isoscalar component, which is purely empirical. Our results obtained with the NLO+N2LO DIχEFT calculation of the isovector radii are in agreement with the experimental values. In the magnetic radii the estimated uncertainty is likewise dominated by the isoscalar component. Note that the DIχEFT calculation of the isovector magnetic spectral function does not include the N2LO corrections and strongly overestimates the empirical result (cf. Table I); this discrepancy is not reflected in the uncertainty estimate.

C. Higher derivatives of EM form factors
Higher derivatives of the nucleon EM FFs at t = 0 are of interest for several reasons. In the experimental analysis, the values of higher derivatives allowed (or assumed) in fits of FF data at t < 0 directly affect the extrapolation to t = 0 and extraction of the nucleon charge radii; see Refs. [26][27][28][29][30] for details. In the theoretical studies reported here, higher derivatives of the FFs represent clean chiral observables that can be predicted almost model-independently with minimal uncertainties. The comparison of low-and high-order derivatives reveals the presence of two dynamical scales in the nucleon FFs, which implies a surprisingly rich structure and should be incorporated into the experimental analysis.
In the context of the traditional representation of the FFs as Fourier transforms of 3-dimensional spatial densi- ties, the higher derivatives of the FFs at t = 0 correspond to the higher r 2 -weighted moments of the densities. The connection is given by [29] G E (t) = 1 + r 2 (for either p or n).
Note that for the proton and neutron the magnetic radii are defined as the derivatives of the FFs divided by the magnetic moments; this is not the case for the isovector and isoscalar components. In the following we quote results for the moments; they can be converted to FF derivatives through Eqs. (17) and (18). The higher moments of the isovector FFs are given by the dispersion integrals The factors 1/t n+1 strongly suppress contributions from large t and render the integrals well convergent. The integrals can therefore be evaluated accurately using the DIχEFT spectral functions. Figure 7 shows the integrands for the isovector electric FF derivatives with π (the peak of the spectral function), and the twopion threshold 4M 2 π (the start of the spectral integral). The integrals receive contributions from both regions of t , and their relative importance changes with n. For a rough assessment we can take M 2 ρ /2 = 0.3 GeV 2 as the boundary between the two regions. For n = 1 approximately 2/3 of the integral comes from the region t > M 2 ρ /2, and 1/3 from 4M 2 π < t < M 2 ρ /2. For n = 2, each region contributes about 1/2. For n = 3 and higher, the near-threshold region dominates.
The presence of two dynamical scales in the isovector moments can also be demonstrated by considering the ratios of successive moments, If the dispersion integral were dominated by a certain region of t , the value of the ratio Eq. (20) would be given by the average of 1/t over that region. The ratios thus directly reveal the effective values of 1/t in the integral. Figure 7 shows the ratios of the isovector electric FF moments obtained with the DIχEFT spectral functions. One sees that the ratios start with a value ∼ 1/m 2 ρ at n = 1 and increase to values ∼ 1/(4M 2 π ) at large n. The presence of two dynamical scales implies that the higher FF moments are of "unnatural" size, i.e.  representation and is insensitive to the details of the dynamical calculation presented here. It has has numerous consequences for the interpretation of the FF moments and the analysis of low-t elastic scattering experiments, which will be elaborated below. Table III shows the DIχEFT results for the higher moments of the isovector FFs, G V E and G V M . Because the dispersion integrals with n ≥ 2 sample the spectral functions near threshold, the higher moments can be computed accurately and represent genuine predictions of our approach. This is seen in the intrinsic uncertainty estimates of Table III: with increasing n, the derivatives become less sensitive to higher-order chiral corrections. We emphasize that one should be careful in interpreting the numerical values of the individual moments in Table III, as they contain large factorial factors. The unnatural behavior of the higher moments should be demonstrated by taking ratios (see above) or comparing the moments to a reference FF (see below).
It is worth noting that the ππ rescattering effects included through the timelike pion FF play an important role even in the higher FF moments. The moments with n ≥ 3 receive most of their contributions from the region 4M 2 π < t 10 M 2 π of the spectral integral Eq. (19) (see Fig. 7). The value of |F π (t )| 2 is ∼ 1.3 at t = 4M 2 π , and ∼ 2 at 10M 2 π (see Fig. 3). The enhancement compared to traditional direct χEFT calculations of the spectral functions is therefore quite substantial (note that without the factor |F π | 2 our results at LO and NLO would be identical to those of the direct calculation). Nevertheless our results for the higher moments agree with those of the direct χEFT calculation of Ref. [34] within errors.
In the isoscalar FFs the strength of the spectral functions is located overwhelmingly at the ω meson mass. The higher moments are therefore governed by this single scale and are of natural size. This in turn implies that the higher moments of the proton and neutron FFs are dominated by the isovector component and can be inferred from our DIχEFT results. the empirical parametrization of the isoscalar FFs. We stress that the isoscalar information is used here only to demonstrate that the higher derivatives are dominated by the isovector component, and that the uncertainties associated with the isoscalar parametrization are irrelevant in the higher derivatives. The theoretical results described here have implications for the analysis of electron-proton elastic scattering data at low Q 2 ≡ −t and the extraction of the proton charge radius. The overall behavior of G p E in the region 0 < Q 2 1 GeV 2 is associated with a scale of the order of the vector meson mass M 2 V (V = ρ, ω). The first derivative of G p E at Q 2 = 0 is of the order 1/M 2 V and therefore appears natural, i.e., simple single-scale parametrizations of the finite-Q 2 data give a reasonable estimate of the first derivative. The higher derivatives, however, are governed by the scale 1/(4M 2 π ) and appear unnatural. Single-scale parametrizations or "natural" powers of the first derivative the data give qualitatively wrong estimates of the higher derivatives. To illustrate the point we compare the order-of-magnitude of the higher derivatives obtained from DIχEFT with the ones of the dipole parametrization which provides a good overall description of the FF data at 0 < Q 2 1 GeV 2 with Λ 2 ≈ 0.71 GeV 2 . Figure 9 shows the ratios as obtained with the DIχEFT results. The coefficients c n are defined such that for the dipole FF Eq. (21) the ratio is equal to unity for all n. The ratio Eq. (22) therefore indicates how strongly the actual higher derivatives deviate from the single-scale estimate based on the dipole form. One sees that the ratio is ∼10 2 for n = 4, and reaches values ∼10 5 for n = 8. It shows the striking consequences of the two dynamical scales in the higher FF derivatives, as implied by their dispersive representation. A similar observation regarding the ratio of the actual FF moments to the dipole was made in the context of an empirical analysis in Ref. [58]. The values of the higher derivatives of G p E and their impact on the Q 2 → 0 extrapolation are presently the subject of intense discussions [26][27][28][29][30]. Fits to the low-Q 2 FF data with different classes of functions (polynomials, rational functions) give widely different values of the second and higher derivatives; see Table V for a compilation of recent results. The DIχEFT results are broadly consistent with the range of empirical values. An analysis of FF data incorporating theoretical constraints from DIχEFT will be the subject of a future study [59]. For reference we quote in Appendix B the numerical values of the DIχEFT moments of G p E up to n = 20. While the individual values have little physical significance, their order-of-magnitude and collective behavior could be compared with the pattern and observed in higher-order polynomial fits.
The unnatural behavior of the higher FF derivatives is a consequence of analyticity and the singularities of the πN Born amplitudes, which govern the behavior of the spectral function in the near-threshold region. Dispersive fits to the low-Q 2 FF data correctly implement these features and can provide reliable results for the higher FF derivatives [25,36].

D. Spacelike EM form factors
The DIχEFT approach also allows us to calculate the nucleon FFs at finite t < 0, where they are measured in eN elastic scattering experiments. For the isovector FFs we use the twice-subtracted dispersion relations Here the FFs at t = 0 (charge and magnetic moment) and the first derivatives (electric and magnetic radii) are taken as input, and the dispersion relation predicts the tdependence starting from the second order. The integrals in Eq. (23) are well convergent and can be evaluated with the DIχEFT spectral functions. For the isoscalar FF we use the empirical parametrization Eq. (14), which imposes the correct value of the FF at t = 0. Combining the two we predict the individual proton and neutron FFs. Figure 10 summarizes the results. The estimated uncertainties are dominated by those of the isoscalar component, for which we have only the empirical parametrization. One observes: (a) In the electric FFs G p E and G n E , the LO and NLO approximations describe the experimental data only up to Q 2 ∼ 0.1 GeV 2 , while the NLO + N2LO approximations show good agreement with the data up to Q 2 ∼ 0.5 GeV 2 . This reflects the improvement of the isovector electric spectral function due to the partial N2LO corrections; see Figs. 4 and 5. (b) In the magnetic FFs G p M and G n M , our results describe the data up to Q 2 ∼ 0.2 GeV 2 . In this channel the N2LO corrections cannot be estimated using the method of Sec. II C. Altogether we obtain a very satisfactory description of the nucleon EM FFs with our dynamical approach.

IV. SUMMARY
This work reports a study of the nucleon EM FFs at momentum transfers |t| 1 GeV 2 using a new method combining χEFT and dispersion analysis (DIχEFT). The isovector spectral functions on the two-pion cut are constructed through the elastic unitarity condition. The N/D method is used to separate effects of ππ rescattering from the coupling of the ππ system to the nucleon.
χEFT is employed to calculate the real functions J 1 ± (t) describing the ππ coupling to the nucleon, which are free of ππ rescattering effects, resulting in good convergence. ππ rescattering effects are included through the timelike pion FF |F π (t)| 2 , which can be taken from e + e − annihilation data or LQCD calculations. The new organization is consistent with basic principles of χEFT and represents a major improvement over traditional direct calculations of the spectral functions. It allows us to calculate the isovector spectral functions up to t ∼ 1 GeV 2 (including the ρ meson region) with controlled accuracy. With these spectral functions we are able to evaluate elements of the FFs (radii, higher derivatives, t-dependence in spacelike region) using well-convergent subtracted dispersion relations.
The new method permits a realistic description of the low-t nucleon FFs and their derivatives. While the basic features of the FFs are rooted in analyticity and have been studied earlier in empirical dispersion theory, the new element is that the spectral functions can now be computed in a χEFT framework with controlled accuracy. It makes it possible to represent the information content of the nucleon FFs in the form of a few physical masses and LECs, resulting in a significant reduction of complexity. It also enables new interpretations of FFs in terms of spatial densities [13,14] and a space-time picture of chiral processes [60][61][62].
Our study shows that the derivatives of the EM FFs involve two dynamical scales. The first derivative is governed by the scale 1/M 2 V , while higher derivatives are governed by the scale 1/(4M 2 π ) 1/M 2 V and therefore appear unnaturally large. The rich structure attests to the fact that, through analyticity and the dispersion relations, the FFs of the nucleon are connected to its hadronic couplings and excitation spectrum and reflect the multiple dynamical scales characterizing the latter. The DIχEFT calculations provide an explicit realization of this general feature. While the predicted pattern of the higher moments is likely far beyond what can be extracted from experimental data, our analytic functions can be used in Monte-Carlo simulations to study how well the first few moments can be extracted under realistic conditions. Our numerical estimates of the higher derivatives can be used in an empirical analysis based on bounded least-squares regression. An analysis incorporating theoretical constraints from DIχEFT is in progress [59]. The DIχEFT FF calculations described here could be extended in several directions. The method could be applied to the N -∆ transition FFs as well as the EM FFs of the ∆ itself, which are defined rigorously in the context of S-matrix theory (as poles in the N → πN and πN → πN EM transition amplitudes) and have been studied in relativistic χEFT [63,64]. The method could also be applied to nucleon FFs of other G-parity-even operators, such as the energy-momentum tensor or higher moments of the GPDs. Finally, one might contemplate extending the DIχEFT approach to nucleon FFs of Gparity odd operators with a 3-pion cut, using methods of 3-body elastic unitarity that are presently being developed for the analysis of meson decays and LQCD calculations [65,66].
The present study demonstrates the potential of χEFT to yield fully predictive results for conventional nucleon structure observables. The same approach can be applied to hadronic structure elements appearing in searches for Physics Beyond the Standard Model; see e.g. Refs. [67][68][69].
are, respectively, the physical pion CM momentum and the unphysical nucleon CM momentum in the ππ → NN process. The functions resulting from the Weinberg-Tomozawa contact term (Fig. 1b) and the N Born term (Fig. 1a) are (A7) The contributions of the ∆ Born term (Fig. 1c) are x ∆ ≡ 2k cm p cm The functions F and G appearing in the first terms of Eqs. (A8) and (A9) are [15] F ≡ α(m ∆ + m N ) + β 3 (m ∆ − m N ), (A12) (A15) they are the invariant amplitudes of πN scattering at t > 4M 2 π and s = m 2 ∆ in the conventions of Ref. [62]. The functions D ∆+ and D ∆− appearing in the second terms in Eqs. (A8) and (A9) are The J functions resulting from the N and ∆ Born terms have logarithmic left-hand cuts starting at The singularity results from the intermediate baryon lines going on mass shell and corresponds to the left-hand cut of the ππ → NN PWA. The singularity is contained in the inverse tangent functions in Eqs. (A4), (A5), (A8), and (A5), which have logarithmic branch points at x 2 N,∆ = ±i. The J functions do not have a righthand cut at t > 4M 2 π , in accordance with their definition within the N/D method, Eqs. (9) and (10). While the expressions in Eqs. (A4), (A5), (A8), and (A9) contain prefactors with inverse powers of k cm , they are in fact regular in the limit k cm → 0, because the expressions in parentheses/brackets depending on x N or x ∆ vanish in the limit, x N,∆ = O(k cm ). Further properties of the J functions are discussed in Ref. [15].
The masses and coupling constants used in evaluating the LO expressions are the standard values for the SU (2) flavor group [15]: M π = 139 MeV, f π = 93 MeV, m N = 939 MeV, g A = 1.27, and m ∆ = 1232 MeV, h A = 2.85.
The contributions of the NLO contact term in the πN amplitude are (A20) The value of the LEC c 4 , determined by the procedure described in Sec. II C, is given in Eq. (11). For reference we present in Table VI our numerical estimates of the higher moments of the proton electric FF, obtained by combining the DIχEFT calculation of the isovector moments with the empirical estimate of the isoscalar moments based on Eq. (14). While individual higher moments have little physical significance and cannot realistically be extracted from the data, the order-ofmagnitude and collective behavior of our results could be compared with the patterns observed in fits to FF data [1,2,[25][26][27][28][29][30]57].