Computing the nucleon charge and axial radii directly at $Q^2=0$ in lattice QCD

We describe a procedure for extracting momentum derivatives of nucleon matrix elements on the lattice directly at $Q^2=0$. This is based on the Rome method for computing momentum derivatives of quark propagators. We apply this procedure to extract the nucleon isovector magnetic moment and charge radius as well as the isovector induced pseudoscalar form factor at $Q^2=0$ and the axial radius. For comparison, we also determine these quantities with the traditional approach of computing the corresponding form factors, i.e. $G^v_E(Q^2)$ and $G_M^v(Q^2)$ for the case of the vector current and $G_P^v(Q^2)$ and $G_A^v(Q^2)$ for the axial current, at multiple $Q^2$ values followed by $z$-expansion fits. We perform our calculations at the physical pion mass using a 2HEX-smeared Wilson-clover action. To control the effects of excited-state contamination, the calculations were done at three source-sink separations and the summation method was used. The derivative method produces results consistent with those from the traditional approach but with larger statistical uncertainties especially for the isovector charge and axial radii.


I. INTRODUCTION
The experimental determinations of the proton (electric) charge radius r p E have a discrepancy greater than 5-sigma between the value determined from spectroscopy of muonic hydrogen [1,2] and the CODATA average [3] of experimental results obtained from hydrogen spectroscopy and electron-proton scattering.This presently unresolved "proton radius puzzle" is the focus of various theoretical and experimental efforts 1 .Last year, the CREMA collaboration reported on their study of muonic deuterium [5].Their experiment corroborates the muonic hydrogen result for the proton charge radius, while finding a similar 6-sigma discrepancy for the deuteron charge radius with the CODATA values, and a 3.5-sigma discrepancy to electronic deuterium spectroscopy results [6].Thus, having a reliable ab-initio calculation of the proton charge radius is a highly attractive goal for practitioners of lattice QCD.
The conventional approach for determining quantities like the charge radius on the lattice involves the computation of form factors at several different discrete values of the initial and final momenta, p and p , that are allowed by the periodic boundary conditions, followed by a large extrapolation to zero momentum transfer Q 2 = 0.This introduces a source of systematic uncertainty, analogous to the systematic uncertainty associated with the choices of the fit ansatz and range of Q 2 in extracting the proton charge radius from electron-proton scattering data.Systematic errors of this kind have in fact been proposed as a possible explanation of the radius puzzle [7][8][9].Given that the smallest nonzero value of Q 2 accessible on the largest available lattices is still an order of magnitude higher than in scattering experiments [10], a lattice method for computing r p E and similar observables directly at Q 2 = 0 without the need of a shape fit is highly desirable.
The Rome method, presented in Ref. [11], provides a way to calculate the momentum derivatives of quark propagators on the lattice at zero momentum.This enables calculating the momentum derivatives of the correlation functions at zero momentum and obtaining the form factors and their momentum derivatives at vanishing momenta.To this end, one introduces twisted boundary conditions and takes the symbolic derivative(s) with respect to the twist angle (at zero twist angle) before the numerical evaluation of the path integral over the gauge fields.
For the case of a pion, it was shown in Ref. [12] that the Rome method for momentum derivatives could be used to calculate the pion charge radius with finite-volume effects that are exponentially suppressed, with asymptotic behavior ∼ √ m π L e −mπL .We employ the Rome method for extracting the proton isovector charge radius (r 2 E ) v and the isovector magnetic moment µ v = G v M (0), from matrix elements of the vector current.We also extract the proton axial radius r 2 A and the induced pseudoscalar form factor at zero momentum, G P (0), using nucleon matrix elements of the axial current.We compare the results from the derivative method with those from the traditional approach.
The outline of the paper is as follows.We start by reviewing the electromagnetic and axial form factors in Sec.II.Section III is devoted to describing the traditional approach for isolating the nucleon ground state, extracting the nucleon electromagnetic and axial form factors and the fits to the Q 2 -dependence of the form factors using the z expansion to determine the corresponding radii and form factors at Q 2 = 0. Section IV explains in detail the derivative method for computing the momentum derivatives of matrix elements at Q 2 = 0 using the Rome method, which we use to determine the charge and axial radii in addition to the magnetic and induced pseudoscalar form factors directly at Q 2 = 0.In Sec.V, we describe the lattice methodology and the ensemble of configurations that are used.Finally in Sec.VI, we present our numerical results computed directly at Q 2 = 0, and compare them with the traditional approach.We give our conclusions in Sec.VII.

