Isovector Charges of the Nucleon from 2+1+1-flavor Lattice QCD

We present high statistics results for the isovector charges $g^{u-d}_A$, $g^{u-d}_S$ and $g^{u-d}_T$ of the nucleon. Calculations were carried out on eleven ensembles of gauge configurations generated by the MILC collaboration using highly improved staggered quarks (HISQ) action with 2+1+1 dynamical flavors. These ensembles span four lattice spacings $a \approx$ 0.06, 0.09, 0.12 and 0.15 fm and light-quark masses corresponding to $M_\pi \approx$ 135, 225 and 315 MeV. Excited-state contamination in the nucleon 3-point correlation functions is controlled by including up to three-states in the spectral decomposition. Remaining systematic uncertainties associated with lattice discretization, lattice volume and light-quark masses are controlled using a simultaneous fit in these three variables. Our final estimates of the isovector charges in the $\overline{\text{MS}}$ scheme at 2 GeV are $g_A^{u-d} = 1.218(25)(30)$, $g_S^{u-d} = 1.022(80)(60) $ and $g_T^{u-d} = 0.989(32)(10)$. The first error includes statistical and all systematic uncertainties except that due to the extrapolation ansatz, which is given by the second error estimate. We provide a detailed comparison with the recent result of $g_A^{u-d} = 1.271(13)$ by the CalLat collaboration and argue that our error estimate is more realistic. Combining our estimate for $g_S^{u-d}$ with the difference of light quarks masses $(m_d-m_u)^{\rm QCD}=2.572(66)$ MeV given by the MILC/Fermilab/TUMQCD collaboration for 2+1+1-flavor theory, we obtain $(M_N-M_P)^{\rm QCD} = 2.63(27)$ MeV. We update the low-energy constraints on novel scalar and tensor interactions, $\epsilon_{S}$ and $\epsilon_{T}$, at the TeV scale by combining our new estimates for $g^{u-d}_S$ and $g^{u-d}_T$ with precision low-energy nuclear experiments, and find them comparable to those from the ATLAS and the CMS experiments at the LHC.


I. INTRODUCTION
The axial, scalar and tensor charges of the nucleon are needed to interpret the results of many experiments and probe new physics. In this paper, we extend the calculations presented in Refs. [1][2][3] by analyzing eleven ensembles of 2 + 1 + 1 flavors of highly improved staggered quarks (HISQ) [4] generated by the MILC Collaboration [5]. These now include a second physical mass ensemble at a = 0.06 fm, and an ensemble with a = 0.15 fm and M π ≈ 310 MeV. We have also increased the statistics significantly on six other ensembles using the truncated solver with bias correction method [6,7]. The resulting high-statistics data provide better control over various sources of systematic errors, in particular the two systematics: (i) excited-state contamination (ESC) in the extraction of the ground-state matrix elements of the various quark bilinear operators and (ii) the re-Also, the rate of proton-proton fusion, the first step in the thermonuclear reaction chains that power low-mass hydrogen-burning stars like the Sun, is sensitive to it. The current best determination of the ratio of the axial to the vector charge, g A /g V , comes from measurement of neutron beta decay using polarized ultracold neutrons by the UCNA Collaboration, 1.2772 (20) [8,9], and by PERKEO II, 1.2761 +14 −17 [10]. Note that, in the SM, g V = 1 up to second order corrections in isospin breaking [11,12] as a result of the conservation of the vector current.
Given the accuracy with which g u−d A has been measured in experiments, our goal is to calculate it directly with O(1%) accuracy using lattice QCD. The result presented in this paper, g u−d A = 1.218 (25) (30), is, however, about 1.5σ (5%) smaller than the experimental value. In Sec. VII, we compare with the result g u−d A = 1.271 (13) by the CalLat collaboration. We show that the data on seven HISQ ensembles analyzed by both collaborations agree within 1σ and the final difference is due to the chiral and continuum extrapolation-the fits are weighted differently by the data points that are not common. Based on the analysis of the size of the various systematics in Sec. VI, and on the comparison with CalLat calculation, we conclude that our analysis of errors is realistic. Our goal, therefore, is to continue to quantify and control the various sources of errors.
The standard model does not contain fundamental scalar or tensor interactions. However, loop effects and new interactions at the TeV scale can generate effective interactions at the hadronic scale that can be probed in decays of neutrons, and at the TeV scale itself at the LHC. Such scalar and tensor interactions contribute to the helicity-flip parameters b and b ν in the neutron decay distribution [13]. Thus, by combining the calculation of the scalar and tensor charges with the measurements of b and b ν in low energy experiments, one can put constraints on novel scalar and tensor interactions at the TeV scale as described in Ref. [13]. To optimally bound such scalar and tensor interactions using measurements of b and b ν parameters in planned experiments targeting 10 −3 precision [14][15][16], the level of precision required in g u−d S and g u−d T is at the 10% level as explained in Refs. [13][14][15][16]. Future higher-precision measurements of b and b ν would require correspondingly higher-precision calculations of the matrix elements to place even more stringent bounds on TeV-scale couplings.
In a recent work [1], we showed that lattice-QCD calculations have reached a level of control over all sources of systematic errors needed to yield the tensor charge with the required precision. The errors in the scalar 3-point functions are about a factor of 2 larger. In this paper we show that by using the truncated solver method with bias correction [6,7], (for brevity called AMA henceforth), to obtain high statistics on all ensembles, we are also able to control the uncertainty in g u−d S to the required 10% level. These higher-statistics results also improve upon our previous estimates of the axial and the tensor charges.
The matrix elements of the flavor-diagonal tensor operators are needed to quantify the contributions of the u, d, s, c quark electric dipole moments (EDM) to the neutron electric dipole moment (nEDM) [1,17]. The nEDM is a very sensitive probe of new sources of T -and CP -violation that arise in most extensions of the Standard Model designed to explain nature at the TeV scale. Planned experiments aim to reduce the current bound on the nEDM of 2.9 × 10 −26 e cm [18] to around 10 −28 e cm. These will put stringent constraints on many BSM theories provided the matrix elements of novel CP -violating interactions, of which the quark EDM is one, are calculated with the required precision. In Refs. [1,3], we showed that the disconnected contributions are negligible so we update the connected contributions to the flavor diagonal tensor charges for the light u and d quarks that are taken to be degenerate.
The tensor charges are also extracted as the zeroth moment of the transversity distributions, These are measured in many experiments including Drell-Yan and semiinclusive deep inelastic scattering (SIDIS) and describe the net transverse polarization of quarks in a transversely polarized nucleon. There exists an active program at Jefferson Lab (JLab) to measure them [19]. It is, however, not straightforward to extract the transversity distributions from the data taken over a limited range of Q 2 and Bjorken x, and additional phenomenological modeling is required. The lattice results of g u−d T are the most accurate at present as already discussed in Ref. [3]. Future experiments at JLab and other experimental facilities worldwide will significantly improve the extraction of the transversity distributions, and together with accurate calculations of the tensor charges using lattice QCD elucidate the structure of the nucleon in terms of quarks and gluons.
The methodology for calculating the isovector charges in an isospin symmetric theory, that is, measuring the contribution to the matrix elements of the insertion of the zero-momentum bilinear quark operators in one of the three valence quarks in the nucleon, is well developed [1][2][3][20][21][22]. Calculation of the flavor-diagonal charges is similar except that it gets additional contributions from contractions of the operator as a vacuum quark loop that interacts with the nucleon propagator through the exchange of gluons. In Ref. [1], we showed that these contributions to g u,d,s T are small, O(0.01), and consistent with zero within errors. Thus, within current error estimates, the connected contributions alone provide reliable estimates for the flavor diagonal charges g u,d T and the isoscalar combination g u+d T . A detailed analysis of disconnected contributions to the axial, scalar and tensor charges will be presented in a separate paper. This paper is organized as follows. In Sec. II, we describe the parameters of the gauge ensembles analyzed and the lattice methodology. The fits used to isolate excited-state contamination are described in Sec. III. The renormalization of the operators is discussed in Sec. IV. Our final results for the isovector charges and the connected parts of the flavor-diagonal charges are presented in Sec. V. Our estimation of errors is revisited in Sec. VI, and a comparison with previous works is given in Sec. VII. In Sec. VIII, we provide constraints on novel scalar and tensor interactions at the TeV scale using our new estimates of the charges and precision beta decay experiments and compare them to those from the LHC. Our final conclusions are presented in Sec. IX.

