The origin of the proton radius puzzle

The dimension of the proton, the basic building block of matter, is still object of controversy. The most precise electron-proton scattering data at low transferred momenta are re-analyzed and the extraction of the proton radius is discussed. A recent experiment from the JLAB-CLAS collaboration gives a small value for the radius (The symbol REα\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_E^\alpha $$\end{document} stands for the root-mean-square charge radius of the proton ⟨rE2⟩\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{\langle r_E^2\rangle }$$\end{document}, obtained by the experimental or theoretical Collaboration α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document}.) RECLAS=(0.831±0.007stat±0.012syst)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_E^\mathrm{CLAS}= (0.831\pm 0.007_\mathrm{stat}\pm 0.012_\mathrm{syst})$$\end{document} fm (Xiong et al. in Nature 575:147, 2019), in contrast with previous electron scattering experiments, in particular with the MAINZ experiment (Bernauer et al. (A1 Collaboration), Phys. Rev. C 90:015206, 2014) that concluded REMAINZ=(0.879±0.005stat±0.004syst±0.002model±0.004group)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_E^\mathrm{MAINZ}= (0.879\pm 0.005_\mathrm{stat}\pm 0.004_\mathrm{syst}\pm 0.002_\mathrm{model}\pm 0.004_\mathrm{group})$$\end{document} fm. The experimental results are re-analyzed in terms of different fits of the cross section and of its discrete derivative with analyticity constraints. The uncertainty on the derivative is two orders of magnitude larger than the error on the measured observable, i.e., the cross section. The systematic error associated with the radius is evaluated taking into account the uncertainties from different sources, as the extrapolation to the static point, the choice of the class of fitting functions, and the range of the data sample.