II. DEFINITIONS OF THE FORM FACTORS
The nucleon matrix elements can be parametrized in terms of nucleon form factors as where p, p are the initial and final nucleon momenta, λ, λ label the different polarization states, and u is the nucleon spinor.We are defining the form factors using a current of flavor q in a proton and | p, λ is a proton state.O q,µ X refers to either the vector (X = V ) or the axial (X = A) current.
For the case of the vector current, O q,µ V = qγ µ q, F q,µ V ( p, p ) can be written in terms of the Dirac and Pauli form factors, F q 1 (Q 2 ) and F q 2 (Q 2 ), in Minkowski space as: where m is the nucleon mass and Q 2 = −(p − p) 2 ≥ 0 is the momentum transfer.These form factors can also be expressed in terms of the nucleon electric G E (Q 2 ) and magnetic G M (Q 2 ) Sachs form factors via: The charge and magnetic radii, r 2 E,M , and the magnetic moment, µ, are defined from the behavior of For the axial vector current, O q,µ A = qγ µ γ 5 q, F q,µ A ( p, p ) can be expressed in terms of the axial and induced pseudoscalar form factors, G q A (Q 2 ) and G q P (Q 2 ), as: The axial form factor admits the following expansion for small momentum transfer where g q A is the axial-vector coupling constant and r q A is the axial radius.
In this work, we are considering the isovector electromagnetic Sachs form factors which parametrize the matrix elements of the u − d flavor combination between proton states and, neglecting the isospin breaking effects, are equivalent to the difference between the form factors of the electromagnetic current The isovector axial form factors G v A,P (Q 2 ) are defined in a similar way.

