Nucleon scalar and tensor charges using lattice QCD simulations at the physical value of the pion mass

TL;DR

This paper reports on lattice QCD calculations of nucleon scalar and tensor charges at the physical pion mass, including connected and disconnected contributions, with non-perturbative renormalization and analysis of excited state effects.

Contribution

First direct lattice QCD determination of isoscalar, strange, and charm nucleon charges at the physical point with comprehensive systematic analysis.

Findings

Scalar charges: g_S^u=5.20(42)(15)(12), g_S^d=4.27(26)(15)(12), g_S^s=0.33(7)(1)(4), g_S^c=0.062(13)(3)(5)

Tensor charges: g_T^u=0.782(16)(2)(13), g_T^d=-0.219(10)(2)(13), g_T^s=-0.00319(69)(2)(22), g_T^c=-0.00263(269)(2)(37)

Includes non-perturbative renormalization and excited state analysis.

Abstract

We present results on the light, strange and charm nucleon scalar and tensor charges from lattice QCD, using simulations with $N_f=2$ flavors of twisted mass Clover-improved fermions with a physical value of the pion mass. Both connected and disconnected contributions are included, enabling us to extract the isoscalar, strange and charm charges for the first time directly at the physical point. Furthermore, the renormalization is computed non-perturbatively for both isovector and isoscalar quantities. We investigate excited state effects by analyzing several sink-source time separations and by employing a set of methods to probe ground state dominance. Our final results for the scalar charges are $g_S^u = 5.20(42)(15)(12)$, $g_S^d = 4.27(26)(15)(12)$, $g_S^s=0.33(7)(1)(4)$, $g_S^c=0.062(13)(3)(5)$ and for the tensor charges $g_T^u = 0.782(16)(2)(13)$, $g_T^d = -0.219(10)(2)(13)$, $g_T^s=-0.00319(69)(2)(22)$, $g_T^c=-0.00263(269)(2)(37)$ in the $\overline{\rm MS}$ scheme at 2~GeV. The first error is statistical, the second is the systematic error due to the renormalization and the third the systematic arising from possible contamination due to the excited states.

Figures

Tables

Table 1. Table 1: Input parameters ( β , L , a ​ μ 𝛽 𝐿 𝑎 𝜇 \beta,L,a\mu ) of our lattice simulation with the corresponding lattice spacing and pion mass. The systematic error on the lattice spacing given in the second parenthesis is due to the interpolation to 135 MeV pion mass. The value of the lattice spacing is determined from the nucleon mass using 140 times the statistics as compared to what was used in Refs.  Abdel-Rehim:2015owa ; Abdel-Rehim:2015pwa , namely using 𝒪 ​ ( 215 , 000 ) 𝒪 215 000 {\cal O}(215,000) statistics.

$\beta =2.10$, $a=0.0938\left(3\right)\left(2\right)$ fm, $r_0/a=5.32\left(5\right)$
${48}^3\times 96$, $L=4.5$ fm	$a\mu$	0.0009
	$m_{\pi }$ (GeV)	0.1305(4)
	$m_{\pi }L$	2.98

Table 2. Table 2: The statistics of our calculation. N conf subscript 𝑁 conf N_{\rm conf} is the number of gauge configurations analyzed and N src subscript 𝑁 src N_{\rm src} is the number of source positions per configuration. With N tot subscript 𝑁 tot N_{\rm tot} we denote the total number of statistics, i.e. N tot = N conf × N src subscript 𝑁 tot subscript 𝑁 conf subscript 𝑁 src N_{\rm tot}=N_{\rm conf}\times N_{\rm src} . In the case of the disconnected contributions, N src subscript 𝑁 src N_{\rm src} refers to the number of two-point functions. Also given here is the number of high-precision stochastic vectors, N HP subscript 𝑁 HP N_{\rm HP} , produced for the loops as well as the number of low-precision vectors produced, N LP subscript 𝑁 LP N_{\rm LP} , in the cases that the TSM is employed.

Connected: three-point				Disconnected
$t_f/a$	$N_{\mathrm{conf}}$	$N_{\mathrm{src}}$	$N_{\mathrm{tot}}$	Flavor	$N_{\mathrm{src}}$	$N_{\mathrm{conf}}$	$N_{\mathrm{tot}}$	$N_{\mathrm{HP}}$	$N_{\mathrm{LP}}$
10,12,14	579	16	9264	light	100	2137	213700	2250	-
16	542	88	47696	strange	100	2153	215300	63	1024
18	793	88	69784	charm	100	2153	215300	5	1250

Table 3. Table 3: Renormalization functions Z P subscript 𝑍 𝑃 Z_{P} and Z T subscript 𝑍 𝑇 Z_{T} for the gauge ensemble analyzed in this work, given in the twisted basis. Z P subscript 𝑍 𝑃 Z_{P} renormalizes the scalar operator in the physical basis, whereas Z T subscript 𝑍 𝑇 Z_{T} is the same in both bases. The first error is statistical and the second error is the systematic due to continuum extrapolation.

	$Z_P^{\overline{\mathrm{MS}}}$	$Z_T^{\overline{\mathrm{MS}}}$
Singlet	0.4997(38)(177)	0.8515(3)(51)
Non-singlet	0.5012(75)(258)	0.8551(2)(15)

Table 4. Table 4: Results for the nucleon scalar charges with their jackknife errors. In columns two to six we give the results using the plateau method at t s = 0.94 subscript 𝑡 𝑠 0.94 t_{s}=0.94 , 1.13 1.13 1.13 , 1.31 1.31 1.31 , 1.50 1.50 1.50 and 1.69 1.69 1.69  fm for the light and strange charges, whereas for g S c superscript subscript 𝑔 𝑆 𝑐 g_{S}^{c} the separations t s = 0.56 subscript 𝑡 𝑠 0.56 t_{s}=0.56 , 0.75 0.75 0.75 , 0.94 0.94 0.94 and 1.13 1.13 1.13  fm are given. In the last two columns we list the results from the two-state fit and the summation method that are shown in Figs.  4 ,  5 and  6 with the bands (brown and gray respectively). The final value we select for each observable is shown in Table 7 .

Observable	Plateau at each $t_s$					Two-state	Summation
	0.94	1.13	1.31	1.50	1.69	fit	method
$g_S^{u-d}$	0.732(81)	1.180(144)	1.656(271)	0.930(252)	0.927(380)	1.134(217)	0.584(691)
$g_S^{u+d}$ (conn.)	6.129(118)	7.158(215)	7.495(416)	8.014(349)	8.221(520)	8.353(404)	8.937(1.045)
$g_S^{u+d}$ (disc.)	0.033(58)	0.399(71)	0.697(103)	1.013(159)	1.249(257)	1.303(245)	2.009(468)
$g_S^s$	0.220(17)	0.271(22)	0.328(29)	0.356(44)	0.329(68)	0.365(45)	0.456(77)
$g_S^c$	Plateau at each $t_s$					Two-state	Summation
	0.56	0.75	0.94	1.13		fit	method
	0.058(14)	0.059(13)	0.062(13)	0.063(14)		0.067(25)	0.062(56)

Table 5. Table 5: Results for the nucleon g T u − d superscript subscript 𝑔 𝑇 𝑢 𝑑 g_{T}^{u-d} and connected g T u + d superscript subscript 𝑔 𝑇 𝑢 𝑑 g_{T}^{u+d} with their jackknife errors. The results using the plateau method are shown at t s = 0.94 subscript 𝑡 𝑠 0.94 t_{s}=0.94 , 1.13 1.13 1.13 and 1.31 1.31 1.31  fm. In the last two columns we list the results from the summation method and the two-state fit that are shown in Fig. 7 with the brown and gray bands respectively. The final value we select for each observable is shown in Table 7 .

Observable	Plateau value			Two-state	Summation
$t_s$ (fm)	0.94	1.13	1.31	fit	method
$g_T^{u-d}$	1.036(6)	1.011(11)	1.004(21)	0.985(27)	0.925(40)
$g_T^{u+d}$ (conn.)	0.611(5)	0.593(10)	0.582(16)	0.565(21)	0.540(37)

