Dispersion-theoretical analysis of the electromagnetic form factors of the nucleon: Past, present and future

We review the dispersion-theoretical analysis of the electromagnetic form factors of the nucleon. We emphasize in particular the role of unitarity and analyticity in the construction of the isoscalar and isovector spectral functions. We present new results on the extraction of the nucleon radii, the electric and magnetic form factors and the extraction of $\omega$-meson couplings. All this is supplemented by a detailed calculation of the theoretical uncertainties, using bootstrap and Bayesian methods to pin down the statistical errors, while systematic errors are determined from variations of the spectral functions. We also discuss the physics of the time-like form factors and point out further issues to be addressed in this framework.


Introduction
Nucleons and electrons are the constituents of everyday matter with nucleons accounting for essentially all of its mass. The mass of the nucleon as a bound state of quarks and gluons, on the other hand, arises from the complicated strong interaction dynamics of quarks and gluons in Quantum Chromodynamics (QCD) [1]. The electromagnetic form factors of the nucleon describe the structure of the nucleon as seen by an electromagnetic probe. As such, they provide a window on strong interaction dynamics over a large range of momentum, for recent reviews see, e.g. Refs. [2,3]. Moreover, they are an important ingredient in the description of a wide range of observables ranging from the Lamb shift in atomic physics to the strangeness content of the nucleon [4,5]. At small momentum transfers, they are sensitive to the gross properties of the nucleon like the charge and magnetic moment as well as the radii. At large momentum transfer, in contrast, they probe the quark substructure of the nucleon as described by QCD.
A new twist was recently added to this picture by measurements of the proton charge radius in muonic hydrogen. The proton charge radius was first indirectly measured in the Nobel prize winning electron scattering experiments by Hofstadter et al. [6,7]. While electron scattering was the method of choice to refine the proton radius in the decades following these pioneering experiments, the Lamb shift in electronic hydrogen and muonic hydrogen is also sensitive to the proton radius [8]. The electronic Lamb shift measurements as well as most electron scattering experiments gave the so-called large radius, r P E 0.88 fm, which was also the value given by CODATA [9]. 1 It then came as a true surprise to most researchers when the first measurement of the muonic hydrogen Lamb shift, which has a larger sensitivity to r P E because of the much larger muon mass, led to the so-called small radius, r P E = 0.84184 (67) fm, differing by 5σ from the CODATA value [10]. At about the same time, a high-precision electronproton scattering experiment performed at the Mainz Microtron (MAMI) reinforced the large radius [11]. This glaring discrepancy in such a fundamental quantity, which was believed to be understood since long, became known as the "proton radius puzzle". It led to much experimental and theoretical activity dedicated to uncover its cause, either physics beyond the standard model, or more mundane reasons, such as an underestimation of the experimental uncertainties. Recent experiments on the electronic Lamb shift [12,13,14] and a novel measurement of electron-proton scattering at unprecedented small momentum transfer [15] now point to the latter reason. With the exception of Ref. [13], all of these new determinations of r p consistently give a small proton radius. Consequently, the newest addition of the CODATA compilation lists the proton charge radius as r p E = 0.8414 (19) fm [16]. A short review of the current status is given in Ref. [17]. The important role of dispersion theory in solving this "puzzle" will be discussed below.
This paper is structured as follows: In Sec. 2, we briefly review earlier dispersion-theoretical analyses of the electromagnetic nucleon from factors. The complete formalism to perform such type of analyses is given in Sec. 3, where all basic definitions are given and the various contributions to the spectral functions, the central objects of the dispersive method, are discussed in detail. Furthermore, constraints on the nucleon form factors and two-photon corrections to the electron-proton scattering cross section are presented. Finally, we display in detail methods to determine the theoretical uncertainties, both the statistical and the systematical ones. In Sec. 4, we display the results of our new dispersion-theoretical analysis of the electromagnetic form factors in the space-like region, including novel determinations of the various radii, form factors as well as the ω-meson couplings. Then, we consider the extension to the form factors in the time-like region and discuss the physics encoded in these. We end with a brief summary and an outlook in Sec. 5. In the appendices, we give further details on the extraction of neutron form factors from light nuclei as well as on the construction of the continuum contributions to the spectral functions. We also collect the various parameters of our best fit discussed in the main text.

Short history of dispersive analyses of the nucleon form factors
Here, we briefly review earlier work using dispersion theory to analyze the electromagnetic structure of the nucleon. To be more precise, we only consider investigations that include explicitly the two-pion continuum, which generates the ρ-meson in the isovector part of the spectral function in addition to a very important uncorrelated twopion contribution as first discussed by Frazer and Fulco [18,19,20]. For other work on dispersion relations applied to the nucleon electromagnetic from factors, we refer the reader to the review Ref. [21].
The first groundbreaking work was done by the Karlsruhe group in 1976 [22]. Here, electron-proton (ep) cross section data supplemented by neutron form factor data from elastic and quasi-elastic electron-deuteron scattering were fitted. Besides the two-pion continuum, the spectral functions contained the ω-meson plus additional isoscalar and isovector poles and normalization constants for the data sets. It should be noted that the ep data base was pruned in the sense that in case of inconsistencies between data sets, only one was retained. A dozen of fits with varying number of vector mesons poles and excluding various subsets of data were performed. The best fit (fit 8.2) featured 8 resonance parameters. Theoretical errors were estimated from the variations in the different fits. The resulting proton radii are tabulated in Tab. 1 and the neutron radii in Tab. 2. The neutron magnetic radius could not be determined very precisely at that time. Also notable were the sizable φN N couplings, where N denotes the nucleon, at odd with expectations from the OZI rule [23,24,25].
In 1995, the Bonn-Mainz group (MMD) rejuvenated the dispersion-theoretical approach to the nucleon electromagnetic form factors, as many new form factor results had become available and perturbative QCD had firmly established the behavior of the form factors at large momentum transfer [26]. Fits were performed to the existing form factor data basis of the Bochum group (updated from Refs. [27,28]). The two-pion continuum was still based on the Karlsruhe-Helsinki pion-nucleon (πN ) partial wave amplitudes f 1 ± (t), but the ρ-ω mixing visible in the pion vector form factor was included for the first time. The best fits where obtained with three additional isovector poles and and one additional isoscalar one (besides the ω and the φ). It was found that the onset of perturbative QCD was not seen in these data and the radii and vector meson couplings were consistent with the findings of the Karlsruhe group, see Tab. 1 and Tab. 2. Remarkably, these dispersive fits could not be made consistent with the then existing best value for r p E from ep scattering, r p E = (0.862 ± 0.012) fm [29]. Further, the large deviation in the OZI rule of the φ couplings was confirmed. One year later, the sparse and not very precise existing data on the proton and neutron form factors in the time-like region were included, which revealed some inconsistencies in the time-like data basis for the neutron [30].
In view of new data on the proton and neutron form factors, in particular the first polarization transfer measurements at Jefferson Lab at few GeV 2 squared momentum transfer [31,32], the MMD work was updated, with a particular emphasis on the magnetic radius of the proton and the neutron in [33]. In this work, no error analysis was performed.
A significant improvement of the dispersion relation (DR) analysis was performed in Ref. [34] (BHM). Not only was the data basis enlarged, but also the description of the isoscalar spectral function was improved by including the KK [35,36] and the πρ [37] continua. Furthermore, the 2π continuum was updated in view of new data for the pion vector form factor [38]. All data from the space-like and the time-like regions were included in the fit. In these fits, besides the mentioned continua, the isoscalar spectral function featured the ω, the φ and two poles, where as the 2π continuum was supplemented by 5 effective poles. The uncertainties were calculated by large scale Monte Carlo samplings of all solution with a χ 2 /dof in the range [χ 2 min , χ 2 min + 1.04], where χ 2 min refers to the best fit, corresponding to the 1σ coincidence in the p-value. Different to all earlier fits, the neutron charge radius squared was not included as a constraint. Nonetheless, the extracted value of (r n E ) 2 came out consistent with determinations from low energy atom-neutron scattering, see Tab. 2
that paper, it was stated that the then accepted proton charge radius determined from the Lamb shift in electronic hydrogen, r p E = 0.88 . . . 0.90 fm, see Ref. [39] (and references therein), was inconsistent with the dispersion analysis of the electron scattering data, thus previewing what was later called the "proton radius puzzle". The various radii came out consistent with earlier DR determinations, see Tab. 1 and Tab. 2. The same spectral functions were also used to extract the strength of two-photon corrections from the difference of data obtained by Rosenbluth separation and direct polarization transfer measurements [40]. The so determined two-photon corrections came out comparable to direct calculations available in the literature, such as Refs. [41,42,43].
The high-precision data with Q 2 ≤ 1 GeV 2 that emerged from MAMI-C in 2010 [11,44] called for a further update of the DR analysis. A first DR analysis in Ref. [45] utilized the same continua as BHM with the ω, the φ and three/five effective isoscalar/isovector poles. The fit was done to the reconstructed MAMI cross section data in the one-photon approximation and simultaneously to the neutron form factor data. The uncertainties were obtained varying the continua within reasonable ranges, namely the 2π continuum by 5% and the KK and πρ continua by 20%. Again, a small proton charge radius, r p E = 0.84 fm, emerged and the other radii were also agreeing with early DR determinations, very different to the values quoted in [11].
This work was further improved in various aspects in Ref. [46]. Here, only proton data were investigated, but two-photon corrections to the cross section were calculated and systematically included to the MAMI-C data, overcoming some inconsistencies in older approaches to this problem. Furthermore, to extract the statistical error due to the fitting procedure, a bootstrap approach was implemented. The spectral function was the same as in [45], but in addition, normalization constants for the various data sets were included (in total 31 new parameters) and the χ 2 definition was augmented by the correlation matrix. This method constituted an improvement about earlier error determinations. The uncertainties in the radii were somewhat increased compared to earlier determinations, see Tab. 1 and Tab. 2. The measured proton form factor ratio data for Q 2 < 1 GeV 2 [47,48] were not included in the fits but well described.
The work of Ref. [46] was extended by including neutron space-like form factor data as well as then existing data for the proton and the neutron in the time-like region in [49]. The emphasis of this work was to understand the structures seen by the BaBar collaboration [50] in the region between threshold up to highest measured momentum transfers. These structures (and similar but less pronounced ones in e + e − → nn) could be explained by including a φ(2170) meson as well as the N∆ and ∆∆ thresholds.
A significant improvement of the isovector spectral functions was reported in Ref. [51], based on the highprecision analysis of pion-nucleon scattering in the framework of the so-called Roy-Steiner equations [52]. This work also featured a detailed investigations of the corresponding isospin breaking effects in the pion form factor and the pion-nucleon P-wave amplitudes. The spectral functions given there serve as input for any DR analysis.
The most recent DR analysis in [53] was triggered by the PRad data [15], that measured ep cross sections at extreme forward angles corresponding to unprecedented small momentum transfers. In Ref. [53], fits to the PRad as well as the PRad and MAMI-C data were performed. The best fit to the combined data featured 5 isoscalar and 5 isovector poles, while the PRad data could be well described with 2+2 poles only. Again, the low-Q 2 data for µ p G p E /G p M were not included in the fit but could be well described. The error analysis was improved compared to earlier DR work, the bootstrap method was used to determine the "statistical" error, while the "systematic" error was obtained from varying the number of effective poles, e.g. in the combined analysis the range from (2+2) to (7+7) isoscalar + isovector poles was covered. This led to the very precise proton radii given in Tab. 1. It was pointed out that the statistical error in the PRad analysis is underestimated, consistent with the earlier findings of Ref. [54]. In Ref. [53], no uncertainties on the proton form factors were given and no neutron data were analyzed. In this review, we will fill this gap and present detailed results on these topics. Also, a Bayesian approach to calculate the statistical errors will be presented and compared to the bootstrap method.
It is remarkable how little fluctuations in the extracted values of the nucleon em radii based on dispersion relations have appeared with time, despite a dramatic improvement in the data base and a number of theoretical improvements, related in particular to the isoscalar and the isovector spectral functions and calculation of the theoretical uncertainties. There has also been some related work in the so-called dispersively improved chiral perturbation theory, see [55,56,57,58]. The extracted proton charge radius is consistent with our result, but as noted in Ref. [59], this approach is subject to uncertainties in the ρ-region, different from the exact representation used in the papers discussed above.
We end this section by noting that the so-called strangeness form factors of the nucleon can also be calculated (under certain assumptions) using the DR results for the isoscalar vector mesons, see e.g. Refs. [60,61,62,36]. For more details on this interesting topic, see the reviews [4,5].