III. COMPUTATION OF MATRIX ELEMENTS USING THE TRADITIONAL METHOD
For determining the nucleon matrix elements in lattice QCD, we compute the nucleon two-point and three-point functions, In this section, we use Minkowski-space gamma matrices.Above, χ = abc (ũ T a Cγ 5 1+γ0 2 db )ũ c is a proton interpolating operator constructed using smeared quark fields q and Γ pol = 1 2 (1 + γ 0 )(1 + γ 3 γ 5 ) is a spin and parity projection matrix.The three-point correlators have contributions from both connected and disconnected quark contractions, but we compute only the connected part since for the isovector flavor combination the disconnected contributions cancel out.
We will be tracing the correlators with Γ pol which contains the projector (1 + γ 0 )/2 so that we can effectively write the overlap of the interpolating operator with the ground-state proton as Ω|χ α (0)| p, λ = Z( p)u( p, λ) α [13,14].At large time separations we obtain where ∆E 10 ( p) is the energy gap between the ground and the lowest excited state with momentum p.By taking τ and T − τ to be large, unwanted contributions from excited states can be eliminated.In order to compute C 3 , we use sequential propagators through the sink [15].This has the advantage of allowing for any operator to be inserted at any time using a fixed set of quark propagators, but new backward propagators must be computed for each source-sink separation T .Increasing T suppresses excited-state contamination, but it also increases the noise; the signal-to-noise ratio is expected to decay asymptotically as e −(E− 3 2 mπ)T [16].In order to cancel the overlap factors and the dependence on Euclidean time, we define the normalization ratio, R X N , and the asymmetry ratio, R S as and compute their product as a function of τ ∈ [0, T ] with fixed T .Above, and ∆E min = min{∆E 10 ( p), ∆E 10 ( p )}.
The ratio in Eq. ( 16) gives an estimate of the nucleon matrix element p , λ |O q,µ X |p, λ and produces at large T a plateau with "tails" at both ends caused by excited states.In practice, for each fixed T , we average over the central two or three points near τ = T /2, which allows for matrix elements to be computed with errors that decay asymptotically as e −∆EminT /2 .
Improved asymptotic behavior of excited-state contributions can be achieved by using the summation method [17,18] which requires performing the calculations with multiple source-sink separations.Taking the sums of ratios for each T yields where we choose τ 0 = 1 and c is an unknown constant.The matrix element can then be extracted from the slope of a linear fit to S q,µ X ( p, p , T ) at several values of T .The leading excited-state contaminations decay now as T e −∆EminT .For calculating the form factors - for the case of the vector current and G A (Q 2 ), G P (Q 2 ) for the case of the axial current -we construct a system of equations parameterizing the corresponding set of matrix elements at each fixed value of Q 2 [19].We combine equivalent matrix elements to improve the condition number [20].We find the solution of the resulting overdetermined system of equations by performing a linear fit.This approach makes use of all available matrix elements in order to minimize the statistical uncertainty in the resulting form factors.
The charge and axial radii can be extracted from the slopes of the electric and axial form factors at Q 2 = 0, respectively.For that we need to fit the Q 2 -dependence of each form factor.In order to avoid the model dependence included in the commonly used fit ansatzes, such as a dipole, we use the model-independent z expansion [21][22][23][24], where each form factor can be described by a convergent Taylor series in z which conformally maps the complex domain of analyticity in Q 2 to |z| < 1.We fix a 0 = 1 for fitting G E (Q 2 ) since G E (0) = 1.We use the particle production threshold t cut = (2m π ) 2 for the vector case and t cut = (3m π ) 2 for the axial case.We apply z-expansion fits following the approach of Ref. [25].The intercept and slope of the form factor at Q 2 = 0 can be obtained from the first two coefficients, a 0 and a 1 .We impose Gaussian priors on the remaining coefficients centered at zero with width equal to 5max{|a 0 |, |a 1 |}.We truncate the series with k max = 5 after verifying that using a larger k max produces identical fit results in our probed range of Q 2 .Furthermore, the isovector G P form factor has an isolated pole at the pion mass below the particle production threshold.We thus remove this pole before fitting and perform the z-expansion fit to (Q 2 + m 2 π )G P (Q 2 ).We perform correlated fits by minimizing with respect to {a k }, where S is an estimator of the covariance matrix and the last term augments the chi-squared with the Gaussian priors.For choosing the estimator of the covariance matrix, we use S = (1 − λ)C + λC diag , where λ = 0.1, C is the bootstrap estimate of the covariance matrix and C diag is the diagonal part of C.

IV. DERIVATIVE METHOD
In this section, we explain the details of our approach for extracting the nucleon charge radius directly at Q 2 = 0. We begin with reviewing the Rome method for computing the momentum derivatives of quark propagators in Subsection IV A. The flavor structure of the correlators constructed from the momentum derivatives of the quark propagators is investigated in IV B. In Subsection IV C, we show how to use the momentum derivatives of the quark propagator in order to obtain the first-and second-order derivatives of the nucleon two-and three-point functions with respect to the initial-state momentum p, and then obtain momentum derivatives of matrix elements in Subsection IV D. From the latter one can then extract the charge radius r 2 E , the magnetic moment µ = G M (0), for the case of the electromagnetic vector current, and the axial radius, r 2 A , and the induced pseudoscalar form factor at zero momentum, G P (0), for the case of the axial current.

