Hadron and light nucleus radii from electron scattering

"How big is it?" and "How heavy is it?" are two questions posed within all fields of life and endeavour, with meanings tuned to the objects under consideration. In physics, they typically address signature Poincaré-invariant characteristics of a particle or bound state. We focus on systems governed by strong interactions within the Standard Model of particle physics (SM), which are supposed to be described by quantum chromodynamics (QCD) [1]; namely, hadrons – baryons (proton, neutron, etc.), mesons (pion, kaon, etc.), and light nuclei.

The QCD Lagrangian is built upon the non-Abelian colour gauge group, SU$ _{{\rm{c}}}(3) $, and expressed in terms of gluon and quark fields. Gluons transform as the adjoint representation of SU$ _{{\rm{c}}}(3) $ and quarks form the fundamental representation. This means that gluons and quarks all carry a colour quantum number, a generalisation of the dynamical charge in quantum electrodynamics (QED). Empirically, colour is confined. Consequently, whilst these gluons and quarks provide the basis for developing a rigorous QCD perturbation theory [2], which is crucial in the analysis of high-energy processes, they are essentially different from the quanta of QED. The gluons and quarks used to express the QCD Lagrangian are not the objects that propagate to be captured in detectors. Instead, only composite bound-states, seeded by colour singlet combinations of QCD quanta, are directly measurable.

The proton is a prime example, as highlighted by Fig. 1. Perceived in terms of QCD's Lagrangian degrees-of-freedom, it is a bound-state seeded by a colour singlet three-quark combination (two u quarks and one d quark) that is embedded in a complex medium produced by infinitely many gluons and quark+antiquark pairs that is also colour neutral overall [3]. With so much stored within, it is natural to ask after the size of this highly nontrivial arrangement of quanta.

Figure 1.  (color online) In terms of QCD's Lagrangian quanta, the proton contains two valence up (u) quarks and one valence down (d) quark, the combination of which is a colour singlet. Owing to the character of quantum field theory and the nature of strong interactions, the proton also contains infinitely many gluons and sea quarks, drawn here as "springs" and closed loops, respectively. (The neutron looks similar, except that it is seeded by two d quarks and one u quark.)

The natural scale for the size of atomic systems is the Angstrom, Å. Using natural units ($ \hbar = 1 = c) $, this scale is set by $ 1/[m_e \alpha_{{\rm{em}}}] $, where $ m_e\approx 0.511 $MeV is the electron mass and $ \alpha_{{\rm{em}}}\approx 1/137 $ is the QED fine structure constant. The quarks in QCD's Lagrangian have current masses, $ m_{u,d} $, generated by the Higgs boson, that are on the same scale as $ m_e $ [4, Sec. 59] and the electromagnetic interactions of quarks are also characterised by $ \alpha_{{\rm{em}}} $. However, if the size of the proton in Fig. 1 were measured on the Å scale, then there would be no notion of confinement because modern detectors have the capacity to directly image targets of this size [5].

Of course, the proton is very much smaller than the hydrogen atom; in fact, it provides the nucleus of $ ^1 $H, whose scale is measured in femtometres: 1 fm = $ 10^{-5} $Å. Hence, all the activity in Fig. 1 is contained within a volume on the order of 1 fm$ ^3 $. It is logical to ask the question: Why is the proton so small?

At this point, another issue arises. The proton mass is known to a high level of precision, with a relative error of $ \sim 10^{-8} $ [4]: $m_p = 0.9382720813(58) {{\rm{GeV}}} = 1/[0.210309 {{\rm{fm}}}]$. This value is commensurate with the associated natural length scale; but, surprisingly, it is more than 100-times larger than the sum of the typically quoted current-masses of the proton's valence quarks [4, Sec. 59]. One must then also ask: Why is the proton so heavy?

The proton is small and the proton is heavy. Within the SM, these features are expected to be explained by QCD. However, this anticipation raises yet another issue. Given that $ m_p $ is so much larger than the quark current-masses generated by the Higgs boson, then the bulk of the proton mass should be understandable in the absence of the Higgs mechanism of mass generation; namely, within QCD defined without quark couplings to the Higgs. The problem with this reasoning is that the associated Lagrangian does not possess any mass-scales: without quark current-masses, the chromodynamics Lagrangian is scale invariant. This is readily verified and raises the key, unifying question [6]: How does QCD itself generate the single mass-scale that characterises the proton mass, the proton size – which may be reckoned as the characteristic confinement length, and, in fact, the natural scale for all nuclear physics? Approved and planned experiments and facilities are targeting this puzzle, as described in Refs. [7–15] and related theory analyses [16–18].

The question can be posed as the need to understand the phenomenon of emergent hadron mass (EHM). To achieve that goal, a solution of QCD will be necessary. (It will only be sufficient if QCD is the true theory of strong interactions.) Nonperturbative continuum methods and numerical simulations of lattice-regularised QCD are together providing plausible solution paths [19], although nothing yet approaches the level of mathematical rigour specified in the formulation of the Yang-Mills Millennium Problem [20, pp. 129–152]. Notwithstanding the progress made, EHM, confinement, and possible connections between them, remain fundamental open questions. Not least because confinement underpins proton stability: unlike all other bound states, no proton has ever been seen to decay in the roughly 14-billion years since the Big Bang. In fact, the current lower limit on the proton lifetime is $ 1.6\times 10^{25} $ billion-years [21].

Particle masses can be measured in a variety of ways, e.g., considering long-lived systems with electric charge, like the proton, a Penning trap can serve as a reliable mass spectrometer [22]. Many hadron masses are measured with precision. However, masses are volume-averaged quantities, sensitive to the underlying dynamics but not a keen discriminator between theories. For instance, in strong interaction phenomenology, many widely differing Hamiltonians can produce hadron spectra that are equivalent so far as modern experiments can tell; e.g., consider and compare Refs. [23–28].

