Data-Driven Extraction of Hadron Radii

We describe the statistical Schlessinger Point Method. Grounded in analytic function theory, this method returns an objective assessment of the information contained in any data under consideration, which is independent of assumptions about underlying dynamics, and thus free from practitioner-induced bias. As a test to its robust nature and versatility, we apply it to a fully data-driven extraction of the proton, pion, and deuteron radii.


Introduction
Hadronic systems like baryons (proton, neutron, etc.), mesons (pion, kaon, etc.), and light nuclei (deuteron, etc.) are governed by Quantum Chromodynamics (QCD) the theory that describes strong interactions within the framework of the Standard Model (SM) of particle physics.And for understanding such systems, a single unifying question needs to be answered [1]: How does QCD dynamically generate the single mass-scale that characterizes the proton mass, its size (and, consequently, the characteristic confinement length), and, in fact, the natural scale for all nuclear physics?Notwithstanding the progress made by non-perturbative continuum methods in understanding the phenomenon of Emergent Hadron Mass (EHM) [2][3][4][5][6][7], its nature, the nature of confinement and the relation between them still remain fundamental open questions. 1 This is also due to the fact that while particle masses are known experimentally with extremely high precision (and many theoretical methods produce hadron spectra equivalent to the experimental ones), it is far more difficult to formulate and complete the reliable calculation of a radius within a given theory (and hadron radii are harder to measure).
The main tool to measure hadron radii is electron+hadron scattering, from which one determines the relevant form factor, G(Q 2 ), with high precision on a domain of squared-momentum transfer, Q 2 , that reaches deeply toward Q 2 = 0.Then, for systems with G(0) = 0, one obtains the radius via . ( There are two intertwined issues with this method.One is experimental: achieving the precision, reach to Q 2 ≈ 0, and density of Q 2 -coverage necessary for a reliable radius determination is very challenging.The other is theoretical: determining the radius through Eq. (1) requires knowledge (through modelling, fitting, etc.) of the form factor G one is measuring.Herein, we describe a new technique known as the statistical Schlessinger point method (SPM) [9][10][11], applicable in the same form to a wide variety of measurements [12][13][14][15][16][17].When applied to high-precision, dense electron-hadron scattering data stretching to very low momentum transfers, it returns a determination of the radius which is completely data-driven.The value obtained makes no reference to any functional form for analyzing the data, being instead based on a set of continued fraction interpolations capable of capturing both local and global features of the curve that the data are supposed to be describing.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, and for a robust estimation of the error by means of a statistical bootstrap procedure.

The (Statistical) Schlessinger Point Method
Let us denote with D N the set of all the N pairs values of a certain form factor G i associated to a given momentum transfer square Q 2 i : Within a subset D M ⊆ D N (typically with 4 < M N /2) one can construct the Schlessinger Point Method (SPM) continued fraction interpolator [9] where the M − 1 coefficients a i are recursively determined from the formulas and are such that The interpolator (3) can be cast in the rational form where P M and R M are polynomials whose degree is determined by the size M of the subset D M chosen: for M odd (respectively, even).While the SPM share the same rational expression (5) with a Padé approximant, the idea behind it is completely different.The latter is in fact defined as an expansion of a function near a specific point: its coefficients are thus constructed from values of higher derivatives at that point, so that the approximant's power series agrees with the one of the original function.In the SPM case, Eq. ( 5) is rather constructed using the original function values at different points. 2Thus, I M is able to capture the behavior of v over a range of values extending beyond the one of D M (and even D N ) without the need to compute any derivatives.
How well the interpolator (3) is capable of reproducing the entire dataset D N -and, correspondingly, the function G-mainly depends on the precision of the starting dataset D N .For example, when exact numerical data are considered as in (2), I M is practically indistinguishable from G independently from the number M of points in D M (provided that the latter set is large enough to capture the basic features of the curve G we wish to describe).
However, in the presence of uncertainties one rather has so that the G i are at most statistically distributed around the true curve G with variance i .A veracious reconstruction of the function G can be obtained also in this case if the SPM is combined with resampling to propagate uncertainties.To this end, we generate from (6) a replica set D r N by randomly drawing from each G i in D N a new one, G r i , from a normal distribution N with a mean value equal to G i and a standard deviation equal to its associated uncertainty i : Next, we randomly choose the subset D r M ⊆ D r N and proceed to construct the corresponding SPM interpolator I r M from Eq. (3).
Repeating these steps for a sufficiently large number of replicas n r , gives rise to a large population of interpolators that can be filtered according to suitable physics criteria.In the specific cases at hand we will only require that the interpolator I r M is smooth (continuous with all its derivatives, I r M ∈ C ∞ ) on the positive real axis3 R >0 .No further restriction is imposed, in particular no normalization constraint on the (electric) form factor at zero is required.
Each interpolating function so constructed defines an extrapolation to Q 2 = 0 from which the hadron radius can be determined through Eq. (1).For each replica, the resulting radius is taken as the mean of a suitable number n I of smooth interpolators, with an estimate of the error obtained by evaluating the variance of the replica's radii.In addition, the fact that M is not fixed leads to a second error source, indicated as σ δM and corresponding to the tandard deviation of the obtained r M j p across all j.Therefore, the value of the SPM extracted radius and the associated (statistical) error can be obtained as As in general the SPM results are independent from the number of input points M one has σ δM σ M j r .

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 [18].As explained in detail in [12], one seeks a function g which mimimizes the sum The first term in Eq. ( 9) 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 least-squares fit is obtained as λ → 0 (maximum penalty).
It can be then shown that: a) g is in fact the natural cubic spline interpolant for the values {y i } at {x i }; and, b) an estimate of the optimal value for λ can be determined by means of a (generalised) cross-validation procedure [19], see again [12] for details.

