Direct Comparison of using a $Z$-Transformation instead of the traditional $Q^{2}$ for Extraction of the Proton Radius from $e-p$ Scattering Data

The Proton Radius Puzzle and the recent increase in available elastic electron-proton scattering data has invigorated discussions into the best methodology for analyzing elastic scattering data to extract the proton electric form factor, $G_{E}^{p}$ [1]. The Proton Radius Puzzle refers to a $\sim 7\sigma$ disagreement in the proton rms charge radius, $\sqrt{\left<r_{p}^{2}\right>}$ (written as $r_{p}$ henceforth), when measured with a new muonic hydrogen spectroscopy technique versus conventional measurement techniques (atomic hydrogen spectroscopy and $e-p$ scattering).

Of particular note is the discussion surrounding the technique for extracting $r_{p}$ from $e-p$ scattering data, as it often requires one to make model-dependent analysis choices that affect the final result [2]. In the case where a theoretical model, such as in Ref. [3], is not used, the most impactful choice is that of the $G_{E}^{p}$ fit function and any constraints on that function.

As has been shown [4, 5] that the consistent definition of the proton’s charge radius for all types of radius measurements is

$$r_{p}\equiv\sqrt{-6\left.\frac{\partial G_{E}^{p}\left(Q^{2}\right)}{\partial Q^{2}}\right\rvert_{Q^{2}=0}}$$

(1)

and thus simply determined by the slope of $G_{E}^{p}$ at a $Q^{2}=0$. For electron scattering, which cannot be measured to $Q^{2}=0$, the implication of this is that $G_{E}^{p}$ must be fit with some functional form and then extrapolated to $Q^{2}=0$ to determine the slope as it approaches the limit. This has proven to be a challenge, as the true functional form of $G_{E}^{p}$ is not known. As such, the form chosen must be sufficiently flexible so as to contain the true value of $r_{p}$ within its parameter-space as well as robust enough to not diverge rapidly outside of the region of fitted data, particularly as $Q^{2}\rightarrow 0$.

One such method that has seen increased use is the so-called $Z$-transformation [6]. In this scheme, the negative four-momentum transfer, $Q^{2}$, range is conformally mapped to a unit circle. This mapping is done with the formalism:

$$Z\left(t,t_{\text{cut}},t_{0}\right)=\frac{\sqrt{t_{\text{cut}}-t}-\sqrt{t_{\text{cut}}-t_{0}}}{\sqrt{t_{\text{cut}}-t}+\sqrt{t_{\text{cut}}-t_{0}}}$$

(2)

where $t=q^{2}=-Q^{2}$, $t_{\text{cut}}$ is the highest mapped $Q^{2}$ value, and $t_{0}$ is a free parameter.

The use of the $Z$-transformation is well-motivated. This method is the standard tool for meson transition form factor studies [6, 7, 8, 9, 10, 11, 12]. By applying this mapping with a cut below the two-pion production threshold, the form factor is restricted to the region of analyticity (i.e. it can be represented by a convergent power series). This implies that with a sufficient number of parameters the uncertainty from truncating a power series fit can be minimized, while still arriving at a function with sufficient predictive power. In all, it is an enticing choice that has been used in a large number of analyses [13, 14, 15, 16, 17].

When fitting $Z$-transformed data, Eq. 1 must be adjusted to accommodate this coordinate transform:

$$r_{p}\equiv\sqrt{-6\left.\frac{\partial G_{E}^{p}\left(Z\right)}{\partial Z}\cdot\frac{\partial Z\left(-Q^{2},t_{\text{cut}},t_{0}\right)}{\partial Q^{2}}\right\rvert_{Q^{2}=0}}$$

(3)

However, this method is not without it’s drawbacks [18]. Namely, by nature of design, while the true form factor exists within this parameter-space, so do many incorrect form factor parameterizations that can create local $\chi^{2}$ minima in the fitting routine. One posited solution to this is the use of physically motivated constraints on the power series coefficients (Ref. [6] suggests using $\left|p_{i}\right|\leq 10$ as the coefficients for a $i$-th order power series). The use of these constraints, while conservative, introduces a model-dependence to the fit. Model-dependence, in and of itself, is not problematic; constraints placed on typical form factor polynomial fit (or any other commonly used functions) are model assumptions on the shape and behavior of the form factor. This is only to say that one must take care and be aware of the ramifications of model inputs.

