Nucleon form factors and root-mean-square radii on a (10.8 fm$)^4$ lattice at the physical point

We present the nucleon form factors and root-mean-square (RMS) radii measured on a (10.8 fm$)^4$ lattice at the physical point. We compute the form factors at small momentum transfer region in $q^2\le 0.102$ GeV$^2$ with the standard plateau method choosing four source-sink separation times $t_{\rm sep}$ from 0.84 to 1.35 fm to examine the possible excited state contamination. We obtain the electric and magnetic form factors and their RMS radii for not only the isovector channel but also the proton and neutron ones without the disconnected diagram. We also obtain the axial-vector coupling and the axial radius from the axial-vector form factor. We find that these three form factors do not show large $t_{\rm sep}$ dependence in our lattice setup, and those RMS radii are consistent with the experimental values. On the other hand, the induced pseudoscalar and pseudoscalar form factors show the clear effects of the excited state contamination, which affect the generalized Goldberger-Treiman relation.

We present the nucleon form factors and root-mean-square (RMS) radii measured on a (10.8 fm) 4 lattice at the physical point. We compute the form factors at small momentum transfer region in q 2 ≤ 0.102 GeV 2 with the standard plateau method choosing four source-sink separation times tsep from 0.84 to 1.35 fm to examine the possible excited state contamination. We obtain the electric and magnetic form factors and their RMS radii for not only the isovector channel but also the proton and neutron ones without the disconnected diagram. We also obtain the axial-vector coupling and the axial radius from the axial-vector form factor. We find that these three form factors do not show large tsep dependence in our lattice setup, and those RMS radii are consistent with the experimental values. On the other hand, the induced pseudoscalar and pseudoscalar form factors show the clear effects of the excited state contamination, which affect the generalized Goldberger-Treiman relation.

I. INTRODUCTION
For the deep understanding of the nucleon and nucleus structures, a precise determination of structure functions is an essential ingredient. Recently an unknown effect for the proton charge radius has been revealed as a significant discrepancy between different approaches [1] in the ep scattering [1] process and the muonic hydrogen spectroscopy [2], in which 5.6-σ deviation appears as the socalled "Proton radius puzzle." The measurement of the atomic spectroscopy [1] has also agreed with the value from the ep scattering, while a recent measurement of the regular hydrogen spectroscopy [3] agrees with the value from the muonic hydrogen spectroscopy. Under such a confusing circumstance, the theoretical estimation is demanded as an independent test.
Similarly, the axial-vector form factor and the axial radius are important inputs for the weak process associated with the neutrino-nucleus scattering. The q 2 dependence of the axial-vector form factor can be used to estimate the neutrino properties such as the neutrino mass and mixing angle [4]. Furthermore, the axial-vector coupling g A obtained from the cross section of the muon-nucleus scattering measured in the muon capture experiment serves as an independent test to check consistency with the highprecision data of g A from the neutron beta decay. Having achived the three times higher precision from the current measurement in the muon capture experiment [5], it would also provide g A comparable with that from the neutron beta decay, which is expected to be less than 1% level, using an input of accurate axial radius [4]. The current experimental value of the axial radius [6] is 3% accuracy from the dipole fit of the neutrino-deuteron scattering and the pion electroproduction experiment, whereas, according to the argument in Refs. [4,7], this error may be underestimated by the model-dependent analysis. It means that the theoretical value of the axial radius is desired to use as an input for the analysis in both muon capture and neutrino scattering experiments.
Lattice QCD (LQCD) is able to nonperturbatively determine the QCD values of the nucleon form factors and RMS radii from the first principles. The recent developments of computational techniques and incredible growth of computational resources make it possible to perform a realistic LQCD computation at the physical point with the light (degenerate up and down) and strange quark flavors, even for the baryonic system in which the systematics involved are more complicated than the mesonic system, e.g., the finite volume effect etc.. Although the results for the nucleon form factors at the physical point have been made available by several LQCD groups [8][9][10][11][12][13][14][15], the precision has not been enough to be comparable with the experimental values. This is due to the exponential growth of the statistical noise as the light quark mass gets close to the physical point, besides the possible systematic uncertainties of the finite volume (FV) effects, the excited state contamination and the lattice cutoff effects we should take into account. Some LQCD groups [8,9,13,14] have tried to subtract the excited state contamination by introducing the "2-state ansatz" [8] and the simultaneous fit of the data off the physical point with the use of the ansatz, e.g., chiral perturbation theory [13] or polynomial functions, to remove the FV effects and the lattice cutoff effects. This approach, instead, is an introduction of the other systematic uncertainties originating from the model dependence. For the purpose of high precision to a few % level and below, much effort to remove the systematic uncertainties in LQCD simulations is needed. We think the most reliable way is the direct simulation at the physical point on sufficiently large volume, which is the critical importance for LQCD computation of the nucleon form factors and RMS radii to directly compare experimental values and theoretically verify the prediction in effective models [16][17][18][19].
This work is an extension of our earlier study [15]. In the previous work, one of the authors analyzed the isovector electric (G E ) and magnetic (G M ) form factors, and obtained the isovector electric RMS charge radius r 2 E in a model-independent way with the z-expansion method. Their results were consistent with two experimental values of ep scattering [1] and µH atom spectroscopy [2] within 1-σ statistical error, although statistical error was much larger than experimental ones. It was then difficult to argue which experimental values can be favored in LQCD. For the magnetic moment µ v , LQCD results were also in agreement with experimental value although its statistical precision was almost 15%. On the other hand, somewhat puzzling situation in the axialvector (F A ) and induced pseudoscalar (F P ) form factors in the nucleon axial-vector matrix element occurred. F A was barely consistent with the experimental results in the low q 2 region of 0 ≤ q 2 ≤ 0.2 GeV 2 and the axial charge g A = F A (q 2 = 0) was slightly underestimated in comparison with the experimental value. Furthermore, F P had an apparent discrepancy from the experimental expectation at very low q 2 , which was a consequence of the violation of the generalized Goldberger-Treiman relation.
The purpose of this paper is further reduction of the statistical and systematic errors for the nucleon form factors and understanding of the issues associated with F A and F P raised in the previous work. We have made several essential improvements from the previous work; (i) The lattice size is enlarged from (8.1 fm) 4 to (10.8 fm) 4 employing the stout-smeared O(a)-improved Wilsonclover quark action and the Iwasaki gauge action at β=1.82 [20], which are exactly same as in the previous work. We expect that the spatial extent 10.8 fm has a strong advantage for both suppression of the finite volume effects and reduction of the minimum value of the momentum transfer to q 2 = 0.013 GeV 2 which is a half of the previous work [15]. (ii) The quark masses are carefully tuned to the physical point [20]. The slight deviation from the physical point with m π = 146 MeV in the previous work [20] is removed by adjusting the pion mass to 135 MeV. (iii) Using the variation of the source-sink separation as t sep /a = t sink /a−t src /a = 10, 12, 14, 16, where the largest one is about 1.35 fm, we can examine the possible excited state contributions, which has not been studied in the previous work [20], where a single choice of t sep /a = 15 was used. (iv) Significant reduction of the computational cost is possible to utilize the all-mode-averaging (AMA) method [21][22][23] optimized by the deflation technique [24]. This paper is organized as follows: Section II presents the definition of nucleon form factors to fix the notations in this paper. The general features of the nucleon form factors have been already explained in Sec. II of Ref. [15]. In Sec. III we first present a brief description of gauge configurations, which are a partial set of "PACS10" configurations generated by the PACS Collaboration [20]. We also explain the error reduction technique employed in this study. The results for the nucleon form factors are presented in Sec. IV. We investigate t sep dependence for the nucleon form factors and the corresponding RMS radii. We also discuss the violation of the generalized GT relation associated with the form factors F A , F P , and G P . Section VI is devoted to summary and outlook.
In this paper the matrix elements are given in the Euclidean metric convention. γ 5 is defined by γ 5 ≡ γ 1 γ 2 γ 3 γ 4 = −γ M 5 , which has the opposite sign relative to that in the Minkowski convention ( γ M = i γ and γ M 0 = γ 4 ) adopted in the particle data group [25]. The sign of all the form factors is chosen to be positive. The Euclidean four-momentum squared q 2 corresponds to the spacelike momentum squared as q 2 M = −q 2 < 0 in Minkowski space.