Table 6. Table 6: Comparison of results for the disconnected contributions to the nucleon tensor charge from the plateau method, the summation method and the two-state fit. The results extracted using the plateau method are presented for source-sink time separations t s = 0.75 subscript 𝑡 𝑠 0.75 t_{s}=0.75  fm to t s = 1.31 subscript 𝑡 𝑠 1.31 t_{s}=1.31  fm. The final value selected for each observable is shown in  Table 7 .

Observable	Plateau value					Two-state	Summation
$t_s$ (fm)	0.75	0.84	0.94	1.13	1.31	fit	method
$g_T^{u+d}$ (disc.)	-0.0118(7)	-0.0152(10)	-0.0194(16)	-0.0213(31)	-0.0206(61)	-0.0257( 30)	-0.0270(35)
$g_T^s$	-0.00176(29)	-0.00200(38)	-0.00319(69)	-0.00340(127)	-0.00426(231)	-0.00341(108)	-0.00416(147)
$g_T^c$	-0.00085(88)	0.00010(137)	-0.00263(269)	-0.00731(439)	-0.00750(666)	-0.00226(398)	-0.00062(502)

Table 7. Table 7: Final results of the nucleon’s scalar and tensor charges, selected from the plateau fits. The error in the first parenthesis is the statistical, in the second parenthesis is a systematic error due to the error in the determination of the renormalization functions Z P MS ¯ superscript subscript 𝑍 𝑃 ¯ MS Z_{P}^{\overline{\rm MS}} and Z T MS ¯ superscript subscript 𝑍 𝑇 ¯ MS Z_{T}^{\overline{\rm MS}} , whereas in the third parenthesis is the systematic error taking into account excited state contamination, and is taken as the difference in the mean values from the plateau and two-state fit methods.

$g_S^{u-d}$	$g_S^{u+d}$ (conn.)	$g_S^{u+d}$ (disc.)	$g_S^s$	$g_S^c$
0.930(252)(48)(204)	8.221(520)(291)(132)	1.249(257)(44)(54)	0.329(68)(12)(36)	0.062(13)(3)(5)
$g_T^{u-d}$	$g_T^{u+d}$ (conn.)	$g_T^{u+d}$ (disc.)	$g_T^s$	$g_T^c$
1.004(21)(2)(0)	0.582(16)(3)(1)	-0.0213(31)(1)(44)	-0.00319(69)(2)(22)	-0.00263(269)(2)(37)

Equations

$$S_F\left[\chi ,\overline{\chi },U\right]=a^4\sum _x\overline{\chi }(x)\left(D_W[U]+m_{\mathrm{cr}}+i{\mu }_l{\gamma }_5{\tau }^3-\frac{1}{4}{c}_{\mathrm{SW}}{\sigma }^{\mu \nu }{\mathcal{F}}^{\mu \nu }[U]\right)\chi (x),$$

$${\mathcal{O}}_{S^a}=\bar{q}\frac{{\tau }^a}{2}q,{\mathcal{O}}_{T^a}^{\mu \nu }=\bar{q}{\sigma }^{\mu \nu }\frac{{\tau }^a}{2}q,$$

$$m_s^R=108.6(2.2)\mathrm{MeV},m_c^R=1392.6(23.5)\mathrm{MeV},$$

$$⟨N(p,s^{\prime })|{\mathcal{O}}_{S^a}|N(p,s)⟩$$

$$⟨N(p,s^{\prime })|{\mathcal{O}}_{T^a|}N(p,s)⟩$$

$$G_{\Gamma }^{3pt}(P,\vec{p},\vec{q}=\vec{0},t_s,t_{\mathrm{ins}})$$

$$G^{2pt}(\vec{p},t_s)$$

$$P^4=\frac{1}{4}\left(1\pm {\gamma }_0\right).$$

$$P=P_k\equiv i{P}^4{\gamma }_5{\gamma }_k,k=1,2,3$$

$$J_{\alpha }(x)={ϵ}^{abc}{u}_{\alpha }^a\left[u_{\beta }^b(x)(C{\gamma }_5{)}_{\beta \gamma }{d}_{\gamma }^c(x)\right].$$

$$q_{\mathrm{smear}}^a(\vec{x},t)=\sum _y{F}^{ab}(\vec{x},\vec{y};U(t))q^b(\vec{y},t),F(\vec{x},\vec{y};U(t))={\left(1+{\alpha }_{G}H\right)}^{n_G}(\vec{x},\vec{y};U(t))$$

$$H(\vec{x},\vec{y};U(t))=\sum _{j=1}^3\left(U_j(\vec{x},t){\delta }_{\vec{x}+a\hat{j},\vec{y}}+U_j^{†}(\vec{x}-a\hat{j},t){\delta }_{\vec{x}-a\hat{j},\vec{y}}\right).$$

$${\mathcal{L}}^{(f)}(\Gamma ,t)=\sum _{x_{\mathrm{ins}}}\mathrm{Tr}\left[G_f(x_{\mathrm{ins}};x_{\mathrm{ins}})\Gamma \right],$$

$$G\equiv M_E^{-1}=\frac{1}{N_r}\sum _{r=1}^{N_r}|s_r⟩⟨{\xi }_r|\approx M^{-1}.$$

$${\mathcal{L}}^{(u+d)}(1,t)=\sum _{x_{\mathrm{ins}}}\mathrm{Tr}\left[G_u(x_{\mathrm{ins}};x_{\mathrm{ins}})+G_d(x_{\mathrm{ins}};x_{\mathrm{ins}})\right]⟶\sum _{x_{\mathrm{ins}}}\mathrm{Tr}\left[i{\gamma }_5\left(G_{{\varphi }_u}(x_{\mathrm{ins}};x_{\mathrm{ins}})-G_{{\varphi }_d}(x_{\mathrm{ins}};x_{\mathrm{ins}})\right)\right],$$

$${\mathcal{L}}^{(u+d)}(1,t)=2{\mu }_l\sum _{x_{\mathrm{ins}},y}\mathrm{Tr}\left[G_{{\varphi }_d}(x_{\mathrm{ins}};y)G_{{\varphi }_d}^{†}(y;x_{\mathrm{ins}})\right].$$

$$\frac{1}{N_r}\sum _r|{\xi }_r⟩⟨{\xi }_r|=1+\mathcal{O}\left(\frac{1}{\sqrt{N_r}}\right)$$

$${\mathcal{L}}^{(u+d)}(1,t)=2{\mu }_l\frac{1}{N_r}\sum _{r=1}^{N_r}⟨s_r|s_r⟩+\mathcal{O}\left(\frac{1}{\sqrt{N_r}}\right),$$

$${\mathcal{L}}^{(u+d)}({\sigma }_{\mu \nu },t)=2{\mu }_l\frac{1}{N_r}\sum _{r=1}^{N_r}⟨s_r|{\sigma }_{\mu \nu }|s_r⟩+\mathcal{O}\left(\frac{1}{\sqrt{N_r}}\right)$$

$$M_{E_{\mathrm{TSM}}}^{-1}=\frac{1}{N_{\mathrm{LP}}}\sum _{j=N_{\mathrm{HP}}+1}^{N_{\mathrm{HP}}+N_{\mathrm{LP}}}|s_j{⟩}_{\mathrm{LP}}⟨{\xi }_j|+\frac{1}{N_{\mathrm{HP}}}\sum _{r=1}^{N_{\mathrm{HP}}}\left[|s_r{⟩}_{\mathrm{HP}}-|s_r{⟩}_{\mathrm{LP}}\right]⟨{\xi }_r|,$$

$$R_{\Gamma }(P,t_s,t_{\mathrm{ins}})=\frac{G_{\Gamma }^{3pt}(P,\vec{0},\vec{0},t_s,t_{\mathrm{ins}})}{G^{2pt}(\vec{0},t_s)}.$$

$$R_{\Gamma }(P,t_s,t_{\mathrm{ins}})\propto \frac{{\sum }_{n,n^{\prime }}⟨J|n^{\prime }⟩⟨n|\bar{J}⟩⟨n^{\prime }|{\mathcal{O}}_{\Gamma }|n⟩e^{-E_{n^{\prime }}(t_s-t_{\mathrm{ins}})}{e}^{-E_n(t_{\mathrm{ins}}-t_0)}}{{\sum }_n|⟨J|n⟩{|}^2{e}^{-E_n(t_s-t_0)}},$$

