The nucleon axial, tensor and scalar charges and $\sigma$-terms in lattice QCD

We present results on the nucleon axial, scalar and tensor charges computed within lattice Quantum Chromodynamics. We use three ensembles of gauge configurations generated with physical values of the pion mass to compute these quantities to high accuracy avoiding the need of uncontrolled chiral extrapolations. We determine the values for the axial, scalar and tensor charges for each quark flavor. We include all contributions from valence and sea quarks by using improved methods to compute the disconnected quark loops. For the nucleon axial charge we find $g_A= 1.286(23)$ in agreement with the experimental value. In addition, we extract the nucleon $\sigma$-terms and find $\sigma_{\pi N}=41.6(3.8)$ MeV as well as the strangeness content of the nucleon obtaining for the y-parameter $y =0.0740(59)$.


INTRODUCTION
The fundamental role of the nucleon axial charge in the physics of weak interactions and in beyond the standard model (SM) physics makes its non-perturbative determination of central importance. The nucleon axial charge determines the rate of the weak decay of neutrons into protons and provides a quantitative measure of spontaneous chiral symmetry breaking in hadronic physics. It enters in the analysis of neutrinoless double-beta decay and in the unitarity tests of the Cabibbo-Kobayashi-Maskawa matrix. Equally important are the isovector scalar and tensor charges of the nucleon, which provide essential input for probing novel scalar and tensor interactions at the TeV scale [1].
An ab initio calculation of the axial charge, as precisely as known experimentally from neutron beta decay measurements using polarized ultracold neutrons [2,3], will provide a strong validation of Quantum Chromodynamics (QCD). However, the non-perturbative nature of QCD makes a theoretical calculation of the axial charge, an isovector coupling that we will denote by g u−d A , difficult. Lattice QCD provides a rigorous, non-perturbative formulation of QCD on a Euclidean lattice that allows for a numerical simulation with controlled systematic uncertainties. Numerous past lattice QCD studies [4]. underestimated g u−d A and impeded reliable predictions of the other nucleon charges. Only recently an accurate computation of g u−d A was presented [5] that reproduced the experimental value. It was, however, obtained using chiral extrapolations involving ensembles with heavier than physical pions.
In this work, we evaluate g u−d A using simulations carried out directly at the physical pion mass and including the physical strange and charm quarks in the sea. This avoids chiral extrapolation or any modelling of the pion mass dependence eliminating a systematic error that has been problematic in many studies. Reproducing the value of g u−d A within the lattice QCD framework serves as a most valuable benchmark computation for the prediction of the isovector scalar g u−d S and tensor g u−d T charges, also presented here. Another milestone of our work, is the computation of the axial, scalar and tensor charges for each quark flavor separately, namely g f A , g f S and g f T where f denotes the up, down, strange and charm quarks. In particular, the quark flavor axial charge g f A determines the intrinsic spin carried by the quarks in the nucleon and scalar and tensor charges probe novel interactions   [11,12]. We refer to these ensembles as physical point ensembles. N f is the number of quark flavors in the sea, L (T ) is the spatial (temporal) extent of the lattice in lattice units and a is the lattice spacing determined using the nucleon mass. When two errors are given, the first is statistical and the second is systematic.
that may reveal new physics beyond the SM. The nucleon scalar matrix element is also related to the phenomenologically important π − N and strange σ-terms where a number of phenomenological predictions exist, which need to be tested through non-perturbative lattice calculations.

