A model-independent determination of the nucleon charge radius from lattice QCD

Lattice QCD calculations of nucleon form factors are restricted to discrete values of the Euclidean four-momentum transfer. Therefore, the extraction of radii typically relies on parametrizing and fitting the lattice QCD data to obtain its slope close to zero momentum transfer. We investigate a new method, which allows to compute the nucleon radius directly from existing lattice QCD data, without assuming a functional form for the momentum dependence of the underlying form factor. The method is illustrated for the case of the isovector mean square charge radius of the nucleon $\langle r^2_\mathrm{isov} \rangle$ and the quark-connected contributions to $\langle r^2_p\rangle$ and $\langle r^2_n \rangle$ for the proton and neutron, respectively. Computations are performed using a single gauge ensemble with $N_f=2+1+1$ maximally twisted mass clover-improved fermions at physical quark mass and a lattice spacing of $a=0.08\mathrm{fm}$.

Lattice QCD calculations of nucleon form factors are restricted to discrete values of the Euclidean four-momentum transfer. Therefore, the extraction of radii typically relies on parametrizing and fitting the lattice QCD data to obtain its slope close to zero momentum transfer. We investigate a new method, which allows to compute the nucleon radius directly from existing lattice QCD data, without assuming a functional form for the momentum dependence of the underlying form factor. The method is illustrated for the case of the isovector mean square charge radius of the nucleon r 2 isov and the quark-connected contributions to r 2 p and r 2 n for the proton and neutron, respectively. Computations are performed using a single gauge ensemble with N f = 2 + 1 + 1 maximally twisted mass clover-improved fermions at physical quark mass and a lattice spacing of a = 0.08 fm.

I. INTRODUCTION
The radius of the proton is a fundamental quantity for atomic, nuclear and particle physics. In atomic physics it enters in the determination of the Rydberg constant, the most precisely known constant in nature, as well as in precision tests of Quantum Electrodynamics. In nuclear physics it characterizes the size of the most abundant hadron in nature and in particle physics it is an input for beyond the standard model physics, testing lepton universality and the possible existence of new particles. Electron-proton scattering was traditionally used to determine the proton radius that is extracted from the slope of the electric Sachs form factor G E (q 2 ) extrapolated to zero momentum transfer. The proton radius value extracted from precision measurements of elastic electron scattering cross sections at the Mainz Microtron MAMI found a proton radius of r p = 0.879(5) stat (4) sys fm [1], where the systematic errors are added quadratically. Combining electron scattering data with more accurate data from the hyperfine structure of electronic hydrogen, the value recommended by CODATA was r p = 0.8775(51) fm [2]. In a pioneering experiment in 2010, the radius of the proton was extracted from the Lamb shift measured for muonic hydrogen [3]. The muonic hydrogen determination is much more accurate and the value extracted r p = 0.844087 (39) fm [4] differs by 7 standard deviations from the value recommended by CODATA. This gave rise to the proton radius puzzle. A followup, very accurate spectroscopic measurement using electronic hydrogen found a value that was inconsistent with previous determinations and consistent with the size extracted from the muonic hydrogen experiments, namely, r p = 0.8335(95)fm [5]. However, another spectroscopic measurement found a value of r p = 0.877 (13) fm [6] in agreement with the CODATA value. This discrepancy between the two most recent spectroscopic measurements conducted on electronic hydrogen remains unresolved and is under investigation. However, recently, the proton charge radius experiment at Jefferson Laboratory (PRad), overcoming several limitations of previous electron-proton scattering experiments, enabled measurements at very small forwardscattering angles. Using their results, they extracted a proton radius of r p = 0.831(7) stat (12) sys fm [7] in agreement with the muonic hydrogen extraction, towards a possible resolution of the puzzle. A theoretical determination of the proton radius starting from the fundamental theory of the strong interaction requires a non-perturbative framework. In this work we used the lattice QCD formulation to determine the radius of the proton. The standard procedure in lattice QCD is to compute the electric form factor for finite values of the momentum transfer and then perform a fit to determine the slope at zero momentum transfer. However, on a finite lattice the smallest non-zero momentum is 2π/L where L is the spatial size of the lattice. Therefore, to reach very small momentum transfers one needs very large lattices. In this work we use a direct method to extract the proton radius that does not depend on fitting the form factor. We illustrate the validity of our method by analyzing one ensemble of gauge configuration generated at the physical pion mass. The structure of is this paper is as follows: In section 2 we describe the general lattice setup, in section 3 we explain our new direct approach to extract the radius, in section 4 we discuss our results and in section 5 we compare to other lattice QCD determinations and give our conclusions.