A. Momentum derivatives of quark propagator
On a lattice with finite size and quark fields satisfying periodic boundary conditions, consider a generic correlation function C( p, t) depending on the three-momentum p and Euclidean time t, which after fermionic integration and Wick contractions can be written in terms of quark propagators and operator insertions as, where U are gauge links and P [U ] is the corresponding probabilistic weight in the functional integral.The planewave phase factor e −i p( x− y) can then be absorbed into one of the quark propagators, which results in a momentum dependent quark propagator G[x, y; U, p] = e −i p( x− y) G[x, y; U ]. G[x, y; U, p] can be obtained by solving the lattice Dirac equation with link variables rescaled by a phase factor: Carrying momentum in a propagator with a uniform U (1) background field is the same approach as used in a standard transformation of twisted boundary conditions [26,27].With p restricted to be a Fourier momentum in the finite volume, the above redefinition is exact.However, to obtain a momentum derivative, we must implicitly make use of twisted boundary conditions and allow p to be continuous.We use the expansion of the lattice Dirac operator and D[U, p]G[U, p] = 1 to compute the first-order momentum derivative of the propagator as: where we use the compact notation and similar notation for G(. . .; U, p).Multiplying Eq. ( 25) from the left by G ≡ D −1 leads to: Similarly, we can derive the second-order momentum derivative of the propagator: Using the lattice Dirac operator for the clover-improved Wilson action, the momentum derivatives of the propagators at a fixed gauge background become [11]: We drop U from the propagators for brevity.Γ k V and Γ k T are the point split vector and tadpole currents, respectively.Those are defined using Euclidean gamma matrices, γ k E , as: In the case of a smeared-source smeared-sink propagator (needed in the two-point function), the phase factor can be absorbed into the propagator in the following way: G(x, y; p) = e −i p( x− y) x ,y where K is the smearing kernel.The momentum derivatives can then be calculated using the product rule along with Eq. ( 29) and Eq.(30).Denoting the momentum derivative with for shorter notation, we obtain For the smeared-source point-sink propagator, which is needed for the three-point function and for evaluating Eq. ( 34) and Eq. ( 35), we obtain Organized in this way, we have one additional propagator right-hand-side per derivative.Gaussian Wuppertal smearing [28] is given by K(x, y; p) = x ,x ,... K0(x, x ; p)K0(x , x ; p)...K0(x ... , y; p) , with We use APE-smeared gauge links Ũ [29].The mth derivative of K 0 at zero momentum is equal to Thus, the first-and second-order momentum derivatives of smearing with ) .

B. Flavor structure of correlators constructed from propagator derivatives
In cases where derivatives of nucleon two-point functions need to be evaluated, there is an ambiguity in applying the above procedure: there are three quark propagators, and the momentum could be absorbed into any of them.To resolve this issue, we make explicit use of twisted boundary conditions, with the understanding that before computing any correlation functions we will take the derivative with respect to the twist angle, at vanishing twist angle.
We introduce a third light quark, r, with the same mass as u and d but with twisted boundary conditions, and a corresponding ghost quark that cancels its fermion determinant.The three light quarks {u, d, r} contain an approximate SU(3) flavor symmetry that becomes exact when the twist angle is zero, or in the infinite-volume limit.Under this symmetry group there is a baryon octet that contains the ordinary (untwisted) nucleons, as well as states with one or two r quarks.We are interested in the states with one r quark, and we find that there are two kinds: an isospin singlet and a triplet, the Λ r and Σ r , respectively.This was previously discussed in Ref. [30].
For the states with quark content udr we use interpolating operators where O = ūΓr is a quark bilinear and X is Σ r or Λ r .The initial momentum p is implied in the initial state due to the twisted boundary conditions for the r quark.The ground-state contribution is proportional to the matrix element N ( p )|O|X( p) for which we will evaluate ∂ ∂ p at p = p = 0.In practice, we simply use our already coded expressions for the connected diagrams in the nucleon three-point functions C q 3 with O q = qΓq, q ∈ {u, d}, and replace the propagator connecting the nucleon source and O q with a first-or second-derivative propagator.By comparing the contractions, we find the relations where the r propagator is substituted into the evaluation of the right-hand-side expressions as described above.A similar consideration was made in Ref. [30]; these relations could also be derived from SU(3) symmetry.
When forming ratios, we must use the appropriate two-point functions: taking Eq. ( 16) with the three-point function C X→N   3   , all nucleon two-point functions that take the initial-state momentum p must be replaced by the two-point function for state X.Once we have formed the ratios for the X → N matrix elements, we can invert the relations in Eq. ( 42) to obtain the nucleon matrix elements of O u and O d .

