Nucleon isovector scalar charge from overlap fermions

We calculate the nucleon isovector scalar charge in lattice QCD using overlap fermions on five ensembles of gauge configurations generated by the RBC/UKQCD collaboration using the domain-wall quark action with $2+1$ dynamical flavors. The five ensembles cover five pion masses, $m_\pi \approx$ 139, 171, 302, 337 and 371 MeV, and four lattice spacings, $a \approx $ 0.06, 0.08, 0.11 and 0.14 fm. Three to six valence quark masses are computed on each ensemble to investigate the pion mass dependence. The extrapolation to the physical pion mass, continuum and infinite volume limits is obtained by a global fit of all data to a formula originated from partially quenched chiral perturbation theory. The excited-states contamination is carefully analyzed with 3--5 sink-source separations and multi-state fits. Our final result, in the $\overline{\text{MS}}$ scheme at 2 GeV, is $g_{S}^{u-d}= 0.94 (10)_{stat}(8)_{sys}$, where the first error is the statistical error and the second is the systematic error.


I. INTRODUCTION
The nucleon scalar charge is a fundamental quantity in understanding the internal structure of nucleons and more importantly it is related to the search for new physics beyond the Standard Model (BSM). Together with the tensor charge, it probes novel scalar and tensor interactions at the TeV scale [1]. The nucleon scalar charge is also an important input in the direct search for dark matter [2][3][4][5]. There are numerous ongoing or planned experiments targeted at searching for scalar and tensor interactions [6][7][8][9][10][11][12]. High precision experimental measurements would require precise input of the scalar/tensor charge to put stringent bounds on the existence of new physics.
Unlike the axial charge, the scalar and tensor charges are not well known in experiments. Lattice QCD provides a first-principles non-perturbative formulation for numerical calculation of the fundamental quantities of the QCD theory with controlled uncertainties. In lattice QCD, the nucleon charges are extracted from the matrix elements of the local quark bilinear operators within the nucleon state. For the isovector charges, only the connected insertions are involved and thus are straightforward to compute. As an extension of our previous work on the isovector axial and tensor charges [13][14][15], we compute the isovector scalar charge in this work.
We use a mixed-action approach with overlap fermions in the valence sector and domainwall configurations. Since both domain-wall and overlap fermions are chiral fermions, the nonperturbative renormalization via chiral Ward identities or RI/MOM scheme can be implemented relatively easily and the systematic uncertainty due to the use of an action explicitly breaking chiral symmetry can be avoided at finite lattice spacing. The multi-mass algorithm for overlap fermions allows us to calculate quark propagators for many different quark masses without much additional cost. Five ensembles covering five pion masses including one at the physical value, four lattice spacings in the range 0.06 fm -0.14 fm and five volumes are used in this work, and 3-6 valence pion masses are computed for each ensemble. This enables us to make a reliable extrapolation to the physical pion mass, continuum limit and infinite-volume limit. Excited-states contamination is a main source of systematic uncertainty in lattice calculations of nucleon matrix elements. In order to investigate the excited-states contamination, we explicitly fit up to three states in the correlation functions with 3-5 source-sink separations. This paper is organized as follows. The numerical details about the lattice setup and the computation of correlation functions are presented in Sec. II. The analysis of the correlation functions, in particular the investigation of the exited-state contamination, is presented in Sec. III. In Sec. IV, we describe the renormalization procedure. In Sec. V, we present the renormalized values of the scalar charge and perform the extrapolation to the physical point. A summary is given in Sec. VI.

A. Lattice setup
The results presented in the paper are based on the gauge configurations generated by the RBC/UKQCD collaboration with 2+1 flavor domain-wall fermions [16][17][18]. The relevant parameters of the five ensembles used in this work are collected in Table I. The gauge action is the Iwasaki+DSDR action [19][20][21] for the ensemble 32ID and the Iwasaki action [22] for the rest of the ensembles. See Refs. [16][17][18] for more details about the ensembles.
For the valence quark, we use the overlap fermion action [23]. The overlap Dirac operator is defined as where is the matrix sign function and D w is the Wilson Dirac operator with a negative mass parameter −ρ = 1/2κ − 4. We set κ = 0.2 in our calculation which corresponds to ρ = 1.5. The massive overlap Dirac operator is To accommodate the SU (3) chiral transformation, it is usually convenient to use the chirally regulated fieldψ = (1 − 1 2 D ov )ψ in lieu of ψ in the interpolating field and currents. This is equivalent to leaving the currents unchanged and adopting the effective propagator [24,25] where D c = ρDov 1−Dov/2 is exactly chiral [26], i.e., {γ 5 , D c } = 0.

