First moment of the flavour octet nucleon parton distribution function using lattice QCD

We perform a lattice computation of the flavour octet contribution to the average quark momentum in a nucleon, $\la x\ra^{(8)}_{\mu^2 = 4 \gev^2}$. In particular, we fully take the disconnected contributions into account in our analysis for which we use a generalization of the technique developed in \cite{Dinter:2012tt}. We investigate systematic effects with a particular emphasis on the excited states contamination. We find that in the renormalization free ratio $\frac{\la x \ra^{(3)}}{\la x \ra^{(8)}}$ (with $\la x \ra^{(3)}$ the non-singlet moment) the excited state contributions cancel to a large extend making this ratio a promising candidate for a comparison to phenomenological analyses. Our final result for this ratio is in agreement with the phenomenological value and we find, including systematic errors, $\frac{\la x \ra^{(3)}}{\la x \ra^{(8)}} = 0.39(1)(4)$.


Introduction
Computing nucleon properties from first principles using lattice QCD is a long standing challenge. In particular, moments of parton distribution functions (pdfs) are important benchmark quantities for lattice calculations and provide insight into the structure of hadrons. The computation of various moments of pdfs is therefore a very active research area for lattice calculations, see [2][3][4][5][6][7][8][9][10][11] for recent results. In the work we perform here a first lattice calculation of the octet contribution to the average quark momentum in a nucleon, x (8) µ 2 is presented, where µ 2 denotes the renormalization scale. This quantity is of interest by itself and can be extracted from phenomenological analyses of data from deep inelastic scattering experiments, see below. It therefore can serve as an additional quantity to probe QCD and the structure of hadrons, see [11].
In addition, x (8) µ 2 needs the same renormalization constant as the iso-vector averaged quark momentum in the nucleon x  the renormalization constant cancels. Although we can evaluate non-perturbatively the renormalization constant, its cancelation in the ratio eliminates any uncertainty related to its determination. Comparing the ratio to a phenomenological analysis can thus help to understand, whether renormalization effects can play a role in the presently observed discrepancy between lattice calculations of x (3) µ 2 and phenomenological determinations of this quantity.
Despite the above given motivations to compute x (8) µ 2 there is so far no value of this quantity available from a lattice determination. The reason for this is that x (8) µ 2 is very difficult to calculate since it involves dis-connected, singlet contributions. The definition of the flavour octet moment x (8) µ 2 reads : x (8) where q(x, µ) denotes the sum of the parton distribution function of the quark q and antiquarkq. Using the same convention, we also define the non-singlet contribution x At a value of the renormalization scale of µ 2 = 4GeV 2 the two moments x (3) µ 2 =4 GeV 2 and x (8) µ 2 =4 GeV 2 can be extracted phenomenologically using parton distribution functions determined from deep inelastic scattering data, and read, using the ABM12 pdfs set [12] and the analysis of ref. [13] x (3) (7). (1. 3) The quantities x µ 2 and x (8) µ 2 are related to matrix elements of local operators that can be computed in Euclidean space-time and are hence accessible to lattice QCD calculations. Introducing where λ a are the Gell-Mann matrices acting on a three flavour quark field ψ = (u, d, s).
With ← → D we denote the symmetrized covariant lattice derivative and with the curly brackets the symmetrization and subtraction of the trace. The required matrix element is then given by From the definition in eq. (1.1) it is clear that the octet matrix element N (p, s)|O a=8 {µν} |N (p, s) involves dis-connected contributions which have in general a bad signal to noise ratio, require thus a very high statistics and are consequently very difficult to compute on the lattice. In addition, the dis-connected contribution to x (8) µ 2 is an SU (3) flavour breaking effect and is thus expected to be small. In [1] and [14] we have demonstrated that for the here used (maximally) twisted mass lattice discretization of QCD there are special noise reduction techniques which can help substantially to improve on the signal to noise ratio. Indeed, these techniques allowed us to compute a number of quantities which were very difficult to access before [15,16].
In this work we present a generalization of the particular technique used in the context of the determination of the σ-terms [1] to calculate non-perturbatively the disconnected contribution relevant for x (8) µ 2 . As we will see below, this generalization will indeed provide a statistically significant signal for x (8) µ 2 . Note that contrary to the case of the σ-terms, the formalism to compute dis-connected contributions developed here is not limited to the twisted mass formulation. The technique could also be applied to other flavour octet operators which would allow to determine non perturbatively moments of polarized pdfs or the flavour octet axial coupling of the nucleon g (8) A for instance, see e.g. [14,16]. The paper is organized as follow : after describing the basic ingredients of our computation we give the details of the variance reduction technique used in this work. We then present a study of the systematic effects appearing in our calculation and perform an extrapolation of the results to the physical pion mass in order to compare with the phenomenologically obtained values.