Introduction
The proton, as other fundamental particles, is defined by the mass, the charge, and quantum numbers as spin, isospin, and parity. Many questions are still open about its fundamental properties, as the origin of its mass and spin. The proton, being a composite particle, has a finite spatial dimension, quantified by a radius, more precisely by the root-meansquare charge radius R E ≡ r 2 E . a e-mail: simone.pacetti@unipg.it (corresponding author) b e-mail: egle.tomasi@cea.fr Until 2010 there was a consensus about its value from several experiments that were officially compiled by the wellrecognized Committee on Data for Science and Technology (CODATA). The experiments were of two kinds: the Lamb shift measurements with hydrogen spectroscopy and electron-proton elastic scattering. Most of the recent literature on the subject, as regularly reported in this journal, arose since Lamb shift measurements could be performed on muonic atoms. The electric field felt by the electron or the muon inside the atom is modified by the proton size, and consequently the level splitting. The proton radius appears in a correction term to the measured Lamb shift. Since the muon mass is about 200 times larger than the electron one, the muonic hydrogen atom is more compact, enhancing the sensitivity to this term. The consequence is a smaller error on the proton radius, compared to standard hydrogen spectroscopy. What came as a surprise was that the extracted value of the radius R H μ E = 0.84087(39) fm [3] resulted smaller by 7 standard deviations than the 2010-CODATA number. After several experiments and much discussion, in 2014 the CODATA Group gave the updated result R CO 14 E = 0.8414 (19) [4], stressing that the discrepancy among the different experiments is not understood, even within the atomic physics field. As a matter of fact, experiments with hydrogen atoms give different conclusions about the radius. A study of the 2S-4P transition [5] gives R H e 2S E = 0.8335(95) fm, compatible with muonic hydrogen, wheres a more recent measurement of the 1S-3S transition [6] concludes in the value R value for the radius R CLAS E 0.83 [1]. On the contrary, most results from electron-proton scattering converge, in general, towards a large value of the radius (R el E 0.88) with associated errors much larger than in atomic physics.
Assuming that the interaction occurs through the exchange of one photon of four momentum Q 2 = −q 2 , and neglecting the electron mass, i.e., in the so-called Born approximation, the expression of the reduced cross section in the proton rest frame is: where M is the proton mass, E and θ e are the initial energy and the scattering angle of the electron, dσ/d is the differential cross section and G E and G M are the proton electric and magnetic Sachs form factors (FFs). The reduced cross section σ red is linear in the variable = [1 + 2(1 + τ ) tan 2 (θ e /2)] −1 , being τ = Q 2 /(4M 2 ), with Q 2 = −4E E sin 2 (θ e /2), where E is the scattered electron energy. The kinematics of the scattering process is fully defined by two variables only. They can be chosen as the initial energy and the scattering angle of the electron or, more often, as and Q 2 .
Since the electromagnetic FFs G E and G M are functions of one variable only, Q 2 , the measurement of the reduced cross section at different scattering angles and fixed Q 2 allows to extract their squared values, as the slope and the intercept multiplied by τ , of a distribution which is linear in the variable . This method is called Rosenbluth separation [12]). At small (large) Q 2 values, the electric (magnetic) FF dominates because of the presence of the factor τ multiplying G 2 M . Another method may be used to extract the FFs: it consists in choosing pre-defined functions with FFs as parameters that are determined by a two-dimensional fit of the cross section. Such a method unavoidably pre-determines and constrains the Q 2 dependence of the FFs.
The root-mean-square charge radius of the proton is defined as being proportional to the limit for Q 2 → 0 + of the normalized first derivative of the electric FF, or equivalently, of its logarithmic derivative. For a recent discussion, see Ref. [13]. From the experimental point of view, recent experiments on electron-proton elastic scattering focussed to the largest precision on the cross section at the lowest transferred   [1]. Such a value was obtained by extracting the data on G E from the cross section under the hypothesis G M = 0. The extrapolation for Q 2 → 0 + was based on a pre-defined rational function. The electric FF G E is, indeed, the dominant contribution to the elastic cross section at such small transferred momenta squared. A previous experiment at Mainz [2] obtained R Mainz E = 0.879 ± 0.005 stat ± 0.004 syst fm, based on different analyses. The electric FF G E was extracted by means of two methods: pre-defined spline or polynomial functions, and the Rosenbluth separation. Even though it entails higher errors for the FFs, this second method can be considered as a modelindependent extraction of FFs from the cross section data since it is based on the one-photon exchange assumption only. Further evidence of the difficulty inherent to the extraction of the radius from the FF data, extracted in their turn from the cross section data, is represented by the significantly different values of R Mainz E and R CLAS E obtained by two sets of FF data that, especially in the low-Q 2 region, are indeed compatible, as it is shown in Fig. 1.
From the theoretical point of view, much discussion has been devoted to corrections to the low-Q 2 scattering as radiative corrections, Coulomb corrections, or effects due to nonvanishing electron mass. In some kinematical regions, they are of the order of the claimed precision on the cross section. Such corrections are, in general, model dependent (for a review see Ref. [14]). More interesting for the content of this work are the considerations about the validity of Eq. (2), and the extraction of the radius from the experimental data.
Reducing the transferred momentum squared, i.e., probing the region where Q 2 → 0 + , should lead to a smaller error in the extrapolation. However, if the extrapolation is done with the help of a predefined fit function, the form of the derivative is highly constrained.
Calculating instead the numerical derivative from the cross section data and fitting directly the derivative, Q 2 appears in the denominator: smaller values induce a larger error in the derivative itself. Precise experiments at smaller momentum transfer lead to a kind of fractal behavior of the electric FF derivative, not bringing any new information on the dimension of the proton. Whereas most of the works stress the importance of the precision of the data at low Q 2 , other works (see e.g. Refs. [15][16][17]) show that investigating the intermediate range of Q 2 and then focusing on theorydriven extrapolations would give more severe constraint to the radius. Note that in Refs. [15][16][17] the authors reached different conclusions about the size of the proton. Extending the method of Ref. [18] to the precise data recently obtained from JLab-PRad [1] and Mainz [2] experiments, we suggest to determine the radius directly from the derivative of the electric FF, which is numerically extracted from the cross section data. We show that, by construction, the error associated to the derivative is much larger (two orders of magnitude, at least) than the experimental error on the cross section. The choice of these two experiments is obvious as their data are the most precise and extend to the lowest Q 2 values ever reached.
Before going in the detail of the spreading of the values of the proton radius extracted from different experiments (as well as from the same data with different analysis), we would like to drive the attention of the reader to a physical issue. The definition given in Eq. (2) implies the limit Q 2 → 0 + , at which the elastic cross section diverges, being O Q 2 −2 . This limit is reached when, either the scattering angle θ e vanishes, i.e., the incident particle does not interact and hence is not deflected, or the scattered electron energy vanishes. In this case, the Coulomb interaction induces a capture process and the scattering formalism does not apply anymore. It follows that the determination of the proton radius requires an extrapolation of the scattering process and of the related formalism in the kinematical region where the electron is captured and an atom is formed. The description of the trapping process is highly model dependent, as it involves different types of corrections, as the Coulomb ones, that cannot be calculated exactly. Note that, in general, photo and electro-production reactions are studied separately, i.e., one can not describe experimental observables measured with a real photon as the static limit of electro-production models, the number of reaction amplitudes and the spin structure of the matrix elements being different. The question arises about the extrapolation of a measured quantity in order to obtain the derivative of this quantity. In other words, is it reliable to extrapolate data in regions where the approximation under which those data are obtained, fails?