B. Computation of the correlation functions
The nucleon isovector scalar charge is defined through the nucleon matrix element where u N is the nucleon spinor at zero momentum with the normalization s u N (s)ū N (s) = 1+γ 4 2 , and u and d are up and down quark fields. This charge can be obtained from the ratio of the   three-point function to the two-point function where χ is the proton interpolating operator χ = abc [u aT Cγ 5 d b ]u c and Γ e = 1 2 (1 + γ 4 ) is the positive parity projector for the nucleon propagating in the forward direction. The scalar current is O S (t) = xq (t, x)q(t, x) with q representing the u/d quark fields. The current insertion time t varies between the source and sink time locations. Smeared grid sources with Z3 noise [27][28][29] are used to compute the high-mode part of the quark propagators, while the low-mode part is constructed exactly using the eigenvectors for each point of the source grid G. The correlation functions are calculated by combining the high-and low-mode parts using the low-mode substitution technique [28], which helps to reduce the corresponding statistical uncertainties significantly.
To suppress the excited-states contamination, Gaussian smearing is applied to all quarks at sink and source. The calculation strategy is the same as in our previous works, e.g., [13,30,31], and please refer to these references for further details.
One of the advantages of overlap fermions is that one can compute the quark propagators with multiple quark masses at a small additional cost compared to the cost for the lightest quark mass.
We employ 3-6 quark masses for each of the ensembles [32]. The corresponding valence pion masses m val π are listed in Table I. In order to reduce the excited-states contamination, the source-sink separation t sink has to be large enough and the insertion time t should be far away from both the source and the sink. For each ensemble, we calculate the correlation functions for a number of values of t sink and at least two of them are above 1 fm. The values of t sink are listed in Table I. The data at all t sink will be fitted simultaneously to extract the scalar charge.
In this work we are interested in the isovector scalar charge, which means only the connected insertion, as illustrated in Fig. 1, needs to be considered when calculating the three-point functions.

III. ANALYSIS OF THE CORRELATION FUNCTIONS
Nucleon charges are given by the matrix elements of the quark bilinear operators between the nucleon ground state. The nucleon interpolating operator used in the calculation contains contributions from excited states. We use the so-called two-state fit, in which the contribution of the first excited state is taken into account, to extract the desired matrix element. Keeping the ground state and the first excited state in the spectral decomposition, the ratio of the three-point function to the two-point function can be written as where g S is the scalar charge and ∆M 1 is the mass difference between the first excited state and the ground state. We use this empirical form to describe the contamination from the excited states. In practice, ∆M 1 can be considered as the effective weighted average of the mass differences between several excited states and the ground state. It is usually higher than the mass difference between the first excited state and the ground state. In Eq. 6, the higher powers of e −∆M 1 t sink are dropped since they are negligibly small. We compute the ratio for the u and d flavors separately, denoted by R u (t sink , t) and R d (t sink , t), respectively, and perform a joint fit of R u (t sink , t) and R d (t sink , t) at all t sink simultaneously using the above formula. The parameter ∆M 1 is common for R u (t, t s ) and R d (t, t s ). Thus there are 7 parameters to be determined in the fit: The isovector scalar charge is then given by g u−d The data points with t close to the sink or source suffer large excited-states contamination.
Those data points should not be included in the fits. In our fits, three points at the source and sink ends were dropped for the ensembles 24I, 32I, 32ID, and 48I, while four points were dropped for the ensemble 32Ifine. This guarantees that the distance between the inserted current and the sink/source is around or larger than 0.25 fm, and also that the χ 2 values of all fits are in reasonable range.
In order to check the contributions from the higher excited states, we performed a three-state fit which retains three states in the spectral decomposition. The ratio of the three-point function to the two-point function takes the form where ∆M 2 is the mass difference between the second excited state and the ground state. The terms arsing from the transition matrix elements between the first and second excited states are dropped in the above equation. They are found to be insignificant in the fits and thus are ignored for better stability in the fits. As for the two-state fit case, we perform a joint fit of R u (t sink , t) and R d (t sink , t) at all t sink simultaneously, keeping the parameters ∆M 1 and ∆M 2 common for Comparing the results of the two-and three-state fits, we found some discrepancies in ensembles 24I and 48I. Note that for these two ensembles the smallest value of t sink is ∼ 0.88 fm, which is rather small and may suffer from large excited-states contamination. We drop the data points with the smallest value of t sink and redo the two-state fit for all ensembles. In Fig. 2, we compare the results of the three types of fit: 1) two-state fit with all data points, labeled as "2state-fit-1" in the figure and the following text; 2) two-state fit excluding the data points with the smallest t sink , labeled as "2state-fit-2"; 3) three-state fit with all data points, labeled as "3state-fit". One can see that the results of the three types of fits agree with each other very well for the ensembles 32I, 32ID and 32Ifine, while for the ensembles 24I and 48I, "2state-fit-2" agrees better with "3state-fit" than the "2state-fit-1" does. Therefore, we take the results of "2state-fit-2" for the ensembles 24I and 48I and "2state-fit-1" results for the ensembles 32I, 32ID and 32Ifine as our two-state fit results.
The difference between the two-state fit and the three-state fit results will be taken as an estimation of the systematic uncertainty due to excited-states contamination. The two-and three-state fit results of unrenormalized g u−d S for all ensembles are collected in Table II.   the values of unrenormalized g u S and g d S . The width of the bands indicates one-sigma statistical uncertainty. For each ensemble we present one valence pion mass m val π as a representative case. Fig. 4 is the same as Fig. 3 except that the data points are fitted to the three-state fit formula Eq. (7).