Final Procedure
Summarizing, the statistical SPM extraction of a hadron radius from a given dataset datasets, requires to perform the following steps: (i) Generate n r replicas for the given experimental central values and uncertainties; (ii) Smooth each replica with the associated optimal parameter λ r opt ; (iii) For each number of input points M j , determine the distribution of hadron radii r M j and its associated σ r considering the first n I interpolators which are smooth in the entire positive real axis R >0 ; (iv) Determine the overall σ δM by considering a new M j , repeating step (iii) for this new M-value, and evaluating the standard deviation of the distribution of r M j ; and (iv) Combine this information to obtain the final result for the proton radius and (statistical) error through Eq. ( 8).

Proton Case
The first example of the SPM extraction of a hadron radius is the determination of the proton size.
In 2010, the Particle Data Group (PDG) listed the value of the proton charge radius to be r p = 0.8768(79) [fm].This value was the average between measures obtained with two different methods: ep scattering; and the Lamb shift between different levels of hydrogen atoms.
In fact, with respect to the zero size case, the finite size of the nucleus induces corrections which are proportional to the probability of finding the electron inside the proton; and the latter are in turn proportional to the third power of the proton radius.However, the electron is so light with respect to the proton that it spends very few time inside it; thus, it is not the best of the probes for determining finite size corrections.A muon is much better: being 200 times heavier than the electron, it orbits much more closely to the proton; and this translates into a 10 million-fold increase in the sensitivity to the proton size.
The result obtained from muonic hydrogen (μH ) was found to be r p = 0.84184 (67) [fm] [21], that is 5σ away from the accepted value.This disagreement stimulated many new experiments, with conflicting results between the ep and μH measurements; and as shown in Fig. 1 (rows A-I), today we are dealing with the so-called "proton radius puzzle" [34].The only exception to this pattern has been the results from the most recent ep scattering experiment by the PRad collaboration [28,PRad], reporting the first ep scattering analysis to obtain a radius in agreement with that extracted from μH measurements: While implementing some significant improvements over earlier efforts, e.g., covering an extensive low-Q 2 domain reaching down to the deepest value yet achieved, viz.2.1 × 10 −4 ≤ Q 2 /GeV 2 ≤ 6 × 10 −2 , the PRad collaboration paid also special attention to the effect of data fitting function when extracting the charge radius 4 through Eq. (1).In particular, the optimal function form was selected by using a bootstrap procedure applied to pseudodata produced with fluctuations mimicking the Q 2 -binning and statistical uncertainty of their experimental setup.As no knowledge of actual PRad data was introduced, the PRad extraction is robust; but it also means that, nevertheless, the specific function form was chosen, namely [44] Rational(1, 1) Fig. 1 (color online) Proton electric radius, r p , extractions, various techniques: [30]; and [N] [12].The green-hashed rectangles display the anticipated precision of two forthcoming muon+proton scattering experiments: [J] [31,32]; and [K] [33].The two vertical grey bands indicate the separate error-weighted averages of the small-radius and large-radius extractions.From [7] Therefore a comparative analysis of the PRad and A1 data [20, A1], chosen as a benchmark owing to their reach in Q 2 (with 3.8 , high-precision, and sheer quantity (1422 cross sections measured recorded at beam energies of 0.18, 0.315, 0.45, 0.585, 0.72, 0.855 [GeV]), provides an ideal playground by which to illustrate the strengths of the SPM and see what it can reveal about the discrepancy between the ep and μH scattering results for r p .

Validation
For the SPM extraction of the proton radius from the PRad and A1 datasets, we have chosen [12] n r = 1 000, {M j = 5 + j | j = 1, . . ., n M ; n M = 12} and n I = 5 000.However, before we can trust the procedure when applied to real data, we need to prove that the SPM extraction method outlined above is robust, i.e. that it can reliably extract the proton radius in a variety of different situations. 5e can do this by showing that the SPM can reliably extract the proton radius from a variety of controlled cases; and, to this end, we use replica datasets 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 supplemented with fluctuations drawn according to a normal distribution.The chosen models are the ones discussed in [44].
For each of the 9 models described there, we proceed to generate pseudodata from the following experimental datasets: PRad at 1.1 and 2.2 GeV beam energy and their combinations; and the first 40 low-Q 2 points of the A1 dataset We then first start by showing that the distribution of the extracted SPM radii r p for a given M j is a Gaussian that is in fact nearly M j independent.This is shown in Fig. 2 for the particular case of the replicas generated from the dipole functional form (a.2) with r * p = 0.85 [fm] and the PRad kinematics at 1.1 [GeV] beam energy.A similar behavior is found in all other cases considered.For this generator and dataset one finds r p = 0.8501 [fm] and which shows the anticipated result that σ δM σ M j r .beam energy (upper-left corner), 2.2 [GeV] (lower-left corner) and combuined datasets (upper-right corner) [28] and the low-Q 2 A1 dataset (lower-right corner) [20].The names refer to the models used in [44], with: r p = 0.85 [fm] for the monopole [45], dipole and Gaussian forms; r p = 0.8277 [fm] for Kelly [46]; r p = 0.8681 [fm] for Arringthon I [47] and r p = 0.8965 [fm] for Arrington II [48]; r p = 0.8871 [fm] for Bernauer [49]; r p = 0.879 [fm] for Ye in [50]; and, finally, r p = 0.8765 [fm] for Alarcón [51].In all cases the SPM extrapolation is robust (|δr p | < σ r ) except for the PRad data at energy beam 2.2 [GeV] and combined kinematics, where the SPM for the Bernauer and Ye generators respectively is only marginally robust (|δr p | ∼ σ r ).
Adapted from [12] Fig. 4 (color online) Bias δr p , standard error σ r and RMSE of the SPM extrapolation of the proton radius for 10 3 replicas generated from the nine models described in the text.Notice the near-independence, within a given experimental dataset, of the RMSE from the generator used.Adapted from [12] In Fig. 3 we show the SPM bias, defined as δr p = r p − r * p and standard error σ r obtained when extracting the radius from all nine generators just described for the four different Q 2 -ranges and statistical errors: PRad at 1.1 [GeV] beam energy, at 2.2 [GeV] and their combination [28], as well as the low-Q 2 A1 data [20].In almost the totality of cases our procedure is robust, leading to a bias which is less than the standard error: |δr p | < σ r .However, for the PRad 2.2 [GeV] and combined kinematics, the SPM extrapolation for the Bernauer and Ye generators is only marginally robust (|δr p | ∼ σ r ).This is consistent with the findings in [44] where these two generators show the highest bias for all fitters independently from the Q 2 -binning (see in particular Fig. 11 therein).Notice also that the PRad dataset at 2.2 [GeV] leads to an error which is much smaller (roughly one third) than the one obtained from the 1.1 [GeV] data, and drives the error in the PRad combined binning (reducing it to roughly one half of the value obtained at 1.1 [GeV]).
Finally, following [44], Fig. 4 quantifies the goodness of the SPM proton radius extraction for all kinematics employed in terms of three indicators: the bias δr p ; the standard deviation σ r ; and the Root Mean Square Error (RMSE) As it can be seen, for any given experimental dataset, the RMSE values are approximatively independent from the generator used to obtain the replicas, which indicates that the SPM procedure for the extraction of the proton radius is robust.

Results
Let us begin then with a comparison between the value of r p quoted by the PRad Collaboration [I] and the one extracted through the SPM.
As noticed when validating the SPM extraction, the 2.2 [GeV] data lead to a lower value of r p and a smaller uncertainty (one-third the size) than the 1.1 [GeV] set; consequently the smaller error determined in the 2.2 GeV beam-energy dataset drives the overall uncertainty of the combined result Eq. ( 15) down to roughly 50% of that in Eq. (14a). 6Thus, within mutual uncertainties the SPM extraction reproduces the published PRad result.
Next, considering the first the low-Q 2 region of the A1 data, consisting of N = 40 data in the interval 3.8 × 10 −3 ≤ Q 2 /[GeV 2 ] ≤ 1.4 × 10 −2 , the SPM gives [12]: whereas when all the A1 data are included we get that is the SPM yields the same central value but a larger uncertainty.Evidently, extending the range up to Q 2 ∼ 1 GeV 2 introduces a sensitivity to the number of input points M j , as in this case we find σ δM ∼ σ The SPM results can be then combined, to produce the lowest entry in Fig. 1 -[N]: an ep value which is now compatible with that inferred from the Lamb shift in muonic hydrogenr p = 0.84136 (39) [fm] [21,23], the modern measurement of the 2S → 4P transition-frequency in regular hydrogen r p = 0.8335(95) [fm] [24], the Lamb shift in atomic hydrogenr p = 0.833 (10) [fm] [27], as well as matching the combination of the latest measurements of the 1S → 3S and 1S → 2S transition frequencies in atomic hydrogenr p = 0.8482 (38) [fm] [52], the muonic deuterium determination r p = 0.8356 (20) [fm] [53], and, last but not least, the latest PDG value r p = 0.8409 ± 0.0004 [fm] [54].This illustrates a possible solution to the proton radius puzzle relying on a sound extraction of r p from ep scattering experiments grounded on an analysis scheme that eliminates systematic error introduced by the use of specific, limiting choices for the functions employed to interpolate and extrapolate the available high-precision data.

Pion Case
Similarly to what happens in the proton case, the radius of the charged pion can be extracted once the charged pion elastic form factor, G π Q 2 , is measured.This has been done in pion+electron (πe) elastic scattering at low-Q 2 [55][56][57][58]; and, extractions of r π therefrom are broadly compatible even though within large uncertainties.Given the pion's fundamental role in binding nuclei and its direct links with EHM, this is therefore particularly unsatisfactory, 7 .
Typically, the quoted pion radius is [54]: a value obtained from the analysis of data collected in three distinct πe → πe experiments [56][57][58], supplemented by information from e + e − → π + π − reactions [60,61].This value has dropped 1.5 σ from that quoted two years earlier [62], r π = 0.672 (8) [fm]; thus, it is worth reconsidering available G π Q 2 data, using the SPM to obtain a fresh r π determination.However, reviewing available low-Q 2 πe → πe data [13], only the NA7 sets (Fig. 5) were judged worth considering for such a reanalysis, as they contain: N = 30 data on 0.015 ≤ Q 2 /GeV 2 ≤ 0.119 in the case of the NA7 1984 dataset [55]; and N = 46 data on 0.015 ≤ Q 2 /GeV 2 ≤ 0.25 in the case of the NA7 Fig. 5 (color online) NA7 pion form factor data: upper panel - [55]; lower panel - [56].Although collected thirty-five years ago, these remain the most dense and extensive sets of low-Q 2 data available.Concerning the 1986 set [56] lower panel, the SPM analysis in ref. [13] uses only the square-within-square data.The half-shaded data are too noisy to add information.From [13] 1986 dataset8 [56].The Q 2 -coverage, density, and precision of the data in Ref. [57] are instead insufficient to contribute anything beyond what is already available in the NA7 sets. 9

Validation
The models chosen for the validation of the SPM procedure in the pion case were a monopole, a dipole and a Gaussian: (20) from which, similarly to what we have done in the proton case, replica data sets are built from values of the pion elastic electromagnetic form factors (20) evaluated at the Q 2 points sampled in each NA7 experiment.Then, for a given value of M ∈ S M = {M j = 2 + 4 j | j = 1, 2, 3, 4}, all forms in Eq. ( 20), and each NA7 data set, give rise to a Gaussian distribution of SPM-extracted radii centred on the input value and practically independent of M, and with σ δM σ 3}, but in these instances one neither recovers the input radius nor obtains a Gaussian distribution.Such exceptional M-values were not encountered in our proton's analysis; so their appearance is probably related to the character of the NA7 data.Herein, therefore, we allow the validation procedure to impose a subset of M-values on our SPM analysis, viz.M ∈ S M .
Next, Fig. 6, displays the bias δr π as obtained using the SPM to extract r π from the three generators, the standard error σ r of Eq. ( 8), and the root mean square error, defined as in Eq. (13).As in all cases, |δr π | σ r and the RMSE values are approximatively independent of the model used to obtain the replicas, the SPM extraction is robust also in the case of the pion's datasets considered.Fig. 6 Bias, δr π , standard error, σ r , and RMSE, for SPM extrapolations of the pion radius from 10 3 replicas generated using the three models in Eq. (20).Notably, the RMSE is approximately independent of the model used.Left panels -data from Ref. [55]; and right panels -N = 35 data in Fig. 5, drawn from Ref. [56].From [13]