On the other hand, as highlighted above, radii reveal more about the underlying theory; but relative to masses, it is far more difficult to formulate and complete the reliable calculation of a radius within a given theory and particle radii are harder to measure. Again using the proton as an example, discovery of its anomalous magnetic moment, $ \kappa_p\approx 2 $, in 1933 [29] showed that the proton is not an elementary Dirac fermion. However, the first measurement of its electric radius was not made until the mid-1950s using electron+proton ($ ep $) elastic scattering [30, Table II]: $ r_p=0.77\pm 0.10 $fm.

Electron+hadron scattering remains an excellent radius measurement tool because QED and the electron (probe) are well understood and analyses involving isolated targets avoid the need to account for binding corrections as in bound-state calculations [31]. However, it does have issues. In order to extract a radius, one must measure the relevant form factor, $ F(Q^2) $, with high precision on a material domain of squared-momentum transfer, $ Q^2 $, that reaches deeply toward $ Q^2=0 $ because for systems with $ F(0)\neq 0 $, the radius is obtained via

Achieving the precision, reach to $ Q^2 \approx 0 $, and density of $ Q^2 $-coverage necessary for a reliable radius determination is very challenging.

Naturally, an object's mass and its radius are correlated; and any (putative) SM solution will deliver values for both. Hence, precise knowledge of an array of hadron and light-nuclei masses and radii are necessary to set rigorous benchmarks for theory. Herein, therefore, we describe a mathematical technique, applicable in the same form to a wide variety of measurements, and also calculations, that when applied to high-precision, dense data stretching to very low momentum transfers will return an objective result for the associated charge radius. Commonly described as the statistical Schlessinger point method (SPM) [32–35], but equally well as a statistical multi-point Padé approximant scheme, the approach avoids any specific choice of fitting function in analysing data.

In 2010, the Particle Data Group (PDG) listed the proton charge radius as [36]: $ r_p = 0.8768(79) $ fm. However, controversy was brewing following an extraction of the radius from measurements of the Lamb shift in muonic hydrogen ($ \mu H $) [37]: the result $ r_p= $ 0.84184(67) fm is 5.0 standard deviations ($ 5.0\sigma $) smaller than the then accepted value. As highlighted by Fig. 2 – Rows [A – I], this disagreement stimulated many new experiments, with conflicting results; so, today, we are dealing with the "proton radius puzzle" [51].

Figure 2.  (color online) Proton electric radius, $ r_p $, extractions, various techniques: [A] [50]; [B] [37]; [C] [38]; [D] [39]; [E] [40]; [F] [41]; [G] [42]; [H] [43]; [I] [44] [L] [45]; [M] [46]; and [N] [47]. The green-hashed rectangles display the anticipated precision of two forthcoming muon+proton scattering experiments: [J] [48,49]; and [K] [7]. The two vertical grey bands indicate the separate error-weighted averages of the small-radius and large-radius extractions.

Although, the proton's mass can now be calculated with fair accuracy using modern theoretical tools [25, 52, 53], the same is not true of $ r_p $. This deficiency has left room for many possible solutions of the puzzle to be offered, e.g. [54, Sec. 5]: unknown QCD-related corrections may have been omitted in the $ \mu H $ analysis, whose inclusion might restore agreement with the electron-based experiments that give a larger value; new interaction(s) or particle(s) outside the SM may 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 $ ep $ scattering.

Results from the most recent $ ep $ scattering experiment were released in 2019 [44, PRad]:

This was the first published $ ep $ scattering analysis to obtain a radius in agreement with that extracted from $ \mu H $ measurements: see Fig. 2 – Row I. Notably, the PRad experiment implemented some significant improvements over earlier efforts, inter alia: stretching to $Q^2 = 2.1\times 10^{-4}$GeV$ ^2 $, the deepest reach yet achieved; and simultaneously covering an extensive low-$ Q^2 $ domain, viz. $ 2.1\times 10^{-4} \leq Q^2/{{\rm{GeV}}}^2\leq 6 \times 10^{-2} $. Furthermore, since the charge radius is obtained using Eq. (1), with $ F=G_E^p $, the proton elastic electromagnetic form factor, the PRad analysis paid careful attention to the effect of data fitting function on the extracted charge radius. The prejudice such choices may introduce has been highlighted [55–63].

PRad selected an optimal function form by using a bootstrap procedure applied to pseudodata produced with fluctuations mimicking the $ Q^2 $-binning and statistical uncertainty of their experimental setup. No knowledge of actual PRad data was introduced. This procedure renders the PRad extraction robust; but it also means that, ultimately, a specific function form was chosen [64].

Of particular note is the conflict between the PRad result for $ r_p $ and the values obtained in analyses of all other $ ep $ scattering data, Fig. 2 – Rows [A, C, G]. In this connection, the data described in [50, A1] may be chosen as the benchmark, owing to their reach toward $ Q^2=0 $, high-precision, and sheer quantity. A comparative analysis of the PRad and A1 data therefore provides an ideal means by which to illustrate the strengths of the SPM.

Suppose one has N pairs, ${\mathsf D} = \{(x_i,y_i=f(x_i)), ~ i= 1,\ldots,N \}$, being samples of some analytic function, $ f(x) $, at a given set of enumerable points. In applying the SPM, one constructs a continued-fraction interpolation:

whose coefficients, $ \{a_i|i=1,\ldots, N-1\} $, are computed recursively and guarantee $ {\cal{C}}_N(x_i) = f(x_i) $, $ i=1 \ldots,N $. Mathematically, the method provides an accurate reconstruction of $ f(x) $ within a radius of convergence determined by that one of the function's branch points which lies closest to the real-axis domain containing the data sample. For instance, consider a monopole form factor sampled precisely at $ N>0 $ points, then using any one of those points, the SPM will truly reproduce the function. The continued fraction interpolation encapsulates both local and global aspects of the underlying function. Here, the global quality is crucial: it validates use of the constructed interpolation outside the domain of available samples and thus legitimises use of Eq. (1), i.e., evaluation of the curves' slope at the origin.