$$R_{\Gamma }(P,t_s,t_{\mathrm{ins}})\underset{\Delta (t_s-t_0)\gg 1}{\overset{\Delta (t_{\mathrm{ins}}-t_0)\gg 1}{\to }}{\Pi }_{\Gamma }(P)$$

$${\Pi }_S(P^4)$$

$${\Pi }_T^{ij}(P_k)$$

$$R_{\Gamma }^{\mathrm{sum}}(P,t_s)=\sum _{t_{\mathrm{ins}}=t_0+a}^{t_s-a}{R}_{\Gamma }(P,t_s,t_{\mathrm{ins}}).$$

$$R_{\Gamma }^{\mathrm{sum}}(P,t_s)=C+(t_s-t_0)\mathcal{M}+\mathcal{O}(e^{-\Delta (t_s-t_0)})+\cdots ,$$

$$G_{\Gamma }^{3pt}$$

$$G^{2pt}=c_0{e}^{-E_0(t_s-t_0)}+c_1{e}^{-E_1(t_s-t_0)},$$

$$\mathcal{M}\equiv ⟨0|{\mathcal{O}}_{\Gamma }|0⟩=\frac{A_{00}}{c_0}.$$

$$G^{ns}(p)=\frac{a^{12}}{V}\sum _{x,y,z}⟨u(x)\bar{u}(z)\Gamma d(z)\bar{d}(y)⟩,G^s(p)=\frac{a^{12}}{V}\sum _{x,y,z}⟨u(x)\bar{u}(z)\Gamma u(z)\bar{u}(y)⟩$$

$${\Lambda }_{\Gamma }(p)=(S(p){)}^{-1}{G}_{\Gamma }(p)(S(p){)}^{-1}$$

Full text

I Introduction

The nucleon scalar and tensor charges are fundamental properties of hadron structure but most importantly they are related to the ongoing search for new physics beyond the Standard Model (BSM). The nucleon isovector scalar and tensor charges probe novel scalar and tensor interactions at the TeV scale. Planned neutron $\beta$-decay experiments with higher accuracy would require input on the scalar and tensor charges. Furthermore, the nucleon matrix element of the light, strange and charm scalar quark operator, from which the scalar charge is extracted, is directly related to the nucleon scalar contents or $\sigma$-terms. These quantities are crucial input in experimental dark matter searches Cushman:2013zza that are seeking to directly detect dark matter by measuring the recoil energy of scattering between nuclei and dark matter candidates. These candidates are weakly-interacting massive particles (WIMPs) and according to a number of BSM theories Ellis:2010kf ; Servant:2002hb ; Bertone:2010ww ; Bertone:2004pz they interact with normal matter via elastic scattering. During the scattering process, a WIMP produces a Higgs boson, which then interacts with a nucleon through scalar density operators. For spin-independent elastic scattering, the theoretical expression of the cross section depends quadratically on the nucleon scalar matrix element. This contribution, in fact, brings the largest uncertainty on the nucleon dark matter cross section Ellis:2008hf .

The nucleon tensor charge plays an important role in BSM physics connected to novel CP-violating interactions. Such interactions will lead to a non-zero neutron electric dipole moment (nEDM) and planned experiments to reduce the current bound by two orders of magnitude will constrain many BSM theories. An accurate measurement of the flavor-diagonal tensor charges will be needed in order to translate the new bounds on the nEDM into CP-violating terms in BSM theories and set bounds on new sources of CP-violation Bhattacharya:2016zcn . Accurate values of both scalar and tensor charges are needed also in evaluating the contribution of the CP-odd electron-nucleon interaction contributing to the atomic EDM Yamanaka:2017mef .

Unlike the axial charge, $g_{A}$, the scalar and tensor charges are not well known experimentally. The $0^{+}\rightarrow 0^{+}$ nuclear decays and the radiative pion decay $\pi\rightarrow e\nu\gamma$, respectively, provide limits on their values. Experiments using ultra-cold neutrons are expected to improve these values Bhattacharya:2011qm . In addition, there is a rich experimental program to study the transverse spin structure of the nucleon at Jefferson Lab. A coincidence experiment in Hall A will employ a newly proposed solenoid spectrometer (SoLID) to perform precision measurements from semi-inclusive electro-production of charged pions from transversely polarized 3He target in Deep-Inelastic-Scattering kinematics using 11 and 8.8 GeV electron beams Gao:2010av . SoLID is expected to increase the experimental accuracy of the tensor coupling by an order of magnitude Dudek:2012vr ; Ye:2016prn . On-going experiments at LHC are also probing scalar and tensor interactions for BSM physics at the TeV scale, and they are expected to increase the limits to contributions arising from such interactions by an order of magnitude. This experimental activity makes a precise lattice QCD calculation of the scalar and tensor couplings well-timed: It provides valuable input in the ongoing searches for BSM physics, and sheds light on our understanding of nucleon structure.

Lattice QCD has progressed noticeably in the last few years, due to new algorithmic improvements and the increase in the available computational power. These ongoing advancements allow for lattice QCD simulations at physical values of the pion mass and at increasingly larger volumes. Such simulations eliminate the need for chiral extrapolations, thus reducing a significant source of systematic uncertainties. However, calculations of baryon observables close to or at the physical point have a worse signal-to-noise ratio and larger effects due to excited-state contaminations, rendering such calculations more challenging. Typically one needs one order of magnitude larger statistics as compared to using simulations with heavier pion masses for the same setup. To eliminate excited states one needs larger Euclidean time propagation with exponentially increasing statistical noise and thus large statistics.

In this work we study the light, strange and charm scalar and tensor nucleon charges using a gauge ensemble with two degenerate light flavors ($N_{f}=2$) of twisted mass clover-improved fermions with pion mass fixed to its physical value Abdel-Rehim:2015pwa . Since we are analyzing a single gauge ensemble, cut-off and finite volume effects cannot be evaluated using directly lattice results.

The isovector scalar and tensor charges, $g_{S}^{u-d}$ and $g_{T}^{u-d}$, are straight forward to calculate since they receive only connected contributions arising from the coupling of the operator to valence quarks, as depicted in Fig 1. Several lattice QCD results have been obtained recently including direct evaluation at the physical point Abdel-Rehim:2015owa ; Bhattacharya:2016zcn . The isoscalar charges $g_{S}^{u+d}$ and $g_{T}^{u+d}$ receive additional contributions coming from the coupling of the operator to vacuum quarks, forming disconnected quark loops. The strange and charm charges $g_{S}^{s,c}$ and $g_{T}^{s,c}$ receive purely disconnected contributions which are notoriously difficult to evaluate being computationally very demanding. It is only recently that disconnected diagrams were included in lattice QCD calculations of the scalar and tensor matrix elements Bali:2011ks ; Freeman:2012ry ; Gong:2013vja ; Abdel-Rehim:2013wlz ; Bhattacharya:2015wna ; Yamanaka:2015lfk ; Abdel-Rehim:2016won eliminating an uncontrolled systematic uncertainty. In this work we employ improved stochastic methods to include all the disconnected diagrams with satisfactory accuracy and by applying non-perturbative renormalization to obtain results on the strange and charm scalar and tensor charges with all contributions taken into account directly at the physical point mass.

The paper is organised as follows: In section II we summarize the characteristics of the gauge configurations used and in section III we describe the extraction of the appropriate matrix elements, which for zero momentum transfer yield the charges. In section IV we discuss in detail the lattice QCD computation of both connected and disconnected contributions and their renormalization, and in section V we present our results. In section VI we compare our results with those obtained recently by other lattice QCD groups and in section VII we conclude.

II Simulation details

We analyze an $N_{f}=2$ gauge ensemble produced by the European Twisted Mass Collaboration (ETMC) Abdel-Rehim:2015pwa ; Abdel-Rehim:2015owa at the physical pion mass. The “Iwasaki” improved gauge action Iwasaki:1983ck ; Abdel-Rehim:2013yaa is employed. The lattice volume is $48^{3}\times 96$ and the lattice spacing determined from the nucleon mass is $a=0.0938(2)$ fm. The rest of the parameters regarding this ensemble are listed in Table 1. We shall refer to this ensemble as the “physical ensemble” from now on. In the fermion sector, the twisted mass fermion (TMF) action at maximal twist is employed Frezzotti:2000nk ; Frezzotti:2003ni , including a Clover term Sheikholeslami:1985ij

