Calculation of derivative of nucleon form factors in $N_f = 2+1$ lattice QCD at $M_\pi = 138$ MeV on a (5.5 fm)$^3$ volume

We present a direct calculation for the first derivative of the isovector nucleon form factors with respect to the momentum transfer $q^2$ using the lower moments of the nucleon 3-point function in the coordinate space. Our numerical simulations are performed using the $N_f = 2 + 1$ nonperturbatively $O(a)$-improved Wilson quark action and Iwasaki gauge action near the physical point, corresponding to the pion mass $M_\pi =138$ MeV, on a (5.5 fm)$^4$ lattice at a single lattice spacing of $a = 0.085$ fm. In the momentum derivative approach, we can directly evaluate the mean square radii for the electric, magnetic, and axial-vector form factors, and also the magnetic moment without the $q^2$ extrapolation to the zero momentum point. These results are compared with the ones determined by the standard method, where the $q^2$ extrapolations of the corresponding form factors are carried out by fitting models. We find that the new results from the momentum derivative method are obtained with a larger statistical error than the standard method, but with a smaller systematic error associated with the data analysis. Within the total error range of the statistical and systematic errors combined, the two results are in good agreement. On the other hand, two variations of the momentum derivative of the induced pseudoscalar form factor at the zero momentum point show some discrepancy. It seems to be caused by a finite volume effect, since a similar trend is not observed on a large volume, but seen on a small volume in our pilot calculations at a heavier pion mass of $M_{\pi}= 510$ MeV. Furthermore, we discuss an equivalence between the momentum derivative method and the similar approach with the point splitting vector current.


I. INTRODUCTION
A discrepancy of experimental measurements of the proton charge radius, called the proton radius puzzle, has not been solved yet since the muonic hydrogen measurement was reported in 2010 [1]. The values of the charge radius measured from both elastic electronproton scattering and hydrogen spectroscopy agree with each other [2], while they differ from the one measured from the muonic hydrogen. Several experiments are carried out and also proposed to understand this discrepancy, see Ref. [3] for a review of this puzzle.
Lattice QCD calculation, which is a unique computational experiment to investigate the complicated strong interaction dynamics, based on the first principles of QCD, can tackle the problem as an alternative to actual experiments.
In lattice QCD calculation, the mean square (MS) charge radius can be determined from the slope of the electric nucleon form factor at q 2 = 0. Most calculations including our previous works [4,5] and recent works [6][7][8][9][10][11][12][13][14][15][16][17][18] measure the form factor in q 2 = 0 with discrete lattice momenta, and fit the data with appropriate functional forms, such as dipole form, in order to determine the form-factor slope at the zero momentum point. The choice of fit functions, however, gives rise to relatively large systematic error. Even in a fit using data set including at tiny q 2 data point, its systematic error still remains as large as the statistical error of 2% as presented in Ref. [5]. Such systematic uncertainty needs to be reduced to draw any conclusion on a maximum discrepancy of about 4% observed in the three experiments.
Apart from the systematic error from the choice of fit functions we observed that there is a discrepancy of more than 10% between the experimental value and our result of the isovector root MS charge radius obtained from lattice QCD calculation at the physical point (M π = 135 MeV) on a (10.9 fm) 3 volume [5]. In our lattice QCD simulation, some systematic uncertainties stemming from the chiral extrapolation and finite volume effect are considered to be negligible. Furthermore, excited state contamination in the electric form factor is well controlled and not significant in our simulation, because strong dependence on the time separation between the source and sink operators was not observed. A possible source of this discrepancy could be the effect of finite lattice spacing though it seems too large for O(a 2 ) effect in our non-perturbative improved Wilson quark calculation. Our future calculations performed at the finer lattice spacing will reveal the presence of the systematic error due to the lattice discretization effect. Nevertheless, we are still pursuing other reasons. We are interested in recent development of another approach, called the momentum derivative method, which can directly calculate the slope of the form factor. In comparison to the standard approach, this method should be useful to pin down the source of the current discrepancy between our lattice result and the experimental value.
The momentum derivative method was proposed in Ref. [19], where the slope is determined without assuming fit functions of the form factor. This method employs the moments of the 3-point function in the coordinate space, which can access the derivatives of the form factor with respect to the square of four-momentum transfer q 2 at vanishing q 2 . This method and its variation were applied to the nucleon form factor with the vector current at the pion mass of M π = 0.4 GeV [20], and the physical M π [21], respectively. Another method of the direct derivative calculation using the point splitting vector current was proposed in Ref. [22], and it was applied for the electric, magnetic, and axial-vector form factors at the physical M π in Ref. [23].
In this study we adopt the former method to calculate physical quantities, including the form-factor slopes determined at q 2 = 0 for both the vector and axial-vector channels, in (2+1)-flavor lattice QCD at very close to the physical M π on the (5.5 fm) 3 spatial volume.
The physical quantities obtained from the momentum derivative method are compared with those from the standard analysis for the form factors. Using similar simulation setup as described in our previous works [4,5], possible systematic errors involved in the momentum derivative method are discussed by examining the effect of excited state contaminations with three source-sink separations and also by the finite volume study with two different volumes in our pilot calculations at a heavier pion mass of M π = 0.51 GeV. In this paper we also elucidate the equivalence between the above-mentioned two direct methods through the discussion of an infinitesimal transformation on the correlation functions. This study is regarded as a feasibility study towards more realistic calculation with the PACS10 configurations on the (10.9 fm) 3 volume [5,24,25]. This paper is organized as follows. Section II explains definitions for the nucleon correlation functions and their derivative calculated by moments of the correlation functions used in this study. We also discuss the equality between two types of the direct methods, that were proposed in Ref. [19] and Ref. [23], to calculate the derivative of the form factors in this section. The simulation parameters are described in Sec. III. The results from the derivative of nucleon correlation functions are presented in Sec. IV. Section V is devoted to summary of this study. In two appendices, we firstly describe how the momentum derivative method is associated with a partially quenched approximation and secondly summarize the results obtained from the standard analysis of the form factors.
All dimensionful quantities are expressed in units of the lattice spacing throughout this paper, unless otherwise explicitly specified. A bold-faced variable represents a threedimensional vector.