A. Definition of nucleon form factors
We present our convention of the nucleon form factors. We measure the electric and magnetic Sachs form factors, G E (q 2 ) and G M (q 2 ), which are relevant for the experimental data of elastic electron-nucleon scattering. They are linear combinations of the Dirac and Pauli form factors, F N 1 , F N 2 (N = p, n), as extracted from the nucleon vector matrix elements, The momentum transfer between the nucleon initial state (P ) to the final state (P ′ ) is defined as q = P − P ′ . The electromagnetic current j em α is expressed in terms of the flavor-diagonal vector current: where Q q denotes the charge (in units of proton charge e) for quark flavor q, and the ellipsis denote terms for strange and heavier quarks. Here we rewrite Eq. (4) into the following form for the later discussion, with the isoscalar and isovector currents: j s α =ūγ α u + dγ α d and j v α =ūγ α u −dγ α d. Recall that in the case of m u = m d (SU(2) isospin symmetry), the nucleon threepoint correlation function for the isovector current does not receive any contribution from the disconnected diagram of all quark flavors since they are canceled each other. Therefore, the isovector part of nucleon electromagnetic form factors can be determined only by the connected-type contribution, whose numerical evaluation is much easier in lattice simulations.
On the other hand, the isoscalar component receives the full contribution including the disconnected diagrams even under the exact isospin symmetry. Nevertheless, all disconnected-type contributions from the light and heavier quark flavors are known to be relatively small in comparison to the connected-type contributions especially in the proton (see for example Ref. [26]), whereas it will not in the neutron. Here we simply evaluate individual proton (neutron) form factors in the vector channel, which is extracted from the matrix element with the isoscalar and isovector parts of the electromagnetic current, only by the connected-type contributions, since a computation of disconnected diagram is much costly and beyond our scope of this study. Later we will show the numerical evidence of the appearance of missing effect of disconnectedtype contribution to the electromagnetic form factor in both proton and neutron.
The isovector nucleon form factors can be related to the isovector combination of the proton's and neutron's form factors assuming the SU(2) isospin symmetry where the normalization of the above form factors at q 2 = 0 are given by the proton/neutron electric charge and the magnetic moment: G p E (0) = 1 and G p M (0) = µ p = +2.79285 for the proton and G n E (0) = 0 and G n M (0) = µ n = −1.91304 for the neutron. Therefore, one finds The nucleon axial-vector matrix element is represented with the axial-vector form factor F A (q 2 ) and the induced pseudoscalar form factor F P (q 2 ) as (8) with the axial-vector current, A + α =ūγ 5 γ α d. The axialvector coupling g A = F A (q 2 = 0), which governs the life time of the neutron beta decay, is experimentally determined as g A = 1.2724(23) [27], and the q 2 dependence of F A (q 2 ) is well described by the dipole form [6,28]. On the other hand, the properties of the induced pseudoscalar form factor F P (q 2 ) is less clear in the experiments [29,30].
On the theoretical side, F A (q 2 ) and F P (q 2 ) are related through the generalized Goldberger-Treiman (GT) rela-tion [31,32]: with a degenerate up and down quark massm = m u = m d , derived from the axial Ward-Takahashi (AWT) identity: ∂ α A + α (x) = 2mP + (x). Here, the additional form factor G P (q 2 ) is defined in the nucleon pseudoscalar matrix element with the local pseudoscalar density P + =ūγ 5 d. In order to test the GT relation in LQCD, the simultaneous investigation of three form factors F A (q 2 ), F P (q 2 ) and G P (q 2 ) is essential. The RMS radius R l = r 2 l is defined from the expansion of the normalized form factor G l (q 2 ) in the powers of q 2 : which measures a typical size in the coordinate space.
Here we define G A ≡ F A , and the RMS radius is computed from the derivative of nucleon form factor with respect to q 2 at q 2 = 0, When the form factors can be described by the dipole form, with the dipole mass parameter Λ, the RMS radius R is obtained as R l = √ 12/Λ l . In the z-expansion method [33,34], whose convergence has been carefully studied in Ref. [15], the form factors can be described by a convergent Taylor series in a new variable z which is conformally mapped from q 2 : where the expansion is truncated at the k max -th order with t cut = 4m 2 π for G E and G M , or with t cut = 9m 2 π for F A , and the RMS radius is then determined by r 2 l = −6(c 1 /c 0 )/(4t cut ) for l = {E, M, A}. For a precise determination of the RMS radii r 2 l in LQCD, the data of low q 2 variation is important. The systematic uncertainty due to fitting LQCD data at low q 2 with extrapolation ansatz into q 2 = 0 can be reduced by using small q 2 data we can obtain. It means that large spatial extent is advantageous to this study since the accessible minimum value of q 2 is essentially determined by the spatial extent of the lattice. Our spatial lattice volume of (10.8 fm) 3 allows q 2 min = 0.013 GeV 2 , and it is then 1.8 times smaller q 2 data we can use than (8.1 fm) 3 box employed in the previous work [15]. Furthermore, LQCD analysis concentrating on the small q 2 region allows us to avoid the additional systematic uncertainty of lattice cutoff effect stemming from O((aq) 2 ).

