Proton charge radius extraction from electron scattering data using dispersively improved chiral effective field theory

We extract the proton charge radius from the elastic form factor (FF) data using a novel theoretical framework combining chiral effective field theory and dispersion analysis. Complex analyticity in the momentum transfer correlates the behavior of the spacelike FF at finite $Q^2$ with the derivative at $Q^2 = 0$. The FF calculated in the predictive theory contains the radius as a free parameter. We determine its value by comparing the predictions with a descriptive global fit of the spacelike FF data, taking into account the theoretical and experimental uncertainties. Our method allows us to use the finite-$Q^2$ FF data for constraining the radius (up to $Q^2\sim$ 0.5 GeV$^2$ and larger) and avoids the difficulties arising in methods relying on the $Q^2 \rightarrow 0$ extrapolation. We obtain a radius of 0.844(7) fm, consistent with the high-precision muonic hydrogen results.


I. INTRODUCTION
The proton charge radius is a fundamental quantity of nuclear physics and attests to the hadron's finite spatial extent and composite internal structure. It is defined as the derivative of the proton electric form factor (FF) at zero momentum transfer, (r p E ) 2 ≡ −6 dG p E /dQ 2 (Q 2 = 0), and describes the leading finite-size effect in the interaction with long-wavelength electric fields; see Ref. [1] for a critical discussion of its interpretation. The electric and magnetic FFs at Q 2 > 0 are measured in elastic electron-proton scattering experiments; see Refs. [2,3] for a review. The radius is also extracted from nuclear corrections to atomic energy levels measured in precision spectroscopy experiments. Measurements of muonic hydrogen transitions have obtained a value r p E = 0.84087(39) fm [4,5], significantly smaller than the value of 0.875 fm by the Committee on Data for Science and Technology (CODATA), obtained from electronic hydrogen transitions and some information from electron scattering data [6]. The discrepancy, known as the "proton radius puzzle," is the subject of a lively debate and has stimulated extensive theoretical and experimental research [7,8], including new dedicated low-Q 2 electronproton and muon-proton scattering experiments [9,10].
Determining the charge radius from electron scattering data amounts to inferring the derivative of the FF at Q 2 = 0 from the data at finite Q 2 . From an empirical point of view, the problem presents itself as one of "extrapolation" of the measured FF to Q 2 → 0. Two approaches have been taken in most studies so far; see Ref. [11] for a review and [12] for the general concepts. Descriptive fits (e.g. higher-order polynomial fits) provide excellent descriptions of the data over a wide range of Q 2 , but the functions are generally not well-behaved outside the fitted region [13][14][15]. Predictive models (e.g. fits with low-order polynomials or other smoothly varying functions) permit stable extrapolation but are constrained by either the selected functional form or tightly bounded parameters [16][17][18][19][20]. In both approaches the question arises over what Q 2 range the extrapolation should optimally be performed, and what uncertainties are associated with this choice.
Complex analyticity plays an essential role in the behavior of the proton FF at low Q 2 . The FF is an analytic function of Q 2 , with singularities at Q 2 < 0, starting with the two-pion cut at Q 2 < −4M 2 π . The behavior of the FF at Q 2 > 0, where it is measured in elastic scattering, is governed by the position of these singularities and by their strength (spectral function), which can be calculated using theoretical methods. Analyticity thus implies correlations between the behavior of the FF in different regions of the Q 2 > 0 domain, which are not apparent in purely descriptive fits. It should therefore inform the analysis of low-Q 2 FF data and extraction of the radius [21][22][23]. In this way one can go beyond the method of extrapolation and use data in a wider Q 2 domain to constrain the derivative at zero.
Here we report an extraction of the proton charge radius using a novel predictive theoretical framework that implements analyticity -dispersively improved chiral effective field theory (DIχEFT) [24][25][26]. We express the spacelike proton FF predicted by the theory in a form such it contains the radius as a free parameter, which directly exhibits the correlation between the finite-Q 2 behavior and the derivative at Q 2 = 0 implied by analyticity. We determine the value of the radius by comparing the theoretical predictions with a descriptive global fit of the spacelike FF data [27] over a broad range of Q 2 (optimally up to ∼ 0.5 GeV 2 ), taking into account the theoretical and experimental uncertainties. At the "best" radius the theory describes the data with the same accuracy as the global fit. Our approach thus combines the best features of descriptive and predictive modeling. It recruits the finite-Q 2 FF data for constraining the radius and overcomes the theoretical and experimental limitations of the Q 2 → 0 extrapolation. We obtain a radius of 0.844(7) fm, which reconciles the electron scattering data with the muonic hydrogen value. We comment on possible improvements of the radius extraction and the arXiv:1809.06373v2 [hep-ph] 7 Apr 2019 relation to other approaches.