II. LATTICE SETUP
Our lattice calculations are performed on a single ensemble with N f = 2 + 1 + 1 flavors of twisted mass [8] Wilson clover-improved fermions at maximal twist [9,10]. The quark masses on this ensemble are tuned to their physical masses, however, the simulations do not account for any isospin splitting for the light quarks. Furthermore, the ensemble has been tuned to maximal twist to implement automatic O(a) improvement [11,12]. The gauge ensemble is simulated with a spatial volume of (L/a) 3 = 64 3 and a temporal extent of T /a = 128. The lattice spacing and pion mass have been determined in Ref. [9] with values of a = 0.0809(2) stat (4) sys fm and M π = 138.0(4) stat (6) sys MeV, respectively. The value of M π is very close to its physical value in the isospin limit [13]. For further details on the simulation we refer to Ref. [9]. Measurements are performed on 750 gauge configurations, measuring on only every fourth configuration which corresponds to a separation of 4 HMC trajectories between consequent measurements. All statistical errors in our analysis are computed from a binned jackknife procedure taking into account correlations and possible residual affects of autocorrelations.

A. Electric form factor of the nucleon
The aim of this study is an extraction of the nucleon charge radius, hence we consider the electromagnetic matrix element of the nucleon where on the left-hand side N (p i , s i ) (N (p f , s f )) label nucleon states with initial (final) state momentum p i (p f ) and spin s i (s f ), and j µ is a vector current insertion which will be discussed below. On the right-hand side u(p i , s i ),ū(p f , s f ) denote Dirac spinors, m N the nucleon mass and we have introduced the Dirac and Pauli form factors F 1 (Q 2 ) and F 2 (Q 2 ), which depend on the Euclidean four-momentum transfer Q 2 = −q 2 with q = p f − p i . We will work in Euclidean spacetime throughout this study and in the following discussion we will always replace F 1 (Q 2 ) and F 2 (Q 2 ) by the electromagnetic Sachs form factors G E (Q 2 ) and G M (Q 2 ) which are more convenient for our purposes and related to F 1 (Q 2 ) and F 2 (Q 2 ) via The electromagnetic vector current j em µ (x) can be expressed in terms of vector currents j f µ (x) for the individual quark flavors f through where e f denotes the electrical charge for a quark of flavor f . On the lattice we employ the symmetrized conserved vector current which for a quark field χ f (x) in the twisted mass basis is given by where U µ (x) denotes a gauge link andμ a unit vector in the µ-direction. We restrict our calculation to light valence quarks and the relation of the light quark doublet χ(x) = (χ u (x), χ d (x)) T to the physical up and down quark at maximal twist is given by the chiral rotation In the SU (2) isospin symmetry limit, the matrix elements of the electromagnetic current satisfy where the isovector current j u−d µ (x) has been introduced and |p , |n denote proton and neutron states, respectively. Unlike the electromagnetic combination, the isovector current does not give rise to quark-disconnected diagrams in our lattice simulations. However, in this work we neglect any quark-disconnected contributions as our main goal is to present and benchmark a new method to extract nucleon radii. We remark that these contributions have been shown to be small compared to the quark-connected contribution in Refs. [14,15]. In the twisted mass lattice regularization the flavor symmetry is broken from SU (2) to the subgroup U (1) 3 and the isospin symmetry-based relations Eq. (8) are valid at non-zero lattice spacing up to lattice artifacts. In particular, the quark-disconnected contribution from the insertion of the isovector current j u−d will be discarded from our calculation, since it vanishes in the continuum limit. The point-split vector current in Eq. (5) is the Noether current associated with the residual flavor symmetry group U (1) × U (1) 3 and thus remains conserved. Therefore no additional multiplicative renormalization of the matrix element is required.