Definitions
The electromagnetic (em) structure of the nucleon is determined by the matrix element of the vector current operator for the light quarks q = (u, d, s) T with the charges Q = diag(2, −1, −1)/3 (in terms of the elementary charge), sandwiched between nucleon states as depicted in Fig. 1. Denoting a nucleon state with four-momentum p as |p (for ease of notation, we do not display the corresponding spin or helicity index), with the help of Lorentz and gauge invariance and assuming CP invariance, this matrix element can be expressed in terms of two form factors, where m is the nucleon mass (which can be either the neutron, the proton or the isospin averaged mass) and t = (p − p) 2 the four-momentum transfer squared. For the analysis of data in the space-like region, it is convenient to use the variable Q 2 = −t > 0. The scalar functions F 1 (t) and F 2 (t) are the Dirac and Pauli form factors, respectively. They are normalized at t = 0 as with κ p = 1.793 and κ n = −1.913 the anomalous magnetic moment of the proton and the neutron, respectively, in units of the nuclear magneton, µ N = e/(2m p ). The magnetic moment of the proton and the neutron is thus given by µ p = 1 + κ p and µ n = κ n , respectively. For the theoretical analysis, it is often convenient to work in the isospin basis and to decompose the form factors into isoscalar (s) and isovector (v) parts, where i = 1, 2 . The experimental data are usually given in terms of the Sachs form factors where τ = −t/(4m 2 ). In the Breit frame, G E and G M may be interpreted as the Fourier transforms of the charge and magnetization distributions, respectively. The nucleon radii are defined via the low-t expansion of the form factors, where F (t) is a generic form factor. In the case of the electric and Dirac form factors of the neutron, G n E and F n 1 , the expansion starts with the term linear in t and the normalization factor F (0) is dropped. Note that the slopes of G n E and F n 1 are related via dG n E (Q 2 ) dQ 2 with m n the neutron mass. We remark that alternative information on the proton charge radius can be obtained from Lamb shift measurements in electronic as well as muonic hydrogen, see e.g. the reviews [63,64].
In the space-like region with t < 0, the form factors are real valued quantities. This is different in the timelike region with t > 0. By their very definition, at the nucleon-antinucleon threshold, t thr = 4m 2 , they fulfill the relation for both the proton and the neutron. In the physical region t > 4m 2 , the FFs are complex valued quantities. In Fig. 2, we sketch an exemplary form factor (here: G p M (t)) for all values of t. More precisely, the modulus of the form factor is depicted. For the space-like region, the threshold is located at t = 0, whereas the corresponding threshold in the time-like region is t = 4m 2 . In between these two thresholds, the various vector meson poles (plus continua) build up the spectral function to be discussed in detail below. This region cannot be observed. We note that for the form factors in the time-like region, an additional complication arises due to the strong near-threshold nucleon-antinucleon interactions, which will be considered in Sec. 4.4.  [30]. The colored area between the two dashed lines at t = 0 and t = 4m 2 is the unphysical region where the form factor cannot be observed.

Elementary cross section and polarization transfer
The form factors (FFs) can not be measured directly but are encoded in observables related to electron scattering. Consider for definiteness electron-proton (ep) scattering, where the four-momenta p i are subject to the constraint p 1 +p 2 = p 3 +p 4 . At first order in the electromagnetic finestructure constant α, the Born-approximation, the differential cross section can be expressed through the Sachs FFs as is the virtual photon polarization, θ is the electron scattering angle in the laboratory frame, and (dσ/dΩ) Mott is the Mott cross section, which corresponds to scattering off a point-like particle, where E 1 (E 3 ) is the energy of the incoming (outgoing) electron. Two quantities out of the energies, momenta and angles suffice to determine this cross section and are related for such an elastic process. Specifically, in the laboratory frame with the initial nucleon at rest and neglecting the electron mass, we can write In experiment, the differential cross section is usually given for a fixed total energy as a function of the scattering angle, so that a small scattering angle corresponds to a small momentum transfer. This is exactly the reason why a precise determination of the em radii is so difficult. At large momentum transfer, the contribution from the magnetic FF dominates the cross section. The contribution from the electric and the magnetic form factor can be read off form the reduced cross section σ R defined in Eq. (11). The reduced cross section σ R depends linearly on for a given Q 2 , with slope G E (Q 2 ) and intercept τ G 2 M (Q 2 ). This is called the Rosenbluth separation [65]. Two-photon corrections to this cross section will be discussed in Sect. 3.11. Also, to investigate the neutron FFs, one measures electron scattering of a light nucleus like deuterium or 3 He. This requires, however, some accurate few-body technique to disentangle the neutron contribution from the scattering cross section, as discussed briefly in App. A.
In early ep scattering experiments, it was found that the form factors could be well approximated by the dipole form, G dip (Q 2 ), with G n E (Q 2 ) = 0 in this approximation. Employing these dipole FFs in the integrated cross section Eq. (11) defines the so-called dipole cross section, σ dip . Often, the form factors or the measured cross sections are given relative to G dip (Q 2 ) and σ dip , respectively.
A method to directly measure the form factor ratio G E /G M in polarized electron scattering off the proton, − → e p → − → e p (or similarly for scattering off the deuteron or 3 He), has been proposed in Refs. [66,67]. A simultaneous measurement of the two recoil polarizations (longitudinal, P l , and transverse, P t ) allows one to measure directly the ratio While this only determines the form factor ratio (and not the individual FFs), many systematic uncertainties cancel out and make this observable an important benchmark for any theoretical form factor calculation. Let us briefly discuss the determination of the form factors in the time-like region. They can be extracted from the cross section data e + e − ↔pp and e + e − →nn for the proton and the neutron, respectively. As only very few differential cross section data exist in the time-like region, a separation of G E and G M is often not possible and one either makes an assumption like e.g. G E = G M in the analysis of the data or one extracts the effective form factor |G eff |, discussed below. For a review on the nucleon em form factors in the time-like region, see Ref. [2].
We now consider the process, in more detail. It is convenient to choose the center-ofmass (CM) frame, i.e., p 1,2 = (E, ±k e ) and p 3,4 = (E, ±k p ).
The photon momentum q then determines the center-ofmass energy by q 2 = (p 1 + p 2 ) 2 = t = E 2 CM = (2E) 2 . In the metric used here, time-like q implies positive q 2 . The three-momenta k e , k p appear in the phase-space factor β = k p /k e , which in the limit of neglecting the electron mass yields the velocity of the proton, and m p is the proton mass. We denote the emission angle of the proton by θ. The differential cross section in the one-photon-exchange approximation then is in terms of the electric and magnetic Sachs form factors and C(q 2 ) is the Sommerfeld-Gamow factor that accounts for the Coulomb interaction between the final-state particles Integrating over the full angular distribution gives the total cross section This defines the effective form factor G eff For neutrons, the formulas are equivalent except for the Sommerfeld-Gamow factor which is not present in that case. Beyond the Coulomb final-state interactions, higher order QED corrections are usually neglected. For the timereversed process, the phase space factor is inverted, yielding Taking into account the angular dependence of pp production, one can express the differential cross section via the angular asymmetry A, see Ref. [68], This can be determined from the FF ratio R = |G E /G M |. The dispersion integral is calculated for the space-like value of t also shown by a cross.

