Nucleon form factors on a large volume lattice near the physical point in 2+1 flavor QCD

We present results for the isovector nucleon form factors measured on a $96^4$ lattice at almost the physical pion mass with a lattice spacing of 0.085 fm in 2+1 flavor QCD. The configurations are generated with the stout-smeared $O(a)$-improved Wilson quark action and the Iwasaki gauge action at $\beta$=1.82. The pion mass at the simulation point is about 146 MeV. A large spatial volume of $(8.1~{\rm fm})^3$ allows us to investigate the form factors in the small momentum transfer region. We determine the isovector electric radius and magnetic moment from nucleon electric ($G_E$) and magnetic ($G_M$) form factors as well as the axial-vector coupling $g_A$. We also report on the results of the axial-vector ($F_A$), induced pseudoscalar ($F_P$) and pseudoscalar ($G_P$) form factors in order to verify the axial Ward-Takahashi identity in terms of the nucleon matrix elements, which may be called the generalized Goldberger-Treiman relation.


I. INTRODUCTION
The nucleon vector and axial elastic form factors are good probes to investigate the internal structure of the nucleon [1]. Although great theoretical and experimental efforts have been devoted to improving our knowledge of the nucleon structure, there are several unsolved problems associated with fundamental properties of the proton and neutron. The proton radius puzzle, where highprecision measurements of the proton's electric charge radius from the muonic hydrogen Lamb shift [2,3] disagree with well established results of both electron-proton scattering and hydrogen spectroscopy [4], is currently one of the most intriguing problems in this field [5]. The neutron lifetime puzzle, where the discrepancy between the results of beam experiments and storage experiments remains unsolved, is another open question that deserves further investigation in terms of the nucleon axial-vector coupling g A [6].
Much effort has been devoted to calculate the nucleon form factors with lattice QCD since 1980's. Unfortunately, satisfactory results, that are consistent with the experiments, have not yet been achieved for the nucleon structure in the previous lattice QCD simulations. The current situation is that we are still struggling for reproducing well-known experimental results, e.g., the axialvector coupling and the electric charge radius. This means that we have not yet achieved our final goal that a single-nucleon state is properly generated in lattice QCD calculation. The discrepancy between existing lattice calculations and experimental values could be attributed to unresolved systematic errors in numerical simulations. A major part of the systematic errors should be due to the fact that the simulated quark mass used in simulations is heavier than the physical one. However, the particular quantities like the axial-vector coupling [59] and the electric charge radius still show some discrepancies with respect to the experimental values in the previous lattice QCD simulations at almost physical point [8][9][10].
In this paper, we aim to calculate the nucleon form factors on the very large spatial volume in realistic lattice QCD, where the light quark masses are down to their physical values. For this purpose, we use the 2+1 flavor QCD gauge configurations generated on an (8.1 fm) 4 lattice at almost physical point (the pion mass m π at the simulation point is about 146 MeV) by the PACS Collaboration [11]. There are three reasons for paying special attention to its large spatial volume of (8.1 fm) 3 .
First of all, the large finite volume effect, represented by the m π L scaling (L denotes the spatial extent of the lattice volume), on measurement of the axial-vector coupling g A was reported in Ref. [12,13]. According to their conclusion, one can estimate that spatial sizes of 7.7−9.4 fm (m π L ≈ 5.7 − 6.9) for m π ≈ 146 MeV is necessary for keeping the finite volume effect on g A at or below 1%.
Secondly, thanks to the large spatial volume, we can access to the small momentum transfer region up to 152 MeV. The momentum squared (q 2 ) dependence of the nucleon form factors can be carefully examined on the very large spatial volume. Indeed, the values of a given form factor at very low q 2 help to more accurately determine the slope of the form factor at q 2 = 0, which is associated with the root-mean-square (RMS) radius R. In this context, uncertainties in the determination of R are sensitive to how large is the spatial extent L in phys-ical units. For the case of the electric form factor (G E ), R corresponds to the charge radius.
Finally we revisit the recent claim in the literature that the charge density of the proton, which has the shape of an exponential in a classical argument, is widely distributed in space [5]. The q 2 dependence of the Sachs form factors G E (q 2 ) and G M (q 2 ) are roughly described by the dipole shape: G D (q 2 ) = Λ 2 /(q 2 + Λ 2 ) 2 with a single parameter Λ called as the dipole mass parameter around 0.84 GeV [1]. The dipole form assumes that the spatial charge distribution is falling exponentially at large r as ρ(r) ∝ e −rΛ [1]. For a sake of simplicity, we adopt the spatial charge distribution as ρ(r) = ρ 0 e −rΛ , where ρ 0 is determined by the normalization condition. The RMS radius R is defined in terms of the charge density as which gives a relation of R = √ 12/Λ. We next introduce the truncated RMS radius R(r cut ), where the integral (1) is stopped at r = r cut , and then calculate the following ratio as a function of r cut [5]: where X = r cut Λ = √ 12r cut /R. To get more than 98% of R, the value of r cut must be greater than 2.75R, which is remarkably large value [5]. If we take it seriously, this intuitive argument suggests that the spatial extent of 2.75R × 4 [60], which is roughly 10 fm, is required for precise determination of the proton charge radius within a few % accuracy on a periodic hyper-cubic lattice.
This paper is organized as follows: In Sec. II, we present a brief introduction of general features of nucleon form factors. In Sec. III, we first summarize simulation parameters in 2+1 flavor ensembles generated by the PACS Collaboration [11] and then give some basic results from the nucleon two-point function. We also describe the lattice method (standard ratio method) for calculating the isovector form factors of the nucleon from relevant three-point correlation functions. The results of our lattice calculations for the nucleon form factors are presented in Sec. IV, which is divided into three major subsections. We first determine four kinds of nucleon couplings: the vector coupling g V , the axial-vector coupling g A , the scalar coupling g S and the tensor coupling g T , which are directly accessible from the three-point correlations function at zero momentum transfer in Sec. IV A. Section IV B presents the results of the isovector electric (G E ) and magnetic (G M ) form factors, including detailed study of the momentum transfer dependence and determination of the isovector electric and magnetic RMS radii and also the isovector magnetic moment. The last subsection (Sec. IV C) is devoted to a discussion of the results of three form factors obtained in the axial-vector and pseudoscalar channels, which are related by the axial Ward-Takahashi identity in terms of the nucleon matrix elements. Finally, we close with a brief summary and our conclusions in Sec. V.