In physical cases, one is presented with N data that are statistically distributed around some curve which must be accurately reconstructed. So long as N is large, then one can introduce a powerful statistical aspect to the SPM. To wit, one selects at random $ M < N $ points from the set $ \mathsf{D} $. In typical applications, $ 6\lesssim M \lesssim N/2 $. Then, in theory, one can build $ C(N,M)=N!/[M!(N-M)!] $ – a huge number – of distinct interpolating functions. In practice, the useable number is reduced somewhat by imposing physical constraints on their behaviour. Herein, we consider electric and magnetic form factors at spacelike momenta. Sampled on this domain, such functions must be continuously differentiable. This is the restriction we implement. It is a weak constraint, so the useable number of interpolants remains very large.

Using PRad data obtained at the $ 1.1 $GeV electron beam energy as an example: $ N=33 $; so $ M\in [6,17] $ may be chosen, giving $ {\cal{O}}(10^6 - 10^9) $ possible interpolators; and one can productively work with the first $ 5 000 $ that are continuously differentiable monotonic functions on the entire range of $ Q^2 $ data. Importantly, no other constraints need be applied; especially, no unity constraint on $ G_E^p(Q^2=0) $.

Rather than present a dry mathematical account, we continue with this concrete descriptive illustration. An extrapolation to $ Q^2=0 $ is defined by each interpolating function, using which $ r_p $ can be calculated via Eq. (1). Then, for every M, a value of $ r_p $ follows as the average of all results obtained from the $ 5 000 $ curves.

In completing the analysis, one must also account for the experimental errors in the given data set. This can be done by using a statistical bootstrap process. Namely, one generates, e.g., $ 1 000 $ replicas of each set by substituting each datum by a new point, Gaussian distributed at random around a mean specified by the datum itself with variance equal to its error. This leads to a probability distribution function, $ {\cal{N}}(\mu,\sigma) $, characterising the $ r_p $ values that is, to an excellent approximation, normal, with average $ r^{{ } M}_p $ and standard deviation $ \sigma=\sigma_{r}^{{ } M} $. (If a constant distribution is used instead, the mean is unchanged and the uncertainty is reduced. The Gaussian distribution therefore provides a more conservative uncertainty estimate.)

Further, the fact that M is deliberately varied leads to a second uncertainty source, $ \sigma_{{ }{\delta M}} $, that can be estimated by changing $ M\to M' $, repeating the procedure just described for this new $ M' $-value, and evaluating the standard deviation of the $ r^{{ } M}_p $ distribution for different M values.

Combining these considerations, then for a given data set, the SPM result is

For the PRad $ 1.1 $GeV beam example, one computes radii for each of the values $\{M_j=5+j \vert\ j = 1,\dots,n_{{ }{M}};~ n_{{ }{M}} = 12\}$, so that for any given data set there are 60-million values of $ r_p $, every one calculated from an independent interpolation. In this case, $ \sigma_{{ }{\delta M}}\ll \sigma_{r}^{{ } M_j} $ for all js in the range specified above [47, Supplement].

Before applying the SPM to any set of real data, one final issue must be addressed. Namely, in ideal situations, experimental data are statistically distributed around that curve which properly represents the observable. Since they do not lie on the curve, empirical data should not be directly interpolated. This feature can be handled by smoothing with a roughness penalty, that may be implemented following the ALGOL procedure detailed in Ref. [65] and made concrete in Ref. [47, Sec. 3]. The approach is characterised by a rigorously defined optimal roughness penalty, $ \epsilon $. Its value is a self-consistent output of the smoothing procedure: $ \epsilon =0 $ means the data are untouched by smoothing, whereas $ \epsilon =1 $ returns a linear least-squares realignment of the data. In the PRad example, $ \epsilon \simeq 0 $.

Distilling the preceding remarks, the SPM extraction of a radius from any high-quality N-element set of electron+hadron scattering data usually proceeds as enumerated in Table 1.

(i)	Generate 1 000 replicas for the given experimental central values and uncertainties.
(ii)	Smooth each replica using the associated optimal parameter, $ \epsilon $;
(iii)	For each number of input points M, chosen in the range $ 6\lesssim M \lesssim N/2 $, determine the distribution of proton radii $ r^{{ }{M}}_p $, associated $ \sigma^{{ }{M}}_r $, and overall $ \sigma_{{ }{\delta M}} $;
(iv)	Combine this material to calculate the final result for the radius and (statistical) uncertainty through Eq. (4).

Table 1.  The SPM extraction of a radius from any high-quality N-element set of electron+hadron scattering data usually proceeds in these four steps.

With any such method, it is fair to ask whether the scheme is robust, i.e., whether it can reliably extract a known value in a diverse array of relevant cases. For the proton radius case, that was checked using a wide variety of models employed over time by various teams to fit the world's $ ep $ scattering data [66–72] and therefrom generating a proton electromagnetic form factor $ G_E^p $ with a known value for the radius. Using those models, Ref. [47] generated replicas with the binning in $ Q^2 $ and errors of the A1 [50] and PRad [44] data. In all cases, irrespective of the $ G_E^p $ model employed, the SPM was shown to return the radius value used to generate the pseudodata, with results being practically independent of M, the number of initial input points. In Sec. VII, this standard validation procedure is detailed in connection with a proposal to extract the deuteron charge radius [73].

With the method defined, one may ask what the SPM can reveal about the discrepancy between the $ ep $ scattering results for $ r_p $: Fig. 2 – [A, C] vs. [I]. The first test is a comparison between the value of $ r_p $ quoted by the PRad Collaboration [I] and the SPM.

There are $ N=33 $ points of PRad $ 1.1 $ GeV beam energy data. Analysed as described in Sec. III, one finds [47]:

With the $ 2.2 $ GeV beam energy, PRad collected $ N=38 $ data, using which the SPM yields [47]:

Evidently, treated independently using the SPM, 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.

Combining the data from both beam energies [47]:

which is displayed in Fig. 2 – [N]. In this result, the smaller error determined in analysing the beam-energy $ 2.2 $ GeV data has driven the overall uncertainty down to roughly 50% of that in Eq. (5a). These observations accord with those made by the PRad Collaboration [44, Supplement – Fig. S16]. As highlighted by Fig. 2 – [I] cf. [N], within mutual uncertainties, Eq. (6) reproduces the published PRad result.