Data samples
We consider three sets of data on the electric FF, labelled as P 1 , P 2 and A 1 1 , obtained by the JLab-PRad experiment [1] at E beam = 1.1 GeV and E beam = 2.2 GeV, and by the Mainz experiment [2] with the Rosenbluth method. The P 1 and P 2 data are more precise than those of the set A 1 extending also to a minimum Q 2 which is lower by an order of magnitude, as it is shown in the upper panel of Fig. 1.
In the light of their precision and density in Q 2 , one may think a priori that P 1 and P 2 data would constrain better the proton radius. As it can be seen in the left panel of Fig. 1, they show that apparently a plateau is reached in the low-Q 2 region, helping for the extrapolation at Q 2 = 0. Indeed, the plateau is a visual effect of the logarithmic scale in the Q 2 axis, as it is easily verified by considering the linear representation of the same data shown in right panel of Fig. 1. The black solid and a red dashed curves represent the expected behavior of the electric FF when assumed linear in Q 2 , with slopes corresponding to R black E = 0.83 fm and R red E = 0.87 fm, respectively. As pointed out in Ref. [17], it appears that data at larger Q 2 would better disentangle the two values of the radius.

Data on the discrete derivative
The data set of the N − 1 experimental values of the discrete derivative of the electric FF G E is obtained from the data set on the electric FF itself, through the following formulae [19] with k = 1, 2, . . . , N − 1. Only the first-order central finite difference has been considered for two reasons: 1 As will be discussed in Sect. 3, a set of N data point is an N -tuple of , i.e., the form factor and its derivative at each Q j value.
(a) depending only on two measured values of the FF, it gives the lowest errors on the discrete derivatives; (b) the non-homogeneous distribution of the experimental points in Q 2 makes the adoption of higher-order estimators unreliable.
The wide spreading and the large errors of the data on the discrete derivative of G E , shown in the left panels of Fig. 2, are not surprising. They can be understood by considering the function G E (Q 2 ), which approximates the FF in the limit Q 2 → 0 + up to terms of the order O (Q 2 ) 2 . Such an approximation is defined as: is the first derivative of the FF and can be interpreted as a kind of reconstructed FF from the derivative. In fact, by considering the Taylor series of G E (t) and G which, being the FFs analytic functions at t = 0, 2 converge in the same disk D R = {t : |t| < R} of radius R < 4M 2 π , we obtain the series of G E (t) as so that the difference between the two functions G E (t) and In the light of that, we can reconstruct the FF data from those on its discrete derivative and compare them to the original ones in order to visualize how much the errors on G E do influence its derivative, especially in the limit t → 0. Hence, starting from the data set D of Eq. (4), we compute the discrete derivative D of Eq. (3), through the formulae of Eq. (5), and then we define the set 2 The FFs are analytic in the whole q 2 = −Q 2 complex plane with a discontinuity cut (4M 2 π , ∞) over the positive real axis. with As it can be inferred from Eq. (6), in the limit Q 2 = t → 0 + , the points of the two data sets D andD, Eqs. (4) and (7), should describe the same quantity, i.e., the electric FF G E (Q 2 ). The crucial difference between the data points of these two sets is that in the first case the FF is the measured observable, while in the second, the measured quantity is its discrete derivative. The sets of data on G E and G E are shown in the right panels of Fig. 2 for the three considered experiments [1,2]. Two key points have to be highlighted.
• The wide spreading and large errors obtained for the data on the discrete derivative, G E,k with k = 1, 2, . . . , N − 1, shown on the left panels of Fig. 2, are drastically reduced in the observable G E , because the quantities G E,k appear in the expression for G E,k of Eq. (8), multiplied by the corresponding transferred momenta squared Q 2 k , whose values are very small.
• The fact that the data on the reconstructed electric FF G E have the same trend as those on G E but with orders of magnitude larger errors, does confirm that, in the limit Q 2 → 0 + , the accuracy of the first derivative extracted from the data is largely downgraded with respect to that of the original quantity, i.e., G E .
The latter statement is clearly proven from the results obtained when the root-mean-square charge radius is extracted directly from the data on the discrete derivative of the electric FF.