[TABLE]

where $D_{W}[U]$ denotes the massless Wilson-Dirac operator, $\tau^{3}$ is the third Pauli matrix acting in flavour space, $m_{\rm cr}$ is the bare untwisted mass tuned to its critical value and $\mu_{l}$ is the bare twisted light quark mass. The last term in Eq. (1) is the clover-term, with $c_{\rm SW}$ the so-called Sheikoleslami-Wohlert improvement coefficient, which is fixed to $c_{\rm SW}=1.57551$ Aoki:2005et , $\mathcal{F}^{\mu\nu}[U]$ is the field strength tensor and $\sigma^{\mu\nu}=(1/2)[\gamma_{\mu},\gamma_{\nu}]$. With $\chi(x)$ we denote the light quark doublet in the twisted basis, $\chi=(u,d)$.

The TMF action is particularly attractive for hadron structure calculations as it provides an automatic $\mathcal{O}(a)$ improvement without requiring further operator improvement. Additional advantages are the infrared regularization of small eigenvalues that makes dynamical simulations faster and the simplified renormalization of operators Frezzotti:2003ni ; Jansen:2005cg ; Farchioni:2004us ; Farchioni:2005bh . However, due to $\mathcal{O}(a^{2})$ lattice artefacts that lead to instabilities in the numerical simulations, particularly at quark masses close to their physical values, the addition of a clover term was required. The latter reduces isospin symmetry breaking effects, while preserving the automatic $\mathcal{O}(a)$ improvement. Another advantage of TMF is that scalar matrix elements are multiplicatively renormalizable Dinter:2012tt , hence a mixing between the bare light and strange scalar matrix elements, seen in other Wilson-type fermion actions, does not occur. TMF also obey a powerful property, which allows an effective increase of the signal-to-noise ratio of the disconnected quark loops, known as the one-end trick Michael:2007vn ; Foster:1998vw ; McNeile:2006bz . In section IV we give more details on the techniques used for the computation of the disconnected contributions. The reader interested in more technical details regarding the twisted mass action and the simulations of the gauge ensemble used in this work is referred to Refs. Abdel-Rehim:2015pwa ; Abdel-Rehim:2014nka ; Boucaud:2007uk ; Boucaud:2008xu ; Baron:2009wt .

III Matrix element decomposition

The quantity of interest is the forward nucleon matrix element $\langle N(p)|\mathcal{O}_{\Gamma}|N(p)\rangle$, where $|N(p)\rangle$ is a nucleon state with momentum $p$ and $\mathcal{O}_{\Gamma}$ is either the local scalar or tensor operator. In the physical basis, these operators read

[TABLE]

respectively, where $q=u,d$ and $\sigma^{\mu\nu}=(1/2)[\gamma_{\mu},\gamma_{\nu}]$. The $\tau^{a}$ matrix acts in flavour space. We consider both isovector and isoscalar quantities, for which we take $\tau^{a}=\tau^{3}$ and $\tau^{a}=\mathbb{1}$, respectively. The individual contributions for $g_{S}^{u,d}$ and $g_{T}^{u,d}$ can then be extracted from the isovector and isoscalar combinations. This is equivalent to calculating directly these contributions by substituting $\tau^{a}$ with the corresponding projectors onto the up- and down-quarks in Eq. (2). Unless otherwise specified, all expressions are given in Euclidean space. For the strange and charm quarks we use Osterwalder-Seiler fermions, that is, they are introduced as heavy doublets similar to the light quark doublet, $\chi(f)=(f^{+},f^{-})$, where $f=s,c$. The action for these doublets is the same as Eq. (1), but with the light twisted mass $\mu_{l}$ replaced by the corresponding mass of the given heavy quark, $\mu_{f}$, and $f^{\pm}$ refers to choosing $\pm\mu_{f}$. We have tuned the bare twisted mass of the strange and charm quarks to reproduce the physical $\Omega^{-}$ and $\Lambda_{c}^{+}$ mass, respectively. The values we obtain are $a\mu_{s}=0.0259(3)$ and $a\mu_{c}=0.3319(15)$ Alexandrou:2017xwd . At our fixed lattice spacing we get for the renormalized quark masses at 2 GeV in the $\overline{\rm MS}$-scheme are

[TABLE]

where only statistical errors are quoted. A more complete analysis, including systematic errors will follow in the future. The nucleon scalar and tensor charges can be extracted from the corresponding matrix elements of the operators of Eq. (2) at zero momentum transfer, which are decomposed as

[TABLE]

From the above matrix elements, the scalar and tensor charges can be obtained from $G_{S}(0)\equiv g_{S}$ and $A_{T10}(0)\equiv g_{T}$. Depending on whether the operators are either the individual up- or down-quark contributions or the isovector or isoscalar combinations, the corresponding charge is obtained. We note here that for non-zero momentum transfer, the form factors $B_{T10}(Q^{2})$ and $\tilde{A}_{T10}(Q^{2})$ appear in Eq. (5), where $Q^{2}$ is the momentum transfer in Euclidean space. We do not consider these form factors in this work.

IV Lattice evaluation

IV.1 Correlation functions

In lattice QCD, the matrix elements of Eqs. (4) and (5) are computed by constructing appropriate three-point correlation functions. Since we are interested in extracting the charges we need the matrix elements with zero momentum transfer. We thus give here the corresponding expression for the three-point function for the case $\vec{q}=\vec{0}$ as well as the nucleon two-point function needed for canceling the Euclidean time evolution and unknown overlaps of the interpolating field with the nucleon state:

[TABLE]

where $t_{0}$, $t_{\rm ins}$ and $t_{s}$ are the source, insertion and sink time coordinates, respectively. The projector matrix $P^{4}$ is given by

[TABLE]

For the scalar charge, the unpolarized projector $P=P^{4}$ is used in the three-point function, whereas for the tensor charge, the polarized projector

[TABLE]

is required. We work in the rest frame, i.e. the source and sink carry zero momentum, therefore we also set $\vec{p}=0$. We use the proton interpolating operator given by

[TABLE]

In order to increase the overlap with the proton ground state we apply Gaussian smearing Alexandrou:1992ti ; Gusken:1989ad at the source and sink. The smeared quark fields read

[TABLE]

and $H$ is the hopping term realized as a matrix in coordinate and color space,

[TABLE]

We also apply APE-smearing to the gauge fields that enter the hopping matrix. For the parameters $\alpha_{\rm G}$ and $n_{\rm G}$ of the Gaussian smearing we use the values $\alpha_{\rm G}=4.0$ and $n_{\rm G}=50$, optimized such as to yield a proton root mean square radius of about 0.5 fm. The APE-smearing parameters are $N_{\rm APE}=50$ and $\alpha_{\rm APE}=0.5$. In our calculations we choose the source-positions $(\vec{x}_{0},t_{0})$ randomly in order to decrease correlations among measurements.

As already mentioned, for isovector quantities the disconnected contributions cancel in the isospin limit up to lattice artefacts, which we expect to be small for the twisted mass clover-improved action used here. In order to evaluate the connected three-point function, shown diagrammatically on the left panel of Fig. 1, we use the sequential inversion approach through the sink Dolgov:2002zm . Within this method, the sum over the sink spatial coordinates, $\vec{x}_{s}$, in Eq. (6) is carried out through an inversion of the Dirac operator with an appropriately constructed source that combines the two forward propagators with the projector and the quantum numbers of the interpolating field at the sink. This so-called sequential propagator is thus required per choice of the sink time coordinate $t_{s}$ and sink projector, whereas all insertion times as well as any insertion operator can be obtained practically without additional computational cost. We perform inversions for five sink-time slices, namely $t_{s}/a=10,12,14,16$ and $18$, which correspond to $t\approx 0.9-1.7$ fm in physical units, for the scalar operator where excited states contributions were found to be significantly large and three sink-times for the tensor operator, $t_{s}/a=10,12$ and $14$. We use four separate projectors, namely $P^{4}$ and $P_{k},\,k=1,2,3$ as given in Eqs. (8) and (9).