IV. RENORMALIZATION
We use the regularization independent momentum subtraction (RI/MOM) scheme [33,34] [34]. Thanks to the statistical enhancement of using the volume source propagator [35], the statistical uncertainty at a given RI/MOM scale p 2 can be smaller than 0.1% and the major uncertainty of Z S comes from the systematic one of the estimated 4-loop effect in the perturbative matching between the RI/MOM and MS schemes, and also the value of Λ QCD , scale running, lattice spacing, and fit range. The overall uncertainty on most of the ensembles is about 1.5% [36], but it is slightly larger on the 32ID ensemble since the minimum p 2 we used is smaller and thus the matching uncertainty is larger.
We also investigated the RI/SMOM scheme [37,38] which has a better perturbative matching convergence at least up to 2-loop level [38]. Based on the calculation on the 64I ensemble which has the same setup as the physical point ensemble 48I but smaller lattice spacing (0.084 fm) and also that on 48I (0.114 fm), we obtained the Z MS S (2GeV) values using different RI/MOM and RI/SMOM scales, as shown in Fig. 5. The RI/MOM cases use the momenta along the body diagonal direction with cutoff µ p 4 µ /( µ p 2 µ ) 2 < 0.26 (blue and purple crosses for the results on the coarser and finer lattice spacings), and the RI/SMOM cases use the surface-diagonal momenta such as p = (k, k, 0, 0) and p = (k, 0, k, 0) (blue and purple dots). All the data are normalized   behavior of Z S in the p 2 ∈[4-10] GeV 2 region observed in the 48I case (blue dots, the same as we obtained in Ref. [34]) using the RI/SMOM scheme remains in the 64I case (purple dots), and thus is not a discretization effect. At the same time, the non-linear p 2 dependence for p 2 > 10 GeV 2 becomes weaker at smaller lattice spacing, and thus it would be a discretization effect at O(a 4 p 4 ).
Generally speaking, the p 2 dependence of the SMOM case is highly non-linear and it is hard to find a linear window to eliminate the discretization error. On the other hand, the result through the MOM scheme shows perfect linear p 2 dependence with the slope decreasing in a 2 . Thus we use the RI/MOM scheme in this work and leave further comparisons between MOM and SMOM schemes to a separate work.
Our results of the scalar current renormalization constants are listed in Table III, with two uncertainties from the statistics and systematics.