II. LATTICE METHODOLOGY
The parameters of the eleven ensembles used in the analysis are summarized in Table I. They cover a range of lattice spacings (0.06 < ∼ a < ∼ 0.15 fm), pion masses (135 < ∼ M π < ∼ 320 MeV) and lattice sizes (3.3 < ∼ M π L < ∼ 5.5) and were generated using 2+1+1flavors of HISQ fermions [4] by the MILC Collaboration [5]. Most of the details of the methodology, and the calculational and the analyses strategies are the same as described in Refs. [1,3]. Here we will summarize the key points to keep the paper self-contained and highlight the new features and analysis.
We construct the correlation functions needed to calculate the matrix elements using Wilson-clover fermions on these HISQ ensembles. Such mixed-actions, cloveron-HISQ, are a nonunitary formulation and suffer from the problem of exceptional configurations at small, but a priori unknown, quark masses. We monitor all correlation functions for such exceptional configurations in our statistical samples. For example, evidence of exceptional configurations on three a15m310 lattices prevents us from analyzing ensembles with smaller M π at a = 0.15 fm using the clover-on-HISQ approach. The same holds for the physical mass ensemble a12m130.
The parameters used in the construction of the 2-and 3-point functions with clover fermions are given in Table II. The Sheikholeslami-Wohlert coefficient [24] used in the clover action is fixed to its tree-level value with tadpole improvement, c sw = 1/u 3 0 , where u 0 is the fourth root of the plaquette expectation value calculated on the hypercubic (HYP) smeared [25] HISQ lattices.
The masses of light clover quarks were tuned so that the clover-on-HISQ pion masses, M val π , match the HISQon-HISQ Goldstone ones, M sea π . Both estimates are given in Table I. All fits in M 2 π to study the chiral behavior are made using the clover-on-HISQ M val π since the correlation functions, and thus the chiral behavior of the charges, have a greater sensitivity to it. Henceforth, for brevity, we drop the superscript and denote the cloveron-HISQ pion mass as M π . Performing fits using the HISQ-on-HISQ values, M sea π , does not change the estimates significantly.
The highlights of the current work, compared to the results presented in Ref. [3], are • The addition of a second physical pion mass ensemble a06m135 and the coarse a15m320 ensemble.
• The new a12m220L simulations replace the older a12m220L O data. In the a12m220L O calculation, the HP analysis had only been done for τ = 10, while in the new a12m220L data the HP calculation has been done for all values of source-sink separation τ , and the bias correction applied. We have also increased the number of LP measurements on each configurations and both HP and LP source points are chosen randomly within and between configurations. Even though the results from the two calculations are consistent, as shown in Tables XII, XIII and XIV, nevertheless, for the two reasons stated above, we will, henceforth, only use the a12m220L data in the analysis of the charges and other quantities in this and future papers.
• All ensembles are analyzed using the AMA method with much higher statistics as listed in Table I. Our implementation of the AMA method is described in Refs. [1,26].
• The new high statistics data for ensembles a09m310, a09m220 and a09m130W were generated using the smearing parameter σ = 7. This corresponds to a r.m.s. radius of ≈ 7.5 in lattice units or roughly 0.66 fm. As discussed in Sec. III and shown in Figs. 9-17, increasing σ from 5.5 to 7.0 reduces the ESC at a given source-sink separation τ .
• The two-point correlation functions are analyzed keeping up to four states in the spectral decomposition. Previous work was based on keeping two states.
• The three-point functions are analyzed keeping up to three states in the spectral decomposition of the spectral functions. Previous work was based on keeping two states.
We find that the new higher precision data significantly improved the ESC fits and the final combined CCFV fit used to obtain results in the limits a → 0, the light-quark mass M π → 135 MeV and the lattice volume M π L → ∞.

A. Correlation Functions
We use the following interpolating operator χ to create/annihilate the nucleon state: with {a, b, c} the color indices, C = γ 0 γ 2 the charge conjugation matrix, and q 1 and q 2 denoting the two different flavors of light quarks. The nonrelativistic projection (1 ± γ 4 )/2 is inserted to improve the signal, with the plus and minus signs applied to the forward and backward propagation in Euclidean time, respectively [27]. , of the eleven 2+1+1-flavor HISQ ensembles generated by the MILC collaboration and analyzed in this study are quoted from Ref. [5]. All fits are made versus M val π and finite-size effects are analyzed in terms of M val π L. Estimates of M val π , the clover-on-HISQ pion mass, are the same as given in Ref. [1] and the error is governed mainly by the uncertainty in the lattice scale. In the last four columns, we give, for each ensemble, the values of the source-sink separation tsep used in the calculation of the three-point functions, the number of configurations analyzed, and the number of measurements made using the high precision (HP) and the low precision (LP) truncation of the inversion of the clover operator. The second set of calculations, a09m130W , a06m310W and a06m220W , have been done with the larger smearing size σ that is given in Table II The parameters used in the calculation of the clover propagators. The hopping parameter for the light quarks,κ l , in the clover action is given by 2κ l = 1/(m l + 4). m l is tuned to achieve M val π ≈ M sea π . The parameters used to construct Gaussian smeared sources, {σ, NKG}, are given in the fourth column where NKG is the number of applications of the Klein-Gordon operator and the width of the smearing is controlled by the coefficient σ, both in Chroma convention [23]. The resulting root-mean-square radius of the smearing, defined as r 2 √ S † Sdr/ √ S † Sdr, is given in the last column.
At zero momentum, this operator couples only to the spin-1 2 state. The zero momentum 2-point and 3-point nucleon correlation functions are defined as where α and β are spinor indices. The source is placed at time slice 0, τ is the sink time slice, and t is an intermediate time slice at which the local quark bilinear The Dirac matrix Γ is 1, γ 4 , γ i γ 5 and γ i γ j for scalar (S), vector (V), axial (A) and tensor (T) operators, respectively. In this work, subscripts i and j on gamma matrices run over {1, 2, 3}, with i < j.
The nucleon charges g q Γ are obtained from the ground state matrix element N (p, s)|O q Γ |N (p, s) , that, in turn, are extracted using the spectral decomposition of the 2and 3-point correlation functions. They are related as with spinors satisfying To extract the charges, we construct the projected 2and 3-point correlation functions The operator P 2pt = (1 ± γ 4 )/2 is used to project on to the positive parity contribution for the nucleon propagating in the forward (backward) direction. For the connected 3-point contributions, P 3pt = P 2pt (1 + iγ 5 γ 3 ) is used. Note that the C 3pt Γ (t, τ ) defined in Eq. (7) becomes zero if Γ anticommutes with γ 4 , so only Γ = 1, γ 4 , γ i γ 5 and γ i γ j elements of the Clifford algebra survive. The fits used to extract the masses, amplitudes and matrix elements from the 2-and 3-point functions, defined in Eqs. (6) and (7), are discussed in Sec. III.

