Fresh extraction of the proton charge radius from electron scattering

We present a novel method for extracting the proton radius from elastic electron-proton ($ep$) scattering data. The approach is based on interpolation via continued fractions augmented by statistical sampling and avoids any assumptions on the form of function used for the representation of data and subsequent extrapolation onto $Q^2\simeq 0$. Applying the method to extant modern $e p$ data sets, we find that all results are mutually consistent and, combining them, arrive at $r_p=0.847(8)\,$fm. This result compares favourably with values obtained from contemporary measurements of the Lamb shift in muonic hydrogen, transitions in electronic hydrogen, and muonic deuterium spectroscopy.

1. Introduction -The proton is Nature's most fundamental bound state. Composed of three valence constituents, two u-quarks and one d-quark, it seems to be absolutely stable: in the ∼ 14-billion years since the Big Bang, proton decay has not been observed. The proton's extraordinarily long lifetime is basic to the existence of all known matter. Yet, the forces responsible for this remarkable feature are not understood.
Proton structure is supposed to be described by quantum chromodynamics (QCD), the Standard Model quantum field theory intended to explain the character and interactions of the proton (and all related objects) in terms of gluons (gauge fields) and quarks (matter fields) [1]. Today, the proton's mass, m p , can be calculated with good accuracy using modern theoretical tools [2][3][4]; but that is not the case for its radius, r p .
The proton's radius is of particular importance because it relates to the question of confinement, viz. the empirical fact that no isolated gluon or quark has ever been detected. The value of r p characterises the size of the domain within which the current-quarks in QCD's Lagrangian may rigorously be considered to represent the relevant degrees of freedom. (A clearer notion of confinement may appear in a proof that quantum SU c (3) gauge field theory is mathematically well-defined, i.e. a solution to the Yang-Mills "Millennium Problem" [5].) Moreover, it is not just strong interactions which feel the size of r p . An accurate value of the proton's charge radius is also crucial to a precise determination of quantities in atomic physics, such as the Rydberg constant and Lamb shift.
Naturally, m p and r p are correlated. A solution to the Standard Model will deliver values for both. Hence, precise measurements are necessary to set rigorous benchmarks for theory. The problem is that whilst the relative error on m p is ∼ 10 −10 , measurements of r p now disagree amongst themselves by as much as eight standard * binosi@ectstar.eu † cdroberts@nju.edu.cn deviations, 8σ, as illustrated in Fig. 1 -upper panel. This conflict, which emerged following extraction of the proton radius from measurements of the Lamb shift in muonic hydrogen (µH) [6], has come to be known as the "proton radius puzzle" [7,8].
Many solutions of this puzzle have been offered, e.g. some unknown QCD-related corrections may have been omitted in the muonic hydrogen analysis and their inclusion might restore agreement with the electron-based experiments that give a larger value; the discrepancy could signal some new interaction(s) or particle(s) outside the Standard Model, which lead to a violation of universality between electron (e) and muon (µ) electromagnetic interactions; or some systematic error(s) has (have) hitherto been neglected in the analysis of electron scattering. Empirically, novel experiments have been proposed in order to test various possibilities, including µp elastic scattering (MUSE) [21] and ep scattering at very low momentum transfer (PRad) [22]. PRad recently released its result [14]: Significantly, this is the first published analysis of an ep scattering experiment to obtain a result in agreement with the radius extracted from µH measurements. In performing and analysing the ep scattering experiment, the PRad collaboration implemented a number of improvements over previous efforts, which included: reaching the lowest yet achieved momentum-transfersquared, Q 2 = 2.1 × 10 −4 GeV 2 ; and covering an extensive domain of low Q 2 : 2.1×10 −4 ≤ Q 2 /GeV 2 ≤ 6×10 −2 . Moreover, since the charge radius is obtained as where G p E (Q 2 ) is the proton's elastic electromagnetic form factor, PRad paid careful attention to the impact of the choice of fitting form on the extracted charge radius, an issue highlighted previously [23][24][25][26][27][28][29][30][31]. Notably, their functional form was predetermined through a bootstrap procedure applied to pseudodata generated with fluctuations mimicking the Q 2 -binning and statistical uncertainty of the experimental setup, i.e. without knowledge of the actual PRad data [32]. While this procedure renders the PRad extraction robust, it also means that, ultimately, a specific functional form was chosen [32].
We reanalyse the PRad data [14] and also data from the A1 Collaboration [13] using a statistical Schlessinger Point Method (SPM) [18,19]. Following Ref. [20], the SPM has been used widely and effectively to solve numerous problems in hadron physics, especially those which demand model-independent interpolation and extrapolation, e.g. Refs. [33][34][35][36][37]. In this approach, no functional form is assumed. Instead, one arrives at a set of continued fraction interpolations capable of capturing both local and global features of the curve that the data are supposed to be measuring. This latter aspect is crucial because it ensures that the validity of the constructed curves extends outside the data range limits, ultimately allowing for the evaluation of the curves' first derivative at the origin. A robust estimation of the error is also obtained by means of a statistical bootstrap procedure [38]. 2. Theory for interpolation and extrapolation of smooth functions -The foundation for our fresh analysis of G p E (Q 2 ) data, obtained from ep scattering and available on Q 2 min ≤ Q 2 ≤ Q 2 max , is the SPM. In general, given N pairs, D = {(x i , y i = f (x i ))}, being the values of some smooth function, f (x), at a given set of discrete points, a basic SPM application constructs a continued-fraction interpolation: in which the coefficients {a i |i = 1, . . . , N − 1} are constructed recursively and ensure C N ( The SPM is related to the Padé approximant; and the procedure accurately reconstructs any analytic function within a radius of convergence fixed by that one of the function's branch points which lies closest to the domain of real-axis points containing the data sample.
For example, suppose one considers a monopole form factor represented by N > 0 points, each one lying on the curve; then using any one of those points, the SPM will exactly reproduce the function.
In the physical cases of interest herein, one deals with data that are distributed statistically around a curve for which the SPM must deliver an accurate reconstruction. Given that all sets considered are large, N is big enough to enable the introduction of a powerful statistical aspect to the SPM. Namely, one randomly selects M < N points from the set D, typically with 4 < M N/2 [33,35]. In theory, one can then obtain C(N, M ) different interpolating functions; in practice, this number is reduced by introducing physical constraints on their behavior. The minimal N we consider is N = 33, i.e. the PRad data set at a beam energy of 1.1 GeV; thus, choosing M ∈ [6,17] gives O(10 6 − 10 9 ) possible interpolators, out of which we select the first 5 × 10 3 corresponding to smooth monotonic functions on the entire Q 2 domain. No further restriction is imposed; specifically, no unity constraint on G p E (Q 2 = 0) is required. Each interpolating function defines an extrapolation to Q 2 = 0, from which r p can be extracted using Eq. (2). For a given value of M , the value of the radius is then obtained as the average of all results obtained from the 5 000 curves.
To estimate the error associated with the SPM determined proton radius, one needs to account for the experimental errors in each of the data sets. This can be achieved by using a statistical bootstrap procedure. To wit, we generate 1 000 replicas for each set by replacing each datum by a new point, randomly distributed around a mean defined by the datum itself with variance equal to its associated error. The probability distribution function, N (µ, σ), characterising the r p values extracted via the procedure described above on each replica is, to a very good approximation, normal, with average r M p , standard deviation σ = σ M r , skewness β M 3 ≈ 0 and kurtosis β M 4 ≈ 3 (see Supplemental Material, Fig. S1). In addition, the fact that M is not fixed leads to a second error source σ δM , which can be estimated by changing M → M , repeating the aforementioned procedure for this new M -value, and evaluating the standard deviation of the distribution of r M p for different M values.
Consequently, for each kinematics, the SPM result is Herein, we compute results for each one of the values {M j = 5 + j | j = 1, . . . , n M ; n M = 12}, so that for any given data set we have 60-million values of r p , each calculated from an independent interpolation; and, typically, we find σ δM σ M j r for all js in the range specified above. (See Supplemental Material, Eq. II.1).