A. Correlation function with momentum
The exponentially smeared quark operator q S (t, x) with the Coulomb gauge fixing is employed in this study to calculate the nucleon 2-and 3-point functions as where q(t, x) presents a local quark operator, and the color and Dirac indexes are omitted.
A smearing function φ(r) is given in a spatial extent of L as with two parameters A and B. The nucleon 2-point function with the local sink operator is defined as where P + = (1 + γ 4 )/2, and the nucleon operator is given for the proton state by with C = γ 4 γ 2 , the up and down quark operators u, d, and a, b, c being the color indexes.
The smeared source operator N S (t, x) is the same as the local one N L (t, x), but all the quark operators u, d are replaced by the smeared ones defined in Eq. (1). We also calculate the smeared sink 2-point function C SS (t, x), where N L (t, x) is replaced by N S (t, x). The momentum projected 2-point function is then given by with X = L, S and a three-dimensional momentum p. In a large t region, the 2-point function behaves as a single exponential function, where M N and E N (p) are the nucleon mass and energy with the momentum p ≡ |p|, respectively. The overlap of the nucleon operator to the nucleon state is defined by We evaluate the nucleon 3-point functions as where P k is a projection operator for P t = P + and P 5j = P + γ 5 γ j for j = x, y, z, and J O α is an isovector local current operator as J O α = uO α u − dO α d with O α = γ α , γ α γ 5 for the vector (V ) and axial-vector (A) currents, respectively.
The form factors in non-zero momentum transfers are calculated by the momentum projected 3-point function as We evaluate three types of 3-point function, C t Vt (t; p), C 5j V i (t; p), and C 5j A i (t; p) for i = x, y, z to obtain the electric and magnetic form factors, G E (q 2 ), G M (q 2 ), and the axial-vector and induced pseudoscalar form factors, F A (q 2 ), F P (q 2 ).
The asymptotic forms in t ≫ t src and t sink ≫ t, where t src and t sink are defined in Eq. (7), for each 3-point function are given by where the squared momentum transfer is given by The renormalization factors Z V and Z A are defined through the renormalization of the local vector and axial-vector currents on the lattice, respectively. All 3-point functions share a common part of C 2 (t; p), which is similar to the asymptotic form of the 2-point function in Eq. (6), However, C 2 (t; p) is simply eliminated by considering the following ratio [26,27] R k which is constructed from a given 3-point function defined in Eqs. (9)-(11) with appropriate combination of 2-point functions (6). The ratios for each 3-point function give the following asymptotic values in the asymptotic region: which contain the respective form factors.
The MS radius of a form factor G O (q 2 ) is defined by with In the standard way to determine the MS radius, G O (q 2 ) is at first fitted by dipole, quadratic, and z-expansion forms [28,29] given by where the z-expansion makes use of a conformal mapping from q 2 to a new variable z defined as with t cut = 4M 2 π for G E and G M , or with t cut = 9M 2 π for F A , where M π corresponds to the simulated pion mass. Thanks to its rapid convergence of Taylor's series expansion in terms of z, we employ a cubic z-expansion form as a model independent fit as described in our previous work [5]. Using the resulting fit parameter given in each fit, the MS radius can be determined as B. Momentum derivatives of the 2-and 3-point functions As proposed in Ref. [19], the second-order momentum derivative of the 2-point function with respective to p i at the zero momentum point is calculated by where the summation is calculated over −L/2 + 1 ≤ r i ≤ L/2. The superscript (2) in C (2) LS denotes the second-order derivative. In a large t region, a ratio of the derivative function at the zero momentum point to the zero momentum 2-point function becomes where the first term A represents a constant contribution corresponding to the derivative of the amplitude in Eq. (6). It should be noted that the constant A does not contain both the overlap with the local operator, Z L , and its derivative. This is simply because Z X (p) becomes independent of p for the local operator (X = L).
with k = i = j and l = x, y, z. Respective ratios associated with the first to third-order derivatives are defined as with the vector three-point function with the zero momentum, C t Vt (t; 0). The superscript (n) denotes the n-th order derivative.
In this study, we will later determine the magnetic moment, MS charge radius, and MS magnetic radius from the ratios associated with the first-order derivative (25), second-order derivative (26) and third-order derivative (27), respectively, without the q 2 extrapolations of the corresponding form factors toward the zero momentum transfer.
For the axial-vector current, the following three types of the second-order derivatives are considered in the DFF method, C 5j,(2) which are defined with i = j. Just as in the vector cases, the ratios associated with three types of the second-order derivatives are evaluated as below, to directly access the MS axial-vector radius and the value of F P (0)/g A with the axial-vector coupling g A , with l = i, j, ij.
In the asymptotic region of t ≫ t src and t sink ≫ t, the ratios, which are associated with the n-th order derivatives defined in Eqs. (28) and (32), for n ≥ 2, exhibit the following asymptotic behavior where the first term C represents a constant contribution that contains the MS radius of the form factor. The second and third terms can be identified with the contributions from the second-order derivative of − C 2 (t; p)/ C 2 (t; 0) in Eq. (12), whose values coincide with the ones evaluated from R 2 (t) defined in Eq. (24). Thus, the constant C can be isolated using R 2 (t) to subtract the other two contributions from Eq. (33).

