Zemach and Friar radii of the proton and neutron from lattice QCD

We present the first lattice-QCD result for the Zemach and Friar radii of the proton and neutron. Our calculation includes both quark-connected and -disconnected diagrams and assesses all sources of systematic uncertainties arising from excited-state contributions, finite-volume effects and the continuum extrapolation. At the physical point, we obtain for the proton $r_Z^p = ( 1.013 \pm 0.010\ (\mathrm{stat}) \pm 0.012\ (\mathrm{syst}) )~\mathrm{fm}$ and $r_F^p = ( 1.301 \pm 0.012\ (\mathrm{stat}) \pm 0.014\ (\mathrm{syst}) )~\mathrm{fm}$. These numbers suggest small values of the Zemach and Friar radii of the proton, but are compatible with most of the experimental studies.


I. INTRODUCTION
The most accurate determination of the proton's electric (charge) radius is derived from the measurement of the Lamb shift in muonic hydrogen spectroscopy [1,2].This result exhibits a large tension with some ep-scattering experiments [3,4], which is known as the "proton radius puzzle".
To infer the electric radius from the observed Lamb shift, higher-order nuclear structurecontributions need to be subtracted.The leading contribution is the two-photon exchange [5], the dominant, elastic part of which depends on the third Zemach moment of the proton [6][7][8], The associated radius is known as the Friar radius of the proton, A very large Friar radius was once suggested [9] as a possible solution to the proton radius puzzle.For this purpose, however, the Friar radius would need to be so large that the expansion in radii would break down [10,11].
While the traditional proton radius puzzle awaits its final resolution, the goal of reaching a consistent picture of all the fundamental electromagnetic properties of the nucleon has attained a new prominence.Historically, data-driven dispersive approaches had found values of the electric radii of the proton consistent with the lower value of muonic-atom spectroscopy measurements [12,13].For the magnetic properties, a tension between dispersive approaches [14] and z-expansion results [15] appeared, i.e., a separate puzzle beclouds the magnetic properties of the proton.Underlining the importance of the magnetic properties of the proton, several experiments are under way to measure these from spectroscopy on (muonic) hydrogen [16][17][18][19].This can be achieved by measuring, in addition to the Lamb shift, the hyperfine splitting (HFS) in either electronic or muonic hydrogen, which is caused by the magnetic spin-spin interaction between the nucleus and the orbiting lepton.The influence of the electromagnetic structure of the nucleus on the HFS is particularly pronounced for the S-states, since the S-state wavefunction has a large overlap with the nucleus.
The leading-order proton-structure contribution to the S-state HFS of hydrogen depends on the Zemach radius of the proton [7,20], Having a first-principles prediction of the Zemach radius prior to the experimental measurement of the ground-state (1S) HFS in muonic hydrogen with ppm precision [16][17][18][19], from which the Zemach radius could be extracted with sub-percent uncertainty, is highly desirable.
Beyond helping in narrowing down the frequency search range, such a prediction allows for a crucial consistency check.We note that the interpretation of the experimental HFS results relies on theoretical input for the proton-polarizability effect, where a discrepancy has emerged between data-driven approaches and baryon chiral perturbation theory [21].
Eventually, combining precise HFS measurements in electronic and muonic hydrogen, the proton polarizability can be determined from those as well and compared to theory [22].
In this letter we present the first lattice-QCD calculation of the Zemach and Friar radii, building on our results for the electromagnetic radii of the proton and neutron [23,24].Our results for the Zemach and Friar radii of the proton have a total precision of 1.5 %, and are well compatible with most of the experimental determinations [2,10,14,15,21,25].

II. LATTICE SETUP
In order to compute the Zemach and Friar radii of the proton and neutron, we need, according to eqs. ( 1) and ( 3), information on their electric and magnetic form factors.For our lattice determination of the latter, we employ a set of lattice ensembles with N f = 2 + 1 dynamical flavors of non-perturbatively O(a)-improved Wilson fermions [26,27], using the tree-level improved Lüscher-Weisz gluon action [28] and twisted-mass reweighting [29,30], which have been generated as part of the Coordinated Lattice Simulations (CLS) effort [31].The ensembles entering our analysis are listed in table I and cover four lattice spacings a ∈ [0.050, 0.086] fm as well as several pion masses down to slightly below the physical one (E250).The calculation of the reweighting factors and the correction of the strange-quark determinant are described in Refs.[32] and [33], respectively.We include the contributions from quark-connected as well as -disconnected diagrams.For further details concerning the setup of the simulations, the calculation of our raw lattice observables, the extraction of the form factors, and the treatment of excited states, we refer to Ref. [23].All dimensionful quantities are expressed in units of the gradient flow time t 0 [34].To this end, we use the numerical determination at the flavor-symmetric point, t sym 0 /a 2 , from Ref. [35].Only our final results for the radii are converted to physical units using the FLAG estimate [36] t 0,phys = 0.14464(87) fm (4) for N f = 2 + 1.