3.
Smoothing with roughness penalty -Before implementing the statistical SPM, however, one issue must be addressed. Namely, as highlighted above, sound experimental data are statistically scattered around that curve which truly represents the observable. They do not lie on the curve; hence, empirical data should not be directly interpolated.
A solution to this problem is smoothing with a roughness penalty, an approach we have implemented following the ALGOL procedure detailed in [39] and which we now sketch. One begins by assuming the data are good, viz. they are a true measurement of an underlying smooth function. The next step is to identify the correct basis functions for the smoothing operation. These are provided by cubic splines, defined as follows. Consider a sequence of increasing numbers x 1 < x 2 < · · · < x in some interval I = [a, b], a < b; and call these numbers knots. A function g defined on I is a cubic spline with-respect-to (wrt) the knots {x i } if the following two conditions are satisfied. (i ) g is a cubic polynomial on each of the m + 1 subintervals: g is a C 2 function, viz. continuous with two continuous derivatives. All cubic splines (wrt knots {x i }) form a vector space of functions with + 4 degrees-of-freedom. A set of basis functions for this space is: is called "natural" if its second and third derivatives vanish at the interval's endpoints. Now consider the Sobolev space S[I ] of C 2 functions on I . In the roughness-penalty approach to smoothing, one seeks that function g ∈ S[I ] which minimises The first term in Eq. (5) quantifies the data-fidelity of g. The second term introduces the roughness penalty via the smoothing parameter λ ∈ [0, 1]: the original data are recovered for λ → 1 (no penalty) and a linear leastsquares fit is obtained as λ → 0 (maximum penalty).