For the isoscalar as well as the purely disconnected strange and charm quantities one needs to compute the disconnected quark loop and appropriately combine it with the two-point function in order to construct the disconnected three-point function, depicted on the right panel of Fig. 1. The disconnected quark loop for a general $\gamma$-structure $\Gamma$ is of the form

[TABLE]

where $G_{f}(x;y)$ is the propagator of the quark flavor $f$. The sum over all the spatial coordinates $x_{\rm ins}$ implies that one needs to evaluate the so-called all-to-all propagator. It is apparent that it is prohibitively expensive to calculate all-to-all propagators in an exact manner, as this would require $L^{3}$ inversions of the Dirac matrix per quark flavor, compared to two inversions per quark flavor required for the connected three-point function. A feasible alternative is to use stochastic techniques Bitar:1989dn in order to obtain an unbiased estimate of $G_{f}(x_{\rm ins};x_{\rm ins})$, at the cost of introducing stochastic error in the calculation. Briefly, this is usually done by generating a set of $N_{r}$ stochastic sources $|\xi_{r}\rangle$ randomly filled with $\mathbb{Z}_{4}$-noise. Then one solves $M|s_{r}\rangle=|\xi_{r}\rangle$ for $|s_{r}\rangle$ and calculates

[TABLE]

For $N_{r}\rightarrow\infty$, Eq. (14) provides an unbiased estimate of the all-to-all propagator. The number $N_{r}$ required in order to sufficiently suppress the stochastic noise depends on the observable, but in general $N_{r}\sim\mathcal{O}(10^{3})$, which is much smaller than $L^{3}$, hence this calculation is computationally attainable.

As already mentioned, TMFs have a property that allows to reduce the gauge noise of disconnected quark loops. At first we remark that the isoscalar combination of a flavor doublet of the scalar operator transforms into an isovector of the pseudo-scalar operator in the twisted basis, i.e. $\bar{u}u+\bar{d}d=i(\bar{\phi}_{u}\gamma_{5}\phi_{u}-\bar{\phi}_{d}\gamma_{5}\phi_{d})$, where $u,d$ are the quark fields in the physical basis and $\phi_{u},\phi_{d}$ are the quark fields in the twisted basis. The disconnected quark loop for the scalar operator in the twisted basis is then given by

[TABLE]

which, when utilizing the TMF properties becomes

[TABLE]

From this transformation, known as the one-end trick Michael:2007vn ; Foster:1998vw ; McNeile:2006bz , two main advantages emerge. The first is that the gauge fluctuations are significantly reduced due to the $\mu_{l}$ factor, which is $\mathcal{O}(10^{-3})$. The second one is that the sum of $V$ terms (sum over $y$ in Eq. (16)) that appears in transforming the subtraction of propagators into a multiplication, which increases the signal-to-noise ratio from $1/\sqrt{V}$ to $V/\sqrt{V^{2}}$. The two benefits emerging from the one-end trick yield a large reduction in the errors for the same computational cost.

The final expression for the disconnected quark loop of the isoscalar combination of the scalar operator is obtained by using the property

[TABLE]

of the noise sources. Inserting this into Eq. (16), and noting that $|s_{r}\rangle=G|\xi_{r}\rangle$ we get

[TABLE]

where with the bra-ket notation, a trace over spatial volume, spin and color indices must be realized. Similarly, the isoscalar tensor operator in the twisted basis transforms as $\sigma_{\mu\nu}\mathbb{1}\rightarrow i\gamma_{5}\sigma_{\mu\nu}\tau^{3}$, where $\mathbb{1}$ and $\tau^{3}$ act in flavor space. Following the same procedure for the tensor operator we obtain the expression

[TABLE]

with the same noise reduction benefits.

Regarding the scalar matrix element of the strange and charm quarks, we use the heavy doublets in the twisted basis to construct the pseudoscalar current $\frac{i}{2}(\bar{\phi}_{f^{+}}\gamma_{5}\phi_{f^{+}}-\bar{\phi}_{f^{-}}\gamma_{5}\phi_{f^{-}})$, where as already mentioned $f=s,c$ and $f^{\pm}$ refers to taking $\pm\mu_{f}$. Considering both $f^{+}$ and $f^{-}$ to construct these isovector-like combinations in the twisted basis, allows us to take full advantage of the one-end trick and increase the signal-to-noise ratio of the disconnected quark loops in order to obtain the loops $\mathcal{L}^{(s)}$ and $\mathcal{L}^{(c)}$. Similar procedure is followed for the tensor matrix element, as in the light quark loops. Namely, apart from a factor of $1/2$ the same expressions as in Eqs. (18) and (19) are derived for the heavy quarks, where $|s_{r}\rangle$ are obtained by inverting the twisted-mass Wilson-clover operator with the corresponding heavy quark mass.

In addition, for the strange and charm loops we use the Truncated Solver Method (TSM), which provides a way to increase $N_{r}$ at a reduced computational cost. Within this method, a large number of stochastic sources inverted to low-precision and a small number inverted to high-precision are combined to estimate the all-to-all propagator Bali:2009hu ; Alexandrou:2012zz according to

[TABLE]

where the first term is similar to Eq. (14) and the second term corrects for the bias introduced by using low-precision stochastic propagators. The parameters that need to be tuned are the exact number of low- and high-precision sources as well as the low-precision criterion, such that the all-to-all propagator remains unbiased. The latter can be set either by a relaxed stop condition for the residual of the CG algorithm, for instance $|\hat{r}|<10^{-2}$, or equivalently, by fixing the number of iterations. The goal is to choose the ratio $N_{\rm LP}/N_{\rm HP}$ as large as possible, while still ensuring that the final result is unbiased and that $r_{c}\simeq 1$, where $r_{c}$ is the correlation between the $N_{\rm HP}$ propagators in low and high precision. It is customary to set $N_{\rm LP}$ as the number of sources that would be used if the standard stochastic method was to be employed instead of TSM, and then increase the number of $N_{\rm HP}$ until the bias is corrected, see e.g. Refs. Alexandrou:2012zz ; Abdel-Rehim:2013wlz ; Alexandrou:2013wca . We set the low-precision criterion such that $r_{c}\simeq 0.99$, which is sufficiently large for the purposes of our calculation. We show the dependence of $r_{c}$ as a function of the CG iterations on the left panel of Fig. 2 for various values of the twisted mass parameter $\mu$. To determine the exact number of iterations for each $\mu$-value, we interpolate our data as demonstrated on the right panel of Fig. 2 for the case of the strange quark. From this procedure we find $n_{\rm iter}^{\rm LP}=126$ for the strange quark. Following a similar procedure for the charm quark we find $n_{\rm iter}^{\rm LP}=9$. However, by fixing the number of iterations the exact residual might differ amongst the stochastic sources inverted. In order to avoid that, we equivalently opt to set $|\hat{r}_{\rm LP}|=10^{-3}$ as the low-precision criterion in both cases, which still satisfies the condition $r_{c}\simeq 0.99$ and yields iteration numbers very close to the values obtained from the TSM tuning procedure. The values of $N_{\rm HP}$ and $N_{\rm LP}$ for the strange and charm quarks are listed in Table 2. We remark here that applying the TSM method for the light quarks is not as beneficial since, as one can see from the left plot of Fig. 2, the number of iterations required to achieve a high correlation is much larger than for the heavy quarks. In fact, in an attempt to tune the TSM parameters following the discussion of Ref. Bali:2009hu we found that the resulting optimal values reported minimal benefits.

After calculating the two- and three-point functions, we then form the ratio

[TABLE]

In the large time limit, the unknown overlaps of the nucleon interpolating field with the nucleon ground state cancel and the ratio becomes time independent, thus the desired matrix element can be extracted from a fit to a constant. This can be realized by writing Eq. (21) on the hadron level

[TABLE]