Procedure to extract the root-mean-square charge radius
We extract the root-mean-square charge radius of the proton, R E , from the three data sets P 1 , P 2 [1] and A 1 [2] on the discrete derivative of the electric FF G E , i.e., that have been obtained from the original measurements of G E using the formulae of Eq. (5). The procedure consists in fitting the data by exploiting two power series for the FF that differ by the expansion variable.
Electromagnetic FFs are functions of only one complex variable, q 2 , where q is the four-momentum of the virtual photon, and are analytic in the q 2 complex plane with the discontinuity cut (4M 2 π , ∞) along the positive real axis, namely the time-like region. The branch point at q 2 th = 4M 2 π , the socalled theoretical threshold, corresponds to the squared mass of the π + π − system that is the lightest hadronic state that  2 to have a direct interpretation as in terms of mean-square charge radius. Right panels: data on the original (solid circles) and reconstructed (empty circles) electric FF G E and G E . Data of the first, second and third rows are from the data sets P 1 , P 2 and A 1 , respectively couples to the photon, having the same quantum numbers. It follows that the FFs admit Taylor expansions in the variable q 2 centered at the origin q 2 = 0 and converging in the disk of radius R q 2 = q 2 th . On the other hand, since the discontinuity, represented by the branch cut, concerns only the imaginary part of the FF, the real part is a regular function of q 2 , that can be expanded in a Taylor series centered at the origin with a convergence radius larger than R q 2 . In the space-like region, q 2 < 0, where the FFs are real, such a Taylor expansion is particularly suitable to describe the data on both the FF and its first derivative.
We consider the power series with the sequences of coefficients {a k } ∞ k=0 , {b m } ∞ m=0 ⊂ R and where t = Q 2 = −q 2 = −s is the natural expansion variable, i.e., the opposite of the transferred momentum squared, while z(s) is a function that maps the edge of the discontinuity branch cut, lying on the s-complex plane, into the unit circle, see Appendix A for a detailed discussion. The definition of the k-th coefficient a k , with k = 0, 1, . . ., in terms of the 2k-mean-square radius, as well as its relationship with the root-mean-square charge radius R E , ruled by Schwarz's inequality, are given in Ref. [19]. The most suitable form of the function z(s) for studying the behavior of the FF in a neighborhood of the origin s = 0, where its derivative has to be extracted, follows from Eq. (A.1) and reads as a function of s = q 2 and t = Q 2 = −s. The power series for the first derivative of the FF with respect to the t variable, obtained from Eqs. (9) and (A.3), are Finally, from the definition given in Eq.
(2) and the FF power series of Eq. (9), the coefficients a 1 and b 1 can be related to the root-mean-square charge radius as It is worth noting that the extrapolation to zero of the numerical derivative, and hence the extraction of the radius, is independent of the global normalization that indeed entails the 0-th coefficients, that in this case are not present.