III. FITS TO BARYONIC χPT
In Refs.[23,24] we have combined the parametrization of the Q 2 -dependence of the form factors with the extrapolation to the physical point (M π = M π,phys , a = 0, L = ∞).For this purpose, we have fitted our form factor data to the next-to-leading-order expressions resulting from covariant baryon chiral perturbation theory (BχPT) [37].While explicit ∆ degrees of freedom are not considered in the fit, we include the contributions from the relevant vector mesons, as discussed in detail in Ref. [23].For the physical pion mass we use the value in the isospin limit [38], so that in units of t 0 , we employ √ t 0,phys M π,phys = 0.09881 (59).Here, the uncertainty of M π,iso in physical units is neglected since it is entirely subdominant compared to the uncertainty of the scale √ t 0,phys .
We perform several such fits, applying different cuts in the pion mass (M π ≤ 0.23 GeV and M π ≤ 0.27 GeV) and the momentum transfer (Q 2 ≤ 0.3, . . ., 0.6 GeV 2 ), and, at the same time, varying our model for the lattice-spacing and/or finite-volume dependence, in order to estimate the corresponding systematic uncertainties.The aforementioned relatively strict cuts in Q 2 are required because the BχPT expansion, from which our fit formulae are derived, is only applicable for low momentum transfers.By including the contributions from vector mesons, the range of validity of the resulting expressions can be extended [37,39,40].Nevertheless, as the heaviest vector meson we consider in the isovector channel is the ρ, momentum transfers larger than M 2 ρ ≈ 0.6 GeV 2 cannot safely be described in this way.For further technical details on our implementation of the BχPT fits, we refer to Ref. [23].
We have extensively crosschecked our excited-state analysis as well as the parametrization of the Q 2 -dependence and the extrapolation to the physical point; for details, see Ref. [23] and its appendices.