A natural comparison for PRad data is provided by the measurements reported in [50, A1]. They encompass $ 1 400 $ cross-sections, recorded at beam energies $ 0.18 $, $ 0.315 $, $ 0.45 $, $ 0.585 $, $ 0.72 $, $ 0.855 $ GeV, and stretch to $ Q^2=3.8 \times 10^{-3} $ GeV$ ^2 $, which is deep but still roughly 20-times larger than the smallest PRad value.

Considering first the low-$ Q^2 $ region, consisting of $ N=40 $ data in the interval $ 0.0038 \leq Q^2/{{\rm{GeV}}}^2 \leq 0.014 $, the SPM returns [47]:

Expanding the selection to include all available A1 data, the SPM yields the same central value but a larger uncertainty: $ r^{{\mathrm{A1}}}_p/{{\rm{fm}}}=0.857\pm0.021_{\mathrm{stat}} $. In this instance, $ \sigma_{{ }{\delta M}}\sim \sigma_{r}^{{ } M_j} $. Thus, extending the range up to $ Q^2\sim 1 $GeV$ ^2 $ merely introduces a perceptible sensitivity to the number of data sampling points, M. Now, whereas the original $ r_p $ estimate from the A1 Collaboration is [50]: $ r^{\mathrm{A1}-\mathrm{coll.}}_p/{{\rm{fm}}} = 0.879\pm0.005_\mathrm{stat}\pm0.006_\mathrm{syst} $, the function-form independent SPM reanalysis of the A1 data, Eq. (7), reconciles the PRad and A1 values and shows both to be consistent with the $ \mu H $ experiments.

Since they are now consistent, it is sensible to combine the PRad and A1 results [Eqs. (6) and (7)]:

This value is the lowest entry in Fig. 2 – [N]. Evidently, the SPM, which produces form-unbiased interpolations of data as the basis for a well-constrained extrapolation, reveals that there is no discrepancy between the proton radius obtained from $ ep $ scattering and that inferred from the Lamb shift in muonic hydrogen – $ r_p = $ 0.84136(39) fm [37, 39], the modern measurement of the $ 2S \rightarrow 4P $ transition-frequency in regular hydrogen – $ r_p = 0.8335(95) $fm [40], and the Lamb shift in atomic hydrogen – $ r_p = 0.833(10) $fm [43]. These values also, therefore, match the combination of the latest measurements of the $ 1S\rightarrow 3S $ and $ 1S\rightarrow 2S $ transition frequencies in atomic hydrogen – $ r_p=0.8482(38) $fm [74] and even the muonic deuterium determination $ r_p=0.8356(20) $fm [75]. A recent assessment of proton radii extractions settles on the following value [4]:

This illustration shows that the proton radius is not a puzzle. Instead, it highlights that a reliable extraction of $ r_p $ from $ ep $ scattering experiments is possible, so long as one works with precise data, densely packed at very-low-$ Q^2 $, and adopts an analysis scheme that eliminates systematic error introduced by the use of specific, limiting choices for the functions employed to interpolate and extrapolate the data.

The proton is a relativistic composite object, so its electromagnetic interactions involve two form factors, Dirac – $ F_1 $ and Pauli – $ F_2 $ [30], in terms of which the electric and magnetic form factors are

Three-dimensional Fourier transforms of $ G_{E,M} $ were long interpreted as measures of the spatial distributions of electric charge and magnetisation inside the proton [76]. Today, their interpretation has changed [77]; but as emphasised above, the importance of measuring these key dynamical characteristics of the proton is undiminished.

Notably, for a structureless, noninteracting fermion, the Dirac equation entails $ F_2\equiv 0 \Rightarrow G_E\equiv G_M $. Perturbative corrections in QED only introduce small modifications for elementary fermions [78, 79]; hence, any material differences between $ G_E $ and $ G_M $ are marks of compositeness. As remarked in Sec. I, the first signal of this for the proton was the discovery of its anomalous magnetic moment $ \kappa_p :=\mu_p-1:= F_2(0) \approx 2 $ [29]. Nevertheless, as recently as twenty years ago, available $ ep $ scattering data were consistent with $ \mu_p G_E (Q^2)/G_M(Q^2) = 1 $. Only with the operation of an electron accelerator combining high energy, luminosity, and beam polarisation was it revealed that $ \mu_p G_E(Q^2)/G_M(Q^2) \neq 1 $ [79]. It follows that proton electric and magnetic radii must be different; unless some mechanism leads dynamically to

Eq. (11) is curious because if true, then the proton's Pauli and Dirac radii are not independent observables. Instead, their difference is a positive number fixed by the proton Foldy term [80]. No symmetry constraints entail Eq. (11). Thus, the result may only emerge as a consequence of SM dynamics; but the precision of today's theory is insufficient to answer the question of whether it does, or not.

As emphasised by Fig. 3, the empirical situation is also unclear. The most precise data suitable for extraction of $ r_p^{{\rm{mag}}} $ were obtained by the A1 Collaboration [50, 70]. However, entries [B, C] in Fig. 3 highlight that when using precisely the same methods to analyse the world's data, markedly different values for $ r_p^{{\rm{mag}}} $ can be obtained, depending on whether the A1 set is included or excluded. Facing this conflict, the PDG [4] quotes an average of [B, C] in Fig. 3:

Figure 3.  (color online) Upper frame. Proton magnetic radius extractions, various techniques: [A] [50, 70] — A1 Collaboration; [B] [81] — A1 data; [C] [81] — world average omitting A1 data; [D] [82] — dispersion theory; [E] [83] — dispersion theory; [F] [84] — lattice QCD. Light-grey band — Eq. (12); within mutual uncertainties, this value agrees with $ r_p $ in Eq. (9), which is indicated by the thin vertical blue band. Lower frame. SPM result in Eq. (13), obtained from A1 data [85].

and given this controversy, it is reasonable to ask what the statistical SPM can objectively reveal about the proton magnetic radius from extant data [85].