B. Two-and three-point functions
The lattice QCD evaluation of the nucleon matrix element in Eq. (1) requires the computation of two-and three-point functions where t i , t op and t f label initial (source), operator insertion and final (sink) timeslices, respectively. The spin projectors Γ ν are given by Γ 0 = 1 2 (1 + γ 0 ) for ν = 0 and Γ k = Γ 0 iγ 5 γ k for ν = k = 1, 2, 3. By time translation invariance only time differences are of physical relevance and it is convenient to introduce the source-sink time separation t sep = t f − t i and the shorthand t = t op − t i . Our kinematic setup is chosen such that the final state is produced at rest, i.e. p f = 0, p i = − q. Finally, the proton interpolating field is given in the physical basis by where C denotes the charge conjugation matrix and is transformed into the twisted basis in terms of χ u/d by using the chiral rotation in Eq. (7). Since nucleon structure calculations are known to be hampered by a severe signal-to-noise problem it is crucial to increase the overlap of the interpolating operator with the desired nucleon ground state, effectively suppressing excited states and allowing to extract a signal at smaller euclidean time separations. To this end we apply Gaussian smearing to the quark fields [16,17] with a choice of α G = 0.2 and N G = 125 smearing steps. The hopping matrix H ab ( x, y; U (t)) is given by where a, b are color indices. Furthermore, we use APE-smeared [18] gauge links U in the construction of H with a smearing parameter of α APE = 0.5 and N APE = 50 smearing steps. For the computation of three-point functions we use sequential inversions through the sink [19]. Therefore, we need to perform separate inversions for each choice of the source-sink time separation, the projector index and -in principle -the momentum at sink. However, the latter is fixed to be zero, as mentioned before. Since we are interested in the nucleon charge radius it is sufficient to consider three-point functions projected with Γ 0 . Three-point functions are computed for five values of t sep as listed in Table I. In order to obtain comparable effective statistics at each value of t sep we add additional source positions for increasing values of t sep . The number of source positions N s per source-sink time separation have also been included in Table I together with the total number of measurements N meas . The source positions are randomly and independently chosen on each gauge configuration. Two-point functions are computed with matching statistics for each value of t sep from the forward-propagators obtained in the calculation of three-point functions, hence there is a possible choice of either using matching statistics for two-and three-point functions at each value of t sep or always using the full available statistics for the two-point functions. However, we found that fully preserving correlation by using exactly matching statistics yields a slight advantage with respect to the resulting statistical fluctuations.
Calculations are performed using an appropriately tuned multigrid algorithm [20][21][22] for the efficient inversion of the Dirac operator that is required for the computation of the quark connected diagrams.

C. Ratio method
Extracting the physical matrix elements in Eq. 1 requires the cancellation of unknown overlap factors in the three-point function, which can be achieved by forming an optimized ratio involving two-point functions   [23-25] At large Euclidean time separations t and t f − t the ground state contribution is expected to dominate asymptotically and the ratio approaches a plateau The electromagnetic Sachs form factors can be extracted from Π µ (Γ ν , q) for appropriate choices of insertion and projector indices, . For the computation of the nucleon electric charge radius we use the first relation, which gives by far the best signal for G E (Q 2 ) and in addition allows to obtain a result directly at zero momentum transfer. The determination from the second relation involving G E (Q 2 ) is impeded due to the momentum factor appearing on the right-hand side of the equation and similarly in the last equation for G M (Q 2 ). A direct method for derivatives of form factors at zero momentum has been put forward in [26], based on the algebraic definition of the momentum derivative and quark propagator expansion, which promotes the n-point correlation function by one point per derivative and is thus very costly. In Ref. [27] we have explored model-independent position space methods to remedy this issue for G M (Q 2 ) without resorting to fits. We remark that one of these methods called momentum elimination in the plateauregion is similar to the approach we will introduce in the next section to allow for a direct, model-independent computation of the nucleon electric charge radius, which is otherwise hindered by the discrete nature of momenta on the lattice.