where $|n\rangle$ is the $n^{\rm th}$ eigenstate of the QCD Hamiltonian with the quantum numbers of the nucleon, and $E_{n}$ is the rest frame energy of that state. We note that when $|n\rangle=|N\rangle$ and $|n^{\prime}\rangle=|N\rangle$, where $|N\rangle$ is the nucleon ground state, the desired matrix element $\langle N|\mathcal{O}_{\Gamma}|N\rangle$ appears in Eq. (22). The exponential terms containing energies of excited states become small compared to the matrix element and compared to unity when $t_{\rm ins}-t_{0}\gg 1$ and $t_{s}-t_{0}\gg 1$, in which case the ratio reduces to the desired ground state matrix element.

IV.2 Excited states investigation

In order to assure that the extracted matrix element corresponds to the nucleon ground state, we employ three methods to assess whether contributions due to excited states to the ratio of Eq. (22) are sufficiently suppressed.

The first method, known as the plateau method, is commonly used in extracting matrix elements. One computes the three-point function for several sink-source time separations and examines the time dependence of the ratio given by Eq. (22). If $\Delta(t_{\rm ins}-t_{0})\gg 1$ and $\Delta(t_{s}-t_{0})\gg 1$, where $\Delta=E_{1}-E_{0}$ then contributions from excited states are expected to be small and ground dominance leads to a time-independent region (plateau). However, due to the fact that the approach to the plateau value is not monotonic, identifying the plateau region can become a difficult task. For this reason, several $t_{s}$ need to be computed for each matrix element, which should demonstrate a convergence towards a single constant plateau value as $t_{s}$ is increased. As already mentioned, within the sequential inversion through the sink, the sink time-slice is fixed, thus a new set of inversions is required for each new value of $t_{s}$. An additional issue arising as the sink-source time separation increases is that the error to signal increases exponentially, as shown in Fig. 3. Therefore, the statistics required for constant error increase exponentially making the calculation computationally very demanding.

In the plateau region one fits the ratio

[TABLE]

over $t_{\rm ins}$ to a constant to obtain the desired matrix element. The procedure is repeated for several increasing values of $t_{s}$ until the plateau value does not change, in order to ensure that the contaminations from excited states are suppressed. The scalar and tensor charges of the nucleon ground state, at zero momentum transfer in the large Euclidean time, are then extracted from the corresponding ratios

[TABLE]

The second technique we employ is the summation method Maiani:1987by , which has been applied in a number of recent calculations concerning nucleon charges Capitani:2012gj ; Abdel-Rehim:2015owa ; Abdel-Rehim:2016won . In this approach, a sum of the ratio over the insertion time $t_{\rm ins}$ is performed

[TABLE]

From Eq. (22) and keeping terms up to $E_{1}$ one can see that the sum over the exponentials is a geometric series, thus it can be easily carried out and reads

[TABLE]

where $C$ is a constant independent of $t_{s}$ and $\mathcal{M}\equiv\langle N|\mathcal{O}_{\Gamma}|N\rangle$ is the desired matrix element. The matrix element $\mathcal{M}$ is then the slope of a straight line fit of $R_{\Gamma}^{\rm sum}(P,t_{s})$ w.r.t. $(t_{s}-t_{0})$. In general, since we now fit over two parameters $C$ and $\mathcal{M}$, the summation method has larger errors on the matrix element.

The third approach to extract the desired matrix element is the so-called two- or three-state fit. Within this method the contributions due to the first or second excited state are taken into account. In this analysis we consider terms involving the first excited state i.e. we perform a two-state fit. We take into account several values of the sink-source separation and perform a simultaneous combined fit with respect to $t_{\rm ins}$ and $t_{s}$, taking into account all terms that involve the ground state and the first excited state. From Eq. (22), one can see that considering all terms up to the first excited state gives the following expression for the three-point function

[TABLE]

where $E_{0}=m_{N}$, the mass of the nucleon, $A_{nm}=\langle J|n\rangle\langle m|\bar{J}\rangle\langle n|\mathcal{O}_{\Gamma}|m\rangle\;,A_{nm}\in\mathbb{R}$, and we note that $A_{01}=A_{10}$. Similarly, the expression for the two-point function is given by

[TABLE]

where $c_{n}=|\langle J|n\rangle|^{2}$. In the above expressions, $|0\rangle$ denotes the ground state of the nucleon and $|1\rangle$ the first excited state. We perform a simultaneous fit to the three- and two-point functions given in Eqs. (27) and (28) that includes seven fit parameters, namely $A_{00},A_{01},A_{11},c_{0},c_{1},E_{0}$ and $E_{1}$. The desired matrix element $\mathcal{M}$ is then obtained from

[TABLE]

These fits are very robust and enable us to extract the excited state contribution accurately. We note that for a consistency we also perform a direct fit to the ratio of Eq. (22), which includes five fit parameters instead of seven. We find consistent results of the matrix element between the two fits, albeit the errors from the latter fit are larger. We thus use the results extracted from the seven-parameter fit.

We consider that excited states are sufficiently suppressed when the three methods mentioned above yield consistent values for $\mathcal{M}$, and take the plateau fit at the $t_{s}$ for which this agreement holds as our final value. We give as a systematic error due to excited state effects the difference between the mean values extracted from the plateau and two-state fit.

IV.3 Finite lattice spacing and volume effects

Since we are using a single ensemble we cannot directly evaluate lattice artifacts due to the finite lattice spacing and volume. However, we have computed similar matrix elements using $N_{f}=2$ Alexandrou:2010hf ; Alexandrou:2011db ; Alexandrou:2011nr and $N_{f}=2+1+1$ Alexandrou:2013joa ensembles for higher than physical pion masses, where we compared results using three different lattice spacings and different volumes. Within our statistical errors, results on for example the axial charge were found to be consistent for lattice spacings between $a\sim 0.9$ fm and $a\sim 0.6$ fm. A continuum extrapolation yielded a value consistent with that determined with the ensemble at $a\sim 0.9$ fm Alexandrou:2010hf . We thus expect cut-off effects to be small for our improved action. Similarly, comparing the axial charge extracted for $Lm_{\pi}=3.3$ and $Lm_{\pi}=4.3$ we observed no detectable volume effects. A further indication that volume effects are under control is based on a study of the nucleon $\sigma$-terms Abdel-Rehim:2016won extracted from the same scalar matrix element using the same ensemble employed here. Our value for the $\sigma_{\pi N}$ is in agreement with the result of Ref. Durr:2015dna obtained from an extensive analysis using the Feynman-Hellmann theorem which was corrected for finite volume effects. This leads us to expect that finite volume effects are reasonably small, although an investigation of volume effects at high accuracy is called for at the physical point. We are currently investigating volume effects on these quantities using a spatial lattice size of $L_{s}=64$ at the same pion mass as the one of this work.

V Results

In this section we present our results for the light, strange and charm scalar and tensor charges. Connected contributions are computed for five source-sink separations for the unpolarized projector (Eq. (8)) used for the scalar charges, namely $t_{s}/a=10,12,14,16$ and $18$, which correspond to about $0.9,1.1,1.3,1.5$ and $1.7$ fm. For the polarized projector (Eq. (9)) applied for the tensor charges we use $t_{s}/a=10,12$ and $14$, which as becomes clear in the following, prove sufficient for our analysis. The disconnected quark loops are calculated for all time slices. In Table 2 we summarize the statistics of our calculation for both the connected and disconnected contributions. Also listed is the number of high- and low-precision stochastic sources used in the TSM for the strange and charm charges. The statistics for $t_{s}/a=10,12$ and $14$ analyzed in this study are about six times more than Ref. Abdel-Rehim:2015owa . The two larger time separations, namely $t_{s}/a=16$ and $18$ that we have introduced in this work serve as a further check for excited state effects.

In order to renormalize the scalar and tensor matrix elements it is sufficient to evaluate the renormalization functions of the non-singlet and singlet scalar and tensor quark bilinears, which for maximally twisted fermions are given by $Z_{P}$ and $Z_{T}$, for the gauge ensemble we use in this work. We determine them non-perturbatively in the $\overline{\rm MS}$-scheme at a scale of 2 GeV. Details on the computation of the non-singlet renormalization functions are given in Ref. Alexandrou:2015sea . For flavor singlet operators, disconnected fermion lines lead to significant increase in the computational effort. In order to calculate the renormalization coefficients non-perturbatively, we consider the bare vertex functions Gockeler:1998ye

[TABLE]

where $G^{ns}$ and $G^{s}$ are the non-singlet and singlet cases, respectively, and $V$ is the lattice volume. The amputated vertex function can be derived from the vertex function as

