Precise determination of proton magnetic radius from electron scattering data

We extract the proton magnetic radius from the high-precision electron-proton elastic scattering cross section data. Our theoretical framework combines dispersion analysis and chiral effective field theory and implements the dynamics governing the shape of the low-$Q^2$ form factors. It allows us to use data up to $Q^2\sim$ 0.5 GeV$^2$ for constraining the radii and overcomes the difficulties of empirical fits and $Q^2 \rightarrow 0$ extrapolation. We obtain a magnetic radius $r_M^p$ = 0.850 $\pm$0.001 (fit 68%) $\pm$0.010 (theory full range) fm, significantly different from earlier results obtained from the same data, and close to the extracted electric radius $r_E^p$ = 0.842 $\pm$0.002 (fit) $\pm$0.010 (theory) fm.


INTRODUCTION
The electromagnetic form factors (EM FFs) are the most basic expressions of the nucleon's finite spatial extent and composite internal structure. They describe the elastic response to external electric and magnetic fields as a function of the 4-momentum transfer Q 2 and can be associated with the spatial distributions of charge and current in nucleon. The traditional representation of FFs in terms of 3-dimensional spatial densities at fixed instant time x 0 = const. is appropriate only for nonrelativistic systems such as nuclei [1]. For relativistic systems such as hadrons, the spatial structure is expressed through 2dimensional transverse densities at fixed light-front time x + = x 0 + x 3 = const. In the context of QCD these transverse densities can be regarded as projections of the nucleon's partonic structure (generalized parton distributions) [2][3][4]. The EM FFs thus reveal aspects of the spatial distribution of quarks and their orbital motion and spin and have become objects of great interest in nucleon structure studies [5,6].
The value of the electric and magnetic proton FFs at Q 2 = 0 is given by the total charge and magnetic moment of the proton, G p E (0) = 1, G p M (0) = µ p = 2.793. The leading information about the spatial structure is in the first derivatives of the FFs at Q 2 = 0. They are conventionally expressed in terms of the equivalent electric and magnetic 3-dimensional root-mean-square radii, (1) this does not imply an actual physical interpretation in terms of 3-dimensional densities; the proper interpretation in terms of 2-dimensional densities is discussed below [1]. Besides their importance for nucleon structure, the FF derivatives are needed in tests of atomic bound-state calculations in quantum electrodynamics and in precision measurements of the Rydberg constant [7,8]. The proton electric (or charge) radius is extracted from the proton FFs measured in electron-proton elastic scattering, and from the nuclear corrections to atomic energy levels (electronic and muonic hydrogen) measured in precision spectroscopy experiments; see Refs. [9,10] for a review. Apparent discrepancies between the different extraction methods ("proton radius puzzle") have engendered intense experimental and theoretical efforts, including dedicated new FF measurements at low Q 2 using electron and muon beams [11,12]. Recent results seem to converge around r p E = 0.84 fm [11,[13][14][15]. The magnetic radius can be extracted only from elastic scattering measurements. Recent determinations based on the Mainz A1 data [16,17], using methods developed in the context of the charge radius extraction, have resulted in a range of values that disagree with each other, r p M = 0.78(2) fm [16], 0.914(35)fm and 0.776(38)fm [18], and significantly depart from older results ∼0.85 fm [19]. It is necessary to resolve these discrepancies and determine the proton magnetic radius with an overall accuracy and consistency commensurate with those achieved in the electric radius.
Here we report an extraction of the proton magnetic radius from electron scattering data using a novel theoretical framework based on dispersion analysis and chiral effective field theory (DIχEFT) [20][21][22][23]. It implements analyticity and the dynamics governing the shape of the low-Q 2 FFs and allows us to use data up to Q 2 ∼ 0.5 GeV 2 for constraining the radii, increasing the sensitivity to the magnetic FF. It overcomes the difficulties in extraction methods based on empirical fits and Q 2 → 0 extrapolation (functional form bias, unstable extrapolation), particularly the issues related to the normalization of data sets taken at different incident energies. DIχEFT was used in Ref. [20] to extract r p E from an empirical FF parameterization [24] and delivered a value of r p E = 0.844(7) fm, as accepted in the CODATA 2018 update and confirmed by more recent measurements [11,14,15]. In this work we use the method to extract both r p M and r p E from a direct analysis of the cross section data, dominated by the Mainz A1 data [16,17]. We obtain r p M = 0.850 ±0.001 (fit 68%) ±0.01 (theory full range) fm, significantly different from the values extracted from the same data using other methods [16,18], and surprisingly close to r p E . In the course we also improve our extraction of r p E and verify the robustness of the results. DIχEFT is a method for calculating nucleon FFs combining dispersion analysis and chiral effective field theory. The theoretical foundations are described in detail in Refs. [21][22][23]; applications to FF fits are discussed in Ref. [20]. The FFs are represented as dispersion integrals over t ≡ −Q 2 . The spectral functions on the two-pion cut at t > 4M 2 π are calculated using (i) the elastic unitarity relation; (ii) πN amplitudes computed in χEFT at leading order, next-to-leading order, and partial nextto-next-to-leading order accuracy; (iii) the timelike pion FF measured in e + e − annihilation experiments. The approach includes ππ rescattering effects and the ρ resonance and generates accurate spectral functions up to t ∼ 1 GeV 2 . Higher-mass t-channel states are described by effective poles. The parameters specifying the dynamical input (the low-energy constants of the χEFT calculation, and the strength of the effective poles) are related by the sum rules of dispersion theory and can be expressed in terms of the nucleon charges, magnetic moments, and radii. For each (assumed) value of r p E and r p M the theory thus generates a unique prediction for G p E (Q 2 ) and G p M (Q 2 ) with controlled theoretical uncertainties; see Ref. [20] for a summary plot. It predicts the "shape" of the spacelike FF as determined by analyticity (position of singularities) and dynamics (strength of singularities). In this way the values of the radii are correlated with the predicted behavior of the FFs at finite momentum transfers Q 2 ∼ 1 GeV 2 , allowing the use of such data for radius extraction. A computer code generating the DIχEFT FF predictions and further information are available [25].
For our radius extraction we use the high-precision data in electron-proton elastic scattering from the Mainz A1 experiment, which dominate the world data [16,17]. The experiment measured the elastic scattering cross section at momentum transfers 0.003 Q 2 1 GeV 2 and incident electron energies E Beam from 0.18 to 0.855 GeV. The 2-dimensional data cover the cross section at a fixed Q 2 at various values of the virtual photon polarization parameter, , and allow for separation of the contributions of G p E (Q 2 ) and G p M (Q 2 ) through global fits, generalizing the traditional Rosenbluth method [16,17].
An important issue in global fits is the normalization of the data sets taken at different energies. The normalization of the cross section data (both absolute and relative, between different energies) is limited by the knowledge of the absolute luminosity in the different settings and subject to considerable uncertainties. The combination of data taken at different energies therefore requires rescaling of the data sets, which depends on the functional form of the FFs or on other assumptions. In the context of empirical fits the effect of the rescaling on the random uncertainties of the data was studied in detail in Refs. [17,18]. In the context of our approach this  [16,17], with the normalization of sets determined by our fit. Band: Theoretical model (DIχEFT) with parameters (r p E , r p M ) obtained from our best fit. The band shows the range of the model predictions obtained by varying the parameters in the 68% confidence interval of the fit; it does not include the intrinsic theoretical uncertainty of the model [20]. Both data and model are divided by the cross section evaluated with the dipole FFs (Λ 2 = 0.71 GeV 2 ). problem is naturally solved by the fact that the theoretical model predicts the shape of the FFs at finite Q 2 (in dependence of the radii). We can therefore perform a global fit with floating normalizations of the data sets, which can adjust themselves to the theoretical model; the physical information is in the variation of the data with Q 2 , which tests the theoretical predictions for the shape and fixes the radius parameters through the best fit.
As the figure-of-merit for the global fit with floating normalizations we use a χ 2 function of the form The summation is over the N dat data points labeled by i. σ thy,i is the theoretical electron-proton elastic scattering cross section at the kinematic point (E i , Q 2 i ), evaluated with the DIχEFT FFs G p E and G p M with the parameters (r p E , r p M ) (the expression of the elastic scattering cross section in terms of the FFs is given in Ref. [17]). σ exp,i is the measured cross section and ∆σ exp,i is the random uncertainty. The data points are grouped in N set sets measured under the same running conditions; the normalization is assumed to be constant inside each set, but its value is unknown. The N set parameters (Λ 1 , ..., Λ Nset ) represent the floating normalizations in each set; k(i) denotes the index k of the set to which data point i belongs. (A detailed discussion of how the experimental normalizations were defined and obtained can be found in Ref. [26].) The χ 2 defined by Eq. (2) is thus a function of the theory parameters (r p E , r p M ) and the normalization parameters (Λ 1 , ..., Λ Nset ). Minimization is performed with respect to all the parameters simultaneously. The values of the Λ k (k = 1...N set ) at the minimum are found to be equal to unity within 1%; this indicates that the normalization determined in the original analysis of Ref. [17] is reproduced reasonably by our fit; the values themselves have no physical significance otherwise (nuisance parameters). The values of (r p E , r p M ) at the minimum correspond to the best fit to the data and represent the proton radii extracted with our method.
To estimate the uncertainties of the extracted radii, we use the criterion ∆χ 2 = 2.7 to determine the 68% confidence interval, corresponding to the simultaneous estimation of two independent parameters. The uncertainties of the physical parameters r p E and r p M are affected also by the statistical fluctuations of the normalization parameters Λ k ; we have estimated the total statistical uncertainties using a bootstrap method and found them to be very close to the uncertainties of r p E and r p M one would obtain from the variation of the reduced χ 2 , obtained after minimization with respect to (Λ 1 , ..., Λ Nset ). In the final uncertainty we also include the theoretical uncertainty of the DIχEFT FFs (see below) [20].

RESULTS
In the fit we include the cross section data up to a maximum momentum transfer, Q 2 < Q 2 max . Suitable values of Q 2 max are determined by considering the balance of experimental and theoretical uncertainties and the sensitivity of the cross sections to the model parameters r p E and r p M [20]. Our standard fit uses Q 2 max = 0.5 GeV 2 and includes 1285 of the 1422 Mainz A1 data points. The overall quality of the description of the experimental cross sections is shown in Fig. 1 (the plots show also the data at Q 2 > Q 2 max , which were not included in the fit). One sees that all features of the kinematic dependence of the data (with the floating normalization determined by the fit) are reproduced by the theoretical model. The reduced χ 2 profile in the physical parameters r p E and r p M , obtained after minimization with respect to the normalization parameters, is shown in Fig. 2. One observes that the variations of χ 2 in r p E and r p M are approximately independent, and that clear minima are obtained in both parameters. Minimizing with respect to the radii, we extract r p E = 0.842 ± 0.002 fm and r p M = 0.850 ± 0.001 fm with a reduced χ 2 of 1.39.
To test the robustness of the results we have performed fits with different values of Q 2 max and found little effect on the extracted radii. In fact, using the entire Mainz A1 data set up to Q 2 max = 1 GeV 2 gives r p E = 0.843 ± 0.002 fm and r p M = 0.850 ± 0.001 fm with a reduced χ 2 of 1.43. This shows that the theoretical model (evaluated with the "best" value of the radii) accurately describes the Q 2 dependence of the data over the entire range considered here.
As a further test we have performed fits to the rebinned version of the Mainz A1 data of Ref. [18], where the original data sets are rescaled to a common normalization using empirical functional forms, including the effects on the random uncertainties. The fit with Q 2 max = 0.5 GeV 2 uses 569 of the 658 rebinned data points and gives radii r p E = 0.840 ± 0.002 fm and r p M = 0.849 ± 0.001 fm with a reduced χ 2 of 1.07, in good agreement with our fit to the original data. Extending Q 2 max to include all the rebinned data we obtain r p E = 0.841 ± 0.003 fm and r p M = 0.849±0.001 fm with a reduced χ 2 of 1.10, showing similar stability the fit to the original data. Overall, the tests show that the extracted radii are not sensitive to the choice of data sets used in the fits.
The Jefferson Lab PRad experiment has reported a new measurement of the electron-proton elastic cross section down to Q 2 ∼ 10 −4 GeV 2 , significantly extending the reach of earlier measurements [11]. We have performed a fit including the PRad data in addition to the Mainz A1 data and found no change in the extracted r p E and r p M within uncertainties. This happens because the DIχEFT model naturally describes the Q 2 dependence of the low-Q 2 data, with the same value of r p E as favored by the higher-Q 2 data [20]; this was also observed in the analysis of Ref. [27]. Note that the low-Q 2 data are sensitive mostly to r p E , and that our present extraction of r p M requires us to include data up to Q 2 ∼ 0.5 GeV 2 .
In our assessment of the errors of the extracted radii we must include also the intrinsic theoretical uncertainty of the DIχEFT model. This refers to the uncertainty in the predictions for G p E and G p M for given values of r p E and r p M , which results from the modeling of the highmass states in the dispersion integral (effective poles) and was estimated in Ref. [20]. Performing fits with different values of the effective pole mass we estimate the effect on the extracted radii as ∼ ±0.010 fm for both r p E and r p M ; the interval should be regarded as the "full plausible range" of the theoretical uncertainty. Our final results for the extracted radii are thus r p M = 0.850 ±0.001 (fit 68%) ±0.010 (theory full range) fm and r p E = 0.842 ±0.002 (fit 68%) ±0.010 (theory full range) fm.

DISCUSSION
Several aspects of our method and results merit further discussion. In our theory-based extraction method the main impact on the radii comes from the data at "higher" Q 2 ∼ 0.1-0.5 GeV 2 (see Fig. 2 and Ref. [20]). One observes that the magnetic radius is actually better determined than the electric one (see Fig. 2), because the cross section at the "higher" Q 2 is dominated by the contribution of the magnetic FF. In contrast, in methods based on the Q 2 → 0 extrapolation the cross section is always dominated by the electric FF, rendering extraction of the magnetic radius extremely difficult. Our method therefore offers principal advantages for the analysis of the proton's magnetic structure.
The results of our analysis validate previous results for the proton magnetic radius ∼0.85 fm, obtained using dispersive fits of the earlier world data; see Ref. [19] and references therein. They disagree with the results obtained from various empirical fits of the Mainz A1 data [16,18]. This indicates that the observed discrepancies are due to the extraction methods (analyticity, correlations be-tween Q 2 regions from dispersion relations) rather than the different data sets.
The values of the electric and magnetic radii extracted from the data are very close. While this may be accidental, it is qualitatively consistent with the nonrelativistic quark model picture (independent particle motion in an L = 0 orbital, no spin-orbit interactions). Using our method we can also determine the proton's transverse charge and magnetization radii, which refer to the relativistic representation of the FFs in terms of transverse densities and can be related to the generalized parton distributions [1][2][3][4]. The derivatives at Q 2 = 0 of the Dirac and Pauli FFs, F p 1 and F p 2 , are related to those of the electric and magnetic FFs by [κ p ≡ F p 2 (0) = µ p − 1 is the anomalous magnetic moment, m is the proton mass] Equations (1) and (4)-(6) linearly relate b 2 p 1 and b 2 p 2 to (r p E ) 2 and (r p M ) 2 . Using the results of our fit, we obtain b 2 1 = 0.394±0.002 (fit 68%) ±0.011 (theory full range) fm 2 and b 2 2 = 0.531 ± 0.002 (fit) ±0.019 (theory) fm 2 . It is interesting to note that, if one neglected the small difference between the extracted electric and magnetic radii and set (r p E ) 2 = (r p M ) 2 , the transverse charge and magnetization radii would be related as i.e., the difference would be entirely proportional to the proton magnetic moment. Equation (7) is the partonic expression of the approximate equality of the electric and magnetic radii. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contract DE-AC05-06OR23177. J.M.A. acknowledges support from the Community of Madrid through the Programa de atracción de talento investigador 2017 (Modalidad 1), the Spanish MECD grants FPA2016-77313-P.