D. Summation method
In nucleon structure calculations it is notoriously difficult to reach ground state dominance due to an exponential signal-to-noise problem, hence a careful analysis of the corresponding systematics is required. While we have lattice data available for five values of t sep as listed in Table I it is not a priori clear that this is sufficient to control excited state effects. Therefore, we use the so-called summation method [28][29][30] in addition to the direct plateau method, which allows for a stronger suppression of excited states where t ex = 2a for the conserved vector current insertion. The contribution from the next higher state with a mass gap ∆ is now suppressed by an additional factor e −∆tsep in contrast to the plateau method for which the suppression is only ∼ e −∆t .
The mean-square charge radius of the nucleon r 2 E is defined from the electric Sachs form factor in Eq. (2) through expansion around small Q 2 Since for the proton and the isovector combination one has G p,u−d E (0) = 1, the mean square charge radius can be extracted from the expansion via computing For the neutron the leading term in Eq. (20) is absent and the normalization factor G n E (0) = 0 is dropped in the definition to obtain a finite result, hence r 2 E n can be computed in the same way from Eq. (21). Moreover, for the proton and the isovector combination it is possible to define a root mean square charge radius, i.e.
which is not meaningful for the neutron as its mean-square charge radius r 2 E n is negative. Computing r 2 E requires an evaluation of the derivative w.r.t. Q 2 which is not directly possible in a finite box as only finite, discrete momenta are accessible. In the literature several methods have been used to circumvent this issue, e.g. dipole fits. However, any such method introduces model-dependence and potentially large, systematic uncertainties considering the physical values of momenta that are typically available in lattice simulations and the steepness of G E (Q 2 ) close to zero momentum transfer. For a recent review on common methods to parametrize electromagnetic form factors we refer to Ref. [31]. In this study we aim to evaluate the pertinent derivative in Eq. (21) by an integral method for a suitably defined form factor in the small momentum region given below in Eqs. (28), (29). It allows for systematic probing of the dependence on the available supporting lattice data and has a well defined infinite volume as well as continuum limit without any further adjustments of parametrization. We present the calculation in its simplest 1-dimensional form based on (numerical) form factor data for on-axis 3-momenta ( q ∝ e k ). It can be extended to incorporate arbitrary momentum directions, which, however, requires taking into account anisotropy effects and analytic continuation. In the continuum we have Q 2 = −2m N (E N (q) − m N ) given our momentum setup for the nucleon at source and sink, and we can straightforwardly replace the Q 2 derivative by for on-shell E N = m 2 N + q 2 . In particular we can use   Thus on the lattice we define the form factor In Eq. (25) q max denotes the maximal value of on-axis momentum, for which numerical data for the ratio R µ and thus the form factor Π 0 can be obtained given the achieved statistical precision for two-and threepoint correlation functions. Apart from the dependence on the source-sink time separation t sep , checking the saturation of the integral defining r 2 E under variation of q max will be a major point of the systematic error analysis. Together with Π(q) in Eq. (25) we obtain its counterpart in position spaceΠ(n) from the discrete Fourier transform for −N/2 ≤ n ≤ N/2 with n = y/a, N = L/a andΠ(−y) =Π(y). Sample data for the position space form factorΠ(y, q max ) with varying cut-off q max in the discrete Fourier transform is shown in Figure 1.
where T 2n are Chebyshev polynomials of the first kind, this gives and with respect to Eq. (24) where P n (k 2 ) = (−1) n T 2n (k/2) − 1 /k 2 are polynomials ink 2 of degree n − 1. The difference quotient in Eq. (29) is an analytic function with a well-defined infinite volume and continuum limit and gives in particular Since for the Chebyshev polynomials T 2n (x) = (−1) n+1 (2n) 2 2 x 2 + O x 4 we also have the familiar integral / summation form or in terms of the original data Π(q) in Eq. (25) The kernel weight w in Eq. (32) for the lattice setup used in this study and exemplary data for the integrand in Eq. (33) are shown in Figure 2. Inspection of the integrands in the right-hand panel of Fig. 2 shows, that the contribution from the form factor at on-axis momenta beyond q/(2π/L) = 5 ∼ 7 will not be significant given our currently achieved statistical accuracy of the form factor.