A1 data relate simultaneously to proton electric and magnetic radii. The SPM result for the electric radius is discussed in connection with Eq. (7). Regarding $ r_p^{{\rm{mag}}} $, Ref. [70, Supplement] lists $ 77 $ points obtained via a Rosenbluth separation [86] of the $ ep $ cross-sections into values for $ G_{E,M} $ on $ 0.015 \leq Q^2/{{\rm{GeV}}}^2 \leq 0.55 $. The nine lowest-$ Q^2 ~ G_M $ data are widely scattered. This is an issue with the Rosenbluth technique when applied to $ G_M $ at low-$ Q^2 $. Included in a SPM analysis, these nine points add noise without affecting the central value. Thus, Ref. [85] reported a study based on the $ N=68~ G_M $ points on $ {\cal{D}}_{A1} = \{Q^2 | 0.027 \leq Q^2/{{\rm{GeV}}}^2 \leq 0.55 \}$.

The SPM analysis of $ r_P^{{\rm{mag}}} $ follows the pattern employed for $ r_p $. In this case: $ N=68 $ and $ M\in [7,19] $. Consequently, there are O$ (10^9 - 10^{15}) $ potential interpolators; and steps (i)–(iv) in Table 1 are implemented within this space [85]. As with $ r_p $, the optimal smoothing parameter $ \epsilon \simeq 0 $ in all $ r_p^{{\rm{mag}}} $-related cases.

An extensive set of SPM robustness tests in connection with A1 $ G_M^p $ data revealed [85] that for $ M \in {\cal{O}}_M= \{ M_{ji} = 3 + 4 j + i \vert j=1,2,3,i=1,2,3\}$ one does not recover a satisfactory Gaussian distribution of radii. (This contrasts with the $ r_p $ case, where all values of M are associated with true normal distributions.) Consequently, only $ M\in {\cal{S}}_M=\{M_j=3+4j \vert\ j = 1,2,3,4\} $ were employed, leading to $ 20 $–million results for $ r_p^{{\rm{mag}}} $, each computed from an independent interpolation; hence,

As apparent in Fig. 3, within mutual uncertainties, this result agrees with that reported in Ref. [50], although the central value is 5.1% larger. In comparison, the analogous SPM extraction of $ r_p $ from A1 data, Eq. (7), returns a central result that is 2.6% lower than originally reported [50].

Now write

Eq. (11) is recovered with $ \delta_P=0 $. Rigorous contact may be made between quantum field theory and quantum mechanics when one views the Poincaré-invariant Dirac and Pauli form factors from a light-front perspective [87, 88]. Then [77], $ F_1^\prime(0) $ relates to the light-front transverse mean-square proton electric radius, $ b_E^\perp $, and $ F_2^\prime(0) $ to the analogous magnetic radius, $ b_M^\perp $. From this standpoint, Eq. (14) states that $ b_M^\perp>b_E^\perp $ so long as $ \delta_P>-1 $.

Using the radii discussed above, the pairings Eqs. (7, 13), Eqs. (9,13), Eqs. (9,12) yield the results drawn in Fig. 4:

Figure 4.  (color online) Experiment-based results for deviations ($ \delta_P \neq 0 $) from Eq. (11) as listed in Eqs. (15). The vertical grey line and associated band is the uncertainty-weighted average.

The uncertainty-weighted mean of these values is $ \delta_P= 0.218 \pm 0.193 $. At low-$ Q^2 $, $[1-\mu_p G_E(Q^2)/G_M(Q^2)] \propto\delta_P$; so, improved accuracy might be attained through new high-precision low-$ Q^2 $ polarisation-transfer measurements of $ \mu_p G_E(Q^2)/G_M(Q^2) $ [85, 89].

Regarding Fig. 4, it is clear that analyses of existing data are consistent with Eq. (11) at a 26% confidence level, i.e., with this probability, the proton's Dirac and Pauli radii are not truly independent observables; so, Eq. (11) is established as a useful approximation. Moreover, the analyses reveal that the proton's Pauli radius is very probably larger than its Dirac radius. These conclusions place valuable constraints on pictures of proton structure [85]. This is important because there are many possible contributions to $ F_2(Q^2\simeq 0) $, including [11, 17, 24, 25, 90–96]: quark anomalous magnetic moments, orbital angular momentum within the proton, meson cloud contributions, etc.; and there are related influences on $ F_1(Q^2\simeq 0) $. An insightful understanding of the accuracy of Eq. (11) will require an explanation of how these factors combine to yield a simple (approximate) algebraic outcome; and, tied to that, experiments and analyses that supply a precise result for $ r_p^{{\rm{mag}}} $.

Comparing Fig. 5 with Fig. 1, it is plain that mesons and baryons are distinct forms of hadron matter. Indeed, the pion (π) is Nature's lightest hadron [4]; and with a mass similar to that of the μ-lepton, the π can seem unnaturally light. Modern theory explains the pion as simultaneously both a pseudoscalar light-quark+antiquark bound-state and, in the absence of Higgs-boson couplings into QCD, the massless Nambu-Goldstone boson which emerges as a consequence of dynamical chiral symmetry breaking. This dichotomy is readily explained so long as one employs a symmetry-preserving treatment of QCD's quantum field equations [3, 27, 97–101]. Such proofs reveal that pion observables provide a very clean window onto EHM [9, 10, 12–14, 19].

Figure 5.  (color online) In contrast to the proton in Fig. 1, the $ \pi^+ $ contains one valence u-quark, one valence $ \bar d $-quark; yet, owing to the strong-interaction, it also holds infinitely many gluons and sea quarks. ($ \pi^- $ is $ d\bar u $ and $ \pi^0 $ is $ u \bar u - d\bar d $.)

Although the π was predicted and discovered more than seventy years ago [102, 103], very little has empirically been revealed about its structure [104–106]. This is primarily because pions are unstable. For instance, $ \pi^\pm $ decay via weak interactions characterised by a leptonic decay constant $ f_\pi = 0.0921(8) $GeV with a lifetime of just $ 26.033(5) $ns [4]. Herein, our focus is the electric radius of the charged pion.