II. GENERAL FEATURES OF NUCLEON FORM FACTORS
In general, four form factors appear in the nucleon matrix elements of the weak current. Here, for example, we consider the matrix element for neutron beta decay. In this case, the vector and axial-vector currents are given by , and then the general form of the matrix element for neutron beta decay is expressed by both the vector and axialvector transitions: where q = P n − P p is the momentum transfer between the neutron (n) and proton (p). The vector (F V ) and induced tensor (F T ) form factors are introduced for the vector matrix element: and also the axial vector (F A ) and induced pseudoscalar (F P ) form factors for the axial-vector matrix element: These matrix elements are here given in the Euclidean metric convention (we have defined σ αβ = 1 2i [γ α , γ β ]) [61]. Thus, q 2 denoted in this paper, which stands for Euclidean four-momentum squared, corresponds to the spacelike momentum squared as q 2 M = −q 2 < 0 in Minkowski space.
The vector part of weak processes is related to the nucleon's electromagnetic matrix element through an isospin rotation if heavy flavor contributions are ignored under the exact isospin symmetry. A simple exercise in SU (2) Lie algebra leads to the following relation between the vector part of the weak matrix elements of neutron beta decay and the difference of proton and neutron electromagnetic matrix elements where j em α = 2 3ū γ α u − 1 3d γ α d. Therefore, this relation gives a connection between the weak vector and induced tensor form factors and the isovector part of electromagnetic nucleon form factors where M N denotes the nucleon mass, which is defined as the average of neutron and proton masses. given by the isovector combination of the Dirac (Pauli) form factors of the proton and neutron as where individual form factors F N 1,2 (N = p, n) are defined by Experimental data from elastic electron-nucleon scattering is usually described in terms of the electric G E (q 2 ) and magnetic G M (q 2 ) Sachs form factors, which are related to the Dirac and Pauli form factors: Their normalizations at q 2 = 0 are given by the proton's (neutron) electric charge and magnetic moment [62]. Therefore, one finds represents the isovector combination of the electric (magnetic) form factors of the proton and neutron.
The q 2 dependence of the Sachs form factors G N E (q 2 ) and G N M (q 2 ) are experimentally known that the standard dipole parameterization G D (q 2 ) = Λ 2 /(Λ 2 + q 2 ) 2 with Λ 2 = 0.71 [(GeV) 2 ] (or Λ = 0.84 [GeV]) describes well the magnetic form factors of both the proton and neutron and also the electric form factor of the proton, at least, in the low q 2 region. In general, if there is no singularity around q 2 = 0 for a given form factor G(q 2 ), the normalized form factor can be expanded in powers of q 2 .
G(q 2 ) = G(0) 1 − 1 6 r 2 q 2 + 1 120 where the first coefficient determines the mean squared radius r 2 , which is a typical size in the coordinate space. For the dipole form, the root-mean-square (RMS) radius R is given as R = r 2 = √ 12 Λ by the dipole mass parameter Λ.
For the axial-vector part of weak processes, the axialvector form factor at zero momentum transfer, namely, the axial-vector coupling g A = F A (0), is precisely determined by measurements of the beta asymmetry in neutron decay. The value of g A = 1.2723 (23) is quoted in the 2016 PDG [4]. The reason why g A deviates from the corresponding vector coupling g V = F V (0) = 1 is that the axial-vector current is strongly affected by the spontaneous chiral symmetry breaking in the strong interaction [14,15]. In this sense, this particular quantity allows us to perform a precision test of lattice QCD in the baryon sector.
The q 2 dependence of the axial-vector form factor F A (q 2 ) has also been studied in experiments, where the dipole form F A (q 2 ) = F A (0)/(q 2 + M 2 A ) 2 is a good description for low and moderate momentum transfer [16,17]. Recently, the axial mass parameter M A is much paid attention to by neutrino oscillation studies [18,19].
Although the induced pseudoscalar form factor F P (q 2 ) is less well known experimentally [20,21], it is theoretically known that two form factors F A (q 2 ) and F P (q 2 ) in the axial-vector channel are not fully independent. It is because the axial Ward-Takahashi identity: ∂ α A + α (x) = 2mP + (x) leads to the generalized Goldberger-Treiman (GT) relation among three form factors [22,23]: wherem = m u = m d is a degenerate up and down quark mass and the pseudoscalar nucleon form factor G P (q 2 ) is defined in the pseudoscalar nucleon matrix element with a local pseudoscalar density, P + (x) ≡ū(x)γ 5 d(x). Therefore, the q 2 dependences of three form factors, F A (q 2 ), F P (q 2 ) and G P (q 2 ) are constrained by Eq. (17) as a consequence of the axial Ward-Takahashi identity. Therefore, three form factors, F A (q 2 ), F P (q 2 ) and G P (q 2 ) are very important to verify the axial Ward-Takahashi identity in terms of the nucleon matrix elements.
The latest lattice calculations of the nucleon structure have been carried out with increasing accuracy [8-10, 12, 13, 24-27]. It seems that there is still a gap between experimentally known values and the lattice results, especially for the electric-magnetic nucleon radii and the magnetic moment. Our preliminary results computed with almost physical pion mass on very large volume shows agreement with experimental results [28][29][30].
In this paper, we present an update of our previous studies [28][29][30], including the three the axial-vector, induced pseudoscalar and pseudoscalar form factors.

III. SIMULATION DETAILS
We use the 2+1 flavor QCD gauge configurations generated with the 6-APE stout-smeared O(a)-improved Wilson-clover quark action and the Iwasaki gauge action [31] on a lattice L 3 × T = 96 3 × 96 at β = 1.82, which corresponds to a lattice cutoff of a −1 ≈ 2.3 GeV (a ≈ 0.085 fm) [11]. Periodic boundary conditions are used for the gauge and quark fields in all four directions. The stout smearing parameter is set to ρ = 0.1 [32]. The improvement coefficient, c SW = 1.11, is determined nonperturbatively by the Schrödinger functional (SF) scheme [33]. The improved factor c A for the axial-vector current becomes very small at the nonperturbative value of c SW and is consistent with zero within the statistical error [33]. Therefore, we do not consider O(a) improvement of quark bilinear current. The hopping parameters for the light sea quarks (κ ud , κ s ) = (0.126117, 0.124790) are chosen to be near the physical point. For the first time, a simulated pion mass reaches down to m π ≈ 146 MeV on a large spatial volume of (8.1 fm) 3 in 2+1 flavor QCD.
The degenerated up-down quarks are simulated with the DDHMC algorithm [34,35] using the even-odd preconditioning and the twofold mass preconditioning [36,37]. The strange quark is simulated with the UVPHMC algorithm [38][39][40][41][42][43] where the action is UV-filtered [44,45] after the even-odd preconditioning without domain decomposition. The total number of gauge configurations reaches 200 which corresponds to 2000 trajectories after thermalization. Each measurement is separated by 10 trajectories. The results for the hadron spectrum and other physical quantities are already presented in Ref. [11]. We use the jackknife method with a bin size of 50 trajectories for estimating the statistical errors. Our preliminary results of the nucleon vector form factors with less number of measurements were first reported in Ref. [28,29].