Simulation Details
The lattice action used in our simulations includes as active degrees of freedom, besides the gluon field, a mass-degenerate light up and down quark doublet as well as a strange-charm quark pair in the sea, a situation which we refer to as the N f = 2 + 1 + 1 setup. We use the Iwasaki action [17] for the pure gauge action and the twisted mass fermion formulation for the Dirac action. In particular, we make use of the formulation of refs. [18,19] for the light mass degenerate u-d sector, while the action introduced in refs. [20,21] 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 approximately their experimental values. More information about the N f = 2 + 1 + 1 setup scheme and further simulation details can be found in ref. [22,23]. For the results we will show here, we have employed two values of the lattice spacing determined using the nucleon mass in [2]. They read a ≈ 0.082 fm (β = 1.95) and a ≈ 0.064 fm (β = 2.1). In addition, we will use a number of quark masses corresponding to pion masses in the range of 300 MeV − 500 MeV. The parameters of the ensembles used in this work are summarized in Table 1 In order to fix the notation, we introduce the twisted mass lattice Dirac operator D f,tm for a doublet of mass degenerate quarks : (2.1) Here D W [U ] is the Wilson Dirac operator, µ f denotes the bare twisted mass and τ 3 is the third Pauli matrix. For further needs we also introduce the operators D f,± denoting the upper and lower flavour components of D f,tm [U ], referred to as the Osterwalder-Seiler (OS) lattice Dirac operator: where tr denotes the trace in flavour space. We will call D f,± [U ] the lattice Dirac operator of an Osterwalder-Seiler 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, where the twist angle ω f = π/2 at maximal twist. In addition, ψ f with index f = l, s will denote quark field doublets of light (l) or strange (s) 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 . Following the notation of Eq. (2.2) we will denote with s ± the two components of ψ s . Employing the OS Dirac operator in the valence sector for the strange quark leads to a mixed action where the strange OS quark mass has been tuned to match within errors the unitary Kaon mass.

Nucleon matrix elements
The nucleon two-point function is defined in the physical basis by where the source position is fixed to 0 in order to lighten notations and t thus denotes the source-sink separation. We also introduced the parity projectors Γ ± = (1 ± γ 0 )/2. The subscript N refers to the proton or to the neutron state for which the standard interpolating fields are given by the formulae: Note that using discrete symmetries and anti-periodic boundary conditions in the time direction for the quark fields, one can show that C + N,2pts (t) = −C − N,2pts (T − t). Let us also recall that an exact symmetry of the action leads to the following relation at a finite value of lattice spacing, C ± n,2pts (t) = C ± p,2pts (t) [24]. In order to improve the overlap between the ground state and the interpolating fields of the nucleon 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.
The nucleon three-point functions is given by {44} is one of the twist-2 operators introduced in Eq. (1.4), t op is the time of insertion of the operator, and t s denotes the so-called source-sink separation. Note that the precise definition of the strange quark field entering in the operator O a=8 {44} is postponed to the next subsection.
Using the two-and three-point correlators of Eqs. (2.4) and (2.6), we construct the following ratio: where am N is the nucleon mass in lattice units and ∆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 R a=8 (t s , t op ) the bare quantity x (8) bare and correspondingly x bare .

