Nucleon mass and isovector couplings in 2+1-flavor dynamical domain-wall lattice QCD near physical mass

Michael Abramczyk, Thomas Blum, 2 Taku Izubuchi (出渕卓), 2 Chulwoo Jung, 2 Meifeng Lin, Andrew Lytle, Shigemi Ohta (太田滋生), 7, 2 and Eigo Shintani (新谷栄悟) (for RBC and UKQCD collaborations) Physics Department, University of Connecticut, Storrs, CT 06269-3046 RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973 Physics Department, Brookhaven National Laboratory, Upton, NY 11973 Computational Science Initiative, Brookhaven National Laboratory, Upton, NY 11973 INFN, Sezione di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Roma RM, Italy Theory Center, The High Energy Accelerator Research Organization (KEK), Oho 1-1, Tsukuba, Ibaraki 3050801, Japan Deaprtment of Particle and Nuclear Physics, The Graduate University of Advanced Studies (SOKENDAI), Hayama, Kanagawa 2400193, Japan RIKEN Advanced Institute of Computational Sciences (AICS), Kobe, Hyogo 650-0047, Japan (Dated: November 12, 2019) KEK-TH-2167, RBRC-1320

In our earlier works calculated with degenerate up-and down-quark mass set at considerably heavier than physical values [4][5][6], we observed the vector-current form factors behaving reasonably well in trending toward experiments: both Dirac and Pauli mean-squared charge radii and the isovector anomalous magnetic moment appeared to linearly depend on the pion mass squared. The radii extrapolated to the physical mass undershoot the experimental value by about 25 percent [17,18]. It would have been interesting to see if the present calculation confirmed this earlier trend, or if it could help resolve the discrepancy seen between muon Lamb shift experiment [17] and older electron scattering ones. However, our current numerical precision from relatively small statistics and large momentum transfer is yet to be competitive with Lamb shift experiments [17,18] which now are prevailing [19,20]. So we would like to defer reporting our form factors at finite momentum transfers until a future date when we will have better precision.
In our earlier calculation of axial-vector current form factors, we saw a significant deficit in the calculated axialvector charge, g A , and form factors in general appear more susceptible to finite-size effects than the vectorcurrent ones [4,5]. We find in our present calculations at lighter pion mass that this deficit persists, and we investigate potential causes for this in detail in Section V.
In contrast, the isovector tensor coupling showed interesting downward departure at the lightest mass to about 1.0 away from the flat higher-mass values of about 1.1 [6]. Whether this trending continues in our present calculations at considerably lighter mass is obviously an interesting question. The tensor and scalar couplings are also relevant to the search for new physics beyond the standard model such as neutron electric dipole moment [21,22]. Some preliminary analyses of these nucleon observables have been reported at recent Lattice conferences [23][24][25][26][27]. In addition, the LHP collaboration also calculated some nucleon observables [28] using a RBC+UKQCD 2+1flavor dynamical DWF ensemble [9].