A. Nucleon two-point functions
Let us first examine the nucleon rest mass and the dispersion relation, which are obtained from the nucleon two-point functions. In order to compute nucleon energies or matrix elements, we define the nucleon (proton) operator as where the superscript T denotes a transposition and C is the charge conjugation matrix defined as C = γ 4 γ 2 . The indices abc and ud label color and flavor, respectively. The subscript X of the nucleon operator specifies what type of the smearing for the quark propagators. In this study, we use two types of smearing function φ: the local function as φ(x i − x) = δ(x i − x) (denoted as X = L) and the exponential smearing function: with A = 1.2 and B = 0.11 (denoted as X = S). For a simplicity, x 1 = x 2 = x 3 is chosen. We then construct two types of the two-point function with the exponential smearing source (the source-time location denoted as t src ) as where X = L (local) or S (smear) stands for a type of smearing at the sink operator (the sink-time location denoted as t sink ). A projection operator P + = 1+γ4 2 can eliminate contributions from the opposite-parity state for |p| = 0 [46,47]. In this study, for nonzero spatial momentum, we use the nine lowest momenta: p = 2π/L × n with a vector of integers n ∈ Z 3 . All choices of n are listed in Table II.

B. Nucleon spectra and dispersion relation
In Fig. 1 we plot the nucleon effective mass with |p| = 0 for two cases: Smear-local denotes the nucleon two-point function with the smeared source and the local sink operators and smear-smear for the smeared source and the smeared sink ones. We observe that the smear-local case shows good plateau for t/a ≥ 6. On the other hand, the signal becomes noisier for the smear-smear case. We also measure the nucleon energies E N (p) from the smear-local case of the nucleon two-point function with the choice of 9 cases of nonzero spatial momenta. All fit results of E N (p) obtained from the single exponential fit are tabulated in Table III. Figure 2 shows a check of the dispersion relation for the nucleon. The vertical axis shows the momentum squared defined through the relativistic continuum dispersion relation as p 2 con = E 2 N (p)−M 2 N , while the horizontal axis is the momentum squared given by the lattice momentum p 2 lat = (2π/L) 2 × |n| 2 in lattice units. As can be seen in Fig. 2, the relativistic continuum dispersion relation is well satisfied up to |n| 2 = 9.

C. Three-point correlation functions for nucleon form factors
The nucleon form factors are extracted from the threepoint correlation functions consisting of the nucleon source and sink operators with a given local current J O α (O = S, P, V, A and T ) located at the time-slice t: where P k denotes an appropriate projection operator to extract the form factors and q = p − p ′ represents the three-dimensional momentum transfer. The local current is given by We then calculate the following ratio constructed from the three-point correlation function C P O,α with nucleon two-point functions C XS : which is a function of the current operator insertion time t at the given values of momenta p and p ′ for the initial and final states of the nucleon.
We consider only the rest frame of the final state with p ′ = 0, which leads to q 2 = 2M N (E N (q) − M N ) for the squared four-momentum transfer. Nucleon energy E N (q) simply abbreviated as E N , hereafter. In this kinematics, R O,α (t, p ′ , p) is represented by a simple nota-tion R O,α (t, q). We calculate only the connected diagrams to concentrate on the isovector part of the nucleon form factors. The current insertion points t are moved between the nucleon source and sink operators, both of which are exponentially smeared in the Coulomb gauge, separated by 15 time slices. In the physical units, the source-sink separation of t sep /a = 15 (denoted as t sep = t sink − t src ) corresponds to 1.26 fm. In previous calculations [9,10,13,24,25,27], the maximum values of the source-sink separation are typically as large as 1.3-1.4 fm, where the excited-state effect is marginally insignificant.
The three-point correlation functions are calculated by the sequential source method with a fixed source location [48,49]. This method requires the sequential quark propagator for each choice of a projection operator P regardless of the types of the current J O α . To minimize the numerical cost, we consider two types of the projection operator P t = P + γ 4 and P 5z = P + γ 5 γ 3 in this study.
The ratio (22) with appropriate combinations of the projection operator P k (k = t, 5z) and the α component of the current j O α gives the following asymptotic values [23]: and for the vector current (O = V ) in the limit when the Euclidean time separation between all operators is large, t sink ≫ t ≫ t src with fixed t src and t sink . Similarly, we get for the axial-vector current (O = A), and for the pseudoscalar (O = P ), respectively. Here, we recall that the lattice operators receive finite renormalizations relative to their continuum counterparts in general. Thus, all above form factors G v E , G v M , F A , F P and G P are not renormalized yet. Their renormalized values require the renormalization factor Z O (O = V, A, P ), which is defined through the renormalization of the quark bilinear currents The renormalization factors Z V , Z A and Z P are independently obtained by the Schrödinger functional (SF) scheme at the vanishing quark mass defined by the partially conserved axial-vector current (PCAC) relation [50]. The isovector electric G v E and magnetic G v M form factors are simply abbreviated as G E and G M , hereafter.