Lattice evaluation
While the light quark fields used in the operator O a=3,8 {44} are the unitary fields we use, as mentioned already above, a different action for the valence strange quark. In practice we introduce a doublet of mass degenerate quarks with a mass aµ s tuned to reproduce the unitary Kaon mass. This procedure introduces an error due to the uncertainty on the determination of the matching mass that we will discuss later on but will allow us to use an efficient noise reduction technique that will be explained in the next section.
Consider the following operator in terms of the field in the twisted mass basis : Performing the rotation to the physical basis, we obtain : ) Note that J 8 keeps the same form in the two bases and that J 8 is only one possible choice for a discretization of the operator O a=8 {44} . While the two-point nucleon correlators of Eq. (2.4) give only rise to quark-connected Wick contractions, in general the three-point functions of Eq. (2.6) yield both quark-connected (illustrated in Fig. 1a) and quark-disconnected (illustrated in Fig. 1b) 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 ±,a N,3pt (resp. D ±,a 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. Note that for a = 3, the disconnected part is a O(a 2 ) effect which vanishes in the continuum limit and can thus be neglected. Introducing the contribution of the disconnected fermion loop to D ±,a N,3pt on a given gauge configuration U in our setup reads 14) The connected contributions C ±,a N,3pt (t s , t op ) have been evaluated using standard techniques for three-point functions (sequential inversions through the sink), see e.g. ref. [2].

Estimation of disconnected loops
We describe here the generalization of the variance reduction method for twisted mass fermions introduced and discussed in [1,14,25,26].
Consider the identity For the generation of the random sources we have used a Z 2 noise taking all field components randomly from the set {1, −1}.
Note that δ (µ,µs) ± , and by construction its variance, is proportional to the mass difference µ l − µ s and vanishes on each configuration in the limit µ s → µ l . Our approach, thus, exactly encodes the fact that the disconnected contributions we are interested in vanish in the SU (3) flavour limit.

Renormalization
The renormalization of the operator O a=3 is known to be multiplicative from our previous work [2] and have been obtained non perturbatively using the methodology developed in [27]. The renormalization factor Z µµ DV (β) read : Note that an independent calculation performed in [28,29] gives compatible results. In the limit µ s = µ l = 0, SU (3) flavour is an exact symmetry of the action and the operators O a=3,8 {µν} belong to the same flavour multiplet. They thus share the same renormalization pattern in a mass independent scheme. The ratio x (3) x (8) is thus renormalization free.

Results
As a first step, we have investigated the magnitude of the stochastic noise introduced by the method described in section 2.3. To this end, we used a fixed number of gauge configurations for a gauge field ensemble at coupling β = 1.95 and twisted mass parameter aµ l = 0.0055. We show in Fig. 2, the ratio R a=8 disc. for a fixed source-sink separation of t s ∼ 1 fm as a function of t op for N ξ = 6 and N ξ = 12. As can be seen, the signal is compatible with zero within the errors when using N ξ = 12. Furthermore we observe that the error does not decrease much when the number of stochastic sources is doubled. This means presumably that the error is already dominated by the intrinsic noise of the gauge field fluctuations. We nevertheless decided to use N ξ = 12 throughout this benchmark study of the octet moment.
In principle, it would therefore be possible to reduce the error further by using more stochastic noise vectors. However, the contribution of the disconnected 3-point function is small compared to the value of the connected part as demonstrated in Fig. 3 where we show the connected contribution R a=8 connected (t s = 12a, t op ), the disconnected contribution R a=8 disconnected (t s = 12a, t op ) and the full correlator R a=8 full (t s = 12a, t op ). In particular, the statistical errors on the connected part are about 0.016 (∼ 2.5%) which is significantly larger than the size of the disconnected contributions. Despite this fact, which would make neglecting the dis-connected contributions tempting, we always include them in the following analysis. Note that the R a=8 disconnected is proportional to the difference between the light and strange quark mass. Since we are using a mixed action setup, our results depends on an approximate tuning of the strange quark mass. In all the figures we use aµ s (β = 1.95) = 0.018 for β = 1.95 and aµ s (β = 2.10) = 0.015 for β = 2.10.f Those values can be compared to the values obtained for instance in [30] where the strange quark mass has been determined using the Ω − mass and correspond to aµ s (β = 1.95) = 0.0194 and aµ s (β = 2.10) = 0.0154. However, by using data on the β = 1.95 ensemble we have explicitly checked that changing the bare strange quark mass by more than ∼ 10% does not lead to any significant change in the value of the disconnected contribution.
connected (black triangles) and of their sum, R