METHODOLOGY
Lattice QCD formulation and ensembles. In this work we use three gauge ensembles with the parameters listed in Table I. These ensembles use the twisted mass fermion discretization scheme [6,7] and include a clover-term [8]. Twisted mass fermions (TMF) provide an attractive formulation for lattice QCD allowing for automatic O(a) improvement [7], an important property for evaluating the quantities considered here. A clover-term is added to the TMF action to allow for smaller O(a 2 ) breaking effects between the neutral and charged pions that lead to the stabilization of simulations with light quark masses close to the physical pion mass. For more details on the TMF formulation see Refs. [9,10] and for the simulation strategy Refs. [11,12].
Two ensembles have been generated with two mass degenerate light quarks (N f = 2) that reproduce the physical pion mass [11] and lattice sizes of 48 3 x 96 and 64 3 x 128 allowing for examining finite volume dependence. The third ensemble has been generated using a lattice of size 64 3 x 128 with two degenerate light quarks and the strange and charm quarks (N f = 2+1+1) in the sea with masses tuned to produce, respectively, the physical mass of the pion, kaon and D s -meson, keeping the ratio of charm to strange quark mass m c /m s 11.8 [13]. For the valence strange and charm quarks we use Osterwalder-Seiler fermions [14] with mass tuned to reproduce the mass of the Ω − and the Λ + c baryons [15], respectively.
Nucleon matrix elements. The axial, tensor and scalar flavor charges g f A,T,S are obtained from the nucleon matrix elements N |j f Γ |N of the current j f Γ (x) = ψ f (x)Γψ f (x) at zero momentum transfer, where f denotes the quark flavor, and Γ = γ µ γ 5 for the axial-vector current, Γ = 1 for the scalar current and Γ = σ µν for the tensor current. The renormalization group invariant Computation of correlators. The nucleon matrix element of the current j f Γ is extracted by computing an appropriately defined three-point correlator G 3pt , as well as the nucleon two-point correlator, G 2pt , at zero momentum. These correlation functions are constructed by creating a state from the vacuum with the quantum number of the nucleon at some initial time (source) that is annihilated at a later time t s (sink), where we take the source time to be zero. Such an operator creates a tower of states and, in order to improve overlap with the ground state, we employ Gaussian smearing [16,17] to the quark fields at the source and the sink. In addition, we apply APE smearing [18] to the gauge links entering the hopping matrix of the Gaussian smearing in order to reduce unphysical ultra-violet fluctuations. We employ a multi-grid solver for the production of quark propagators [19] that has been extended to the case of the twisted mass operator, where it was shown to yield a speed-up of more than one order of magnitude at the physical point as compared to the conjugate gradient method [20].
The resulting propagators are used to construct the two-and three-point correlators. The latter consist of two contributions, one arising when the current couples to a valence quark and one when coupled to a sea quark. The former is referred to as giving rise to a connected and the latter to a disconnected contribution. The connected contributions are evaluated using sequential inversions through the sink. Since in this method t s and the four spin projection matrices needed for the extraction of the charges are fixed, four sets of sequential inversions are performed for each value of t s in the rest frame of the nucleon. For the disconnected contributions, we made a number of improvements as compared to our previous analysis of the cA2.09.48 ensemble where a stochastic estimation of the loops combined with sloppy inversions was employed [21][22][23]. For the analysis of the cB211.072.64 ensemble we instead use a combination of methods optimized for physical point ensembles [24], namely, we employ full dilution in spin-color space and a partial dilution in space-time using Hierarchical Probing (HP) [25] in order to compute those off-diagonal elements exactly. The approach takes advantage of the exponential decay of off-diagonal elements with the distance to reduce stochastic noise. For the light quarks we use HP up to a distance of 2 3 lattice units. Since for the light quarks, this exponential decay is slow we combine with deflation using up to 200 low modes. For the strange and charm quarks deflation is not necessary and we use HP up a distance of 2 2 lattice units. As for the analysis of cA2.09.64 ensemble we also employ the so-called oneend trick [26,27] allowing for an improved signal-to-noise ratio. Renormalization. Lattice QCD matrix elements must be renormalized to extract physical observables. We use the RI MOM scheme [28] to compute non-perturbatively the renormalization functions. We implement it using the momentum source method [29] and remove lattice spacing artifacts by subtracting O(g 2 a ∞ ) terms computed in perturbation theory [30,31]. We distinguish between non-singlet and singlet renormalization functions, where for the latter we compute, in addition to the connected, the disconnected contributions. The non-singlet and singlet renormalization functions Z A , Z P and Z T for the N f = 2 + 1 + 1 cB211.072.64 case are computed using N f = 4 ensembles simulated at the same β value and at five values of the pion mass so the chiral limit can be taken. The values for the renormalization functions are listed in Table II. We estimate the systematic error by varying the fit ranges used for the extrapolation of the RI MOM scale µ 0 → 0 removing any residual dependence on the scale. The non-singlet Z A is scheme and scale independent, while all other renormalization functions are converted to the MS scheme at a scale of 2 GeV. For Z singlet A we use the conversion factor calculated to 2loops in perturbation theory [32]. The conversion factor for Z singlet S and Z singlet T is the same as in the corresponding non-singlet case.  II: Non-singlet (Z ns ) and singlet (Z s ) renormalization constants computed using N f = 2 and N f = 4 ensembles and used for the renormalization of the matrix elements computed for the cA2.09.48 and cA2.09.64 [21,22], and cB211.072.64 ensembles, respectively.
Analysis of correlators. The nucleon matrix element can be extracted by taking a ratio of G 3pt (t s , t ins ) and G 2pt (t s ), where t ins is the time of the current insertion and ∆E is the energy gap between the ground and first excited states. This ratio becomes time independent for large values of t s and t ins yielding a plateau the value of which gives the desired nucleon matrix element, M. In practice, t s cannot be chosen arbitrarily large because the statistical errors grow exponentially with t s . Thus, we need to use the smallest t s that ensures convergence to the nucleon state. In this work we use several values of t s and increase the statistics as we increase t s to keep the statistical error approximately constant, which is essential  to reliably assess excited states [21,33]. In Table III we give the values of t s used for the connected contribution and the associated statistics. A careful analysis is then performed, employing different methods to study ground state convergence [34]: i) We fit the ratio of Eq. (1) to a constant for various values of t s and seek convergence of the fitted value as t s increases. We refer to this approach as the plateau method. ii) We perform two-and threestate fits that take into account the contributions of the first and second excite states, respectively, using multiple values of t s simultaneously. iii) We sum over the insertion time t ins , omitting the source and sink time slices, of the ratio in Eq. (1) to obtain a linear dependence on t s where M is given by the slope. We refer to this approach as the summation method, which is expected to suppress excited states faster than the plateau method. In Fig. 1 we show the analysis from which we extract the nucleon isovector charges. For all isovector charges the two-state fit yields a robust value, which is stable as we increase the lower value of the time separation t low s used in the fits. In addition, the value extracted from the two-state fit is confirmed by the three-state fits. We take as our final result the value extracted from the two-state fit when it agrees with the summation method. For the analysis we use correlated fits. We repeat the same analysis for the cA2.09.48 1 and cA2.09.64 ensembles reaching similar conclusions.
We perform the same analysis for the disconnected contribution with the results for the isoscalar u + d and strange charges given in Fig. 2. The excited state effects are milder and the plateau values converge rapidly yielding agreement with the two-state and summation fits. A similar behaviour is observed for the charm charges. We thus take the plateau value that is in agreement with the two-state fit as our final value. The disconnected contributions are non-zero for the light and strange quarks. The charm contribution is non-zero for the axial and scalar cases. These results clearly demonstrate that disconnected contributions cannot be neglected at the physi- cal point. They are enhanced in comparison to the values obtained at heavier pion masses where for example the disconnected part of g u+d A using a N f = 2 + 1 + 1 ensemble simulated at a pion mass of m π = 370 MeV was -0.07(1) [36] as compared to -0.15(2) for the cB211.072.64 ensemble.