Dispersion relations and spectral decomposition
Dispersion relations (DRs) are based on unitarity and analyticity. Here, DRs relate the real and imaginary parts of the electromagnetic nucleon form factors. Let F (t) be a generic symbol for any one of the four independent nucleon form factors. We write down an unsubtracted dispersion relation of the form where t 0 is the threshold of the lowest cut of F (t) (see below) and the i defines the integral for values of t on the cut. The convergence of an unsubtracted dispersion relation for the form factors has been assumed. For proofs of such a representation in perturbation theory, see Ref. [69] (and references therein). One could also use a once subtracted dispersion relation, since the normalization of the form factors at t = 0 is known. However, in what follows, we will only employ the unsubtracted form give in Eq. (26). Most importantly, by Eq. (26) the electromagnetic structure of the nucleon can be related to its absorptive behavior. In Fig. 3 we display the analytic structure underlying the dispersion integral in Eq. (26). The various ingredients (continuum cuts, vector meson poles) will be discussed in detail below. The imaginary part Im F entering Eq. (26) can be obtained from a spectral decomposition [70,71]. For this purpose, consider the electromagnetic current matrix element in the time-like region (t > 0), which is related to the space-like region (t < 0) via crossing symmetry. This matrix element is given by where p andp are the momenta of the nucleon and antinucleon created by the current j em µ , respectively. The four-momentum transfer squared in the time-like region is t = (p 3 + p 4 ) 2 .
Using the LSZ reduction formalism, the imaginary part of the form factors is obtained by inserting a complete set of intermediate states as [70,71] (28) where N is a nucleon spinor normalization factor, Z is the nucleon wave function renormalization, andJ N (x) = J † (x)γ 0 with J N (x) a nucleon source. This decomposition is illustrated in Fig. 4. It relates the spectral function to on-shell matrix elements of other processes, as detailed below.
The states |λ are asymptotic (observable) states of momentum p λ . They carry the same quantum numbers as the current j em µ : for the isovector component (29) of the current j em µ . Here, I and J denote the isospin I = 0, 1 and the angular momentum J = 1 of the photon, whereas G, P and C give the G-parity, parity and charge conjugation quantum number, respectively. Furthermore, these currents have zero net baryon number. Because of Gparity, states with an odd number of pions only contribute to the isoscalar part, while states with an even number contribute to the isovector part. For the isoscalar part the lowest mass states are: 3π, 5π, . . . , KK, KKπ, . . . , (30) and for the isovector part they are: 2π, 4π, . . . .
Associated with each intermediate state is a cut starting at the corresponding threshold in t and running to infinity. As a consequence, the spectral function Im F (t) is different from zero along the cut from t 0 to ∞, with t 0 = 4 (9) M 2 π for the isovector (isoscalar) case. The spectral functions are the central quantities in the dispersion-theoretical approach. Using Eqs. (27,28), they A B Fig. 5. Two-pion cut contribution to the isovector form factors, given in terms of the pion vector form factor F V π (represented by A) and the ππ →N N P-waves f 1 ± (represented by B). Solid, dashed, and wiggly lines denote nucleons, pions, and the external photon, respectively, while the dash-dotted line indicates the cutting of particle propagators.
can in principle be constructed from experimental data. In practice, this program can only be carried out for the lightest two-particle intermediate states.
The longest-range, and therefore at low momentum transfer most important continuum contribution comes from the 2π intermediate state which contributes to the isovector form factors [72]. A novel and very precise calculation of this contribution has recently been performed in Ref. [51] including the state-of-the-art pion-nucleon scattering amplitudes from dispersion theory, as detailed below. In the isoscalar channel, the inclusion of the KK [35,36] and ρπ continua [37] was first introduced in Ref. [34] in the dispersive analysis of the em form factors. These important ingredients are discussed in more detail below. Apart from the continua, there are also single vectormeson pole contributions. As will become clear in the following, the contributions from the continua and the poles are sometimes strongly intertwined, e.g. the ρ-meson pole is indeed generated as part of the 2π-continuum, as known since long [18,19,20].

Two-pion continuum
In the isovector channel, the lowest continuum contribution is given by the two-pion exchange as depicted in Fig. 5. Therefore, the unitarity relations for the nucleon form factors reads [20] is the vector form factor of the pion and the f 1 ± (t) are the P-wave πN partial waves in the t-channel. Watson's theorem ensures that the lefthand side of the equations stays real, as long as the same ππ phase shift is used in the calculation of the pion form factor and the ππ →N N partial waves. Therefore, in the most recent determination of the two-pion continuum, the same three variants of the phase shift δ 1 1 in the data fits for F V π (t) were used as in the Roy-Steiner analysis of pionnucleon scattering [52]. The full consistency among all ingredients entering the unitarity relation that was achieved in Ref. [51] was a key improvement over earlier calculations, and thus the representation of the two-pion continuum given there will be discussed in what follows.
It is important to discuss the range of validity of the 2π approximation to the unitarity relation. Strictly speaking, the 4π threshold opens at √ t = 4M π = 0.56 GeV, but it is well known from phenomenology that the 4π contribution is completely negligible below the ωπ threshold at √ t = 0.92 GeV, see e.g. [73], and only becomes sizable once the ρ , ρ resonances are excited. This can also be understood from chiral perturbation theory, where the 4π contribution appears first at three loop order [74]. For this reason, the two-pion cut contribution to the isovector form factors is Let us now discuss in more detail the various ingredients entering the Eqs. (32). We start with the vector (em) form factor of the pion. It is given by The ππ intermediate states produce the unitarity relation with the ππ P-wave phase shift δ 1 1 . Eq. (34) reflects Watson's final-state theorem [75], which states that the phase of F V π has to coincide with the ππ scattering phase shift (up to multiple integers of π). Neglecting higher intermediate states unitarity determines F V π (t) up to a polynomial P (t) in terms of the Omnès factor Ω 1 .
(35) In fact, the representation (35) provides a very efficient and accurate parameterization of the experimental data, up to the distortions due to ρ-ω mixing. This isospinviolating effect can be included via a modification of F V π (t), with the ω mass M ω and width Γ ω . The parameters α and , where parameterizes the strength of the ω-ρ mixing, are fit to recent form factor data, see Refs. [77,78,79] below √ t = 1 GeV, using the same ππ phase shifts as in the RS analysis [52]. The latter has been determined from Roy and Roy-like equations by the Bern [80] and the Madrid-Cracow group [81]. To get a better handle on the uncertainty estimate for the final spectral functions from the pion vector FF, in Ref. [51] a variant of the Bernese phase shift was also considered. It includes effects from the ρ and the ρ in an elastic approximation [82].
Next, we discuss the the t-channel partial waves f 1 ± (t), given by [19] f J with the t-channel scattering angle z t = (s − u)/(4p t q t ), the P J are the Legendre polynomials, and the momenta are q t and p t = t/4 − m 2 . Further, the standard decomposition of the πN scattering amplitude T (π a (q)+N (p) → π b (q ) + N (p )) in the isospin limit has been used, where a, b are isospin indices, the τ a are isospin Pauli matrices, I = ± refers to isoscalar/isovector amplitudes and s = (p+q) 2 , t = (p −p) 2 , u = (p−q ) 2 are the Mandelstam variables subject to the constraint s+t+u = 2(m 2 +M 2 π ). The best way to determine the pion-nucleon scattering amplitudes are undoubtedly dispersion relations, as they allow for a systematic continuation from the physical region into the unphysical ones and further make best use of the existing scattering data. The most modern and accurate investigations are based on the Roy-Steiner (RS) equation analysis of the Bonn group, developed and performed in Refs. [83,84,85,52,86] (for earlier work by the Karlsruhe-Helsinki group, see e.g. Refs. [87,88]). The RS equations, originally developed in [89,90] (and references therein), are hyperbolic DR of the form which have a number of advantages compared to other formulations (like e.g. fixed-t DR). They combine all physical regions, display an explicit s ↔ u crossing, require the absorptive parts only in regions where the corresponding partial wave expansions converge, and, further, a judicious choice of the parameter a allows to increase the range of convergence. The RS equations have a limited range of validity, where √ s m , √ t m denotes the so-called matching point for the s-and t-channel partial waves, respectively. The required inputs to solve the RS equations are the S-and Pwaves above the matching point, the higher partial waves (D-, F-, . . .) and the inelasticities. An important constraint are the pion-nucleon scattering lengths deduced from pionic hydrogen and pionic deuterium [91] (for a recent update, see [92]). The output of the RS equations are the so-called subthreshold parameters, which allow one to reconstruct the scattering amplitude in the unphysical region, such as the f 1 ± (t) in the pseudophysical region required for the isovector spectral functions. Some basic definitions of the πN scattering amplitude in the unphysical region are given in App. B. It should also be mentioned that the results of the RS analysis were given with theoretical uncertainties, which to our knowledge has been the first time that a dispersive analysis of pion-nucleon scattering provided these, for details see [52].
In Ref. [51], the isospin-violating effects beyond the ρω mixing in the pion form factor where also worked out, leading to an improved representation of the unitarity relations, Eq. (32), namely For a more detailed discussion of this representation, the reader is referred to Ref. [51].
Putting pieces together, the isovector spectral functions divided by t 2 based on Eq. (40) are shown in Fig. 6. These nicely exhibit the ρ-resonance at √ t = 0.77 GeV as well as a remarkable enhancement on the left shoulder of the resonance. This shows that the ρ is indeed generated by unitarity [18] and thus no explicit ρ-meson is required in the isovector spectral function. The visible enhancement on the left shoulder of the ρ can be traced back to the fact that the partial wave amplitudes f 1 ± (t) have a singularity on the second Riemann sheet [88] (originating from the projection of the nucleon pole terms in the invariant pion-nucleon scattering amplitudes) located at very close to the physical threshold at t 0 = 4M 2 π . The isovector form factors inherit this singularity (on the second sheet) and the closeness to the physical threshold leads to the pronounced enhancement between √ t = 0.3− 0.6 GeV shown in Fig. 6. This issue will be taken up below.
isospin limit with ρ-ω mixing isospin limit with ρ-ω mixing The uncertainties displayed in Fig. 6 originate from three different sources: 1) the subthreshold parameters b − 00 , b − 01 , a − 00 and a − 01 (as defined in App. B), 2) the pion-pion phase shift δ 1 1 (t) and 3) the data for the pion form factor F V π (t). In fact, the uncertainty of the subthreshold parameters from the RS analysis is in fact the dominating effect below 1 GeV. We note that the effect of the ρ-ω mixing is small, as the comparison of the black and red dashed lines in Fig. 6 shows. Note also that this consistent inclusion of isospin-breaking effects in the pion em form factor and the πN partial waves constitutes a major achievement compared to earlier analyses. The two-pion continuum contribution to the isovector form factors is displayed below in Fig. 8 .
Based on the DR, Eq. (26), it is is straightforward to derive sum rules for the normalizations and radii of the isovector form factors. These were first considered in Ref. [72] for the various nucleon radii, see also [51], where 353 is the isovector magnetic moment of the nucleon. Note that the sum rules for the radii remain unchanged if a once-subtracted dispersion relation is used instead of the unsubtracted one. Cutting the integrals at Λ = 2m, one finds It is remarkable that just using a simple ρ-exchange using e.g. a Breit-Wigner or a Gounaris-Sakurai form [93], the corresponding isovector radii would be sizeably underestimated (by about 40%), as inspection of Fig. 6 reveals. Thus, any dispersive analysis that does not include the full two-pion continuum but only the ρ-resonance in the isovector spectral function below 1 GeV will simply miss important physics. We will come back later to these sum rules.