B. High Statistics Using the AMA Method
We have carried out high-statistics calculation on all the ensembles using the truncated solver method with bias correction (called AMA for brevity) [6,7]. In this method, correlation functions are constructed using quark propagators inverted with high precision (HP) and low precision (LP). The bias corrected correlators on each configuration are then given by where C LP and C HP are the 2-and 3-point correlation functions constructed using LP and HP quark propagators, respectively, and x LP i and x HP i are the source positions for the two kinds of propagator inversion. The LP stopping criteria, defined as r LP ≡ |residue| LP /|source| varied between 10 −3 and 5 × 10 −4 , while that for the HP calculations between 10 −7 and 10 −8 .
As discussed in Ref. [26], to reduce statistical correlations between measurements, N HP maximally separated time slices were selected randomly on each configuration and on each of these time slices, N LP /N HP LP source positions were again selected randomly. The number of sources, N LP and N HP , used are given in Table I. An important conclusion based on all our calculations with O(10 5 ) measurements of nucleon charges and form factors carried out so far (see Refs. [1,3,26,28,29]), is that the difference between the LP and the bias corrected estimates (or the HP) is smaller than the statistical errors.
To further reduce the computational cost, we also used the coherent sequential source method discussed in Ref. [26]. Typically, we constructed four HP or LP sequential sources on four sink time slices, and added them to obtain the coherent source. A single inversion was then performed to construct the coherent sequential propagator. This was then contracted with the four original propagators to construct the four three-point functions. All of these propagators were held in the computer memory to remove the I/O overhead.
Our final errors are obtained using a single elimination jackknife analysis over the configurations, that is, we first construct the average defined in Eq. (8) on each configuration. Because of this "binning" of the data, we do not need to correct the jackknife estimate of the error for correlations between the N LP LP measurements per configuration.

III. EXCITED-STATE CONTAMINATION
To extract the nucleon charges we need to evaluate the matrix elements of the currents between ground-state nucleons. The lattice nucleon interpolating operator given in Eq. (1), however, couples to the nucleon, all its excitations and multiparticle states with the same quantum numbers. Previous lattice calculations have shown that the ESC can be large. In our earlier works [1,3,26,28], we have shown that this can be controlled to within a few percent using the strategy summarized below.
The overlap between the nucleon operator and the excited states in the construction of the two-and threepoint functions is reduced by using tuned smeared sources when calculating the quark propagators on the HYP smeared HISQ lattices. We construct gaugeinvariant Gaussian smeared sources by applying the three-dimensional Laplacian operator, ∇ 2 , N GS number of times, i.e., (1 + σ 2 ∇ 2 /(4N GS )) NGS on a delta function source.
The input smearing parameters {σ, N GS } for each ensemble are given in Table II along with the resulting root-mean-square radius defined as r 2 √ S † Sdr/ √ S † Sdr. We find that, as a function of distance r, the modulus of the sum of the values of the twelve spin-color components at each site, √ S † S, is well described by a Gaussian, and we use this ansatz to fit the data. The results for the root-mean-square radius given in Table II show weak dependence on the lattice spacing or the pion mass for fixed σ, and are roughly equal to the input σ. Throughout this work, the same smearing is used at the source and sink points.
The analysis of the two-point functions, C 2pt , was carried out keeping four states in the spectral decomposition: where the amplitudes and the masses of the four states are denoted by A i and M i , respectively. In fits including more than two states, the estimates of M i and the A i for i ≥ 2 were sensitive to the choice of the starting time slice t min , and the fits were not always stable. The fits were stabilized using the empirical Bayesian procedure described in Ref. [28]. Examples of the quality of the fits are shown in Figs. 22-29 in Ref. [29]. The new results for masses and amplitudes obtained from 2-, 3and 4-state fits are given in Table XII. In Fig. 1, we compare the efficacy of different smearing sizes in controlling excited states in the 2-point data on the three ensembles a09m130, a0bm310 and a06m220. In each case, the onset of the plateau with the larger smearing size occurs at earlier Euclidean time t, however, the statistical errors at larger t are larger. The more critical observation is that, while M 0 overlap, the mass gaps a∆M i are significantly different in two cases. Thus the excited state parameters are not well-determined even with our high statistics, O(10 5 ) measurements, data. More importantly, except for the a06m310 case, the mass gap a∆M 1 obtained is much larger than 2aM π , the value expected if N ππ is the lowest excitation. Based on these observation, we conclude that to resolve the excited state spectrum will require a coupled channel analysis with much higher statistics data.
The results of different fits for the bare charges extracted from the three-point data, given in Table XIII, indicate that these differences in the mass gaps do not significantly effect the extraction of the charges. At current level of precision, the variations in the values of the mass gaps and the corresponding values for the amplitudes compensate each other in fits to the 2-and 3-point data.
The analysis of the zero-momentum three-point functions, C (3pt) Γ (t; τ ) was carried out using 3-states in its spectral decomposition: where the source point is at t i , the operator is inserted at time t, and the nucleon state is annihilated at the sink time slice t f . The source-sink separation is τ ≡ t f − t i . The state |0 represents the ground state and |n with n > 0 the higher states. The A i are the amplitudes for the creation of state i with zero momentum by the nucleon interpolating operator χ. To extract the matrix elements, the amplitudes A i and the masses M i are obtained from the 4-state fits to the two-point functions.
Note that the insertion of the nucleon at the sink time slice t f and the insertion of the current at time t are both at zero momentum. Thus, by momentum conservation, only the zero momentum projections of the states created at the source time slice contribute to the threepoint function.
We calculate the three-point correlation functions for a number of values of the source-sink separation τ that are listed in Table I. To extract the desired matrix element 0|O Γ |0 , we fit the data at all τ and t simultaneously using the 3-state ansatz given in Eq. (10).
In the simultaneous fit to the data versus t and multiple τ , we skip t skip points adjacent to the source and the sink for each τ to remove points with the largest ESC. The value of t skip is tuned for each ensemble and given in Table XIII. The above procedure was followed for each matrix element on each ensemble. The simultaneous fits in t and τ give the τ → ∞ estimate of the desired matrix element 0|O Γ |0 .
To visualize the ESC, we plot the data for the following ratio of correlation functions in Figs. 9-17 and show the various fits corresponding to the results in Table XIII. In the limit t → ∞ and τ − t → ∞, this ratio converges to the charge g Γ . At short times, the ESC is manifest in all cases. For sufficiently large τ , the data should exhibit a flat region about τ /2, and the value should become independent of τ . The current data for g A , g S and g T , with τ up to about 1.4 fm, do not provide convincing evidence of this desired asymptotic behavior.
Comparing the data between the pair of calculations with different smearing size carried out on the three ensembles, a09m130, a06m310 and a06m220, we find a significant reduction in the ESC on increasing the smearing size. Nevertheless, for the axial and tensor charges, the 2and 3 * -state fits give consistent estimates for the ground state matrix elements, and for the two calculations. The agreement between these four estimates has increased our confidence in the control over ESC. The results for g u−d S , obtained using 2-state fits, have larger uncertainty as discussed in Sec. III A, but are again consistent except those from the a06m220 ensemble.
This higher statistics study of the ESC confirms many features discussed in Ref. [3]: , and the convergence to the τ → ∞ value is monotonic and from below.
10% for τ > 1 fm, and the convergence to the τ → ∞ value is also monotonic but from above.  • For a given number of measurements at the same τ and t, the statistical precision of g u−d T is slightly better than that of g u−d A . The data for g u−d S is noisy, especially at the larger values of τ . On many ensembles, it does not exhibit a monotonic increase with τ . Using the current estimates as the criteria, to get g u−d S with the same precision as g u−d A will require ≈ 5 times more measurements.
• The data for each charge and for each sourcesink separation τ becomes symmetric about τ /2 with increasing statistical precision. This is consistent with the cosh(t − τ /2) behavior predicted by Eq. (10) for each transition matrix element.
• The variations in the results with the fit ranges selected for fits to the two-point functions and the number, t skip , of points skipped in the fits to the three-point data with multiple values of τ decrease with the increased statistical precision.
• Estimates from the 2-and the 3 * -state fits overlap for all fourteen measurements of g u−d A and g u−d T .
• The 3 * -state fits for g u−d S are unreliable, with too namy poorly determined parameters. To extract our best estimates, we use 2-state fits.
• The largest excited-state contribution comes from the 0|O Γ |1 transition matrix elements. We, therefore discuss a poor person's recipe to get estimates based on the 2 * fits in Sec. III A from data at a single value of τ .
Our conclusion on ESC is that with O(10 5 ) measurements, 3 * fits, the choice of smearing parameters used and the values of τ simulated, the excited-state contamination in g u−d A and g u−d T has been controlled to within a couple of percent, i.e., the size of the quoted errors. The data for g u−d S are more noisy, and we take results from the 2-state fit as our best estimates. In general, for calculations by other groups, where data with reasonable precision are only available at a single value of τ , we show that the 2 * fit gives a much better estimate than the plateau value.
A. A poor person's recipe and g u−d S Our high statistics calculations allow us to develop the following poor person's recipe for estimating the ground state matrix element when data are available only at a single value of τ . To illustrate this, we picked two values of τ ≈ 1 fm (τ = {6, 7}, {8, 10}, {10, 12}, {16, 20} in lattice units for the a ≈ 0.15, 0.12, 0.09, 0.06 ensembles) for which we have reasonably precise data at all values of t and for all three isovector charges. We then compared the estimates of the charges from the 2 * fit to data at these values of τ with our best estimate from the 3 * fit (2-state for g u−d S ) to the data at multiple τ and t. Fits for all ensembles are shown in Figs. 9-17 and the results collected in Table XIII. In the case of g u−d A and g u−d T we get overlapping results. This suggests that, within our statistical precision, all the excited-state terms that behave as cosh ∆M (t − τ /2) in the spectral decomposition are well-approximated by the single term proportional to 0|O|1 in the 2 * fit. Isolating this ESC is, therefore, essential. Also, the remainder, the sum of all the terms independent of t is small. This explains why the values of the excited state matrix elements 1|O|1 and 0|O|2 , given in Table IV, are poorly determined.
We further observe that in our implementation of the lattice calculations-HYP smoothing of the lattices plus the Gaussian smearing of the quark sources-the prod- Since this condition holds for the physical nucleon spectrum, it is therefore reasonable to expect that the charges extracted from a 2 * fit to data with τ 1 fm are a good approximation to the τ → ∞ value, whereas the value at the midpoint t = τ /2 (called the plateau value) is not. This is supported by the data for g u−d A and g u−d T shown in Table XIII; there is much better consistency between the 3 * results and 2 * fits to data with a single values of τ 1 fm versus the plateau value.
In this work, the reason for considering such a recipe is that estimates of g u−d S have much larger statistical errors, because of which the data at the larger values of τ do not in all cases exhibit the expected monotonic convergence in τ and have large errors. As a result, on increasing n in an n-state fit to data with multiple values of increasing τ does not always give a better or converged value. We, therefore, argue that to obtain the best estimates of g u−d S one can make judicious use of this recipe, i.e., use 2 * fits to the data with the largest value of τ that conforms with the expectation of monotonic convergence from below. In our case, based on such analyses we conclude that the 2-state fits are more reliable than 3 * fits for g u−d S . These fourteen values of g u−d S used in the final analysis are marked with the superscript † in Table XIII. The same strategy is followed for obtaining the connected contribution to the isoscalar charges, g u+d A,S,T , that are given in Table XIV.