IV. RESULTS
The lattice data for the electric Sachs form factor that are used as input in our calculation of the nucleon radii are shown in Figure 3 as a function of Q 2 . We have included data from the ratio method in Eq. (15) at all five available source-sink time separations as well as from the summation in Eq. (19). Results are shown for the proton, neutron and the isovector combination. In any case, we observe that the summation method yields results compatible with the ratio method for the largest value of t sep .
In Figure 4 we show our lattice data for the difference −6(G u−d E (Q 2 ) − G u−d E (0))/Q 2 as a function of Q 2 in physical units for all five values of t sep /a and the summation method. The extrapolation bands are computed by evaluating Eq. (29) for any given value of and multiplying the expression byk 2 /Q 2 leading to where Q 4 4m 2 + Q 2 = q. Note that the factor related to substitutingk 2 → Q 2 is one for the nucleon radius i.e. at Q 2 =k 2 = 0. However, at non-zero Q 2 , the change of variablek 2 → Q 2 has to be taken into account explicitly to make contact with the discrete lattice data, which is achieved by the above expression. Here we have consistently chosen a rather small value of q max = 4 · (2π/L) for the extrapolations in all six plots. Still, the resulting bands describe well the lattice data for Q 2 0.5 GeV. Data points entering the construction of the extrapolation are shown in red and fall onto the curve by definition. At larger momenta we observe the expected deviation from the lattice data. This behavior is particularly visible for the data at smaller values of t sep which are statistically more precise.

A. Study of the large-Q 2 tail contribution
In order to further study the systematics associated with the truncation of the Fourier transform for Π(q) we have explored ways of generating synthetic data for the tail of the form factor to avoid this truncation altogether. To this end we use several fit models that are applied to the lattice data at low values of Q 2 . In a second step we use the jackknife samples for the fit parameters to generate synthetic data samples for values of Q 2 corresponding to on-axis momenta q = (q, 0, 0) T with q > q max . This synthetic data are then used to perform the Fourier transform of Π(q) for all q > q max in addition to the actual lattice data at q ≤ q max . A simple choice for the proton and isovector form factor is the dipole fit model where A, B are free parameters of the fit. For the neutron we use a Galster-like parametrization [32,33] instead of the dipole fit model. This model takes into account that G n E (0) = 0 The free parameters of this model are denoted by C and D, while Λ 2 = 0.71 GeV 2 is treated as a constant.
Since the dipole fit model approaches zero quadratically for increasing values of Q 2 one might anticipate that using synthetic data from this model yields a result not very different from the one obtained by a simple truncation in the Fourier transform of Π(q) at a reasonable value of q max . Therefore, we consider an 1/ log fit model with an unphysically slower decrease at large values of Q 2  with two fit parameters E and F , which allows us to further investigate the dependence on the choice of the model. This model can be applied to all three isospin combinations, but it requires a careful choice of the lower bound of the fit range. In Figure 5 results are shown for fitting the different models to the proton and neutron data for G E (Q 2 ). Corresponding results for the isovector insertion are very similar to the ones for the proton. As expected, the bands from the dipole and Galster-like fits exhibit a much sharper drop-off in the large-Q 2 tail of the form factor than the corresponding bands from the 1/ log-fit. In general, fit ranges have been chosen such that Q 2 ∈ 0 GeV 2 , 1.8 GeV 2 for the dipole and Galster-like fits and Q 2 ∈ 0.7 GeV 2 , 1.8 GeV 2 for the 1/ log-fit. We found this choice to yield good values of χ 2 /d.o.f. regardless of the source-sink time separation as well as for data from the summation method. As mentioned before, the results for the fit parameters can be used to generate synthetic data for any desired value of Q 2 in the large-Q 2 tail of G E (Q 2 ). This allows us to compute the position space form factorΠ(y, q max ) for any value of q max substituting missing lattice data in the Fourier transform by synthetic data from a given fit model. Typical results from this procedure are shown in Figure 6. Regardless of the model, the signal becomes increasingly peaked at y = 0 for larger values of q max while the large-y tail becomes flatted out, which is the expected behavior. Moreover, including data generated from the 1/ log fit model enhances this effect further at large values of q max compared to the dipole or Galster-like fit. Some results for the final extrapolations together with lattice and model data as a function of Q 2 are shown in Figure 7 for all three isospin combinations and different fit models. We find that there is hardly any difference visible between the extrapolations obtained from using synthetic data from the dipole or Galster-like fit model in the left panels and the 1/ log fit model in the right panels of Figure 7. In the presented examples we have employed model data for either q ≥ 6 · (2π/L) or q ≥ 5 · (2π/L) depending on whether results have been obtained at finite values of t sep or from the summation method, respectively. At any rate, this confirms the expectation from the kernel weight function in Figure 2 that the final result is indeed saturated by the first few values of on-axis momentum q. Therefore, residual systematic effects related to the truncation of the Fourier transform must be small. (0))/Q 2 together with the continuous extrapolation band computed from Eq. (35). The first, second and last row show results for the isovector data at tsep = 1.29 fm, the proton data from the summation method and the neutron data at tsep = 0.97 fm, respectively. In the left column the dipole fit and Galster-like fit models have been used in the generation of synthetic data, while for the right column synthetic data from a 1/ log-fit have been used. For the isovector combination and the neutron synthetic data has been used for q ≥ 6 · (2π/L) in the construction of the extrapolation, while for the more noisy proton data from the summation method the switch to model data occurs at q ≥ 5 · (2π/L).