Three-pion continuum
The lowest isoscalar continuum is given by 3-pion exchange as depicted in Fig. 7. There, A refers to the γ → 3π transition amplitude, that is given at low energies by the anomalous Wess-Zumino-Witten Lagrangian [94,95] and B corresponds to the 3π →N N amplitude. An analysis based on unitarity alone of this contribution does not exist, but it has been shown in chiral perturbation theory at leading [96] and subleading [97] orders, that there is no enhancement on the left wing of the ω resonance. This argument is outlined in App. C. Thus, the usual inclusion of the ω as a vector meson pole is justified. In case of the φ, the situation is, however, more complicated as discussed next.

KK continuum
The first important continuum contribution to the isoscalar spectral function is the one from KK states, as evaluated in Refs. [35,36] from an analytic continuation of KN scattering data. The KK contribution to the imaginary part A B Fig. 7. Three-pion cut contribution to the isoscalar form factors, given in terms of the γ → 3π (represented by A) and the 3π →N N (represented by B) transition amplitudes. Solid, dashed, and wiggly lines denote nucleons, pions, and the external photon, respectively, while the dash-dotted line indicates the cutting of particle propagators.
of the isocalar form factors is given by [35,36] Im F where p t = t/4 − m 2 and q t = t/4 − M 2 K , with M K the charged kaon mass. Further, F V K (t) is the kaon form factor, defined via whereas the b 1/2, ±1/2 1 (t) are the J = 1 partial wave amplitudes for KK → NN [35,36]. Having determined these imaginary parts, the contribution of the KK-continuum to the form factors is obtained from the dispersion relation Eq. (26).
The b (t) in the above equations are the kaon-nucleon partial wave amplitudes with total angular momentum J = 1 (for definitions, see App. D). For t ≥ 4m 2 the partial waves are bounded by unitarity, In the unphysical region 4M 2 K ≤ t ≤ 4m 2 , however, they are not constrained by unitarity. In Ref. [35], the amplitudes b calculated. Strictly speaking this calculation provides an upper bound on the spectral function since one replaces the amplitudes and the form factor in Eqs. (44,45) by their absolute values. The striking feature in the spectral function is a clear φ resonance structure just above the KK threshold. The resonance emerges in the partial wave amplitude b 1/2, 1/2 1 as well as in the kaon form factor F K . In contrast to the 2π continuum, there is no strong enhancement on the left wing of the φ resonance which sits directly at the KK threshold. This is completely analogous to the lowest-lying isoscalar resonance, the ω(782), which also does not exhibit any enhancement on its left shoulder.
The resulting contribution to the nucleon form factors can be parameterized by a pole term at the φ mass [34]: with a KK 1 = 0.1054 GeV 2 and a KK 2 = 0.2284 GeV 2 . As a consequence, the contribution of the KK continuum to the electromagnetic nucleon form factors can conveniently be included in the analysis via Eq. (48). The form factor contributions from Eq. (48) are also shown in Fig. 8.

ρπ continuum
Another important contribution to the isoscalar spectral function is the correlated ρπ exchange, that was investigated in the Bonn-Jülich nucleon-nucleon interaction model in Ref. [98]. Since in that work cancellations between φ-exchange contribution and this correlated πρ-exchange was found, the ρπ contribution to the isoscalar spectral function was worked out in Ref. [37]. This continuum contribution was evaluated in terms of a dispersion integral which in turn can be represented by an effective pole term γ V N N Fig. 9. Vector meson dominance: The photon only couples through vector mesons, V = ρ, ω, φ, . . ., to the nucleon.
for a fictitious ω meson with a mass M ω = 1.12 GeV [37]: with a ρπ 1 = −1.01 GeV 2 and a ρπ 2 = −0.04 GeV 2 . In the form factor analysis, one uses this effective pole instead of the full spectral function.
There is very little sensitivity in the dispersive fits to a ρπ 2 , which can vary between −0.04 and −0.4 without affecting the outcome of the fit. If the ω pole is treated as a real resonance, the latter value is consistent with f ω ∼ 10 for a ρπ 1 = −1.01 if the coupling constants g i ω N N (i = 1, 2) from Ref. [37] are used as input (for a precise definition of these couplings, see Sec. 3.8).
In Fig. 8, we show the contribution of the 2π, KK, and ρπ continua to the electromagnetic nucleon form factors F 1 and F 2 . The 2π contributes to the isovector form factors while the KK and ρπ continua contribute to the isoscalar form factors. The KK and ρπ contributions have opposite sign and partially cancel each other. The dominant contribution to F s 1 comes from the ρπ continuum while for F s 2 the KK contribution is larger. While the KK and ρπ contributions can be represented by simple pole terms, the expressions for the 2π continuum Eq. (40) are more complicated. This is related to the strong enhancement close to the 2π threshold on the left wing of the ρ resonance discussed above. Finally, note that these continuum contributions enter as an independent input in the dispersive analysis. They are not fitted to cross section or form factor data.

Vector meson poles
In the most simple picture, the photon couples to the nucleon through vector mesons only (i.e. there is no direct photon-nucleon coupling), the so-called vector meson dominance (VMD) picture, see e.g. [99,100,101,102], as depicted in Fig. 9. In this picture, the form factors take the simple form with and the couplings f V can be deduced from the leptonic decay widths V → e + e − , Also, we have identified s 1 , s 2 with the ω, φ and v 1 with the ρ. Each such vector meson comes with two couplings, the vector coupling a V 1 and the tensor coupling a V 2 . One also employs the ratio of the tensor to the vector coupling, κ V , defined via While for some resonances these couplings can be deduced from nucleon-nucleon scattering data, in the dispersive analysis, they are considered as fit parameters (with the exception of the ρ, which is completely determined from the 2π continuum). In the pure VMD picture with only ρ and ω vector mesons contributing, one can relate the tensor-to-vector coupling ratio to the isovector and isoscalar anomalous magnetic moments of the nucleon, such that However, extracting κ ρ from the two-pion continuum leads to a larger value, κ ρ 6, consistent with extractions from nucleon-nucleon scattering, see e.g. the discussion in [26]. The corresponding imaginary part, i.e. the contribution to the spectral function for any vector meson reads: As already discussed in Sec. 3.4, the ρ-meson is entirely generated by the two-pion continuum, so that an explicit ρ will never appear in the spectral function. Different to that, the lowest isoscalar mesons are the ω and the φ, which are explicitly taken into account. As noted before, the related 3π continuum has a very small nonresonant contribution, that can be safely neglected [96], see also App. C. Also, in the isoscalar region around 1 GeV, we consider the KK and ρπ continua, which tend to cancel, and an additional residual φ pole. Because of the complicated structure of the isoscalar spectral function around 1 GeV, it is no longer possible to extract useful φN N couplings, as it was done in earlier works, where one just had the φ-pole in this region. The large φ-couplings found in these earlier studies are clearly an artifact of the simplified isoscalar spectral function assumed in this region.

Structure of the spectral functions
As discussed above, the spectral function can at present only be obtained from unitarity arguments and experimental data for the lightest two-particle intermediate states (2π and KK). Furthermore, the ρπ continuum contribution has been calculated in the Bonn-Jülich N N model.
The remaining contributions to the spectral function can be parameterized by vector meson poles. On the one hand, the lower mass poles can be identified with physical vector mesons such as the ω and the φ. The higher mass poles on the other hand, are simply an effective way to parameterize higher mass strength in the spectral function. These effective poles at higher momentum transfers appear in the isoscalar (s 1 , s 2 , ...) and isovector channels (v 1 , v 2 , ...) It should be noted that we are dealing with an ill-posed problem here [103,104], that means increasing the number of poles will from some point on not improve the description of the data. Therefore, the strategy has always been to use as few poles as possible. We come back to this issue in Sec. 3.12.
Putting all pieces together, the spectral function has the general structure For the light isoscalar vector mesons, the residua in the pole terms can be related to their couplings. Only rough estimates exist for these: [37]. Note that the dominant vector ωN N coupling is taken to be positive, consistent with one-boson-exchange in the nucleon-nucleon interaction. These ranges are used as constraints in the fits. The masses of the effective poles (s 1 , s 2 ,-. . . , v 1 , v 2 , . . .) are fitted to the data. We remark that to ensure the stability of the fit [104], we demand that the residua of the vector meson poles are bounded, |a V i | < 5 GeV 2 (this can also be considered a naturalness argument for the couplings), and that no effective poles with masses below 1 GeV appear. Furthermore, the masses of these effective poles should also be smaller than 5 GeV. We generally do not include widths for the effective poles. However, if one wants to mimic the imaginary part of the form factors in the time-like region, one can e.g. allow for a large width for the highest mass effective pole (see, e.g., Ref. [34]). A cartoon of the resulting (isoscalar and isovector) spectral functions is shown in Fig. 10. The vertical dashed line separates the phenomonologically wellconstrained low-mass region from the effective vector meson poles at higher masses.