B. Transition and excited state matrix elements
The only transition matrix element that has been estimated with some degree of confidence is 0|O Γ |1 as can be inferred from the results given in Table IV. Including information from Figs. 9-17, our qualitative conclusions on it are: • Estimates of 0|O A |1 vary between −0.1 and −0. 3 and account for the negative curvature evident in the figures. All ground-state estimate of g u−d A converge from below.
• Estimates of 0|O S |1 vary between −0.2 and −0. 5 and account for the larger negative curvature observed in the figures. All ground-state estimates of g u−d S also converge from below. Our long term goal is to improve the precision of these calculations to understand and extract an infinite volume continuum limit value for the transition matrix elements. C. A caveat in the analysis of the isoscalar charges g u+d A,S,T keeping only the connected contribution In this paper, we have analyzed only the connected contributions to the isoscalar charges g u+d A,S,T . The disconnected contributions are not included as they are not available for all the ensembles, and are analyzed for different, typically smaller, values of source-sink separation τ because of the difference in the quality of the statistical signal. Since the proper way to extract the isoscalar charges is to first add the connected and disconnected contributions and then perform the fits using the lattice QCD spectral decomposition to remove excited state contamination, analyzing only the connected contribution introduces an approximation. Isoscalar charges can be defined without having a disconnected part in a partially quenched theory with an additional quark with flavor u . However, in this theory the Pauli exclusion principle does not apply between u and u quarks. The upshot of this is that the spectrum of states in the partially quenched theory is larger, for example, an intermediate u ud state would be the analogue of a Λ baryon 1 . Thus, the spectral decomposition for this partially quenched theory and QCD is different. The problem arises because our nstate fits assume the QCD spectrum since we take the amplitudes and masses of states from the QCD 2-point function when fitting the 3-point function using Eq. (10). One could make fits to 3-point functions leaving all the parameters in Eq. (10) free, but then even 2-state fits become poorly constrained with current data.
We assume that, in practice, the effect due to using the QCD rather than the partially quenched QCD spectra to fit the connected contribution versus t and τ to remove ESC is smaller than the quoted errors. First, the difference between the plateau value in our largest τ data and the τ → ∞ value is a few percent effect, so that any additional systematic is well within the quoted uncertainty. Furthermore, for the tensor charges the disconnected contribution is tiny and consistent with zero, so for the tensor charges one can ignore this caveat. For the axial and scalar charges, the disconnected contribution is between 10-20% of the connected, and possible systematic effect due to extrapolating the connected and disconnected contributions separately is neglected.