Studies similar to this one have been performed previously, as in Refs. [19, 20], which we have used as a framework for our study. Particularly, we aim to address claims made in Ref. [21] that reanalyzes the PRad [22] proton electric form factor data using the $Z$-transformation technique. A claim is made that the PRad proton radius uncertainties are underestimated due to the use of the “Rational(1,1)” fit function instead of the $Z$-transformation technique. This claim, if justified, naturally leads to a conclusion that the results from the $Z$-transformation technique are a lower limit on precision and accuracy in the extraction of the proton radius when using realistic data. This work aims to test that claim.

Here, we aim study the power of this technique to accurately measure the proton charge radius relative to that of other techniques. We acknowledge that the Ref. [6] posit that trial fits of model data are problematic, as models have to make assumptions on the behavior of the form factor. However, we argue that by studying the technique with several models the affects of these assumptions can be mitigated and a better understanding of the technique can be realized.

In this study we closely follow the methodology of Ref. [19]. To test the robustness of the $Z$-transformation, we generate $G_{E}^{p}$ pseudodata from parameterizations with a known radius. These pseudodata are then analyzed with several fit functions both in $Q^{2}$ and transformed to $Z$. This technique allow us to assess the sensitivity of the $Z$ conformal mapping to the input and statistical fluctuations of data. For this study, we use 6 different parameterizations of $G_{E}^{p}$ as generators:

• •

  Alarcón, Higinbotham, Weiss, and Ye (AW) fit [3], a parameterization in terms of radius-independent and -dependent parts ($r_{p}=0.844$fm)

• •

  Arrington, Melnitchouk, and Tjon (AMT) fit [23], a rational(3,5) parameterization ($r_{p}=0.878$fm)

• •

  Arrington fit [24], an inverse-polynomial fit ($r_{p}=0.829$fm)

• •

  Bernauer fit [25, 26, 27], a 10th-order polynomial fit ($r_{p}=0.887$fm)

• •

  Standard dipole fit [28] ($r_{p}=0.811$fm)

• •

  Kelly fit [29], a rational(1,3) fit ($r_{p}=0.863$fm)

Each of these parameterizations are used to generate equally spaced (in $Q^{2}$) values of $G_{E}^{p}$ and apply an uncertainty of $0.2\%$ to each point, corresponding to a $0.4\%$ uncertainty on the measured cross section. These pseudodata begin at $Q^{2}_{\text{min}}=0.004$ GeV² and are spaced every $0.001$ GeV² up to a variable cutoff $Q^{2}_{\text{max}}$. $Q^{2}_{\text{max}}$ is varied by adding a single data point at a time up to a total of 500 points (that is, $0.005$ to $0.504$ GeV²). At each $Q^{2}_{\text{max}}$ value, 5000 data sets are generated for each parameterization where each data point is smeared by a Gaussian distribution with a standard deviation equal to the $0.2\%$ uncertainty on the point.

A keen eye may note that some of these parameterizations incorporate corrections for two-photon exchange (TPE), while others do not. We make no attempt to change this as, particularly at high $\epsilon$, the corrections for TPE are quite small often at the sub-percent level [30]. While TPE has been suggested as a contribution to the Proton Radius and Proton Form Factor Puzzles, there is no clear consensus as to whether it can account for these differences [31]. While the inclusion or exclusion of TPE will change the underlying form of the data, our work here assumes no knowledge of the underlying form and is testing different fit functions ability to extrapolate the slope to $Q^{2}=0$ without this knowledge. We do not need the fit function to explain any other physical quantities or to be a true representation of the underlying function of the form factor.

Each data set is fit with polynomials of increasing order, $n=1\ldots 4$, as defined by

$$f\left(Q^{2}\right)=p_{0}\times\left[1+\sum\limits_{i=1}^{n}p_{i}\left(Q^{2}\right)^{i}\right]$$

(4)