Constraints
The number of parameters in the spectral function (i.e. the various meson couplings a V i (i = 1, 2) and the masses of the effective poles) is reduced by enforcing various constraints.
The first set of constraints concerns the low-t behavior of the form factors. We enforce the correct normal-  ization of the form factors as given in Eq. (3). The nucleon radii, however, are not included as a constraint. The exception to this is the squared neutron charge radius, which in some dispersive fits has been constrained to the value from low-energy neutron-atom scattering experiments [106,107]. In the new fits discussed later, we implement this constraint using the high-precision determination of the neutron charge radius squared based on a chiral effective field theory analysis of electron-deuteron scattering [108,109], Another set of constraints arises at large momentum transfers. Perturbative QCD (pQCD) constrains the behavior of the nucleon electromagnetic form factors for large momentum transfer. Brodsky and Lepage [110] worked out the behavior for Q 2 → ∞, in terms of the leading order QCD β-function. The anomalous dimension γ ≈ 2 depends weakly on the number of flavors, N f [110]. The power behavior of the form factors at large Q 2 can be easily understood from perturbative gluon exchange. In order to distribute the momentum transfer from the virtual photon to all three quarks in the nucleon, at least two massless gluons have to be exchanged. Since each of the gluons has a propagator ∼ 1/Q 2 , the form factor has to fall off as 1/Q 4 . In the case of F 2 , there is additional suppression by 1/Q 2 since a quark spin has to be flipped. The analytic continuation of the logarithm in Eq. (59) to time-like momentum transfers −Q 2 ≡ t > 0 yields an additional term, ln(−t/Λ 2 ) = ln(t/Λ 2 ) − iπ for t > Λ 2 . Employing the Phragmen-Lindeloef theorem [88], it follows that the imaginary part has to vanish in the asymptotic limit. Taking these facts into account, the proton effective FF can be described for large time-like mo-mentum transfer t by [111] |G p with the parameters from a fit to data prior to the 2013 measurement by the BaBar collaboration [50], given as A = 72 GeV −4 and Λ = 0.52 GeV. The power behavior of the form factors leads to superconvergence relations of the form with n = 0 for F 1 and n = 0, 1 for F 2 . These will be employed in the current analysis. In earlier DR analyses, modifications of the superconvergence relations were used including e.g. some higher order corrections. These should be, however, abandoned as the data are simply not sensitive to such corrections. We note that these superconvergence relations have already been used in Ref. [22], i.e. before the pQCD analysis. Consequently, the number of effective poles in Eqs. (56,57) is determined by the stability criterion mentioned before, that is, we take the minimum number of poles necessary to fit the data. The number of free parameters is then strongly reduced by the various constraints (unitarity, normalizations, superconvergence relations). These constraints can be implemented as what is called "hard constraints" or "soft constraints", respectively. In the former case, one solves a system of algebraic equations relating the various parameters (couplings, masses), thus reducing the number of free parameters in the fit (for an explicit representation, see e.g. [26]). In the latter case, the χ 2 is augmented by a Lagrange multiplier enforcing the corresponding constraints, see Sec. 3.12. Both options are viable and have been used.
It is straightforward to enumerate the nunber of fit parameters, which is given by the couplings and masses of the vector meson, N V = 4 + 3N s + 3N v , with N s (N v ) the number of the effective isoscalar (isovector) poles and the 4 represents the ω and φ couplings, minus the number of constraints, given by N C = 4 + 6 + 1, referring to the low-t, the high-t constraints and the neutron charge radius squared, respectively. If the latter in not included, N C = 10. Putting pieces together, we have in total

Two-photon effects
The interest in two-photon corrections was triggered by the high precision measurements of the form factor ratio G E /G M using the polarization transfer method reported in Refs. [31,32]. These results were found to be in striking disagreement with the world data based on the Rosenbluth separation. Even after removing inconsistencies from the Rosenbluth data base [112], this discrepancy remained, triggering a flurry of works on two-photon corrections beyond the work of Refs. [113,114,115], which neglected, however, the effects of the structure of the nucleon in the calculation of the two-photon box and crossed-box diagrams, see Fig. 11. These diagrams were calculated in various approaches, like in hadronic models, using generalized parton distributions or using dispersive methods, see e.g. Refs. [116,117] for reviews and very recent work in Ref. [118]. Here, we concentrate on the work presented in [46], because the two-photon corrections given there are applied in the DR analyses since then. The corrections to the electron-proton cross sections at order α 3 are given by the interference of the one-photonexchange amplitude M 1γ and the amplitudes from vacuum polarization, vertex corrections, self-energy corrections and the two-photon-exchange amplitude M 2γ and additionally the contribution from Bremsstrahlung. The main data set that are considered in the DR analyses already contains a set of calculations of such corrections by Maximon and Tjon [115]. This calculation contains improvements towards earlier works by Mo and Tsai [114] but still uses a soft-photon approximation, particularly relevant for the two-photon exchange (TPE) contribution. This contribution to the corrected cross section can be expressed through a factor of (1 + δ 2γ ) as so that We briefly discuss the soft-photon approximation by Maximon and Tjon since only the difference between any new evaluation of the 2γ corrections and this approximation is required for the purification of the ep scattering data.
Ref. [115] separates the IR-divergent part of the TPEamplitude by considering the poles in the photon propagators, i.e. one vanishing photon momentum. The resulting factor is where λ is an infinitesimal photon mass and E 1 (E 3 ) the incoming (outgoing) electron energy. The logarithmic infrared singularity in λ is canceled by a term in the Bremsstrahlung correction, so that the complete cross section is λ-independent. The same cancellation takes place, if both δ 2γ,IR and the Bremsstrahlung correction are calculated in the older approximation scheme by Mo and Tsai. In Ref. [46], the interference between the 1γ-amplitude and the 2γ-amplitude was calculated. In this notation, the metric tensor from the photon propagator has already been contracted. Then, M 1γ is given in terms of the conventional lepton spinors u(p) and the elastic nucleon-vertex Γ ν (q) from Eq. (2). The 2γ-amplitude contains the lepton tensor whereas the hadronic tensor for nucleon or ∆ intermediate states are respectively. Here, Γ µα γ∆→N (p, k) is the transition vertex in terms of the Raita-Schwinger spinor field Ψ (a) µ (p) [119]. Various parameterizations of this transition matrix element exist, see e.g. Refs. [120,121]. The corresponding electric G E , magnetic G M and Coulomb G C transition form factors can be related to the helicity amplitudes measured in pion electroproduction off the nucleon, see e.g. Ref. [122]. Accounting for this momentum dependence also in the ∆N γ vertices in the box and crossed box diagrams was the main improvement in Ref. [46] compared to some earlier calculations. Further, in the denominator of the photon propagator for the pure nucleon graph, one includes an infinitesimal photon mass λ to regulate the infrared divergences. The loop containing the ∆ is not IR divergent because of the mass of the ∆.
The S F and S αβ in Eqs. (68,69) are the conventional spin-1/2 and spin-3/2 propagators, respectively. The calculation of the crossed box graph proceeds accordingly. In Fig. 12, the -dependence at Q 2 = 3 GeV 2 is shown, which allows for a comparison to previous calculations Ref. [42] (red solid line). Displayed is the difference of the calculation in [46] to the soft-photon approximation by Maximon and Tjon [115]. such as [42]. In this case, the dependence on the nucleon FF parameterization largely cancels out. The use of the pole fit parameterization from Ref. [42] indeed reproduces their result. Lowering the Q 2 -value in the calculation decreases the nucleon-TPE correction. For the intermediate ∆, the situation is different, one finds a stronger dependence on the FF parameterizations. In Fig. 13 the results of the calculation that employs the helicity amplitudes obtained from data on electroproduction of nucleon resonances [123] are displayed. These corrections can be parameterized conveniently by a set of FFs, determined in Ref. [122] and used in Ref. [124] for a similar calculation albeit without realistic NFFs. This form of the γN ∆-vertex does not deviate significantly from recent data and is numerically well-treatable. These results are similar to the ones of Ref. [125], where different transition form factors are employed. As stated before, the sum of δ 2γ,N and δ 2γ,∆ from Ref. [46] constitute the two-photon corrections employed in the DR analysis of the Bonn-Darmstadt group. Their effect on the high-precision data from Mainz [11] is displayed in Fig. 14. Note that the original data contain an approximation of the two-photon correction given by [126].
where Z is the nuclear charge (here, Z = 1). As pointed out in Ref. [127], this approximation is only valid as Q 2 → 0 and has the wrong sign for some kinematical regions. Thus, this contribution is subtracted from the data and the two-photon corrections from Ref. [46] are added. The differences are quite visible. The corrections from Ref. [46] will be also employed in the calculations presented in the next section. Nevertheless, an updated calculation of these corrections would be welcome.