IV. RESULTS OF NUCLEON FORM FACTORS
All the results presented in this work are obtained with 200 configurations. To reduce the statistical uncertainties, we perform 64 measurements of the two-point and three-point correlation functions on each configuration. For 64 measurements, we use 8 different source locations, 4 choices for a temporal direction on a 96 4 lattice, and 2 choices of both forward and backward directions in time.
In the analysis, all 64 sets of three-point correlation functions and nucleon two-point function are folded together to create the single-correlation functions, respectively. It can reduce possible autocorrelation among measurements. The statistical errors are evaluated by the jackknife analysis with a bin size of 5 configurations.
We employ 9 cases of spacial momentum transfer q = π/48 × n with |n| 2 ≤ 9 as described in Table II. The minimum momentum transfer is about 152 MeV, which is so small as the pion mass. A difference between the results for Q8 and Q9 cases could probe the possible lattice discretization error. For zero momentum transfer q = 0, the ratio (22) becomes which vanishes unless Γ O = 1(S), γ 4 (V ), γ i γ 5 (A), and σ ij (T ) with i, j = 1, 2, 3 [49]. The non-vanishing ratio gives an asymptotic plateau corresponding to the bare value of the coupling g O relevant for the O channel in the middle region between the source and sink points, when the condition t sink ≫ t ≫ t src is satisfied. With our choice of the projection operators k = t, 5z, we can determine 4 different couplings: the vector coupling g V , the axial-vector coupling g A , the scalar coupling g S and the tensor coupling g T , from the following four ratios, whose values are not yet renormalized. In Fig. 3, we plot the above four ratios as a function of the current insertion time slice t for the vector (left upper panel), axial-vector (right upper panel), scalar (left lower panel) and tensor (right lower panel) channels. The best plateau is clearly observed in the vector channel since the vector coupling g V corresponds to the conserved charge associated to the isospin symmetry. The exact symmetry would suppress the statistical fluctuations and hinder parts of the excited-state contamination. Therefore, a very long plateau indeed appears in the ratio of R t V,4 (t, 0). It is also worth recalling that the time-reversal average was implemented in our averaging procedure on each configuration with multiple measurements, which include both forward and backward propagations in time, as described previously.
In the vector channel (left upper panel), t-dependence of R t V,4 (t, 0) is symmetric between the source (t/a = 0) and sink (t/a = 15) points. Although this symmetric behavior with respect to t is hindered by larger statistical fluctuations in both the scalar (left lower panel) and the axial-vector (right upper panel) channels, its behavior is clearly seen in the tensor channel (right lower panel) with relatively small errors. In the case of the tensor, good plateau is observed only in the middle region between the source and sink points. Therefore, we choose the range of 6 ≤ t/a ≤ 9, where an asymptotic plateau behavior reveals in the ratio of R k O,α (t, 0) with our choice of a source-sink separation.
In each plot, a solid line indicates the average value over range of 6 ≤ t/a ≤ 9 and a shaded band displays one standard deviations estimated by the jackknife method.
The obtained values of the bare couplings g O , which are not yet renormalized, are summarized in Table IV.
We next evaluate the renormalization factor of the local vector current Z V through a property of g ren V = 1 for the renormalized value of the vector charge g ren V under the exact isospin symmetry. Therefore, the value of Z V can be evaluated by 1/g bare Figure 4 plots the result of Z V , which shows a good plateau in the time range of 2 ≤ t/a ≤ 14. The constant fit result with one-standard-error band denoted by three solid lines shows a good consistency with the value of Z V (red line) obtained by the Schrödinger functional (SF) scheme at the vanishing PCAC quark mass [50]. This consistency may suggest that the excited-state contamination in three-point functions is not serious with our choice of a source-sink separation. We also draw the value of Z A with the SF scheme for comparison in the same figure. The difference between Z V and Z A is 1.5%, which indicates a fairly small chiral symmetry breaking effects in our simulations.
The local axial-vector current is renormalized with the value of Z A = 0.9650(68)(95) obtained with the SF scheme [50] shown in Fig. 4. We then plot the renormalized value of the axial charge g ren as a function of the current insertion time slice t in Fig. 5. Three solid lines represent the fit result with one-standard error band in the time region of 6 ≤ t/a ≤ 9, while the uncertainty in determination of the value of Z A is take into account by shaded band. We finally obtain the renormalized value of the axial charge which slightly underestimates the experimental value of 1.2723(23) by less than 10%. However, the dominant statistical error is of the same order. We also recall that t-dependence of R 5z A,3 (t, 0) in Fig. 3 shows large wiggles, which seem to break time reversal between the source and sink points. However the time-reversal feature should eventually be restored as observed in the vector and tensor channels. This suggests that the ratio of R 5z A,3 (t, 0) still has large statistical fluctuations, which are not well under control. A new class of statistical error reduction techniques such as low-eigenmode deflation and all-modeaveraging (AMA) [51] should be useful and efficient to improve the statistical accuracy of R 5z A,3 (t, 0) with the limited number of configurations. We intend to examine this possibility in future research. The isovector electric form factor G E and magnetic form factor G M are separately obtained by the ratios (23) and (24). Figure 6 shows ratios of three-and two-point functions (23) and (24) as a function of the current insertion time slice t for the isovector electric form factor G E (left) and magnetic form factor G M (right) at the nine lowest nonzero momentum transfers after extracting the relevant kinematical factors. Although the plateau region is not clearly defined at the higher momentum transfers, the time region of 6 ≤ t/a ≤ 9 certainly exhibits an asymptotic plateau at low momentum transfers with our choice of a source-sink separation. Therefore, we evaluate the values of both electric and magnetic form factors by constant fits to the data points in the range of 6 ≤ t/a ≤ 9 denoted by the gray shaded area as same as the cases of the coupling g O (O = V, A, S, T ).
In Table V, we compile the results of G E and G M form factors together with the results of F 1 and F 2 form factors, which are evaluated by linear combinations of G E and G M with appropriate coefficients determined by Eqs. (23) and (24) with measured values of the momentum squared q 2 and the nucleon mass M N . The q 2 dependence of G E (left panel) and G M (right panel) are plotted in Fig. 7, respectively. For the finite three momentum q, we use the nine lowest nonzero momenta: q = π/48 × n with |n| 2 ≤ 9. The lowest nonzero momentum transfer q 2 in present work reaches the value of 0.024(1) [(GeV) 2 ], which is much smaller than 4m 2 π even at m π ≈ 146 MeV. In each panel, we also plot Kelly's fit [52] (red dashed curves) as their experimental data. The simulated electric form factor G E is fairly consistent with the experiment, especially at low q 2 . The magnetic form factor G M agrees with the Kelly's fit albeit with rather large errors.
This feature suggests that our results successfully reproduce experimental value for the electric RMS radius. On the other hand, noisier signals for G M would hinder us from the precise determination of the magnetic RMS radius. The electric (magnetic) RMS radius, r 2 E ( r 2 M ), can be evaluated from the slope of the respective form factor at zero momentum transfer with l = E (M ). Recall that G M at zero momentum transfer, whose value corresponds to the isovector magnetic moment µ v = G M (0), cannot be directly measured for the kinematical reason [23]. Therefore, q 2extrapolation is necessary to evaluate the value of G M (0). In principle, low q 2 data points are crucial for determination of both the RMS radii and magnetic moment. However, the accessible momentum transfer is limited on the lattice since the components of three momentum are discrete in units of 2π/L. In this sense, a large spatial size L is required for precise determination of the RMS radii ( r 2 E and r 2 M ) and magnetic moment (µ v ). The determination of the slope (or q 2 -extrapolation) is highly sensitive to how we fit the q 2 -dependence of the form factors. In the previous studies [9,13,[24][25][26], the dipole form G(q 2 ) = a 0 /(1 + a 1 q 2 ) 2 , and the polynomial form G(q 2 ) = k=0 d k q 2k have been often adopted for We also plot Kelly's fit [52] as their experimental data.
G E and G M . However, the fitting form ansatz may tend to constrain an interpolation (or extrapolation) and introduce a model dependence into the final result of the RMS radius (or the value of G(0)). In order to reduce systematic errors associated with the interpolation or extrapolation of the form factor in momentum transfer, we use the model-independent z-expansion method [53,54].

Analysis with z-expansion method
Let us recall the analytic structure of the form factors in the complex q 2 -plane. There is a non-analytic region for G(q 2 ) due to threshold of two or more particles, e.g. 2π continuum. Hence the q 2 -expansion, G(q 2 ) = k=0 d k q 2k , does not converge if |q 2 | > 4m 2 π where q 2 = −4m 2 π is a branch point associated with the pion pair creation for G = G E or G M [54].
The z-expansion (denoted as z-Exp) makes use of a conformal mapping from q 2 to a new variable z. Since this transformation makes the analytic domain mapped inside a unit-circle |z| < 1 in the z-plane, the form factors can be described by a convergent Taylor series in terms of z: where k max truncates an infinite series expansion in z [63]. Here, t cut = 4m 2 π , where m π corresponds to the simulated pion mass, is chosen for the vector channel (G = G E or G M ), while t cut = 9m 2 π , which is associated to 3π continuum, will be later chosen for the axial-vector channel (G = F A ) [18]. We here note that k max must ensure that terms c k z k become numerically negligible for k > k max for a modelindependent fit. Although |c k /c k−1 | < 1 is expected for sufficiently large k, the range of possible values of k max is limited by the condition k max ≤ 8 (7) for the G E (G M ) form factor that are calculated at totally ten (nine) data points in this study. We then first check the stability of the fit results with a given k max . Figure 8 shows the fit results for all ten data points of G E (q 2 ) with k max = 2, 3 and 8 in the z-expansion. Similarly, we also plot the fit results for all nine data  points of G M (q 2 ) with k max = 2, 3 and 7 in Fig. 9. At a glance, one can see from Fig. 8 and Fig. 9 that the curvature becomes smaller in z variable than q 2 variable for both cases of G E and G M . This indicates that a power series in terms of z is clearly more efficient than that of q 2 . It is also observed that both fit results are not sensitive to the choice of k max . This suggests that the z-expansion gives a rapid convergence series for both For GE, the ratios of |c k /c k−1 | with k ≥ 2 become zero within the statistical error, while the ratios of |c k /c k−1 | for GM reach a convergence value less than unity at k ≈ 5 or 6. cases. Indeed, as shown in Fig. 10, the values of |c 1 /c 0 | are insensitive to the choice of k max if k max ≥ 3. For G E , the ratios of |c k /c k−1 | with k ≥ 2 become zero within the statistical error, while the ratios of |c k /c k−1 | for G M reach a convergence value less than unity at k ≈ 5 or 6. Thus, the value of k max ≤ 7 is large enough to guarantee that the z-Exp analysis makes a model independent fit and reduces systematic uncertainties to determine the RMS radii and the value of G M (0). For these reasons, k max = 8 (7) is hereafter chosen for the G E (G M ) form factor in the z-Exp method.
In Fig. 11, we show the fit results of G E (q 2 ) (left panel) and G M (q 2 ) (right panel) with three types of fits: dipole (green dashed curve), quadratic (blue dot-dashed curve) and z-Exp (red solid curve) fits. All the fits reasonably describe all ten (nine) data points for G E (G M ) with χ 2 /dof < 1. However, if one takes a closer look at low q 2 , the fit curve given by the z-expansion fit goes nicely through the data points in low q 2 region. This implies that the z-Exp fit is relatively insensitive on the higher q 2 data points as was expected. The RMS radius is determined to be r RMS = −6(c 1 /c 0 )/(4t cut ) (z-Exp fit), √ −12a 1 (dipole fit) and −6d 2 /d 0 (quadratic fit), while the coefficients of c 0 , a 0 and d 0 for G M correspond to the value of the magnetic moment, respectively. All results obtained by various fits are compiled in Table VI. In Fig. 12, we compare the results of r 2 E (left panel) and µ v (right panel) from the z-Exp fit (red square) with those of quadratic (blue circle) and dipole (green diamond) fits. In the left panel, two horizontal bands represent experimental results from ep scattering (upper) and muonic hydrogen (µ-H) atom (lower). The dipole results are underestimated in comparison to the corresponding experimental values. Although each z-Exp fit result has relatively larger error than the other results, its error may include both statistical and systematic uncertainties without any model dependence. Moreover each result of r 2 E and µ v from the z-Exp fit is most con-sistent with respective experiment. As our final results, we then quote the value of the (isovector) electric RMS radius: r 2 E = 0.915 ± 0.099 fm (37) and the value of the (isovector) magnetic moment: which are obtained from the z-Exp fits. The former is consistent with both the experimental values from ep scattering and µ-H atom spectroscopy, though statistical uncertainties should be reduced down to a few percentages so as to resolve the experimental puzzle on the proton size.

Comparison with previous results
We finally compare our results of r 2 E and µ v with recent lattice results from LHPC [9], PNDME [24], the Mainz group [25], ETMC (denoted as ETMC'13 [26] and ETMC'17 [10]) and RBC-UKQCD (denoted as RBC-UKQCD'09 [13] and RBC-UKQCD'17 [27]) as shown in Fig. 13. See also Table VIII for a summary of previous lattice calculations in comparison with this study. Although our results are compatible to both experiments albeit with large errors, many results of both the electric RMS radius and magnetic momentum are often underestimated relative to respective experimental values in the previous calculations as can be seen in Fig. 13.
The following are major points of concern for this issue: 1) either quantities are sensitive to the simulated pion mass, 2) their finite size effects become significant with the pion mass decreasing, 3) the longer q 2 interpolation or extrapolation to estimate them causes larger systematic uncertainties. The last point is connected to the second point since the larger spatial volume enables  us to access the nucleon form factors at smaller momentum transfers. In this context, the LHPC calculation, where the simulations are performed on the largest spatial size of 5.6 fm with almost physical quark masses, is a favorite among the previous calculations. Indeed, the LHPC results show better agreements with the experiments, though the electric RMS radius is slightly smaller than the PDG value. The simulated pion mass in the LHPC calculation is comparable to that of ours. Our lattice setup is, however, superior to the LHPC calculation with respect to our spatial size of 8.1 fm. The larger spatial volume provides rich information about momentum squared dependence of the nucleon form factors in the low q 2 region. Raw data of F 1 and F 2 are available in tables of Ref. [9]. Therefore, it is worth comparing our results of F 1 and F 2 form factors with their results in the same plots.
In Fig. 14, we show both the LHPC results (cross symbols) and our results (open circles) for the Dirac F 1 (left panel) and Pauli F 2 (right panel) form factors as a function of momentum squared q 2 . In each panel, the dashed curves correspond to Kelly's fit [52] as their experimental data. In the LHPC results, they use two types of kinematics regarding momentum p ′ for the final state of the nucleon state when constructing the three-point correlation functions (22). The smallest momentum transfer q 2 min = 0.044 [(GeV) 2 ] in the LHPC calculation is given by the choice of p ′ = 2π/L × (−1, 0, 0), while we consider only the rest frame of the final state with p ′ = 0.
Our results and the LHPC results are consistent with each other in both F 1 and F 2 form factors within the statistical errors. We have found that the size of our statistical errors is similar to that of the LHPC data points which are calculated in the rest frame of the final state with p ′ = 0. It is clear that q 2 -dependence of the form factors at low q 2 are much constrained by our results that contain the smallest nonzero q 2 data point. We primary focus on the values of r 2 1 and κ v since the data of the F 2 form factor exhibits large statistical fluctuations in the low q 2 region.
We then extract isovector Dirac RMS radius ( r 2 1 ) and anomalous magnetic moment (κ v = F v 2 (0)) together with the Pauli RMS radius ( r 2 2 ) from the F 1 and F 2 form factors with the various fits as summarized in Table VII. Our results of r 2 1 and κ v obtained from the z-Exp fits are more consistent with both experiments albeit with large errors. The dipole fits tend to yield smaller errors in comparison with the other fits. For comparison, the LHPC results of r 2 1 , F v 2 (0) and r 2 2 are also listed in Table. VII. Their quoted errors on both r 2 1 and κ v are, however, small by a factor of two or three in comparison with our results obtained by analysis with the z-Exp fit. This is partly because the LHPC results are given by the dipole fits with large number of q 2 data points. Indeed, if we adopt the dipole fit rather than the z-Exp fit, the statistical uncertainties on the obtained results become small as those of the LHPC results. Al-though the LHPC results are also roughly consistent with experiments, the shorter q 2 interpolation (or extrapolation) that was achieved in our study by the use of larger volume reduces the systematic uncertainties on the determination of r 2 1 (or κ v ).
C. Results of nucleon form factors in the axial-vector and pseudoscalar channel

