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

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

Ken-Ichi Ishikawa Graduate School of Science, Hiroshima University, Higashi-Hiroshima, Hiroshima 739-8526, Japan RIKEN Advanced Institute for Computational Science, Kobe, Hyogo 650-0047, Japan RIKEN Advanced Institute for Computational Science (AICS) was renamed to RIKEN Center for Computational Science (R-CCS) in April 2018.    Yoshinobu Kuramashi Center for Computational Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan RIKEN Advanced Institute for Computational Science, Kobe, Hyogo 650-0047, Japan    Shoichi Sasaki Department of Physics, Tohoku University, Sendai 980-8578, Japan RIKEN Advanced Institute for Computational Science, Kobe, Hyogo 650-0047, Japan    Natsuki Tsukamoto Department of Physics, Tohoku University, Sendai 980-8578, Japan    Akira Ukawa Center for World Premier International Research Center Initiative (WPI), Japan Society for the Promotion of Science, Tokyo 102-0083, Japan RIKEN Advanced Institute for Computational Science, Kobe, Hyogo 650-0047, Japan    Takeshi Yamazaki Faculty of Pure and Applied Sciences, University of Tsukuba, Tsukuba, Ibaraki, 305-8571, Japan Center for Computational Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan RIKEN Advanced Institute for Computational Science, Kobe, Hyogo 650-0047, Japan
July 14, 2019
###### Abstract

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

###### pacs:
11.15.Ha, 12.38.-t 12.38.Gc
preprint: UTHEP-721, UTCCS-P-113, HUPD-1806

PACS Collaboration

## I Introduction

The nucleon vector and axial elastic form factors are good probes to investigate the internal structure of the nucleon Thomas:2001kw (). 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 high-precision measurements of the proton’s electric charge radius from the muonic hydrogen Lamb shift Pohl:2010zza (); Pohl:2013yb () disagree with well established results of both electron-proton scattering and hydrogen spectroscopy Patrignani:2016xqp (), is currently one of the most intriguing problems in this field Sick:2018fzn (). 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  Czarnecki:2018okw ().

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 axial-vector 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 111 After we complete this work, the new calculation method of the axial-vector coupling is developed based on the Feynman-Hellmann method in Ref. Chang:2018uxx (). Their high-precision result is quite consistent with the experimental value. 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 Bhattacharya:2016zcn (); Green:2014xba (); Alexandrou:2017ypw ().

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 lattice at almost physical point (the pion mass at the simulation point is about 146 MeV) by the PACS Collaboration Ishikawa:2015rho (). There are three reasons for paying special attention to its large spatial volume of .

First of all, the large finite volume effect, represented by the scaling ( denotes the spatial extent of the lattice volume), on measurement of the axial-vector coupling was reported in Ref. Yamazaki:2008py (); Yamazaki:2009zq (). According to their conclusion, one can estimate that spatial sizes of fm () for MeV is necessary for keeping the finite volume effect on 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 () 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 help to more accurately determine the slope of the form factor at , which is associated with the root-mean-square (RMS) radius . In this context, uncertainties in the determination of are sensitive to how large is the spatial extent in physical units. For the case of the electric form factor (), 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 Sick:2018fzn (). The dependence of the Sachs form factors and are roughly described by the dipole shape: with a single parameter called as the dipole mass parameter around 0.84 GeV Thomas:2001kw (). The dipole form assumes that the spatial charge distribution is falling exponentially at large as  Thomas:2001kw (). For a sake of simplicity, we adopt the spatial charge distribution as , where is determined by the normalization condition. The RMS radius is defined in terms of the charge density as

 R2≡4π∫∞0ρ(r)r4dr, (1)

which gives a relation of . We next introduce the truncated RMS radius , where the integral (1) is stopped at , and then calculate the following ratio as a function of  Sick:2018fzn ():

 R(rcut)R=√1−e−X[1+X+12X2+16X3+124X4], (2)