Extraction of the root-mean-square charge radius
The root-mean-square charge radius of the proton has been extracted from the data on the electric FF by considering various fitting schemes, which differ by the data sets and the fitting functions. Two data sets have been considered [1,2]: • the set P 1 ∪ P 2 ∪ A 1 containing the data on the electric FF G E ; • the set P 1 ∪ P 2 ∪ A 1 containing the data on the discrete derivative of the electric FF, obtained from those on G E through Eq. (3).
The fit functions are polynomials either in the variable t = Q 2 , the transferred momentum squared, or z(s) = z(−t), following the expression given in Eq. (10). They are obtained by truncating the power series of Eqs. (9) and (11) for the FF and its first derivative, respectively. The polynomials in the variable t of degrees n and n − 1, for the FF and its first derivative are and depend on the same sets of parameters apart from the 0 th coefficient a 0 , that is present only in G (n) The polynomials in the variable z(−t), which is explicitly indicated as the argument, are Due to the non-linear relation between t and z(−t), see Eq. (10), the functions G are not so in the variable z. Nevertheless, comparisons of the results obtained by using fit functions having the same degree n can be done, since they depend on the same number n of parameters.
By considering polynomials in both variables t and z of 12 different degrees, from n = 2 up to n = 13, and two sets of data, for the FF and its discrete derivative, we obtained: 2×12×2 = 48 values for the root-mean-square charge radius of the proton. The results obtained using polynomials in the t and z variable are reported in Tables 1 and 2, respectively. We indicate with R n,t (z) and R n,t (z) the radii extracted by the fit to the data on the FF and its discrete derivative performed with the fit functions G (n) As an example, Fig. 3 shows the parameterizations obtained by performing the fit only to the FF data, i.e., by minimizing the χ 2 's where the index j runs over all the data points Q 2 j , G E, j , δG E, j ∈ P 1 ∪P 2 ∪A 1 and the fit functions are polynomials of the 7 th degree in the variable t and z(−t), depending on the parameter sets {a k } 7 k=0 and {b m } 7 m=0 , respectively.   Figure 4 shows the parameterizations for the FF (left panel) and its discrete derivative (right panel) obtained by fitting only the data corresponding to the latter quantity, where the χ 2 's are In this case, the index j runs over all the data points Q 2 j , G E, j , δ G E, j ∈ P 1 ∪ P 2 ∪ A 1 . The fit function G (7) E (t), defined in the second expression of Eq. (12), is still Fig. 3 Left panel: data sets P 1 (red circles), P 2 (black squares) [1] and A 1 (pink diamonds) [2] on the electric FF, the blue and red bands represent the parameterizations G E,7 (t) and G E,7 (z), respectively. Right panel: data sets P 1 (red circles), P 2 (black squares) [1] and A 1 (pink diamonds) [2] on the first discrete derivative multiplied by −6 and in units of fm 2 , the blue and red dashed bands represent the parameterizations G E,7 (t) and G E,7 (z), respectively. The fit has been performed to the data on the electric FF G E corresponding to all sets P 1 , P 2 and A 1 a polynomial, even though it is of the 6th degree, because of the derivative in t. On the other hands, G (7) E (z), given in the second expression of Eq. (13), due to the derivative with respect to the variable t, looses its polynomial structure in z. Nevertheless, both functions G (7) E (t) and G (7) E (z) depend on 7 parameters, namely those of the sets {a k } 7 k=1 and {b m } 7 m=1 , respectively. Figure 5 shows the values of the root-mean-square charge radius and the normalized χ 2 from the fit of FF data only, with parameterizations in form of t and z(−t) polynomials, G  12) and (13), as functions of the polynomial degree n. Figure 6 shows our main result, i.e., the values of the root-mean-square charge radius, as well as the normalized χ 2 obtained from the fit of the discrete derivative of the FF, namely the data of the set P 1 ∪ P 2 ∪ A 1 , with the functions G  12) and (13).
In both figures the green band shows the root-mean-square charge radius R CLAS E = 0.831 ± 0.007 ± 0.012 fm, extracted by the JLab-PRad experiment [1] from the data sets P 1 and P 2 .
Four facts have to be noted. , are more stable with respect to the degree n, showing slightly larger errors. Fact 3. Above the degree n = 5, the results stabilize, as also proven by the convergence of the χ 2 's, and the root- Fig. 4 The graphs are the same of those of Fig. 3, but in this case the fit has been performed to the data on the discrete derivative of the electric FF, namely on those of the sets P 1 , P 2 and A 1  mean-square radii extracted by fitting to the data on the discrete derivative agree with those obtained by fitting to the data of the the FF, shown in Fig. 7 as dashed and solid bands, respectively. This is a clear confirmation that extracting the root-mean-square charge radius from the derivative of the fit function instead of the fit of the discrete derivative does necessarily entail an important underestimate of the uncertainty. Fact 4. The result shown in Fig. 5 gives an unequivocal demonstration of the pivotal role played by the higher-Q 2 precise data set A 1 from the Mainz experiment [2]. This is proved by the disagreement between the value that we obtained by considering the data set P 1 ∪ P 2 ∪ A 1 , blue and red dashed bands, and the one from the  JLab-PRad experiment [1], solid green band, using the data set P 1 ∪ P 2 (discarding the higher-Q 2 measurement from the Mainz experiment [2]).

