Transverse charge and current densities in the nucleon from dispersively improved chiral effective field theory

Background: The transverse densities $\rho_{1, 2}(b)$ describe the distributions of electric charge and magnetic moment at fixed light-front time and connect the nucleon's elastic form factors with its partonic structure. The dispersive representation of the form factors $F_{1, 2}(t)$ expresses the densities in terms of exchanges of hadronic states in the $t$-channel and permits their analysis using hadronic physics methods. Purpose: Compute the densities at peripheral distances $b = \mathcal{O}(M_\pi^{-1})$, where they are generated predominantly by the two-pion states in the dispersive representation. Quantify the uncertainties. Methods: Dispersively improved chiral effective field theory (DI$\chi$EFT) is used to calculate the isovector spectral functions $\textrm{Im}\, F_{1, 2}(t)$ on the two-pion cut. The method includes $\pi\pi$ interactions ($\rho$ resonance) through elastic unitarity and provides realistic spectral functions up to $t \approx$ 1 GeV$^2$. Higher-mass states are parametrized by effective poles and constrained by sum rules (charges, radii, superconvergence relations). The densities $\rho_{1, 2}(b)$ are obtained from their dispersive representation. Uncertainties are quantified by varying the spectral functions. The method respects analyticity and ensures the correct $b \rightarrow \infty$ asymptotic behavior of the densities. Results: Accurate densities are obtained at all distances $b \gtrsim 0.5$ fm, with correct behavior down to $b \rightarrow 0$. The region of distances is quantified where transverse nucleon structure is governed by the two-pion state. The light-front current distributions in the polarized nucleon are computed and discussed. Conclusions: Peripheral nucleon structure can be computed from first principles using DI$\chi$EFT. The method can be extended to generalized parton distributions and other nucleon form factors.