With results in hand for the charged-pion elastic form factor, $ F_\pi(Q^2) $, then Eq. (1) yields the pion radius, $ r_\pi $. Many model and QCD theory calculations of $ r_\pi $, $ f_\pi $ are available. They show, inter alia, that on a large domain of meson mass, the product of radius, $ r_{0^-} $, and decay constant, $ f_{0^-} $, for ground-state pseudoscalar mesons is practically constant [107, Fig. 2A]:

This feature is a manifestation of EHM, expressed in pseudoscalar meson Bethe-Salpeter wave functions.

The pion form factor has been measured in pion+electron $ (\pi e) $ elastic scattering at low-$ Q^2 $ [108–111]. Extractions of $ r_\pi $ therefrom are broadly compatible; but the uncertainty is large, e.g., $ m_\pi $ is known with more than $ 10 000 $-times better precision. Given the pion's fundamental role in binding nuclei and its direct links with EHM, this is unsatisfactory; especially since $ m_\pi $ and $ r_\pi $ are also correlated [10, Fig. 2], entailing that precise values of both are necessary to set meaningful benchmarks for theory.

Today, one typically sees [4]:

a value obtained from the analysis of data in Refs. [111–113], supplemented by information from $e^+ e^- \to \pi^+ \pi^-$ reactions [112, 113]. This value has dropped $ 1.5 \sigma $ from that quoted two years earlier [114]: $ r_\pi = 0.672(8) $fm. So, it is worth reconsidering available $ F_\pi(Q^2) $ data, using the SPM to obtain a fresh $ r_\pi $ determination.

The densest, most extensive sets of low-$ Q^2 $ data were collected by the NA7 Collaboration [108, 109]. They are depicted in Fig. 6. Within mutual uncertainties, the $ r_\pi $ value published in Ref. [111] agrees with those reported in Refs. [108, 109]. However, referred to the results quoted in Ref. [109], the central value in Ref. [111] is over $ 2\sigma $ smaller.

Figure 6.  (color online) NA7 pion form factor data: upper panel — [108]; lower panel — [109]. Although collected thirty-five years ago, these remain the most dense and extensive sets of low-$ Q^2 $ data available. Concerning the 1986 set [109] — lower panel, the SPM analysis in [115] uses only the square-within-square data. The half-shaded data are too noisy to add information.

Reviewing available low-$ Q^2 ~\pi e \to \pi e $ data [115], only the NA7 sets were judged worth considering for reanalysis using the SPM. The $ Q^2 $-coverage, density, and precision of the data in Ref. [110] are insufficient to contribute anything beyond what is already available in the NA7 sets; and Ref. [111] reports only a radius obtained by assuming a monopole form factor in a fit to unpublished data.

As usual, the SPM analysis follows the pattern established by the $ r_p $ study, Secs. III, IV, using [115] the entire body of $ N=30 $ 1984 NA7 data on $0.015\leq Q^2/ {{\rm{GeV}}}^2\leq 0.119$ [108] and the $ N=35 $ lowest-$ Q^2 $ points from the 1986 NA7 set [109], $ 0.015\leq Q^2/{{\rm{GeV}}}^2\leq 0.144 $. Here, $ 6 \leq M \leq N/2 $, delivering $ {\cal{O}}(10^5 - 10^7) $ or $ {\cal{O}}(10^6 - 10^9) $ possible interpolators for the 1984 and 1986 data, respectively. Steps (i)–(iv) in Table 1 are implemented within these spaces [115]; and as with $ r_p $, $ r_p^{{\rm{mag}}} $, the optimal smoothing parameter $ \epsilon \simeq 0 $ in all $ r_\pi $-related cases.

Subjecting both sets of NA7 data to an array of robustness tests, like those detailed in Sec. VII below, it was found [115] that with $M \in {\cal{O}}_M=\{ M_{ji} = 2 + 4 j + i \vert j= 1,2,3, i=1,2,3\}$, one obtains neither a Gaussian distribution nor recovers the input radius. Given experience with $ r_p $ and $ r_p^{{\rm{mag}}} $, the appearance of such exceptional cases can be linked to the character of the NA7 data. On the other hand, all robustness tests were satisfied for $M\in {\cal{S}}_M= \{M_j=2+4j \vert\ j = 1,2,3,4\}$. Hence, these values were used to extract $ r_\pi $:

which are mutually consistent. Notably, the SPM result is $ \approx 1 \sigma $ below the error-weighted average of each of the original function-choice dependent determinations: $ \bar r_\pi^{84}= $ 0.650(8) fm and $ \bar r_\pi^{86}=0.659(5) $ fm. The uncertainty-weighted average of the radii in Eqs. (18) is

The results in Eqs. (18), (19) are compared with other extractions in Fig. 7: the SPM reanalysis supports the recent downward shift of the PDG average.

Figure 7.  (color online) Pion charge radius: purple circle (PDG, 2020) [4]; blue square (PDG 2018) [114]; cyan diamond [111]; green up-triangle [109]; light-green down-triangle [110]. The vertical grey band marks the uncertainty in Eq. (17). SPM results [115] — gold stars: SPM A-84 from [108]; SPM A-86 from [109]; and orange star — SPM combined, the average in Eq. (19).

Since extant $ \pi e $ elastic scattering data do not match analogous modern $ ep $ data in precision, density of coverage, or low-$ Q^2 $ reach, it is today sensible to refrain from judging the true size of $ r_\pi $. A final determination must await improved $ \pi e $ data.

Suppressing Higgs couplings into QCD, the kaon and pion would be identical; hence, differences between the K and π express Higgs-boson modulations of EHM, viz. interference between Nature's two known sources of mass. One of the simplest examples is the ratio of kaon and pion charge radii. However, the charged-kaon radius, $ r_K $, is very poorly constrained. The commonly quoted value is obtained [4] by averaging the results of separate monopole-squared fits to the distinct data sets [116, 117] depicted in Fig. 8. Employing the SPM to analyse these same data, one finds [115] that the largest possible roughness penalty is returned, $ \epsilon = 1 $. This means that the data contain insufficient information to yield an objective result for the radius: the value extracted is strongly influenced by the practitioner's choice of fitting form. Hence, $ r_K $ should be considered known with a precision that is not better than 10% [115].