B. Further systematics and final values
The smallness of systematic effects due to truncation of the Fourier transform is further corroborated by a direct comparison of the results for r u−d,p E and r 2 E n at different cutoffs q max as well as using synthetic data from the fit models. An overview of results for all available values of t sep together with the summation method is shown in Figure 8. The corresponding numerical values are listed in Table II together with the mean squared radii for the proton and isovector combination. In general, the results using data from the two fit models are in excellent agreement at any given value of t sep . A significant difference is only observed between the two smallest values of the cutoff q max = 4·(2π/L) and q max = 5·(2π/L) for the statistically most precise data at the smallest three source-sink time separations. However, the results obtained with a cutoff of q max = 5 · (2π/L) in the Fourier transform are compatible with the results obtained from modeling the large-Q 2 tail of the form factor.   Regarding excited state contamination there is a weak trend visible for the first few values of t sep . We find that for the proton and isovector radius the results at the largest available value of t sep = 1.62 fm are in good agreement with the summation method. Therefore, we quote the value obtained at the largest available value of t sep and for q max = 5 · (2π/L) as our final results for the isovector and proton radius, i.e.
where the first error is statistical and the second error refers to the scale setting uncertainty. At the current level of precision the statistical error clearly dominates and systematics due to the scale settings are negligible. The data for the neutron combination is more noisy and the signal for r 2 E n is essentially lost after the third source-sink time separation. For that reason we quote the value obtained at t sep = 1.29 fm and q max = 5·(2π/L)  as the final result