A. Predictive theoretical framework
DIχEFT is a method for calculating the nucleon FFs combining chiral effective field theory (χEFT) -a systematic description of strong interactions at distances O(M −1 π ); and dispersion analysis -the use of complex analyticity for connecting the behavior of hadronic amplitudes in different kinematic regions. The method is described in detail in Refs. [24][25][26]; the essential elements used in the present calculation are summarized for reference in Appendix A. The FFs are represented as dispersive 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 next-to-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 allows one to calculate the twopion spectral functions up to t ≈ 1 GeV 2 . Higher-mass t-channel states are described by effective poles, whose strength is fixed by the dispersive integrals for the nucleon charges, magnetic moments, and radii (sum rules) [26]. The nucleon radii thus enter as explicit parameters in the DIχEFT predictions of the spectral functions. Evaluating the finite-Q 2 dispersive integrals with these spectral functions, we obtain an analytic parametrization of the spacelike FFs in which the nucleon radii appear as explicit parameters [see Appendix A, specifically Eqs. (A15) and (A16)]; all other dynamical input is determined independently by the πN scattering data and the pion timelike FF data. We emphasize that in this approach the correlation between the radius and the finite-Q 2 behavior appears through the global analytic properties of the FF (dispersive representation, sum rules), not through a power series expansion at Q 2 = 0; the correlation therefore extends beyond the range of convergence of the power series expansion (the connection with the power series expansion is discussed further in Sec. IV B).
In the application here we take the proton charge radius as a free parameter, to be varied over a range covering the presently discussed values. The neutron charge radius, which enters indirectly through the separate sum rules for the isovector and isoscalar FFs, is fixed at its Particle Data Group (PDG) value [28]; its influence on the proton FF is negligible in the Q 2 region considered here (see Appendix A, specifically Fig. 5). In this way DIχEFT provides us with a family of theoretical predictions of G p E (Q 2 ), with each function respecting analyticity in Q 2 and corresponding to definite value of the proton charge radius. Figure 1 shows the predictions for a set of radii r p E = (0.80, 0.82, 0.84, 0.86, 0.88, 0.90) fm. The dominant uncertainties in the DIχEFT predic-tions of G p E (Q 2 ) (for a given r p E ) arise from the the effective description of high-mass states in the isovector spectral function. We have estimated them conservatively, by varying the position of the effective pole over a range M 2 1 = (0.5-2)×M 2 1 (nom), where M 2 1 (nom) is the nominal value determined in Ref. [26] [see Appendix A, specifically Eq. (A5)]. The results are shown by the bands in Fig. 1. One sees that the uncertainties of the FF predictions are small at low Q 2 (where the dispersive integral is dominated by the two-pion cut and constrained by the given value of the radius), but increase at larger Q 2 (where the dispersive integral becomes sensitive to high-mass states in the spectral function).