Excited states contamination
In order to investigate the contamination of excited states due to the second and third term in Eq. (2.7), we used the same procedure ("open-sink" method) as in [31], namely we study the source-sink dependence for a fixed source to operator separation t op ∼ 0.9 fm ( t op = 11a). Details on the technical implementation of the "open-sink" method can be found in [31]. To this end, a large statistics for the connected part (∼ 23000 measurements) has been used. We plot in Fig. 4 the resulting renormalized ratios R a=3,8 (t s , t op = 11a). The results in the isovector case are represented by orange squares and the results in the flavour octet case (including the disconnected piece) are depicted by red dots. We also represent using blue triangles the disconnected contribution to R a=8 (t s , t op = 11a). As can be seen, the noise for the dis-connected part dominates for large source-sink separation. Nevertheless, we can obtain a reasonable signal up to source-sink separation of about 16a (≈ 1.3 fm). The gray bands indicate the values of x (a=3,8) obtained using a fixed source-sink separation calculation with t s ∼ 1 fm. Note that the fitting ranges have been determined choosing the fit with the longest plateau with a confidence level of least 90%. In the flavour octet case we observe that results obtained for t s /a > 15 and a fixed sourceoperator separation t op = 11a are compatible with the result obtained at a fixed source- sink separation of t s /a = 12. We mention that at the here used value of the pion mass of m PS ≈ 370 MeV, volume (L = 2.6 fm) and lattice spacing (a ≈ 0.082 fm) we found for the isovector channel a relative shift of about 10% for a corresponding comparison at different source sink separations [31]. We performed a similar study of the ratio of the octet to the singlet 3-point functions of Eq. (2.10), The ratio R(t s , t op = 11a) is shown in Fig. 5 as a function of the source-sink separation t s .The constant fit for this ratio obtained with a fixed source to operator time t op = 11a is depicted by a gray band. We used the same criterion as in Fig. 4 to determine the fitting range. In this case we do not observe any significant source-sink dependence of our results and thus there seem to be only a small excited states contamination, at least within the size of our errors. Note that the statistical errors stemming from the disconnected contribution are responsible for most of the total statistical error at large source-sink separation. In the next section we will use the difference between the central value at t s ∼ 1.3fm with the open think method and the central value obtained with the fixed source-sink method at t s ∼ 1fm as an estimate of the systematical error for x (a=3,8) and x (3) x (8) .