Axial-vector FA and induced pseudoscalar FP form factors
In the axial-vector channel, two independent form factors F A and F P can be extracted from kinematically different types of the three-point functions (25). Recall that the z-direction is chosen as the polarized direction through the definition of the projection operator P 5z . Therefore, the longitudinal momentum (q 3 ) dependence explicitly appears in Eq. (25). The three-point functions are classified with the transverse (i = 1 or 2) and longitudinal (i = 3) components. Furthermore, for the case of i = 3, there are two types of kinematics, either q 3 = 0 or q 3 = 0.
We take advantage of the choice of either transverse or longitudinal directions together with explicit q 3dependence so as to separately obtain F A and F P from Eq. (25) in line with Ref. [23]. For convenience, using the ratio of Eq. (25) we define and In the plateau region of Λ A L,T (t, q) we determine the matrix elements of the axial-vector current which has the following relation to the form factors: The most simple way is that F A is obtained from Λ A L (q) with q 3 = 0, while F P is evaluated from Λ A T (q)/M N . However, due to the kinematics, q 3 = 0 is not available for Q3, Q6 and Q9 (labels defined in Table II), while Λ A T (q) is not calculable for Q1, Q4, Q8. We then determine the two form factors through the following equations which explicitly depend on the longitudinal momentum q 3 : Here, N f denotes the number of dynamical quark flavors. In the fourth column, TM-Wilson (TM-Clover) stands for the twisted mass Wilson (clover) Dirac operator, DWF for the domain-wall fermions. In the last column,"R", "S", "TSF" and "GPoF" stand for the standard ratio method, the summation method, the two-state fit method and the generalized pencil-of-function method.
Although both ways are available for Q2, Q5 and Q7, we choose the former in this study. Figure 15 shows ratios of three-and two-point functions (43) and (44) as a function of the current insertion time slice t for the axial-vector form factor (F A ) (left) and induced-pseudoscalar form factor (F P ) (right) at the nine lowest nonzero momentum transfers. The latter is multiplied by 2M N for dimensionless quantity, while both ratios are renormalized with Z A = 0.9650(68), which is given in the SF scheme [50]. We evaluate the values of both axial-vector and induced-pseudoscalar form factors by constant fits to the data points in the time region of 6 ≤ t/a ≤ 9 denoted by the gray shaded area.
We next show the q 2 -dependence of the renormalized form factors, F ren A (left panel) and F ren P (right panel), in Fig. 16. The values of F ren A (q 2 ) and 2M N F ren P (q 2 ) are also summarized in Table IX. In the left panel, the experimental curve is given by a dipole form with the RMS radius of 0.67(1) fm [16,17] for a comparison. At a glance, the F ren A form factor is barely consistent with the experimental curve in the whole region of measured momentum transfers within the current statistics, except for two points at the smallest and second smallest momentum transfers including the axial-vector coupling g A = F A (0). In the right panel, two experiments of muon capture and pion-electroproduction are marked as blue diamond and green asterisk symbols. Our result of F P (q 2 ) is significantly underestimated in comparison with both experiments especially in the low q 2 region. The F P form factor is expected to have a pion pole that dominates the behavior near zero momentum transfer. The red dashed curve in the right panel is given by the pion-pole dominance (PPD) model, where the induced pseudo-scalar form factor is given as which functional form becomes justified by the generalized GT relation (17) at least in the chiral limit (m = 0). Indeed, two experiments of muon capture and pion-photo production follow a curve given by the PPD model. According to the PPD model (45), the observed quenching effect in the value of F P (q 2 ) might be partly due to the underestimation of the axial-vector coupling g A obtained in this study. We next evaluate the axial RMS radius r 2 A from the slope of F A (q 2 )/F A (0) at zero momentum transfer. Three types of fits are used as same as the cases of r 2 E and r 2 M . First of all, Figure 17 shows the z-Exp   [50]. In the left panel, the experimental curve is given by a dipole form with the RMS radius of 0.67(1) fm [16,17], while the red dashed curve in the right panel is given by the pion-pole dominance model. The experimental data points in the right panel are given by muon-capture [55] and pion-electroproduction [21].  fit results for all ten data points of F A (q 2 )/F A (0) with k max = 2, 3 and 8 in the z-expansion. In Fig. 17, the ratio of F A (q 2 )/F A (0) is plotted as a function of the variable z, which is defined by Eq. (36) with t cut = 9m 2 π . Before discussing the stability of the z-Exp fit, we remark that the statistical uncertainties on F A (q 2 )/F A (0) at the lower momentum transfers are extremely suppressed since statistical fluctuations in the numerator and denominator are highly correlated. This observation indicates that the underestimation of F A (q 2 ) at the second smallest momentum transfer compared with the experimental value found in Fig. 16 can be attributed to the reduction of the axial-vector coupling g A .
As shown in three panels of Fig. 17, we again confirm that the fit results are not sensitive to the choice of k max as same as the cases of G E and G M . Therefore, the value of k max = 8 is large enough to guarantee that the z-Exp analysis makes a model independent fit.
For comparison, we next show the fit results of F A (q 2 )/F A (q 2 ) with three types fists: dipole (green dashed curve), quadratic (blue dot-dashed curve) and z-Exp (red solid curve) fits in Fig. 18. All results of the axial RMS radius obtained by three fits are complied in Table X. Although all three fits are mutually consistent with each other because of their large statistical errors, the z-Exp fit tends to give a higher value as r 2 A = 0.46 (11) fm, that is closer to the experimental value of 0.67(1) fm [64].