IV. RENORMALIZATION OF OPERATORS
The renormalization constants Z A , Z V , Z S and Z T of the isovector quark bilinear operators are calculated in the regularization-independent symmetric momentumsubtraction (RI-sMOM) scheme [30,31]. We use the same method as given in Refs. [1,3] and refer the reader   to it for details. Results from six ensembles, a12m310, a12m220, a09m310, a09m220, a06m310 and a06m220, given in Table V are reproduced from Refs. [1,3], along with the additional calculation on the a15m310 ensemble. We briefly summarize the method below for completeness.
The calculation was done as follows: starting with the lattice results obtained in the RI-sMOM scheme at a given Euclidean four-momentum squared Q 2 , we first convert them to the MS scheme at the same scale (horizontal matching) using two-loop perturbative relations expressed in terms of the coupling constant α MS (Q 2 ) [32]. This estimate at µ 2 = Q 2 , is then run in the continuum in the MS scheme to 2 GeV using the 3-loop anomalous dimension relations for the scalar and tensor bilinears [33,34]. These data are labeled by the Q 2 in the original RI-sMOM scheme and suffer from artifacts due to non-perturbative effects and the breaking of the Euclidean O(4) rotational symmetry down to the hypercu-bic group. To get the final estimate, we fit these data versus Q 2 using an ansatz motivated by the form of possible artifacts as discussed in Refs. [1,3].
We find that the final renormalization factors on ensembles with constant a show no significant dependence versus M π . We, therefore, average the results at different M π to get the mass-independent values at each a.
In Table V, we also give the results for the ratios Z A /Z V , Z S /Z V , and Z T /Z V that show much smaller O(4) breaking, presumably because some of the systematics cancel. From the individual data and the two ratios, Z Γ /Z V and g Γ /g u−d V , we calculate the renormalized charges in two ways: because of the conservation of the vector current. These two sets of renormalized charges are given in Table VI.
We are also interested in extracting flavor diagonal charges which can be written as a sum over isovector (u − d) and isoscalar (u + d) combinations. These com- in the MS scheme at 2 GeV at the four values of the lattice spacing used in our analysis. Results for the a = 0.12, a = 0.09 and a = 0.06 fm ensembles are reproduced from Ref. [3].
The errors are obtained by adding in quadrature the errors in the matrix elements and in the renormalization constants given in Table V. The unrenormalized charges are given in Table III. In the last two columns, we also give the results for the bare, g bare,u−d binations renormalize with the corresponding isovector, Z isovector , and isoscalar, Z isoscalar , factors that are, in general, different [35] 2 . Only the isovector renormalization constants are given in Table V.
In perturbation theory, the difference between Z isovector and Z isoscalar appears at two loops, and is therefore expected to be small. Explicit calculations in Refs. [36][37][38] show that Z isosinglet ≈ Z isovector for the axial and tensor charges. Since the two agree to within a percent, we will assume Z isoscalar in this 2 In general, one considers the singlet and non-singlet combinations in a N f -flavor theory. In this paper, we are only analyzing the insertions on u and d quarks that are taken to be degenerate, so it is convenient to use the 2-flavor labels isosinglet (u + d) and isovector (u − d).
work, and renormalize both isovector (u−d) and isoscalar (u + d) combinations of charges using Z isovector . In the case of the tensor charges, this approximation is even less significant since the contribution of the disconnected diagrams to the charges is consistent with zero within errors [1].
In the case of the scalar charge, the difference between Z isosinglet and Z isovector can be large due to the explicit breaking of the chiral symmetry in the Wilson-clover action which induces mixing between flavors. This has not been fully analyzed for our clover-on-HISQ formulation, so only the bare results for g u−d S and g u+d S , and the renormalized results for g u−d S are presented in this work.
The final errors are obtained by adding in quadrature the errors in estimates of the ratios g bare Γ /g u−d,bare V to the errors in the ratios of the renormalization constants, given in Table V. The disconnected contributions to g u+d T have been neglected as they were shown to be tiny in Ref. [1]. The rest is the same as in Table III. The three fits, 11-point, 10-point and the 10 * -point, are defined in the text.

V. CONTINUUM, CHIRAL AND FINITE VOLUME FIT FOR THE CHARGES gA, gS, gT
To obtain estimates of the renormalized charges given in Tables VI and VII in the continuum limit (a → 0), at the physical pion mass (M π 0 = 135 MeV) and in the infinite volume limit (L → ∞), we need an appropriate physics motivated fit ansatz. To parameterize the dependence on M π and the finite volume parameter M π L, we resort to results from finite volume chiral perturbation theory (χFT) [39][40][41][42][43][44][45]. For the lattice discretization effects, the corrections start with the term linear in a since the action and the operators in our clover-on-HISQ formalism are not fully O(a) improved. Keeping just the leading correction term in each, plus possibly the chiral logarithm term discussed below, our approach is to make a simultaneous fit in the three variables to the data from the eleven ensembles. We call these the CCFV fits. For the isovector charges and the flavor diagonal axial and tensor charges, the ansatz is where M ρ in the chiral logarithm is the renormalization scale.
The coefficients, c log 3 , are known in χPT, and with lattice QCD data at multiple values of M π and at fixed a and M π L one can compare them against values obtained from the fits. As shown in Fig. 2, the M π dependence of all three isovector charges is mild and adequately fit by the lowest order term. Since the c log 3 predicted by χPT are large, including it requires also including still higher order terms in M π to fit the mild dependence. In our case, with data at just three values of M π and the observed mild dependence between 320 and 135 MeV, including more than one free parameter is not justified based on the Akaike Information Criterion (AIC) that requires the reduction of χ 2 by two units for each extra parameter. In short, we cannot test the predictions of χPT. For example, in a fit including the chiral log term and a M 3 π term, the two additional terms essentially negate each other over the range of the data, i.e., between 320-135 MeV. If the large χPT value for the coefficient c log 3 of the chiral log is used as an input, then the fit pushes the coefficient of the M 3 π term to also be large to keep the net variation within the interval of the data small. Furthermore, as can be seen from Table VIII, even the coefficients of the leading order terms are poorly determined for all three charges. This is because the variations between points and the number of points are both small. For these reasons, including the chiral logarithm term to analyze the current data does not add predictive capability, nor does it provide a credible estimate of the uncertainty due to the fit ansatz, nor tests the χPT value of the coefficient c log 3 . Consequently, the purpose of our chiral fit reduces to getting the value at M π = 135 MeV. We emphasize that this is obtained reliably with just the leading chiral correction since the fits are anchored by the data from the two physical pion mass ensembles.
The finite-volume correction, in general, consists of a number of terms, each with different powers of M π L in the denominator and depending on several low-energy constants (LEC) [43]. We have symbolically represented these powers of M π L by X(M π L) and Y (M π L). Since the variation of these factors is small compared to the exponential over the range of M π L investigated, we set X(M π L) = Y (M π L) = constant and retain only the appropriate overall factor M n π e −MπL , common to all the terms in the finite-volume expansion, in our fit ansatz. The, a posteriori, justification for this simplification is that no significant finite volume dependence is observed in the data as shown in Fig. 2.
We have carried out four fits with different selections of the fourteen data points and for the two constructions of the renormalized charges. Starting with the 14 calculations, we first construct a weighted average of the pairs of points from the three a09m130, a06m310 and a06m220 ensembles. For errors, we adopt the Schmelling procedure [46] assuming maximum correlation between the two values from each ensemble. This gives us eleven data points to fit.
• The fit with all the data points is called the 11point fit. This is used to obtain the final results.
• Remove the coarsest a15m310 ensemble point from the analysis. This is called the 10-point fit.
• Remove the a12m220S point as it has the largest errors and the smallest volume. This is called the 10 * -point fit.
• To compare results for g u−d A with those from the CalLat collaboration [47] (see Sec. VII), we perform an 8−point fit that neglects the data from the three a ≈ 0.06 fm ensembles.
The results from these four fits and for the two ways of constructing the renormalized isovector charges are given in Table VI. We find that the six estimates for g u−d A and g u−d T from the 11-point, 10-point and 10 * -point fits with the two ways of renormalization overlap within 1σ. As discussed in Sec. VII, for g u−d A , the a15m310 point plays an important role in the comparison with the CalLat results.
For final results we use the 11-point fit to the isovector charges renormalized using g bare Γ /g bare V × Z Γ /Z V as some of the systematics cancel in the double ratio. These fits are shown in Fig. 2.
The lattice artifact that has the most impact on the final values is the dependence of g u−d A and g u−d S on the lattice spacing a. As shown in Fig. 2, in these cases the CCFV fit coincides with the fit versus just a (pink and grey bands overlap in such cases). On the other hand, one can see from the middle panels, showing the variation versus M 2 π , that had we only analyzed the data versus M 2 π (grey band), we would have gotten a higher value for g u−d A and a lower one for g u−d S , and both with smaller errors. Our conclusion is that, even when the observed variation is small, it is essential to perform a simultaneous CCFV fit to remove the correlated contributions from the three lattice artifacts.
The data for g u−d T continues to show very little sensitivity to the three variables and the extrapolated value is stable [3]. A large part of the error in the individual data points, and thus in the extrapolated value, is now due to the poorly behaved two-loop perturbation theory used to match the RI-sMOM to the MS scheme in the calculation of the renormalization constant Z T . Further precision in g u−d T , therefore, requires developing more precise methods for calculating the renormalization constants.
Overall, compared to the results presented in Ref. [3], our confidence in the CCFV fits for all three charges has improved with the new higher precision data. The final results for the isovector charges in the MS scheme at 2 GeV from the 11-point fit to data given in Table VI and renormalized using g bare These results for g u−d S and g u−d T meet the target ten percent uncertainty needed to leverage precision neutron decay measurements of the helicity flip parameters b and b ν at the 10 −3 level to constrain novel scalar and tensor couplings, S and T , arising at the TeV scale [3,13].
Results of the 11-point, 10-point, and 10 * -point fits to the connected contributions to the flavor-diagonal charges g u,d A,T , using the isovector renormalization factor Z isovector A,T , respectively, are given in Table VII. Their behavior versus the lattice spacing and the pion mass is shown in Figs. 3 and 4 using the 11-point fits, again with c log 3 = 0 in the ansatz given in Eq. (12). The data exhibit the following features: • The noticeable variation in the axial charges is in g u A versus a which caries over to g u−d A .
• Estimates for the neutron are given by the u ↔ d interchange.