II. NUMERICAL METHOD
In this paper we concentrate on our results for mass and four isovector couplings of nucleon: vector and axialvector charges and tensor and scalar couplings. Though we summarize their definitions and computational methods later in this section, we refer the readers to our earlier publication [5] for the full details.
The two 2+1-flavor dynamical domain wall fermion gauge field ensembles we use in this work [11,12] were generated jointly by the RBC and UKQCD collaborations on 32 3 × 64 four-dimensional volume and L s = 32 in the fifth dimension with Iwasaki × DSDR gauge action at the gauge coupling of β = 1.75 and 2+1-flavors of dynamical DWF quarks of the domain-wall height of 1.8, strange-quark mass of 0.045, and degenerate up-and down-quark mass of 0.0042 and 0.001 in lattice units. These parameters result in the inverse lattice spacing, a −1 , of 1.378(7) GeV and the DWF residual quark mass of 0.001842 (7). Note that the inverse lattice spacing has been slightly revised from the original [10,11] using the global chiral and continuum fits in conjunction with new physical-mass ensemble sets [12] with Möbius DWF quarks. Thus the heavier of the two ensembles corresponds to the pion mass, m π , of 249.4(3) MeV and spatial lattice extent L of m π L = 5.79 (6), and the lighter to 172.3(3) MeV and 4.00(6), respectively. Our measurement calculations were made using 165 configurations between the molecular-dynamics (MD) trajectory 608 and 1920 with an 8-trajectory interval for the former, and using 39 configurations between 748 and 1420 with a 16trajectory interval for the latter.
We refer to our earlier publications [3,4] for the details of two-and three-point correlation functions for the nucleon. A conventional nucleon operator, , with color indices a, b, and c, quark flavors u for up and d for down, and charge conjugation operator C = γ 4 γ 2 is used. Additionally, gauge-invariant Gaussian smearing [29,30] is applied to suppress excited-state contamination.
Isospin symmetry is enforced for the up and down quarks, and we calculate only proton matrix elements of the third isospin component of the vector and axialvector currents, All quark-disconnected diagrams cancel in these matrix elements.
We use a source-sink separation of 9 lattice units, or about 1.3 fm. This is sufficiently long for the observables reported in this paper to be free of excited-state contamination, as is demonstrated below with augmentative calculations with source-sink separation of 7 lattice units [24,25]. For the two point correlation functions we use the same Gaussian-smeared sources and point or Gaussian-smeared sinks: We will refer to the former as G-L and the latter as G-G respectively. We also calculate the tensor coupling, g T = 1 δq , and scalar coupling, g S , the same way. We used a conventional measurement strategy for the former with seven source-sink pairs for each configuration, and an "all-mode-averaging (AMA)" strategy [31] for the latter with 112 sloppy and four precise solves for each configuration. Table I summarizes the nucleon energies obtained from correlated, single-exponential fits to the G-G two point correlation function. To improve statistics the correlation function is averaged over forward and backward propagating nucleon and anti-nucleon states. The fit ranges were chosen after inspecting the effective masses and taking the shortest distance from the source consistent with an acceptable χ 2 value from the fit. Increasing this minimum distance by a single time unit does not change the energy beyond the statistical error. From these we estimate the nucleon mass for the present ensembles as summarized in Tab. II and Fig. 1. Whereas our prior calculations at heavier pion masses found nucleon mass extrapolating to values much higher than experiment, with pion masses now sufficiently close to physical that the new data extrapolates, linearly in terms of pion mass squared, m 2 π , to a value m N = 0.6894(7)a −1 = 0.6894(7) × 1.378(7)GeV = 0.950(5) GeV. This extrapolation is only slightly more than two standard deviations away from the average of proton and neutron mass,

III. NUCLEON MASS
FIG. 1. Estimated nucleon mass, mN , plotted against estimated pion mass squared, m 2 π , of the present ensembles (ID32, cyan) and the two lightest of ref. [5,6] (I24, magenta) with the new and more accurate estimate for the inverse lattice spacing [12]. The present results linearly extrapolate to the experiment (2) within the statistical error. 0.938918747(6) GeV [19]. The slope in this linear extrapolation is steeper than that observed in our earlier calculations with heavier mass: the result constrains nonlinear dependence of nucleon mass on pion mass squared.