C. Momentum derivatives of the two-point and three-point functions
Let us consider the two-point function of the isospin singlet operator, χ Λr = 3 2 [udr].This can be written in terms of smeared-source smeared-sink quark propagators, G, as where f ββγδ is the spin tensor determining the quantum numbers of the Λ r and G(x, 0; p) = e −i p x G(x, 0).By using the first-and second-order momentum derivatives of a quark propagator at zero momentum given in Eq. ( 29) and Eq. ( 30), one can straightforwardly calculate the momentum derivatives of the two-point correlators.
For connected diagrams, the three-point function with current O Γ = qΓq and zero sink momentum p = 0 can be written as where G refers to a propagator with smeared source and point sink and G S (y) is the sequential backward propagator, which is independent of p.Only the forward propagator G(y, 0; p) needs to be expanded using Eq. ( 29) and Eq. ( 30).Hence, no additional backward propagators are needed.Figure 1 shows graphically the way we compute the momentum derivatives of the correlation functions on the quark level.The derivative method cannot be applied to disconnected diagrams because those involve a quark propagating from a point to the same point and therefore the momentum transfer can not be absorbed into the propagator.

D. Momentum derivatives of the ratio
Because we do not know how Z( p) depends on the momentum, we need to compute the momentum derivatives of the ratio of three-point and two-point functions given in Eq. ( 16).Here and in the following, we use Minkowski-space gamma matrices.We set p = 0 and p = k e j , where e j is the unit vector in j-direction.For computing the first-and second-order momentum derivatives of the ratio in Eq. ( 16), we start by computing the momentum derivatives of the normalization ratio part, R X N , defined in Eq. ( 14): where, for more readability we suppress the τ, T parameters as well as O µ X from the correlation functions and the ratio.We denote the derivatives with a prime, e.g.C 2 (k) ≡ dC2 (k)  dk .We know that C 2 (0) = 0 in the infinite-statistics limit because of parity symmetry.Hence, we can eliminate this from the ratios.Similarly, we can calculate R S (k) and R S (k) which can be used together with Eq. (45) and Eq. ( 46) to calculate the first-and second-order derivatives of the ratio R X .These derivatives are computed on the lattice directly at k = 0 as discussed earlier in the previous section.
From the ground-state contributions to the correlation functions given in Eq. ( 12) and Eq. ( 13), we find the following ground-state contribution to their ratio We take the derivative with respect to k and obtain: (R X ) (k) can be calculated in a similar way.We use the continuum dispersion relation , and find that at k = 0, the second derivative is needed to obtain the slope of F 1 : The same applies for F 2 , G A , and G P .Furthermore, we have at k = 0 : For the renormalized vector current, we use G E (0) = 1 and find nonzero results for the following combinations of j and µ: and for the axial current, with ∂ j = ∂ ∂p j and From Eq. (54) and Eq.(55), we find the following relations for the nucleon magnetic moment µ = G M (0) and squared charge radius r 2 E : The red and blue bands correspond to the combined fits of C 2 (0, t) Λ,Σ and C 2 (0, t).
0.00 0.05 0.10 0.15 FIG. 3: The derived values for Z( p 2 ) from two-state fits of C 2 ( p, t) (black points) followed by a linear fit (grey band) for extracting Z(0) and Z (0).
where Z ≡ Z(0).We expect C 2 t) to vanish due to parity symmetry and our numerical results shown in the left part of Fig. 2 confirm that, which allows us to set Z (0) = 0 in Eq. ( 70).We apply a combined 1-state fit for C 2 (0, t) and C 2 (0, t) Λ,Σ using Eq.(68) and Eq. ( 70) with Z, Z and m being the fit parameters.The results of these fits are shown in Fig. 2, where the slight differences between the momentum derivatives of Σ r and Λ r two-point functions give an indication of the systematic errors associated with the derivative method and motivate the approach described in Sec.IV B for isolating Σ r → N from Λ r → N three point functions when extracting the momentum derivatives of the matrix elements.We also try another approach for extracting Z(0) and Z (0) where we apply two-state fits to C 2 ( p, t) for different discrete values of p 2 which allows us to extract Z( p 2 ).The extracted values for Z( p 2 ) are consistent with a linear dependence on (a p) 2 .By applying a linear fit to Z( p 2 ) against p 2 , Z(0) can be obtained from the intercept and Z (0) from the slope as Z (0) = 2 ∂Z( p 2 ) ∂ p 2 .This is shown in Fig. 3. Table I reports a comparison between the extracted values for Z(0) and Z (0) using the two different approaches and when using [C 2 (0, t)] Σ and [C 2 (0, t)] Λ in the combined fit.All fit methods lead to consistent values for both Z(0) and Z (0).FIG.4: Isovector magnetic moment (left) and isovector charge radius (right).For both µ v and (r 2 E ) v /a 2 , results from the ratio method are shown using source-sink separations T /a ∈ {10, 13, 16}, as well as the summation method.