VI. ASSESSING ADDITIONAL ERROR DUE TO CCFV FIT ANSATZ
In this section we reassess the estimation of errors from various sources and provide an additional systematic uncertainty in the isovector charges due to using a CCFV ansatz with only the leading order correction terms. We first briefly review the systematics that are already addressed in our analysis:  FIG. 3. The 11-point CCFV fit using Eq. (12) to the connected data for the flavor diagonal charges g u A and g d A renormalized in the MS scheme at 2 GeV. Only the data for g u A show a notable dependence on the lattice spacing a. The rest is the same as in Fig. 2. an estimate of the possible uncertainty presented in Ref. [3] have not changed. These are reproduced in Tables V and IX, respectively. With the increase in statistical precision of the bare charges, the uncertainty in the Z Γ is now a significant fraction of the total uncertainty in g u−d A,S,T . • Residual uncertainties due to the three systematics, extrapolations to a → 0 and M π L → ∞ and the variation with M π . Estimates of errors in the simultaneous CCFV fit using the lowest order corrections (see Eq. (12)) are given in rows 3-5 in Table IX. These are, in most cases, judged to be small because the variation with respect to the three variables, displayed in Fig. 2, is small. With increased statistics and the second physical mass ensemble, a06m135, our confidence in the CCFV fits and the error estimates obtained with keeping only the lowest-order corrections in each variable has increased significantly. The exception is the dependence of g u−d S on a as highlighted by the dependence of the extrapolated value on whether the a15m310 point is included (11-point fit) or excluded (10-point fit).
Adding the guesstimates for these five systematic uncertainties, given in rows 1-5, in quadrature, leads to an error estimate that is consistent with the errors quoted in Eq. (13) and reproduced in the sixth row of Table IX. We, therefore, regard the fits and the error estimates given in Eq. (13) as capturing the uncertainty due to the five systematics discussed above.
The χ 2 /d.o.f. of all four fits for the axial and tensor charges given in Table VI are already very small. Therefore, adding higher order terms to the ansatz is not justified as per the Akaike Information Criterion [48]. Nevertheless, to be conservative, we budget for a possible additional systematic uncertainty due to the truncation of the CCFV fit ansatz at the leading order in each of the three variables, by examining the variation in the data   12) to the connected data for the flavor diagonal charges g u T and g d T renormalized in the MS scheme at 2 GeV. Only the data for g u T show a notable dependence on Mπ. The rest is the same as in Fig. 2. in Fig. 2.
For g u−d A , the key difference between our extrapolated value and the experimental results are the data on the a ≈ 0.06 fm lattices. As discussed in Sec. VII, an extrapolation in a with and without these ensembles gives g u−d A = 1.218 (25) and g u−d A = 1.245 (42), respectively. The difference, 0.03, is roughly half the total spread between the fourteen values of g u−d A given in Table VI. We, therefore, quote 0.03 as the additional uncertainty due to the truncation of the fit ansatz.
The dominant variation in g u−d S is again versus a, and, as stated above, the result depends on whether the a15m310 point is included in the fit. We, therefore, take half the difference, 0.06, between the 11-point and 10point fit values as the additional systematic uncertainty. One gets a similar estimate by taking the difference in the fit value at a = 0.06 and a = 0. For g u−d T , the largest variation is versus M 2 π . Since we have data from two ensembles at M π ≈ 135 MeV that anchor the chiral fit, we take half the difference in the fit values at M π = 135 and 220 MeV as the estimate of the additional systematic uncertainty.
These error estimates, rounded up to two decimal places, are given in the last row of Table IX. Including them as a second systematic error, our final results for the isovector charges in the MS scheme at 2 GeV are: Error From Estimates of the error budget for the three isovector charges due to each of the five systematic effects described in the text. The symbols ⇑ and ⇓ indicate the direction in which a given systematic is observed to drive the central value obtained from the 11-point fit. The second last row gives the errors quoted in Eq. 13 from the 11-point fit. The last row gives the additional systematic uncertainty assigned to account for possible uncertainty due to the using the CCFV fit ansatz with just the lowest order correction terms.
Similar estimates of possible extrapolation uncertainty apply also to results for the connected contributions to the flavor diagonal charges presented in Eq. (14). Their final analysis, including disconnected contributions, will be presented in a separate publication.
Our new estimate g u−d S = 1.022(80)(60) in the MS scheme at 2 GeV is in very good agreement with g u−d S = 1.02(8)(7) obtained by Gonzalez-Alonso and Camalich [49] using the conserved vector current (CVC) relation g S /g V = (M N − M P ) QCD /(m d − m u ) QCD with the FLAG lattice-QCD estimates [50] for the two quantities on the right hand side. Using CVC in reverse, our predictions for (M N − M P ) QCD , using lattice QCD estimates for m u and m d , are given in  [50] and recent results [51,52].