Fit strategies and error analysis
In this section, we briefly describe how the fits of the spectral functions to data are performed and how the statistical and systematic errors can be determined. First, we discuss the quality of the fits, which is measured in terms of the total (traditional) χ 2 , where C i are the cross section data at the points Q 2 i , θ i and C(Q 2 i , θ i , p ) are the cross sections for a given FF parameterization for the parameter values contained in p. Moreover, the n k are normalization coefficients for the various data sets (labeled by the integer k), while σ i and ν i are their statistical and systematical errors, respectively. A more refined definition of the χ 2 is given by [46] in terms of the covariance matrix V ij = σ i σ j δ ij + ν i ν j . This latter definition accounts for the correlation between the various fit parameters. A fit to form factor data uses the same definitions, except for the absence of the normalization factors. One also considers the reduced χ 2 , which is given by:  with N D the number of fitted data points and N F the number of independent fit parameters, see Sec. 3.10.
As noted in Sec. 3.10 the various constraints on the form factors can be implemented algebraically (hard constraints) or by modifying the χ 2 (soft constraints). The latter type of constraints are implemented as additive terms to the total χ 2 in the following form where x is the desired value and p is a strength parameter, which regulates the steepness of the exponential well and helps to stabilize the fits [128,34]. One method to estimate the fit (statistical) errors is the bootstrap procedure, see e.g. Ref. [129]. One simulates a large number of data sets compared to the number of data points by randomly varying the points in the original set within the given errors assuming their normal distribution. Let us consider the radius extraction. In that case, one fits to each of these data sets separately, extracts the radius from each fit and consider the distribution of these radius values, which is sometimes denoted as bootstrap distribution. The artificial data sets represent many real samples. Therefore, this radius distribution represents the probability distribution that one would get from fits to data from a high number of measurements. The precondition for using this method are independent and identically distributed data points which is fulfilled when the χ 2 sum does not depend on the sequential order of the contributing points. For n simulated data sets, the errors thus scale with 1/ √ n. However, to get a more realistic uncertainty, we exclude one percent of the data points from the sample and can so determine the lowest and highest value of the extracted radius. The same procedure can, of course, also be applied to the full form factors. In Fig. 15, we again use the 71 PRad data points to show the bootstrap extraction of the proton charge radius and its statistical uncertainty based on 1000 samples. The extracted error thus reads (a similar plot is obtained for the magnetic radius) [53] δ(r p E ) stat. = ±0.012 fm , δ(r p M ) stat. = ±0.005 fm . (77) We note that the bootstrap error for r p M for the PRad data given in [53] is corrected here.
Another statistical tool to estimate the error intervals of our model parameters is the Bayesian approach, see e.g. Ref. [130] (and references therein). In contrast to the interpretation of probabilities in the classical (also called frequentist) approach, where the probability is the frequency of an event to occur over a large number of repeated trials, the Bayesian method uses probabilities to express the current state of knowledge about the unknown parameters, which allows one to estimate the uncertainty as a statement about the parameters. The key ingredients to a Bayesian analysis are the prior distribution, which quantifies what is known about the model parameters prior to data being observed, and the likelihood function, which describes information about the parameters contained in the data. The prior distribution and likelihood can be combined to derive the posterior distribution by means of Bayes' theorem: where "paras" denotes the parameters. It is the main goal of a Bayesian statistical analysis to obtain the posterior distribution of the model parameters. The posterior distribution contains the total knowledge about the model parameters after the data have been observed. From a Bayesian perspective, any statistical inference of interest can be obtained through an appropriate analysis of the posterior distribution. For example, point estimates of parameters are commonly computed as the mean of the posterior distribution and interval estimates can be calculated by producing the end points of an interval that correspond with specified percentiles of the posterior distribution. A powerful and easyto-implement method to access posterior distribution is Markov Chain Monte Carlo (MCMC) algorithm. A systematic illustration of Bayesian analysis application can be found in Ref. [131].
As an example, we implement a Bayesian analysis for the fit to PRad data where the 2s+2v configuration of the spectral function is used. The likelihood function is given by with the χ 2 objective function defined in Eq. (74). Here, p contains the model parameters and D = {d i } denotes the PRad data points. N is the normalization constant. Two different prior distributions shown in Fig. 16 are considered to test the stability of the obtained statistical outputs from our Bayesian analysis. We apply a particular MCMC sampling algorithm called ParaMonte [132] to acquire a Monte Carlo sample from the posterior distribution. The obtained posteriors of the parameters m V 1 and a ω 1 are taken as an example to show the equivalence of normal and uniform priors we used as shown in Fig. 17. The statistical estimates of form factor and radius errors from our Bayesian analysis are discussed in next section.
Next, we discuss the extraction of the systematic uncertainties, which is always the most difficult task. Our strategy is similar to what was already done in Ref. [22], namely to vary the number of isoscalar and isovector poles around the values corresponding to the best solution, where the total χ 2 does not change by more than 1%. An example of this is given in Tab. 3 taken from Ref. [53]. Here, only the PRad data [15] are considered. The best fit corresponds to 2 isoscalar and 2 isovector poles, so we can read off the systematic errors in this case as [53] δ(r p E ) syst. = ±0.001 fm , δ(r p M ) = +0.018 −0.012 fm .
We note that while the absolute χ 2 does not change, the reduced one worsens as the number of fit parameter increases. As expected, the systematic error is larger for the magnetic radius as at low Q 2 , the electric FF dominates.
More detailed results will be given below.  Table 3. Fit to the PRad data with varying numbers of isoscalar (s) and isovector (v) effective poles. Given are the total and the reduced χ 2 and the resulting values for the proton radii. The * marks the best solution which defines the central values for the radii.

Physics results
In this section, we display a number of physics results, in particular we discuss fits including proton polarization transfer and neutron form factor data, and present new uncertainty analyses, thus extending and deepening the work of Ref. [53]. We also discuss the inclusion of data for the time-like form factors and the related physics. First, however, we want to sharpen and validate our toolbox to pin down the errors on the example of the PRad data.

Detailed analysis of the PRad data
The PRad data [15] are given at two beam energies, E = 1.1, 2.2 GeV, covering squared momentum transfers in the range Q 2 = 2 · 10 −4 − 6 · 10 −2 GeV 2 , in total 71 differential cross section data points. Using this data set, we will make a detailed comparison of the bootstrap and the Bayesian methods to extract the statistical uncertainty. The extraction of the systematic uncertainty for these data was already discussed in Sec. 3.12. Before continuing, it is worth noting that from the proton data alone, the isospin of a given pole is not determined. One can, however, simply assign a given number of isoscalar and isovector poles besides the continuum contributions, which have a given isospin, as well as the ω and φ mesons. This ambiguity will be resolved once neutron data are also fitted, see Sec. 4.2.
We consider first the Bayesian analysis described in Sec. 3.12. We assume two sets of priors, the normal and the uniform distributions depicted in Fig. 16. In both cases, the constraints on the various couplings and masses discussed in Sec. 3.10 are already included. Note further that in case of the uniform distribution, the prior for the unknown mass m v 1 is biased towards smaller values. The resulting proton em radii are equal within 3 significant digits for these two different prior distributions, see Tab. 4. Note that the systematic errors of these data have already discussed in Sec. 3.12.
Next, we compare the radius extraction from the Bayesian and the bootstrap method, which are shown in Fig. 18 for the proton charge radius. As can be read off from this figure and also seen in Tab. 4, the results are very similar,  Table 4. Statistical uncertainty in the proton electromagnetic radii from the PRad data using two different Bayesian distributions and the bootstrap approach.
with the bootstrap having slightly larger errors, which is due to our conservative choice considering the 99% quantile. The resulting normalized form factors G E (Q 2 )/G dip (Q 2 ) G M (Q 2 )/(µ p G dip (Q 2 )) as well as their uncertainties for the two methods are shown in Fig. 19. The differences are  negligible. The form factor ratio µ p G p E /G p M measured below Q 2 = 1 GeV [47,48] is also well described, as already displayed in Fig. 2 of Ref. [53] Note also that the magnetic form factor does not display any bump-dip structure below Q 2 < 1 GeV 2 as found in the MAMI analysis [44].
Having shown the equivalence of both methods here, in what follows we will stick to the bootstrap procedure, which is easier to implement in case of large data sets with a larger number of fit parameters.

Fits to proton and neutron data
We are now in the position to analyze the the full data set.
To be more precise, for the proton we fit to the cross section data from PRad [15] and from MAMI-C [44] as well as to the polarization transfer data from Jefferson Lab [133,134,135,136] above Q 2 = 1 GeV 2 (note that the data from Refs. [31,32] are updated in Refs. [133,136], respectively, and thus do not appear in the data base) together with the neutron form factor world data base already used in [45]. The size of the data base and the Q 2 -ranges we are fitting is provided in Tab. 5. We also include the constraint on the neutron charge radius squared, updated to the latest value given in Eq. (58) from Ref. [108]. Ultimately, we need to reassess the neutron data base by performing chiral EFT analyses of electron scattering of the deuteron and (polarized) 3 He. This, however, goes beyond the scope of the present work.  Table 5. Data base used in the fits.
Before showing the results of the best fit and the corresponding statistical and systematic uncertainties, it is worth pointing out we made extensive searches for solutions with altogther 36 combinations of isoscalar (is) and isovector (iv) poles, ranging from 3 + 3 to 8 + 8 is+iv poles, with the reduced χ 2 varying by less than 5%, in most cases even by less than 1%. We noticed that the fits with a larger number of is than iv poles turned out to be slightly better.
The best solution has 6 + 4 is+iv poles, and the fits to the ep cross section data with Q 2 < 1 GeV 2 , the proton form factor ratio with Q 2 > 1 GeV 2 , the electric form factor of the neutron and the magnetic form factor of the neutron are shown in Figs. 20, 21, 22, 23, respectively. The corresponding central values for the various vector mesons masses, vector meson couplings and the normalization constants of the MAMI and PRad data are collected in Tab. 6 in App. E. It is remarkable that while the isoscalar spectral functions requires a number of high mass effective poles, the effective isovector poles all have masses  below 2.3 GeV. We note that the vector coupling of the residual φ comes out small, consistent with expectations from the OZI rule. The tensor coupling is, however, quite large, but the next effective isoscalar pole has a comparable tensor coupling of opposite sign, see also the discussion in Sec. 4.3. We also note that the various normalization factors are deviating from one by less than 1%. Let us now discuss the predictions of and physics related to these fits. First, we extract the various radii from these fits, where the first error is statistical (based on the bootstrap procedure explained in Sec. 3.12) and the second one is systematic, based on the variations in the spectral functions discussed before. The values for the radii are completely consistent with earlier determination, cf. Ta-0.5 1 Fig. 23. Best fit to the neutron magnetic form factor data. For notations, see Fig. 20. bles 1,2, but with a much improved uncertainty estimation. Clearly, given the data set we fitted, the systematic uncertainty is largest for the neutron magnetic radius. The statistical uncertainty is small for all three radii. It is also interesting to give the radii of the Dirac and Pauli from factors in the isospin basis, in particular for the comparison with lattice QCD results. The reason for this is that the contribution of the so-called disconnected diagrams to the isoscalar FFs, which are notoriously difficult to calculate. These radii are given by: Similarly, the electric and magnetic radii in the isospin basis are (using the conventions given in Eqs. (6,7)) r s E = 0.773 ± 0.002 +0.002 −0.003 fm, r v E = 0.900 ± 0.002 ± 0.002 fm, r s M = 0.801 ± 0.008 +0.010 −0.038 fm, We note that the central values in Eq. (83) lead to the squared isovector radii, of (1/2)(r v E ) 2 = 0.405 fm 2 and µ v (r v M ) 2 = 1.72 fm 2 , which are perfectly consistent with the sum rule estimates in Eq. (43) but have, of course, much smaller uncertainties. It is also interesting to compare the isovector radii We also note that the value for the squared isovector with a recent state-of-the-art lattice QCD calculation at physical pion masses [137], While the value of the isovector charge radius is consistent with ours, the latttice value for the isovector magnetic radius is smaller than ours, that is there is some tension. It remains to be seen what future lattice calculations will give. As in earlier fits [46,53], the data for the proton form factor ratio µ P G p E /G p M for Q 2 < 1 GeV 2 , which do not participate in the fit, are well described, see the inset in Fig. 21. This points towards consistency between the twophoton corrected cross section data and the ratio data, that are not affected by such corrections. The situation is, however, different for larger momentum transfers. In Figs. 24,25 we display G p E (Q 2 ) and G p M (Q 2 ), that did not participate in the fits. Because the proton form factor ratio tends to zero at Q 2 8 GeV 2 , marked deviations from the dipole form are observed. Only at very large momentum transfer, the fall-off required by pQCD is observed. More precisely, we find that Q 4 F p,n 1 (Q 2 ) starts to level off beyond 30 GeV 2 , whereas that is not the case yet for Q 6 F p,n 2 (Q 2 ). Clearly, in this region of momentum transfer, more data are needed to pin down the form factors more precisely and to eventually see the onset of perturbative QCD. This is entirely consistent with earlier findings, see e.g. Refs. [26,34].
Note that the long-range part of the Breit-frame charge and magnetization distributions that follows from the Sachs form factors can be interpreted in terms of a "pion cloud" and some additional short-range contributions from the ρ and other short-ranged physics. However, we emphasize that this separation is scale-dependent and thus not unique. A general discussion of the pion cloud can be found in Refs. [138,139].