IV. VECTOR CHARGE
Signals for the isovector vector charge, g V , are shown in Fig.2 for both m π = 172.3(3) and 249.4(3) MeV ensembles. Robust time-independent plateaux are seen. In the following we average values in the range 3 ≤ t op ≤ 6, and find 1.450(4) for the heavier and 1.447(9) for the lighter ensemble, respectively. The values compare with the inverse of the vector-current renormalization, Z V , computed in the meson sector [11], 0.664(5) −1 = 1.506 (11) and 0.669(4) −1 = 1.495 (9). Alternatively we linearly extrapolate the two calculated vector-charge values of 1.450(4) at m f a = 0.0042 and 1.447(9) at 0.001 to the chiral limit, m f a = −m residual a = −0.0018427 to obtain a value 1.474 (11). This is to be compared with the inverse of the meson-sector vector-current renormalization, Z V , in the chiral limit, Z −1 V = 0.673(8) −1 = 1.49(2). Thus the nucleon vector-charge in the chiral limit renormalizes to unity within statistical errors: Z V g V = 0.992 (14). Had the renormalized vector charge, Z V g V , deviated from unity, a likely cause would be excited-state contamination, excitedstates excitedstate|g V |0 , through violation of vector-current conservation at O(a 2 ) in lattice spacing a. That our vector charge renormalizes to unity within a couple of percent statistical error constrains such excited-state contamination. As is shown in Fig. 3, signals for isovector vector charge, g V , from source-sink separation of 7 and 9 lattice units agree to 0.6 percent, well within statistical errors. If we assume the first excited state is the ground-state nucleon plus a pion of mass m π a = 0.1249(2) with a unit of lattice momentum, 2π/L = 0.1963, then it decays as exp(−2 × (0.1249 + 0.1963)) = 0.5260 in two lattice units of time from 7 to 9. So the relative amplitude of this state in the source times the O(a 2 ) mixing matrix element cannot exceed approximately one percent.