Conclusions
Recent and precise data on elastic electron scattering on the proton at very low transferred momentum were collected for the extraction of the proton radius. The results from the JLab-CLAS collaboration Ref. [1] are critically discussed and compared with the Rosenbluth set of data extracted from the Mainz experiment [2] spanning a large range of low Q 2 values. The choice of the Rosenbluth set among several extractions is justified by the work of Ref. [19], showing that FF data extracted by pre-imposed polynomial fits have a smaller statistic error but the result on the radius is highly constrained leading to a large systematic error, difficult to be evaluated. In the case of Ref. [1], the magnetic FF was neglected. Therefore, the pre-definition of a fitting function was not necessary, and G E was directly related to the measured observable, the cross section. The proton radius is related to the derivative of the cross section. By calculating and plotting the numerical derivative point by point, one can observe that measurement at the smallest Q 2 give the largest error on the derivative. The value and the error on the Fig. 7 Summary of Figs. 5 and 6. Root-mean-square charge radii: R n,t (blue band), R n,t (blue dashed band), R n,z (red band) and R n,z (red dashed band). The green band represents the result from the JLab-PRad experiment [1] extracted radius depend on the Q 2 range and on the fitting function. The Mainz data [2] lead to the smaller error on the radius.
The main finding of this work is that all published values of the proton radius extracted from elastic electron-proton scattering suffer from under-evaluation of the associated error. This fact has two main sources: • Numerical: the extrapolation of the electric form factor based on any pre-defined function, in a definite range, gives a strong constrain to the derivative and its limit to the static point; • Physical: elastic scattering is a dynamical two-body problem. Its extrapolation to the photon point, i.e., a static limit corresponding to a compound object, is outside the range of validity of the formalism. The extrapolation from the scattering cross section becomes hazardous.
We have shown that a sub-selection of the Q 2 range, a choice of possible functions and of the FF data set may lead to extract any value in the range 0.7-0.9 fm. Our method, based on a direct extrapolation of the derivative built from the data themselves, shows that a large error should be attributed to any number for the radius. Such error is sufficiently large to bring consistency among all measurements. We note that a similar situation has been going on for several years concerning the determination of the electric and magnetic form factors of the proton at large trans-ferred momenta, where polarized and unpolarized experiments found inconsistent results. In this case, too, the observable i.e., the electric form factor appears as the derivative of the elastic cross section in the Rosenbluth plot, leading to an underestimation of its error in unpolarized ep experiments [20]. In the present case, the analysis with numerical derivatives give larger errors that make all results consistent among each other and consistent with atomic and electron scattering experiments.
The contradiction among the extracted values of the proton radius is due therefore to the attempt of pushing the information gathered from experiments beyond the limits imposed by the physical observables. In the present case, the observable is the cross section, but the physical quantity is proportional to its derivative. Furthermore, the extrapolation of this quantity outside the validity domain of the method (the scattering formalism) should be taken with caution.
In the light of all these reasons, for all intents and purposes, this conclusion is general and applies to experiments and analyses not explicitly considered in this paper. It depends neither on the functional form of the fitting function (different classes of functions: polynomials in t and z, spline in Ref. [14] and references therein) nor on normalization issues, as the radius by definition is proportional to a logarithmic derivative.