VII. COMPARISON WITH PREVIOUS WORK
A summary of lattice results for the three isovector charges for N f = 2−, 2+1-and 2+1+1-flavors is shown in Figs. 5, 6 and 7. They show the steady improvement in results from lattice QCD. In this section we compare our results with two calculations published after the analysis and comparison presented in Ref. [3], and that include data from physical pion mass ensembles. These are the ETMC [36,37,53] and CalLat results [47].
The ETMC results g u−d A = 1.212 (40), g u−d S = 0.93(33) and g u−d T = 1.004 (28) [36,37,53] were obtained from a single physical mass ensemble generated with 2-flavors of maximally twisted mass fermions with a clover term at a = 0.0938(4) fm, M π = 130.5(4) MeV and M π L = 2.98. Assuming that the number of quark flavors and finite volume corrections do not make a significant difference, one could compare them against our results from the a09m130W ensemble with comparable lattice parameters: g u−d A = 1.249 (21), g u−d S = 0.952 (74) and g u−d T = 1.011 (30). We remind the reader that this comparison is at best qualitative since estimates from different lattice actions are only expected to agree in the continuum limit.
Based on the trends observed in our CCFV fits shown in Figs. 2-4, we speculate where one may expect to see a difference due to the lack of a continuum extrapolation in the ETMC results. The quantities that exhibit a significant slope versus a are g u−d A and g u−d S . Again, under the assumptions stated above, we would expect ETMC values g u−d A = 1.212 (40) to be larger and g u−d S = 0.93 (33) to be smaller than our extrapolated values given in Eq. (13). We find that the scalar charge fits the expected pattern, but the axial charge does not.
To understand why the results can be different, we first review the notable differences between the two calculations. CalLat uses (i) Möbius domain wall versus clover for the valence quark action. This means that their discretization errors start at a 2 versus a for PNDME. They also have no uncertainty due to the renormalization factor since Z A /Z V = 1 for the Möbius domain wall on HISQ formalism. (ii) They use gradient flow smearing with t gf /a = 1 versus one HYP smearing to smooth high , for N f = 2-2+1-and 2+1+1-flavors. The lattice results are from: PNDME'18 (this work); PNDME'16 [3]; LHPC'12 [71]; PNDME'11 [13] and RQCD'14 [58]. The estimates based on the conserved vector current and phenomenology are taken from Gonzalez-Alonso'14 [49] and Adler'75 [72]. , for N f = 2-2+1-and 2+1+1-flavors. The lattice and phenomenology results quoted from: PNDME'18 (this work); PNDME'16 [3]; PNDME'15 [1] LHPC'12 [71]; RBC/UKQCD'10 [73] ETMC'17 [36]; RQCD'14 [58] and RBC'08 [62]. The phenomenological estimates are taken from the following sources: Kang'15 [74]; Goldstein'14 [75]; Pitschmann'14 [76]; Anselmino'13 [77]; Bacchetta'13 [78] and Fuyuto'13 [79]. frequency fluctuations on the gauge configurations. This can impact the size of statistical errors. (iii) Different construction of the sequential propagator. CalLat inserts a zero-momentum projected axial current simultaneously at all time slices on the lattice to construct the sequential propagator. Their data are, therefore, for the sum of contributions from insertions on all time slices on the lattice, i.e., including contact terms and insertion on time slices outside the interval between the source and the sink. CalLat fits this summed three-point function ver-sus only the source-sink separation τ using the 2-state fit ansatz. (iv) The ranges of τ for which the data have the maximum weight in the respective n-state fits are very different in the two calculations. The CalLat results are obtained from data at much smaller values of τ , which accounts for the smaller error estimates in the data for g u−d A . (v) CalLat analyze the coarser a ≈ 0.15, 0.12 and 0.09 fm ensembles. At a ≈ 0.15 fm, we can only analyze the a15m310 ensemble due to the presence of exceptional configurations in the clover-on-HISQ formulation at lighter pion masses. On the other hand, computing resouces have so far limited CalLat from analyzing the three fine a ≈ 0.06 fm and the physical mass a09m130 ensembles so far.
A combination of these factors could easily explain the ≈ 5% difference in the final values. The surprising result, shown in Table XI, is that estimates on the seven ensembles analyzed by both collaborations are consistent and do not show a systematic difference. (Note again that results from two different lattice formulations are not, a priori, expected to agree at finite a.) These data suggest that differences at the 1σ level (see also our analysis in Table IX) are conspiring to produce a 5% difference in the extrapolated value. Thus, one should look for differences in the details of the CCFV fit.
We first examine the extrapolation in a. A CCFV fit keeping our data from only the eight a ≈ 0.15, 0.12 and 0.09 fm ensembles gives a larger value, g u−d A = 1.245 (42), since the sign of the slope versus a changes sign as is apparent from the data shown in the top three panels of Fig. 2. Thus the three a ≈ 0.06 ensembles play an important role in our continuum extrapolation.
Our initial concern was possible underestimation of statistical errors in results from the a ≈ 0.06 lattices. This prompted us to analyze three crucial ensembles, a09m130, a06m310 and a06m220, a second time with different smearing sizes and different random selection of source points. The consistency between the pairs of data points on these ensembles suggests that statistical fluctuations are not a likely explanation for the size of the undershoot in g u−d A . The possibility that these ensembles are not large enough to have adequately explored the phase space of the functional integral, and the results are possibly biased, can only be checked with the generation and analysis of additional lattices.
The chiral fits are also different in detail. In our data, the errors in the points at M π ≈ 310, 220 and 130 MeV are similar, consequently all points contribute with similar weight in the fits. The errors in the CalLat data from the two physical mass ensembles a15m130 and a12m130 are much larger and the fits are predominately weighted by the data at the heavier masses M π ≈ 400, 350 310 and 220 MeV. Also, CalLat finds a significant change in the value between the M π ≈ {400, 350, 310} MeV and M π ≈ 220 MeV points, and this concerted change, well within 1σ errors in individual points, produces a larger dependence on M π . In other words, it is the uniformly smaller values on the M π ≈ {400, 350, 310} MeV en-sembles compared to the data at M π ≈ 220 MeV that makes the CalLat chiral fits different and the final value of g u−d A larger.
To summarize, the difference between our and CalLat results comes from the chiral fit and the continuum extrapolation. The difference in the chiral fit is a consequence of the "jump" in the CalLat data between M π = {400, 350, 310} and the 220 MeV data. The CalLat data at M π ≈ 130 MeV do not contribute much to the fit because of the much larger errors. We do not see a similar jump between our M π ≈ 310 and 220 MeV or between the 220 and the 130 MeV data as is evident from Fig. 2. Also, our four data points at M π ≈ 310 MeV show a larger spread. The difference in the continuum extrapolation is driven by the smaller estimates on all three fine a ≈ 0.06 ensembles that we have analyzed. Unfortunately, neither of these two differences in the fits can be resolved with the current data, especially since the data on 7 ensembles, shown in for the proton on the seven 2+1+1-flavor HISQ ensembles that have been analyzed by us and the CalLat collaboration [47]. The results are consistent within 1σ in most cases.
Even with the high statistics calculation presented here, the statistical and ESC errors in the calculation of the scalar charge are between 5-15% on individual ensembles. As a result, the error after the continuum extrapolation is about 10%. A summary of results for g u−d S , presented in Fig. 6, shows significant reduction in the error over time.
The variation of the tensor charge g u−d T with a or M π or M π L is small. As a result, the lattice estimates have been stable over time as shown in Fig. 7. The first error estimate in our result, g u−d T = 0.989 (32) (10), is now dominated by the error in Z T .