C. Equivalence on two definitions of momentum derivatives
In this subsection we intend to discuss an equivalence of the two direct derivative methods.
In the following discussion, variables x, y, z represent four-dimensional coordinates.
The momentum derivative of the quark propagator at the zero momentum is given by where G(x, y) represents the quark propagator. This definition is employed in our calculation and Refs. [19,21,30,31]. Another definition [22] of the momentum derivative is given by means of the point splitting vector current γ j as where G p (x, y) represents the quark propagator calculated through the phase rotation of the gauge link as U j (x) → e ip j U j (x) with the phase associated with the momentum p j . The definition of γ j appearing in Eq. (35) is given by We will discuss an equivalence of the above two definitions as below.
Let us consider an infinitesimal transformation of the quark field ψ( with an arbitrary infinitesimal parameter α(x) depending on x. Requiring the invariance of the expectation value of the quark propagator under this transformation, one finds the following relation: where V µ (z) = ψ(z) γ µ (z)ψ(z) corresponds to the conserved vector current, and ∂ f µ represents a forward difference. When α(x) = εx j with ε = 0, we obtain Here let us consider the proton 2-point function with the local operators having exact isospin symmetry. Performing the transformation of the u quark fields in the 2-point function, one obtains the following relation where V u ν is the conserved vector current for the u quark. The factor of 2 comes from the fact that the proton operator has two u quark fields. A similar relation is obtained through the same transformation of the d quark fields as where V d ν is the conserved vector current for the d quark. As for the isovector current , the corresponding relation is given by a difference of the above two relations as, where there is no contribution from quark disconnected diagrams in the left hand side, since they can be canceled due to exact isospin symmetry.
In our pilot calculations at a heavier pion mass of M π = 0.51 GeV, we verify the left and right equality of the above equation (41), i.e., the equivalence between the two definitions of the momentum derivative. We numerically confirm that they are reasonably consistent with each other through verification of the following two equations as, The numerical verification of Eq. (42) is essentially the same as that of the renormalization factor of Z V = 1. It should be noted that the equality described in Eq. (43) is valid only if both the 2-point and 3-point functions are constructed by the local nucleon operators N L . If the smeared nucleon operators N S are used, the corresponding relation becomes more complicated due to explicit spatial dependence of the smearing function, i.e., φ(r) appearing in Eq. (1). On the other hand, the equality described in Eq. (42) can be applicable even for the smeared nucleon operators N S , since the coordinate in the direction of the derivative is not spatial but temporal. The smearing function φ(r) is clearly independent of the temporal coordinate.
Finally, recall that the smeared nucleon source operator is adopted in our whole calculation. Thus, the quantity, which we actually calculated as defined in Eq. (23), may exhibit a slight difference from the right hand side of Eq. (43) in the asymptotic region. This difference is absorbed only into the difference of the constant A, which is attributed to the difference in the overlaps to the nucleon, Z X (p), defined in Eq. (6), and also their derivatives, between X = L and S. Therefore, it should not matter for extraction of the physical quantities using the DFF method.
Hereafter, we adopt the definition of Eq. (34) for the momentum derivative in our whole calculation, since the calculation cost of the right hand side in Eqs. (42) and (43) is much smaller than that of the left hand side.

III. CALCULATION PARAMETERS
The configuration used in this study was generated with the Iwasaki gauge action [32] and the six stout-smeared Clover quark action near the physical point for intended purpose where the finite-size study was done for light hadron spectroscopy in N f = 2 + 1 lattice QCD at the physical point [24,25]. The lattice size is L 3 × T = 64 3 × 64 corresponding to (5.5 fm) 4 with a lattice cutoff, a −1 = 2.3162(44) GeV [25]. Details of the parameters for the gauge configuration generation are summarized in Ref. [24].
For the measurements for the nucleon correlation functions, the same quark action as in the gauge configuration generation is employed with the hopping parameter κ = 0.126117 for the light quarks and the improved coefficient c SW = 1.11 [33]. The quark propagator is calculated using the exponential smeared source in Eq. method [34][35][36] with the deflated Schwartz Alternative Procedure (SAP) [37] and Generalized Conjugate Residual (GCR) [38] for the measurements as in our previous work [5]. We compute the combination of correlator with high-precision O org and low-precision O approx as where the superscript f, g denotes the transformation under the lattice symmetry G. In our calculation, it is translational symmetry, e.g., changing the position of the source operator, and rotating the temporal direction using the hypercube symmetry of the configuration.  Table I. We also take the average of the forward and backward 3-point functions, and also three 3-point functions with the projector P 5j in all the three spatial directions j = x, y, z to increase statistics.
In our calculation we obtain M π = 0.1382(11) GeV and M N = 0.9499 (27) GeV, which agree with the previous calculation done with the same configuration [24,25]. The result of M N is shown by the solid lines in Fig. 1, which is obtained from a single exponential fit of C LS (t; 0) in the asymptotic region. Although a little finite volume effect of 3 MeV was observed in the pion mass M π by comparing with results obtained on the L = 64 and L = 128 lattice volumes [24], the finite volume effect was not seen in the nucleon mass M N [25]. For the renormalization factors, we adopt Z V = 1/R t Vt (t; 0) with G E (0) = 1, and Z A = 0.9650(68) evaluated by the Schrödinger functional scheme [39]. The physical quantities obtained from the DFF method are compared with the standard analysis of the form factors evaluated with the same configuration, and also the ones from the larger volume calculation of (10.9 fm) 3 in Ref. [5]. The parameters for the larger volume data used in the comparison are summarized in Table II. For a study of finite volume effect for physical quantities obtained from the DFF method, a small test calculation is carried out at a heavier pion mass of M π = 0.51 GeV using the N f = 2 + 1 configurations with a −1 = 2.194 GeV generated in Ref. [40]. The parameters for this study are tabulated in Table III.
with a fit range of t = 10-14, using the linear fixed slope with the measured M N obtained from the standard nucleon 2-point function. This fit result describes the data well in the large t region, though in the small t region we observe a deviation from the linear behavior which indicates the unwanted excited-state contributions appearing in the ratio. In the following analyses, the fit result R 0 2 is utilized to eliminate the common contribution in R k,(n) Oα,(l) (t) of Eq. (33) for n ≥ 2, so as to extract the physical quantities of interest as discussed in Sec. II B. The t dependence of R t,(2) V 4 ,(l) (t) is presented in Fig. 3 together with the data of R 2 (t). The slope of R t,(2) V 4 ,(l) (t) reasonably agrees with the one of R 2 (t) as expected. At a glance, it is observed that R t,(2) V 4 ,(l) (t) has the smaller effect from excited states than the one appearing in R 2 (t), since the t dependence of R t,(2) V 4 ,(l) (t) exhibits almost linear behavior even in a small t region.
Considering the second derivative of C t Vt (t; p), one finds the asymptotic form of R t,(2) based on the asymptotic form of C t Vt (t; p) given in Eq. (9), where we use the following relation with q 2 = 2M N (E N (p 2 ) − M N ) and the condition of q = p.
As discussed in Sec II B, the last two terms can be removed with R 2 (t) or a set of two quantities: R 0 2 obtained from R 2 (t) (through Eq. (45)) and the measured M N obtained from the standard spectroscopy. Figure 4 shows the result of the effective MS charge radius r 2 E determined from in each t using two measured quantities of R 0 2 and M N (denoted as diamond symbols). For comparison, we also plot the result (denoted as circle symbols) given by a naive determination of r 2 E eff as 3(R t,(2) V 4 ,(l) (t) − R 2 (t)) using the raw data of R 2 (t). A little t dependence is observed in the naive subtraction due to the non-negligible excited state contamination in R 2 (t) as explained earlier. Since the former value of r 2 with the solid lines. The value of r 2 E obtained from the above analysis (denoted as the DFF method) is tabulated in Table IV. In the following analysis, the systematic error of quantities obtained from the DFF method will be quoted in the same manner as described above.
The result of r 2 E is compared with the one determined from the fitting of the q 2 dependence of G E (q 2 ) in Fig. 4. The result of r 2 E obtained from the dipole fit defined in Eq. (18) on the data of G E (q 2 ) has the smaller statistical error (inner error bar) than that of the above mentioned DFF method, while the rather large systematic uncertainty (outer error bar) on the fit result of G E (q 2 ) is inevitable due to the choice of the fit form. A systematic error is estimated by maximum discrepancy of the results obtained with different fit forms, namely the quadratic and z-expansion forms as defined in Eqs. (19) and (20). Therefore, the results obtained from both the standard and DFF methods, are mutually consistent within the total errors whose sizes are comparable. The two results obtained in this study also reasonably agree with the previous result obtained from the standard method in our calculation performed on the larger volume (L = 128) at the physical point [5].
It is worth mentioning that, compared to the values obtained from the standard calculation with the form factor G E (q 2 ), the result of r 2 E in the DFF method is likely closer to both experimental values from the electron-proton scattering [2], 0.882(11) fm 2 , and muonic hydrogen spectroscopy [41], 0.823(2) fm 2 , while there remains a discrepancy of more than 10%. It might not be attributed to excited state contributions. The reason is that the data given with t sep = 12 completely agrees with those with t sep = 14 as shown in Fig. 5 in the DFF method. The data given with t sep = 16 is also included in the figure, though it is too difficult to determine that there is an obvious dependency with respect to t sep because of its large statistical errors. In Ref. [23], a large t sep dependence was reported for the calculation of r 2 E using another derivative method. We, however, consider that it could be caused by their choice of smearing parameters for the quark operators, since our smearing parameters are highly optimized to eliminate the excited state contributions in the nucleon 2-point and 3-point functions.
The finite volume effect is other possible source of the systematic errors in the DFF method. The large effect on the quantity of r 2 E is not observed in our standard analyses with the form factor G E (q 2 ) obtained from the L = 64 and L = 128 lattice data shown in Appendix B. A significant effect, however, is reported in a previous study in momentum derivative calculations of meson correlation functions [30], and also found in our pilot study of r 2 E eff at a heavier pion mass of M π = 0.51 GeV. Figure 6 shows that there is a large discrepancy of r 2 E eff on volumes of the spacial extent of 2.9 and 5.8 fm at M π = 0.51 GeV, while the two results determined from the form factor in the standard way are highly consistent with each other and also the larger volume result from the DFF method.
Thus, for a precise determination of r 2 E , it is an important future task to investigate the systematic uncertainty in the DFF method due to the finite size effect near the physical point, using our large volume configuration of L = 128. Finally, we remark on an attempt of the improved analysis in the DDF method for the case of meson form factors [30,31] in order to reduce the finite volume effect. This improvement can be also extended to the nucleon form factors. results, which are given by the average of results from t sep = 12, 14 and 16 within the standard analysis in our previous calculation [5], are also tabulated.   V i ,(k) (t), using its asymptotic form obtained with Eq. (25). The effective magnetic moment is defined by with k = i = j, whose form can be read off by considering the derivative of the asymptotic form of C 5j V i (t; p) in Eq. (10). It is here worth pointing out that the ratio R   (denoted as the DFF method) is tabulated in Table IV The result in the DFF method is in good agreement with the two fit results of G M (q 2 ) obtained in this study (L = 64) and also from the larger volume calculation (L = 128) [5].
However, the significant reduction of the systematic error in the DFF method may expose an underestimation of µ in comparison with the experimental value, µ = 4.70589, in PDG20 [42]. In order to investigate possible other systematic errors in the DFF method, the direct comparison of the data with t sep = 12, 14, and 16 is first presented in Fig. 8.
The three data statistically agree with each other in the middle t region. We thus consider that the present data set shows no significant uncertainty associated with the excited state contamination in the data of t sep = 14, though the statistical error of t sep = 16 is larger than those for the other two data.
We rather consider that the discrepancy from the experiment might be caused again by a finite volume effect. As shown in Fig. 9, in our pilot study of the DFF method at M π = 0.51 GeV, we observe a large finite volume effect of more than 10% in comparison with the data on the (2.9 fm) 3 and (5.   V i ,(kl) (t) yields the following asymptotic form: with l = k where l and k are the directions of the derivative defined in Eq. (27). When l = k is set in Eq. (27), the right hand side of the above equation is multiplied by a factor of three. In our analysis all the l data are averaged.
Using the asymptotic form of Eq. (50), the effective MS magnetic radius in the DFF method is determined by with the first derivative R ) for the same reason as in the case of r 2 E eff . In the middle t region the two results are mutually consistent, so that we determine the value of r 2 M by a constant fit of r 2 M eff in the region of t = 5-9, which is plotted by the solid lines in Fig. 10, and tabulated in Table IV.
Although the statistical error of r 2 M obtained from the fit of G M (q 2 ) on the L = 64 lattice volume is smaller than that of the DFF method, much larger systematic error regarding the choice of fit function makes their total accuracy worse than the DFF result.
While the result of r 2 M given in the DFF method has smaller systematic error, its central value is much smaller than the experimental value, r 2 M = 0.733(32) fm 2 [42]. It might be caused by a finite volume effect in the DFF calculation, since this quantity is indeed considered to be more sensitive to finite volume effects compared to other quantities discussed earlier. In the DFF method, the quantity of r 2 M requires the third derivative of the 3-point function i.e., the third moment of the 3-point function in coordinate space. Therefore, the larger spatial extent is naturally required for its higher moment compared to the other quantities considered in the previous subsections. Figure 11 shows that in our pilot calculation of r 2 M eff at M π = 0.51 GeV, more significant finite volume effect is observed than the case of r 2 E eff in Fig. 6: the data of r 2 M eff on the smaller volume becomes negative in a large t region in contrast to that of r 2 E eff .
Another possible source of systematic errors comes from the excited state contamination.
Comparing the three data sets using t sep = 12, 14, and 16 in Fig. 12, the t sep dependence is not clearly seen due to their large statistical fluctuations. In a future work, it is important to investigate the above two systematic errors by calculations with the larger volume (L = 128) and the large variation of t sep . V i ,(kl) (t). The cyan band represents the experimental result [42].

E. MS axial radius r 2
A For the MS axial radius r 2 A , the same analysis as in the case of r 2 E with R t,(2) A j ,(i) (t), where i = j in the subscript expresses the direction of the derivative defined in Eq. (29). The asymptotic form of R 5j, (2) A j ,(i) (t) is given by The effective MS axial-vector radius is defined by and its value is plotted as a function of t in Fig. 13. The data is compared with the one determined from a naive subtraction with the raw data of R 2 (t) as 3(R 5j, (2) A j ,(i) (t)−R 2 (t)). Both estimations for r 2 A eff exhibit a reasonably flat behavior in the middle t region, respectively.
We determine the value of r 2 A by a constant fit of the data of r 2 A eff in the region of t = 5-9. The result obtained by the DFF method is tabulated in Table IV. As shown in Fig. 13, the DFF result is fairly consistent with the two dipole-fit results of F A (q 2 ) obtained in this study (L = 64) and also from the larger volume calculation (L = 128) [5]. Again, the total accuracy of the dipole fit results is slightly worse than the DFF result that can avoid any model dependence.
As in the case of the dipole fit results, the result of the DFF method is little smaller than the experiment [43,44], 0.449(13) fm 2 . It can be attributed to excited state contamination.
It is simply because a visible difference between the data from t sep = 12 and 14 is seen in the flat region of r 2 A eff as a function of t as shown in Fig. 14. Furthermore, the data of t sep = 16 statistically agrees with the experiment, though the statistical uncertainties are not small enough to make a firm conclusion. In addition, recently, it was reported that a large effect due to the excited state contamination exists in determination of r 2 A even from the fitting of the q 2 dependence of F A (q 2 ) [13,17]. Since we have only a few variations of t sep , we will need to verify whether or not the discrepancy from the experiment can be explained by the systematic uncertainties from the excited state contamination using the data with the larger variation of t sep . It is worth pointing out that the systematic error associated with the finite volume effect in r 2 A eff might be smaller than the other quantities obtained from the DFF method discussed earlier. This is expected from our pilot calculation at M π = 0.51 GeV as shown in Fig. 15. The difference between the data on the larger and smaller volumes is about 10% in the middle t region, which it is much smaller than those of other quantities as shown in Figs. 6, 9, and 11.
The size of the finite volume effect is supposed to depend on the MS radius of the target form factor. If the MS radius is small, in other words, the form factor has a broad shape in the q 2 space, its moment in the coordinate space is less sensitive to the finite volume. This is because the narrow spatial distribution in the coordinate space is given by the inverse A j ,(i) (t). The cyan band represents the experimental result [43,44].  in Eq. (31). Their asymptotic forms are given by where i = j. Using them, the effective value of 2M N F P (0)/g A can be defined in two ways as where the second term of R 5j, (2) A j ,(i) appearing in Eq. (57) is the one used for the determination of r 2 A eff as discussed in the last subsection. Figure 16 shows that the two different estimations of 2M N F eff P (0)/g A with t sep = 14 (denoted with filled symbols) provide inconsistent results: the result obtained from Eq. (56) is much smaller than that of Eq. (57). The same behavior is also seen in the data from the smaller and larger source-sink separations for t sep = 12 and 16. Note that the expected value of 2M N F P (0)/g A is much larger than these two observed values according to the following reasons: (1) 2M N F P (q 2 1 )/g A ∼ 40 is observed even at the lowest non-zero q 2 = q 2 1 in the standard method with t sep = 14. (2) The q 2 dependence of F P (q 2 ) is expected to rapidly increase in the limit of q 2 → 0 as shown in Appendix B.
This difference between the two estimations of 2M N F eff P (0)/g A can be attributed to a finite volume effect. This is simply because a similar trend is observed in our pilot calculation at M π = 0.51 GeV, where the discrepancy between the two results becomes resolved in a larger volume calculation as shown in Fig. 17. A large systematic effect stemming from the finite volume is easily understood in this case, since the induced pseudoscalar form factor in the coordinate space, which corresponds to the one given by the inverse Fourier transform of F P (q 2 ), has a very broad structure. This is naively expected from the fact that F P (q 2 ) has a sharp peak near the origin at the physical M π corresponding to the large contribution of the pion pole in the pion pole dominance model. The moment of such a broad function in the coordinate space could not avoid the strong dependence of the finite volume.
In addition to the finite volume effect, in our previous studies, we observed other problem that the lattice data of F P (q 2 ) differs from the pion pole dominance model [4,46], which can be explained by large excited state contamination [5]. Similar discussions were reported in Refs. [13-15, 17, 47]. In order to fully resolve the problems, more comprehensive investigations are necessary in this particular quantity. A i ,(ij) (t) (circle symbols) and R 5j, (2) A j ,(j) (t) (triangle symbols). The black, red and blue symbols represent the data obtained with t sep = 12, 14 and 16.

V. SUMMARY
We have calculated the MS radii and magnetic moment for the isovector nucleon form factors using the DFF method, which is a direct calculation method of the derivative of the form factor proposed in Ref. [19], in the N f = 2 + 1 QCD near the physical point on the (5.5 fm) 3 volume. We have also discussed an equivalence of the method used in this study to another derivative method [22].
The results from the DFF method near the physical point are compared with the ones from the standard form factor calculation on the same volume (L = 64) and also on the larger volume (L = 128) in our previous work [5]. For r 2 E , µ, r 2 M , and r 2 A , the statistical uncertainties of the results in the DFF method are relatively larger than those obtained from the fitting of the q 2 dependence of the corresponding form factors. However, the DFF method can avoid the systematic error associated with the model dependence of the fit form.
Such a systematic error is dominant over the statistical error in the standard method, and then it makes the total accuracy worse than the DFF results. We have also confirmed that the results from the DFF method are in good agreement with the ones from the standard analysis with the form factors within the combined errors of the statistical and systematic uncertainties in both methods, except for the quantity of F P (0). Two ways to determine F P (0) in the DFF method provide inconsistent results. A similar discrepancy is observed in our pilot calculation at a heavier pion mass of M π = 0.51 GeV on the smaller volume of (2.9 fm) 3 , while it becomes resolved on a larger volume of (5.8 fm) 3 .
We thus have considered that the discrepancy comes from a finite volume effect, and the significant finite volume effect can be related to a steep behavior of F P (q 2 ) near the origin, in other words, its broad shape in the coordinate space.
In our pilot calculation at a heavier pion mass of M π = 0.51 GeV, we have also found that the DFF method is more sensitive to finite volume effect than the standard method as reported in the previous works [30,31]. Therefore, an undetermined systematic error stemming from the finite volume effect might exist in the DFF results from our numerical simulations on the (5.5 fm) 3 volume near the physical point, although they agree with the results obtained by the standard method. One of the important future works is a comprehensive study of systematic errors in the DFF method, including the one from excited state contamination, at the physical point using the (10.9 fm) 3 volume as is done in our previous study of the form factors in the standard method [5]. In this future direction, some improvements in the DFF method are needed to reduce the finite volume effect in the analysis level. Such an improved analysis in the DFF method was proposed for the case of meson form factors [30,31], so that it can be extended and applied to the nucleon form factors.

ACKNOWLEDGMENTS
We thank members of the PACS collaboration for useful discussions. Numerical calcula- Under a partially quenched approximation of the quark field subjected to uniform phase rotation, the condition of det D p j = 1 is imposed in Eqs. (A2) and (A3). It thus ends up that the second and third terms in Eq. (A4) disappear, since both terms arise from the derivative of det D p j . Therefore, the momentum derivative of the quark propagator is expressed only by the first term in Eq. (A4) under the partially quenched approximation. The results for the renormalized isovector G E (q 2 ), G M (q 2 ), F A (q 2 ), and F P (q 2 ) with t sep = 12, 14, and 16 are shown in Figs. 18-21 together with those obtained in the previous calculation on a 128 4 lattice [5] for comparison. For F A (q 2 ) and F P (q 2 ), the error of Z A = 0.9650(68) [39] is included in their errors. The values of each form factor obtained with t sep = 12, 14, and 16 are tabulated in Table VI, VII, and VIII, respectively. Our results with t sep = 12, 14, and 16 on the L = 64 lattice volume are consistent with each other, and also statistically agree with the data on the L = 128 lattice volume [5] in all the form factors, except for F P (q 2 ). There is a clear discrepancy between the data of F P (q 2 ) obtained with t sep = 12 and 14 on the L = 64 lattice volume. It is considered to be caused by significant effect from excited state contamination as reported in Refs. [5, 13-15, 17, 47].
It is worth remarking that the value of g A obtained with t sep = 12 is slightly smaller than the one obtained with t sep = 14. Since the data with t sep = 16 has a much larger statistical error than the two data, the t sep dependence is unclear as shown in Fig. 22. The figure also shows that any appreciable t sep dependence was not observed in the range of t sep = 10-16 in the larger volume calculation [5], and they are statistically consistent with all the data on the L = 64 lattice. Therefore, it is not clear whether this slight difference between two results from t sep = 12 and 14 is just a statistical fluctuation or related to the systematic errors, e.g., the finite volume effect and the excited state contamination. To clarify this point, a further systematic study using a large set of different t sep with the statistical error as small as the two data is needed.
The fit results with dipole, quadratic, and z-expansion functions in Eqs. (18)- (20) for G E (q 2 ), G M (q 2 ), and F A (q 2 )/g A are summarized in Tables IX-XI. The value of χ 2 /dof is obtained from a correlated fit. The maximum q 2 in each fit denoted by q 2 cut in the tables is chosen to obtain an acceptable χ 2 /dof. In all the fits for G E (q 2 ) and F A (q 2 )/g A , we use the condition of G E (0) = F A (0)/g A = 1 imposed in the respective fit functions.
The central values and their statistical errors presented in Table IV are [5] given after taking average of three data sets calculated with t sep = 12, 14 and 16 is also plotted by the asterisk symbol. The red curve represents Kelly's parametrization of the experiment data [48].   FIG. 20. Same as Fig. 18 for F A (q 2 ). The red curve is given by a dipole form with the dipole mass [43,44] and g A [42]. FIG. 21. Same as Fig. 18 for 2M N F P (q 2 ). The red curve is given by the pion-pole dominance model (PPD) 2M N F PPD P (q 2 ) = 4M 2 N F A (q 2 )/(M 2 π + q 2 ) with a dipole form of F A (q 2 ) using the dipole mass [43,44] and g A [42]. The experimental result of the pion-electroproduction [49] is also plotted by the diamond symbol. The data of L = 128 are slightly shifted to the negative x direction for clarity. The cyan band represents the experimental value [42].    (20). The fit range is q 2 = q 2 1 -q 2 cut . In all the fit we employ G E (0) = 1.  (20). The fit range is q 2 = q 2 1 -q 2 cut . In all the fit we employ F A (0)/g A = 1.