V. AXIAL-VECTOR CHARGE
As is the case for the vector current, we use the localcurrent definition for the axial-vector current. Because the two local currents are connected by a chiral rotation, they share a common renormalization, Z A = Z V , that relates them with the corresponding conserved global currents, up to O(a 2 ) discretization. This is an advantage of the DWF scheme. Thus for the axial-vector charge, g A , it is better to look at its ratio, g A /g V , to the vector charge, for precision, as is demonstrated in Fig. 4. The calculated value of the ratio, g A /g V , underestimates the experimental value of 1.2732(23) [19] by about 10% and does not depend much on the pion mass, m π , in the range from about 418.8(1.2) MeV down to 172.3(3) MeV from four recent RBC+UKQCD 2+1-flavor dynamical DWF ensembles [9][10][11]. The result appears to have been confirmed by several other major collaborations [21,[32][33][34][35] using different actions but with similar lattice spacings and quark masses, though extrapolations to physical mass seem to differ. It is especially important for calculations with Wilson-fermion quarks [21,32] to remove the O(a) systematic errors at the linear order in the lattice spacing, a [33].
We are obviously suffering from some systematics that make our calculations undershoot the experimental value of g A /g V = 1.2732(23) [19]. Indeed we see possible signs of inefficient sampling: First we observe an unusually long-range autocorrelation when we divide the lightest ensemble at m π = 172.3(3) MeV into two halves, earlier and later, in hybrid MD time, as in Fig. 6. Indeed when we further divide into four consecutive quarters in MD time, the axial-vector charge starts at a value consistent with experiment but monotonically decreases to a value below unity, as in Fig. 7.
Importantly, we also note that no such undersampling is seen in any other isovector observables we have looked at, including the vector charge, g V , quark momentum fraction, x u−d and quark helicity fraction, x ∆u−∆d , and that blocked-jackknife analyses with block size of 2 showed strong correlation of two successive gauge configurations for g A and g A /g V but not for the other observables.
A similar but weaker sign of unusually long-range autocorrelation can be seen in the lightest of our earlier ensembles [9] at m π = 331.3(1.4) MeV when we divide it into four consecutive quarters in hybrid MD time as in Fig. 8. However no such sign of undersampling is seen in the other two ensembles, the present one at m π = 249.4(3) MeV and another at 418.8(1.2) MeV from the earlier work [9]. In other words, the strongest sign of undersampling is seen at the smallest finite-size scaling parameter, m π L ∼ 4.00 (6). Another weaker indication of a finite-volume effect is seen at the second smallest but not the second lightest at m π L ∼ 4.569 (15). No effect is seen for larger values at m π L ∼ 5.79(6) or 5.813 (12). This of course does not prove that the problem is caused by the finite lattice spatial volume, but suggests so.  That there is a long-range autocorrelation in this observable is corroborated by blocked-jackknife analysis with block sizes of 2, 3, and 4, as is summarized in Tab. III: the statistical error of the axial-vector charge keeps growing to at least beyond a block size of 3 while those for the other observables do not grow at all except perhaps for tensor coupling which nonetheless stops growing earlier.
If an observable appears to have long-range autocorrelation, it would be interesting to look at its correla-  tion with the topology of the gauge configurations. We explored this possibility by plotting jackknife samples against the topological charge (see Fig. 9), and did not find correlation. We can also look at whether our low-mode deflation affected this, though the available information is limited to about half of the configurations of what we are presenting from the 172-MeV ensemble (see Fig. 10.) Albeit with this limitation we do not find any correlation either: average of the lowest 100 eigenmodes does not differ between the two halves. Some difference emerges as we go to higher eigenmodes but does not appear to be significant.
On the other hand, as pointed out earlier in this paper, similar long-range autocorrelation was seen in the 331.3(1.4)-MeV ensemble [24] that is at the second  It may be also instructive to remember earlier phenomenological analyses such as one done using the MIT bag model that estimates g A /g V = 1.09 without pion [36], and another by the Skyrmion model that gives only a conditionally convergent result of 0.61 that is strongly dependent on pion geometry [37].
To explore such spatial dependence arising from pion geometry, we divided the AMA samples into two spatial halves such as 0 ≤ x < L/2 and L/2 ≤ x < L for each of the three spatial directions in order to check if there is any uneven spatial distribution (see Fig. 11.) We found that the calculation fluctuates in space. Larger spatial volume would stabilize the calculation better.

VI. TENSOR AND SCALAR COUPLINGS
Plateau signals for the bare isovector tensor, g T = 1 δu−δd , and the scalar coupling, g S , are presented in Fig. 12. The tensor-coupling signals are very clean and do not show any mass dependence. As was mentioned earlier, this observable shows weaker but still relevant signs of long-lasting autocorrelation similar to that of the axial-vector charge in the lighter, 172.3(3)-MeV, ensemble [25]. Yet the agreement with the heavier ensemble where there is no such autocorrelation reassures this is less problematic in the tensor coupling than in the axialvector charge. The scalar plateaux are also well defined albeit with larger statistical errors. No mass dependence can be seen here either. Bilinear operators are renormalized using the Rome-Southampton nonperturbative renormalization (NPR) method [13]. This method allows determination of lattice matching factors (Z-factors) to regularization independent (RI) schemes using lattice simulations directly, and without recourse to lattice perturbation theory. The validity of the method requires that the renormalization scale µ is sufficiently separated from the QCD scale and the scale of the lattice cutoff: (3) A perturbative calculation in the continuum is then conventionally used to convert the RI-scheme results to the MS scheme. The relation between bare (lattice) operators and their renormalized counterparts is Here O is one of S, V , or T for scalar, vector, and tensor bilinears respectively. We compute the amputated Landau-gauge fixed Green's functions of the operators of interest between external quark states with momentum p 1 and p 2 . Let these quantities be denoted Λ Γ (p 1 , p 2 ), where Γ is one of {S = 1, V = γ µ , T = σ µν }. Here we use the RI/SMOM scheme defined in [16], for which the momenta satisfy p 2 1 = p 2 2 = (p 1 − p 2 ) 2 . This differs from the original RI/MOM scheme proposal [13] which instead takes p 1 = p 2 . This improved kinematic setup maintains a single renormalization scale, but ensures momentum flows through the vertex thereby suppressing unwanted infrared effects. It was also found in [16] that the RI/SMOM to MS matching factors have smaller O(α s ) coefficients.
There are two principal SMOM variants defined in [16], called RI/SMOM and RI/SMOM γµ , which differ in the definition of wavefunction renormalization. The bilinear Tr σ µν Λ µν T,R (p 1 , p 2 ) = 1 on the amputated Green's functions in the chiral limit.
Here the traces are over both spin and color indices, and the normalization factors are such that Eqs. (5) are satisfied by the bare operators in the free theory.
The propagators used to construct Λ Γ (p 1 , p 2 ) are computed using momentum sources, which results in very low statistical noise using a modest number of configurations. Additionally, twisted boundary conditions on the quark fields allow us to continuously vary the magnitude of p while keeping the orientation fixed, resulting in smooth data as a function of renormalization scale.
Results for the RI/SMOM and RI/SMOM γµ intermediate schemes are presented in Tabs. IV, V, and VI. The wavefunction renormalization factors have been removed using Z V /Z q determined from (5), where Z V relates the local 4-d current to the conserved 5-d current. These results are then converted to the MS scheme using the two-loop perturbative expressions calculated in [38,39].
Taking the average of results from the two intermediate schemes after conversion to modified MS, and taking the full difference as an estimate of the systematic, we find for the renormalization factors of the scalar and tensor couplings: Z S (MS, 2GeV) = 0.642(8) (22) and Z T (MS, 2GeV) = 0.731 (8) (24) From these, we obtain our estimates for the renormalized isovector tensor and scalar couplings as presented in Table VII. Neither is dependent on mass.  .177(1) 1.166(1) 1.155(1) 1.146(1) 1.138(1) 1.131(1) 1.125(1) The tensor coupling is in good agreement with a value, about 1.0, obtained at the lightest mass in our previous calculations [6] and also with later calculations by others [40,41]. Its errors are dominated by a schemedependent systematics in non-perturbative renormalization, at about five percent, due mainly from the relatively low lattice cut off.
The scalar coupling, though noisier, is in broad agreement with other later calculations [40,41]. The scalar errors are still dominated by statistical noise, but will eventually encounter the same scheme-dependent systematics in non-perturbative renormalization.

VII. CONCLUSIONS
The nucleon mass calculated in the two present ensembles extrapolate linearly in pion mass squared, m 2 π , to a value m N = 0.950(5) GeV at the physical point. This is to be compared with the average of proton and neutron mass experimental values, 0.938918747(6) GeV [19]. The slope in this linear extrapolation is steeper than that observed in our earlier calculations with heavier mass: the result constrains non-linear dependence of nucleon mass on pion mass squared.
The isovector vector charge renormalizes to unity in the chiral limit. This narrowly constrains excited-state contamination in the Gaussian smearing.
The ratio of the isovector axial-vector to vector charge shows a deficit of about ten percent. This is in agree-ment with some other major lattice numerical calculations [21,[32][33][34][35] using different actions but with similar lattice spacings and quark masses. The origin of this deficit is still to be understood.
We obtained good signals for isovector tensor coupling. It does not depend on mass and extrapolates to 1.05 (5) at physical mass with MS 2-GeV renormalization. This is in agreement with the value obtained for the lightest pion mass of 331.3(1.4) MeV in our earlier work [6] and also with calculations by others [40,41].
Isovector scalar coupling is noisier but again does not show mass dependence, and is in agreement with other calculations [40,41].

ACKNOWLEGMENT
The ensembles were generated using four QCDOC computers of Columbia University, Ediburgh University, RIKEN-BNL Research Center (RBRC) and USQCD collaboration at Brookhaven National Laboratory, and a Bluegene/P computer of Argonne Leadership Class Facility (ALCF) of Argonne National Laboratory provided under the INCITE Program of US DOE. Calculations of nucleon observables were done using RIKEN Integrated Cluster of Clusters (RICC) at RIKEN, Wako, and various Teragrid and XSEDE clusters of US NSF. SO was partially supported by Japan Society for the Promotion of Sciences, Kakenhi grant 15K05064.   (8)