S[I]
is an infinite-dimensional space; but it can be shown that when > 1, the minimiser lies in the finite dimensional space of natural cubic splines with knots located at the data points, {x i }. In fact, the following theorem holds. Let g be any smooth function on I for which g(z i ) = y i , i = 1, . . . , ; and suppose that s is the natural cubic spline interpolant for the values with equality if and only if g ≡ s.
At this point the smoothing parameter λ is somewhat arbitrary, with a typical value near 1/(1 + h 3 /6) where h is the average spacing of the data sites: h ∼ 5 × 10 −4 , 1.5 · 10 −3 , and 1 × 10 −3 for the PRad data at beam energy 1.1 GeV, 2.2 GeV and their combination [14]; and h ∼ 1.7 × 10 −3 for the Mainz data [13]. On the other hand, an estimate of the optimal value for λ can be determined by means of a (generalised) cross-validation procedure [40], which we now explain. Pretend that observation "k" is lost, so that only the remaining − 1 points are available for constructing a smoothing spline with respect to λ. Denote the solution of this reduced problem byš k ; by definition,š k minimises The quality ofš k as a predictor for a new observation can be judged by the difference y k −š k (x k ); and this leads to the cross-validation procedure, i.e. λ opt is the value of λ for which the following function is minimised: On average, we find λ opt = 1 − , with ∼ 2 × 10 −6 , 6 × 10 −6 , and 4 × 10 −6 for PRad 1.1 GeV beam energy, 2.2 GeV, and combined values, respectively, and ∼ 1.5× 10 −6 for the Mainz data. 4. Final procedure, validation and results -As explained above, the SPM extraction of the proton radius from a set of ep scattering data requires the following steps: (i ) generate 1 000 replicas for the given experimental central values and uncertainties; (ii ) smooth each replica with the associated optimal parameter λ opt ; (iii ) for each number of input points M j ∈ [6,17], determine the distribution of proton radii r M j p , its associated σ M j r , and the overall σ δM ; and (iv ) combine this information to obtain the final result for the proton radius and (statistical) error through Eq. (4).
One might wonder if the proposed SPM extraction method is robust, i.e. whether or not it can reliably extract the proton radius in a diverse array of cases. We checked this by using a wide variety of models that have been employed to fit the world's ep scattering data [41][42][43][44][45][46][47] to generate a proton electromagnetic form factor G E p with a known value for the radius. From these, we generated replicas with the Q 2 -binning and uncertainties of the PRad [14] and A1 [13] data sets. In all cases, regardless of the generator employed, we found that the SPM returns the radius value used to generate the pseudodata; and, furthermore, the result is practically independent of the number of initial input points M j . (See Supplemental Material, Sec. II.2).
The first G E p data from which we extracted the proton radius are those from the PRad experiment [14], which reported data using 1.1 GeV (N = 33) and 2.2 GeV (N = 38) electron beams. Analysed separately, the SPM gives: r 1.1 p /fm = 0.842 ± 0.008 stat for the 1.1 GeV data; and r 2.2 p /fm = 0.824 ± 0.003 stat for the 2.2 GeV data. Treated alone, the PRad data at 2.2 GeV leads to a lower value of the proton radius and a smaller error (one-third the size) than are obtained from the 1.1 GeV data; moreover, it drives the error in the PRad combined binning, reducing it to roughly one-half the value obtained using the 1.1 GeV data alone. These observations accord with those made by the PRad Collaboration [14]: see, in particular, Fig. S16 in the associated Supplemental Material.
Our final result, obtained from a combined analysis of the PRad data, is which is displayed in Fig. 1  Therefore, we also applied our method to the A1 data, first restricting the analysis to the low-Q 2 region, consisting of N = 40 data in the interval 3.8 × 10 −3 ≤ Q 2 /GeV 2 ≤ 1.4 × 10 −2 . In this case we obtained ( Fig. 1 -lower panel): Eliminating the restriction to the low-Q 2 region yields the same central value but a larger error: r A1 p /fm = 0.857 ± 0.021 stat . In this case, σ δM ∼ σ M j r ; so, extending the range of squared momentum transfer up to Q 2 ∼ 1 GeV 2 limits the ability of the SPM to provide an M -independent result.
The original A1 Collaboration estimate is [13]: r A1−coll. p /fm = 0.879 ± 0.005 stat ± 0.006 syst . Thus, whilst the SPM reanalysis of the A1 data, Eq. (10), has a larger statistical uncertainty, it yields a value that agrees with both the PRad estimate and the µH experiments. 5. Conclusions -We calculated the proton charge radius, r p , by analysing high-precision ep scattering data obtained in modern experiments [13,14] using a statistical sampling approach based on the Schlessinger Point Method for the interpolation and extrapolation of smooth functions. An important feature of this scheme is that no specific functional form is assumed for the interpolator, i.e. it produces a form-unbiased interpolation as the basis for a well-constrained extrapolation. All considered ep scattering data sets yielded consistent results [Eqs. (9) and (10)]; and combining them we find: (11) which is indicated by the gold band in the lower panel of Fig. 1.
Consequently, according to this analysis, there is no discrepancy between the proton radius obtained from ep scattering and that determined from the Lamb shift in muonic hydrogen -r p = 0.84136(39) fm [6,10], the modern measurement of the 2S→4P transition-frequency in regular hydrogen -r p = 0.8335(95) fm [11], the Lamb shift in atomic hydrogen -r p = 0.833(10) fm [15], the combination of the latest measurements of the 1S→3S and 1S→2S transition frequencies in atomic hydrogenr p = 0.8482 (38) fm [16], or even the muonic deuterium determination r p = 0.8356 (20) fm [17]. Furthermore, our analysis suggests that the explanation for the mismatch which spawned the "proton radius puzzle" lies in an underestimation of the systematic error introduced by the use of specific, limiting choices for the functions employed to interpolate and extrapolate ep scattering data. SM I. Validation models and procedure. -The validity of the SPM procedure for extracting the proton radius can be checked against replica data sets built from values of the proton elastic electromagnetic form factor, G p E , evaluated at the experimentally available Q 2 and generated using specific models in which the proton radius r p is known a priori. To mimic the variability of real data, these values are then redistributed by introducing fluctuations drawn according to a normal distribution. The models chosen are those discussed in Ref. [32], which range from (a) standard functions to (b) parametrisations of experimental data and, finally, (c) a theoretical calculation. a. Standard functions. In this category, we have considered [41]: a monopole form corresponding to a Yukawa charge distribution of the proton; a dipole form Bin counting FIG. S1. Probability distribution function associated with the SPM extraction of rp from 1 000 replicas generated using a dipole model, Eq. (a.2), with rp = 0.85 fm, from PRad data at 1.1 GeV beam energy kinematics. Plainly, the reconstructed probability distribution is a Gaussian whose characteristics are practically insensitive to the number of input points Mj.
corresponding to an exponential charge distribution; and, finally, a Gaussian form For all these models, r p = 0.85 fm.
b. Parametrisations of experimental data. The first parametrisation considered is that in Ref. [42]: with m p = 0.938272 GeV being the proton mass. The coefficients a i , b i , and the expected value of r p are: 24  Two different fits are employed in Refs. [43,44]. In the first [43], the proton's form factor is expressed via the following series: The second [44] implements r p = 0.8965 fm and is characterised by a fifth-order continuous fraction: Ref. [45] employed a tenth-order polynomial expansion of G p E : with proton radius r p = 0.8871 fm and the following coefficients (in GeV −i ) The final parametrisation is that in Ref. [46], which is a fit to the world's data, yielding r p = 0.879 fm, expressing the proton's electric form factor as with t cut = 4m 2 π , the two-pion particle production threshold (m π = 0.13957 GeV is the pion mass), and t 0 = 0.7 GeV. The coefficients are: (b.10) c. Theoretical calculation. Ref. [29] exploited a form factor analysis method that uses a combination of effective field theory and dispersion relations. A parametrisation of the results, yielding the CODATA value of the proton radius [9]: 0.8765 fm, is given by Observation 1. For a given M j , the distribution of SPM extracted radii is a Gaussian whose characteristics are practically independent of M j . This is illustrated in Fig. S1 for the case of replicas generated using the dipole functional form, Eq. (a.2), with r * p = 0.85 fm and the PRad kinematics at 1.1 GeV beam energy. Qualitatively identical behaviour is found in all other cases considered. For this generator and data set, the SPM yields r p = 0.8501 fm and r . Observation 2. Defining the bias as δr p = r p − r * p , then the SPM extraction of the proton radius is robust: in almost all cases, |δr p | < σ r , where σ r is the standard error in Eq. (4) -main text. This is demonstrated by Fig. S2, which displays δr p as obtained using the SPM to extract r p from all nine generators described   [14]; and the low-Q 2 Mainz data set (lower-right) [13]. In almost all cases, the SPM extrapolations are robust (|δrp| < σr). For the PRad data at 2.2 GeV and combined kinematics, the SPM results for the Eq. (b.6) and (b.8) generators are marginally robust (|δrp| ∼ σr).
in Sec. SM I over four different Q 2 ranges and statistical errors: PRad at 1.1 GeV beam energy, 2.2 GeV beam, 1.1, 2.2 GeV combined [14]; and the low-Q 2 Mainz data [13]. In all but 3 of the 36 cases, |δr p | < σ r . The three exceptions are PRad 2.2 GeV kinematics with the generators in Eqs. (b.6), (b.8), and the PRad-combined with the Eq (b.8) generator. This outcome is consistent with the analysis in Ref. [32], which reported that the Eq. (b.6), (b.8) generators consistently showed the highest bias of all fitters, independent of the Q 2 -binning (see Fig. 11 therein). In this connection, it is also notable that, as already remarked, the PRad 2.2 GeV error is much smaller than that associated with the 1.1 GeV data, driving down the error in the combined set to roughly one-half of that found with the 1.1 GeV data alone. Bias, δrp, standard error, σr, and RMSE, Eq. (II.2), for SPM extrapolations of the proton radius from 10 3 replicas generated using the nine models described in Sec. SM I. Notably, within a given experimental data set, the RMSE is approximately independent of the generator used.
(RMSE) RMSE = (δr p ) 2 + σ 2 r , (II.2) then for any given experimental data set, as shown in Fig. S3, the SPM analysis produces RMSE values that are approximatively independent of the generator used to obtain the replicas. This means that the SPM procedure satisfies a standard "goodness of fit" criterion [32], in consequence of which our SPM extractions of the proton radius can objectively be judged as sound.