IV. EXTRAPOLATION OF THE FORM FACTORS AND INTEGRATION
Given that the Zemach radius and third Zemach moment are defined as integrals over all possible (spacelike) values of Q 2 [cf.eqs.( 1) and ( 3)], an extrapolation of the BχPT fits beyond their range of applicability is required if they are to be employed to parametrize the form factors.For each model, we evaluate the BχPT formula for G p,n E and G p,n M , using the low-energy constants as determined from the corresponding fit, at the physical point and at twenty evenly spaced points in Here, Q 2 cut is the cut in the momentum transfer corresponding to the respective variation of the BχPT fit.
In the next step, we fit a model which obeys the large-Q 2 constraints on the form factors from perturbation theory [41] to these data points and their error estimates.We note that the data points exhibit an extremely high correlation due to the way we generate them.Taking these correlations into account when adjusting the extrapolation model would thus not be meaningful, and also technically challenging because the resulting covariance matrices are extremely badly conditioned.To describe the Q 2 -dependence, we use the model-independent z-expansion [42], with where we employ τ 0 = 0 and τ cut = 4M 2 π,phys .We truncate the z-expansion beyond m = 9, and incorporate the four sum rules from Ref. [43] for each form factor, which ensure the correct asymptotic behavior of the latter for large Q 2 .The normalization of the electric form factor is enforced by fixing a p 0 = 1 and a n 0 = 0, respectively.For the determination of the Zemach radius, we fit G E and G M simultaneously, similar to the crosscheck of our analysis in Ref. [23], so that we have eleven independent fit parameters altogether.For the third Zemach moment, on the other hand, only the electric form factor is required, so that we fit only G E and have five independent fit parameters here.The extrapolation fits are performed for the proton and neutron independently.Using more than twenty data points for each form factor or a higher degree of the z-expansion does not increase the overlap between the original BχPT fit and the extrapolation any further.
For the numerical integration of eqs.( 1) and (3), we smoothly replace the BχPT parametrization of the form factors by the z-expansion-based extrapolation in a narrow window around Q 2 cut .Concretely, we use the following estimate for the form factor term, where for the third Zemach moment, respectively.In eq. ( 9), F χ (Q 2 ) represents our fit to BχPT, while F z (Q 2 ) denotes the z-expansion parametrization of the form factors.For the width of the window in which we switch between the two parametrizations, we choose ∆Q 2 w = 0.0537t −1 0 ≈ 0.1 GeV 2 .We remark that for a consistent calculation of the third Zemach moment, the replacement according to eq. ( 9) has to be applied to all terms in eq. ( 1), i.e., also to the value of ⟨r 2 E ⟩.The cancellation between the different terms of eq. ( 1) at small Q 2 does not occur at the required numerical accuracy on all our bootstrap samples.To facilitate the numerical integration, we therefore regulate the small-Q 2 contribution to the integral for the proton by replacing t 0 Q 2 → t 0 Q 2 + 1 × 10 −7 in the denominator, which changes the result for ⟨r 3 E ⟩ p (2) by less than 10 % of its statistical error.
The two parametrizations and their weighted average according to eq. ( 9) are illustrated in fig. 1 for the case of the Zemach radius of the proton.While the BχPT formula is clearly not reliable for Q 2 ≳ 1.7 GeV 2 ≈ 0.9t −1 0 , the z-expansion behaves well for arbitrarily large momenta, which is due to the sum rules [43] we have included.In the region where we adjust the z-expansion to the BχPT parametrization (0 < Q 2 ≤ 0.6 GeV 2 for the case shown in fig.1), however, the two curves overlap so closely that they are indistinguishable by eye.The blue curve, which is the one we use for the integration, smoothly switches from the orange (BχPT) curve to the green (z-expansion) one in a tight window around Q 2 cut = 0.6 GeV 2 = 0.322t  The orange curve shows one of the BχPT fits to our lattice data with Q 2 cut = 0.6 GeV 2 ≈ 0.322t −1 0 , the green curve the z-expansion-based extrapolation, and the blue curve the weighted average of the two according to eq. ( 9).
Replacing the BχPT parametrization smoothly with a constant zero instead of the zexpansion-based extrapolation [i.e., setting F z (Q 2 ) ≡ 0 in eq. ( 9)] allows one to estimate the contribution of the form factors at Q 2 > Q 2 cut to the resulting Zemach radius and third Zemach moment, respectively.For Q 2 cut = 0.6 GeV 2 (our largest, i.e., least stringent, value for the cut), we find that the relative difference of the thus obtained value for r p Z to the actual result using the corresponding variation of the BχPT fits is less than 0.9 %.In other words, the form factor term at Q 2 > 0.6 GeV 2 contributes less than 0.9 % to the Zemach radius of the proton.For the third Zemach moment, the denominator in the integrand suppresses the large-Q 2 contribution to the integral even more strongly than for the Zemach radius.Accordingly, we find a corresponding relative contribution of less than 0.3 % to the third Zemach moment of the proton.
Due to this smallness of the contribution of the extrapolated form factors, the precise form of the chosen model for the extrapolation only has a marginal influence on the resulting values for the Zemach radius and third Zemach moment.For example, if we replace the z-expansion by a dipole ansatz (which also fulfills the constraints from Ref. [41]), we find that the Zemach radius of the proton derived from any of our fit variations changes by at most 20 % of the entire systematic error quoted in eq. ( 10) below.Thus, adding the variation in r p Z due to the extrapolation model quadratically to the systematic uncertainty in eq. ( 10) would not change the latter significantly.
Finally, we note that the major advantage of our approach based on the BχPT fits over an integration of the form factors on each ensemble is that the Zemach and Friar radii can be computed directly at the physical point, so that an extrapolation of results for the radii to the physical point, which would entail further significant systematic uncertainties, is not required.