Pseudoscalar GP form factor
As described in Sec. II, it is theoretically known that two form factors F A and F P in the axial-vector channel are not fully independent. Theoretically, q 2 dependences of F A and F P should be constrained by the generalized GT relation (17) with the pseudoscalar form factor G P as a consequence of the axial Ward-Takahashi identity. To test whether the correct behavior of the generalized GT relation is well satisfied in our simulations, we also calculate the pseudoscalar matrix element, which is described by the single form factor G P (q 2 ) defined in Eq. (18). Figure 19 shows ratio of three-and two-point functions (26) as a function of the current insertion time slice t for the pseudoscalar form factor G P (q 2 ) at the nine lowest nonzero momentum transfers after extracting the relevant kinematical factors. The plateau region is clearly defined even at the higher momentum transfers with our choice of a source-sink separation. We thus evaluate the values of the pseudoscalar form factors by constant fits to the data points in the time region of 6 ≤ t/a ≤ 9 denoted by the gray shaded area as same as the cases of  other form factors. We next plot the bare pseudo-scalar form factor G P (q 2 ), which is not renormalized, in Fig. 20. The measured q 2 dependence of G P (q 2 ) resembles that of F P (q 2 ), where the relatively strong q 2 dependence appears in the lower q 2 region as was expected from the pion-pole contribution. In the PPD model, the pion-pole dominance holds even in G P (q 2 ). Combined with Eq. (17) and Eq. (45), a naive pion-pole dominance form G PPD P (q 2 ) is given as Thus one may realize that the ratio of the PPD forms, G PPD P and F PPD P does not depend on q 2 and gives the low-energy constant B 0 as with a help of the Gell-Mann-Oakes-Renner (GMOR) relation for the pion mass: m 2 π = 2B 0m . As shown in Fig. 21, the ratio of G P (q 2 )/F ren P (q 2 ) indeed exhibits a flat q 2 dependence at lower q 2 . We then estimate the low-energy constant B 0 by a constant fit to the plateau value using six data points at the lower momentum transfer. We then get the bare value of the lowenergy constant as B 0 = 3.10 (25) [GeV], which is represented by blue solid line with a shaded band in Fig. 21. This value is fairly consistent with the one evaluated by the GMOR relation with the simulated pion mass and a (bare) quark mass (am PCAC = 0.001577(10)) obtained from the pion two-point correlation functions with the PCAC relation [11,50]. This observation strongly indicates that the G P (q 2 ) form factor shares a similar pionpole structure with the F ren P (q 2 ) form factor and the ratio of their residues is highly constrained by the GMOR relation.