Figure 8.  (color online) Existing low-$ Q^2 $ kaon form factor data (red squares [116] and blue circles [117]). Simple monopole-squared fits to this data are used for existing inferences of the charged-kaon electric radius [4]. The shaded region indicates the anticipated expansion of $ Q^2 $ coverage achievable with high-luminosity K beams at a new QCD facility [7, 8].

A new QCD facility at CERN promises to deliver high-luminosity π and K beams [7, 8]. Thus, one may expect the next decade to deliver the first new $ \pi e $ and $ K e $ elastic scattering data in over thirty years, with the anticipated $ Q^2 $ coverage indicated by the shaded area in Fig. 8. Should such data on this expanded domain meet the stringent demands of density and precision that have above been established as necessary, then they will enable the first truly objective determinations of pion and kaon charge radii.

The deuteron (d) is the simplest compound nucleus in Nature: as sketched in Fig. 9, a neutron and proton, somehow bound by SM strong interactions. Detailed nuclear theory studies of the deuteron exist. As illustrated elsewhere [118], they typically treat the neutron and proton as pointlike and bind them by some model nucleon+nucleon potential. No realistic QCD-connected calculations of deuteron structure are available.

Figure 9.  (color online) In terms of QCD's Lagrangian quanta, this image might be a fair sketch of the deuteron; but since no realistic QCD-connected calculations of deuteron structure are available, it is currently impossible to know with any confidence.

Notwithstanding the absence of SM theory, the deuteron mass is known with great precision [45]: the relative uncertainty is $ 3.0 \times 10^{-10} $. On the other hand, repeating the stories told above, $ r_d $, the deuteron radius, is far less well known. Extracted from $ ed $ elastic scattering, the value usually quoted is $ r_d = 2.111(19) $fm [45]. However, the value $ r_d = 2.12562(78) $fm is reported from measurements of $ 2P \to 2S $ transitions in $ \mu D $ atoms [75]. This is 25-times more precise than the $ ed $ scattering value and, referred to its own error, $ 19\sigma $ larger. Further clouding the issue, an analogous measurement using standard atomic deuterium returns $ r_d = 2.1415(45) $fm [119], which is $ 3.5\sigma $ above the $ \mu D $ value.

Within the context established by this level of disagreement, a new experiment has been proposed at Jefferson Lab [73, DRad]. Based upon the PRad experimental configuration, DRad aims for an extraction of $ r_d $ from $ e d $ elastic scattering with precision better than 0.25%. Supporting that effort and following the PRad approach, Ref. [120] reports a study that uses four parametrisations of the deuteron charge form factor, $ G_C^d $, with the aim of identifying fitting functions that can robustly extract $ r_d $ from precise data. As highlighted by the examples drawn above and hereafter shown explicitly, the statistical SPM provides a powerful alternative.

Our illustration proceeds by following the SPM validation approach employed in Refs. [47, 85, 115]. Namely, suppose experimental data on $ G_C^d $ can be obtained with similar precision and $ Q^2 $-coverage as PRad proton data. Then, choose thirteen different models for $ G_C^d $, each of which has its own associated value of $ r_d $. The selected models fall into two classes: (I) four parametrisations of available $ e d $ scattering data; and (II) nine theory-based models, ranging from elementary to elaborate. All correspond to a deuteron with unit electric charge. Class (I) is the collection of parametrising functions used in Ref. [120]: [121, Abbott1]; [122, Abbott2]; [123, Parker]; and [124, 125, SOG]. The simple models in Class (II) are a monopole, a dipole, and a Gaussian, each defined by $ r_d=2.1 $fm; and a quadratic and a cubic polynomial, with their own distinct radii, determined in fits to a model treatment of neutron+proton scattering [126]. The elaborate models in Class (II) are drawn from the analyses in Refs. [127, 128] and correspond to the following treatments of $ e d $ scattering: relativistic impulse approximation based on a one-boson-exchange nucleon+nucleon potential (IA); IA augmented by meson exchange currents (IAMEC); relativistic impulse approximation using a Reid soft core potential (RSC); and RSC augmented by meson exchange currents (RSCMEC). The associated form factors are determined via spline interpolations of the corresponding numerical results, which are available elsewhere [129] and are linked with the following radii (in fm): IA, $ r_{{d}}=2.15798 $; IAMEC, $ r_{{d}}=2.15727 $; RSC, $ r_{{d}}=2.0997 $; and RSCMEC, $ r_{{d}}=2.09912 $.

It is worth noting that all models in Class (II) lead to large $ \chi^2 $ values if used to fit existing experimental data.

As described in Ref. [120], DRad will emulate PRad kinematics, collecting data on $ 2 \times 10^{-4}\leq Q^{2}/{{\rm{GeV}}}^2 \leq 0.05 $ using electron beam energies $ 1.1 $, $ 2.2 $ GeV, with anticipated bin-by-bin statistical uncertainties in the range $ 0.02\%-0.07\% $ and systematic uncertainties of 0.06%– 0.16%. Working with this information, the next phase of the illustration consists in following the steps listed in Table 1.

Analysing the results from this procedure enables one to make the following observations.

[A] Defining $ {\cal{S}}_M^{1.1}= \{M_j=5+4 j | j=1,2,3\} $, $ {\cal{S}}_M^{2.2}= \{M_j= 2+4 j+i | j=1,2,3, i=1,2\} $, $ {\cal{S}}_M^{1.1 \cup 2.2}= \{M_j= 1+4 j+ 2 i | j=1,2,3, i=1,2\} $, and considering the different planned DRad electron beam energy settings, $ E_\gamma $, then for a given value of $ M^{E_\gamma}\in {\cal{S}}_M^{E_\gamma} $ and all thirteen fitting models, the distribution of SPM-extracted radii is a very clear Gaussian, which is centred on the radius input value and whose characteristics are practically independent of M: in these cases, $ \sigma_{{\delta M}} < \sigma^{{M_j}}_r $.