Chiral behaviour
We plot our results for x (3,8) µ 2 =4 GeV 2 as a function of the pseudoscalar meson mass m 2 PS in Fig. 6. The results for x (3) µ 2 =4 GeV 2 (respectively x (8) µ 2 =4 GeV 2 ) are shown using red down triangles (resp. blue up triangle) for a lattice spacing of 0.082 fm and by an orange square (resp. an orange dot) for a lattice spacing of 0.064 fm. Neglecting, in a first step, excited state contaminations and in order to investigate the chiral limit behaviour of our data, all the results have been obtained using a fixed source-sink separation of approximately 1 fm. Of course, in our final result we will add the shift from the excited states contamination as a systematic error. As can be seen in Fig. 6 the results obtained at two different lattice spacings are compatible, and O(a 2 ) effects can be neglected.
The phenomenological estimates [13] of Eq. (1.3) are represented by two black stars. As can be seen in the graph, for both, x (3) and x (8) and for unphysically large values of the pion mass the lattice data lay consistently above the phenomenological values, a phenomenon that is well know from previous investigations of x (3) , see e.g. [11], and is here demonstrated for the first time for x (8) .  As explained before, systematic errors stemming from discretization effects and from the matching mass dependence can be neglected. With the present data set we cannot estimate safely the volume dependence of x (3,8) , however knowing that they are negligible in the isovector case in our N f = 2 calculation [3] we assume them to be small and neglect them, too. Our dominating source of systematic error for the individual matrix elements is then due to the excited states contamination. In addition, as the chiral extrapolation is not performed here, it remains as a presently unquantifiable systematic error.
In fig. 7 we show the ratio x (3) x (8) as a function of the pion mass together with the phenomenological value. Here the situation is somewhat different in that the lattice data agree with the phenomenological analysis even at the here used unphysically large pion masses. This suggests that some of the systematic uncertainties cancel in the ratio. As discussed in section 3.1, the systematic error stemming from the excited states contamination is negligible for the ratio x (3) x (8) . Therefore, we conclude for the individual moments that either the non-perturbative renormalization or the chiral extrapolation is the most probable systematic effect leading to the observed discrepancies in x (3,8) . It would therefore be highly desirable to perform simulations directly at the physical point to eliminate the uncertainty from the chiral extrapolation. Work in this direction is in progress [32,33].
Since the lattice data are flat as function of the pion mass, we perform a simple constant extrapolation to the physical pion mass. As a systematic error, we take the difference between the data points at the smallest and the largest pion mass. We find then for x (a=3) where the first error is statistical and the second one is our estimate of the systematic error. When comparing to phenomenological extractions [13] we find that the ratio x (3) x (8) is, within errors, compatible with phenomenological extractions at least for the here used simple extrapolation to the physical point.

Conclusion
In this work we have performed a benchmark computation of the flavour octet combination of the first moment of parton distribution functions of the nucleon using N f = 2 + 1 + 1 twisted mass fermions tuned to maximal twist. This quantity can provide a first hint of the contribution of the strange quark to the quark momentum in the nucleon. We have utilized the fact that x (8) µ 2 shares the same renormalization property as the standard isovector moment x (3) µ 2 . Using our result for the non-perturbative renormalization constant obtained in [2] we can thus provide an O(a 2 ) improved estimate of x (8) µ 2 . The calculation takes into account the notably noisy dis-connected contribution, which are shown to have a negligible contribution in our setup.
In our investigation of x µ 2 [31], that the excited states contamination can be substantial and has to be taken into account as a systematic error. We also found that, similar to the case of x   However, when the ratio x (a=3) µ 2 =4 GeV 2 / x (a=8) µ 2 =4 GeV 2 is considered, the lattice data show an overall agreement with the corresponding phenomenological ratio, indicating that in the ratio systematic errors cancel. Since in the ratio the excited states contamination is small, the here observed agreement cannot be due to this systematic uncertainty. One reason why the ratio agrees with the phenomenological value can be that the renormalization constants cancel. Another possibility is that both, x (8) µ 2 and x (3) µ 2 have a very similar chiral limit behaviour such that the effects from the chiral extrapolation cancel. It would be highly desirable to have therefore lattice results directly at the physical point for both quantities. Since the lattice data are rather flat, we performed a simple constant extrapolation to the physical point for the ratio x (a=3) µ 2 =4 GeV 2 / x (a=8) µ 2 =4 GeV 2 and find x (3) x (8) = 0.39(1)(4) (4.1) This value is in agreement with the phenomenological extraction of this quantity from deep inelastic scattering data [13].