Test for the axial Ward-Takahashi identity
In order to verify the axial Ward-Takahashi identity in terms of the nucleon matrix elements, we define the following ratio inspired by the generalized GT relation (17) with the simulated nucleon mass M N . If the ratio reveals a q 2 independent plateau in the entire q 2 region, m AWTI should be regarded as an alternative (bare) quark mass. As shown in Fig. 22, the ratio m AWTI has no appreciable q 2 dependence. It indicates that three form factors well satisfy the generalized GT relation. Using data points at the six lowest momentum transfers, we can read off am AWTI = 0.00451 (48), which is roughly 3 times heavier than the PCAC quark mass [11,50]. Since the PCAC relation is also a consequence of the axial Ward-Takahashi identity, a relation m PCAC ≈ m AWTI was expected. However, it is not the case.
This issue seems to be associated with finite lattice spacing artifacts since the above consideration does not take into account O(a) improvement of the axial-vector current. The renormalized O(a)-improved operator takes the form The ratio of GP (q 2 )/F ren P (q 2 ) as a function of momentum squared q 2 . A flat q 2 dependence is observed due to the fact that the GP (q 2 ) form factor shares a similar pionpole structure with the F ren P (q 2 ) form factor. A blue solid line represents the (constant) fit result and a shaded band denotes the fit range and statistical error estimated by the jackknife method, while a red dashed line represents the bare low-energy constant evaluated by the ratio of m 2 π /(2mPCAC).
where A + α (x) =ū(x)γ α γ 5 d(x) and P + (x) =ū(x)γ 5 d(x). Strictly speaking, this improved axial-vector current satisfies the axial Ward-Takahashi identity on the lattice: . Therefore, the generalized GT relation (17) is slightly modified by the presence of O(a) improvement term in Eq. (49). However, the term proportional to c A only contributes to modification of the F P form factor as F imp P (q 2 ) = F P (q 2 ) + ac A G P (q 2 ), and then the modified GT relation is expressed by Starting from the modified GT relation (50), the (bare) quark massm becomes where the first term is nothing but the ratio m AWTI defined in Eq. (48) and the second term corresponds to a correction term due to O(a) improvement of the axialvector current. Although the ratio m AWTI may indeed receive the O(a) correction, which yields linear dependence of q 2 , the effect of O(a) improvement is insignificant in the low q 2 region. Furthermore, the improvement coefficient c A is found to be very small with our choice of c SW = 1.11 [33]. The improvement coefficient c A is expected to be of the order of a few 0.01 in lattice units [33]. In Fig. 23, we plot the ratio m imp AWTI defined in Eq. (51) as a function of momentum squared q 2 . Circle symbols represent the unimproved results obtained in Eq. (51) with c A = 0.0 in lattice units. After adopting the renormalized O(a)-improved operator (49), the central value of the unimproved results at each q 2 point is likely to move in the brown band, which represents the region between the lower (c A = 0.03 in lattice units) and upper (c A = −0.03 in lattice units) limit. Figure 23 indicates that the systematic uncertainties that arise from the O(a) improvement term in the axial-vector current are much smaller than the current statistical errors in this study. The large deviation from a blue dashed line, which denotes the PCAC quark mass, mostly stays the same. We thus conclude that the issue of m AWTI ≫ m PCAC is not directly related to finite lattice spacing artifacts. We rather speculate that this issue is connected with modification of the pion-pole structure in both the F P (q 2 ) and G P (q 2 ) form factors as will be described in the next subsection. The improvement coefficient is set to c A = 0 in our whole analysis, hereafter.  [11,50] in lattice units.

