Flavor decomposition of the nucleon unpolarized, helicity and transversity parton distribution functions from lattice QCD simulations

We present results on the quark unpolarized, helicity and transversity parton distributions functions of the nucleon. We use the quasi-parton distribution approach within the lattice QCD framework and perform the computation using an ensemble of twisted mass fermions with the strange and charm quark masses tuned to approximately their physical values and light quark masses giving pion mass of 260 MeV. We use hierarchical probing to evaluate the disconnected quark loops. We discuss identification of ground state dominance, the Fourier transform procedure and convergence with the momentum boost. We find non-zero results for the disconnected isoscalar and strange quark distributions. The determination of the quark parton distribution and in particular the strange quark contributions that are poorly known provide valuable input to the structure of the nucleon.

The flavor decomposition of proton charges and form factors had been under investigation in the last few years with calculations of disconnected diagrams using ensembles at or near the physical values for the quark masses(physical point) [11][12][13][14][15][16][17][18][19][20][21][22][23][24]. Such a success in lattice calculations of hadron structure is partly due to available computational resources, but also due to novel methodologies to extract the individual quark Mellin moments and form factors. A notable example is the hierarchical probing [25], which improves the signal significantly. In this work, we extend hierarchical probing to non-local operators. Such an approach was shown to be successful in our first calculation on the helicity PDF [26].
Matrix elements of non-local operators are of great interest in recent years. These can be related to light-cone PDFs through a factorization and matching procedure. Methods to access the x dependence of PDFs, such as the quasi-PDFs [27,28], pseudo-ITDs [29], and current-current correlators [30,31] are now well established. Progress in terms of the renormalizability, renormalization prescription, and factorization of light-cone PDFs has been made, alleviating major sources of systematic uncertainties. For their application in lattice QCD see the recent reviews of Ref. [32][33][34]. Results from lattice QCD simulations on the x-dependence of PDFs are very promising, and therefore, their flavor decomposition is the extension of these investigations. The matching for the singlet case, as well as the mixing with the gluon PDFs have been recently addressed [35,36]. In Ref. [26] we presented the first calculation for the helicity PDFs including disconnected contributions and the flavor decomposition for the up-, down-and strangequark PDFs. Here we extend the calculation to the three types of collinear PDFs, that is the unpolarized, helicity and transversity PDFs. While such calculations are becoming feasible, there are a number of computational challenges before taming the statistical uncertainties. To date, calculations of disconnected contributions for matrix elements of non-local operators at the physical point do not exist.
The paper is organized as follows: Section II presents the methodology of extracting the nucleon matrix elements, the renormalization and matching. Section III focuses on the details of the calculation and the techniques employed. The results for the disconnected and connected matrix elements are presented in Sections IV and V, respectively. In Section VI, we extract the vector, axial and tensor charges using the data with zero length for the Wilson line. The PDF reconstruction and flavor decomposition is given in Section VII. Finally, we give out conclusions in Section VIII.

A. Nucleon bare matrix elements
The main component of this study is the calculation of the nucleon matrix elements of non-local operators, that is where |N (P 3 ) is the nucleon state with momentum boost along the z-direction, i.e. P = (0, 0, P 3 ). The fermionic field ψ(x) ≡ ψ( x, t) can be either the light quark doublet ψ ≡ (u, d) T or the strange quark field ψ(x) ≡ s(x). The Wilson line W (z) is constructed in the direction parallel to the nucleon boost P and extends from zero length to up to half of the lattice, L/2, in both positive and negative directions. The Dirac structure of the operator, Γ, acts in spin space and depends on the type of the collinear PDF under study. Without loss of generality, one can take that the momentum boost is in the z-direction (k). Based on this, we employ the following Γ matrices: • Γ = γ 0 for the unpolarized distribution q(x); • Γ = γ 5 γ 3 for the helicity distribution ∆q(x); • Γ = σ 3j with j = 3 for the transversity distribution δq(x).
We calculate both the isovector and isoscalar quark contributions to Eq. (1), and, therefore, we introduce the superscript Υ for the matrix elements. The matrix Υ acts on the light quark sector and takes the value τ 3 (1) for the isovector (isoscalar) distribution, where τ 3 = diag(1, −1), is the third Pauli matrix.
The matrix elements are computed from the ratio of three-and two-point functions defined as where t s is the source-sink separation, and t ins the insertion time of the three-point function. We use the proton interpolating field N α = ε abc u a α (x) d bT (x)Cγ 5 u c (x) with C = γ 0 γ 2 . The three-point function projector depends on the operator under study and can be found in Table I and for the two-point function we use P = (1 ± γ 0 )/2. To increase the number of measurements for the disconnected diagrams, we average the three-and two-point functions over plus and minus parity projectors. The operator O is defined as O( y, t ins ; z) =ψ( y + zẑ, t ins )Υ Γ W ( y + zẑ, y)ψ( y, τ ) , and is inserted in the three-point function of Eq. (2). The Wick contractions lead to two topologically different diagrams, as shown in Fig. 1b and 1c. In Fig. 1a, we also show pictorially the two-point function of Eq. (2). For the case that the fermionic field in Eq. (3) is ψ = (u, d) T and Υ = τ 3 , we obtain the matrix elements for the isovector distribution u − d, which receive contribution from the connected diagram only (Fig. 1b). However, in the case where Υ = 1, the three-point function takes contributions from both connected and disconnected diagrams. For the nucleon, the strange-quark contribution comes exclusively from the disconnected diagram. We emphasize that, disconnected contributions have a considerably smaller signal-to-noise ratio compared to the connected ones and their evaluation requires the use of stochastic and gauge noise reduction techniques described in detail in Sec. III A 2.
The ratio of three-over two-point functions becomes, and is used to obtain the matrix elements of Eq. (1). This relation is meaningful for the ground-state contribution. To isolate the latter, we apply constant (plateau) fits in a range where the operator insertion time is large enough, and away from the source and sink. Besides the plateau fit method, we employ different techniques allowing the extraction of the nucleon matrix element, as described in Sec. IV A.

