Sigma terms and strangeness content of the nucleon with $N_f=2+1+1$ twisted mass fermions

We study the nucleon matrix elements of the quark scalar-density operator using maximally twisted mass fermions with dynamical light ($u$,$d$), strange and charm degrees of freedom. We demonstrate that in this setup the nucleon matrix elements of the light and strange quark densities can be obtained with good statistical accuracy, while for the charm quark counterpart only a bound can be provided. The present calculation which is performed at only one value of the lattice spacing and pion mass serves as a benchmark for a future more systematic computation of the scalar quark content of the nucleon.


Introduction
A number of experiments (see for instance refs. [1][2][3][4][5]) have been designed to investigate the nature of Dark Matter (DM) by detecting and/or measuring the recoil energy of nuclei hit by a hypothetical DM particle. One very popular example of such DM candidates is a weakly interacting massive particle (WIMP) like the ones that are predicted in a large class of models [6]. However, definite evidences for a direct detection of WIMPs have not been observed up till now. Nonetheless, the various ongoing experiments provide rather severe constraints for the parameters of many DM models.
A possible scenario for the detection of a WIMP type of DM particles relies on the idea that the WIMP -due to its assumed large mass -produces a Higgs boson which in turn couples to the various quark flavour scalar density operators taken between nucleon states, as depicted in Fig. 1. In fact, at zero momentum transfer, the cross section for spin independent (SI) elastic WIMPnucleon (χN) scattering reads [7] The functions G q f depend on several parameters of the particular model σ SI,χN is computed in, among which the WIMP mass, m χ . A precise expression of the cross section in the constrained minimal supersymmetric Standard Model (CMSSM) can be inferred from ref. [7]. The expression above for σ SI,χN is in particular valid in the SU(2) isospin limit for u and d quarks, an approximation that is always understood in the following. What we want to compute is the magnitude of the dimensionless and renormalization group invariant (RGI) coupling f T f , which depends on the mass m q f of the quark of flavour f and the nucleon mass, m N . For short in the following we will refer to the matrix element N|q f q f |N as the "scalar f -quark content" of the nucleon.
As one can see from Eq. (1), the cross section, σ SI,χN depends quadratically on f T f , and it is thus very sensitive to the size of the scalar content contributions of different flavours. Already a O(10%) variation of f T f can lead to significant changes in σ SI,χN . It is therefore necessary to compute accurately and with controlled error the hadronic matrix elements N|q f q f |N . The implications of the hadronic uncertainty in comparing models with present data of WIMP direct detection has been for instance presented recently in [8].
One way to calculate for various flavours the scalar quark content of the nucleon is provided by chiral perturbation theory (χPT). The results of such a calculation are usually parametrized in terms of σ πN and σ 0 defined by the where m l is the common mass of the u and d light quarks. A low energy theorem relates σ πN to the pion nucleon scattering amplitude extrapolated to the Cheng-Dashen point [9]. The functional form of the extrapolation formula can be established using dispersion relations and in this way σ πN has been found to be σ πN = 79 (7) MeV [10]. Note that a recent result which uses Lorentz covariant baryon chiral perturbation theory gives σ πN = 59 (7) MeV [11]. The non-singlet matrix element σ 0 can be obtained by looking at the pattern of the SU(3) symmetry breaking visible in the spectrum of the baryon octet. A determination using χPT obtained in [12] has given σ 0 = 36(7)MeV.
A direct measure of the magnitude of the strange quark content of the nucleon relative to the light quark content is represented by the ratio The y N -parameter can be related to σ πN and σ 0 by y N = 1 − σ 0 /σ πN and one finds, using the previously quoted numbers, y N = 0.44 (13) [13]. This result leads to a surprisingly large strange quark content of the nucleon, and consequently to a large strange quark contribution to the σ SI,χN cross section in Eq. (1). On the other hand, the quoted uncertainty on y N is rather large, i.e. of the order of 30%. As pointed out in ref. [14] and discussed above, it is quite important to provide a precise value for the strange quark content of the nucleon in order to be able to interpret the ongoing and planned experimental searches of WIMPs.
Lattice QCD can provide a determination of these nucleon matrix elements from first principles. The difficulty involved in these computations has limited for a long time the possibility to isolate a physical signal for these quantities. The aim of this paper is to show that it is indeed possible to accurately compute y N and test whether the strange quark content of the nucleon is as large as indicated by χPT. To this end, we will use Wilson lattice QCD with maximally twisted mass fermions since this framework is particularly suited for such a calculation, as we will discuss below.
On the lattice, two approaches have been followed for the calculation of the light and strange quark content of the nucleons. The first is based on the Feynman-Hellman theorem [15] that relates the nucleon scalar matrix element to the dependence of the nucleon mass on the quark masses via the equations where we denote by m s the strange quark mass and the derivative are intended to be evaluated at the quark mass values that correspond to the physical pion and kaon masses. This approach is often referred to as the "spectrum method". The other approach, the so-called "direct method", consists in evaluating directly the matrix elements appearing in Eqs. (4) and (5).
Both approaches are numerically very challenging. In the spectrum method, a number of simulations at different values of the strange quark mass are needed in order to evaluate the derivative. In the direct method, the computation of the disconnected diagrams is required, which is a highly demanding task since they often show quite a bad signal to noise ratio, thus requiring very high statistics. In addition, as pointed out in ref. [16], when lattice discretizations that break chiral symmetry are used, a mixing between the bare light and strange scalar quark density matrix elements occurs under renormalization. This mixing is not present for chiral invariant (e.g. overlap) fermions [17]. It can also be avoided up to O(a 2 ) effects when -like we do in this workmaximally twisted mass fermions are used in a mixed action setup as explained below.
There exist already a number of lattice QCD computations of the strange quark content of the nucleon, see the works of refs. [16,[18][19][20][21][22][23][24][25][26][27][28] and the review paper [13]. Although these calculations indicate that the strange quark content of the nucleon is smaller than suggested from χPT, the quoted results are affected by large statistical errors, making it difficult to reach definite conclusions.
In this paper we present a method which allows to compute the strange quark content of the nucleon with small statistical errors. By means of a benchmark calculation at only one value of the lattice spacing, volume and quark mass we demonstrate that it is indeed possible to achieve for the parameters f Ts and y N (see Eqs. (1) and (3)) a signal to noise ratio significantly (more than 4 standard deviations) different from zero. It remains of course open the question of the size of systematic effects which we plan to address in the future by including in the computation more lattice spacings, volumes and quark masses.
As we said, a main ingredient that enables us to obtain accurate values for the strange quark content of the nucleon is the use of Wilson lattice QCD with (maximally) twisted mass fermions for our simulations. Besides its property of automatic O(a)-improvement, it offers the advantage that special techniques for computing disconnected diagrams (see below) can be employed, which significantly reduce the size of the numerical noise typically affecting the computation of such diagrams. In addition, with our choice of twisting for valence fermions the operator matrix elements relevant for the various flavour contents of the nucleon turn out to be all multiplicatively renormalizable. Although the same property is valid in chiral invariant (e.g. overlap) lattice formulations, twisted mass fermions are computationally much less demanding, enabling us to work with a large number of gauge configurations and at fairly big volumes and small lattice spacing.
As a last point, we want to mention that for this computation we employ gauge configurations generated with N f = 2+1+1 dynamical quarks, in which besides a mass degenerate u and d light doublet also a mass non-degenerate strange s and charm c pair is present in the sea. This will allow us for the first time to also study the charm quark content of the nucleon.
2 Computational methods

Lattice action
The lattice action used in our simulations includes as dynamical degrees of freedom, besides the gluon field, a mass-degenerate light up and down quark doublet as well as a strange-charm quark pair, a situation which we refer to as the N f = 2 + 1 + 1 setup. While in the pure gauge sector we use the Iwasaki action [29], for the fermion part twisted mass fermions are used. In particular, concerning sea quarks, we make use of the formulation of refs. [30,31] for the light mass degenerate u-d sector, while the action introduced in refs. [32,33] is employed for the mass non-degenerate c-s sector. The quark mass parameters of the heavy flavour pair have been tuned so that in the unitary lattice setup the Kaon and D-meson masses, take (approximatively) their experimental values. More information about this scheme and further simulation details can be found in ref. [34].
In order to fix the notation, we give here the explicit form of the twisted mass action for a doublet of degenerate quarks : Here µ f > 0 denotes the bare twisted mass. The hopping parameter κ f is an alias for the bare standard mass m 0f = ((2κ f ) −1 − 4)/a. Eq. 6 defines D (f ) tm [U] the flavour degenerate Wilson twisted mass operator. When κ f is tuned to its critical value, κ cr , the lattice QCD formulation known as maximally twisted mass fermions is achieved, which guarantees O(a) improvement of physical observables. The value of κ cr for all valence flavours is taken to be the same as in the sea quark sector [33].
For further needs we also introduce the operators D f,± denoting the upper and lower flavour components of D where tr denotes the trace in flavour space. D f,± [U] then is the Dirac operator of an Osterwalder-Seiler lattice quark with mass ±µ f .
When we discuss below the 2-point and 3-point correlation functions necessary for this work, we will use the so-called physical basis of quark fields denoted as ψ f . The physical field basis is related to the twisted quark field basis, χ f , by the following field rotation 1 where the twist angle ω f = π/2 at maximal twist. In the following, ψ f with index f = l, s, c will denote quark field doublets of light (l), strange (s) or charm (c) quarks depending on the mass µ f chosen in the valence sector.
Since ψ f will always refer to the physical basis we will denote with u and d the two components of ψ l . Staying close to the notation of Eq. (7) we will denote with s ± (resp. c ± ) the two components of ψ s (resp. ψ c ).
In order to have a consistent mixed action setup, the values of the bare quark mass parameters µ s and µ c in Eq. (6) have been tuned such that the Kaon and D-meson masses of the unitary setup [34] are matched.
In particular the matrix element entering σ πN will be calculated in such a way that (Euclidean) unitarity is preserved at finite lattice spacing, while those of interest for the strange and charm content of the nucleon are evaluated in a mixed action setup where unitarity violations represent mere O(a 2 ) artefacts [33]. Such cutoff effects are expected to be numerically small in line with the findings from previous mixed action studies carried out on ETMC N f = 2 [35] and N f = 2 + 1 + 1 [36] gauge ensembles.

Nucleon scalar matrix elements
The nucleon two-point function is defined in the physical quark basis by where ... denotes field average, Γ ± = (1 ± γ 0 )/2 are the parity projectors, x src ≡ (t src , x src ) is the space-time location of the source and τ ≡ t − t src stands for the source-sink separation. The subscript N refers to the proton or to the neutron state for which the interpolating fields are given by the formulae where C is the charge conjugation matrix. Note that due to translational invariance C ± N,2pts (τ ) does not depend on the spatial source location, x src , which we can thus choose freely. Let us also recall that, since here we work in the SU(2) isospin limit approximation, an exact symmetry of the action (i.e. P ×(u ↔ d) where P is parity) leads to the relation C ± n,2pt (τ ) = C ± p,2pt (τ ) [37].
The three-point functions of interest in this paper are defined by where O f for f = l, s, c denotes an appropriate (see below Eq. (13)) lattice regularization of the light, strange or charm quark scalar density and τ op = t op − t src is the operator-to-source time separation. Since we are considering an operator with a non vanishing vacuum expectation value, we also need to introduce the corresponding vacuum subtracted correlator To be specific, the operator O f will be given for our case by depending on the quark flavour of interest. Using these operators, we shall obtain the multiplicatively renormalizable, O(a) improved matrix elements relevant for this paper.
For each flavour f , the bare scalar matrix element at zero momentum transfer N(p)|O f (0)|N(p) can be written in terms of an effective coupling bare constant g S,f and the nucleon spinor u N , in the form Using the two-and three-point correlators of Eqs. (9) and (11), we build for f = l, s, c the ratio where ∆M is the mass gap between the lowest nucleon state and the first excited state with the same quantum numbers. One can thus extract from the asymptotic time behaviour of the various R f (τ, τ op ) the bare effective scalar couplings g S,f , which are in turn simply related to the nucleon sigma terms of interest. For instance, at maximal twist the lattice regulated versions of Eqs. (4) and (5) will read The systematic errors O(e −∆M |τ −τop| ) and O(e −∆M |τop| ) originating from the finiteness of the time separations τ −τ op and τ op will be neglected in this work. However they can have a non-negligible impact on the evaluation of nucleon matrix elements, as shown for instance in [38]. We are therefore planning to address this problem in a forthcoming publication.

Lattice discretization and evaluation of correlators
The main aim of the paper is to study whether the improved methods to compute disconnected diagrams as applicable for twisted mass fermions will indeed lead to a calculation of the quark contents of the nucleon with significantly reduced errors compared to earlier works. The analysis performed in this work concentrates therefore on one ensemble of a 32 3 × 64 lattice volume with a lattice spacing of a = 0.0779(4) fm (β = 1.95) where the error quoted is only statistical [39], and a pion mass of approximately 390 MeV (aµ l = 0.0055).
In order to improve the overlap between the ground state and the interpolating fields we use Gaussian smearing of the quark fields appearing in the interpolating fields. We also use APE smearing of the gauge links involved in the Gaussian smearing, following the same strategy as in [40,41].
In the twisted basis the scalar operators read and are hence given by the pseudo scalar density. There is a one-to-one correspondence between the bare operators O f and the bare operators O f introduced in Eq. (13) which is given at maximal twist by: While the two-point nucleon correlators of Eq. (9) give only rise to quarkconnected Wick contractions, in general the three-point functions of Eq. (11) yield both quark-connected (illustrated in Fig. 2a) and quark-disconnected (illustrated in Fig. 2b) contributions. In the following we will refer to them simply as to connected and disconnected fermionic Wick contractions (or diagrams) and shall write with C ±,f N,3pt (resp. D ±,f N,3pt ) corresponding to the connected (resp. disconnected) quark diagrams, defined as where the symbol [...] is a shorthand for all the connected fermionic Wick contractions. In particular, the contribution of the disconnected fermion loop to D ±,f N,3pt on a given gauge configuration U in our setup reads For the strange and charm content of the nucleon only disconnected diagrams contribute to the three-point correlator, while for the light quark content both kinds of fermionic diagrams matter. The connected contributions to N(p)|O l |N(p) have been evaluated using standard techniques for three-point functions ("sequential inversions through the sink"). In this method one needs to fix the sink-to-source separation τ = t − t src and we choose, as in ref. [40], τ = 12a corresponding in physical units to a separation of τ ≈ 1 fm.
Since, using discrete symmetries and anti-periodic boundary conditions in the time direction for the quark fields, one finds where T denotes the lattice time extent, in order to increase the signal over noise ratio we have averaged contributions related by the symmetry relations (23) and (24). In addition, we have carried out Dirac matrix inversions at a number (denoted by N src in the following) of randomly chosen source points per gauge configuration, with the goal of better exploiting the gauge field information contained in each configuration.
terms (see Eq. (16)). The renormalization pattern is thus as straightforward as for chirally invariant overlap fermions.

Numerical estimate of disconnected loops
Let us shortly sketch the variance reduction method for the evaluation of with twisted mass fermions introduced in [42,43]. The method has already been applied to study the η ′ meson in [44]. It relies on the fact that the difference between the twisted basis Dirac matrices D implying that For the practical calculation, we introduce a set Ξ of N ξ independent random volume sources, {ξ [1] , . . . , ξ [r] , . . . , ξ [N ξ ] }, satisfying where i = 1, ..., 12 refers to the spin and color indices of the source and [. . . ] Ξ denotes the average over the N ξ noise sources in Ξ .
Multiplying Eq. (26) by a Γ matrix and the noisy sources ξ * [r] (y) and ξ [r] (x) and taking the trace over spin and color indices we get For the generation of the random sources we have used a Z 2 noise taking all field components randomly from the set {1, −1}. We note that in the case of Γ = γ 5 , the quantity in Eq. (28), after summation over x ≡ x op , provides an unbiased estimator of the disconnected fermion loops of Eq. (22).

Performance of the method for disconnected loops
As a first step, we performed a study to determine the optimal number, N ξ , of stochastic volume sources to be used for evaluating the disconnected diagrams of Eq. (22). For a given number N conf of gauge field configurations increasing N ξ beyond some value will not improve the signal over noise ratio (SNR) since the noise induced by the fluctuations of the gauge fields will eventually become dominating.
For the present test we have used 677 gauge field configurations and chose fixed time separations, τ = 12a and τ op = 6a. Fig. 3 shows the SNR as a function of the number of stochastic source samples N ξ employed to evaluate the disconnected loops of Eq. (22). As indicated in the figure, for the computation of the two point function we have used N src = 1, 2, 3, 4 randomly chosen point sources per gauge configuration.
The figure demonstrates that the signal over noise ratio reaches a plateau for N ξ > 7, meaning that for larger values of N ξ the error is dominated by the gauge field fluctuations. The finding that for N ξ > 7 a plateau of the SNR is reached holds true for all values of N src we have used. In addition we have checked that the above conclusion remains essentially valid in the quark mass range we intend to explore. Fig. 3 also demonstrates that the SNR increases when more source points per gauge field configuration are used. When we change the number of source points from 1 source per configuration to 4, we find a decrease of the error by a factor of approximately ∼ 1.6. Although this does not correspond to the optimal factor of 2, using N src -values moderately larger than 1 turns out to be a convenient and economic way to increase the signal over noise ratio. In the final analysis, we will always use N ξ = 12 and N src = 4.
In Fig. 4 we compare the efficiency of the method discussed in sect. 2.4, which is based on the peculiar property, see Eq. (25), of twisted mass fermions, with another noise reduction technique relying on the hopping parameter expansion of the Dirac operator. This latter technique is not restricted to twisted mass lattice QCD and has been introduced in [45,46]. We refer the interested reader to the appendix B of [43] for an implementation in the case of twisted mass fermions. As can be seen from the figure the twisted mass specific variance reduction technique improves the signal over noise ratio by a factor ∼ 3. Performing a simple extrapolation in the number of gauge field configurations we estimate that with the hopping parameter expansion technique O(10000) configurations would be needed to reach a result 5σ away from zero, while only O(1000) configurations are necessary to obtain the same accuracy with our twisted mass specific technique.

The Pion-Nucleon σ-term, σ πN
We first concentrate on the determination of σ πN defined in Eq. (2). Since for this quantity only the up and down quarks come into play, we work in a fully unitary setup, where valence and sea quarks are regularized in the same way. In the following we will denote by R conn. (resp. R disc. ) the contribution of C +,f N,3pt (resp. D +,f N,3pt ), see Eqs. (20,21), to the ratio R l defined in Eq. (15). In Fig. 5 we show our results obtained for R l , R conn. and R disc. as functions of τ op = t op − t src for a fixed sink to source separation, τ = t − t src = 12a. R disc. has been computed using measurements over 842 configurations with N ξ = 12 randomly chosen volume sources. R conn. has been computed using 510 configurations by employing the fixed sink method.
The connected part, R conn. , denoted by the black filled circles, shows a pronounced time dependence indicating the contribution of excited states. In this work, we do not attempt to quantify the size of this systematic effect since our goal here is more to investigate whether statistically significant values for the scalar quark contents of the nucleon can be obtained. The disconnected part, R disc. , denoted by blue triangles in Fig. 5, clearly corresponds to a small contribution compared to the connected part R conn. , of the order of ∼ 10% of the full ratio, R l , represented by the red diamonds in the figure.
In Fig. 6 (a zoom of Fig. 5) we show only the disconnected contribution. Note the change in the scale on the vertical axis. It is encouraging that, by employing the techniques described above, we can indeed obtain a non-zero signal at a ∼ 4σ level. In order to determine the "plateau" values of R conn. , R disc. and R l , we performed several fits to a constant through the data varying the fit interval. The results are summarized in Table 1. We find that the disconnected contribution is about ∼ 8% of the connected one. Nevertheless since the error on the connected contributions is smaller than the value of R disc. , the disconnected contribution cannot be neglected when computing the ratio R l . We finally remark that all the statistical errors in this work are computed using the bootstrap method [47].
An estimate of the systematic error on σ πN can be given on the basis of the spread of the results one gets by varying the time interval [τ op 1 , τ op 2 ] over which   Table 1 Plateau values for the ratio R disc. , R conn. and R l relevant for the extraction of σ πN for different time intervals [τ op 1 , τ op 2 ]. We also include the χ 2 by degrees of freedom (χ 2 /ndof) and the confidence level (CL).
the plateau is taken, as displayed in Table 1. For this purpose we construct the distribution of all fit results, weighted by their confidence level, and take the variance of this distribution as our estimate of the systematic errors. Using this procedure we find where the errors correspond to statistical and systematic uncertainties, re-  spectively. Note that the dominant contribution to the systematic error comes from the connected part of the ratio R l . Also, the value obtained here corresponds to only one pion mass of about 390 MeV and an extrapolation to the physical value of the pion mass will be finally needed. Notice that chiral perturbation theory predicts that σ πN vanishes in the zero quark mass limit. While the calculation of σ πN at several quark masses and lattice spacings is beyond the scope of this paper, we remark that simulations in this directions are under way.

Strange content of the nucleon
For each value of the valence quark mass, one can define the quantity: where µ f is the bare quark mass and N|O f |N the bare matrix element corresponding to the quark flavour f . As argued in Appendix A, f T f is a RGI quantity in our mixed action setup. For the ensemble used in this study, the nucleon mass am N has been determined in [41] and is am N = 0.510 (7). In Fig. 7 we show R disc. versus τ op for a quark mass of aµ s = 0.016. We see that the ratio is ∼ 5σ away from zero in the middle of the "plateau". We have performed several fits varying the time interval to extract a "plateau" value. Results for f Ts , σ 0 and y N are summarized in Table 2. While y N does not depend strongly on the fit window and seems thus to be, within our accuracy, free of excited state contaminations, f Ts and σ 0 are affected by an excited state contamination in a similar way as we observed for σ πN . As we already remarked before, a quantitative evaluation of this systematic effect goes beyond the goal of this paper and will be addressed in the future. Our present results at a pion mass about 390 MeV for f Ts , σ 0 and y N are where the first number in parenthesis represents the statistical error and the second the systematic uncertainty. In order to estimate the magnitude of systematic effects the same strategy as in the case of σ πN has been employed. The statistical error on f Ts also includes the error on the nucleon mass determination. Note that the value of y N quoted is obtained directly from the ratio of three point functions C ±,s,sub N,3pt (τ, τ op ) and C ±,l,sub N,3pt (τ, τ op ) and agrees within error with the value of y N estimated from σ 0 and σ πN (y N = 1 − σ 0 /σ πN ≈ 0.092). We stress again that our measurement of the strange content of the nucleon leads to a value of y N which is different from zero at a 5σ level. It is interesting to compare f Ts to the value, f T l , one gets in the light quark sector, for which we obtain a significantly larger value, namely f T l = 0.117 (12)(3) (at m P S ∼ 390 MeV).  Table 2 Results of the fits to f Ts , σ 0 and y N at aµ s = 0.016. The notations are the same as in Table 1.

Charm quark content
Following the same strategy as in sect. 3.2, we have carried out the first study of the charm quark content of the nucleon. This is possible because we have at our disposal N f = 2 + 1 + 1 simulations with a fully dynamical charm quark degree of freedom.
We show in Fig. 8 the dependence of R c on τ op (at τ = 12a), using exactly the same statistics as in the light and in strange quark sectors. Unfortunately, for the charm quark content no hint of a plateau is visible. Signal and noise have equal order of magnitude and our results are compatible with zero. For comparison we also show the results for the strange quark content obtained in the previous section as a grey band. From our data we can only establish the inequality

Conclusion
In this work we have performed a benchmark calculation of the scalar quark contents of the nucleon by directly computing the matrix elements N|O f |N for f = l, s, c. Extending the calculation to strange and the charm quark flavours became possible owing to N f = 2 + 1 + 1 dynamical simulations recently carried out by the ETM Collaboration [34]. Our calculations were In evaluating these nucleon matrix elements, maximally twisted mass fermions are very helpful in two respects. The first is that the twisted mass fermion regularization provides a framework where it is possible to efficiently evaluate quark-disconnected diagrams. The second is that a consistent lattice framework can be set up where the matrix elements of interest are multiplicatively renormalizable and at the same time O(a) improved. As a result of these technical benefits we have been able to control the disconnected contributions and provide statistically significant values for σ πN = m l N|ūu +dd|N and σ 0 = m l N|ūu +dd − 2ss|N .
In the case of the scalar charm content of the nucleon, our statistics was not sufficiently large to yield a signal above the statistical noise. We could thus only give the bound | N|O c |N | | N|O s |N |.
We remark that for phenomenological applications, the relevant quantities are actually the RGI quantities m s N|O s |N and m c N|O c |N . Given our current statistical accuracy, it is therefore unclear at this moment whether the large Yukawa coupling of the charm quark can compensate the smallness of the matrix element.
The most important achievement of this paper is the rather accurate evaluation of the ratio y N = 2 N|ss|N / N|ūu +dd|N between the strange and the light quark content of the nucleon (see Eq. (3)), for which we find y N = 0.082 (16)(2). The value we obtain is small compared to estimates from chiral perturbation theory but is in line with recent lattice results obtained by other groups [21,25,27,48].
Naturally the results presented in this paper need to be further scrutinized. In particular a careful study of the unwanted excited state contamination must be carried out to reduce the magnitude of the systematic errors associated to these effects. Finally data points at various lattice spacings and pion masses are necessary to be able to safely perform an extrapolation to the continuum limit and to the physical pion mass.

A Scalar density renormalization with maximally twisted Wilson quarks
In this Appendix we want to discuss the renormalization properties of scalar quark operators. We will separately discuss the unitary and the mixed action case in the setting offered by maximally twisted Wilson fermions [30].
We work here in the so-called physical quark basis and adopt the notations of refs. [31][32][33]. We recall however that, as it is customary (see e.g. ref. [49]), the operator renormalization constants (RC's), being independent of the twisting angle (in all mass-independent schemes), are named after the form operators take in the twisted quark basis.

A.1 Unitary degenerate doublet
It has been proved in refs. [30,31] that at maximal twist 2 in the case of a degenerate u, d doublet, with µ u = −µ d ≡ µ l , quark mass and scalar density renormalize according to the formulae where Eq. (A.2) is written in the physical quark basis (with r u = −r d as a consequence of the above chosen values of µ u and µ d ) and the last term in its r.h.s. represents the mixing of the quark scalar density operator with the identity.
In this paper this term is of no importance because it will be automatically subtracted out in the computation of the nucleon matrix elements as described in the main text. For this reason, in order not to overload the forthcoming formulae, this mixing will not be indicated anymore.

A.2 Unitary non-degenerate doublet
For a pair of maximally twisted mass non-degenerate quarks, which for concreteness we name s and c, one finds [31,33] the more complicated relations where the bare mass parametersm and ǫ coincide, respectively, with the average mass and the mass difference of c and s quarks in free theory.

A.3 OS valence fermions
The mass RC of each OS valence fermion in the action is Z µ = Z −1 P . This follows from the proof provided in ref. [33] or from the extension of an old argument given in ref. [50] that we reproduce in sect. B for completeness. We thus get We recall that Z P is an even function of the r-Wilson parameter.
In the philosophy of the "mixed action" approach proposed in ref. [33], a pair of mass degenerate OS fermions, denoted (in the so-called physical basis) as s + and s − , with opposite values of the Wilson parameter (r s + = −r s − = 1) is introduced to represent the valence s quark, with the understanding that no Wick contractions between the fermion s + and the fermion s − is allowed. This is done also in sect. 2.2, with µ s > 0 denoting the bare s quark mass.
Consider a correlator where besides the strange quark scalar density only (renormalized) operators containing no strange quark are present. Then the insertion of the renormalized combination (we recall that the divergent mixing with the identity must be subtracted out) is finite, i.e. no new divergences are introduced. The reason for this fact can be traced back to the cancellation of chiral violating effects (coming from "quark disconnected" -i.e. OZI [51][52][53] violating -diagrams) between the two self-contractions ("loops") of the two valence quarks regularized with opposite values of r. Alternatively this result can be ascribed to the fact that, having the members of the s + , s − pair opposite values of the r parameters they look like a mass degenerate (valence) flavour doublet (in the "physical" quark basis), e.g. just as the mass degenerate u and d pair discussed above. Naturally it remains the fact that the theory is not unitary since valence and sea quarks are regularized differently. This lack of unitarity leaves behind only O(a 2 ) effects [33].

B Mass renormalization for OS valence fermions
For completeness in this appendix we want to explicitly prove the relation (A.7) along the lines of ref. [50] in a setting where N v ≥ 2 valence OS fermions are introduced over an SU(N f = 2) (or SU(N f = 2 + 1 + 1)) maximally twisted sea. The renormalization condition (A.7) is valid for anyone of the N v OS fermions, being Z P the renormalization constant of the non-singlet pseudo-scalar quark density.
The key observation to prove Eq. (A.7) is to recall that in order to be entitled to use the techniques that are usually employed to derive WTIs, it is necessary to have a fully local formulation of the theory. This means that for the purpose of dealing with a mixed action case, one has to keep in mind that for each OS valence quark in the action a corresponding ghost with equal mass (and opposite statistics) has to be introduced. Without loss of generality, for the purpose of proving Eq. (A.7), we can assume that all valence quarks (and ghosts) have the same bare (twisted) mass, µ OS . The generalization to the non-degenerate valence quarks is straightforward.

WTIs for OS fermions
Let A a µ , a = 1, 2, . . . , N 2 v − 1 be the (non-singlet) axial vector current constructed in terms of only valence fermions (and no ghosts) and let us assume that we are exactly at maximal twist. For convenience we shall work in the twisted fermion basis where the valence quark mass term has the expression The argument can be split into four parts I -The (bare/lattice) axial WTI between on-shell states reads [54] α|∇ µ A a µ (0)|β = 2µ OS α|S a (0)|β − α|X a (0)|β , and X a is the chiral variation of the OS Wilson term, the explicit expression of which we do not need in this discussion. The only thing we need to know about X a is its mixing pattern (dictated by dimensional argument and the symmetry P χ × (µ OS → −µ OS ), where P χ is the formal parity acting on the twisted fields, see ref. [30]) which reads is immediately obtained by setting Our aim is to prove the relation from which the formula Z −1 P = Z µ (B.10) follows.
II -To this end we need to extend the previous equations to the case where the divergence of the axial current is inserted together with the singlet pseudoscalar densityP 0 = f (χ f iγ 5 χ f + ghosts) (B.11) which results from considering contributions coming from valence quarks as well as the associated ghosts. We note that as the axial current we are considering is only made up of valence quarks, it cannot rotate the ghost fields. We thus get for the WTI where the operatorP 0 is inserted α|∇ µ A a µ (x)P 0 (y)|β = 2 α|S a (x)|β δ(x − y) + +2µ OS α|S a (x)P 0 (y)|β − α|X a (x)P 0 (y)|β , (B.12) with external states such that one does not get identically vanishing matrix elements. Introducing the decomposition (B.4), we first rewrite the previous equation in the form α|Z A ∇ µ A a µ (x)P 0 (y)|β = 2 α|S a (x)|β δ(x − y) + +2(1 − η X )µ OS α|S a (x)P 0 (y)|β − α|X a (x)P 0 (y)|β . (B.13) We remark that, whenX a is inserted with a local operator, it gives rise to localized terms plus genuinely O(a) contributions. Simple symmetry considerations and the fact that in X a only valence fermions (and not ghosts) appear, imply α|X a (x)P 0 (y)|β = 2c X α|S a (x)|β δ(x − y) + O(a) , (B.14) where c X is again a finite function of g 2 0 . Substituting into Eq. (B.13) and neglecting irrelevant (for the present argument) O(a) terms gives where we have multiplied both members by ZP 0 and used Eqs. (B.7) and (B.8). It must be stressed that consistency with continuum WTIs (universality) at vanishing µ OS [54,55] requires the identification which would finally prove Eq. (B.9). The reason for the validity of Eq. (B.21) is that, when one considers the diagrams contributing to ZP 0 , one realizes that the OZI-violating diagrams where a valence quark is self-contracted (closed into a "loop") is exactly cancelled by the contribution where the corresponding ghost, present inP 0 , is closed into a "loop". One is thus only left with diagrams in which neither valence nor ghost self-contractions appear, hence exactly with the diagrams that contribute to the non-singlet pseudo-scalar density RC, Z P , where such self-contractions are forbidden by flavour conservation.