Transverse densities have emerged as a key concept in nucleon structure physics. The functions ρ 1,2 (b) describe the transverse coordinate distributions of charge and current in the nucleon at fixed light-front time x + = x 0 + x 3 and provide a spatial representation appropriate to the relativistic nature of the dynamical system [1][2][3][4]. They are defined as two-dimensional Fourier transforms of the invariant form factors (FFs) F 1,2 (t) parametrizing the current matrix element between nucleon states. At the same time, they represent a projection of the generalized parton distributions (GPDs) describing the distribution of partons in light-front longitudinal momentum and transverse position [2,3]. As such, the transverse densities connect the FFs measured in low-energy electronnucleon elastic scattering with the partonic structure probed in high-energy processes such as deep-inelastic scattering and hard exclusive processes.
The nucleon FFs F 1,2 (t) at spacelike momentum transfers t < 0 can be interpreted in terms of hadronic exchanges between the current and the nucleon. The mathematical framework is provided by the dispersive representation of the FFs based on analyticity in t. The FFs at t < 0 are expressed as integrals over their imaginary parts on the cut at t > t thr > 0, Im F 1,2 (t), the so-called spectral functions, which correspond to t-channel states with definite hadronic composition and quantum numbers. A similar dispersive representation can be derived for the transverse densities [5]. It establishes a correspondence between the densities at a given distance b and the arXiv:2204.11863v1 [hep-ph] 25 Apr 2022 hadronic exchanges at various masses t > 0 [6]. In particular, it connects the densities at large distances b 1 fm with the lowest-mass t-channel states and permits a systematic study of "peripheral" nucleon structure [5,[7][8][9]. Using this framework one can compute the peripheral densities from first principles employing methods of hadronic physics. One can also explore the duality between the hadronic exchanges in the t-channel and the partonic structure in s-channel [6].
The lowest-mass t-channel state in the nucleon electromagnetic FFs is the two-pion (ππ) state. It appears in the isovector channel and saturates the isovector spectral functions up to t ≈ 1 GeV 2 . The ππ system in this mass region interacts strongly and forms the ρ resonance at t ≈ 0.6 GeV 2 . The picture of "vector dominance" abstracted from this situation explains many observations in the phenomenology of the electromagnetic form factors of the nucleon and other hadrons. The spectral functions in the ππ channel have been constructed empirically using methods of hadronic amplitude analysis, such as elastic unitarity with input from πN scattering data [10,11], or Roy-Steiner equations [12]. The transverse densities have been studied using the dispersive representation with such empirical spectral functions [6].
It would be interesting if the transverse densities could be computed using spectral functions derived from chiral effective field theory (chiral EFT). This approach would open up several new possibilities. First, chiral EFT is predictive and allows one to reduce the information content of the spectral functions and densities to a few universal parameters (low-energy constants), which are determined from independent measurements. Second, one can quantify the uncertainties in the peripheral densities using the parametric expansion of chiral EFT. Third, because chiral EFT is a point-particle field theory, one can explore the duality between t-channel exchanges and schannel structure at a microscopic level and derive a partonic representation of the chiral processes. Fourth, because of the universality of chiral EFT one can relate the electromagnetic densities to other elements of peripheral nucleon structure.
Traditional chiral EFT calculations of the spectral functions are limited to the near-threshold region t − 4M 2 π ∼ few M 2 π and cannot describe the ρ resonance region, because the strong ππ interactions amount to large higher-order corrections [13][14][15][16]. When applied to the transverse densities, this only allows one to compute the densities at very large distances b 2 fm, where they are extremely small and of little practical interest [5,[7][8][9]. In order to go to smaller distances, one needs an approach that takes into account the ππ interactions in the ρ region in a different manner. Dispersively improved chiral EFT (DIχEFT) is an approach that incorporates ππ interactions through elastic unitarity and enables EFTbased calculations of the spectral functions in the ρ meson region [17][18][19] (an equivalent alternative formulation is described in Ref. [20]). Using the N/D method, the spectral function is separated into a part containing the nonperturbative interactions in the ππ system, which is taken from the measured pion timelike form factor, and a part describing the coupling of the ππ system to the nucleon, which can be computed with chiral EFT with good convergence. The DIχEFT spectral functions have been used to compute the electromagnetic form factors [18,19] and extract the proton radii from elastic scattering data [21,22]; the approach has also been applied to the nucleon scalar FF [17]. A first study of transverse densities has been performed in leading-order (LO) accuracy [23].
In this work we use DIχEFT to compute the peripheral transverse charge and magnetization densities ρ 1,2 (b) in the nucleon and study their properties. We construct the ππ spectral functions Im F 1,2 (t) in partial next-toleading-order (N2LO) accuracy. High-mass states are described by effective poles, whose parameters are fixed by dispersive sum rules and superconvergence relations. We compute the transverse densities in the dispersive representation and quantify the region of distances where they are dominated by the ππ state. We quantify the uncertainties of the densities resulting from the low-energy constants and the high-mass poles. We also compute the transverse light-front current densities and construct 2dimensional images of the transversely polarized nucleon. We discuss the interpretation of the results and possible applications to physics studies with transverse densities.
Novel aspects of the present study of transverse densities are as follows: (a) The dispersive representation respects the analytic properties of the FFs (position and strength of singularities) and produces densities with the correct asymptotic behavior at b → ∞. This makes it possible to reliably compute the peripheral densities and estimate their uncertainties. Methods based on the Fourier transform of empirical FFs become unstable at large b and are not adequate for peripheral densities [24]. (b) The DIχEFT approach incorporates ππ interactions and the ρ resonance and produces realistic spectral functions up to t ≈ 1 GeV 2 , which allows one to compute the densities down to distances b 0.5 fm, substantially smaller than possible with traditional chiral EFT [5,[7][8][9]. This means that a large fraction of transverse nucleon structure is now amenable to an EFT description and can be deconstructed in terms of effective degrees of freedom and low-energy constants, representing a significant gain of information. It also means that the EFT description at large distances can be matched with a quark modelbased description of the densities at distances b 1 fm, enabling studies of quark-hadron duality in the transverse densities. (c) The uncertainties of the densities are estimated in the context of the dispersive representation, by varying the elements of the spectral functions. The unknown high-mass part of the spectral functions is parametrized by a random ensemble of high-mass poles, whose distribution is constrained by stability criteria imposed on the spacelike form factors. This new formulation minimizes the model dependence in the description of the high-mass states [18,19] and enables robust uncertainty estimates for the densities.
The article is organized as follows. In Sec. II we describe the methods used in the present study, including the properties of the transverse densities, their dispersive representation, the construction of the DIχEFT spectral functions, and the uncertainty estimates in the dispersive representation. In Sec. III we describe the results, including the DIχEFT isovector spectral functions, the isovector transverse densities, the proton and neutron densities, and the light-front current densities in the polarized nucleon. In Sec. IV we discuss the results and outline possible future applications to quark-hadron duality and other structures. In Appendix A we collect the nucleon radii and their uncertainties, which serve as input parameters in the DIχEFT calculation. In Appendix B we summarize the formulas for the N functions appearing in the calculation of the spectral functions with the N/D method. In Appendix C we describe the parametrization of the isoscalar spectral functions used in the calculation of proton and neutron densities.

A. Transverse densities
The transition matrix element of the electromagnetic current operator between nucleon (proton, neutron) states with 4-momentum transfer ∆ ≡ p − p is described by the FFs F 1 (t) and F 2 (t) (Dirac and Pauli FFs); see Ref. [23] for details. They are functions of the invariant momentum transfer t = ∆ 2 , with t < 0 in the physical region of elastic scattering. Their values at t = 0 are given by the nucleon charges and anomalous magnetic moments (in units of nuclear magnetons), The isovector and isoscalar components are defined as the same convention is used for other quantities (densities, radii; see below). The FFs are invariant functions and can be analyzed without specifying the form of relativistic dynamics or choosing a reference frame. Their interpretation in terms of spatial distributions requires specific choices. In the light-front form of relativistic dynamics one follows the evolution of strong interactions in light-front time x + ≡ x 0 + x 3 and describes the structure of systems at fixed x + [25][26][27]. In this context it is natural to consider the FFs in a class of reference frames where the 4-momentum transfer has only transverse components, ∆ 0 = ∆ 3 = 0, ∆ T ≡ (∆ 1 , ∆ 2 ) = 0, with |∆ T | 2 = −t. The FFs can then be represented as 2-dimensional Fourier integrals over a transverse coordinate variable b, with b ≡ |b|, Interpretation of the transverse densities in a proton state localized at the transverse origin, x = y = 0, Eq. (7). ρ1(b) describes the spin-independent J + current at transverse radius b. ρ2(b) describes the spin-dependent distortion in a nucleon polarized in the y-direction.
(the two-dimensional Fourier integrals can also be expressed as radial Fourier-Bessel integrals) [1][2][3][4]. The functions ρ 1,2 (b) are the transverse densities. Their spatial integrals reproduce the total charge and anomalous magnetic moment, The functions ρ 1,2 (b) describe the transverse spatial distributions of electric charge and anomalous magnetic moment in the nucleon at fixed light-front time. The distributions are frame-independent (they are invariant under longitudinal light-front boosts and transform kinematically under transverse boosts) and provide a spatial representation appropriate to the relativistic nature of the dynamical system; see Refs. [2,7,28] for discussion. A simple interpretation of the densities can be provided in nucleon states which are localized in transverse position (see Fig. 1; we use x, y, z to denote the 1, 2, 3 spatial directions) [2,7]. In a state where the nucleon is localized at the transverse origin, and its spin quantized along the y-direction, the expectation value of the current J + at the transverse position b is (same for p → n). Here (...) represents a factor resulting from the normalization of states [7], φ is the angle of the vector b relative to the x-axis, and S y = ±1/2 is the expectation value of the spin in the y-direction. m N is the nucleon mass (same for p and n). Thus ρ 1 (b) describes the spin-independent part of the plus current in the localized nucleon state, and ρ 2 (b) describes the spinand angle-dependent part of the current in a transversely polarized nucleon. Note that ρ 2 (b) satisfies the integral relation which is obtained from Eq.(6) by integration by parts. The derivatives of the FFs at t = 0 are related to the average squared transverse radii of the distributions, The factor 4 results from the 2-dimensional distributions and replaces the well-known factor 6 in the representation of the FF derivatives in terms of conventional 3dimensional radii. While the 3-dimensional radii have a physical interpretation only in nonrelativistic systems, the transverse radii here are averages of spatial distributions that have a well-defined meaning for relativistic system. The relation between the transverse radii and the 3-dimensional Dirac and Pauli radii r 2 1,2 is The nucleon radii are used as parameters in the dynamical calculations of the transverse densities in this work. Because form factor phenomenology usually quotes the 3-dimensional nucleon radii, we present our calculations such that they use the 3-dimensional radii as input, keeping in mind that they are related to the transverse radii by Eq. (12). The empirical values of the 3-dimensional radii and their uncertainties are summarized in Appendix A and will be quoted in the following. In studies of the nucleon's partonic structure in QCD one considers the transverse coordinate distributions of quarks and antiquarks with a given light-cone momentum fraction x in the proton, f a (x, b) andf a (x, b), where a = u, d, ... denotes the quark flavor. They are defined as the Fourier representation of the GPDs H a (x, ξ = 0, t), which describe the form factors of partons with light-cone plus momentum fraction x in the proton, in the situation where the plus momentum difference between the proton states is ξ = 0 and the momentum transfer has only transverse components [2,3], In this context the transverse charge density ρ p 1 (b) represents the integral over x of the difference of the proton's quark and antiquark distributions at transverse distance b, weighted by the quark charges e a , which can be interpreted as the cumulative charge of the partons in the proton at the transverse radius b. Equivalently, ρ p 1 (b) represents the Fourier transform of the first moment of the charge-weighted GPDs, A similar relation connects the density ρ p 2 (b) with the proton helicity-flip GPDs E a (x, ξ = 0, t) [2,3]. The transverse densities are thus directly related to the nucleon's transverse partonic structure in QCD.
The concepts of light-front quantization and partonic structure referenced here are used only for the interpretation of the transverse densities but are not needed for their computation. The densities are simple Fourier transforms of the invariant FFs and can be computed using hadronic physics methods such as dispersion theory and effective field theory. It is this "dual" character that makes the transverse densities so useful for nucleon structure studies.

B. Dispersive representation
The nucleon FFs are analytic functions of t. The physical sheet has a principal cut at positive real t > t thr ; the threshold t thr depends on the isospin channel (see below). The FFs satisfy unsubtracted dispersion relations, which express the functions at complex t as integrals over their imaginary parts on the cut. The real functions Im F i (t ) are known as the spectral functions. They correspond to processes in which the electromagnetic current with timelike momentum transfer t > t thr couples to the nucleon through a hadronic state in the t-channel. These processes occur in the unphysical region below the twonucleon threshold, t thr < t < 4m 2 N , where the spectral functions cannot be measured directly and have to be constructed using theoretical methods; see Ref. [29] for a review. In the isovector FFs the lowest-mass t-channel state is two-pion state with t thr = 4M 2 π ; in the isoscalar FF it is the three-pion state with t thr = 9M 2 π . The transverse densities Eqs. (4) can be computed as the Fourier transform of the dispersive representation of FFs, Eq. (17) [5,6]. One obtains a dispersive (or spectral) representation of the densities as (19) where K 0 and K 1 are the modified Bessel functions of the second kind. Equations (18) and (19) express the densities at a given distance b as a superposition of contributions of t-channel states (or exchanges) with squared mass t. The modified Bessel functions decay exponentially at large arguments, The dispersive integrals for the densities therefore converge exponentially at large t, and the integration effectively extends over masses This provides a mathematical formulation of the connection between the masses of the exchanges and the ranges of the spatial distributions in the nucleon. In particular, the peripheral densities at b = O(M −1 π ) are governed by lowest-mass states in the dispersive representation and can be computed and analyzed accordingly. Note that the actual spectral composition of the densities -how much the states with various √ t contribute to the densities at given b -depends on the distribution of strength in the spectral functions and can be established only with specific models of the latter.
The dispersive representation offers several theoretical and practical advantages for studying the peripheral transverse densities compared to other approaches. (a) The dispersive representation permits efficient calculation of the peripheral densities. The exponential convergence of the dispersion integrals reduces the contribution from the high-mass region, where the spectral functions are poorly known. Calculations can focus on the low-mass region -the two-pion part of the isovector spectral function, for which dedicated theoretical methods are available. The high-mass region can be parametrized by effective poles, whose coefficients are fixed by sum rules; only the overall strength in this region is relevant to the peripheral densities, not the details of the distribution. (b) The dispersive representation automatically generates densities with the correct asymptotic behavior at b → ∞. The asymptotic behavior of the densities at b → ∞ is governed by the analytic properties of the FFs in t (position and strength of singularities), which are explicitly realized in the dispersive representation. The densities exhibit an exponential decay with a range governed by lowest-mass exchanges, modified by a pre-exponential factor resulting from the behavior of the spectral function near the threshold [7]. The spectral integrals Eqs. (18) and (19) permit stable numerical evaluation of the densities in the region where they are exponentially small. 1 In contrast, methods calculating the densities as Fourier transforms of the FFs become numerically unstable at large b, even if the FF parametrization have correct analytic properties. 2 (c) The dispersive representation enables uncertainty estimates of the peripheral densities. The densities generated by Eqs. (18) and (19) depend smoothly on the parameters of spectral function, even at large b where they are exponentially small. Varying the parameters of the low-mass spectral functions one can estimate the uncertainties of the peripheral densities in a manner that respects analyticity and is numerically stable. Methods using the frequency spectrum of the Fourier transform of FFs for estimating the uncertainties of the densities are not appropriate for distances significantly larger than 1 fm [24].
The spectral functions obey certain integral relations (sum rules), which result from the constraints on the nucleon FFs at t = 0 and their derivatives in the dispersive representation Eq. (17), Additional relations follow from the asymptotic behavior of the nucleon FFs at large spacelike |t|, This behavior is predicted by the QCD hard scattering mechanism up to logarithmic corrections (counting rules) [27]. In F 1 the predicted behavior is approximately observed in data at |t| 1 GeV 2 ; in F 2 the predicted behavior is not observed at presently available momentum transfers [31]; see Ref. [32] for a possible theoretical explanation in the context of QCD. When imposed on the dispersive representation of the FFs, Eq. (26) implies the relations (so-called superconvergence relations) At the level of the densities, these relations constrain the behavior in the limit b → 0. This can be derived from the dispersive representation, Eqs. (18) and (19), by using the limiting behavior of the modified Bessel functions at small argument (here z ≡ √ tb), where "analytic" denotes the parts that are analytic at z = 0 (constants or positive powers). Equation (27) implies that The kernel K 0 ( √ tb) in Eq. (18) diverges logarithmically for b → 0 according to Eq. (30); this would cause a logarithmic divergence of ρ 1 (b); however, the coefficient of the divergent term is zero because of Eq. (27). 3 Similarly, Eqs. (28) and (29) imply that

C. Spectral functions
The spectral functions of the isovector nucleon FFs on the two-pion cut have been computed using various theoretical approaches, such as analytic continuation of pion-nucleon amplitudes [10], chiral EFT [13][14][15][16], and Roy-Steiner equations [12]. Here we employ the method of DIχEFT, which combines general methods of dispersion theory (elastic unitarity in the ππ channel, N/D method) with specific dynamical input from chiral EFT. The foundations of the method and its applications to FFs are described in detail in Refs. [17][18][19]. Here we summarize only the main steps and the new features arising in the present application to densities. The new features are: (a) We now construct the spectral functions Im F 1,2 (t) as needed for the transverse densities, whereas in Refs. [18,19] we worked with Im G E,M (t). (b) We impose the superconvergence relations resulting from the asymptotic behavior of F 1,2 (t), which provide additional constraints on the high-mass part of the spectral functions. (c) We implement a more flexible parametrization of the high-mass part of the spectral function to enable more realistic uncertainty estimates.
The isovector spectral functions are organized as where Θ is the step function (1 if the argument is true, otherwise 0). The first term is the contribution of the two-pion cut that starts at 4M 2 π and extends up to t max ≈ 1 GeV 2 (see below); this part is calculated theoretically. The second term represents the contribution of high-mass states above t max of unspecified hadronic composition; this part is parametrized through effective poles.
The two-pion spectral functions are obtained from the elastic unitarity relation in the two-pion channel (see Fig. 2). The relation is written in a manifestly real form by applying the N/D method [35][36][37] k cm is the center-of-mass momentum of the ππ state in t-channel, Γ i (t)(i = 1, 2) are the ππ → NN partial-wave amplitudes, and F π (t) is the pion timelike form factor. The functions Γ i (t) and F π (t) in Eq. (35) are complex for t > 4M 2 π because of ππ rescattering but have the same phase. The ratios N i (t)(i = 1, 2) in Eq. (36) are real for t > 4M 2 π and possess only left-handed singularities; they are free of ππ rescattering effects and describe the coupling of the ππ system to the nucleon. These function can be computed in χEFT with good convergence. The ππ rescattering effects in Eq. (36) are contained in the squared modulus of the pion form factor |F π (t)| 2 . This function is measured in e + e − → π + π − exclusive annihilation experiments and can be taken from a parametrization of the data [38]. In this way Eq. (36) factorizes the rescattering effects in the ππ channel and permits chiral EFT-based computation of the spectral functions. While Eq. (35) and Eq. (36) are strictly valid only up to the 4-pion threshold t < 16 M 2 π , the e + e − exclusive annihilation data indicate that 4-pion and other states are strongly suppressed up to ∼ 1 GeV 2 [38], so that the elastic unitarity relations can practically be used up to t max = 1 GeV 2 [10,35].
We have computed the functions N i (t) in relativistic χEFT with N and ∆ intermediate states at LO and NLO accuracy. Contributions at N2LO accuracy have been estimated assuming that the loop corrections have the same functional form as the tree-level result (partial N2LO, or pN2LO approximation). The χEFT results for N i (t) can be obtained from those of the functions J 1 ± (t), which appear in the N/D representation of the unitarity relation for the G E,M (t) form factors and were calculated in our earlier studies [18,19]. The explicit formulas are given in Appendix B.
The LO and NLO results for the functions N i (t) are given in terms of known low-energy constants in the chiral Lagrangian. The pN2LO estimates involve one unknown parameter, λ i , describing the size of the loop corrections relative to the tree-level result. The total results for the function N i (t) are thus given by (i = 1, 2) The explicit form of N i (t)[N2LO-tree] is given in Eq. (B10). In summary, at pN2LO accuracy the two-pion part of each of the spectral functions in Eq. (34) involves one free parameter, λ i (i = 1, 2). The high-mass part of the isovector spectral functions in Eq. (34) is parametrized by effective poles. It is important to note that for our calculation of peripheral densities we need only a summary description of the high-mass strength of the spectral function, because of the strong numerical suppression of large t in the dispersion integrals Eqs. (18) and (19). (see Sec. II B). We construct an appropriate parametrization by making reasonable assumptions about the positions of the poles, treating variations of the position as part of the theoretical uncertainty, and and fixing the strength of the poles through the dispersive sum rules Eq. (22)  Several observations suggest that the main strength of the isovector spectral functions beyond the two-pion region is located around t ≈ 2 GeV 2 , and that higher values of t are strongly suppressed. (a) The e + e − → hadrons exclusive annihilation data show that the cross section at t > 1 GeV 2 is dominated by the 4π channel and concentrated around t ≈ 2 GeV 2 [38]. This suggests similar behavior of the nucleon spectral functions, even if the connection with the annihilation cross section is at the amplitude level and cannot be made explicit. (b) Dispersive fits to the spacelike nucleon FFs with flexible parametrizations of the high-mass states using multiple effective poles find most of the strength in the region around t ≈ 2 GeV 2 [39,40]. (c) The dual resonance model describes the isovector spectral functions of the pion or nucleon FFs through the exchange of vector resonances with masses M 2 n = M 2 ρ (1 + 2n), with M 2 0 ≡ M 2 ρ . The first resonance after the ρ has mass M 2 1 = 3 M 2 ρ = 1.8 GeV 2 . If in the nucleon FFs the resonance contributions decrease rapidly with n, the dominant contribution beyond the ρ should come from this resonance.
Based on these observations, we parametrize the highmass part of the isovector spectral function Im F V 1 (t) as The first term is a delta function, the second term is the derivative of a delta function. The pole masses have the nominal value Their actual values are considered undetermined and will be allowed to vary randomly in a plausible range; their distribution will be constrained by further physical requirements, and the resulting variation in physical quantities will be regarded as a theoretical uncertainty of the model (see Sec. II D). Note that the sum of the functions in Eq. (41) parametrizes both the "strength" and the "shape" of the high-mass spectral function in region around t ≈ 2 GeV 2 in an effective form. 4 Combining the two-pion part given by Eqs. (36), (39), and (40), and the high-mass part given by Eq. (41), the isovector spectral function Im F V 1 (t) contains three unknown parameters: We fix these parameters through the sum rules for the isovector charge and radius, see Eqs. (22) and (23), and the superconvergence relation for F V 1 , see Eq. (27): where The relations Eqs. (44)-(46) constrain weighted integrals of the total isovector spectral function Im F V 1 (t). The integrals extend over the two-pion and the high-mass parts of the spectral functions (here n = 0, 1, 2), The integral over the two-pion part is a continuous integral and computed numerically; the integral over the high-mass part is a sum of delta function derivative integrals and computed exactly. Since the integrands depend linearly on the parameters Eq. (43), one obtains a system 4 In dispersive fits of the nucleon isovector FFs, the high-mass states are traditionally parametrized as a sum of simple delta functions at different positions. These parametrizations give clusters of poles around ∼ 2 GeV 2 with varying signs and unnaturally large coefficients 1 [39], which can effectively be combined to a sum of a single delta function and delta function derivatives. In this sense our parametrization Eq. (41) is equivalent to the traditional sum of simple delta functions. An advantage of our parametrization is that all coefficients have natural size ∼ 1, and that one can directly implement the feature that all the "structures" arise from the same mass region. of linear equations for the parameters, which can easily be solved.
In the spectral function Im F V 2 (t), we parametrize the high-mass part as which compared to Eq. (41) includes also a term with a second derivative of a delta function. The three pole masses again have the nominal value and will be allowed to vary randomly in an interval around this value (see Sec. II D). The combined spectral function now contains four unknown parameters: They are fixed through the sum rules for the isovector magnetic moment and the magnetic radius, see Eqs. (24) and (25), and by the two superconvergence relations for F 2 , see Eqs. (28) and (29): where In our approach the unknown parameters in the spectral functions are fixed by the dispersive sum rules and expressed in terms of the values and derivatives of the FFs at t = 0. Since values of the FFs (charges and magnetic moments) are known, this leaves the derivatives (radii) as the effective parameters of our model. With the spectral functions determined by the radii, our approach can then predict the spacelike FFs and the densities in terms of the radii. This particular "information flow" is made possible by the analytic properties of the FFs, which relate integrals over the spectral functions to derivatives of the FF at t = 0.  When computing the densities in the dispersive representation, Eqs. (18) and (19), the integrals receive contributions from the two-pion and the high-mass parts of the spectral functions, Eq. (34). An important question is how the relative contributions depend on the distance b, which represents an external parameter in the integrals. Figure 3 shows the two-pion and high-mass contributions to ρ V 1 (b) obtained with our spectral functions (here, for the the nominal parameter values). The top panel shows the absolute contributions; the bottom panel shows the relative contributions, i.e., the fractions of ρ V 1 (b) due to the two-pion and high-mass parts. One observes that in ρ V 1 (b) the two-pion part accounts for > 80% of at b > 1 fm, and > 60% at b > 0.5 fm. Similar relative contributions are found in the dispersion integral for ρ V 2 (b) (not shown in the figure). In ρ V 2 (b) the two-pion part accounts for > 97% of at b > 1 fm, and > 50% at b > 0.5 fm. These findings are central to our approach, as they quantify the dominance of the two-pion state at large distances and justify the summary description of the high-mass states for the purpose of computing the peripheral densities.
In the present study our focus is on the isovector chan-nel, where the two-pion state in the spectral functions generates the dominant contributions to the peripheral nucleon densities. In order to compute the individual proton and neutron densities in the dispersive representation, we need also the isoscalar spectral function. A parametrization of the isoscalar spectral functions, constructed along similar lines as for the isovector spectral function but relying more on empirical information, is described in Appendix C.

D. Uncertainty estimates
Our dispersive approach allows us to estimate the uncertainties of the spectral functions and the densities obtained from them. We consider two sources of uncertainties: (I) Uncertainties due to the parametrization of the high-mass part of the spectral functions. The high-mass part of the isovector spectral function is linked to the lowmass part through the dispersive sum rules, Eqs. (44)- (46) and Eqs. (53)-(56). The parametrization of the high-mass part can therefore influence the low-mass part of spectral functions and indirectly affect observables sensitive to the low-mass part, such as the peripheral densities. We estimate this uncertainty by varying the positions of the high-mass poles in the isovector spectral functions, Eqs. (42) and Eqs. (51). As the plausible range of variation we consider This range allows for variations of the pole masses with a maximum/minimum ratio of 2, which is a very significant change. Eq. (59) covers the entire region of the secondary peak of the e + e − annihilation cross section above the ρ resonance [38]. In the context of the dual resonance model, Eq. (59) corresponds to varying the pole position from the n = 1 resonance at 3M 2 ρ to values that are half way between this one and the n = 0 or 2 resonances. Note that we let the mass parameters in the delta functions vary independently of each other over the given range, so that the parametrization represents a wide range of "shapes" of the spectral function.
We further constrain the set of mass parameters by requiring that the variation of the spacelike form factor generated by the spectral function be within a certain range around the nominal value. This is essentially a stability condition, which eliminates extreme values of the mass parameters that would lead to large excursions of the spacelike form factor and can be ruled out on physical grounds. We implement this by requiring that (here i = 1, 2) where t ref < 0 is a spacelike t value. In the following applications we choose t ref = −1 GeV 2 and = 0.1 for both F V 1 and F V 2 ; the choice is justified in the following; other choices are possible. The parameter variation Eq. (59), supplemented by the stability condition Eq. (60), generates a functional variation in the high-mass part of the spectral function which we regard as the theoretical uncertainty of our model. Note that the stability condition Eq. (60) restricts the variation of the theoretical FF prediction relative to the nominal value of the model, not relative to an experimental value; no fitting to the FF data is performed here. The parameters t ref and are chosen such that the resulting theoretical model uncertainty of the FFs is reasonable and covers the experimental data. In this way the experimental FF data are used only in estimating the theoretical uncertainty of the model, not in determining the nominal model prediction.
To map out the theoretical uncertainty in practice, we generate a random ensemble of mass parameters in the range of Eq. (59) and retain those for which the spacelike FFs satisfy the condition Eq. (60). We then use this restricted ensemble to generate uncertainty bands in the spectral functions and transverse densities (and possibly other quantities derived from the spectral functions). Figure 4 illustrates the procedure in the case of F V 2 . One observes that the procedure generates natural uncertainty bands, which are approximately symmetric around the nominal value. The resulting uncertainty will be quoted as "high-mass uncertainty" in the results below.
We emphasize that the procedure respects analyticity and the dispersive sum rules. Each instance in the ensemble corresponds to a form factor with correct analyticity in t, and a density with correct asymptotic behavior at large b. Each instance represents a spectral function that satisfies the sum rules Eqs. (44)- (46) and Eqs. (53)-(56) and produces form factors and densities with the correct normalization. The only differences between the instances are in the form of the high-mass spectral function, and in the distribution of strength between the low-mass and high-mass regions.
(II) Uncertainties due to the nucleon radii. The nucleon radii determine the spectral function parameters through the dispersive sum rules Eqs. (44)-(46) and Eqs. (53)-(56). We can estimate the resulting uncertainty by varying the value of the radii. The empirical values of the radii and their uncertainties are summarized in Appendix A. For our uncertainty estimate, we vary the the radii in a range corresponding to their empirical uncertainty We then follow the effect of this variation from the spectral function to the form factors and densities calculated as dispersive integrals. The resulting uncertainty will be quoted as "radius uncertainty" in the results below.  The middle and lower panels of Fig. 5 show the uncertainties in the two-pion part of the spectral functions resulting from the parametrization of the high-mass part and from the nucleon radii, estimated with the procedure of Sec. II D. (Note that the figure shows only the variation of the two-pion part of the spectral function; the high-mass part undergoes a corresponding variation with the parameters, so that the sum rule are satisfied; this part is not shown in the figure.) One observes: (a) The uncertainties of the spectral functions are negligible in the region of t from the threshold at 4M 2 π to ∼ 0.4 GeV 2 . In this region the chiral expansion of the functions N i (t) is well convergent, and the spectral functions represent genuine predictions of the theory. Notice that the enhancement of the spectral functions through ππ rescattering is already very significant in this region [23]. (b) In the region of the ρ resonance, the spectral func-tions show significant uncertainties from the high-mass states and from the radii. The behavior is this region is mainly constrained by the sum rules, so that the spectral functions become sensitive to the parameters, as expected. The relative uncertainties are of order unity and approximately the same in Im F V 1 (t) and Im F V 2 (t). In the present study we use the spectral functions to compute the peripheral densities. The dispersive integrals for the peripheral densities converge rapidly and sample mostly the two-pion part of the spectral functions; the contribution of high-mass states is strongly suppressed (see Fig. 3). Our DIχEFT method and uncertainty estimates aim to provide a realistic description of the two-pion part, while parametrizing the high-mass part in summary form. The dispersive integrals for the spacelike FFs converge more slowly and are more sensitive to the high-mass states. The computation of spacelike FFs therefore generally places stronger demands on the description of the high-mass states than are needed in the present study. Still, it is instructive to see how our simple spectral functions perform in the computation of the spacelike form factors, for which experimental data are available. Figure 6 shows the spacelike FFs F 1 (t) and dictions with the nominal values of the high-mass pole position and the nucleon radii. One observes that the nominal predictions agree very well with the empirical FFs extracted from experimental data [41]. This indicates that our assumptions made in parametrizing the high-mass part of the spectral function (in particular, the rapid saturation at low masses t ≈ 3M 2 ρ ) are realistic at the quantitative level. The middle and lower panels show the uncertainties of the predictions due to the position of the high-mass poles and the values of the nucleon radii, estimated with the procedure of Sec. II D. One observes that the procedure gives a reasonable uncertainty estimate of the spacelike FFs that is approximately symmetric around the nominal value (by design of the procedure) and covers the experimental values. We emphasize that our goal here is not to predict or analyze the spacelike FFs, but just to validate that our spectral function results are compatible with the spacelike FF data. It is clear that a much more accurate description of the spacelike FFs could be achieved within our framework if the value of the high-mass poles were used as fit parameters, as is commonly done in dispersive fits [39,40].  Figure 7 shows the densities ρ V 1 (b) and − ρ V 2 (b) on a logarithmic scale and their relative uncertainties. Figure 8 shows the radial densities 2πbρ V 1 (b) and −2πb ρ V 2 (b) on a linear scale and their absolute uncertainties.

B. Isovector densities
One observes: (a) The densities exhibit an exponential decrease at b 0.5 fm, as dictated by the analytic properties of the form factor (see Sec. II B). This behavior is naturally obtained from the dispersive representation and is the principal reason for the use of this method for the computation of peripheral densities. (b) The uncertainties of the densities resulting from the high-mass part of the spectral function and from the nucleon radii show a characteristic dependence on b (nodes, maxima). This dependence is explained by the way in which these parameters influence the low-mass and high-mass parts of the spectral functions through the dispersive sum rules (see Sec. II C). (c) The uncertainty bands are bounded at large distances and permit stable estimates of the uncertainties of the peripheral densities. In both ρ V 1 (b) and ρ V 2 (b), the estimated relative uncertainties from highmass states and radii are 10% at distances b > 0.5 fm.
(d) At b > 0.5 fm the relative uncertainties of ρ V 1 and ρ V 2 are comparable. At b < 0.5 fm, the relative uncertainty of ρ V 2 is larger than that of ρ V 1 . This happens because at small b the dispersion integral samples the high-mass part of the spectral function, and our parametrization of Im F V 2 allows for more variation than that of Im F V 1 . We emphasize that our theoretical calculations are aimed only at the "peripheral" densities, with the boundary in b determined by the uncertainty estimates. Our results for both the densities and uncertainty estimates are to be understood in this sense. At distances b 0.3 fm the densities become sensitive to the details of the high-mass states in spectral function. Our effective description is not expected to be accurate in this region but still allows us to estimate the uncertainty and demonstrate its increase. Figures 7 and 8 also show the empirical densities, computed as Fourier transforms of the FF fit of Ref. [41], see Eqs. (4) and (8). Because it is difficult to quantify the uncertainties of the Fourier transform densities at large b (see discussion in Sec. III A), we show only their central value. One observes that our theoretical results show excellent agreement with the empirical densities at distances b 0.5 fm, where our calculations are accurate according to our intrinsic uncertainty estimates. Our result obtained with the nominal parameters follow the empirical densities even down to b → 0. This strongly indicates that our assumptions made in parametrizing the highmass part of the spectral function are realistic regarding the nominal values, and that the uncertainty of our predictions at small b could be reduced by constraining the variation of the high-mass part through further theoretical arguments or empirical information.

C. Proton and neutron densities
The DIχEFT method allows us to predict the isovector densities, which contain the effect of the two-pion states and dominate peripheral nucleon structure. In order to compute the individual proton and neutron densities in the dispersive representation we need also the isoscalar spectral functions. For the purposes of this study we use the parametrization of Appendix C, which describes the isoscalar spectral function through effective poles and implements the dispersive sum rules in a similar manner as in the isovector sector. The parameters are fixed in terms of the isoscalar nucleon radii and their uncertainties. The proton and neutron densities and their uncertainties are obtained by combining the isovector and isoscalar densities. As the high-mass uncertainty of the proton and neutron densities we take only the one resulting from the isovector densities estimated in Sec. III B, which is expected to be dominant. As the radius uncertainty of the proton and neutron densities we quote the change of the density under variation of the nucleon's "own" radius, i.e., r 2 p 1,2 for the proton and r 2 n 1,2 for the neutron, corresponding to a simultaneous change of the isovector and isoscalar radii. This is the dominant radius uncertainty in most of the densities. The only exception is the neutron charge density, where the uncertainty resulting from the change of the neutron radius is smaller than that resulting from the proton radius, through its combined effect on the isovector and isoscalar densities. Figures 9 and 10 show the densities ρ p,n 1 (b) and ρ p,n 2 (b) densities and their uncertainties obtained in this way. One observes: (a) The calculation predicts the peripheral nucleon densities with good accuracy. In the proton charge density ρ p 1 , and the proton and neutron magnetization densities ρ p,n 2 , the relative uncertainties are estimated at 10% at b > 0.5 fm (these densities are uniformly positive or negative, so that one can sensibly quote the relative uncertainty). (b) In the neutron charge density ρ n 1 , the absolute uncertainty is estimated at > 50% near the positive maximum at b = 0.65 fm, and rapidly decreasing at larger b (this density has different sign in different regions). (c) The nominal DIχEFT predictions agree well with the empirical densities at b 0.5 fm. In particular, the theoretical calculation reproduces the behavior of the neutron density, which changes from negative values at b 1.5 fm to positive values at 1.5 b 0.35 fm to negative values at b 0.35 fm [4,28]. This behavior arises as the result of a delicate cancellation of isovector and isoscalar densities in the different regions of b.

D. Current densities in polarized nucleon
Combining our results for the densities ρ 1 (b) and ρ 2 (b), we can compute the J + current density in the transversely polarized nucleon, Eq. (7). In order to display the theoretical uncertainties it is useful to show a onedimensional projection of the two-dimensional current density. We consider the current density Eq. (7) in the nucleon with S y = +1/2 on the transverse x axis, where b = (b x , 0) with b x < 0 or > 0, which is given by This function describes the current density to the "left" and "right" when looking at the nucleon from z = +∞ (see Fig. 1). Notwithstanding its piecewise definition in Eq. (63), it is a smooth function of b x because ρ 2 (|b x | = 0) = 0. Figure 11 shows the J + current density Eq. shows that the internal motion due to the nucleon spin causes a significant distortion of the plus momentum distribution and attests to the essentially relativistic character of the system. A "mechanical" interpretation of the left-right asymmetry of the peripheral densities in traditional chiral EFT, as arising from the motion of a soft pion in the nucleon's periphery, has been developed in Refs. [8,9,42,43].

IV. DISCUSSION
In the present work we have computed the peripheral transverse charge and magnetization densities in the nucleon using the DIχEFT method and quantified their uncertainties. The main findings are: (a) The dispersive representation permits stable calculation of the peripheral densities. The densities exhibit the exponential decrease implied by analyticity of the form factor and depend smoothly on the parameters of the spectral function. (b) Uncertainties can be estimated by allowing for variation of the spectral function (functional form, parameters) and following its effect on the densities. The procedure makes use of the particular "information flow" implied by analyticity and relates the peripheral densities to the spacelike form factor in a controlled man-ner. (c) Using a minimal parametrization of the highmass part of the isovector spectral function, the isovector densities are computed with an estimated accuracy of ∼ ±10% at b 0.3 fm.
In the present calculation we have not used any spacelike form factor data beyond the nucleon radii (FF derivatives at t = 0) to constrain the isovector spectral functions. In particular, we do not fit the high-mass part of the spectral function to the spacelike form factor data, as is done is dispersive fits. [The stability condition Eq. (60), controlling the variation of the high-mass spectral function, applies to the variation relative to the nominal theoretical prediction, not relative to the FF data.] Our results represent theoretical predictions based on a minimal parametrization of the high-mass spectral function, and our uncertainty estimates should be understood in this sense. It is clear that a much more accurate description could be achieved if spacelike FF data were used to constrain the high-mass part of the isovector spectral functions. Our estimates of the high-mass uncertainty therefore should not be regarded as "final," but rather as showing how far one can go without fitting spacelike FF data.
The methods developed here enable an EFT-based computation of the transverse densities down to distances b 0.5 fm. At such distances the transverse densities can be described in approaches using other effective degrees of freedom, e.g. quark models. This makes it possible to match the EFT-based description with quark model predictions of the transverse densities at "intermediate" distances and explore quark-hadron duality in new ways. While quark models may not be able to accurately reproduce the absolute densities, they can predict qualitative features such as the spin/isospin dependence and flavor decomposition of the densities, which can lead to interesting conclusions when matched with the EFT description. We plan to explore the use of transverse densities for quark-hadron duality studies in a separate work.
Some comments are in order regarding the interpretation of our results in terms of a "pion cloud" of the nucleon. It is true that the isovector densities at distances 0.5 fm are generated mostly by the two-pion states in the dispersive representation (see Fig. 3). One might be tempted to explain this in a picture where a bare nucleon fluctuates into a pion-nucleon state, and the peripheral structure arises from the propagation of that pion. Such a picture is indeed obtained in traditional chiral EFT, where the pion and nucleon are pointlike, and the peripheral densities emerge from the propagation of the pointlike pions. However, in our unitarity-based approach the pion is not pointlike and has an extended structure of the same range as the nucleon, as required by unitarity. The results of our approach should therefore not be interpreted in terms of the traditional pion could picture. The space-time interpretation of the densities in the unitarity-based approach is an interesting question which we plan to investigate in a future study.
In the present study we have applied our unitaritybased approach to the peripheral densities associated with the nucleon matrix element of the electromagnetic current operator. The approach could be extended to compute the peripheral densities of other operators whose form factors possess a two-pion cut, such as the QCD energy momentum tensor (spin-2 operator) or the leading-twist QCD operators whose matrix elements determine the moments of the GPDs (twist-2, spin-n operators, n ≥ 1). This would allow one to "deconstruct" not only the nucleon's electromagnetic current but also its peripheral partonic structure in terms of EFT degrees of freedom. One difference between the electromagnetic and the generalized form factors is that for the latter the "radii" (derivatives at t = 0) are generally not known from independent measurements, so that one has to adjust the procedure of fixing the parameters of the spectral functions and recruit new sources of information.

Appendix A: Nucleon radii
In this appendix we list the values of the nucleon radii used as parameters in the DIχEFT calculation of the spectral functions. The Dirac and Pauli radii of the pro- ton and neutron are defined in terms of the FF derivatives (same for p → n). They are related to the conventional electric and magnetic radii by (same for p → n). Here κ p,n are the anomalous magnetic moments, and µ p,n = Q p,n + κ p,n = (2.793, −1.913) are the ordinary magnetic moments of the nucleons. We estimate the values of r 2 1,2 and their uncertainties from the empirical values of r 2 E,M and their uncertainties, neglecting correlations between the uncertainties of r 2 E and r 2 M . We use the following numbers and sources: r 2 p E = (0.7090 ± 0.0168) fm 2 [22], r 2 n E = (−0.1161 ± 0.0022) fm 2 [44], r 2 p M = (0.7225 ± 0.0170) fm 2 [22], r 2 n M = (0.7465 ± 0.0156) fm 2 [44]. The proton and neutron radii r 2 1,2 thus obtained are summarized in Table I. We also list in Table I the isovector and isoscalar combinations of the radii, defined by Eqs. (47)-(48) and (57)-(58), and Eqs. (C6)-(C7) and (C8)-(C9), which enter in the sum rules for the isovector and isoscalar spectral functions. In calculating the uncertainties we neglect correlations between the uncertainties of the proton and neutron radii. Note that the isovector/isoscalar combinations of r 2 2 defined by Eqs. (57)-(58) and (C8)-(C9) involve the nucleon anomalous magnetic moments and cannot directly be interpreted as nucleon radii.

Appendix B: N functions
In this appendix we give the expressions for the N i (t) functions appearing in the unitarity relations for the isovector spectral functions, Im F V 1,2 (t) Eqs. (36) and (38). We do not compute these functions explicitly but express them in terms of our earlier results for the J 1 ± (t) functions appearing in the unitarity relations for the Im G V E,M (t) spectral functions [18]. The Dirac/Pauli FFs F 1,2 (t) and the electric/magnetic FFs G E,M (t) are related by or, inversely, which hold for any complex t. The elastic unitarity relation for Im G E,M (t), in its manifestly real form analogous to Eq. (36), is written as (here t > 4M 2 π ) where k cm is given in Eq. (37). The relation between the functions N 1,2 (t) and J 1 ± (t) is (B8) Note that where p cm = m 2 N − t/4 is the unphysical nucleon center-of-mass momentum in the ππ → NN process in the t-channel.
The explicit expressions for the J 1 ± (t) functions in chiral EFT at LO and NLO accuracy are given in Ref. [18]. The expressions for the functions N 1,2 (t) at this order can be obtained from those results using Eqs. (B7) and (B8). Note that in this procedure Eqs. (B7) and (B8) are evaluated only in the region of elastic unitarity in the two-pion channel, t < t max ≈ 1 GeV 2 , far away from the singularity at t = 4m 2 N . In our calculation at pN2LO accuracy we approximate the full N2LO corrections to the N 1,2 (t) functions by rescaling the tree-level N2LO result, see Eq. (40). For the N2LO tree-level result we use the simplified form the overall coefficients are irrelevant since they are absorbed by the parameters λ i in Eq. (40). Equation (B10) differs from the exact N2LO tree-level result [as obtained from Eqs. (B7) and (B8) and the N2LO tree level result for the J 1 ± (t)] only by additive terms ∝ M 2 π with coefficient of order unity, which are numerically small at t = 20-50 M 2 π where the pN2LO term Eq. (40) is significant. The simplified form Eq. (B10) ensures that each of the functions N 1,2 (t) gets an pN2LO term with its own parameter λ 1,2 , which can then be fixed through the dispersive sum rules.

Appendix C: Isoscalar parametrization
In this appendix we present a simple parametrization of the isoscalar spectral functions, which is used in the dispersive calculations of the individual proton and neutron densities. The isoscalar spectral functions are constructed along similar lines as the isovector ones in Sec. II C, but relying more on empirical information.
The isoscalar spectral function starts with the 3π channel and can be organized in a similar way as Eq. (34) (here i = 1, 2), The strength in the 3π channel is overwhelmingly concentrated in the ω resonance at t = M 2 ω = 0.61 GeV 2 ; for an estimate of non-resonant 3π contributions in chiral EFT, see Ref. [45]. We parametrize the 3π part of the spectral function as (i = 1, 2) High-mass strength appears at t 1 GeV 2 through the KK channel and the φ resonance, as well as through other hadronic states such as πρ [46,47]. We assume that the high-mass part of the isoscalar spectral function is approximately exhausted by these states at t ∼ 1 GeV 2 and parametrize the spectral functions as  Table I. Dashed lines: Empirical form factors of Ref. [41]. All FFs are shown divided by the standard dipole FF.