Methodology
Our analysis follows these steps: • We consider only published data on the electric form factor. We focus on two sets of data: (1) Ref. [1], recently obtained from the CLAS11 collaboration at JLab (these are the most precise data, at the smallest values of transferred momenta); (2) the set of data from Mainz [2], less precise and extending to larger Q 2 . • We calculate the numerical derivative of the measured electric form factor. • We apply a number of fits for the extrapolation of the numerical derivative (that is, indeed, the variable related to the physical quantity, the radius) to Q 2 → 0, with different functions: namely different order polynomials, functions of the variables t and z in the complex plane.
• the values of the radius and their errors are derived.
• the results show that the central values span a large interval, statistically compatible taking into account the errors.
Funding Open access funding provided by Università degli Studi di Perugia within the CRUI-CARE Agreement.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study based on data that are already published, produced by other independent groups, and hence these data will not be deposited.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. of the time-like region ranging from the origin up to the theoretical threshold, are mapped onto those of the segment (−1, 1) in the z complex plane, as it is shown Fig. 8, where the correspondence s ↔ z is highlighted with colors. In particular, by considering as a reference point s = t 0 , the mapping through the function z(s) of the subsets −∞ < s < t 0 and t 0 < s < q 2 th is s ∈ (−∞, t 0 ] ∪ (t 0 , q 2 th ) ↓ z(s) ∈ z(t 0 ) = 0, z(−∞) = 1 ∪ z(q 2 th ) = −1, z(t 0 ) .
The value of the parameter t 0 ∈ (−∞, q 2 th ) has to be set equal to the four momentum squared under investigation, indeed, being mapped onto the origin z(t 0 ) = 0, such a value guarantees a better convergence rate for the series in the z variable. In our case, aiming to extract the first derivative of the FF at Q 2 = 0, we set t 0 = 0, so that z(s) = Each k-th term of the series, with k ≥ 1, has a pole of order 2k at z = 1, that can be expanded in a further power series corresponding to the (2k − 1)-th derivative of a geometric series, i.e., In the light of that, the original series becomes where the last member is the formal expression of the power series in the variable z. The m-th coefficient b m , with m ∈ N, is then obtained as It is a combination of the m coefficients a j , with j = 1, 2, . . . , m, as an example the first four coefficients are b 0 = a 0 , b 1 = 4q 2 th a 1 , b 2 = 8q 2 th a 1 + 2q 2 th a 2 , b 3 = 4q 2 th 3a 1 + 16q 2 th a 2 + 16(q 2 th ) 2 a 3 . (A. 2) The first derivative of the FF with respect to the t variable is still analytic in a neighbourhood of the origin and hence it is expansible in power series. The series can be obtained from those of Eq.