VIII. CONSTRAINING NEW PHYSICS USING PRECISION BETA DECAY MEASUREMENTS
Nonstandard scalar and tensor charged-current interactions are parametrized by the dimensionless couplings S,T [13,82]: These couplings can be constrained by a combination of low energy precision beta-decay measurements (of the pion, neutron, and nuclei) combined with our results for the isovector charges g u−d S and g u−d T , as well at the Large Hadron Collider (LHC) through the reaction pp → eν + X and pp → e + e − + X. The LHC constraint is valid provided the mediator of the new interaction is heavier than a few TeV.
In Fig. 8 (left) we show current and projected bounds on { S , T } defined at 2 GeV in the M S scheme. The beta decays constraints are obtained from the recent review article Ref. [80]. The current analysis includes all existing neutron and nuclear decay measurements, while the future projection assumes measurements of the various decay correlations with fractional uncertainty of 0.1%, the Fierz interference term at the 10 −3 level, and neutron lifetime with uncertainty δτ n = 0.1s. The current LHC bounds are obtained from the analysis of the pp → e + M ET + X. We have used the ATLAS results [81], at √ s = 13 TeV and integrated luminosity of 36 fb −1 . We find that the strongest bound comes by the cumulative distribution with a cut on the transverse mass at 2 TeV. The projected future LHC bounds are obtained by assuming that no events are observed at transverse mass greater than 3 TeV with an integrated luminosity of 300 fb −1 .
The LHC bounds become tighter on the inclusion of Z-like mediated process pp → e + e − + X. As shown in Fig. 8 (right), including both W -like and Z-like mediated processes, the current LHC bounds are comparable to future low energy ones, motivating more precise low energy experiments. In this analysis we have neglected the NLO QCD corrections [83], which would further strengthen the LHC bounds by O(10%). Similar bounds are obtained using the CMS data [84,85].

IX. CONCLUSIONS
We have presented a high-statistics study of the isovector and flavor-diagonal charges of the nucleon using clover-on-HISQ lattice QCD formulation. By using the truncated solver with bias correction error-reduction technique, we have significantly improved the statistical precision of the data. Also, we show stability in the isolation and mitigation of excited-state contamination by keeping up to three excited states in the analysis of data at multiple values of source-sink separation τ . Together, these two improvements allow us to demonstrate that the excited-state contamination in the axial and the tensor channels has been reduced to the 1 − 2% level. The highstatistics analysis of eleven ensembles covering the range 0.15-0.06 fm in the lattice spacing, M π = 135-320 MeV in the pion mass, and M π L = 3.3-5.5 in the lattice size allowed us to analyze the three systematic uncertainties due to lattice discretization, dependence on the quark mass and finite lattice size, by making a simultaneous fit in the three variables a, M 2 π and M π L. Data from the two physical mass ensembles, a09m130 and a06m135, anchor the improved chiral fit. Our final estimates for the isovector charges are given in Eq. (15).
One of the largest sources of uncertainty now is from the calculation of the renormalization constants for the quark bilinear operators. These are calculated nonperturbatively in the RI-sMOM scheme over a range of values of the scale Q 2 . As discussed in Ref. [3], the dominant systematics in the calculation of the Z's comes from the breaking of the rotational symmetry on the lattice and the 2-loop perturbative matching between the RI-sMOM and the MS schemes.
Our estimate g u−d A = 1.218(25)(30) is about 1.5σ (about 5%) below the experimental value g A /g V = 1.2766 (20). Such low values are typical of most lattice QCD calculations. The recent calculation by the CalLat collaboration, also using the 2+1+1-flavor HISQ ensembles, gives g u−d A = 1.271 (13) [47]. A detailed comparison between the two calculations is presented in Sec VII. We show in Table XI that results from seven ensembles, that have been analyzed by both collaborations, agree within 1σ uncertainty. Our analysis indicates that the majority of the difference comes from the chiral and continuum extrapolations, with 1σ differences in individual points getting amplified. Given that CalLat have not analyzed the fine 0.06 fm ensembles and their data on the two physical pion mass ensembles, a15m130 and a12m130 have much larger errors and do not contribute significantly to their chiral fit, we conclude that our error estimate is more realistic. Further work is, therefore, required to resolve the difference between the two results.
Our results for the isovector scalar and tensor charges, g u−d S = 1.022(80)(60) and g u−d T = 0.989(32)(10), have achieved the target accuracy of 10% needed to put bounds on scalar and tensor interactions, S and T , arising at the TeV scale when combined with experimental measurements of b and b ν parameters in neutron decay experiments with 10 −3 sensitivity [13]. In Sec. VIII, we update the constraints on S and T from both low energy experiments combined with our new lattice results on g u−d S and g u−d T and from Atlas and CMS experiments at the LHC. We find that the constraints from low energy experiments combined with matrix elements from lattice QCD are comparable to those from the LHC.
For the tensor charges, we find that the dependence on the lattice size, the lattice spacing and the light-quark mass is small, and the simultaneous fit in these three variables, keeping just the lowest-order corrections, has improved over that presented in Ref. [1].
We have also updated our estimates for the connected parts of the flavor-diagonal charges. For the tensor charges, the contribution of the disconnected diagram is consistent with zero [1,2], so the connected contribution, g u T = 0.790 (27) and g d T = −0.198(10) for the proton, is a good approximation to the full result that will be discussed elsewhere. The extraction of the scalar charge of the proton has larger uncertainty. The statistical errors in the lattice data for g u−d S (a, M π , M π L) are 3-5 times larger than those in g u−d T (a, M π , M π L), and the data show significant dependence on the lattice spacing a and a weaker dependence on the pion mass M π . Our new estimate, g u−d S = 1.022 (80) (60), is in very good agreement with the estimate g u−d S = 1.02(8)(7) obtained using the CVC relation g S /g V = (M N − M P ) QCD /(m d − m u ) QCD in Ref. [49]. In Table X,  In this Appendix, we first present the masses and amplitudes obtained from fits to the 2-point function using the spectral decomposition Eq. (9) in Table XII. These are used as inputs in the fits to the 3-point functions using Eq. (10). We then give in Tables XIII and XIV the results of 2 * -, 2-and 3 * -state fits used to control the ESC in the extraction of the isovector and the connected contribution to the isoscalar axial, scalar and tensor charges for the fourteen calculations. The data and the 2 * -, 2and 3 * -state fits are shown in Figs. 9-17. In each case, we compare the 2 * fit on data from two source-sink separations with τ ≈ 1 fm with the 2-or 3 * -state fit using data from multiple values of τ .  FIG. 9. Comparison between the 2 * and 3 * fits to the axial charge g u−d A data from the a ≈ 0.15 fm (top row) and a ≈ 0.12 fm (bottom 4 rows) ensembles. The results of the fits are sumamrized in Table XIII along with the number of points t skip skipped. The first two columns show 2 * fits to data versus t at a single value of τ , while the third panel shows the 3 * fit using data at multiple values of τ . The labels give the ensemble ID, and the values of τ used in the fits. The τ → ∞ value is given by the grey band in each panel.    Comparison between the 2 * and 3 * fits to the tensor charge g u−d T data from the a ≈ 0.06 fm ensembles. The rest is the same as in Fig. 9.