Distortion of pion-pole structure
In the previous subsections, we have observed that q 2 dependences of F A , F P and G P are well constrained by the generalized GT relation, while the "pion-pole" structure of F P and G P are likely shared with each other. However, the bare quark mass (m AWTI ) evaluated from the ratio (48) is roughly 3 times heavier than the value of m PCAC . This puzzle would be related to the discrepancy between our result of F P (q 2 ) and experiments.
We first speculate that the measured F P (q 2 ) can be described by the following PPD inspired form: with the pole mass of m pole , which is not necessarily con- strained by our simulated pion mass m π . If this functional form well describes the measured F P (q 2 ), we may define the effective "pion-pole" mass from F P (q 2 ) as which would exhibit a flat q 2 dependence. In Fig. 24, we plot the effective pole-mass plot as a function of q 2 . A horizontal dashed-dotted line represents the value of the simulated pion mass m π in lattice units. Clearly, there is no appreciable q 2 dependence of the effective pole mass m pole . Especially, the six data points at lower q 2 show a fairly good consistency among them within statistical errors. We then get the pole mass value as am pole = 0.1099(74), which is given by a constant fit to data in the shaded region. However, the value of m pole = 256(17) MeV is roughly twice heavier than the simulated pion mass (m π ≈ 146 MeV). This indicates that the "pion-pole" structure in F P (q 2 ) is modified with the larger pole mass. We then define a modification factor α pole as follows and obtain its value as α pole = 3.05(41) that can account for the discrepancy between m AWTI and m PCAC through the GMOR relation. We thus estimate an improved value of the bare quark mass as which should be really consistent with the value of m PCAC .  (17) 3.05(41) 0.00145 (12) We next plot the quantity of m quark in lattice units as a function of q 2 in Fig. 25. Again, there is no appreciable q 2 dependence especially at low q 2 . The horizontal dashed line represents the value of m PCAC , while the solid line indicates the fit result of m quark and shaded band displays the fit range and 1 standard deviation. The value of m quark is obtained as am quark = 0.00145 (12), (56) which is in good agreement with the value of am PCAC = 0.001577 (10). The importance of our findings is twofold: 1) our results of all three form factors, F A (q 2 ), F P (q 2 ) and G P (q 2 ), are subjected to the generalized GT relation (17) as a consequence of the axial Ward-Takahashi identity, 2) the discrepancy between our result of F P (q 2 ) and experiments is mainly caused by distortion of the pion-pole structure in both F P (q 2 ) and G P (q 2 ), though the reason is not yet known. Therefore, in this work, we prefer not to estimate the pseudoscalar coupling g P and pionnucleon coupling g πN N , since both quantities are highly sensitive to the pole structure of F P (q 2 ).  [11,50] in lattice units.
MeV) on a large spatial volume (8.1 fm) 3 . We first analyze the vector nucleon matrix element, which is described by the electric (G E ) and magnetic (G M ) form factors. We carefully examine both G E and G M form factor shapes with a model-independent analysis based on the z-expansion method. As a result, we obtain the isovector electric RMS radius r 2 E and magnetic moment µ v , which are consistent with experimental values. We would like to emphasize that the former quantity is, for the first time, consistent with both the experimental values from ep scattering and µ-H atom spectroscopy, though further reduction of statistical and systematic errors (including electromagnetic and isospin breaking effects) down to less than 1% is required to resolve the experimental puzzle on the proton size.
We have also analyzed axial-vector nucleon matrix element, which is described by the axial-vector (F A ) and induced pseudoscalar (F P ) form factors. We have found that although the axial charge g A = F A (0) is slightly underestimated in comparison with the experimental value, the axial vector form factor F A is barely consistent with experiments in the low q 2 region (0 ≤ q 2 ≤ 0.2 [(GeV) 2 ]). However, the pseudoscalar form factor F P is considerably underestimated at very low q 2 .
To understand this issue, we have, in addition, calculated the pseudoscalar matrix element, which is described by a single form factor of G P called as the pseudoscalar form factor. We examine q 2 dependences of three form factors, F A , F P and G P , which should be constrained by the generalized GT relation as a consequence of the axial Ward-Takahashi identity. We have observed that the G P form factor shares a similar "pion-pole" structure with the F P form factor. The "pion-pole" structure is, however, found to be distorted due to the larger pole mass than the simulated pion mass. If this effect is taken into account for an estimation of the bare quark mass by using three form factors through the generalized GT relation, we can fully verify the axial Ward-Takahashi identity in terms of the nucleon matrix elements. The bare quark mass obtained in this study is consistent with an alternative quark mass obtained from the pion twopoint correlation functions with the PCAC relation. We, however, have not yet fully understood the origin of the "pion-pole" distortion found in the F P and G P form factors. The similar issue in the axial-vector channel was reported in Ref. [57].
The PACS collaboration is generating new ensembles on a significantly large spatial volume of (10.8 fm) 3 at the physical point (m π ≈ 135 MeV) [58]. We then plan to develop the present calculation for more precise determination of the nucleon form factors and also address all of unsolved issues described in this paper. Such planning is now underway.