RESULTS
In Table IV and Fig. 3 we present our final values for the isovector charges for the three ensembles. Comparing the values extracted from the two N f = 2 ensembles with Lm π ∼ 3 and Lm π ∼ 4 we see no indication of volume effects. This corroborates our previous results at heavier than physical pion masses where no volume effects were detected for g u−d A [37]. Similarly, a study of cut-off effects using three N f = 2+1+1 ensembles with lattice spacings a = 0.089(5), 0.070(4), 0.056(4) fm, revealed that cut-off effects are negligible for a range of pion masses spanning 260 MeV to 450 MeV [37].
In Fig. 3 we compare with recent results from other lattice collaborations considering only results computed using simulations with approximately physical pion mass The same description of Fig. 1   i.e. excluding chiral extrapolations. Good agreement is seen for g u−d A using simulations over a range of lattice spacings demonstrating that lattice spacing effects remain indeed small. Having found no detectable finite volume and lattice spacings effects, the agreement between results from the N f = 2 + 1 + 1 cB211.072.64 ensemble and the N f = 2 ensembles means that any unquenching effects are smaller than the errors. This conclusion is in agreement with an older study using TMF ensembles at a similar lattice spacing at a pion mass of m π ∼ 370 MeV, where we found g A = 1.141 (18) for N f = 2 + 1 + 1 and g A = 1.140 (27) for N f = 2 [38].
For the case of g u−d S , we find a value that is larger as compared to other lattice QCD determinations. This is because g u−d S increases with t s and we use larger time separations combined with increased statistics allowing us to better control excited states [33], in particular for the cB211.072.64 ensemble, where seven values t s were used. Similarly, our value for g u−d T tends to be smaller since this quantity decreases with increasing values of t s . Given that our analysis for the cB211.072.64 ensemble has the largest statistics and the biggest number of t s , we consider the values extracted using this ensemble as the best determination of these quantities.  The values extracted for the single flavor charges are tabulated in Table V for the cB211.072.64 ensemble, considered as our best determination. We also include the values for the σ-terms. The π − N σ-term is defined as σ πN = 1/2(σ u + σ d ) and it constitutes one of the fundamental low-energy parameters paying a significant role for phenomenological studies of low energy scattering and dark matter searches.