where $p_{0}$ is a normalization term. By examining Eqs. 1 and 4, the radius is then extracted as $r_{p}=\sqrt{-6p_{1}}$. We also explore applying bounds to the fitting parameters and the use of Rational(N,M)-type functions as used in the recent PRad analysis [22]. For each of the 5000 fits done with each polynomial at each $Q^{2}_{\text{max}}$ point, the extracted radius, $r_{p}$, and $\chi^{2}$ are recorded. The final radius from each polynomial at each $Q^{2}_{\text{max}}$ point is treated to be the radius extracted from the mean slope of all fits and the uncertainty of the slope is taken to be the standard deviation of the slopes fitted and then propagated to the mean extracted radius. Additionally, to assess the power of the fitting techniques, at each $Q^{2}_{\text{max}}$ point we calculate the Mean Squared Error as

$$\text{MSE}=\text{bias}^{2}+\sigma^{2},$$

(5)

where bias is the difference between the mean extracted radius and the input radius, and $\sigma$ is the rms-spread of extracted radii.

Each of these data sets is then transformed from $Q^{2}$ to $Z$ using Eq. 2. For this work, we set $t_{\text{cut}}=4m_{\pi}^{2}$, where $m_{\pi}$ is the charged pion mass, to restrict data below the two-pion production threshold and $t_{0}=0$. These values were chosen to match the values that were used in the proof-of-concept performed in Ref. [6]. The $Z$-transformed data is then fit with polynomials of increasing order, $n=1\ldots 4$ (analogous to Eq. 4 with $Q^{2}$ replaced by $Z$). This functional form, along with this choice of $t_{0}$, allows us to write Eq. 3 as $r_{p}=\sqrt{-1.5p_{1}/t_{\text{cut}}}$.

These $Q^{2}$ and $Z$-transformed data are then refit with a polynomial but with added bounds to constrain the fits. In each of these, prescribed bounds from literature are followed for the $Q^{2}$ [2] and $Z$ [6] data. As the $Z$-transformation dramatically changes the shape of the $G_{E}^{p}$ data so that bounds for data in $Q^{2}$ are not applicable to $Z$-transformed data and vice versa. This means that the bounds applied to fits in $Q^{2}$ and fits in $Z$ are not equivalent to each other. However, given that the techniques allow for different formulation of bounds on the fit parameters, the authors assert that the most fair technique to compare (un)bounded fits in $Q^{2}$ to (un)bounded fits in $Z$.

When assessing the quality of these fits in this study, there are two qualities to keep in mind:

1. 1.

   How well does the fit reflect the data?

2. 2.

   How well does the fit extract the input radius?

The first question is easily addressed with a simple ‘goodness of fit’ test. In this work, we report the mean and standard deviation of the reduced chi-squared, $\chi^{2}_{r}$, calculated for each fit at each $Q^{2}_{\text{max}}$. A $\chi^{2}_{r}$ of approximately 1 indicates that the fit represents the data well, whereas $\chi^{2}_{r}\gg 1$ indicates underfitting and $\chi^{2}_{r}\ll 1$ indicates overfitting. An important note is that while a satisfactory $\chi^{2}_{r}$ value is critical for extracting the proton radius, as the data must be properly fit, it does not indicate that the fit will accurately extrapolate to guarantee a correct radius extraction [32]. This statement is even more pertinent as additional terms are added to the fit. Additional terms bring in additional moments that can cause the fit function to vary wildly and unpredictably when extended beyond the range of the measured (or simulated in this case) data.

The second question is addressed by studying the bias and the variance of the extracted radii. The bias is defined as the difference of the mean radius extraction at each $Q^{2}_{\text{max}}$ and the input model radius. The variance is the standard deviation of the mean extracted radius. In this context, both of these quantities are best assessed when compared to the magnitude of the proton radius puzzle; an unreasonably large bias could yield a result that suggests the incorrect solution to the puzzle, whereas an unreasonably large variance would be unable to discern between the purported solutions. To aid in this assessment, lines representing the magnitude of the proton radius puzzle are added to the bias-variance plots in this work. As these two values can often be inversely correlated when optimizing a fit, we choose to define the error using the Mean Squared Error as in Eq. 5. A smaller MSE indicates an overall improvement in the radius extraction. It should be noted that just as having $\chi^{2}_{r}\approx 1$ does not indicate that the extracted radius will be correct, having a low MSE is also insufficient if the fit does not well represent the data. The two properties must be taken together to show that the fit can both reproduce the data and extract the correct radius.