B. Descriptive global fit
In order to confront the DIχEFT predictions with the experimental proton FF we use the results of a descriptive global fit [27]. It employs a bounded polynomial zexpansion [29] and determines G p E and G p M directly from the cross section and polarization data. Sum rules are imposed to ensure the correct normalization at Q 2 = 0 and the asymptotic scaling behavior at large Q 2 [15]. The treatment of uncertainties includes the covariance matrix of the fit itself, the systematic errors arising from the tension between data sets, and the uncertainty from two-photon exchange corrections at high Q 2 . In the original work of Ref. [27] the proton radii in the fit were fixed at their presumed empirical values (r p E = 0.879 fm from CODATA, r p M = 0.851 fm from PDG); the radius constraints were then removed when evaluating the fit uncertainty. In the present study we have used the same fitting method to generate a family of fits with different values of the proton charge radius, covering the range 0.8-0.9 fm (see Fig. 1). The magnetic radius is kept fixed at its PDG value; the uncertainty resulting from this simplification is negligible and covered by the quoted overall fit uncertainty. In these fits we only include the uncertainties from the covariance matrix. One observes that in the global fits there is little correlation between the proton charge radius and the value of the FFs at finite Q 2 (the present data cover the range 0.01 GeV 2 ), as expected for this descriptive approach. predictions provide a description of the data with the same accuracy as the global fit up to Q 2 ∼ 1 GeV 2 (actually up to even larger values Q 2 ∼ 2 GeV 2 , which are not included in Fig. 1).

C. Method for radius extraction
The observations suggest a simple method for extracting the proton radius from the FF data: Compare the DIχEFT predictions with the global fit in the region where both descriptions are valid, and determine the radius by best agreement. The method combines the advantages of the descriptive fit (reliable uncertainty estimates in the region where there are data) and the predictive theory (correlation of the radius and the finite-Q 2 behavior through analyticity).
The method can be optimized by choosing the "best" Q 2 region for the comparison. Figure 2 shows the theoretical uncertainty of the DIχEFT FF prediction for a fixed proton radius, and the experimental uncertainties of the FF obtained from the global fit. It also shows the variation of the DIχEFT FF predictions under a certain change of the proton radius (here, ∆r = ±0.01 fm, ±0.02 fm, etc.), which quantifies the sensitivity of the FF to the proton radius as a function of Q 2 . One observes: (a) At "low" Q 2 (∼0.1 GeV 2 ) the theoretical uncertainty of the DIχEFT FF is much smaller than the experimental uncertainty of the global fit, which is mainly due to normalization errors (inconsistencies between data sets). The variation of the FF under the change ∆r = ±0.01 is larger than the theoretical uncertainty but smaller than the fit uncertainty. (b) At "high" Q 2 (∼1 GeV 2 ) the theoretical uncertainty is larger than the fit uncertainty. The variation of the FF under ∆r = ±0.01 is smaller than the theoretical uncertainty and comparable to the fit uncertainty. Based on Fig. 2 we choose the upper limit of the Q 2 -region for radius extraction as Q 2 max ∼ 0.5 GeV 2 ; with this upper limit the theoretical error remains smaller than, or at most comparable to, the fit error. The lower limit we choose as Q 2 min ∼ 0.01 GeV 2 ; this represents the lowest value for which FF data are presently available, and the results are not sensitive to this choice. (The performance of our method with smaller values of Q 2 max and the relation to low-Q 2 fits using a power series expansion are discussed in Sec. IV B.) . The variation quantifies the sensitivity of the theoretical FF predictions to the proton radius parameter (cf. Fig. 1).
To quantify the agreement of the theoretical model with the global fit and extract the radius, we use a figure of merit in the form of a reduced χ 2 , The sum runs over N bins in Q 2 covering the range (Q 2 min , Q 2 max ); we use N = 50; the results are not sensitive to the binning. thy i denotes the DIχEFT prediction of G p E (Q 2 i ) in the i'th Q 2 bin for the given r p E ; fit i denotes the global fit result for G p E (Q 2 i ) in the same bin. The theoretical and fit uncertainties, ∆thy i and ∆fit i , are added in quadrature. Figure 3 shows the reduced χ 2 as a function of r p E . One observes: (a) The dependence is approximately quadratic, indicating a natural best agreement without tension.

III. RESULTS
Using the above method and Q 2 -range, we have extracted the proton charge radius and its uncertainty, obtaining a value r p E = 0.844 (7) fm. The uncertainty estimate is based on the combined theoretical and global fit uncertainties entering in our figure of merit Eq. (1) and corresponds to a confidence interval with ∆χ 2 ≤ 1. The extracted radius is consistent with the high-precision muonic hydrogen results and clearly disfavors the CO-DATA result. Our result therefore suggests that the finite-Q 2 electron scattering data agree well with the muonic hydrogen results, and that the disagreement is rather between the electronic and muonic hydrogen results. We note that some recent measurements with electronic hydrogen have yielded a value consistent with the muonic result [30], while others agree with the current CODATA value [31].