Results
Once applied to the NA7 datasets, the SPM returns the following values for the pion radius: (21a) which are mutually consistent.Notably, the SPM result is ≈ 1 σ below the error-weighted average of each of the original function-choice dependent determinations: r 84 π = 0.650(8) fm and r 86 π = 0.659(5) fm.The uncertainty-weighted average of the radii in Eq. ( 21) is then which implies that the SPM extraction supports the recent downward shift of the PDG average (see Fig. 7).
Clearly, extant πe elastic scattering data do not match analogous modern ep data in precision, density of coverage, or low-Q 2 reach: thus, at the moment, it is sensible to refrain from judging the true size of r π , as a veracious determination must await improved πe data 10 .

Deuteron Case
The last hadron radius we are going to consider is the one of the deuteron, the simplest compound nucleus in Nature (a neutron and proton, somehow bound by strong interactions).As was already the case for the proton and the pion, the deuteron size r D is not well known, and varies according to the method used: extracted from eD elastic scattering, the value usually quoted is r D = 2.111 (19) [fm] [29]; however, the value r D = 2.12562(78) 10 The situation for the kaon is even more dire.The commonly quoted value of r K = ± 0.031 [fm] [54], is obtained by averaging the results of separate monopole-squared fits to the distinct data sets in Refs.[63,64].For these datasets, the SPM interpolants relax to a linear least-squares fit to the data which returns r K = 0.53 [fm]; being consistent with the value obtained with a monopole fit, this outcome further highlights the inadequacy of available data so far as a precise determination of the kaon radius is concerned.Fig. 7 (color online) Pion charge radius: purple circle (PDG, 2020) [54]; blue square (PDG 2018) [62]; cyan diamond [58]; green up-triangle [56]; light-green down-triangle [57].The vertical grey band marks the uncertainty in Eq. ( 19).SPM results [13] -gold stars: SPM A-84 from [55]; SPM A-86 from [56]; and orange star -SPM combined, the average in Eq. (22).From [13] [fm] is reported from measurements of 2S transitions in μD atoms [53].This is 25-times more precise than the eD scattering value and, referred to its own error, 19σ larger.An analogous measurement using standard atomic deuterium yields r D = 2.1415 (45) [fm] [65], which is 3.5σ above the μD value.
To clarify this "deuteron radius puzzle", a new experiment has been proposed at Jefferson Lab [66,DRad], which is based upon the PRad experimental configuration, and aims for an extraction of r D from eD elastic scattering with precision better than 0.25%.