CONCLUSIONS
Results on the nucleon axial, tensor and scalar charges are presented for three ensembles of twisted mass cloverimproved fermions tuned to reproduce the physical value of the pion mass. A notable result of this work is the accurate computation of g u−d A that agrees with the experimental value of 1.2732(23) [35] for all the three ensembles. An additional milestone is the evaluation to an unprecedented accuracy of the flavor charges directly at the physical point taking into account the disconnected contributions. We show that the charm axial charge is nonzero and obtain a value for g s A that is more accurate than recent phenomenological determinations that confirms the smaller values recently suggested by the NNPDF [42] and JAM17 [43] analyses both of which, however, carry a large error. Using the N f = 2 + 1 + 1 ensemble that allows for the most robust determination we find that the intrinsic quark spin contribution in the nucleon is 1 2 ∆Σ = 1 2 f =u,d,s,c g f A = 0.188 (16). We also obtain for the non-singlet combination g u+d−2s A = 0.520 (20).
The accurate evaluation of the isovector scalar and tensor charges to about 10% and 3% error respectively put constraints on the possible allowed scalar and tensor effective S and T and new physics searches. Using the scalar matrix element we extract the nucleon σterms that are important for direct dark matter searches and for phenomenological studies of π − N scattering processes. For the N f = 2 + 1 + 1 ensemble we find σ πN = 41.6(3.8) MeV, in agreement with the pioneering chiral perturbation theory analysis that yielded σ πN ∼ 45 MeV [44]. It confirms a smaller value already suggested from previous lattice QCD studies [45][46][47] as compared to the one recently extracted using data from pionic atoms [48]. The y-parameter, defined as y = 2 N |ss|N N |ūu+dd|N , gives a measure of the strangeness content of the nucleon. We find a value of y = 0.0740(59).
For the scalar couplings f f N = m f m N N |ψ f ψ f |N we find f u N = 0.0169 (18) and f d N = 0.0257 (25) taking the ratio of m u /m d = 0.513(31) from Ref. [13], f s N = 0.0425(60) and f c N = 0.115 (24). These couplings are needed for estimates of Higgs-induced lepton flavor violation [49].