A. Possible improvements
The proton radius extraction reported here could be improved in several aspects: (a) by reducing the theoretical uncertainty of the DIχEFT FF predictions through an improved description of the high-mass spectral functions (t > 1 GeV 2 ), which would be possible with a more flexible parametrization and further theoretical constraints; (b) by reducing the experimental uncertainties of the FF data, especially in the region Q 2 0.2 GeV 2 , which would allow us to limit the theory-experiment comparison to a smaller Q 2 -interval (Q 2 max ∼ 0.2 GeV 2 ), where the theoretical uncertainties are smaller.
In our analyticity-based framework the main impact on the proton radius comes from FF data at moderate Q 2 (∼ 0.1-0.5 GeV 2 ) rather than at the lowest available Q 2 . The forthcoming FF data at very low Q 2 from the Jefferson Lab PRad experiment [9] (down to a few ×10 −4 GeV 2 ) can complement the results of our study by reducing the normalization errors of the data (cf. the discussion below) and enabling an independent radius extraction using traditional extrapolation methods. They can also validate our theoretical framework, e.g. by extracting higher FF derivatives, which enable sensitive tests of the DIχEFT spectral functions [25].

B. Relation to low-Q 2 FF fits
In Ref. [18] the proton radius was extracted from fits to the low-Q 2 FF data (Q 2 max 0.02 GeV 2 ) with a truncated power series in Q 2 , in which the coefficient of the first-order term in Q 2 (the proton radius) was determined by the data, and the coefficients of the higher-order terms Q 4 , Q 6 , etc. (the higher moments) were calculated in standard χEFT and supplied as a theoretical input. It is worth explaining how our dispersive method relates to that approach when the fit is restricted to the low-Q 2 region. Figure 4 shows the DIχEFT FF parametrization (with the radius determined by our above best fit), as well as its first-order and second-order expansion in Q 2 ; the coefficients are given by the derivatives of the FF, evaluated using the dispersive integral with the DIχEFT spectral functions. One sees that the full DIχEFT FF is well approximated by the second-order expansion, with the second-order term giving a correction of 1% (5%) at Q 2 = 0.03 GeV 2 (0.06 GeV 2 ). The second derivative of the DIχEFT FF changes only by a few percent when the radius is varied within the range r p E = (0.80, 0.90) fm, so that the coefficient of the second-order term may be regarded as a fixed theoretical input. Our method therefore effectively reduces to that of Ref. [18] in the region where the second-order approximation is valid (Q 2 max < 0.05 GeV 2 ). The advantage of our method is that it is not limited by the power series expansion and allows us to include FF data at significantly higher Q 2 into the fit, making the radius extraction more robust. This can be seen in the trend shown in Fig. 3: increasing the value of Q 2 max decreases the uncertainty in the radius, as the fit is constrained by more data, while at the same time the theoretical uncertainties are still under control. For reference we note that, if we restricted our fit to Q 2 max = 0.1 GeV 2 , we would obtain a radius r p E = 0.849(10) fm, which agrees well with the result of Ref. [18] in the central value and the uncertainty.
We also point out that the second derivative (or moment) of the DIχEFT FF is significantly larger than the standard χEFT results used in Ref. [18], because the ππ rescattering effects included in the DIχEFT calculation increase the spectral function on the two-pion cut in the near-threshold region; see the discussion in Ref. [25]. This shows that in the approach of Ref. [18] the theoretical uncertainty would deteriorate quickly if it were applied at Q 2 values where the second-order (and higher-order) terms become sizable. We note that the ππ rescattering effects are noticeable even in the higher FF derivatives (n ≥ 3); nevertheless, for these quantities the DIχEFT results agree with the standard χEFT predictions within uncertainties; see Ref. [25] for details. Explicit expressions for the FF derivatives in standard χEFT can be found in Ref. [32].
Some comments are in order regarding the experimental uncertainties in FF fits at low Q 2 (∼ 0.01 GeV 2 ). The dominant uncertainty in this case results from the normalization errors of the FF data. In the fits the normalization of different data sets is adjusted, using the given fit function or theoretical model, and requiring that G p E (0) = 1. Excellent fits can be achieved after this rescaling of the data sets, and one may be tempted to conclude that the radius could be extracted with great precision from them. However, such reasoning would be circular: The rescaling of the data sets depends on the fit function, and at low Q 2 the behavior of that function is governed by the radius, so that the procedure essentially forces the rescaled data to reproduce the assumed radius, resulting in a loss of sensitivity to the radius. This effect can be seen in the global fit uncertainties shown in Fig. 2. The dark shaded band shows the uncertainty of the global fit for a certain fixed value of the radius, i.e., constraining the first-order term in the fit function. This uncertainty vanishes rapidly in the limit Q 2 → 0, reflecting the trivial fact that the constrained fit with an assumed radius reproduces that value of the radius. The dashed area above the dark shaded band shows the variation in the global fit result under a change of the radius by ∆r = ±(0.01, 0.02, 0.03, 0.04) fm. The combined area represents the actual experimental uncertainty in the radius extraction in low-Q 2 fits. It is much larger than the fixed-radius uncertainty at low Q 2 and vanishes much more slowly in the Q 2 → 0 limit. Altogether, this points to a principal limitation of radius extraction from fits to low-Q 2 FF data. Our DIχEFT method overcomes this limitation by recruiting higher-Q 2 data for the radius extraction (up to Q 2 max ∼ 0.5 GeV 2 ), whose impact is only minimally affected by normalization errors.