[B] Defining the bias: $ \delta r_M=r_M-r^{{\rm{input}}}_M $, then an assessment of the robustness of the SPM $ r_d $ extraction is enabled by Fig. 10. With beam-energy $ 1.1 $ GeV (left panel), the SPM delivers a robust radius result, marked by $ |\delta r_{{d}}|<\sigma_r $, for all data generators except the Gaussian and quadratic models. Similarly, for the combined $ 1.1 $ and $ 2.2 $ GeV data (right panel), the radius determined via the SPM is robust for all data generators except the Gaussian, quadratic and cubic models. On the other hand, with $ 2.2 $ GeV beam energy kinematics (central panel, different scale for $ \delta r_d $), one finds large biases and errors associated with all five elementary generators; and although the SPM leads to a robust radius determination for the generators based on experimental data, the exposed contrast suggests that analysis of the $ 2.2 $ GeV kinematic setup alone should be avoided.

Figure 10.  (color online) Bias, $ \delta r_{{d}} $, and associated standard error, $ \sigma_r $, for the SPM extraction of $ r_d $ from $ 1 000 $ data replicas generated using the thirteen models discussed in Sec. VII on the kinematic domain and with uncertainties anticipated in [120, DRad].

[C] Fig. 11 expands on [B] by quantifying the reliability of SPM analyses of anticipated DRad data using three measures: bias, $ \delta r_{d} $; standard deviation, $ \sigma_r $; and root mean square error (RMSE)

Figure 11.  (color online) Bias $ \delta r_{{d}} $, standard error $ \sigma_r $ and RMSE of the SPM extrapolation of the deuteron radius for $ 1 000 $ replicas generated from the thirteen models described in Sec. VII. 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.

Evidently, 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. Consequently, the SPM analyses satisfy a standard "goodness of fit" criterion [64]; hence, may objectively be judged to return a reliable expression of the information actually contained 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.

It is now worth comparing these remarks and Figs. 10, 11 with Ref. [120, DRad – Figs. 2, 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 $ r_p $ are unsuitable for the $ r_d $ case. So, a new method, tuned specifically to $ e d $ scattering, needed to be developed. It yielded two new best fitter functions, neither of which coincides with any of the forms used in connection with the proton radius extraction.

This contrasts markedly with the SPM validation tests. As shown above, for practical purposes, all procedures used in connection with the reliable extraction of $ r_p $ are equally effective for a SPM extraction of $ r_d $: 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. Moreover, the SPM analysis supplies a mathematical demonstration that, used alone, the $ 2.2 $GeV beam data cannot deliver an objective result for $ r_d $. Finally, considering the middle row of Fig. 11, the SPM delivers a mean variance $ \bar\sigma_r = $ 0.0055, 0.0075 for the $ 1.1 $ GeV and combined $ 1.1 $+$ 2.2 $ GeV beam data, respectively, which compares favourably with the mean variance ($ 0.0035 $) obtained using the $ ed $-tuned best fitters in Ref. [120, Fig. 6]. In applying the SPM to anticipated DRad data and subject to the considerations discussed in this section, one would be using 75-million (45-million with $ E_\gamma=2.2 $ GeV excluded) independent deuteron radius values in the determination of $ r_d $ via Eq. (4).

As a solution is sought for the strong interaction problem within the Standard Model (SM), precise knowledge of the Poincaré-invariant radii of hadrons (e.g., proton, pion, kaon), and light nuclei (e.g., deuteron) has become a high priority. Yet, recent measurements and new analyses of older data have revealed uncertainties and imprecisions that severely limit the use of these radii as benchmarks for SM theory. The last decade has revealed that reliable extraction of such radii from electron+hadron scattering requires: (i) precise data, stretching as close as possible to $ Q^2=0 $, and densely covering a domain $ Q^2\lesssim 0.1 $GeV$ ^2 $; and (ii) elimination of all subjectivity/bias associated with practitioner-dependent choices made in fitting the data and subsequently extrapolating to estimate a $ Q^2=0 $ derivative. Herein, we have not addressed requirement (i), which is a challenge for facility hardware, and experiment design and conduct. Regarding requirement (ii), different approaches have been suggested, each typically tuned to a given reaction. We have reviewed an alternative that is applicable in the same form to all cases; namely, the statistical Schlessinger point method (SPM).

The SPM is a tool for interpolating data (broadly defined); and, significantly, from the interpolations, delivering an estimate, with known uncertainty, for the curve underlying that data. This information may then be used to predict the values of quantities properly determined by data outside the domain of existing measurements. The SPM's basic strengths are a rigorous mathematical foundation in analytic function theory; the absence of assumptions about theories which may underly the data; the elimination of practitioner-induced bias; and mechanical, uniform implementation in diverse, unrelated contexts, viz. the absence of a need for system-specific tuning. With these foundations, the SPM returns an objective expression of the information contained in the data under consideration; and because it makes no reference to theories of physical phenomena, that expression serves as a test for both the experiment and potentially relevant theory.

Whilst it will now be clear, it is nevertheless worth emphasising that the SPM can be used for much more than extracting radii. The range of existing applications includes the extrapolation of theory predictions for hadron form factors to large values of $ Q^2 $, e.g., Refs. [35, 94, 96]; continuation from Euclidean to Minkowski space in order to obtain light-front specific quantities from Poincaré-covariant wave functions, e.g., Refs. [130, 131]; 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. [132]; and the analysis of light-nuclei deep inelastic scattering data [133] and its interpretation in terms of parton distribution functions [134]. There is no limit to the scope of the SPM. Given precise data, densely covering an appropriate domain, it 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.

In preparing this article we benefited from constructive interactions with L. Chang, O. Denisov, J. Friedrich, H. Gao, D. W. Higinbotham, R. J. Holt, V. Mokeev, W.-D. Nowak, C. Quintans, P. E. Reimer, W. Xiong and J. Zhou. Use of the computer clusters at the Institute for Nonperturbative Physics at Nanjing University is gratefully acknowledged.