V. RESULTS
The renormalized values of g u−d S extracted from the two(three)-state fit for all ensembles are shown in the left(right) panel of Fig. 6 as a function of valence pion mass. In order to obtain the result at the physical point, we perform a joint fit of all data points to the following form, g u−d S (m 2 π,val , m 2 π,sea , a, L) = C 0 + C 1 m 2 π,val + C 2 m 2 π,sea + C 3 a 2 + C 4 e −m π,val L .
Notice that the gauge action of the ensemble 32ID is different from that of the other ensembles; the coefficient of the a 2 term for 32ID should not be the same with the others. We denote it by In order to investigate the systematics in the extrapolation, we performed the following alternative fits. 1) We performed the extrapolation with different formulas by adding a log term m 2 π log m 2 π and using a different volume dependent term m 2 π √ mπL e −mπL , and the extrapolated results did not change. 2) For each ensemble, there are 1-2 valences pion masses that are very close to the sea pion mass, i.e., m val π = 383MeV for 32Ifine, m val π = 295MeV and 316MeV for 32I, m val π = 321MeV and 348MeV for 24I, m val π = 149MeV for 48I and m val π = 174MeV for 32ID. We fit these 7 data points to the formula g u−d S (m 2 π,val , m 2 π,sea , a) = C 0 + C 1 (m 2 π,val + m 2 π,sea ) + C 3 a 2 . The volume dependent term is ignored since we found this term is not important in the fit. The extrapolated value of g u−d S is 0.93(0.21). The agreement between the results of this fit using the nearly unitary data points and the original fit using all data points supports the validity of our partially quenched scheme.
3) We dropped the data point of 32Ifine in fit 2) since it has the largest pion mass, and redid the fit. The extrapolated result is g u−d S = 0.99(0.21). The systematic error in the extrapolation is then estimated by the differences between the original fit and these alternative fits.
Another source of systematics comes from the uncertainties of the renormalization factors, which are ∼ 1.5% as shown in Table III. This systematic error is estimated by 1.5% of the central value of g u−d S and is added quadratically to the systematic error due to the excited-states contamination and the extrapolation to obtain the total systematic error.
Our final result of the isovector scalar charge is where the first error is the statistical error and the second error is the systematic error. The width of the bands indicates one-sigma statistical error. The black diamond is the extrapolated value of g u−d S at the physical pion mass, continuum limit and infinite volume.  In Fig. 7, we compare our result with a number of other lattice calculations: ETMC'20 [39], Mainz'19 [40], JLQCD'18 [41], PNDME'18 [42], RQCD'14 [43] and LHPC'12 [44]. ETMC'20 [39] presented the results from a N f = 2+1+1 twisted mass ensemble with physical pion mass and lattice spacing a ∼ 0.08 fm. Mainz'19 [40] [45] and the PNDME [46,47] collaborations are not included in the comparison.
As shown in Refs. [13,48], there is no mixing from the glue and quark disconnected insertions to the connected insertions. Thus, it is meaningful to define u and d contributions separately in the connected insertions. Such separation can be compared to those from the DIS and Drell-Yan experiments for x and g A for example, when such separation is accommodated in the global analysis of the parton distribution functions [48]. Here we present the extrapolated values of the scalar charges for u and d from the connected insertions: g u S (CI) = 4.04 (19) stat (54) sys , g d S (CI) = 3.11 (14) stat (64) sys , where the systematic errors are estimated by the same method as for g u−d S . By the same token, the isoscalar matrix element for the connected insertion in the MS scheme at 2 GeV is g u+d S (CI) = g u S (CI) + g d S (CI) = 7.18(32) stat (80) sys .
The systematic uncertainties of g u S (CI), g d S (CI) and g u+d S (CI) mainly come from the alternative fit 2) as described above. There is a large cancellation in this source of systematic uncertainty when taking the difference of g u S (CI) and g d S (CI), therefore it dose not contribute much to the systematic uncertainty of g u−d S .

VI. SUMMARY
We have presented the result of the nucleon scalar charge from a lattice calculation using overlap fermions on domain-wall configurations. The calculation is performed on five ensembles with various values of the pion mass, lattice spacing and volume, covering the pion mass from the physical value to 371 MeV and lattice spacing from 0.06 fm to 0.14 fm. Using the multimass algorithm for overlap fermions, 3-6 valence quark masses are obtained for each ensemble. Mainz'19 [40], JLQCD'18 [41], PNDME'18 [42], ETMC'20 [39], RQCD'14 [43] and LHPC'12 [44]. The error bars are the statistical and systematic errors added quadratically.