C. Comparison with other studies
Although electromagnetic form factors have been studied within lattice QCD since many years, it is only recently that they have been extracted using simulations with physical values of the light quark masses, referred to as physical point ensembles. Also, while there are a number of studies for the isovector combination, namely, the difference between proton and neutron form factors, for the proton and neutron themselves the results are scarce. This is due to the complexity of evaluating accurately contributions from disconnected quark loops. In fact, only ETMC has computed quark-disconnected contributions using physical point ensembles [14]. Therefore, we compare the results of our direct method with recent results extracted using simulations with nearly physical pion masses neglecting the quark-disconnected contributions. Such contributions yield a correction at the 15% level for the neutron electric form factors, while for the proton the correction is only at the one percent level. We summarize below recent simulations used by various collaborations for the extraction of the nucleon electromagnetic form factors.
• ETMC analyzed three physical point ensembles [14,34] of twisted mass fermions: a N f = 2 + 1 + 1 ensemble with m π L = 3.62 which has also been used in the present study, and two N f = 2 ensembles with m π L = 2.98 and m π L = 3.97 and the same lattice spacing of a = 0.0938(3)(1) fm. In any case, results are automatically O(a)-improved, hence cut-off effects are of O(a 2 ).
• LHPC analyzed two N f = 2 + 1 ensembles simulated with two step HEX-smeared clover fermions: an ensemble with m π = 149 MeV, lattice spacing a = 0.116 fm and m π L = 4.21 [35] and one ensemble at m π = 135 MeV and m π L = 4 with a finer lattice spacing of a = 0.093 fm [36]. In the latter study they employed a momentum derivative method to extract directly the radii. They use improved currents, therefore their cut-off effects are O(a 2 ).
• The PNDME collaboration [39] analyzed eleven ensembles of N f = 2 + 1 + 1 highly improved staggered quarks (HISQ) simulated by MILC. They used a mixed action setup with clover fermions in the valence sector. The ensembles have lattice spacings a 0.06, 0.09, 0.12, 0.15 fm and the pion masses are m π 135, 225, 315 MeV. Using a combined fit analysis, they performed a chiral, continuum and infinite volume extrapolation. We limit ourselves here to their results obtained using the m π = 135 MeV ensemble with a = 0.0570(1) fm and m π L = 3.  [14], and N f = 2 and mπL 3 (orange upwards-pointing triangles) [34]; LHPC using an N f = 2 + 1 ensemble with mπL = 4.21 [35] (brown downwards-pointing triangles) and an N f = 2 + 1 ensemble with mπL = 4 [36] (square magenta); PACS using an N f = 2 + 1 ensemble with mπL = 6 [37] (cyan rhombus) and an ensemble with mπL = 7.4 [38] (purple pentagon); PNDME [39] using an N f = 2 + 1 + 1 mixed action and their physical point ensemble with mπL = 3.7 (black crosses). The experimental result extracted from the muonic hydrogen [3] is given by the vertical dashed-dotted line and the value referenced by CODATA [40] by the dotted vertical line. The PDG value [41] is given by the dashed vertical line. We note that the PRad [7] experiment finds compatible value for the proton radius with the muonic hydrogen experiment albeit with much larger error, thus has not been included in the figure.
In Fig. 9 we compare the results for the radii as extracted from this study with the aforementioned studies. For the isovector combination we find that the value computed within our direct method is compatible with all other studies. We note that one of the three ETMC values is extracted by analyzing the same ensemble but using a dipole fit to extract the radius instead. For the proton we find very good agreement with all other determinations except the one from PACS, using m π L = 7.4 [38] ensemble, finding a higher value than the rest. In the PACS study the analysis on the large volume is carried out using only small time separations. Whether this larger value, which is consistent with the experimental determinations, is due to the usage of larger volume needs further investigation. Clearly, however, all lattice QCD results have errors that are compatible with the difference between the experiment result from the µH experiment and that from CODATA. Therefore it is not yet possible to make a distinction between these two values based on lattice results. For the neutron, the relative errors on lattice data are much larger and there is overall agreement among lattice results. The value extracted in this work is in fact very close to the PDG value.

V. SUMMARY & OUTLOOK
In this work we extract the r.m.s charge radii avoiding a fit ansatz to the electric form factors of the nucleon, that may introduce a model error. The method has been shown to be insensitive towards the large-Q 2 tail of the form factor by testing different ansätze to model the large Q 2 -dependence of the form factors where lattice QCD data are not available. In particular, even using an unphysical ansatz that falls very slowly with Q 2 the results are unaffected, demonstrating that lattice QCD results up to about Q 2 ∼ 2 GeV 2 determine the radii for the currently available lattice volume and statistical precision. Comparing the results for r p−n E , r p E and r n E extracted by applying this approach to the values when the traditional approach of fitting the electric form factors is used for the same ensemble we find consistent values. Moreover, our values agree with the model-independent determination using the approach of algebraic derivative and quark propagator expansion for the direct calculation of form factor derivatives [26] employed by LHPC [36]. In all approaches the errors on the proton charge radius are still large and cannot distinguish between the values extracted from muonic hydrogen and older electron scattering experiments. We plan to implement this approach to study the magnetic radii and moments where an increased precision in the lattice QCD data will be required. For the neutron, one has also to include the disconnected contribution which may bring the value closer to the experimental one. Similarly, disconnected diagrams will have to be included if aiming for 2% precision on the proton radius. Finally, non-zero lattice spacing and finite volume artifacts need to be evaluated. This can only take place when ensembles at physical quark mass values are simulated for smaller lattice spacings and larger volumes. Such a program will be possible as these ensembles are simulated and analyzed.