B. Electromagnetic form factors
The "plateau plots" in Fig. 4 show the results we obtain using the momentum derivative approach for both G v M (0) (left), computed using Eq.(60), and (r 2 E ) v /a 2 (right), extracted from Eq. (61).In each case, we show results from both the ratio method and the summation method.G v M (0) increases for increased source-sink separations, indicating that the excited-state contributions are significant in this case.The relative statistical uncertainty is much larger for (r 2 E ) v /a 2 , and therefore we are unable to resolve any significant excited-state effects.
Figure 5 shows a comparison between our results using the derivative method and the traditional approach for both the isovector magnetic moment µ v = G v M (0) (bottom row) and the isovector charge radius (r 2 E ) v (top row).In Fig 5, we present the results extracted using both the ratio method with the smallest source-sink separation T /a = 10 and the summation method.When going to the summation method, G v E (Q 2 ) decreases significantly whereas G v M (Q 2 ) increases (especially for small Q 2 ) towards the corresponding phenomenological curve from Kelly [36].This shows the non-trivial contribution from excited states associated with the ratio method using T /a = 10.The summation points for G v E (Q 2 ) lie slightly above the corresponding Kelly curve while those for G v M (Q 2 ) show a good agreement with the Kelly curve.The derivative method's results for both G v M (0) and (r 2 E ) v using the summation method are consistent with both the traditional method's results and the experiment but with statistical errors roughly twice as large as the traditional approach for the isovector magnetic moment and three times as large for the isovector charge radius, as reported in Table II.