V. MODEL AVERAGE AND FINAL RESULT
As in Refs.[23,24], we do not have a strong a priori preference for one specific setup of the BχPT fits.Thus, we determine our final results as well as the statistical and systematic error estimates from an average over the different fit models and kinematic cuts, using weights derived from the Akaike Information Criterion (AIC) [44][45][46][47][48][49].All values for the Zemach radii and third Zemach moments entering the average are listed in the Supplemental Material, together with the associated weights.More details on our procedure can be found in section V of Ref. [23].As our final results, we obtain r n Z = (−0.0411± 0.0056 (stat) ± 0.0040 (syst)) fm, (12) This corresponds to Friar radii of r p F = (1.301± 0.012 (stat) ± 0.014 (syst)) fm and r n F = (0.198 ± 0.017 (stat) ± 0.010 (syst)) fm, respectively.
In fig.2, our numbers for the proton are compared to other determinations based on experimental data.There are three main types of experiments which have been employed in the literature to compute the Zemach radius of the proton: muonic hydrogen HFS [2], electronic hydrogen HFS [50], and ep scattering.In order to extract the proton Zemach radius from an HFS measurement, input on the proton-polarizability effect is required.This can be either taken from BχPT [21] or evaluated in a data-driven fashion, i.e., using information on the spin structure functions [51][52][53] (as was done in Refs.[2,25]).The form factors measured in ep-scattering experiments, on the other hand, can be analyzed with many different fit models, e.g., by employing a (modified) power series [10], a z-expansion [15], or dispersion theory [14].
This work disp., ep data [14] z-exp., ep data [15] A1 ep scatt.[10] H HFS [25] BχPT+H HFS [21] BχPT+µH HFS [21] µH HFS [2] 2.25 2.50 2.75  While our result for r p Z agrees within one combined standard deviation with the extractions based on BχPT [21] and the z-expansion-based analysis of world ep-scattering data [15], and still within two combined standard deviations with the data-driven HFS extractions [2,25] and the analysis of the A1 ep-scattering experiment [10], we observe a 2.6 σ tension with the dispersive analysis of world ep-scattering data [14].We also note that our estimate is smaller than all of the above experimental determinations except the one combining BχPT and electronic hydrogen HFS, which is slightly smaller than ours.
The proton's third Zemach moment can be extracted from ep-scattering experiments in the same way as its Zemach radius, and we also compare to these results in fig. 2. Again, our value is comparatively small, but this time in good agreement with both the z-expansion-based [15] and dispersive analyses [14].Against the analysis of the A1 ep-scattering experiment [10], on the other hand, we observe a clear tension of 5.3 σ in ⟨r 3 E ⟩ p (2) .In interpreting the aforementioned discrepancies, one must take into account that our results for the Zemach radii and third Zemach moments are not independent from those for the electromagnetic radii [23,24] because they are based on the same lattice data for the form factors and the same BχPT fits.Indeed, we observe a correlation of around 80 % both between ⟨r 2 E ⟩ p and r p Z and between ⟨r 2 M ⟩ p and r p Z , while our correlation between ⟨r 2 E ⟩ p and r p F is even larger, about 95 %.A large positive correlation between the proton's electric and Zemach radii has also been reported in the experimental literature [22,54].Hence, our small results for ⟨r 2 E ⟩ p and ⟨r 2 M ⟩ p in Refs.[23,24] naturally imply small values for r p Z and r p F .By contrast, the dispersive analysis [14] arrives at a significantly larger magnetic radius than the A1-data analyses [3,43] and our lattice-QCD-based extraction [23,24].This may explain why we observe a larger tension in the Zemach radius (which equally probes electric and magnetic properties) with Ref. [14] than with Ref. [10], even though the situation is exactly the opposite for the third Zemach moment / Friar radius (which only probes the electric properties).For a deeper understanding of the underlying differences, a comparison of the full Q 2 -dependence of the form factors would be required, rather than merely of the radii.Furthermore, the role of higher-order electromagnetic corrections should be clarified.
The Zemach radius of the proton can also be computed in the framework of heavy-baryon chiral perturbation theory [55], which yields a much larger value of r p Z = 1.35 fm.However, the authors of Ref. [55] do not quote an error estimate on this number and claim it to be in good agreement with the experimental results, so that the uncertainty is presumably rather large.
Our results for the neutron are very well compatible with the z-expansion-based analysis of world en-scattering data [15], albeit with a more than 2 times larger error.

VI. CONCLUSIONS
We have performed the first lattice-QCD calculation of the Zemach and Friar radii of the proton and neutron, which includes the contributions from quark-connected and -disconnected diagrams and presents a full error budget.The overall precision of our results for the proton is sufficient to make a meaningful comparison to data-driven evaluations.Our final estimates, which are given in eqs.(10) to (13), point to small values for the Zemach and Friar radii of the proton, but are consistent with most of the previous determinations within two standard deviations.We agree rather well with the dispersive analysis of Ref. [14] regarding the electric properties of the proton (i.e., the Friar radius), but to a much lesser degree on its magnetic properties (i.e., the Zemach radius).
We stress that our results are highly correlated with those for the electromagnetic radii [23,24].Thus, our relatively low values for the Zemach and Friar radii of the proton are not unexpected, and they do not give rise to an independent puzzle from the lattice perspective.
Table II

Figure 1 .
Figure1.Product of the electric and normalized magnetic form factors of the proton at the physical point evaluated with different parametrizations.The orange curve shows one of the BχPT fits to our lattice data with Q 2 cut = 0.6 GeV 2 ≈ 0.322t −1 0 , the green curve the z-expansion-based extrapolation, and the blue curve the weighted average of the two according to eq. (9).

Table I .
[23]view of the ensembles used in this study.For further details, see table I of Ref.[23].