where . To get more than 98% of , the value of must be greater than , which is remarkably large value Sick:2018fzn (). If we take it seriously, this intuitive argument suggests that the spatial extent of  222In a spatial box with periodic boundary conditions, one may consider that a half size of corresponds to ., 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 Ishikawa:2015rho () 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 , the axial-vector coupling , the scalar coupling and the tensor coupling , which are directly accessible from the three-point correlations function at zero momentum transfer in Sec. IV.1. Section IV.2 presents the results of the isovector electric () and magnetic () 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.3) 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 , and then the general form of the matrix element for neutron beta decay is expressed by both the vector and axial-vector transitions:

 ⟨p|V+α(x)+A+α(x)|n⟩=¯up(OVα(q)+OAα(q))uneiq⋅x, (3)

where is the momentum transfer between the neutron () and proton (). The vector () and induced tensor () form factors are introduced for the vector matrix element:

 OVα(q)=γαFV(q2)+σαβqβFT(q2) (4)

and also the axial vector () and induced pseudoscalar () form factors for the axial-vector matrix element:

 OAα(q)=γαγ5FA(q2)+iqαγ5FP(q2). (5)

These matrix elements are here given in the Euclidean metric convention (we have defined 333The sign of all form factors is chosen to be positive. Remark that our definition has the opposite sign relative to that in the Minkowski convention ( and ) adopted in the particle data group Patrignani:2016xqp ().. Thus, denoted in this paper, which stands for Euclidean four-momentum squared, corresponds to the spacelike momentum squared as 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 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

 ⟨p|¯uγαd|n⟩ = ⟨p|¯uγαu−¯dγαd|p⟩ (6) = ⟨p|jemα|p⟩−⟨n|jemα|n⟩ (7)

where . Therefore, this relation gives a connection between the weak vector and induced tensor form factors and the isovector part of electromagnetic nucleon form factors

 Fv1(q2) = FV(q2), (8) Fv2(q2) = 2MNFT(q2), (9)

where denotes the nucleon mass, which is defined as the average of neutron and proton masses. () is given by the isovector combination of the Dirac (Pauli) form factors of the proton and neutron as

 Fv1,2(q2)=Fp1,2(q2)−Fn1,2(q2), (10)

where individual form factors () are defined by

 (11)

Experimental data from elastic electron-nucleon scattering is usually described in terms of the electric and magnetic Sachs form factors, which are related to the Dirac and Pauli form factors:

 GNE(q2)=FN1(q2)−q24M2NFN2(q2), (12)
 GNM(q2)=FN1(q2)+FN2(q2). (13)

Their normalizations at are given by the proton’s (neutron) electric charge and magnetic moment 444For the proton, and , while and for the neutron.. Therefore, one finds

 FV(0)=Fv1(0) = GvE(0)=1, (14) 2MNFT(0)=Fv2(0) = GvM(0)−1=3.70589, (15)

where () represents the isovector combination of the electric (magnetic) form factors of the proton and neutron.

The dependence of the Sachs form factors and are experimentally known that the standard dipole parameterization with (or [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 region. In general, if there is no singularity around for a given form factor , the normalized form factor can be expanded in powers of .

 G(q2)=G(0)(1−16⟨r2⟩q2+1120⟨r4⟩q4+⋯), (16)

where the first coefficient determines the mean squared radius , which is a typical size in the coordinate space. For the dipole form, the root-mean-square (RMS) radius is given as by the dipole mass parameter .

For the axial-vector part of weak processes, the axial-vector form factor at zero momentum transfer, namely, the axial-vector coupling , is precisely determined by measurements of the beta asymmetry in neutron decay. The value of is quoted in the 2016 PDG Patrignani:2016xqp (). The reason why deviates from the corresponding vector coupling is that the axial-vector current is strongly affected by the spontaneous chiral symmetry breaking in the strong interaction Nambu:1961tp (); Nambu:1961fr (). In this sense, this particular quantity allows us to perform a precision test of lattice QCD in the baryon sector.

The dependence of the axial-vector form factor has also been studied in experiments, where the dipole form is a good description for low and moderate momentum transfer  Bernard:2001rs (); Bodek:2007ym (). Recently, the axial mass parameter is much paid attention to by neutrino oscillation studies Bhattacharya:2011ah (); Bhattacharya:2015mpa ().

Although the induced pseudoscalar form factor is less well known experimentally Choi:1993vt (); Gorringe:2002xx (), it is theoretically known that two form factors and in the axial-vector channel are not fully independent. It is because the axial Ward-Takahashi identity: leads to the generalized Goldberger-Treiman (GT) relation among three form factors Weisberger:1966ip (); Sasaki:2007gw ():

 2MNFA(q2)−q2FP(q2)=2^mGP(q2), (17)

where is a degenerate up and down quark mass and the pseudoscalar nucleon form factor is defined in the pseudoscalar nucleon matrix element

 ⟨p|P+(x)|n⟩=¯up(γ5GP(q2))uneiq⋅x (18)

with a local pseudoscalar density, . Therefore, the dependences of three form factors, , and are constrained by Eq. (17) as a consequence of the axial Ward-Takahashi identity. Therefore, three form factors, , and 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 Bhattacharya:2016zcn (); Bhattacharya:2013ehc (); Green:2014xba (); Yamazaki:2008py (); Yamazaki:2009zq (); Capitani:2015sba (); Alexandrou:2013joa (); Alexandrou:2017ypw (); Ohta:2017gzg (). 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 Yamazaki:2015vjn (); Kuramashi:2016lql (); Tsukamoto:2017fnm (). In this paper, we present an update of our previous studies Yamazaki:2015vjn (); Kuramashi:2016lql (); Tsukamoto:2017fnm (), 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 -improved Wilson-clover quark action and the Iwasaki gauge action Iwasaki:2011np () on a lattice at , which corresponds to a lattice cutoff of GeV ( fm) Ishikawa:2015rho (). Periodic boundary conditions are used for the gauge and quark fields in all four directions. The stout smearing parameter is set to  Morningstar:2003gk (). The improvement coefficient, , is determined nonperturbatively by the Schrödinger functional (SF) scheme Taniguchi:2012kk (). The improved factor for the axial-vector current becomes very small at the nonperturbative value of and is consistent with zero within the statistical error Taniguchi:2012kk (). Therefore, we do not consider improvement of quark bilinear current. The hopping parameters for the light sea quarks = (0.126117, 0.124790) are chosen to be near the physical point. For the first time, a simulated pion mass reaches down to MeV on a large spatial volume of in 2+1 flavor QCD.

The degenerated up-down quarks are simulated with the DDHMC algorithm Luscher:2003vf (); Luscher:2005rx () using the even-odd preconditioning and the twofold mass preconditioning Hasenbusch:2001ne (); Hasenbusch:2002ai (). The strange quark is simulated with the UVPHMC algorithm deForcrand:1996ck (); Frezzotti:1997ym (); Frezzotti:1998eu (); Frezzotti:1998yp (); Aoki:2001pt (); Ishikawa:2006pb () where the action is UV-filtered deForcrand:1998sv (); Alexandrou:1999ii () 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. Ishikawa:2015rho (). 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. Yamazaki:2015vjn (); Kuramashi:2016lql ().

### iii.1 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

 (19)

where the superscript denotes a transposition and is the charge conjugation matrix defined as . The indices and label color and flavor, respectively. The subscript 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 (denoted as ) and the exponential smearing function: with and (denoted as ). For a simplicity, is chosen.

We then construct two types of the two-point function with the exponential smearing source (the source-time location denoted as ) as

 CXS(tsink−tsrc,p)=14Tr{P+⟨NX(tsink,p)¯NS(tsrc,−p)⟩}, (20)

where (local) or (smear) stands for a type of smearing at the sink operator (the sink-time location denoted as ). A projection operator can eliminate contributions from the opposite-parity state for  Sasaki:2001nf (); Sasaki:2005ug (). In this study, for nonzero spatial momentum, we use the nine lowest momenta: with a vector of integers . All choices of are listed in Table 2.

### iii.2 Nucleon spectra and dispersion relation

In Fig. 1 we plot the nucleon effective mass with 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 . On the other hand, the signal becomes noisier for the smear-smear case. We also measure the nucleon energies 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 obtained from the single exponential fit are tabulated in Table 3.

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 , while the horizontal axis is the momentum squared given by the lattice momentum in lattice units. As can be seen in Fig. 2, the relativistic continuum dispersion relation is well satisfied up to .

### iii.3 Three-point correlation functions for nucleon form factors

The nucleon form factors are extracted from the three-point correlation functions consisting of the nucleon source and sink operators with a given local current ( and ) located at the time-slice :

 (21)

where denotes an appropriate projection operator to extract the form factors and represents the three-dimensional momentum transfer. The local current is given by where is a Dirac matrix appropriate for the channel .

We then calculate the following ratio constructed from the three-point correlation function with nucleon two-point functions :

 RkO,α(t,p′,p)=CPkO,α(t,p′,p)CSS(tsink−tsrc,p′)√CLS(tsink−t,p)CSS(t−tsrc,p′)CLS(tsink−tsrc,p′)CLS(tsink−t,p′)CSS(t−tsrc,p)CLS(tsink−tsrc,p), (22)

which is a function of the current operator insertion time at the given values of momenta and for the initial and final states of the nucleon.

We consider only the rest frame of the final state with , which leads to for the squared four-momentum transfer. Nucleon energy simply abbreviated as , hereafter. In this kinematics, is represented by a simple notation . We calculate only the connected diagrams to concentrate on the isovector part of the nucleon form factors. The current insertion points 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 (denoted as ) corresponds to 1.26 fm. In previous calculations Bhattacharya:2013ehc (); Green:2014xba (); Capitani:2015sba (); Alexandrou:2017ypw (); Yamazaki:2009zq (); Ohta:2017gzg (), 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 Martinelli:1988rr (); Sasaki:2003jh (). This method requires the sequential quark propagator for each choice of a projection operator regardless of the types of the current . To minimize the numerical cost, we consider two types of the projection operator and in this study.

The ratio (22) with appropriate combinations of the projection operator () and the component of the current gives the following asymptotic values Sasaki:2007gw ():

 RtV,4(t,q)→√EN+MN2ENGvE(q2) (23)

and

 R5zV,i(t,q)→−iεij3qj√2EN(EN+MN)GvM(q2) (24)

for the vector current () in the limit when the Euclidean time separation between all operators is large, with fixed and . Similarly, we get

 (25)

for the axial-vector current (), and

 R5zP(t,q)→iq3√2EN(EN+MN)GP(q2) (26)

for the pseudoscalar (), respectively. Here, we recall that the lattice operators receive finite renormalizations relative to their continuum counterparts in general. Thus, all above form factors , , , and are not renormalized yet. Their renormalized values require the renormalization factor (), which is defined through the renormalization of the quark bilinear currents

 [¯uΓOd]ren=ZO[¯uΓOd]lattice. (27)

The renormalization factors , and 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 Ishikawa:2015fzw (). The isovector electric and magnetic form factors are simply abbreviated as and , 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 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 with as described in Table 2. 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.

### iv.1 Nucleon three-point function with zero momentum transfer

For zero momentum transfer , the ratio (22) becomes

 RkO,α(t,0)=CPkO,α(t,0)CSS(tsink−tsrc,0), (28)

which vanishes unless , and with  Sasaki:2003jh (). The non-vanishing ratio gives an asymptotic plateau corresponding to the bare value of the coupling relevant for the channel in the middle region between the source and sink points, when the condition is satisfied.

With our choice of the projection operators , we can determine 4 different couplings: the vector coupling , the axial-vector coupling , the scalar coupling and the tensor coupling , from the following four ratios,

 RtV,4(t,0) → GE(0)=FV(0)=gV, (29) R5zA,3(t,0) → FA(0)=gA, (30) RtS(t,0) → gS, (31) R5zT,12(t,0) → gT, (32)

whose values are not yet renormalized.

In Fig. 3, we plot the above four ratios as a function of the current insertion time slice 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 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 . 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), -dependence of is symmetric between the source () and sink () points. Although this symmetric behavior with respect to 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 , where an asymptotic plateau behavior reveals in the ratio of with our choice of a source-sink separation.

In each plot, a solid line indicates the average value over range of and a shaded band displays one standard deviations estimated by the jackknife method. The obtained values of the bare couplings , which are not yet renormalized, are summarized in Table 4.

We next evaluate the renormalization factor of the local vector current through a property of for the renormalized value of the vector charge under the exact isospin symmetry. Therefore, the value of can be evaluated by since .

Figure 4 plots the result of , which shows a good plateau in the time range of . The constant fit result with one-standard-error band denoted by three solid lines shows a good consistency with the value of (red line) obtained by the Schrödinger functional (SF) scheme at the vanishing PCAC quark mass Ishikawa:2015fzw (). 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 with the SF scheme for comparison in the same figure. The difference between and 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 obtained with the SF scheme Ishikawa:2015fzw () shown in Fig. 4. We then plot the renormalized value of the axial charge as a function of the current insertion time slice in Fig. 5. Three solid lines represent the fit result with one-standard error band in the time region of , while the uncertainty in determination of the value of is take into account by shaded band. We finally obtain the renormalized value of the axial charge

 grenA=1.163(75)stat(14)ZA, (33)

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 -dependence of 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 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-mode-averaging (AMA) Blum:2012uh () should be useful and efficient to improve the statistical accuracy of with the limited number of configurations. We intend to examine this possibility in future research.

### iv.2 Results of nucleon form factors in the vector channel

#### iv.2.1 Electric Ge and magnetic Gm form factors

The isovector electric form factor and magnetic form factor 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 for the isovector electric form factor (left) and magnetic form factor (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 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 denoted by the gray shaded area as same as the cases of the coupling ().

In Table 5, we compile the results of and form factors together with the results of and form factors, which are evaluated by linear combinations of and with appropriate coefficients determined by Eqs. (23) and (24) with measured values of the momentum squared and the nucleon mass . The dependence of (left panel) and (right panel) are plotted in Fig. 7, respectively. For the finite three momentum , we use the nine lowest nonzero momenta: with . The lowest nonzero momentum transfer in present work reaches the value of 0.024(1) [], which is much smaller than even at MeV. In each panel, we also plot Kelly’s fit Kelly:2004hm () (red dashed curves) as their experimental data. The simulated electric form factor is fairly consistent with the experiment, especially at low . The magnetic form factor 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 would hinder us from the precise determination of the magnetic RMS radius. The electric (magnetic) RMS radius, (), can be evaluated from the slope of the respective form factor at zero momentum transfer

 ⟨r2l⟩ = −6Gl(0)dGl(q2)dq2∣∣∣q2=0 (34)

with (). Recall that at zero momentum transfer, whose value corresponds to the isovector magnetic moment , cannot be directly measured for the kinematical reason Sasaki:2007gw (). Therefore, -extrapolation is necessary to evaluate the value of .

In principle, low 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 . In this sense, a large spatial size is required for precise determination of the RMS radii ( and ) and magnetic moment ().

The determination of the slope (or -extrapolation) is highly sensitive to how we fit the -dependence of the form factors. In the previous studies Bhattacharya:2013ehc (); Green:2014xba (); Yamazaki:2009zq (); Capitani:2015sba (); Alexandrou:2013joa (), the dipole form , and the polynomial form have been often adopted for and . 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 ). In order to reduce systematic errors associated with the interpolation or extrapolation of the form factor in momentum transfer, we use the model-independent -expansion method Boyd:1995cf (); Hill:2010yb ().

#### iv.2.2 Analysis with z-expansion method

Let us recall the analytic structure of the form factors in the complex -plane. There is a non-analytic region for due to threshold of two or more particles, e.g. continuum. Hence the -expansion, , does not converge if where is a branch point associated with the pion pair creation for or  Hill:2010yb ().

The -expansion (denoted as z-Exp) makes use of a conformal mapping from to a new variable . Since this transformation makes the analytic domain mapped inside a unit-circle in the -plane, the form factors can be described by a convergent Taylor series in terms of :

 G(z)=kmax∑k=0ckz(q2)k (35)

with

 z(q2)=√tcut+q2−√tcut√tcut+q2+√tcut, (36)

where truncates an infinite series expansion in  555We use the singular value decomposition (SVD) algorithm to solve the least squares problem for high degree polynomials.. Here, , where corresponds to the simulated pion mass, is chosen for the vector channel ( or ), while , which is associated to continuum, will be later chosen for the axial-vector channel (Bhattacharya:2011ah ().

We here note that must ensure that terms become numerically negligible for for a model-independent fit. Although is expected for sufficiently large , the range of possible values of is limited by the condition (7) for the () 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 .

Figure 8 shows the fit results for all ten data points of with and 8 in the -expansion. Similarly, we also plot the fit results for all nine data points of with and 7 in Fig. 9. At a glance, one can see from Fig. 8 and Fig. 9 that the curvature becomes smaller in variable than variable for both cases of and . This indicates that a power series in terms of is clearly more efficient than that of . It is also observed that both fit results are not sensitive to the choice of . This suggests that the -expansion gives a rapid convergence series for both cases. Indeed, as shown in Fig. 10, the values of are insensitive to the choice of if . For , the ratios of with become zero within the statistical error, while the ratios of for reach a convergence value less than unity at or 6. Thus, the value of 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 . For these reasons, (7) is hereafter chosen for the () form factor in the z-Exp method.

In Fig. 11, we show the fit results of (left panel) and (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 () with . However, if one takes a closer look at low , the fit curve given by the -expansion fit goes nicely through the data points in low region. This implies that the z-Exp fit is relatively insensitive on the higher data points as was expected. The RMS radius is determined to be (z-Exp fit), (dipole fit) and (quadratic fit), while the coefficients of , and for correspond to the value of the magnetic moment, respectively. All results obtained by various fits are compiled in Table 6.

In Fig. 12, we compare the results of (left panel) and (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 scattering (upper) and muonic hydrogen (-) 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 and from the z-Exp fit is most consistent with respective experiment. As our final results, we then quote the value of the (isovector) electric RMS radius:

 √⟨r2E⟩=0.915±0.099fm (37)

and the value of the (isovector) magnetic moment:

 μv=4.81±0.79, (38)

which are obtained from the z-Exp fits. The former is consistent with both the experimental values from 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.

#### iv.2.3 Comparison with previous results

We finally compare our results of and with recent lattice results from LHPC Green:2014xba (), PNDME Bhattacharya:2013ehc (), the Mainz group Capitani:2015sba (), ETMC (denoted as ETMC’13 Alexandrou:2013joa () and ETMC’17 Alexandrou:2017ypw ()) and RBC-UKQCD (denoted as RBC-UKQCD’09 Yamazaki:2009zq () and RBC-UKQCD’17 Ohta:2017gzg ()) as shown in Fig. 13. See also Table 8 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 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 region. Raw data of and are available in tables of Ref. Green:2014xba (). Therefore, it is worth comparing our results of and 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 (left panel) and Pauli (right panel) form factors as a function of momentum squared . In each panel, the dashed curves correspond to Kelly’s fit Kelly:2004hm () as their experimental data. In the LHPC results, they use two types of kinematics regarding momentum for the final state of the nucleon state when constructing the three-point correlation functions (22). The smallest momentum transfer [] in the LHPC calculation is given by the choice of , while we consider only the rest frame of the final state with .

Our results and the LHPC results are consistent with each other in both and 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 . It is clear that -dependence of the form factors at low are much constrained by our results that contain the smallest nonzero data point. We primary focus on the values of and since the data of the form factor exhibits large statistical fluctuations in the low