B. Non-perturbative renormalization
The bare matrix elements of Eq. (1) must be renormalized in coordinate space prior to obtaining the quasi-PDFs, which are defined in the momentum space (see Sec. II C). The renormalization of both the non-singlet and singlet quantities is multiplicative. In this work, we use the non-singlet renormalization function, as the difference with the singlet is expected to be small [37], which was demonstrated numerically for local operators [18,20]. The mixing between the unpolarized and helicity singlet-quark PDFs with the gluon PDFs arises at the matching level because there is no additional non-local ultraviolet divergence in the quasi-PDF [35,36,38].
To renormalize the matrix elements we apply the regularization independent (RI ) scheme, and use the momentum source method [39] that offers high statistical accuracy. More details on the setup can be found in Refs. [40,41]. In Refs. [42,43], we proposed an extension of the renormalization prescription to include non-local operators, which we also follow in this work. The conditions for the renormalization functions of the non-local operator, Z Γ , and the quark field, Z q , are V(p, z) (S(p)) is the amputated vertex function of the operator (fermion propagator) and S Born (p) is the tree-level of the propagator. These conditions are applied at each value of z separately. We improve Z q , and consequently Z Γ , by subtracting lattice artifacts calculated to one-loop level in perturbation theory and to all orders in the lattice spacing, dZ ∞ q (p). The details of the calculation can be found in Ref. [40]. The RI-type schemes are mass-independent, and therefore Z Γ is calculated at several ensembles with different values for the quark masses. Eventually, a chiral extrapolation is applied to remove residual artifacts related to the quark mass. Here we use five ensembles that are generated at different pion masses in the range of 350 MeV -520 MeV with a lattice volume of 24 3 × 48. For the chiral extrapolation to be meaningful, all quark masses should be degenerate. Therefore, we use N f = 4 ensembles generated by the Extended Twisted Mass collaboration (ETMC) that are dedicated to the renormalization program. These ensembles have the same lattice spacing and action parameters as the N f = 2 + 1 + 1 ensemble used for the production of the nucleon matrix elements.
Z Γ depends on the RI renormalization scale µ 0 , and it will be converted to MS and evolved at a scale of choice. To reliably perform this procedure we use several values of µ 0 , and the conversion and evolution formulas is applied on Z Γ obtained at each scale. We choose the initial scale µ 0 such that discretization effects are small [40]. In particular, the 4-vector momentum p, which is set equal to µ 0 , has the same spatial components: p = (p 0 , p 1 , p 1 , p 1 ). The values of p 0 and p 1 are chosen such that the ratio p 4 (p 2 ) 2 is less than 0.28, as suggested in Ref. [44]. The values of a µ 0 cover the range [1,5]. For each µ 0 value, we apply a chiral extrapolation using the fit to extract the mass-independent Z RI Γ,0 (z, µ 0 ) at each value of the initial scale. Z RI Γ,0 (z, µ 0 ) is converted to the MS scheme and evolved to µ=2 GeV (µ= √ 2 GeV) using the results of Ref. [42] for the unpolarized and helicity (transversity) PDFs. The conversion and evolution depends on both the initial scale µ 0 and the final scale in the MS. The appropriate expressions have been obtained to one-loop perturbation theory in dimensional regularization. Therefore, there is residual dependence on the initial scale µ 0 , which is eliminated by taking the limit (a µ 0 ) 2 → 0 using a linear fit on the data in the region (a µ 0 ) 2 ∈ [1 − 2.6].
The last step of the renormalization program is the conversion to a modified MS scheme (MMS), developed in Ref. [41], and is given by where This scheme was introduced to satisfy particle number conservation. In Ref. [41] we showed that the difference between MS and MMS is numerically very small, but brings the PDFs closer to the phenomenological ones. In Eq. (9) µ F is the factorization scale set equal to the MS scale. Ci, Si, Ei and sgn are the special functions cosine integral, sine integral, exponential integral, and sign function, respectively. The coefficients e (i) Γ depend on the operator: {e Γ , e Due to the presence of the Wilson line, both the matrix elements and renormalization functions are complex functions. As a consequence, a complex multiplication is required to extract the renormalized matrix element, that is For simplicity in the notation, the dependence on z, P 3 , scheme and scale is implied. As can be seen, the real (imaginary) part of the renormalized matrix elements are not simple multiples of the real (imaginary) part of the bare matrix element. Therefore, controlling systematic uncertainties in the renormalization is an important aspect of the calculation.

C. Quasi-PDFs and matching to light-cone PDFs
Quasi-PDFs are defined as the Fourier transform of the renormalized nucleon matrix elements in Eq. (4) with respect to the Wilson line length zq Note that the renormalized matrix elements, h Υ Γ , depend on the scheme and scale, which also propagates toq(x, P 3 ). For simplicity in the notation, this dependence is implied. As mentioned in the previous paragraph, the matrix elements are renormalized in the MMS scheme and evolved to 2 GeV ( √ 2 GeV) for the unpolarized and helicity (transversity) PDFs.
On the lattice, we can only evaluate the matrix elements for discrete and finite values of z. Therefore, the integral of Eq. (11) is replaced by a discrete sum over a finite number of Wilson line lengths where dz/a = 1. For the summation in Eq. (12) to accurately reproduce Eq. (11), both the real and imaginary parts of the matrix element should be zero beyond z max . Practically, this is not always possible due to the finite momentum boost and limited volume in the lattice formulation. The choice of the cutoff z max , which is anyway limited up to L/2, requires an extensive study. We note that systematic effects related to the reconstruction of the PDFs is operator dependent, as each matrix element may have different large-z behavior. We will show the results of such analysis in Sec. VII B.
From the finite-momentum quasi-PDF, it is possible to obtain the light-cone parton distribution function (infinite momentum) through the so-called matching procedure. This is accomplished through a convolution of the quasi distribution with a kernel evaluated in continuum perturbation theory within the large momentum effective theory (LaMET) [33,45]. The matching formula reads and the factorization scale µ is chosen to be the same as the renormalization scale. The matching kernel C contains information on P 3 which, in principle, is eliminated in q(x, µ). However, there is residual P 3 due to the limitations in accessing large values of momentum from lattice QCD (see, e.g., Ref. [41]) and the matching kernel being available to limited order in perturbation theory. Most of the calculations of the matching kernel have been performed to one-loop level (see, e.g., Refs. [30,[46][47][48][49][50][51][52]). Recently, the computation of the kernel C was extended to two loops [52][53][54][55]. In this study, we employ the kernel in the MMS-scheme which is known at one-loop level [41]. This matching kernel relates the quasi-PDFs defined in the MMS scheme at some scale, to the light-cone PDFs in the MS at the same scale.
For the unpolarized and helicity we choose a scale of 2 GeV, while for the transversity we choose √ 2 GeV.
To calculate the anti-quark distributions from q(x), we exploit the crossing relations [56], that is

III. LATTICE SETUP
The computation is performed using one gauge ensemble of N f = 2 + 1 + 1 clover-improved twisted mass fermions and the Iwasaki improved gluonic action [57] generated by ETMC [58]. The fermionic action of the light quarks in the "twisted basis" takes the form Here, χ T l (x) = (u, d) is the light quark doublet in the twisted basis, D W [U ] is the massless Wilson-Dirac operator and F µν [U ] is the field strength tensor. The last term is weighted by c SW , the Sheikoleslami-Wohlert [59] clover coefficient. The heavy quark twisted mass action is similar to the light quark action in Eq. (15). However, it contains an additional term, proportional to the parameter µ δ , due to the non-degeneracy of the heavy quarks and reads where χ T h (x) = (s, c). Moreover, µ l and µ σ are the twisted mass parameter. The mass terms m l and m h are the (untwisted) Wilson quark masses tuned to the critical value m crit (i.e. at maximal twist), which ensures automatic O(a) improvement [60] for parity even operators. The equivalent discussion for non-local operators can be found in Refs. [38,[61][62][63]. Fields in the "physical basis" can be obtained from the twisted basis through the transformation with α = π/2 at maximal twist. From now on we will use fields in the physical basis.
The ensemble we use has lattice volume V = 32 3 × 64, with a lattice spacing of a = 0.0938 fm. The pion mass is approximately equal to m π = 260 MeV and m π L ≈ 3. In Table II, one can find the summary of the main parameters characterizing the ensemble. For further details see Ref. [58].  II: Parameters of the ensemble used in this work. The nucleon mass (mN ), the pion mass (mπ) and the lattice spacing (a) are determined in Ref. [64].
A. Numerical methods

Connected diagrams
To improve the overlap between the states generated by the interpolating field N α (x) = ε abc u a α (x) d bT (x)Cγ 5 u c (x) and the proton ground state we employ Gaussian smearing [65,66]. In addition, we use APE smearing for the gauge links that enter the Gaussian smearing. The optimal parameters for the Gaussian and APE smearing techniques, determined in Ref. [67], are (α G , N G ) = (4.0, 50) and (α APE , N APE ) = (0.5, 50), respectively. Moreover, to further improve the overlap with the boosted proton ground state, we use momentum smearing [68], as it has been proven to drastically reduce the statistical noise in the matrix elements of boosted hadrons [69]. In Ref. [70] and in this work, the momentum smearing parameter has been tuned to ξ = 0.6, which minimizes the statistical errors. The momentum smearing operator S on a quark field ψ(x) reads where ξ is the momentum smearing parameter and j runs over the spatial directions, with U j (x) being the link in the spatial j-direction.
To evaluate the connected contributions to the three-point functions, we employ the sequential method [71] through the sink. Moreover, to further increase the number of measurements, we compute the three-point functions with N src different source positions on each configuration and we boost the nucleon along all the spatial directions and orientations, i.e. ±x, ±y, ±z. Indeed, in Ref. [41] it was found that the statistical uncertainty decreases as 1/ √ N src N dirs for all the operators under consideration, with N dirs = 6 being the number of directions of the nucleon boost. The number of source positions employed depends on the nucleon boost, and is N src = 8 for the two lowest values of the momentum and N src = 14 at P 3 = 1.24 GeV for the third. The source-sink separation is t s = 0.94 fm for the lowest momentum value, and t s = 1.13 fm for the two highest ones. The value employed for t s at P 3 = 1.24 GeV is expected to be large enough to suppress excited-states contamination [41]. In Table III we report the number of measurements for the connected contributions at each value of P 3 .

Disconnected contribution
The evaluation of the disconnected quark loops with a Wilson line in the boosted frame constitutes the most computationally demanding aspect of this work. The isoscalar three-point function of Eq. (2) (τ = 1 and ψ( The three-point function C 3pt contains a connected and disconnected part. The latter is given by where x = (t s , x), y = (t ins , y) and z = (0, 0, 0, z). The quantity G f ( x, t x ; y, t y ) is the all-to-all propagator with quark flavor f = u, d, s, from each lattice point x to any point y. Eq. (19) is a correlation of two parts: the nucleon two-point function and the quark loop with a Wilson line. The latter can be written as where the trace in the second line is intended over volume, spin and flavor indices. Due to the presence of the all-to-all propagator G f ( x, t x ; y, t y ), the exact evaluation of the disconnected contribution in Eq. (21) would require ≈ 10 7 inversions of the Dirac operator per configuration for the lattice that we are considering. In contrast, stochastic techniques allow to carry out the computation of the disconnected loops with a feasible but yet high computational cost compared to the currently available resources. Stochastic techniques employ noise sources ξ r (x) that obey two properties where the product between source vector has to be intended as a tensor product in volume, spin and color subspaces. Given the set of vectors ξ r (x), the all-to-all propagator can be constructed by solving the equation M (x, y) being the Dirac twisted-mass operator. Having the set of solutions φ r (x), the all-to-all propagator can be estimated via The number of stochastic vectors required so that the stochastic error becomes comparable to the gauge error should be much smaller than calculating the all-to-all propagator exactly. In addition, exploiting a property of the twisted mass operator, it is possible to design a stochastic algorithm that further reduce the computational cost. Recalling the transformation of Eq. (17), the insertion operator Γ in twisted basis reads with α = π/2 at maximal-twist. Depending on the operator Γ it is possible to exploit two properties of the twistedmass operator to evaluate the loop of Eq.
In this case, we apply the standard one-end trick, that exploit the following property of the twisted mass operator The transversity operator Γ = σ 3j belongs in this category.
2. If {Γ tm , γ 5 } = 0, then the loop in twisted basis possess the same analytical form as in the physical basis The quantity of interest can be computed with the generalized one-end trick, exploiting the following property where D W is the massless Wilson-Dirac clover operator. Note that if the twisted mass parameter becomes very small (close to the physical point) this type of one-end trick is approaching the standard definition in Eq. (24). The helicity operator Γ = γ 5 γ 3 and the unpolarized Γ = γ 0 belong to this category.
One of the technical aspects of the calculation is the evaluation of the traces in Eqs. (26) - (28). With small quark masses, the contribution to the loops coming from the low modes of the spectrum of the Dirac operator may be sizeable, and contributes significantly to the stochastic noise [74]. Therefore, we compute the first N ev = 200 eigenpairs λ j , |v j of the squared Dirac operator M u M † u , that allow to reconstruct exactly the low-mode contribution to the disconnected quark loops. At this stage, stochastic techniques can be employed with the deflated operator to evaluate the high-modes contribution to the traces. To reduce the stochastic noise, we use the hierarchical probing algorithm [25], that allows to reduce the contamination to the trace coming from off-diagonal terms up to a distance 2 k . This improvement is achieved by partitioning the lattice with 2 d(k−1)+1 Hadamard vectors, where d = 4 is the number of dimensions of the lattice. Finally, to remove the contamination from off-diagonal terms in spin-color subspaces, we apply full dilution [75]. The algorithm employed in the present work has been successfully used in other studies involving the evaluation of disconnected contributions [17][18][19][20].
For each value of the proton boost P 3 = 0.41, 0.83, 1.24 GeV, we evaluated the two-point functions contributing to the disconnected diagram of Eq. (19) with N srcs = 200 source positions (see Sec. III B). In addition, apart from averaging over all possible directions and orientations of the nucleon boost, we also average over forward and backward projections. In Table IV, we report the total statistics collected for the disconnected three point correlators. We note that, for the second largest value of the momentum used, namely P 3 = 1.24 GeV we use ≈ 10 6 measurements. In addition, we also compute the matrix elements for P 3 = 1.65 GeV and all Γ using approximately the same statistics as the previous smaller boost. While this number of statistics is not sufficient to obtain the same statistical accuracy as lower momenta, it allows us to check whether convergence with P 3 is reached.

Loops
Two-point functions   The two-point functions enter the calculation through the ratio of Eq. (4), but also contribute to the evaluation of the disconnected diagram, as shown in Eq. (19). For this reason, to obtain a significant amount of measurements for the disconnected contributions we compute the two-point functions with a large number of source positions, N 2pt src = 200, and consider boosts of the nucleon along the different spatial directions and orientations. This procedure allowed us to considerably reduce the statistical error in the disconnected contributions at small computational cost because the same loops is combined with all 200 two-point functions on the same configurations. Given that the computational cost of the two-points function is considerably lower compared to the one required to evaluate the disconnected quark loops, using multiple source positions is highly beneficial. In Table IV we report the number of measurements of the two-point function performed at each nucleon boost P 3 .
The two-point function can be written as with |n being the n th energy state of the interpolator N α (x) and E n (P ) its energy. We performed the analysis by keeping up to two terms in the expansion of Eq. (30). In particular, the two-state fit function of the two-point correlator consists of while the effective energy reads with ∆E = (E 1 − E 0 ) and B = c 1 /c 0 . In Fig. 2, we show the results of the two-state fits of the correlator and the effective energy for the parameters E 0 , ∆E, c 0 and c 1 /c 0 , varying the low-end of the fit interval t min . The results show that the fits on the correlator or the effective energy lead to the same ground state energy. Furthermore, the plateau and two-state fits converge at t s /a = 9 for momentum 1.24 GeV.
In Table V we report the parameters extracted using one-and two-state fits. The results for E 0 are obtained with the plateau fit of the effective energy, and they are compatible with the values extracted using two-state fit results. In Fig. 3 we show the effective energy for the second largest momentum P 3 = 1.24 GeV, together with the plateau fit   and two-state fit results. By iterating the fit procedure described above over the data for the different nucleon boosts, we reconstructed the dispersion relation with m N being the nucleon mass. In Fig. 3 we show the observed trend of the energy with the nucleon boost P 3 together with a linear fit performed with the function of Eq. (33), giving a 2 m 2 N c 4 = 0.2678(8) and c 2 = 1.003 (7). As can be seen, the lattice data are fully compatible with the dispersion relation for all values of P 3 .

IV. DISCONNECTED MATRIX ELEMENTS
Obtaining the disconnected contributions to the up-, down-and strange-quark PDFs is the central goal, and most laborious aspect of this work. We use the techniques outlined in Section III A 2 to extract the matrix elements and study systematic uncertainties, such as excited-states contamination.

A. Excited-states contamination
To extract reliably the ground-state contribution to the matrix elements, we evaluate the ratio between the threeand two-point functions at seven source-sink separations, ranging from t s = 0.563 fm to t s = 1.126 fm in steps of a = 0.0938 fm. For disconnected contributions, the evaluation of different source-sink separations does not require new inversions. This allowed us to study the excited-states contamination to the matrix elements using several t s values and three analysis methods: plateau fit, two-state fit and summation method. We briefly summarize these methods.
1. Two-state fit. In Sec. III B, the two-point correlator is expanded up to the first excited state. Likewise, we can expand the three-point correlator keeping terms up to the first excited state. This gives four terms, that is Being interested in the forward kinematic limit allows one to reduce the number of independent parameters, since A 0,1 = A 1,0 . We performed a fit of the ratio of Eq. (4) with the function where the parameters c 1 /c 0 and ∆E are determined through the effective energy fit and the results are reported in Table V. Thus, the parameters determined by fitting the ratio of three-and two-point functions are A 0,0 /c 0 , A 0,1 /A 0,0 and A 1,1 /A 0,0 . A 0,0 /c 0 corresponds to the matrix element we are interested in. Such fits are weighted by the statistical errors, and therefore the fit is driven by the most accurate data. Since the statistics for the disconnected contributions are independent of th t s value, we repeat the two-state fits modifying, each time, the starting value of t s entering the fit (t low s ). The results from the two-state fit method, allows us to verify ground-state dominance by comparing the matrix elements from the individual t s values.
2. Plateau fit. For 0 τ t s and ∆E t s 0, the first term in the ratio of Eq. (35) dominates. Thus, the matrix elements can be extracted by performing a constant fit on the ratio of Eq. (4) in the region defined 0 τ t s , with large enough source-sink separation. We exclude from the fit range three points from left and right, i.e. we evaluate the weighted average in the interval τ ∈ [3, t s − 3]. While the excited-states contamination decreases with t s , at the same time the statistical uncertainty exponentially increases. For this reason, the determination of the ground state of the matrix elements with the plateau fit method is a challenging task, and the results need to be compared with other analysis techniques.
3. Summation method. Summing over the insertion time τ of the ratio of the three-and two-point functions we find [76,77] S(t s ) = Thus, the matrix elements corresponds with the slope of the straight line S(t s ), and can be measured by performing a linear regression.
Using the three aforementioned approaches we analyze the excited-states effects on the matrix elements for the unpolarized, helicity and transversity PDFs. In the next three subsections, we present the analysis of the real and imaginary parts of the matrix elements at P 3 = 1.24 GeV, as a representative example.

Unpolarized
We start by discussing the analysis of the unpolarized isoscalar u + d disconnected matrix elements. In Fig. 4 we show the ratio of three-and two-point functions at P 3 = 1.24 GeV for the unpolarized operator for z/a = 3 and we compare the results obtained with the three analysis methods reported in Sec. IV A. The real part of the matrix elements shows no substantial dependence on the source-sink separation, and the plateau fit results obtained at different t s give all compatible results. The dependence of the two-state fit on the lowest source-sink separation t low s , included in the fit shows a constant trend, which is also compatible with the results obtained with the summation method. The reduced chi-square χ 2 /d.o.f = 0.96 suggests that the function of Eq. (35) provides a good description of the data. To extract the matrix elements, we compute the constant correlated fit of the plateau fit results starting from t low s /a = 9. In contrast to the real part, the imaginary part shows a large effect due to the excited-states contamination. However, the two-state fit is compatible with the plateau value using t s /a = 11. Also, the results obtained at different t low s using the summation method are compatible with the other methods within uncertainties. As final results for the unpolarized matrix element, we report the ones from the plateau fit for t s /a = 11, which is compatible with the results obtained with the two-state fit and summation method.

Helicity
The disconnected contributions to the helicity isoscalar matrix elements is purely real and exhibits a non-negligible dependence on the source-sink separation. We observe a decreasing behavior as t s increases, which is a behavior also observed in the axial charge [18]. In addition, the plateau fit for t s /a > 10 are compatible with the two-state fit. Therefore, we use the plateau fit for t s /a = 10 as our final results, so that statistical uncertainties are controlled.

Transversity
The ratio of the three-and two-point functions for the disconnected contributions to the isoscalar transversity matrix elements does not shows dependence on the source-sink separation for both the real and imaginary parts (see Fig. 6). Thus, the matrix elements are computed from the plateau fit for t s /a = 9, both for the real and the imaginary parts. Using the criterion adopted for selecting the final results for each PDF case we compared the extracted matrix elements in Fig. 7 for all values of z. In summary, the plateau fits are evaluated at t s /a = 9, 10, 9 (t s /a = 11, 10,9) for the real (imaginary) part of the unpolarized, helicity and transversity, respectively. For the summation method, we employ all t s values available, except for the imaginary part of the unpolarized and the real part of the helicity matrix elements, where t s /a = 6 is excluded as explained above. The two-state fit is performed in the range t s /a ∈ [6,12] in all cases, except for the imaginary part of the unpolarized and transversity matrix elements, where the t s /a = 6 value, as explained, is not included in the regression.
Our conclusions for the isoscalar disconnected matrix elements apply also to the strange matrix elements, with the excited-states contamination showing similar effects.

B. Momentum dependence
As explained in Sec. II C, the matrix elements and the quasi-PDFs have a dependence on the nucleon boost P 3 , which also enters the matching formula leading to the PDFs. An important aspect of the study is the investigation of the momentum dependence of the matrix elements, which affects the convergence to the light-cone PDFs. In Fig. 8, we present the results for the renormalized strange and isoscalar disconnected matrix elements as a function of the momentum boost. For the unpolarized case, the real part decreases in magnitude as the P 3 increases, and becomes compatible with zero. In contrast, its imaginary part is non-zero and shows convergence for the two largest values of P 3 . We find that the isoscalar disconnected matrix elements share the same qualitative behavior as the strange-quark ones.
First results on the helicity distribution appeared in Ref. [26]. Here we show results with increased statistics, and with the addition of P 3 = 1.65 GeV. The matrix elements show a mild residual dependence on the momentum. The imaginary part of the renormalized matrix elements arises entirely from the complex multiplication with the renormalization function and the bare matrix elements (see Eq. (10)). Indeed, as mentioned already in Sec. IV A 2, the disconnected contribution to the bare matrix element for the helicity distribution is purely real.
The real part of the matrix elements for the transversity distribution exhibit a strong dependence on the nucleon boost, changing dramatically as we increase P 3 from 0.83 GeV to 1.24 GeV. However, results obtained for P 3 = 1.65 GeV show agreement with those for P 3 = 1.24 GeV albeit the large uncertainty. From the current results it is still unclear if convergence is reached. However, in order to fully check this would require a larger momentum and much more measurements to reach the required accuracy. This is beyond the current study and will be tested in a followup work. We thus, construct the PDFs using the results for P 3 = 1.24 GeV. In Sec. VII D 2, we comment on the region of x affected by the gap observed in the real part of the matrix elements as momentum increases. In contrast, the imaginary part is fully compatible with zero for the two lowest momenta, while it is slightly non-zero at large z at the highest momentum.

V. CONNECTED MATRIX ELEMENTS
The evaluation of the connected matrix elements contributing to the three types of PDFs has been studied in our previous works. In particular, we refer the reader to the study of Ref. [41], where several sources of systematic uncertainties were discussed in great detail. For completeness, we briefly discuss here the connected contributions which are needed for the flavor decomposition. In Fig. 9, we show the momentum dependence of the bare connected contributions to the isoscalar and isovector matrix elements for the three types of PDFs. In all cases, as the nucleon boost increases, the real part of the matrix elements decay to zero faster, and the magnitude for the imaginary part increases in the region z/a < ∼ 9. The unpolarized matrix elements show convergence with P 3 while the imaginary parts of the helicity and transversity distributions increase in magnitude. We note that in order to compute the connected contributions for a fourth larger boost would require new inversions and large number of measurements to reduce the errors sufficiently enough to check convergence. Thus for the current work we opt to use the results for P 3 = 1.24 GeV for the connected parts since the focus of this work is the evaluation of the disconnected contributions. The uncertainty on the unpolarized distribution is smaller as compared to the other two distributions. This behavior is due to the fact that both the three-and two-point functions share the same projector, (1 + γ 0 )/2, which increases the correlation between the two quantities and, as a result, drastically decreases the noise-to-signal ratio. We note also that the transversity distribution reported here is the average over the two insertions σ 3j with j = 1, 2. show respectively the isoscalar connected contribution and the isovector unpolarized matrix elements. The left column shows the real part and the right the imaginary part. The same flavor combinations are reported respectively for the helicity and transversity distributions in the 3 rd and 4 th rows, and in the last two rows. We show the matrix elements at P 3 = 0.41 GeV (blue), 0.83 GeV (green) and 1.24 GeV (red). Data points are slightly shifted to to improve readability.

VI. NUCLEON CHARGES
The nucleon charges are usually extracted from the nucleon matrix elements of local operators. This limit is obtained from the matrix elements of non-local operators at z = 0. Since the charges are frame independent, any value of P 3 may be used. Indeed, in Sec. V we demonstrate that the z = 0 have little dependence on P 3 For the disconnected contributions, we have the matrix elements at P 3 = 0 (rest frame). For the connected contributions to the charges, we use the lowest momentum, so we control statistical uncertainties. In what follows, we will show the results obtained for the isovector u − d, isoscalar u + d and strange-quark vector, axial and tensor charges, g V , g A and g T . First, we describe our results for the disconnected contributions, whose total number of measurements in the rest frame is N meas = 66 · 10 3 . In Fig. 10 we show our results for the renormalized disconnected vector, axial and tensor charges. The integral over the volume (i.e. Fourier transform at zero momentum transfer) of the trace of the vector currentψ(x)γ 0 ψ(x) is zero because quark and antiquark loops contributions cancel each other. Thus, the disconnected contribution to the unpolarized isoscalar and strange matrix elements in the absence of the Wilson line and at P 3 = 0 expected to be zero is verified and this constitutes a consistency check of our computations. Due to excited-states contamination, the disconnected isoscalar and strange vector charges g u+d (disc) V and g s V obtained with the two-state and plateau fits are not compatible with zero at small source-sink separation. In particular, the two-state fit results become compatible with zero when t low ≥ 10. The results obtained with the summation method have the largest uncertainties and are compatible with zero. The axial charge shows the larger contamination from excited states. In particular, the plateau fit results show a decreasing trend with the t low , converging to a constant value for t low /a ≥ 11 for the isoscalar and t low /a ≥ 8 for the strange charges, which are selected as our final values. In contrast, both the results obtained with two-state fit and summation are constant and compatible with the selected plateau fit results. We find g u+d (disc) A = −0.104(10), g s A = −0.0320 (28) .
The results on the disconnected contributions of the tensor charge show very mild excited states effects. We use the plateau value extracted by fitting the ratio to a constant for t s /a ≥ 9 for the isoscalar connected tensor charge g u+d(disc) T and for t s /a ≥ 8 for g s T . Our final results for these two quantities are g u+d (disc) T = −0.00818(91), g s T = −0.00265 (60) .
We stress that despite the agreement of the strange tensor charge with the value extracted using local operators [12,78], a direct comparison is not meaningful since we are using gauge ensembles simulated with heavier than physical pion mass. It thus comes with no surprise that the disconnected isoscalar tensor charge differs from the value obtained at the physical pion mass.
The connected contributions to the nuclear charges are computed for the smallest momentum P 3 = 0.41 GeV. Using these results we can extract the values for each quark flavor for the vector, axial an tensor charges. The details on the computation of the connected isoscalar and isovector contributions are given in Sec. V. The connected contributions used to extract the charges are obtained from plateau fits with t s /a = 12. The nucleon axial and tensor charges are given in Table VI. We note that for the vector charge g V we find results that are consistent with charge conservation.  VI: Results for the isovector (first column), isoscalar connected (second column) and disconnected (third column) and for the up (fourth column), down (fifth column) and strange (sixth column). We show our results on the axial (second row) and tensor (third row) nucleon charges.

VII. PARTON DISTRIBUTION FUNCTIONS
A. Isoscalar and isovector renormalized matrix elements In Fig. 12 we show the momentum dependence of the total renormalized isoscalar and isovector matrix elements, including disconnected contributions. The renormalized matrix elements are reported as a function of zP 3 . We renormalize and apply the matching procedure independently for the isoscalar and isovector distributions, allowing us to obtain the individual up and down quark PDFs.
The source-sink separation used is t s = 0.94 fm for the lowest momentum and t s = 1.13 fm for P 3 = 0.83 and 1.24 GeV. From previous studies of the isovector distributions (see, e.g. Ref. [41]) and nucleon charges [18], we expect that excited-states contamination is more significant for the nucleon three-point correlators of the axial and tensor currents as compared to the vector. However, the statistical uncertainty is larger for these quantities and within our current errors the source-sink separation employed at the highest momentum is sufficient to suppress excited states to this level of accuracy [41].
A summary of the results for the connected matrix elements in the absence of the Wilson line are reported in Table VII. The momentum dependence of all the matrix elements analyzed is negligible for z = 0 as expected For example, the isoscalar connected matrix elements for the unpolarized distribution, h u+d (z = 0), for the largest momentum differs from the others by less then 1%. The isovector unpolarized matrix elements at z = 0 is independent of the momentum boost, and equal to 1, as expected from charge conservation. Regarding the isoscalar helicity case, we still find agreement for different P 3 within uncertainties, but with larger fluctuations of the mean values, as the disconnected contribution is about ∼ 17% of the connected part. We note that the ∆h u−d (z = 0) is compatible with the experimental value g u−d A = 1.27641(56) [79]. Insensitivity to the momentum boost is also observed in the transversity case.

B. Truncation of the Fourier transform
In order to construct the x-dependence of PDFs we need to take the Fourier transform. Since the matrix elements are determined for discrete finite number of z values, we study the dependence on the cutoff z max to understand systematic effects related to the reconstruction. In particular, for all types of distributions we verify that a z max exists such that, addition of information for z > z max in the Fourier transform leaves the PDF unchanged within statistical uncertainties. This value of z is selected as the maximum value z max included in the Fourier transform and, typically, the matrix elements at this value has a vanishing real part. Note that the latter is just a qualitative criterion and in practice we always check by increasing z for convergence. In Fig. 11 we show the dependence on the cutoff z max for the isoscalar and isovector distributions at P 3 = 1.24 GeV. We compare the results obtained with the discrete Fourier transform of Eq. (12) with the results from the Bayes-Gauss-Fourier transform (BGFT) [80]. The latter is an advanced reconstruction technique based on Gaussian process regression, which allows to obtain an improved estimate of quasi-PDF for continuous values of x, starting from a discrete set of data obtained with lattice QCD computations. The chosen values of z max for each quark flavor and each operator are given in Table. VIII. From Fig. 11 it is clear that increase of the cutoff beyond the reported values of

C. Isoscalar and isovector distributions
The isoscalar and isovector PDFs are extracted from the corresponding renormalized matrix elements shown in Fig. 12. For the isoscalar combination, we add both the connected and disconnected contributions. We plot the matrix elements against zP 3 , which is the argument of the exponential in the Fourier transform.
In all cases, we find that the matrix elements for the lowest momentum P 3 = 0.41 GeV do not decayed to zero for large z, demonstrating, as expected, that the momentum is not large enough. By increasing the momentum to P 3 = 1.24 GeV, the matrix elements become consistent with zero within their uncertainties. While the imaginary parts show a residual momentum dependence, the convergence must be checked at the level of the reconstructed PDF. This is due to the fact that P 3 enters the matching kernel and affects the convergence. Therefore, to address the momentum convergence as we increase P 3 , we show in Fig. 13 the momentum dependence of the isoscalar and isovector PDFs. We use the standard Fourier transform, with the values of z max given in Table VIII, as discussed in Sec. VII B. As can be seen in Fig. 13, the overall dependence on the two largest values of the momentum is relatively small. Dependence on P 3 is observed in the unpolarized isoscalar PDF. In general, the PDFs for the smallest momentum, do not show convergence, and exhibit non-physical oscillations due to the presence of systematic effects in the reconstruction of the x-dependence. However, such oscillations are suppressed for the higher values of P 3 . The isoscalar and isovector helicity distributions have a similar magnitude and exhibit milder dependence on the boost as compared to the unpolarized. In particular, both isoscalar and isovector helicity distributions are consistent for P 3 = 0.83 GeV and P 3 = 1.24 GeV. Finally, the isoscalar and isovector transversity distributions also show nice convergence with P 3 for the two largest values. These distributions will be used for the flavor decomposition presented in Sec. VII D together with comparison of our data with phenomenology.

Light quark distributions
Our results on the isoscalar and isovector distributions presented in Sec. VII C allow us to extract the up and down quark contributions for the unpolarized, helicity and transversity distributions. The disconnected contributions are taken into account in all cases. We stress that the comparison with phenomenology can only be qualitative for a number of reasons: i) We use an ensemble with larger than physical pion mass. We know from previous studies that there is a non-negligible pion mass dependence on the PDFs; ii) lattice systematics, such as cut-off effects, are not taken into account; iii) the renormalization ignores mixing present in the case of the unpolarized and helicity singlet PDFs; and iv) errors are still sizable and may hide systematics, such as convergence with the boost. However, it is still interesting to compare with phenomenology keeping these caveats in mind. The results for the unpolarized PDF at the largest momentum are compared with data by NNPDF3.1 [81], while the helicity distribution is compared with JAM17 [82] and NNPDF POL1.1 [2]. Finally, the quark transversity distribution obtained in this study is compared against the SIDIS data [83] and SIDIS data constrained by the value of tensor charge g T computed in lattice QCD [83]. For the anti-quark region for the NNPDF3.1 data, we include the crossing relations of Eq. (14), such that we show the antiquark distributions in the negative-x region.
The light-quark contributions to the unpolarized PDF show good agreement with phenomenology in the region x > ∼ 0.2. Also, the region x < ∼ −0.2 both estimates are compatible with zero. Note that lattice results for the small-x region (|x| < ∼ 0.15) suffer from uncontrolled uncertainties due to the reconstruction of the PDFs and the values of the lattice spacing used. The case of the helicity distributions is very interesting, as it has non-negligible contribution from the disconnected diagram. Our results for the up quark helicity show similar features as the NNPDF data, but are have higher values. The down quark distribution gives compatible results both with NNPDF POL1.1 and JAM17 data for all x in the physical region [−1, 1]. The transversity distribution is the least known collinear PDF and it (bottom panels) distributions at P 3 = 1.24 GeV (red band). We also show the NNPDF results [2,81,84] (blue band) and JAM17 [82] (orange band) phenomenological results. For the transversity PDF we compare against the SIDIS data [83] (green band) and SIDIS data constrained by the value of tensor charge g T computed in lattice QCD [83] (gray band).
is not well-constrained by SIDIS data. As a result, global fits for the light quark δq(x) carry large relative error of ≈ 50−100% [83]. A more precise phenomenological estimate of the transversity PDFs can be obtained by constraining the distributions with the value of the tensor charge g T computed within lattice QCD [83]. A comparison with the latter, reveals a similar agreement as for the helicity PDFs. We would like to stress that the overall qualitative agreement is very promising, as this computation is done using simulations with heavier than physical pions.

Strange quark distributions
The strange distributions presented here are computed using the renormalized matrix elements shown in Fig. 8. The values of z max employed in the Fourier transform defining the quasi-PDF are reported in Table VIII. The criterion adopted to select z max is to analyze the dependence of the PDF as z max is increased, as discussed in the previous section. In Fig. 15 we show the unpolarized, helicity and transversity PDFs. The antiquark distribution reported here takes into account the crossing relations in Eq. (14), showing the anti-quark distributions in the negative x region. Although the unpolarized PDFs extracted from the matrix element using the two largest momenta tend towards the phenomenological result, there is still some residual dependence, which points to the need to increase the momentum boost to check the independence on P 3 . Due to the simultaneous suppression of the real part of the matrix elements and the enhancement of the imaginary part,s(x) becomes symmetrical with respect to x = 0 as the momentum boost increases. This symmetry feature is exploited in the global fits. The results for the helicity distribution are approximately symmetric in the quark and antiquark regions, and are compatible with the results from the NNPDF POL1.1 [2] and with JAM17 global fits analysis both of which have larger uncertainties. Our results, thus, provide valuable input for phenomenological studies. In fact, this is more evident for the strange transversity distribution where experimental results are lacking. We obtained results on the transversity PDF with  Besides the individual s(x) ands(x) distributions, there is also an interest on the strange-quark asymmetry. This is partly due to the fact that there is no symmetry to suggest that the two distributions have to be the same.
The strange and anti-strange asymmetry has been discussed within chiral effective theory [85,86], perturbative evolution of QCD [87], and a physical model for parton momenta [88]. Here, we study the asymmetry using our data for P 3 = 0.41, 0.83, 1.24 GeV, and the results are shown in Fig. 16. In contrast to the individual s(x) ands(x) distributions, here we find that there is no momentum dependence in the strange-quark asymmetry. We also note that the difference between s(x) ands(x) is a non-singlet combination and, thus, does not mix with the gluon PDFs. Focusing on the most accurate results at P 3 , we find that the asymmetry vanishes at x > ∼ 0.2 and is small but non negligible in the small-x region. This conclusion is, at present stage, qualitative, and an investigation of systematic effects is needed before drawing quantitative conclusions.

E. Moments of nucleon PDFs
In this section, we calculate the moments x n of the three PDFs considering n = 0, ..., 3. The n−th moment of the unpolarized, helicity and transversity distributions are defined as where we employed the crossing relations of Eq. (14), to write the moments as a function of the quark distributions only. In Table IX we report the results for the isovector, isoscalar and flavor diagonal moments. The zero-th moments are compatible with the nucleon charges reported in Tab VI. This is a non-trivial check, as the calculation of the charges follows a totally different procedure and undergoes a Fourier transform and matching.

VIII. CONCLUSIONS
In this work we present a study of the x-dependence of proton collinear quark PDFs from lattice QCD considering both connected and disconnected diagrams. These contributions are necessary to determine the individual-flavor contributions to PDFs. We present results for the up, down and strange quark for the unpolarized, helicity and transversity PDFs. This work extends our first calculation for the flavor decomposition of the helicity PDFs [26]; here we increase the statistics by about a factor of two, and include results on the unpolarized and transversity PDFs.
The main goal of this work is to explore the feasibility of the calculation of disconnected quark loops with non-local operators. To this end, we provide the necessary details on technical and theoretical aspects, as well as the examination of some sources of systematic uncertainties. The calculation is carried out using one ensemble of N f = 2+1+1 twisted mass fermions simulated with quark mass value that produces a pion mass of 260 MeV. Using a single ensemble, we can address excited-states contamination, reconstruction of the x dependence, and the convergence with increasing the momentum boost in the final PDFs.
The matrix elements contain non-local operators with the length of the Wilson line extending up to half the spatial extend of the lattice. The proton states are boosted with momentum using three values, namely P 3 = 0.41, 0.83, 1.24 GeV. Several values of the source-sink time separation are considered. As we increase the boost we also increase the source-sink time separation in order to investigate the effect of excited states; we employ up to t s = 1.13 fm for the highest momentum (see Figs. 4 -7). Both the isovector and isoscalar flavor combinations are calculated, with the latter receiving contributions from the connected and disconnected diagrams. All matrix elements are renormalized multiplicatively using the RI scheme and evolved to the modified-MS scheme at a scale of 2 GeV (unpolarized and helicity) or √ 2 GeV (transversity). The renormalization is followed by the transform to the momentum space, x, which produces the quasi-PDFs. We apply the standard Fourier transform that is confirmed using the Bayes-Gauss-Fourier transform, finding compatible results (see Fig. 11). The quasi-PDFs are matched to the light-cone PDFs using one-loop perturbation theory. The matching kernel contains information on the renormalization scheme and scale for the quasi-PDFs (modified-MS scheme at 2 or √ 2 GeV), which are then matched to the light-cone PDFs in the MS scheme at the same scale. In this proof-of-principle study we neglect the mixing with the gluon PDFs for the unpolarized and helicity case. The extraction of the latter has its own challenges and will be included in our future studies. To test the influence of systematic uncertainties, we calculate the charges by integrating the PDFs (Table VI), and compare with the values obtained directly from the matrix elements (Table IX). We find good consistence among the results.
We find that the light-quark disconnected contributions have the most impact for the helicity PDF, while the transversity disconnected contribution is very small. Regardless, a clear non-zero signal is found in all cases. The strange-quark PDFs are nonzero up to x ∼ 0.5, with the unpolarized and helicity having a similar magnitude, and the transversity being an order of magnitude smaller, as can be seen in Fig. 15. These distributions are very challenging to extract from experimental data due to the lack of sensitivity to the strange-quark. In a qualitative comparison of our results with phenomelogically extracted PDFs we find: i) our results on the unpolarized have a statistical precision which is similar to the NNPDF data; ii) the helicity strange-quark PDF is significantly more accurate than the JAM and NNPDF results and iii) our results for the strange-quark transversity PDF serve as a prediction.
There are a number of improvements that one can do to quantify and eliminate systematic uncertainties. These include, but not limited to, pion mass dependence, mixing under matching with the gluon PDFs for the unpolarized and helicity case, and address the inverse problem using, e.g., Bayesian reconstruction methods. One can also address finite-volume and discretization effects, which require extracting the PDFs using multiple ensembles. However, this work clearly demonstrates the great potential in the extraction of the x-dependence of individual quark PDFs from lattice QCD.