All figures in this article show the fitting techniques applied to the AMT parameterization as an illustrative set. See the Supplemental Information for this work for the full set of figures for all parameterizations as well as textual descriptions of all figures.

The first test we perform on using this methodology is with unbounded polynomials. The authors note that while Ref. [6] proposes that bounds of $\left|p_{i}\right|\leq 10$ be used, these bounds, while conservative, are based on a vector-dominance model. In order to initially assess the robustness in the least model-dependent way, we opt to start by omitting any bounds.

Comparing the outputs of the unbounded polynomial fits in Figs. 1 and 2, we can see that there is a trade-off between choosing to fit in $Q^{2}$ or $Z$. The unbounded $Q^{2}$ fit has a larger mean bias, that decreases with increasing polynomial order. However, this fit has very little variance in the fit value once a moderate $Q^{2}$ cutoff is reached. The unbounded $Z$ fit, on the other hand, has a much smaller mean bias for most models checked (though for some models it diverges rapidly once past a moderate $Q^{2}$ cutoff). The trade-off is that the variance in results from the $Z$ fit is quite large.

Given these results, it is ill-advised to use either of these methods blindly. There may exist some data set where one of these is the optimal choice; however that data set does not exist within this exercise. When comparing the fits in $Q^{2}$ to $Z$ the choice is, in essence, between a highly repeatable incorrect answer or a potentially correct answer that is not repeatable.

The next test involves placing bounds on the fit parameters. As discussed in the previous section, Ref. [6] proposes that, for the fit of $Z$-transformed data, bounds of $\left|p_{i}\right|\leq 10$ be used. The bounds are determined by beginning with a model-independent dispersion relation analysis, but then applies a model-dependent vector-dominance ansatz. The calculations could prescribe a more rigid set of bounds than this, however Ref. [6] decided that it was most appropriate to suggest a conservative implementation to minimize the model dependence.

For the $Q^{2}$ fit we follow the prescription of Ref. [2] and require that the fit parameters alternate signs. This bounded regression forces the result to approximate a monotonic function (in this case, non-increasing). This technique reduces the risk that the fit parameters will have dramatic fluctuations when extrapolated beyond the fit region.

The results of these fits, seen in Figs. 3 and 4, fall largely in line with the results of the unbounded polynomial fits. The main effect is a reduction of the variance, most notably seen below $Q^{2}_{\text{max}}{}\approx 0.1$ GeV². In this region, the variance is still too large to effectively resolve the proton radius puzzle. At $Q^{2}_{\text{max}}$ above this, the difference is negligible. The fit in $Q^{2}$ still has a very large bias and the fit in $Z$ with the vector-dominance bounds has very large variance.

The final functional form used is a Rational(N,M) function. A Rational(N,M) function is defined as the form

$$f\left(Q^{2}\right)=p_{0}\frac{1+\sum\limits_{i=1}^{N}n_{i}Q^{2}}{1+\sum\limits_{j=1}^{M}m_{j}Q^{2}}$$

(6)

with $p_{0}$ as a normalization factor. As an example, the Rational(1,1) function then parameterizes the data as

$$f\left(Q^{2}\right)=p_{0}G_{E}^{p}\left(Q^{2}\right)=p_{0}\frac{1+n_{1}Q^{2}}{1+m_{1}Q^{2}}$$

(7)

where the extracted radius is $\sqrt{-6\left(n_{1}-m_{1}\right)}$. This functional form was used for the extraction of the proton radius from the PRad experiment [22]. The decision was made by studying various fitters on generated pseudodata (much like this study) and selecting the most robust fitter for the $Q^{2}$ range of the measurement [20].

In this study, we vary N and M within the range of $0-2$ while restricting them to be within $\pm 1$ of each other. The results of these fits are seen in Fig. 5. An important note here is that the Rational(1,0) function is simply a linear polynomial which yields identical behavior to the linear fit in Fig. 1. Increasing the orders N and M simultaneously increases the $Q^{2}_{\text{max}}$ at which the bias is small compared to the proton radius puzzle, the minimum $Q^{2}_{\text{max}}$ at which the variance is small compared to the proton radius puzzle, and the time taken for the fits to converge. The reason for the increase in the time to convergence is that Rational(N,M) functions cover a very large parameter-space, such that increasing the order introduces additional local minima that slow down the fit. It is noted that it was attempted to include the Rational(2,2) function here, but convergence of the illustrative fits took a prohibitively lengthy amount of time.