Vector meson couplings
As noted before, in the earlier DR analyses, in the isoscalar spectral function below 1 GeV, only the ω and the φ mesons were retained and thus using Eq. (51), one was able to extract the vector and the tensor couplings of these mesons. However, we have shown that in the region of the φ, the isoscalar spectral function is much more complicated, preventing one from extracting φ-meson couplings.
In what follows, we will thus only consider the ω-meson. In this case, we have The earlier fits, which had no restrictions on the residua led to a large vector and a small tensor coupling, from Ref. [22] and Ref. [26], respectively. This vector coupling is sizeably larger than from the determination us-ing forward dispersion relations in nucleon-nucleon scattering, g ωN N 1 = 10.1±0.9 [105]. This smaller value is, however, inconsistent with the approximate dipole behavior of F s 1 (Q 2 ) [140]. Note, however, that in one-boson-exchange models of the NN interaction, one typically finds values of g ωN N 1 (M 2 ω ) 20 which for typical strong form factors translates to g ωN N 1 (0) 10 [141]. Starting with the work of BHM [34], the isoscalar spectral function was considerably improved. In that work, the vector coupling was still large, but the tensor coupling could not be pinned down so precisely, Values within this range where also found in the analysis of the MAMI-C data combined with the proton form factor data for Q 2 > 1 GeV 2 and the neutron FF data base [45] where only central values were given. Finally, the analyses that concentrated mostly on the high-precision ep data from MAMI-C and PRad, the ω couplings took the values from Ref. [46] and Ref. [53], respectively. Finally, we present the results of our new fits, including the statistical and the systematical uncertainty: For the central values, the tensor-to-vector coupling ratio is small, κ ω = 0.06. We note that the uncertainties on the vector coupling are modest, they are much larger for the suppressed tensor coupling. Similar to the findings in BHM, the sign of the tensor coupling is not determined and the range of allowed values is sizeable.

Time-like form factors and final-state interactions
Before discussing the DR fits including the data from the time-like region, it is worth noting two very intriguing experimental findings related to the cross section σ(e + e − → pp) (and the reversed reaction) and the corresponding effective form factor |G p eff |. First, as shown in the upper panel of Fig. 26, there is a strong enhancement in the close-to-threshold region, as comparison with the phase space behavior (normalized to the data at about 50 MeV excess energy) clearly reveals. Note also that due to the Coulomb interaction between the proton and the antiproton, the cross section does not start at zero. No such effects are seen in σ(e + e − → nn). We note that such threshold enhancements are also observed in other processes like e.g. J/Ψ → xpp, Ψ (3686) → xpp with (x = γ, ω, ρ, π, η) and B + → K + pp, see e.g. Refs. [146,147,148]. Second, extending further out in momentum transfer, the BaBar data [50]  and also the BESIII data [142] exhibit some oscillating structures, most pronounced for invariant masses M pp below 2.5 GeV, see the lower panel in Fig. 26. The corresponding neutron data from FENICE [145] and SND [149] are less precise than the proton data, but show a similar behavior for q 2 4 GeV 2 . For recent fits to the timelike proton effective FF accounting for these structures, see [150].
DR fits including space-and time-like data were performed in Refs. [30,34,49]. Here, we focus on the work done in the latest paper. Though that work investigated some issues related to FSI in an exploratory way, it provides the most detailed information on the physics contained in the time-like FFs. In this work, the spectral functions was enlarged to account for the coupling to the newly established φ(2170) vector meson, as well as baryonic triangle graphs with virtual N N π, N ∆π and ∆∆π particles, the first of these giving a simple representation of the strong final-state interactions (FSI), see the discussion below.
sections and the ratio G E /G M from polarization observables on the scattering side in addition to the effective FF and |G E /G M | on the production side were included for the proton. In case of the neutron, G E and G M from scattering data and again the effective FF on the production side were considered. The spectral function included the the 2π, KK and ρπ-continuum, the ω-and φ-contribution and three/five isoscalar/isovector poles restricted in the mass range M V = 1 . . . 1.8 GeV. In addition, the new vector resonance was taken at a mass of 2.125 GeV and its width was determined in the fit, which turned out to be Γ = 0.088 GeV. Good agreement with the existing data, as shown in Fig. 27 for the proton and neutron effective form factors and in Fig. 28 for the proton form factor ratio in the space-like and the time-like region, was obtained. In particular, the SND data for the neutron effective FF show a very similar behavior to the proton effective FF over a large range, which can be well described in this approach. However, the range around the φ(2170) calls for further neutron measurements to allow for a determination of the isospin of the structures in G p eff . Let us now discuss the various threshold effects in the time-like data, in particular the strong enhancement at the pp threshold. This was first observed at LEAR in the inverse reaction pp → e + e − [152] and substantiated by the BaBar collaboration [153], which provided data for e + e − → pp down to the threshold region. As noted before, this threshold enhancement was also observed in a number of other production reactions such as J/ψ and B decays. Several explanations involving unobserved me- son resonances or scenarios that involve NN bound states (baryonium) have been put forward, see e.g. Ref. [146]. More conventional but plausible interpretations of this phenomenon were given in terms of the FSI between the proton and the antiproton, employing either a Migdal-Watson approximation or meson-exchange models of various levels of sophistication to describe the pp interaction, see e.g. the earlier works Refs. [154,155,156,157,158,159]. The latest and arguably most sophisticated approach to this phenomenon employs simple point-like form factors, whose energy dependence is entirely given by the protonantiproton FSI (or the pp initial state interactions in the annihilation process) [160]. The nucleon-antinucleon interaction is based on chiral effective field theory at NLO [161] and NNLO [162]. In this approach, the steep rise of the effective FF for energies close to the pp threshold is explained solely in terms of the pp interactions, cf. Fig. 29, consistent with the findings in Refs. [159,163,164,165,166,167,168,169,170,171]. Also existing experimental information (differential cross sections, form factor ratio) is quantitatively described in this approach. In addition, predictions for various spin-dependent observables, that can be tested in the future with PANDA at FAIR, are also given in that work. Note, however, that this framework is only applicable to the threshold region, that is up excess energies of about 100 MeV. Triangle graphs with virtual N∆π and ∆∆π states were also considered in Ref. [49], picking up an idea of Rosner [172], in order to approximate possible cusp ef- fects. However, the vertices are not well known for these kinematics, so the calculation was reduced to the scalar loop integral parameterized in terms of fitted strength parameter f N ∆/∆∆ ∼ O(1) multiplying each loop structure. The explicit momentum dependence of the vertices was accounted for in terms of overall form factors, with Λ N ∆/∆∆ the respective fitted cut-off parameter. Interestingly, performing fits with these structures instead of the explicit poles discussed before leads to an equally good description of the proton effective FF, very similar to what is shown in Fig. 27 (upper panel). Furthermore, in Ref. [49] the nucleon FFs were also discussed in region of t 0 = 4M 2 π < t < t thr = 4m 2 p which is not accessible by direct measurements, but by analytic continuation in t = q 2 = −Q 2 . In fact, an additional particle emission from the initial state proton can lower the energy of the (virtual proton) to reach below the threshold, as discussed in Ref. [173] for the process pp → e + e − π 0 . To get insight into the unphysical region, it is instructive to use a DR for the logarithm, see e.g. Refs. [174,175,176,177,178]. In principle, this also allows for a separation of the FF phase δ(t) and modulus in the representation G(t) = |G(t)|e iδ(t) . A once subtracted DR for the function ln where the first term vanishes due to the normalization G E (0) = G M (0)/µ p = 1. Experimental information on this integral equation (93) is available in the space-like region t < 0 on G(t) and in the time-like region for t > t thr on the modulus |G(t)|. The solution of this integral equation is not straightforward, it requires additional information to be included (as the problem is ill-conditioned). One possible solution was proposed in Refs. [176,177], namely to consider the integral contributions to the logarithm ln |G(t)| in the space-like region, using definite values for the known part above t thr In Ref. [49], as input for the left-hand-side of Eq. (94), the discretized result of a simultaneous fit to data in all accessible regions were used. A typical result for the modulus of G p E (t) is shown in Fig. 30. The enhancement just below the production threshold could signal the appearance of a broad baryonium state, but clearly more precise data on the time-like nucleon FFs are required to come to a definite conclusion here.