C. Axial form factors
The left-hand side of Fig. 6 shows the isovector induced pseudoscalar form factor G v P (0) extracted using the derivative method, Eq. ( 64).The right-hand side of the same figure shows the extracted r 2 A using Eq.(63). Figure 6 shows the plateau plots for both quantities corresponding to the three available source-sink separations in addition to the summation points.For G v P (0), we see a large increase with the source-sink separation, indicating substantial excited-state effects, and that leads us to conclude that the summation point may not be free from excited-state effects.For r 2 A , the statistical errors are too large to detect any excited-state effects.FIG.5: Isovector electric (top row) and magnetic (bottom row) form factors using both the ratio method with T = 10 a (left column) and the summation method (right column).The blue points show results from the standard method and the red bands show a z-expansion fit to those points.The green band (top) and point (bottom) show the slope and value of the respective form factor at Q 2 = 0, computed using the momentum derivative method.The black curves result from a phenomenological fit to experimental data by Kelly [36].A comparison between our results using the derivative method and the traditional method for both r 2 A and G v P (0) is shown in Fig. 7, top and bottom row, respectively.Shown are results from both the ratio method with T /a = 10 and the summation method.Both G v A (Q 2 ) and G v P (Q 2 ) increase when going to the summation method indicating the significant excited-state contributions for the ratio method with T /a = 10.The extracted value for the axial radius using the derivative method has a much larger statistical error compared to its value from the traditional approach.For G v P in Fig. 7, before fitting we remove the pion pole that is present in the form factor, and then restore it in the final fit curve as was discussed in Sec.III.At T /a = 10, there is a significant disagreement between G P (0) from the traditional and the derivative approaches which is likely due to excited-state effects.The value for G v P (0) using the summation method and the derivative approach seems to be in good agreement with its value from the traditional approach despite the large extrapolation caused by the inclusion of the pion pole in the fit.However, G v P (0) obtained from the derivative method has statistical uncertainties roughly twice as large as the traditional approach.Our results for the axial form factors are reported in Table II.

VII. SUMMARY AND OUTLOOK
In this paper, we presented a derivative method for computing nucleon observables at zero momentum transfer.Using this method helps to avoid the model dependence and the large extrapolation needed in the traditional approach for computing such quantities.We applied the derivative method to the nucleon isovector magnetic moment and electric charge radius as well as the isovector induced pseudoscalar form factor at Q 2 = 0 and the axial radius.
We confirm that our approach produces results consistent with those obtained using the traditional method.For G M (0) and G P (0), there is excellent agreement between the two approaches.This is particularly remarkable in the latter case, since the pion pole produces a very large effect in the extrapolation of G P (Q 2 ) to Q 2 = 0.However, we found that this approach suffers from large statistical uncertainties, especially for the isovector charge and axial radii.This may be connected with the fact that these quantities require a second momentum derivative.However, G P (0) also requires two derivatives and is not as noisy.Our quoted errors are statistical; we still need to estimate and improve control over systematic uncertainties in order to have a reliable calculation.The difference between the CODATA value of (r 2 E ) v and its muonic hydrogen measurement is ∼ 0.06 fm 2 , so it will be a challenge to calculate the charge radius with a total uncertainty significantly less than that.
Our present setup of the derivative method includes computing the momentum derivatives of the nucleon correlators with respect to only the initial nucleon momentum.As suggested originally for the pion charge radius in Ref. [12], one can alternatively obtain the radius by computing the mixed-momentum derivatives of three-point functions i.e., first-order momentum derivatives with respect to both initial-and final-state momenta.A calculation including this alternative approach is currently underway; preliminary results suggest that the statistical uncertainty for the radii is significantly reduced [37].

FIG. 1 :
FIG.1: Left: Nucleon two-point (top) and three-point (bottom) functions .The solid black circles represent the nucleon source and sink, the black square in the three-point function represents the current insertion.The red line refers to the propagator which we use for computing the momentum derivatives of the correlators which carry therefore the derivative vertex (solid red circle).The right panel shows the representation of the derivative vertex for the simplified case of unsmeared propagators.

TABLE I :
Resulting values for Z(0) and Z (0) using either the combined fit of C 2 (0, t) and C 2 (0, t) Λ,Σ or the fit to Z( p 2 ).

TABLE II :
Numerical results for the four different nucleon observables at Q 2 = 0, computed with the traditional method (via z expansion fit to the form factor shape) and with the derivative method.