Results from the SPM Procedure
To support the DRad effort and following the PRad approach, in [67] a study was made to identifying fitting functions that can robustly extract r D from the projected DRad measurements.The models chosen are of two different kind: (i) four parametrizations of the available experimental data; and (ii) nine theoretically based models ranging from the elementary to the extremely elaborated. 11Again, after describing in certain detail the chosen models, we will show that the statistical SPM can provide a powerful alternative to [67].
(a) Abbott-I model [68].In this parametrization, the deuteron charge form factor is represented by where G 0 C = 1 is a normalizing factor fixed by the deuteron charge, Q 0 C = 0.830747 [GeV] and the coefficients a i C read (in [GeV −2i ]) [68]: The radius associated to this model is r D = 2.09345 [fm].
(b) Abbott-II model [68,69].In this second parametrization the deuteron charge form factor is given by where and, finally, The 24 parameters and The corresponding deuteron charge radius is in this case r D = 2.08823 [fm].(c) Parker model [70].This parametrization is essentially is a refit of the first two parametrizations with imposed constraints to prevent singularities in the functional forms of the form factors.One has where G 0 C = 1 and Q 0 C = 0.830747 [GeV] as before, whereas are the a i parameters are given by (in a 2 a 3 a 4 a 5 0.695719 13.9994 0.695719 0.695502 0.695719 The associated deuteron charge radius is in this case r D = 2.06380 [fm] (d) Sum of Gaussian (SOG) model [68,71,72].The fourth parametrization utilizes the SOG method, by which where a dipole form and, finally, a Gaussian form For all these models we set r D = 2.1 [fm].(g) Quadratic and cubic models.These models read and represents respectively a second and third order polynomial fit to the low-Q 2 theoretical points calculated using the covariant spectator theory model of np scattering presented in [74].In the following we will work with the predicted DRad kinematic range described in [67] with beam energy at 1.1 and 2.2 [GeV], using bin-by-bin statistical uncertainties from 0.02 to 0.07% and systematic uncertainties from 0.06 to 0.16%.Let us start by observing that considering the different planned DRad electron beam energy settings, E γ , and defining accordingly S 1.1 M and all thirteen fitting models, the distribution of SPM-extracted radii is a Gaussian centred on the radius input value r * D and whose characteristics are practically independent of M, with σ δM < σ It is now worth comparing the SPM results and Figs. 8, 9 with Ref. [67,5,9] and the associated discussion.The primary conclusion of the latter is that the methods used by the PRad Collaboration to identify "best fitter" functions for use in extracting the proton radius are unsuitable for the deuteron case.So, a new method, tuned specifically to eD scattering, needed to be developed.It yielded two new best fitter Fig. 8 (color online) Bias δr D and associated standard error σ r for the SPM extrapolation of the deuteron radius from 10 3 replicas generated from the thirteen models discussed in the text, for the different Q 2 ranges corresponding to the projected DRad kinematics and bin-by-bin uncertainties.For the 1.1 GeV kinematics (left panel) the SPM provides a robust extrapolation for all models (|δr D | < σ r ), except for the Gaussian and the quadratic ones, with the IA model being marginally robust; similarly for the combined 1.1 and 2.2 GeV kinematics (right panel) the extrapolation is robust for all models but the dipole, Gaussian, quadratic and cubic.The 2.2 GeV kinematics (center) leads instead to large biases and errors (notice the change of scale in this panel); while the SPM leads to a robust extrapolation for models based on the experimental data, analysis of this kinematics alone should be avoided.Adapted from [16] functions modRational (1,1) neither of which coincides with Eq. ( 11), the form used in connection with the proton radius extraction.This contrasts markedly with the SPM results presented.As shown above, for practical purposes, all procedures used in connection with the reliable extraction of the proton and pion radii are equally effective for a SPM extraction of the deuteron radius: one obtains a robust extraction of the deuteron radius in 85% (69%) of cases using DRad 1.1 GeV beam (1.1 + 2.2 combined) kinematics. 12Finally, the SPM delivers an average variance σr between 5.5 × 10 −3 and 7.5 × 10 −3 for the 1.1 [GeV] and combined 1.1 and 2.2 [GeV] kinematics respectively, which is to be compared with the value 3.5 × 10 −3 obtained using the eD-tuned best fitters in Ref. [67,Fig. 6].

Conclusion
The SPM provides a method for delivering an estimate, with known uncertainty, to the curve underlying a set of data; and, by constructing continued fractions interpolants to those data, it allows the extrapolation of quantities outside the domain of existing measurements.Its rigorous mathematical foundation in analytic Fig. 9 (color online) Bias δr π , standard error σ r and RMSE of the SPM extrapolation of the pion radius for 10 3 replicas generated from the three models described in the text.Notice the limited variation of the RMSE with the generator chosen for the 1.1 GeV and combined 1.1 and 2.2 GeV kinematics.Adapted from [16] function theory eliminates the need for assumptions about theories which may underly the data; and, making no reference to those theories, it removes practitioner-induced bias.It thus provides an objective expression of the information contained in the data under consideration, serving as a test for both the experiment as well as potentially relevant theory.
The range of existing SPM applications goes well beyond the extraction of hadron radii illustrated here, including: the extrapolation of theory predictions for hadron form factors to large values of Q 2 , e.g., Refs.[75][76][77]; continuation from Euclidean to Minkowski space in order to obtain light-front specific quantities from Poincaré-covariant wave functions, e.g., Refs.[78,79]; overcoming the challenge of large mass imbalances between valence degrees-of-freedom in delivering predictions for the semileptonic decays of heavy+light mesons, e.g., Ref. [80]; the analysis of light-nuclei deep inelastic scattering data [81] and its interpretation in terms of parton distribution functions [14]; the identification of resonance poles in particle decays [17]; and the identification of crossing odd states (the odderon) in the t-channel of high-energy elastic hadron+hadron scattering [82].
Given a set of precise data, densely covering an appropriate domain, the SPM will return an unbiased estimate of a desired observable with an uncertainty that quantitatively reflects the quality of the data and the length of the extrapolation.

(see Fig. 1 -
Row I, and notice in particular the tension between the PRad result and the values obtained in analyses of all other ep scattering data, Rows [A, C, G]).

Fig. 2 (Fig. 3 (
Fig.2(color online) Probability Distribution Function of the estimated SPM proton radius r p extracted from 10 3 replicas generated from a dipole functional form (with r p = 0.85 [fm]) applied to the PRad 1.1 data kinematics.The reconstructed PDF is a Gaussian that is nearly independent from the number of input points M j .From[12] to the A1 Collaboration original estimate r A1−coll.p = 0.879 ± 0.005 stat ± 0.006 syst[fm], the SPM extraction reconciles the PRad and A1 values and, perhaps more importantly, shows both to be consistent with the μH experiments.
. 8 provides evidence for the robustness of the SPM r D extraction in terms of the bias δr D = r D −r * D and standard deviation σ r .More specifically we have that: • With beam-energy 1.1 [GeV] (left panel), the SPM extraction is robust, |δr d | < σ r , for all data generators except the Gaussian and quadratic models (with the IA model being marginally robust, |δr d | ∼ σ r ); • For the combined 1.1 and 2.2 [GeV] data (right panel), the SPM radius extraction is robust for all data generators except the Gaussian, quadratic and cubic models; • Finally, with 2.2 [GeV] beam energy kinematics (central panel), one finds large biases and errors for many generators (IA, IAMEC, dipole, Gaussian, quadratic and cubic); and although the SPM leads to a robust extraction for generators based on experimental data, these results suggest that analysis of the 2.2 [GeV] kinematic setup alone should be avoided.Finally, Fig. 9 quantifies the goodness of the SPM proton radius extraction for all kinematics employed in terms of the bias δr D , the standard deviation σ r ; and the RMSE = (δr D ) 2 + σ 2 r .As expected, for the 1.1 [GeV] (left panel) and combined 1.1 and 2.2 [GeV] (right panel) kinematics, the RMSE values are practically independent of the generator used to produce the data replicas.As a result, the SPM returns a reliable expression of the deuteron radius actually encoded in the data.Considering the middle panel, i.e., 2.2 [GeV] beam kinematics alone, the elementary generators are again seen to present a challenge.
γ 2 i are not all independent being constrained by a set of 12 relations, see Ref. [67, Eq. (A5)]; they are given by [73]tional theoretical models[73].The first four theory based models we consider are the following: IA (relativistic impulse approximation); IAMEC (relativistic IA plus meson exchange current), RSC (relativistic IA with Reid Soft Core), and RSCMEC (relativistic IA plus meson exchange current with Reid Soft Core).The corresponding form factors are obtained as spline interpolators of the data appearing in Table1which are the results of the theoretical calculations discussed in[73].The associated deuteron radii are: IA, r D = 2.15798 [fm]; IAMEC, r D = 2.15727 [fm]; RSC, r D = 2.0997 [fm]; RSCMEC, r D = 2.09912 [fm].(f) Simple functional forms.In this category we have considered naive models that imitate possible extrapolations at Q 2 → 0. Namely, a monopole form

Table 1
[73]eron charge form factor G D C obtained from the theoretical calculations discussed in[73]for four different models: IA (relativistic impulse approximation); IAMEC (relativistic IA plus meson exchange current), RSC (relativistic IA with Reid Soft Core), and RSCMEC (relativistic IA plus meson exchange current with Reid Soft Core)