B. Nucleon two-and three-point functions for form factors
We first define the nucleon (proton) interpolating operator as where the superscript T denotes a transposition and C is the charge conjugation matrix defined as C = γ 4 γ 2 .
The indices a, b, c and u, d label the color and the flavor, respectively. The function φ X (X = L, S) represents two types of the smearing functions employed in this study: local type (L) given by and exponentially smeared one (S) by φ S ( The nucleon two-point functions are constructed with the source and sink operators located at t src and t sink , respectively: where the smeared operator is employed at the source and at the sink we use both the smeared (X = S) and local (X = L) operators. The lattice momentum is defined as p = 2π/(N s a) × n with a vector of integers n ∈ Z 3 and N s the number of the spatial lattice sites. A projection operator P + = 1+γ4 2 is applied to extract the contributions from the positive-parity state for |p| = 0 [35,36].
In our study, the nucleon form factor is extracted from the nucleon three-point function, using the local currents with the renormalization factor Z O . In the above equation, q = p − p ′ represents the three-dimensional momentum transfer, and P k denotes the projection operator to extract the form factors for unpolarized case, P k = P t ≡ P + γ 4 and polarized case (in z direction) P k = P 5z ≡ P + γ 5 γ 3 . In a conventional way to remove the unwanted nucleon wavefunction, we use the following ratio, as a function of initial and final nucleon momenta, p and p ′ , and the temporal position of local current t in the fixed temporal position of source and sink nucleon interpolation operator t src and t sink . In this study we use various nucleon source-sink separations as t sep /a = t sink /a − t src /a = 10, 12, 14, 16 to examine a possible excited state contamination. Here, by restricting the kinematics of the nucleon final state at rest, where q 2 = 2M N (E N (q) − M N ) with p ′ = 0, the above ratio can be represented as R O,α (t, q).
In the electromagnetic vector channel, O = j em , the ratio of Eq. (19) is supposed to give the following asymp-totic form [32]: in the limit of t sink ≫ t ≫ t src with N = p, n.
Similarly in the axial-vector current and pseudoscalar cases, O = A + , P + , we obtain The three-point correlation functions of Eq. (18) are calculated by the sequential source method with t sink and t src fixed [37,38], which requires to solve the sequential quark propagators individually in the four choices of t sep and the projection operators P = P t , P 5z .

A. PACS10 configurations
In this work we have used a partial set of PACS10 configurations [20]. We briefly present relevant points to make the paper self-contained (see in Ref. [20] for the detailed description).
The gauge configurations in 2+1 flavor QCD with the stout-smeared O(a)-improved Wilson-clover quark action and the Iwasaki gauge action [39] on an N 3 s × N t = 128 3 × 128 lattice at β = 1.82, which corresponds to a (10.8 fm) 4 physical space-time with a lattice cutoff of a −1 = 2.333(18) GeV (a = 0.08457(67) fm) [40] have been generated by PACS Collaboration. The Schrödinger functional (SF) scheme is employed to determine the nonperturbative improvement coefficient c SW = 1.11 [41]. Since the improvement factor for the axial-vector current c A is consistent with zero within the statistical error [41], we do not take account of the O(a) improvement of the quark bilinear currents. The hopping parameters of (κ ud , κ s ) = (0.126117, 0.124902) are carefully chosen to be at the physical point. We use 20 gauge configurations separated by 10 trajectories. The statistical error is estimated by the single elimination jackknife method.

B. Utilization of all-mode-averaging technique
Here we employ the AMA technique to efficiently implement the LQCD computation of two-and three-point functions. For the implementation of AMA, 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, for instance translational symmetry. N org and N G are the number of such a transformed observable for O (org) and O (approx) respectively. To achieve the high performance of AMA, we need to set N org ≪ N G satisfying the strong correlation r between O org and O (approx) , as 2(1 − r) < 1/N G [22,23]. Following Refs. [23,42] we employ the optimized AMA which adopts the deflated Schwartz Alternative Procedure (SAP) [43] and Generalized Conjugate Residual (GCR) [24] in the computation of both high-precision O (org) and low-precision O (approx) . As demonstrated by the performance test in Ref. [42], the utilization of deflated SAP-GCR can significantly save the computational cost compared to the low-mode deflation originally suggested in Refs. [21,22].

C. LQCD parameters
First we tune the parameters for the source and sink smearing function as A = 1.2 and B = 0.16 in Eq. (16). The smearing parameters are slightly different from previous work [15] to gain a better overlap with the ground state in three-point function. As mentioned in Sec. II A, to avoid the considerable lattice cutoff effect, we choose the eight lowest variations of q 2 listed in Table I, up to q 2 = 0.11 GeV 2 , in our analysis.
The renormalization factors Z O (O = V, A) are obtained by the SF scheme at the vanishing quark mass as Z V = 0.95153(76)(1487), Z A = 0.9650(68)(95) [44], where the first error is statistical one and the second is a systematic error coming from the difference of two volumes. In our analysis this systematic error is regarded as negligible since we here choose the larger volume.
We compute the three-point function of Eq. (18) with four variations of t sep /a = 10, 12, 14, 16 to examine the excited state contamination. Since the LQCD calculation with large t sep suffers from the large statistical noise, we increase the N G as shown in Table II to keep the signalto-noise ratio as t sep becomes large. As for the AMA tuning parameter, the approximation is obtained by 5 GCR iteration using 8 4 SAP domain size with 50 deflation fields. O (org) is given in the stopping criteria of residual norm less than 10 −8 .

D. Nucleon effective mass
In Fig. 1 we show the nucleon effective mass plot with the smeared-smeared and smeared-local operators. The single exponential function is used in the correlated fit with the range of t/a =15-20 for the smeared-local and 13-20 for the smeared-smeared cases, and those are consistent with each other. The he nucleon mass from the smeared-local case is obtained as aM N = 0.4041 (47), M N = 0.9416(110) GeV, (25) where the error is statistical. This is consistent with the experimental value M exp N = 0.93891874 GeV obtained by averaging the proton and neutron masses. For the extraction of the form factors in Eqs. (20), (21), (22) and (23), we use the central values of the nucleon mass and the energy with finite momenta determined from the smeared-local case. In fact, even if we input the statistically fluctuating mass and energy onto those equations, the variation of extracted form factors is negligibly small TABLE I. Choices for the nonzero spatial momenta: q = π/(64a) × n, and corresponding nucleon energy EN measured by global fitting of two-point function with the same range as in Q0 (also see in a text). The degeneracy of |n| 2 due to the permutation symmetry between ±x, ±y, ±z directions is listed in the bottom raw. compared to statistical fluctuation of three-point function. We have summarized the measured nucleon mass and energy in Table I together with values of q 2 .

IV. RESULTS FOR NUCLEON ELECTROMAGNETIC FORM FACTORS AND AXIAL FORM FACTOR
A. Electric form factor and electric charge radius Figure 2 shows t dependence of the ratio R t,N j em ,4 (t, q) of Eq. (20) for the isovector electric form factor G v E (q 2 ) at t sep /a = 10, 12, 14, 16 in the smallest four nonzero momenta corresponding to Q1, Q2, Q3, and Q4 (see Table I). We observe clear plateau for all the cases of t sep and |n| 2 thanks to our elaborate tuning of the smearing parameter. The G v E (q 2 ) is extracted by the constant fit with the fit range listed in Table II. In Fig. 3 we plot the t sep dependence of G v E (q 2 ) for the smallest five values of |n| 2 . One can see the data at t sep /a = 12, 14, 16 is statistically consistent within the error, which means there is negligibly small t sep dependence, while the data at t sep /a = 10 differs from others at smaller nonzero q 2 . This observation allows us to use two possible combined values with t sep /a = {12, 14, 16} and t sep /a = {14, 16} to obtain G v E (q 2 ) without considerable excited state contamination. Figure 4 shows the q 2 dependence of G v E (q 2 ) together with the Kelly's fit [45]. The combined results with t sep /a = {12, 14, 16} and t sep /a = {14, 16} are consistent with each other in all q 2 . In the small q 2 region the LQCD data closely follow the Kelly's fit, while in the large q 2 region it is slightly but systematically above that. The figure also shows that our error is much smaller than the one of our previous results at m π = 0.146 GeV on (8.1 fm) 4 in Ref. [15] thanks to the AMA technique described in Sec. III B and tuning the smearing parameters.
We also calculate G p E (q 2 ) and G n E (q 2 ) separately without the disconnected diagram. They have similar properties to G v E (q 2 ): tiny t sep dependence and good plateau in R t,p jem,4 (t, q) and R t,n jem,4 (t, q). The combined results with t sep /a = {12, 14, 16} and {14, 16} for G p E (q 2 ) and G n E (q 2 ) are summarized in Appendix A together with the ones of G v E (q 2 ). The results for G p E (q 2 ) and G n E (q 2 ) are compared with the Kelly's fit in Fig. 5. We observe that G p E (q 2 ) is closer to Kelly's fit rather than G v E (q 2 ). On the other hand, G n E (q 2 ) is much smaller than the Kelly's fit. One possible reason is an uncertainty due to  r 2 E in the isovector, proton and neutron channels. In the row of "This work" we present our best estimates, where the first error is statistical and the second one is systematic as explained in the text. Results for the proton and neutron are obtained without the disconnected diagram. Our previous work was performed on a 96 4 lattice at mπ = 146 MeV in Ref. [15], where only the statistical errors are presented.
the missing disconnected diagram in the isoscalar channel, which could affect G p E (q 2 ) and G n E (q 2 ). Actually, a recent work using the mixed actions with the domainwall sea quarks and the overlap valence ones implies the contribution of the disconnected diagram to G p E (q 2 ) and G n E (q 2 ) amounts to ∼ 0.005 with rather large statistical errors 1 at q 2 ≈ 0.05 GeV 2 in m π = 135 MeV [26], whose magnitude is comparable to the difference between our G n E (q 2 ) and the experimental value. To completely re- Diamond symbols, which are obtained with tsep/a = 15 on a 96 4 lattice at mπ = 146 MeV in Ref. [15], are also plotted for comparison.
solve the problem we need to evaluate the isoscalar vector form factor in the future.
With the use of the correlated fit procedure, we compare four types of fitting functions to examine the uncertainty in the extrapolation of the slope to q 2 = 0: linear function G E (q 2 ) = d 0 + d 1 q 2 , dipole form of Eq. (13), quadratic function G E (q 2 ) = d 0 + d 1 q 2 + d 2 q 4 and the model-independent z-expansion method with Eq. (15) with k max = 3. In Figs. 4 and 5 we find that the dipole form well describes the LQCD results for G v E (q 2 ) and G p E (q 2 ) up to the maximum fitting range of q 2 cut = 0.102 GeV 2 . We plot the fit form dependence of r 2 E in Fig. 6, where the upper and lower bands denote the experimental results of the ep scattering and the spec-troscopy of the muonic hydrogen (µH) atom, respectively. The numerical results are summarized in Table III together with the experimental values. We observe that all the fit procedures show good consistency within the error bars both for the isovector and proton channels with a reasonable χ 2 /dof, which is evaluated by jackknife estimator in correlated fit. We also find that the combined results with t sep /a = {12, 14, 16} are consistent with those with t sep /a = {14, 16} within the error bars, which indicates that the excited state contamination in G E (q 2 ) is under control. Note that, in the case of neutron, one can find a clear deviation from the experimental value due to the lack of the disconnected diagram as already mentioned above.
As shown in this section, the LQCD calculation at the low q 2 region up to 0.11 GeV 2 allows us to successfully reduce the uncertainties stemming from the choice of the fitting procedures. Their central values, however, slightly fluctuate depending on each fitting procedure and choice of t sep range. We then take a result of r 2 E with the dipole form at t sep /a = {12, 14, 16} as our best estimate of the central value and its statistical error. The maximum difference between the central value in the dipole fit with t sep /a = {12, 14, 16} and those in other fitting procedures with two choices of the combined t sep ranges is taken as the systematic error (see Table III).
Although our result of the isovector channel stays around the value of the µH experiment rather than that of the ep scattering experiment, it may be too early to conclude its preference at this stage because of relatively large error bars. For the proton case, on the other hand, LQCD value stays amid those experimental values. For the definite conclusion we need more precise calculation including the disconnected diagram. This is, however, an encouraging situation indicating a possibility that LQCD can distinguish the two experimental results in near future.

B. Magnetic form factor and magnetic RMS radius
The isovector magnetic form factor G v M (q 2 ) is extracted from the ratio R 5z,N jem,i (t, q) of Eq. (21). The anal- . We first plot the t dependence of the ratio with |n| 2 = 1, 2, 3, 4 for t sep /a = 10, 12, 14, 16 in Fig. 7, which show good plateau for all the cases of |n| 2 and t sep /a. We extract G v M (q 2 ) with the constant fit employing the same fitting range as in the G v E (q 2 ) case. Figure 8 shows that the results for t sep /a = 10, 12, 14, 16 agree with each other within 1-σ error bars, though the statistical fluctuation is much larger than the G v E (q 2 ) case. We evaluate G p M (q 2 ) and G n M (q 2 ) separately from each R 5z,N j em ,i (t, q) for N = p, n, where we omit the disconnected diagram. As in G v M (q 2 ), all the ratios of R 5z,N j em ,i (t, q) have reasonable plateaus, and those values are consistent in the four t sep cases. At each q 2 we take two combined val-    . Those values are summarized in Appendix A. Figure 9 shows that the results from the two combined t sep ranges are consistent with each other. These results are compared with that of our previous calculation [15]. Our current result has much smaller error than the previous one, and closer to the Kelly's fit. In Figs. 9 and 10 we observe that the q 2 dependence of G v M (q 2 ) and G p M (q 2 ) is consistent with the Kelly's fit within the 1.5-σ error, though G n M (q 2 ) for t sep = {12, 14, 16} in the smaller q 2 region shows slight deviation from the Kelly's fit. This could be due to the lack of the disconnected contribution as well as G n E (q 2 ) case. Note that the negative disconnected contribution of ∼ −0.03 at q 2 ≈ 0.05 GeV 2 in m π = 135 MeV implied in Ref. [26] could make our result of G n M (q 2 ) closer to the experimental value. We obtain the magnetic RMS radius r 2 M together with the magnetic moment µ = G M (0) with four types of fitting functions as in the electric case. The numerical values of µ and r 2 M for the isovector, proton and neutron channels are summarized in Table IV with the linear, dipole, quadratic forms and the z-expansion method with k max = 3. In Figs. 9 and 10, one can see that the dipole form up to q 2 cut = 0.102 GeV 2 can well describe the LQCD data for both choices of the combined t sep ranges, and they show good consistency with each other. Figure 11 illustrates a comparison between four types of the fit procedures for the magnetic moment µ, which shows good consistency within 1-σ error as well as the case of the electric RMS radius. The situation is similar to the magnetic RMS radius r 2 M in Fig. 12, though the z-expansion method at t sep /a = {14, 16} gives the deviation beyond 1-σ error between two choices of the combined t sep ranges. We take the result of the dipole form with t sep /a = {12, 14, 16} as our best estimate of the central value and the statistical error for µ and r 2 M , and take the maximum difference between the central value in the dipole fit with t sep /a = {12, 14, 16} and those in other fitting procedures with two choices of the combined t sep range as the systematic error (see Table IV). Compared to the previous work [15], the statistical precision is significantly improved with less discrepancies between four types of the fit procedures. The results of the magnetic moment and RMS radius show consistency with the experimental values within 1-σ error bars of the isovector, proton, and neutron channels, though the systematic uncertainty in a choice of the combined t sep range is relatively large due to the excited state contamination. C. Axial-vector coupling, axial-vector form factor, and axial radius

Axial-vector coupling
The axial-vector coupling g A = F A (0) has been extensively calculated with LQCD by various groups (see Ref. [15] and references therein). We first show the t dependence of g A extracted from Eq. (22) with Q0 for t sep /a = 10, 12, 14, 16 in Fig. 13. We observe reasonable plateau for all the cases of t sep . Figure 14 shows that our results of g A in all the t sep cases agree with the experimental value, 1.2724(23) [27].
We determine the central value of g A from the com- regarded as the systematic error. The result is summarized in Table V. Our best estimate of g A in this work is where we also include a systematic error stemming from the error of Z A 2 as the third one. This result entirely agrees with the experiment.
Here it may be useful to present the up-and downquark spin component in the nucleon spin, which can be obtained by decomposing the axial-vector coupling into the up-and down-quark contributions: g u A = 0.967(30)(16)(7) and g d A = −0.306(19)(21)(2), in which the first error is statistical one and, the second and third ones are systematic errors due to the excited state contamination and uncertainty of Z A , respectively. Again note that these results are obtained without the disconnected diagram.

Axial-vector form factor and axial radius
We next show the t dependence of F A extracted from Eq. (22) with Q1, Q2, Q3, and Q4 for t sep /a = 10, 12, 14, 16 in Fig. 15. We observe reasonable plateau for all the cases of t sep and finite n. As in the case of G M (q), we do not observe a significant t sep dependence in Fig. 16. We employ two combined values obtained by applying the constant fit to the data in the two ranges of t sep /a = {12, 14, 16} and t sep /a = {14, 16}, which are used for the investigation of q 2 dependence of F A (q 2 ). The two combined values for F A (q 2 ) are summarized in Appendix A.
We plot the q 2 dependence of F A (q 2 ) in Fig. 17, where any strong curvature is not observed in terms of q 2 . We also find that our two combined results show good agreement with each other and both of them are consistent with the experimental values [6] within 1-σ error bars. The isovector axial-vector coupling F A (0) and the axial RMS radius r 2 A are obtained from several types of fitting procedure with the linear, dipole, quadratic forms and the z-expansion method with k max = 3. Note that we employ t cut = 9m 2 π in Eq. (15) for the z-expansion method.
The dipole form fits for two combined data with different choices of t sep range are presented in Fig. 17. The fit results up to q 2 cut = 0.102 GeV 2 well describe our data. As shown in Fig. 17, the fitted curve for t sep /a = {14, 16} appears slightly below their respective data points. This could be due to a poor determination of the covariance matrix in the correlated fit for the highly correlated data among different q 2 points. The results from four types of fit form are compared graphically in Fig. 18 and numerically in Table V. The fit results for both F A (0) and r 2 A show good consistency among all four types of fitting procedures, and they are also in agreement with the experimental values.
Following the analysis in G E (q 2 ) and G M (q 2 ), we take the result of the dipole form with t sep /a = {12, 14, 16} as our best estimate of the central value and the statistical error for r 2 A , and take the maximum difference between the dipole fit with t sep /a = {12, 14, 16} and other fitting with two t sep ranges as the systematic error (see Table V).

D. Induced pseudoscalar form factor
In Fig. 19 we plot the t dependence of F P (q 2 ) extracted from Eq. (22) with Q1, Q2, Q3, and Q4 for t sep /a = 10, 12, 14, 16. One can observe that its depen- IV. Results for the magnetic moments µ and magnetic RMS radius r 2 M for the isovector, proton and neutron channels. In the row of "This work" we present our best estimates, where the first error is statistical and the second one is systematic as explained in the text. Results for the proton and neutron are obtained without the disconnected diagram. Previous work was performed on a 96 4 lattice at mπ = 146 MeV in Ref. [15], where only the statistical errors are presented. dence has slight convex shape for all the cases in contrast to the form factors G E (q 2 ), G M (q 2 ), and F A (q 2 ) discussed above. Inside the fitting range of t, however, the data points are overlapping within 1σ statistical error, so that employing a constant fit to obtain F P (q 2 ) is appropriate. Figure 20 shows the t sep dependence of F P (q 2 ) at the smallest three values of q 2 . We find that F P (q 2 ) clearly increases as t sep increases. This indicates that the significant contribution from the excited states is involved in the F P (q 2 ) case. In fact, the previous work [15] on a 96 4 lattice at the 146 MeV pion with t sep /a = 15 gives F P (q 2 ) close to t sep /a = 14 in our case.
We plot the q 2 dependence of the normalized induced pseudoscalar form factor 2M N F P (q 2 ) for t sep /a = 14 and 16 compared to the previous work [15] in Fig. 21. The colored curve denotes a prediction of the pion-pole dominance (PPD) model with the measured values of m π , M N and the global fit result of F A (q 2 ) in the dipole form: which successfully describes two experimental results of the muon capture [46] and the pion-electroproduction [29]. The t sep dependence of 2M N F P (q 2 ) found in Fig. 20, where the significant change of form factor as t sep increases appears, gives an important hint to explain a discrepancy from experimental values and PPD model Results for the axial-vector coupling gA = FA(0) and axial-vector RMS radius r 2 A . In the row of "This work" we present our best estimates, where the first error is statistical and the second one is systematic as explained in the text. Results for the proton and neutron are obtained without the disconnected diagram. Previous work was performed on a 96 4 lattice at mπ = 146 MeV in Ref. [15], where only the statistical errors are presented. prediction.
For more careful verification of the excited state contamination in the induced pseudoscalar form factor, we need more accurate data at t sep /a = 16 and an additional calculation with at least one more large t sep , which may help to extrapolate F P in the infinite t sep limit. We note that the baryon chiral perturbation theory suggests the aforementioned discrepancy as the contamination of the π-N excited states in the standard plateau method [47,48]. It is also noted that this contamination is recently investigated using a proper projection [49]. More detailed comparison is also interesting for the future work.

V. PSEUDOSCALAR FORM FACTOR AND GOLDBERGER-TREIMAN RELATION
In the previous section, we have found the relatively large excited state contamination in F P (q 2 ) compared to F A (q 2 ), which may one of the reasons for the considerable discrepancy between the LQCD result of 2M N F P (q 2 ) and the experimental values. We expect that the generalized GT relation of Eq. (9), which is associated with the AWT identity, may also suffer from the serious effects of the excited state contamination.
The pseudoscalar form factor G P (q 2 ) is defined by Eq. (10) and extracted from the ratio R 5z P (t, q) of Eq. (23). Figure 22 shows the t dependence of the ratio with |n| 2 = 1, 2, 3, 4 for t sep /a = 10, 12, 14, 16. We observe that the convex shape is much clearer than the F P case and its top value increases for larger t sep . One can see that the pseudoscalar form factor is also strongly affected by the excited state contributions. In Fig. 23, we plot our values of G P (q 2 ) as a function of q 2 , which are obtained by the constant fit of data for both t sep /a = 14 and 16 choosing the same fit range as the other form factors. The stronger curvature appears around q 2 = 0 as t sep increases. This behavior is found to be similar to F P (q 2 ). Although data points of G P (q 2 ) for t sep /a = 14 are comparable with the previous result [15], where t sep /a = 15 was chosen, the magnitude of G P (q 2 ) for t sep /a = 16 becomes about 10 percent larger than that of t sep /a = 14 at all the simulated q 2 points. The t sep dependence of G P (q 2 ) is much more prominent than that of F P (q 2 ).
According to the PPD model or the generalized GT relation associated with the AWT identity, F P (q 2 ) and G P (q 2 ) are supposed to share the same pion-pole structure, i.e., ∝ 1/(q 2 + m 2 π ), at lower q 2 . In the previous work [15], it was indeed observed that the ratio of G P (q 2 )/F P (q 2 ) exhibited a flat q 2 dependence at lower q 2 and was in good agreement with the bare value of the low-energy constant B 0 = m 2 π /(2m) with the simulated pion mass m π and the PCAC quark massm = m PCAC AWTI . In order to test whether this feature holds against the  11. Magnetic moment µ for the isovector (top), proton (middle) and neutron (bottom) channels obtained by the fitting with the linear, dipole, quadratic forms and the zexpansion method for the combined data. Horizontal bands represent the experimental results. Two types of symbols denote the results with two choices of the combined tsep ranges. Results for the proton and neutron channels are obtained without the disconnected diagram.
variation of t sep , we plot the ratios of G P (q 2 )/F P (q 2 ) for all the case of t sep = 10, 12, 14, 16 in Fig. 24. Each ratio of G P (q 2 )/F P (q 2 ) does not depend on q 2 and those are in good agreement with the bare value of the lowenergy constant B 0 as illustrated by the green band. This strongly indicates that the individual effects of the excited state contamination on the G P (q 2 ) and F P (q 2 ) form factors are canceled in the ratio of G P (q 2 )/F P (q 2 ). In order to test the generalized GT relation of Eq. (9), it is convenient to define the following (bare) quark mass as in Ref. [15]: Since the generalized GT relation is an expression of the AWT identity in terms of the nucleon matrix elements, the above quark mass should coincide with the PCAC quark mass m PCAC AWTI extracted from the pseudoscalar matrix elements. In Fig. 25 we plot the quark mass m GT AWTI as a function of q 2 for all the cases of t sep /a = 10, 12, 14, 16. The results do not show any strong q 2 dependence but they are systematically decreased when t sep increases. Compared to the measured m PCAC AWTI [20] from the pion propagator on the same gauge ensembles, m GT AWTI approaches to m PCAC AWTI as t sep /a is increased: m GT AWTI /m PCAC AWTI = 3.0(1) at t sep /a = 12 to 2.3(1) at t sep /a = 16, quoted from the value of Fig. 25 at the minimum q 2 = 0.013 GeV 2 , with only statistical error. This tendency provides a hint to resolve an issue of "distortion of pion-pole structure" discussed in Ref. [15].

VI. SUMMARY AND DISCUSSION
We have calculated the nucleon electric and magnetic form factors, G E (q 2 ) and G M (q 2 ), for not only the isovector channel but also the individual form factors of proton and neutron without the disconnected diagram on a (10.8 fm) 4 lattice at the physical point in 2+1 flavor QCD. We have also measured the axial-vector form fac-   tor F A (q 2 ) together with the axial-vector coupling g A and the axial radius and the induced pseudoscalar form factor F P (q 2 ). Utilizing the optimized all-mode-averaging (AMA) technique with the Wilson-clover fermion, we have investigated the effects of the excited state contamination by varying t sep from 0.85 fm to 1.35 fm with t sep /a = 10, 12, 14, 16 in the plateau method, which has not been studied in the previous work [15]. After elaborate tuning of the sink and source functions, we can obtain clear signal of the nucleon asymptotic state for G E (q 2 ), G M (q 2 ), and F A (q 2 ) without significantly large excited state contamination. Taking account of the uncertainties with the extrapolation onto q 2 = 0 and the excited state contamination, our best estimates for the RMS radii and magnetic moments are obtained as follows: , where the first error is statistical, the second one is systematic error of which the uncertainty of possible excited state contamination and the fit dependence for extrapolation to q 2 = 0 (see Secs. IV A and IV B for the details) are included. They are comparable with the experimental values of r 2 E , which are given by 0.939(6) fm (isovector) and 0.875(6) fm (proton) for the ep scattering, and 0.907(1) fm (isovector) and 0.8409(4) fm (proton) for the µH spectroscopy. The experimental values of r 2 M is given by 0.862 (14) fm (isovector), 0.776(38) fm (proton) and 0.864(9) fm (neutron), and those for the magnetic moment are µ v = 4.70589, µ p = 2.79285 and µ n = −1.91304 quoted from PDG'18 [27]. Although our results for the electric RMS charge radius in the isovector channel seem to favor the experimental result of the µH spectroscopy within 1-σ error, it is still too early to draw any definitive conclusion because of rather large error bars of 4% level. For the proton and neutron channels we leave the inclusion of the disconnected diagram to future work.
In Fig. 26 we compare our results with those obtained by previous LQCD calculations for the isovector channel (see Table VI for the simulation parameters). Those errors are combined with statistical and systematic errors in the quadrature, except that g A for PACS'18 has only statistical error. The electric RMS charge radius given by the PNDME [50] and ETM [10] Collaborations are below the experimental values beyond 1-σ error. In spite of that the other LQCD calculation, e.g., the ETM Collaboration, has also been at the physical pion, their values have differed from our results and experimental values. This may be due to some sorts of finite volume effect on their results [10]. The spatial extent of 10.8 fm in our case allows q 2 = 0.013 GeV 2 as the minimum momentum transfer, which is 6 times smaller than that of the ETM Collaboration [10] 3 , who have employed gauge configurations with N f = 2 twisted mass fermion on a (4.5 fm) 3   Fig. 3 for FA(q 2 ) with four lowest nonzero momentum transfers.
box at the physical pion mass. It means that our simulation on large volume is a strong advantage for the determination of a slope at q 2 = 0 to correctly obtain RMS ETM Collaboration appear in Ref. [51]. Their results of electric RMS radius still has a large discrepancy from experimental values. As we have argued in this paper, this may be due to finite volume effect on their relatively small spatial size, which is up to 6 fm (N f = 2), compared to our calculation on 10.8 fm. FIG. 17. Same as Fig. 4 for the axial-vector form factor FA(q 2 ). Experimental line is obtained by the dipole form with the value of dipole mass [6,28] and gA [27].
radius. Actually the deficit of the charge radius observed in the ETMC results is consistent with the theoretical expectation discussed in Refs. [18,20]. For the magnetic RMS radius and the magnetic moment, our results are consistent with the experimental values, though we find relatively large error bars compared to other LQCD results. Our results indicate that G v M is sensitive to the source-sink separation t sep rather than G v E . In Eq. (30) our large systematic error takes into account such an uncertainty due to the excited state contamination. We find similar stories for the proton and neutron charge radii and the magnetic moment as shown in Fig. 27, even though our result is obtained from only the connected diagram. As discussed in Secs. IV A and IV B we expect that the disconnected contribution may compensate the difference between our results and the experimental val-

ues.
Our best estimates for the axial vector coupling and the axial radius are obtained as in which the first error is statistical and the second one is systematic explained in Secs. IV C 1 and IV C 2. They are comparable with experimental values, g A = 1.2724(23) [27] and r 2 A = 0.67(1) [6,28]. Our result of the axial RMS radius is consistent with the experimental value, while the 7% precision is 4.5 times larger than the experimental one. For the axial-vector coupling, our value is also consistent with the experimental one, though the 2% precision of the former is an order-of-magnitude larger than that of the latter. We also present the results for g u A and g d A neglecting the disconnected contribution, which are in agreement with other LQCD results [11,54] within 2-σ error. In comparison with other LQCD results as shown in Fig. 26, our results show consistency with experimental values within comparable magnitude of error bars to other groups. For the axial radius and the axial-vector coupling our results are significantly improved from the ETMC's results [11]. Here again the spatial lattice size may play a crucial role. Our results on a (10.8 fm) 3 spatial box, which is about 14 times larger than a ∼(4.5 fm) 3 spatial box employed by the other groups e.g., the ETM Collaboration, clearly show consistency of G A (q 2 ) with the Kelly's fit (see Fig. 17), and we observe less excited state contamination effects.
On the other hand, the induced pseudoscalar form factor F P in the axial-vector channel shows clear t sep dependence and considerable deficit from the experimental value in very low q 2 region. Investigation of the generalized GT relation associated with F A , F P , and G P strongly suggests a sizable amount of the excited state contributions to the determination of F A and F P in the plateau method. More dominant excited state contamination compared to the other form factors could be a res- 2M N F P PPD (q) experiment (muon capture) experiment (pion-electroproduction) m π =0.146 GeV, 96c, t sep /a=15 m π =0.135 GeV, 128c, t sep /a=16 m π =0.135 GeV, 128c, t sep /a=14 FIG. 21. q 2 dependence of the normalized induced pseudoscalar form factor 2MN FP (q 2 ). We also plot the results in the previous work [15] for comparison. The colored band shows the function of the pion pole dominance model using the fit result of FA(q 2 ) with the dipole form in Fig. 17. olution of "distortion of pion-pole structure" [15], and it would be solvable once the high-precision data at larger t sep is available. Note that, thanks to our spatial size more than 10 fm, we can first obtain the low-q 2 LQCD data of induced pseudoscalar form factor, which is close to q 2 in MuCap experiment [5], at the physical point. This is the first lattice QCD calculation that succeeds in simultaneously reproducing the experimental values for r 2 E , µ, r 2 M , g A , and r 2 A , and makes an important step for the LQCD calculations to successfully improve its precision to be comparable with the experimental results. In order to assure the reliability of the results, a next step would be further reduction of both statistical and systematic errors such as the cutoff effects and the isospin breaking effects including quark disconnected diagrams.   [52], PNDME [50], ETMC [10], Hasan et al. [12] and PACS [15]. (Bottom) Same as top panels for the axial RMS radius and the axial-vector coupling obtained by CLS-Mainz [53], Green et al. [54], PNDME [9,14], ETMC [11], CalLat [13], and PACS [15]. Those errors are total one combined with statistical and systematic errors in the quadrature, except that gA for PACS '  VI. Summary of simulation parameters in recent LQCD calculations for the nucleon form factors. N f denotes the number of dynamical quark flavors. In the column "Fermion," "TM-Clover" stands for the twisted mass clover-improved Wilson-Dirac operator, "ST-Clover" denotes the stout smeared Wilson-clover fermion, and "HEX-Clover" denotes the HEX smeared Wilson-clover fermion. "MDWF" denotes Möbius domain-wall fermion. In the column of "Method," "R," "S," "TSF," "D," and "FH" stand for the standard plateau (ratio) method, the summation method, the two-state fit method, the derivative method and the method based on the Feynman-Hellmann theorem. Observables