[TABLE]

where $S(p)$ is the propagator in momentum-space. For the singlet vertex function the disconnected contribution is amputated using one inverse propagator because the closed quark loop does not have an open leg. More details will be given in a fore-coming publication renormNf2_paper . The values used to renormalize the lattice scalar and tensor matrix elements are given in Table 3 for both isovector and isoscalar quantities. Perturbatively, the difference between the singlet and non-singlet renormalization functions for both $Z_{P}$ and $Z_{T}$ is zero up to two loops, as presented in Ref. Constantinou:2016ieh . A non-zero difference is present only for the scalar and axial operators. In particular, for the scalar operator which breaks chirality similarly to the pseudoscalar and tensor ones, the difference is extremely small for the Iwasaki gauge action combined with the value of $c_{SW}$ used in our simulations. This behavior partly explains the small difference we find non-perturbatively renormNf2_paper as presented in Table 3.

We illustrate our results in what follows using a common format, namely we show two plots per observable: In the first plot we display for each $t_{s}$ the ratio of Eq. (21) as a function of $t_{\rm ins}-t_{s}/2$, such that the midpoint time of the ratio coincides for all source-sink separations at $t_{\rm ins}-t_{s}/2=0$. We also include in the same plot the horizontal bands resulting from fitting the ratio, the summation method and the two-state fit. In the second plot we summarize the three methods we employ in the calculations, by demonstrating the convergence of our results extracted from fitting the ratio in the plateau region as a function of $t_{s}$, as well as, by including the results of the summation method and the two-state fit as we vary the lowest value of $t_{s}$ considered in the fits, denoted by $t_{s}^{\rm low}$. Throughout, all errors including the error bands in the fits are calculated using jackknife resampling.

V.1 Nucleon scalar charge

In this section we discuss our results on the scalar charge of the nucleon. In Fig. 4 we show the two plots as discussed above for both the isovector scalar charge, $g_{S}^{u-d}$, and the connected contributions to the isoscalar scalar charge, $g_{S}^{u+d}$. We note that the isovector scalar charge is noisy since it results from subtracting two large values. This explains the fact that at $t_{s}=1.31$ fm the plateau is consistently higher than all the rest. The statistics for this sink-source separation is the same as for the two smaller, whereas for $t_{s}=1.50$ fm and $t_{s}=1.69$ fm the statistics is about five and seven times more. Thus, we interpret the higher value from the plateau as a statistical effect. The rest of the time separations yield consistent values in the plateau region. We take the value of the plateau at $t_{s}=1.50$ fm, which is in agreement with the results from the summation method and two-state fit, as our final value for $g_{S}^{u-d}$. The difference between the value extracted from the plateau and the two-state fit is taken as a systematic error due to excited state contamination.

Our previous studies of the isoscalar scalar charge, $g_{S}^{u+d}$, have shown large contamination due to excited states Abdel-Rehim:2015owa ; Abdel-Rehim:2013wlz . As can be seen from the results shown in Fig. 4, the apparent curvature and the increasing trend in the plateau regions of the ratio as $t_{s}$ becomes larger confirm this observation. Both values at $t_{s}=1.50$ fm and $t_{s}=1.69$ fm are consistent. The accuracy obtained for $g_{S}^{u+d}$ allows an accurate determination at $t_{s}=1.69$ fm ensuring ground state dominance and it is the value we adopt. The results regarding the isovector and isoscalar scalar charges in this study corroborate the findings from Refs. Abdel-Rehim:2015owa ; Abdel-Rehim:2016won that large source-sink separations and high statistics are required for a reasonable agreement of all three methods. We take the difference between the plateau value at $t_{s}=1.69$ fm and the one extracted from the two-state fit starting at $t_{s}^{\rm low}=1.13$ fm as the systematic error due to residual excited states for both $g_{S}^{u-d}$ and the connected $g_{S}^{u+d}$.

We show the ratio for the disconnected contributions to $g_{S}^{u+d}$ in Fig. 5. As with the connected contribution to $g_{S}^{u+d}$, we need to go to large time separation $t_{s}$ in order to suppress sufficiently the excited states, with the $t_{s}$-dependence being more pronounced. We note that for the disconnected quantities one can produce the ratio for all source-sink separations because the two-point function and the quark loop are produced for all time slices. We plot the dependence on $t_{s}$ on the right panel of Fig. 5 where we show the plateau fits for all $t_{s}$ between about $1.0-2.2$ fm. We select the plateau fit value at $t_{s}=1.69$ fm as our final result for $g_{S}^{u+d}$, which is in agreement with the values obtained for larger values of $t_{s}$. This yields a contribution, which is about 15% of the connected contribution to $g_{S}^{u+d}$ and is approximately the same as the percentage found for a twisted mass ensemble with a pion mass of $m_{\pi}=373$ MeV Abdel-Rehim:2013wlz .

The results on the strange and charm scalar charges are purely disconnected and the ratios from which $g_{S}^{s}$ and $g_{S}^{c}$ are extracted are shown in Fig. 6. A first observation is that both $g_{S}^{s}$ and $g_{S}^{c}$ are non-zero and smaller than the disconnected contributions coming from the light quark loops, as expected. We find that the strange scalar charge is about $5\%$ of the isoscalar connected scalar charge, whereas the charm contribution is two orders of magnitude smaller. The source-sink separation $t_{s}=1.69$ fm at which the ratio for $g_{S}^{s}$ converges is comparable to the light charges. For $g_{S}^{c}$ the relative errors are larger and the plateau values overlap at smaller sink-source time separations. We choose the plateau value $t_{s}=0.94$ fm as our final estimate for $g_{S}^{c}$ with the systematic error due to excited state effects as the difference between this and the value from the two-state fit shown in Fig. 6.

In Table 4 we collect the values for the light, strange and charm scalar charges from the plateau fits for all $t_{s}$. Also listed are the results from the summation method and the two-state fit that are shown in Figs. 4, 5 and 6 with the brown and gray bands, respectively. The final values for the scalar charge that we select from the plateau fits are listed in Table 7, where apart from the statistical error shown in the first parenthesis, we show in the second parenthesis the systematic error estimated from the error in the determination of the renormalization functions $Z_{P}^{\overline{\rm MS}}$, as well a systematic error due to excited state effects in the third parenthesis, taken as the difference in the mean values between the plateau and two-state fit methods.

V.2 Nucleon tensor charge

The results for the tensor charge are illustrated in this section in a similar manner to the discussion of the scalar charge . In Fig. 7 we show the isovector and connected isoscalar tensor charges. From our previous study Abdel-Rehim:2016won we know that excited states effects are less severe for $g_{T}$. Indeed, a milder dependence on $t_{s}$ is confirmed also for this twisted mass ensemble and the ratio appears to converge at $t_{s}=1.31$ fm for both the isovector and isoscalar quantities. This can be inferred also from the two-state fit and the summation method results where consistent values are extracted already for the lowest $t_{s}^{\rm low}=0.94$ fm. Therefore, for the tensor charge the analysis is carried out only for three sink-source time separations.

The disconnected contributions to $g_{T}^{u+d}$ are shown in Fig. 8. As can be seen, we obtain a clear non-zero negative value. The values extracted from fitting the plateau for different sink-source separation show convergence for $t_{s}=1.31$ fm as well as agreement with the summation and two-state fit methods. We, thus, select the plateau value at this source-sink time separation as our final value for $g_{T}^{u+d}$. For the strange and charm tensor charges we obtain very small negative values, which are about 10% of the disconnected light contributions, as demonstrated in Fig. 9. We select the plateau value at $t_{s}=1.13$ fm for both $g_{T}^{s}$ and $g_{T}^{c}$. Combining all disconnected contributions gives an upper bound of about 0.1% when compared with the connected $g_{T}^{u+d}$.

We tabulate the results of the connected tensor charge in Table 5 and of the disconnected contributions in Table 6. As with the scalar charges, the error in the first parenthesis is statistical, in the second parenthesis the systematic error due to the error in the determination of the renormalization functions $Z_{T}^{\overline{\rm MS}}$ and in the third parenthesis the systematic error due to excited state effects, taken as the difference in the mean values between the plateau and two-state fit methods.