For the fits in $Q^{2}$, the functional form of a Rational(N,M) function (provided both N and M at least 1) is by far the most precise and has the least bias. Unsurprisingly, given the study done for the PRad experiment, the Rational(1,1) function has minimal bias around $Q^{2}_{\text{max}}{}\approx 0.07\text{GeV}^{2}$ for nearly all models (though perfect agreement with the PRad study should not be expected as this analysis uses a higher $Q^{2}_{\text{min}}$, different point density, and different point-to-point uncertainties).

This functional form, when applied to $Q^{2}$ data, shows by far the most robust fitting of the parameterizations considered. However, as with all of the fitting forms studied, it will not always be the best choice in all scenarios.

From these illustrative fits, there are a number of conclusions that can be drawn:

• •

  When comparing fits in $Q^{2}$ and $Z$, both unbound and bound, no technique is clearly superior for extracting the proton radius. The technique that yields lower MSE depends not only on the order of polynomial and the $Q^{2}_{\text{max}}$ range, but also on the $G_{E}^{p}$ parameterization fit.

• •

  In contrast, the Rational(N,M) fitting procedure, when N and M are both at least 1 and a minimum $Q^{2}_{\text{max}}$ is reached, consistently has an MSE value of nearly an order of magnitude lower than the other fits, with the exception of when fitting the Bernauer parameterization. It should also be noted that the Rational(N,M) fits behave very poorly for low $Q^{2}_{\text{max}}$ values, even giving unphysical results for higher (N,M) orders. A caveat to the smaller MSE values is that it is partially driven by the substantially smaller variance of the fits as there appear to be fewer local minima in the parameter-space.

• •

  Applying bounds to the fitting procedure decreases the variance of the extracted radius. This is expected behavior, as bounds limit the phase space of possible solutions. However, bounds appear to do very little to improve the extracted bias. While physically motivated bounds (e.g. the slope must be negative to ensure a positive charge radius) are a useful tool, additional bounds may induce a false belief in the accuracy of the extraction.

• •

  As the true functional form of $G_{E}^{p}$ is not known, no single fitting function will be ideal for every possible data set. As such, it is imperative that any analysis that seeks to extract the proton charge radius does their due diligence to determine the best fitting function for their data. Ref. [20] details a procedure for an analysis of the most robust fitter for a given data range. Highlighting the care that must be taken, this reference also found that Rational(N,M) functions were inadequate for particularly low $Q^{2}$ data.

• •

  The true form factor is guaranteed to lie within the phase space of a polynomial fit in $Z$-transformed data. This technique aims not only to accurately extract the proton radius, but to gain insight into the true functional form of the proton form factor. This additional insight comes with a trade-off of increased uncertainties on the extracted proton radius. As thoroughly described in Ref. [33], modeling to explain data and modeling to predict data are two different goals that are not always achieved with identical techniques. The $Z$-transformation technique aims to achieve a deeper explanation of the data than the other techniques investigated, which in turn can negatively impact its predictive power.

• •

  As discussed in Sec. 2, it must be stressed that $\chi^{2}$ or MSE alone cannot be used to assess the extracted radius from a fit. The $\chi^{2}$ of a fit is only related to the ability of the fit to reproduce the fitted data. The MSE of a fit is only related to the extracted and input radius, with no regard to data points that were fit. Both values must be taken together to assess the strength of a particular fit.

In our work, there was no single data set where applying a $Z$ conformal mapping to the $Q^{2}$ range improved the extraction of the proton charge radius. This result should not imply that this is not a useful technique. Rather, the $Z$ transformation should be considered a useful tool in the physicists toolbox after they have assessed the correct tool for the job such as Ref. [20] finds that a 2^nd order polynomial fit in $Z$ is a robust fit for the PRad data points.

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under contracts DE-AC05-06OR23177 and DE-AC02-05CH11231. This material is based upon work supported in part by NSF award PHY-1812421.