Summary and outlook
This paper served two purposes. First, we have reviewed the dispersion-theoretical approach to the electromagnetic form factors of the nucleon, with particular emphasis on the constraints posed by unitarity and analyticity on the spectral functions. Second, we have performed new fits including recent high-precision data on electron-proton scattering for squared momentum transfers Q 2 < 1 GeV 2 , the proton form factor ratio in the range Q 2 1 · · · 8.5 GeV 2 and the world data basis on the neutron electric and magnetic form factors, including the recent accurate extraction of the neutron charge radius squared from chiral nuclear EFT. We also have sharpened the toolbox to determine the statistical and systematical uncertainties. This led to a number of new results concerning the various nucleon electromagnetic radii, the form factors and the ωN N couplings. We would like to stress again that DRs have always found a small proton charge radius, r p E 0.84 fm, with a slightly larger proton magnetic radius, r p M 0.85 fm. As before, we find that the neutron magnetic radius is the largest, r n M 0.87 fm. Consistent with earlier analyses, the onset of pQCD is not seen in the existing form factor data. We have also discussed our present understanding of the physics in the time-like region, where a strong enhancement of the cross section for e + e − → pp (and in its reversed process pp → e + e − ) is observed, that can be understood in terms of proton-antiproton final-state interactions (or initial-state interactions for the reserved process). Furthermore, there are interesting oscillating structures in the cross section that require additional poles and/or thresholds.
Clearly, there are a number of issues that require more data and/or further investigations: -For the neutron data basis, a thorough analysis of the existing electron-deuteron and electron-3 He scattering data based on chiral effective field theory and including two-photon corrections should be performed. This would allow to consistently analyze the proton and neutron form factors based on the dispersive approach applied directly to cross section data. -A new combined analysis of the space-like data (as done here) with the time-like data should be performed, including the improved knowledge of the pp final-state interactions obtained in the last decade. This would also sharpen the predictions for future measurements with PANDA at the FAIR facility. -Data on ep scattering or the polarization transfer at Q 2 10 GeV 2 are urgently needed to investigate the onset of perturbative QCD. It will also be interesting to find out whether the form factor ratio really crosses zero as the present data seem to indicate.
-It would also be useful to improve our understanding of the spectral functions in the unphysical region, see Figs. 2 and 30. This requires more work based on logarithmic dispersion relations such as Eq. (93). Finally, let us point out the the dispersion-theoretical approach to the nucleons electromagnetic form factors has matured and become a precision tool to analyze electron scattering and form factor ratio data. In the future, it will also be extended to analyze the upcoming muon-proton scattering data from the MUSE [179] and AMBER [180].

A Neutron form factors from light nuclei
As there are no free neutron targets, the neutron form factors must be extracted from electron scattering off light nuclei. In this appendix, we briefly outline how this can be achieved in case of the simplest nucleus, the deuteron. The extension to systems with 3 or 4 nucleons is similar but more complicated.
The deuteron is a spin-1 particle. Using Lorentz invariance, time-reversal invariance as well as parity and current conservation, the matrix element of the deuteron em currents takes the form, see e.g. [67,181] (as usual, in units of the elementary charge e) where p , p are the deuteron four-momenta, λ , λ the corresponding helicities, m d = 1.8756 GeV is the deuteron mass and Q 2 ≡ −k µ k µ = −(p − p) 2 ≥ 0. Furthermore, the polarization four-vectors ξ µ are subject to the constraints ξ µ (p, λ)p µ = 0 and ξ * µ (p , λ )p µ = 0. Instead of the scalar functions G i (Q 2 ) (i = 1, 2, 3), one often uses the deuteron charge (G C ), magnetic (G M ) and quadrupole (G Q ) form factors given by with η = Q 2 /(4m 2 d ). These form factors are subject to the normalizations: in terms of the deuteron magnetic moment, µ d = 0.857 µ N , and the deuteron quadrupole moment, Q d = 0.286 fm 2 .
In the one-photon exchange approximation, the unpolarized elastic electron-deuteron (ed) scattering cross sec-a) b) Fig. 31. Electromagnetic response of the deuteron. Diagram a) depicts the IA, that is sensitive to the nucleon em FFs (black circle) whereas diagram b) denotes the so-called two-body corrections (depicted by the shaded box) as explained in the text. The hatched triangles denote a deuteron wave function and solid (wiggly) lines represent nucleons (photons).
tion in the laboratory (lab) frame is given by with E the energy of the incoming electron, θ the electron scattering angle in the lab frame and Q 2 ≥ 0 is the squared momentum transfer. The structure functions A(Q 2 ) and B(Q 2 ) are related to the three form factors of the deuteron via As can be seen, from the unpolarized cross section one can not disentangle the charge and the quadrupole form factors. This can be achieved by considering polarization data, e.g. the tensor analyzing power T 20 (Q 2 , θ) is sensitive to a different combination of the three deuteron FFs.
Using a chiral expansion (or a meson-exchange model in older times), these deuteron FFs can now be related to the single nucleon form factors as depicted in Fig. 31. The so-called impulse approximation (IA), where the photon couples to one nucleon, is depicted in graph a). This contribution is evidently sensitive to the single nucleon em FFs. To be precise, one takes the proton FFs as given and is then entirely sensitive to G n E and G n M , see below. There are, however, a number of corrections as displayed by graph b) in Fig. 31. This include one-and two-pion exchange currents as well as photon-four-nucleon vertices (or heavy meson exchanges in the older language). The status of the nuclear currents based on a systematic evaluation of the few-body wave functions and the current operators is given in Ref. [182]. Consider for example the charge density that generates the charge form factor. As the deuteron is an isoscalar, its em response is entirely sensitive to the isoscalar form factors. Modulo higher order corrections, at leading order in the chiral expansion, the charge density gives rise to the form factors: ∞ 0 (u 2 (r) + w 2 (r)) j 0 kr 2 dr , where one works in the Breit frame with k 2 = Q 2 , the direction of the photon momentum is taken along the positive z-axis and k = |k|. Also, u(r) and w(r) are the deuteron S-and D-wave wave functions, normalized to one, see e.g. [141], and j 0 (x), j 2 (x) are the conventional Bessel functions. Corrections to this leading order results can be worked out straightforwardly, these are given for a chiral EFT approach in Ref. [109] (and references therein).

B Pion-nucleon scattering in the unphysical region
Here, we briefly discuss the subthreshold expansion of the πN amplitudes, which proceeds in terms of the variables ν = (s − u)/(4m) and t around ν = t = 0 , a ± mn ν 2m t n , where the upper/lower entry corresponds to I = ±, and the a ± mn , b ± mn are the subthreshold parameters. TheĀ,B are the Born-term-subtracted amplitudes, defined as X ± (ν, t) = X ± (ν, t) − X ± pv (ν, t), X ∈ {A, B}, (102) with Here, the subscript 'pv' refers to the pseudovector πN coupling as required from chiral symmetry and g denotes the πN coupling constant, to be identified later with g c for the charged-pion vertex. In the RS analysis of πN scattering, the value g 2 c /(4π) = 13.7(0.2) [91] has been used. This value is in line with the most recent determination from nucleon-nucleon scattering [183,184]. More details on the subthreshold expansion of the πN scattering amplitude is given in Refs. [88,52].

C Analysis of the three-pion contribution
Here, we briefly review the arguments of Ref. [96] that there is no strong enhancement on the left shoulder of the ω resonance, the lowest vector meson in the three-pion channel (for a recent update, see Ref. [97]). The imaginary parts of the isoscalar electromagnetic form factors open at the three-pion threshold t 0 = 9M 2 π . The three-pion cut contribution is given by where the symbol 'A' refers to the γ → 3π and 'B' to the 3π →N N transition, respectively, and dΓ 3 is the measure on the invariant three-body phase space. This can be explicitly worked out in baryon chiral perturbation theory [185], where to leading order in the small parameter p, namely at order O(p 7 ), the two-loop diagrams shown in Fig. 32 can contribute to the isoscalar imaginary parts, however, graph (d) vanishes because of an isospin factor zero. The isoscalar imaginary parts in the heavy nucleon limit m → ∞ can be given in compact form (note that this corresponds to switching off all higher order corrections starting at O(p 8 )) Im G S M (t) = g A m (8π) 4 F 6 π L(t) 3t 2 − 10tM 2 π + 2M 4 π +g 2 A 3t 2 − 2tM 2 π − 2M 4 π +W (t) t 3 + 2t 5/2 M π − 39t 2 M 2 π − 12t 3/2 M 3 π +65tM 4 π − 50 √ tM 5 π − 27M 6 π +g 2 A 5t 3 + 10t 5/2 M π − 147t 2 M 2 π + 36t 3/2 M 3 π +277tM 4 π − 58 √ tM 5 π − 135M 6 π , with and in terms of the kinematical variables of the two independent pions (the third one can be expressed in terms of these and the nucleon momenta), The behavior near threshold t 0 = 9M 2 π of the imaginary parts for finite pion mass, Eqs. (105,106), is given by which corresponds to a stronger growth than pure phase space This feature indicates (as in the isovector case) that in the heavy nucleon mass limit m → ∞ normal and anomalous thresholds coincide. In order to find these singularities for finite nucleon mass m an investigation of the Landau equations is necessary [186]. By using standard techniques [186] one finds one anomalous threshold of diagrams (a) and (b) at √ t c = M π 4 − M 2 π /m 2 + 1 − M 2 π /m 2 , which is very near to the (normal) threshold t 0 = 9M 2 π and indeed coalesces with t 0 in the infinite nucleon mass limit. Note that diagram (d) does not possess this anomalous threshold, but only the normal one.
The resulting spectral distributions weighted with 1/t 2 are shown in Fig. 33. Very different to the isovector spectral functions discussed in Sec. 3.4, they show a smooth rise and are two orders of magnitude smaller than the corresponding isovector ones, cf. Fig. 6. There is indeed no enhancement of the isoscalar electromagnetic spectral function near threshold. Even though the isoscalar and isovector electromagnetic form factors behave formally very similar concerning the existence of anomalous thresholds t c very close to the normal thresholds t 0 , the influence of these on the physical spectral functions is rather different for the two cases. Only in the isovector case a strong enhancement is visible. This is presumably due to the different phase space factors, which are (t − t 0 ) for the isovector and isoscalar case, respectively. In the latter case, the anomalous threshold at t c = 8.9M 2 π is thus effectively masked. This justifies the standard DR approach of only taking the ω-meson as the lowest pole in isoscalar spectral function.