VI Comparison with other calculations

In this section we compare the values of the nucleon scalar and tensor charges that we extract from the analysis presented here with a set of other recent lattice calculations as well as with results from phenomenology, when available.

As already mentioned, the scalar and tensor charges have received particular attention recently due to their implication in BSM physics. The isovector charges are the most studied in lattice QCD due to the absence of disconnected contributions. Whereas not all calculations are carried out using simulations directly at the physical pion mass, a number of these results are extrapolated to the physical point, providing a comparison with our results that are computed directly at the physical point. The PNDME collaboration has recently presented results on the isovector scalar and tensor charges from nine $N_{f}=2+1+1$ ensembles using the highly-improved staggered quarks (HISQ) action produced by the MILC collaboration, at three values of the lattice spacing and a pion mass range of about $m_{\pi}=140-315$ MeV Bhattacharya:2016zcn . Chiral extrapolations were performed and systematic uncertainties were studied. We also compare with the RQCD collaboration that obtained results from clover-improved fermions on eleven $N_{f}=2$ ensembles at three lattice spacings and several volumes, with a lowest pion mass of $m_{\pi}=150$ MeV Bali:2014nma . The LHPC Green:2012ej has analyzed a number of $N_{f}=2+1$ ensembles of clover-improved fermions produced by the BMW collaboration, domain-wall-fermions (DWF) by the RBC and UKQCD collaborations, as well as a mixed action scheme which uses DWF valence quarks on Asqtad staggered sea quarks generated by the MILC collaboration spanning a pion mass in the range of $m_{\pi}=149-356$ MeV. Both RQCD and LHPC performed chiral extrapolations and examined lattice artifacts.

In  Fig. 10 we compare our results on the isovector scalar charge with the values extrapolated at the physical point from the aforementioned collaborations, as well as with their value at their smallest pion mass if this is at or below $150$ MeV. In the same plot two results from phenomenological analyses are included. One is obtained by employing a quark model with spherically symmetric quark wave functions Adler:1975he to obtain an estimate for $g_{S}^{u-d}$. The second used a conserved vector current (CVC) relation $g_{S}/g_{V}=\delta M_{N}^{\rm QCD}/\delta m_{q}^{\rm QCD}$, where $\delta M_{N}^{\rm QCD}$ and $\delta m_{q}^{\rm QCD}$ are the mass differences of the proton and neutron and the up- and down-quarks, respectively, in pure QCD Gonzalez-Alonso:2013ura ; Abdel-Rehim:2016won . As can be seen, there is a very good agreement among all lattice calculations as well as with the phenomenology results.

In Fig. 10 we also compare the value we extract for the isovector tensor charge with the lattice calculations from PNDME, RQCD and LHPC, and additionally from the RBC/UKQCD collaboration using $N_{f}=2+1$ domain-wall-fermions at a pion mass range of $m_{\pi}=330-670$ MeV Aoki:2010xg . We furthermore include a number of phenomenology results from Refs. Kang:2015msa ; Goldstein:2014aja ; Pitschmann:2014jxa ; Anselmino:2013vqa ; Radici:2015mwa ; Fuyuto:2013gla . The lattice QCD results are very accurate and show an excellent agreement among them. Their errors are noticeably smaller as compared to the phenomenological results, illustrating the important input that lattice QCD is currently providing.

In Fig. 11 we compare our results for the connected parts of the isoscalar scalar and tensor charges with selected results from the PNDME and the LHPC collaborations at various pion masses, using the lattice ensembles described previously, as well as with TMF from a gauge ensemble at $m_{\pi}=373$ MeV Abdel-Rehim:2013wlz . Regarding $g_{S}^{u+d}$, PNDME obtained results at two pion masses, $m_{\pi}=220$ MeV and $310$ MeV using clover valence quarks on $N_{f}=2+1+1$ HISQ sea fermions Bhattacharya:2013ehc and performed a linear extrapolation to the physical point to obtain $g_{S}^{u+d}=7.15(65)$, which agrees with our value. The same group calculated $g_{T}^{u+d}$ on the nine gauge ensembles for which they obtained the isovector charges Bhattacharya:2016zcn and after performing a chiral extrapolation they obtain $g_{T}^{u+d}=0.598(33)$, which is in good agreement with the value extracted in this work. In general, there is agreement among lattice QCD for $g_{S}^{u+d}$ and $g_{T}^{u+d}$ over a range of pion masses. The tendency for lower values regarding $g_{S}^{u+d}$ and higher values regarding $g_{T}^{u+d}$ at heavier pion masses can be explained by the fact that older results have typically used smaller sink-source time separations. Since these quantities are affected by excited state contaminations that tend to decrease and increase their values, respectively, contributions from excited states might explain the higher and lower values, respectively, obtained in more recent calculations.

Besides our computation, only the PNDME has evaluated the disconnected contributions Bhattacharya:2015gma ; Bhattacharya:2015wna at pion masses around $m_{\pi}=220$ MeV and $m_{\pi}=310$ MeV. Disconnected contributions, besides the physical ensemble, were also computed for a gauge ensemble of $N_{f}=2+1+1$ twisted mass fermions at $m_{\pi}=373$ MeV Abdel-Rehim:2013wlz . We compare our results for the scalar and tensor charges regarding the disconnected isoscalar as well as the strange contributions in Figs. 12 and 13, respectively. Concerning the scalar charges, both ETMC and PNDME find that the disconnected contribution to $g_{S}^{u+d}$ and $g_{S}^{s}$ are non-zero and positive. For $g_{S}^{s}$ the PNDME values tend to be larger compared to our values and further study of the systematics is called for. Concerning the tensor charges, our result for the disconnected contribution to $g_{T}^{u+d}$ is clearly non-zero and negative. While PNDME also find a negative value their error on the value obtained at the smallest pion mass that they analyzed is much larger making it consistent with zero. We note that the preliminary value of Ref. Gambhir:2016jul agrees with our result. Similarly the value of $g_{T}^{s}$ is clearly nonzero and negative.

VII Conclusions

The nucleon scalar and tensor charges are computed within lattice QCD using simulations generated with two dynamical degenerate light quarks with mass fixed to reproduce approximately the physical pion mass. Both isoscalar and isovector combinations are obtained including the disconnected contributions. We also compute the nucleon strange and charm scalar and tensor charges for the first time. After a careful investigation of excited states contributions we obtain in the $\overline{\rm MS}$ at 2 GeV the following values

[TABLE]

[TABLE]

where the first error is the statistical error, the second is the systematic error due to the determination of the renormalization functions and the third error is the systematic error due to the excited states, estimated by taking the difference between the mean value obtained from the plateau and two-state fit methods. We stress that both isovector and isoscalar charges are renormalized non-perturbatively with the non-singlet and singlet renormalization functions, respectively. We find that the disconnected contributions to the tensor charge are negligible whereas for the scalar they make about 15% of the total value in the case of the isoscalar combination. In addition, excited states are found to be more severe in the case of the scalar as compared to the tensor. We note that since these results were produced using one ensemble of twisted mass clover-improved fermions we cannot provide systematics errors due to finite lattice spacing. Results from other lattice QCD groups close to the physical point are only reported for the isovector and the connected isoscalar combinations. Overall lattice QCD results are in agreement with a couple of exceptions and produce non-zero values for the scalar strange and charm charges, whereas the strange and charm tensor charge are consistent with zero.

Acknowledgments

We would like to thank the members of the ETMC for a most enjoyable collaboration. We acknowledge funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 642069. This work was partly supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project IDs s540 and s625 on the Piz Daint system, by a Gauss allocation on SuperMUC with ID 44060 and in addition used computational resources from the John von Neumann-Institute for Computing on the Jureca and the BlueGene/Q Juqueen systems at the research center in Jülich. We also acknowledge PRACE for awarding us access to the Tier-0 computing resources Curie, Fermi and SuperMUC based in CEA, France, Cineca, Italy and LRZ, Germany, respectively. We thank the staff members at all sites for their kind and sustained support. K.H. and Ch. K. acknowledge support from the Cyprus Research Promotion Foundation under contract T$\Pi$E/$\Pi\Lambda$HPO/0311(BIE)/09.