C. Relation to empirical dispersive fits
The proton radius was extracted previously from dispersive FF fits in which the two-pion spectral functions were constructed by analytic continuation of empirical πN amplitudes [22,23]. In these approaches the twopion spectral functions are completely determined before they are placed in the dispersive integrals and used to evaluate FFs and radii. Our method is different in that the two-pion spectral functions are computed in DIχEFT and contain an unknown low-energy constant (related to the nucleon radii, cf. Appendix A), which can vary and adjust the strength of the spectral functions in the ρ meson peak and above [26]. This increases the flexibility of the FF description and enables a more robust radius extraction. We point out that the DIχEFT spectral functions at partial N2LO accuracy, evaluated with a realistic range of proton radii, agree very well with those of the Roy-Steiner analysis of Ref. [33], but differ significantly from those of Refs. [22,23] in the ρ meson mass region; see Ref. [26] for a detailed comparison. Even so, the empirical dispersive fits have consistently obtained proton radii ∼ 0.84 fm [22,23], in agreement with our result.

Appendix A: Form factor parametrization
In this appendix we summarize the DIχEFT calculation of the nucleon FFs and describe the parametrization used in the radius extraction in the present work. Further information about the method and other applications can be found in Refs. [24][25][26].
The nucleon electric FFs are separated into isovector and isoscalar components, , (A1) and represented as dispersive integrals over t ≡ −Q 2 , The integrands involve the imaginary parts of the FFs on the cuts at t > t thr > 0, which are known as the spectral functions. In the isovector component t thr = 4M 2 π , and the spectral function is organized as (A3) The first term accounts for the contribution of the twopion cut from t thr up to t max ∼ 1 GeV 2 and is calculated theoretically using the elastic unitarity relation in the ππ channel, the πN amplitudes computed in χEFT, and the empirical timelike pion FF; the explicit expressions are given in Refs. [25,26]. The couplings entering in the χEFT amplitudes at LO and NLO accuracy are determined by pion-nucleon scattering data. At partial N2LO accuracy the χEFT result involves one unknown low-energy constant, λ, which represents a free parameter; schematically The second term in Eq. (A3) accounts for the highmass states in the spectral function above t max and is parametrized by an effective pole, where the pole position M 2 1 = 2.1 GeV 2 is inferred from the e + e − annihilation data and the pole strength represents a free parameter; the justification for this approximation and its accuracy are discussed in Ref. [26]. The values of the parameters λ and a Since the value of the isovector charge is known, this leaves the isovector radius as the only free parameter of the isovector spectral function Eq. (A3). Note that the relation between the original parameters and the charge and radius implied by Eqs. (A6) and (A7) is linear, Expressing the parameters in terms of the radius, substituting the spectral function in Eq. (A2), and performing the dispersive integral, we obtain a parametrization of the spacelike isovector FF (t < 0) in the form The functions A V E (t) and B V E (t) represent, respectively, the dispersive integrals over the parts of the spectral function that are independent of (r V E ) 2 , and proportional to (r V E ) 2 . The functions have the same analytic structure as the full FF and embody the full complex t dependence in the spacelike region, as dictated by the dispersive representation. Note that the linear decomposition in Eq. (A9) results from the linear relation of (r V E ) 2 to the original theoretical parameters; it does not imply an expansion in t or other approximation to the complex t-dependence of the FF.
In the isoscalar component of Eq. (A2) t thr = 9M 2 π , and the spectral function is parametrized as the sum of two effective poles, The first pole is at the ω resonance mass and accounts for the 3π strength; the second pole is at the φ mass and effectively accounts for the φ resonance and other hadronic contributions at higher masses. The parameters a ω E and a φ E are fixed by the sum rules for the nucleon's isoscalar electric charge and radius, which imply a linear relation similar to Eq. (A8). Altogether we obtain a representation of the spacelike isoscalar FF analogous to the isovector case, Eq. (A9), Combining the isovector and isoscalar parametrizations, Eqs. (A9) and Eqs. (A13), we obtain a parametrization of the proton an neutron FFs in terms of the isovector and isoscalar radii, (A14) It may also be expressed directly in terms of the individual nucleon radii, The functions A p,n (t) describe the radius-independent part of the nucleon FFs (in the context of the dispersive parametrization defined above); the functions B E (t) andB E (t) describe, respectively, the parts proportional to the charge radii of the "same" and the "other" nucleon, under the constraints of isospin symmetry. Several properties of the functions follow immediately from their definitions: A p,n E (0) = Q p,n E , dA p,n E /dt(0) = 0, B E (0) = 0, 6 dB E /dt(0) = 1, where Q p,n E = {1, 0} is the nucleon electric charge. In particular, these conditions ensure that the nucleon FFs described by Eqs. (A15) and (A16) have the correct values at zero momentum transfer, G p,n E (0) = Q p,n E , and that the first derivatives of the nucleon FFs are controlled exclusively by the radii of the "same" nucleon, 6 dG p,n E /dt(0) = 6 dB E /dt(0) (r p,n E ) 2 = 6 (r p,n E ) 2 . (A21)  Figure 5 shows the contributions of the individual terms in the proton FF parametrization Eq. (A15) as functions of Q 2 = −t. One observes: (a) The radiusindependent term A p E starts with value 1 and derivative 0 at Q 2 = 0, cf. Eq. (A20); its value remains close to 1 for all Q 2 < 1 GeV 2 . (b) The proton-radius-dependent term (r p E ) 2 B E starts with value 0 at Q 2 = 0 and accounts for the derivative of G p E at Q 2 = 0; it causes most of the Q 2dependence of the FF (i.e., the deviation of G p E from 1) at Q 2 < 1 GeV 2 . (c) The neutron-radius-dependent term (r n E ) 2B E starts with value 0 and derivative 0 at Q 2 = 0; its contribution to G E remains < 0.01 for Q 2 < 1 GeV 2 . Altogether, the sum A p E + (r p E ) 2 B E practically accounts for the entire value of the FF and its Q 2 -dependence all at Q 2 < 1 GeV 2 ; the (r n E ) 2B E term represents a percent-level correction. This makes the parametrization Eq. (A15) particularly useful for extracting the proton radius from the FF data. The parametrization of the neutron FF Eq. (A16) follows the same pattern, with the role of the radii reversed.
The numerical evaluation of the FFs with the radiusdependent dispersive FF parametrizations, Eqs. (A15) and (A15), is extremely simple. The functions A p,n E (t), B E (t) andB E (t) can be pre-computed and tabulated as functions of Q 2 = −t. The FFs are then generated by multiplying these functions with the radii and combining the terms. This method was used to produce the plots of Fig. 1 and perform the radius extraction summarized in Fig. 3. The parametrization can be used also in other studies of low-Q 2 proton and neutron FFs. Tables of the functions are are available upon request.