Flavor decomposition for the proton helicity parton distribution functions

We present, for the first time, an \textit{ab initio} calculation of the individual up, down and strange quark helicity parton distribution functions for the proton. The calculation is performed within the twisted mass clover-improved fermion formulation of lattice QCD using one ensemble of dynamical up, down, strange and charm quarks with a pion mass of 260 MeV. The lattice matrix elements are non-perturbatively renormalized and the final results are presented in the $\overline{ \rm MS}$ scheme at a scale of 2 GeV. We give results on the $\Delta u^+(x)$ and $\Delta d^+(x)$, including disconnected quark loop contributions, as well as on the $\Delta s^+(x)$. For the latter we achieve unprecedented precision compared to the phenomenological estimates.

Introduction. The theory of the strong interaction, Quantum Chromodynamics (QCD), successfully explains the structure of hadrons and their interactions. The fundamental degrees of freedom in QCD are the quarks and gluons, collectively referred to as partons. Partons are responsible for the rich internal structure of hadrons. Most of the knowledge of the complex hadronic structure comes from parton distribution functions (PDFs), a set of number densities describing the non-perturbative QCD dynamics. Distribution functions are universal quantities and, therefore, can be accessed by a variety of high-energy scattering processes. The cross section of such processes can be factorized into a perturbative component calculable in perturbative QCD, and a non-perturbative part expressed in terms of the partonic densities. The generalized parton distributions (GPDs) and transverse momentum distribution functions (TMDs) complement PDFs, and are necessary for the 3dimensional mapping of the hadrons.
At leading order within the parton model, the PDFs have a simple interpretation. The unpolarized PDFs are interpreted as the probability to find an unpolarized parton with a longitudinal momentum fraction x within an unpolarized nucleon. The helicity PDFs can be interpreted as the difference between finding quarks with spins aligned and opposite to that of a longitudinally polarized nucleon. The colinear PDFs are completed with the transversity PDFs, which have the interpretation of finding quarks polarized in the same or opposite direction as a transversely polarized nucleon. PDFs play a central role in the on-going experimental program of major facilities, such as, BNL, CERN, DESY, Fermilab, JLab and SLAC (see, e.g., Refs. [1][2][3]). These experiments provide a wealth of measurements that are jointly analyzed within the framework of phenomenological fits. Based on the available experimental data, the most well-studied colinear distributions are the unpolarized, followed by the helicity with an order of magnitude less experimental data sets, namely a few hundred data sets [4,5]. The transversity PDFs are even less-known [6]. The accessible kinematical region is more limited for the helicity and transversity PDFs as compared to the unpolarized PDFs, and therefore, the reconstruction of PDFs uses input from models. Such input introduces dependence on the functional forms employed. As a consequence, the extraction of the helicity and transversity PDFs are, to some extent, driven by the fit functions, due to the lack of experimental data (see, e.g., Ref. [6] for a discussion). The dependence on the analysis procedure is evidence by the tension among some of the global analyses [5,[7][8][9].
The focus of this work are the helicity PDFs, which are typically accessed experimentally in deep-inelastic scattering (DIS), semi-inclusive DIS, Drell-Yan, and protonproton scattering processes. Currently, the global analyses use next-to-leading order (NLO) corrections in perturbative QCD (NNPDF POL 1.1, DSSV14, JAM17) [5,[7][8][9]. In these analyses, the up and down contributions, ∆u(x), ∆d(x) are better constrained in the valence sector, with ∆u(x) being more precise. On the other hand, constraining ∆s(x) is not successful, as the kinematic regions of some of the data sets (e.g., the W -boson production data) are not sensitive to the strangeness [5]. The situation somewhat improves with the inclusion of kaon production SIDIS data, but it is still unsatisfactory, and influenced by theoretical assumptions. In the recent work of the JAM Collaboration [7] the authors used inclusive and semi-inclusive data, and find, for both sets of data, the strange polarization to be very small and consistent with zero. More details on the global analyses can be found in the recent reports [6,10].
Based on the current status of phenomenological analyses, an extraction of the PDFs from first principles is highly desirable. Here we present the first extraction of the up, down and strange helicity PDFs for the proton using lattice QCD, the only known ab initio formulation of QCD. We study both the valence and sea quark contributions that allow one to perform a controlled decomposition of the various distributions. To obtain the x-dependence of PDFs, we implement the quasi-PDF method [11]. This approach is based on correlation functions that are calculable on a Euclidean lattice. The matrix elements are between include hadron state with a finite momentum P = (0, 0, P 3 ). A non-local operator with fermion fields separated by a distance z connected by a Wilson line, is inserted between the proton states. Note that the Wilson line is in the same spatial direction as P . Naturally, the matrix elements are defined in coordinate space, with z varying from zero up to half the spatial extent of the lattice. To extract physical quantities, a Fourier transform is applied on the matrix elements to obtain the so-called quasi-PDFs, which are defined in momentum space, x. For large values of P 3 , the momentum boost in the quasi-PDFs can be interpreted as a Lorentz boost, recovering the light-cone PDF. The difference between quasi-PDFs and light-cone PDFs is O Λ 2 QCD /P 2 3 , m 2 N /P 2 3 and is calculable in continuum perturbation theory within the Large Momentum Effective Theory (LaMET) [12]. A successful research program on obtaining the PDFs using the quasi-PDFs method was developed since Ji's proposal, leading to theoretical and numerical advances . Recently, an exploratory study appeared on the strange and charm unpolarized PDFs [51] using ensembles with pion mass 310 and 690 MeV. However, the work only presents matrix elements in coordinate space. Other methods on extracting the x-dependence of distribution functions have been discussed . For an extensive review of the lattice calculations using the quasi-PDFs method, as well as other approaches to extract PDFs, see Refs. [81,82].
Lattice implementation. Based on the quasi-PDFs approach, the light-cone PDFs are obtained by the convolution of quasi-PDFs and the corresponding analytic expression for the matching kernel calculated in continuum perturbation theory. The quasi-PDFs are defined in momentum space as and are Fourier transform of hadronic matrix elements These matrix elements are calculable on a Euclidean lattice. In Eq. (3), the proton initial and final states, |N (P ) , carry the same momentum P = (P 0 , 0, 0, P 3 ), as the PDFs are obtained in the forward kinematic limit.
Here we focus on the helicity PDFs, ∆q ≡ g q 1 (x), and therefore we use the axial non-local operator, which contains a Wilson line W (0, z) of length z that guarantees gauge invariance. The bare matrix elements M(z, P 3 ) must be renormalized with an appropriate renormalization function, Z(z, µ), to remove divergences. We calculate Z(z, µ) using the RI -type prescription proposed in Refs. [16,17] which is applied at each value of z separately. We refer the Reader to Ref. [34] for notation. Due to the presence of the Wilson line in the operator, extracting the singlet renormalization functions is very challenging, as it involves a disconnected diagram. Here we use the non-singlet function indicated by Z(z, µ). We note that the difference between the singlet and non-singlet renormalization functions is expected to be small, as is the case of the local axial-vector operator [83]. This small difference has its origin to the fact that the difference between singlet and non-singlet arises to two loops in perturbation theory [84]. In addition to the logarithmic divergences and finite renormalization, the definition of Eq. (4) also removes the power-law divergence of the Wilson line. Z(z, µ) is obtained at an RI scale µ 0 . In our analysis, we convert to the MS scheme at a scale µ = 2 GeV. An additional conversion factor is used to bring Z(z, µ) in the modified MS-scheme [34]. Therefore, the scale dependence appears in the renormalized matrix element M R (z, P 3 , µ). While the matrix elements of local operators mix under renormalization [85], the non-local operators under study do not mix in the renormalization process, as discussed in Refs. [32,33,35]. This is because there is no additional non-local ultraviolet divergence in the quasi-PDF, an argument that holds to all orders in perturbation theory. However, the mixing occurs at the matching level and should be eliminated. To disentangle the singlet PDFs requires the matrix elements of the gluon PDFs, which is beyond the scope of this work. The nature of the mixing was also discussed earlier in Ref. [21] using the auxiliary field approach.
The most widely-used method to obtain the quasi-PDFs is via the discretized Fourier transform of Eq. (1). More recently, alternative reconstruction techniques are being explored [44,69,80,86,87]. In this work, we compare the standard Fourier transform, with the Bayes-Gauss-Fourier transform [87]. We find agreement between the two approaches, indicating that the behavior of the lattice results at the large-x region are not due to the discretization of the Fourier transform. We thus present results using the discretized Fourier transform.
As can be seen in Eq. (1), the quasi-PDFs depend on the nucleon momentum P 3 , which should be finite but large. This dependence is expected to be removed by the matching kernel which is calculated to a given order in continuum perturbation theory. The matching kernel for the quasi-PDFs approach has been extensively studied (see, e.g., Refs. [22][23][24][45][46][47][48][49]. In this work we employ the oneloop matching kernel in the modified MS-scheme, as defined in Ref. [34]. Note that we choose the factorization scale to be the same as the renormalization scale µ. The final step in extracting the light-cone PDFs is the application of the nucleon mass corrections, that have been calculated analytically in Ref. [13]. Numerical Methods. Obtaining M(z, P 3 ) is the most computationally demanding part of the calculation, as it contains a non-local operator, and must be calculated in the boosted frame. We perform the calculation including, for the first time, connected and quark-disconnected diagrams, for both the light and strange quark. In the light sector, we extract the isovector and isoscalar combinations, which are decomposed into the up and down quark helicity PDFs. The calculation is performed using an ensemble of two light, a strange and a charm quark (N f = 2 + 1 + 1) within the twisted mass fermion formulation with clover term. The lattice spacing is a = 0.093 fm and the lattice volume is 32 3 ×64 (L ≈ 3 fm). The pion mass is about 260 MeV and m π L ≈ 4.
The evaluation of the connected diagram uses the techniques outlined in Ref. [34], including the implementation of the momentum smearing method [88], and five stout smearing steps with parameter ρ = 0.15, on the Wilson line entering the operator. Both smearing methods contribute to the reduction of the statistical noise. We refer to Ref. [34] for the details. We use a total number of measurements N meas = 392, 1552 and 6320, for momenta P 3 = 0.41, 0.83 and 1.24 GeV, respectively. The sourcesink separation is t s = 0.94 fm for the lowest momentum and t s = 1.13 fm for the other two.
The evaluation of the quark-disconnected diagrams involves the computation of disconnected quark loops that have to be combined with the nucleon two-point correlators. The disconnected quark loop with Wilson line reads where D −1 q (x ins ; x ins + z) is the quark propagator, whose endpoints are connected by a Wilson line of length z. To reduce the stochastic noise coming from the lowmodes [89], we computed the first N ev = 200 eigen-pairs of the squared Dirac twisted-mass operator. From the eigen-pairs, the low-modes contribution to the all-to-all propagator can be exactly reconstructed and the highmodes contribution can then be evaluated with stochastic techniques. In particular, the stochastic evaluation of the disconnected loops is based on well-established techniques developed for local operators, such as hierarchical probing [90]. The latter allows for reduction of the contamination of the off-diagonal terms in the evaluation of the trace of Eq. (6), up to a distance 2 k . This is done using Hadamard vectors as basis vectors for the partitioning of the lattice. Here, the hierarchical probing algorithm has been implemented with k = 8 in 4-dimensions, leading to 512 Hadamard vectors. In addition, for the stochastic evaluation of the disconnected loops in this work we make use of the one-end trick [91,92] and fully dilute spin and color sub-spaces. We have recently employed successfully such methods in other studies of disconnected contributions. For more details see Refs. [83,[93][94][95].

Results for the connected and disconnected contributions.
For each value of the proton momentum, P 3 = 0.41, 0.83 and 1.24 GeV, we compute the two-point correlator for 200 source positions to reach a good statistical accuracy. We also take all spatial orientations of P 3 and W , that is, ±x, ±y, ±z. Moreover, both in the two-point and disconnected three-point functions we average over the forward and backward contributions. The total number of configurations analyzed is 330 for the two smallest momenta, and 480 for the largest one, bringing the total statistics to 66,000 and 96,000, respectively. Momentum smearing is applied for the two largest values To properly take into account the contamination of the excited states occurring at small source-sink separations t s , we compute the disconnected three-point correlators at t s = 0.75, 0.84, 0.94, 1.03, 1.13 fm, and perform a twostate fit analysis, following the procedure described in Ref. [93]. We find that the two-state fit gives results that are in agreement with those obtained form the plateau method analysis for t s = 1.13 fm. We will use the results from the plateau method for what follows. In Fig. 1  we show the real part of the bare matrix elements using t s = 1.13 fm for the disconnected contributions. The disconnected part of the isoscalar combination u + d, is smaller than the connected isoscalar contribution as expected. The real strange matrix element is about half as compared to the disconnected u + d. However, in both u + d and strange we clearly obtain a non-zero signal with the statistical uncertainties under control. The imaginary part of the bare disconnected matrix element is compatible with zero at each z, and is not shown. The matrix elements smoothly decay to zero and for z/a > 8 become compatible with zero. We note that the increase in the error for z = 8a in the disconnected part of the matrix element, i due to using hierarchical probing with length 2 k and k = 3. This is verified by repeating the evaluation of the disconnected diagrams with k = 2, and confirm that the same behavior occurs at z = 4a and its multiples, reflecting the limitation of the hierarchical probing technique when dealing with large lengths of the Wilson line. In taking the Fourier transform in Eq.
(1), we choose the cutoff z max such that the renormalized matrix element is compatible with zero. Since for the isoscalar and isovector matrix elements this occurs at different values of the Wilson line length z, we use different cutoffs z max for the two quantities. In particular, for the isoscalar case (the sum of connected and disconnected contributions) at P 3 = 1.24 GeV, we use z max = 7a, which is below the hierarchical probing length of 8a. While, for the isovector case, the matrix element is compatible with zero at z max /a = 9.
Two additional important issues need to be addressed in order to extract the PDFs, namely the dependence of the results on the momentum boost and the accuracy of the discrete Fourier transform. We examine these issues b considering the x∆d + (x) ≡ x ∆d + ∆d distribution, since the behavior is similar for the other two. To extract the ∆d + (x) distribution we apply renormalization and matching procedures separately on the isovector, isoscalar and strange quasi-PDFs. As mentioned above, we neglect the mixing with the gluon PDFs at the matching level.
In Fig. 2 we show the momentum dependence of x∆d + (x). We observe that, while when increasing the momentum from 0.41 GeV to 0.83 GeV there is a discrepancy in particular for large values of x, when we further increase the momentum to P 3 = 1.24 GeV, the results become compatible. This suggests that convergence has been reached within the limits of our current precision. In Fig. 2 we also show the dependence of the x∆d + (x) distribution on the cutoff z max adopted in the computation of the isoscalar and isovector quasi-PDFs. Despite the fact that when increasing z max , the resulting distribution tends to show more pronounced oscillations, the results for differrent z max all agree within uncertainties. In order to estimate the extent of the systematic effect due to the discretization and truncation of the Fourier Transform (FT), we employ the Bayes-Gauss-Fourier Transform (BGFT) [87]. As can be seen in Fig.  2, the distribution obtained with the BGFT technique is compatible with the standard reconstruction based on the discrete FT.
Flavor decomposition and comparison with phenomenology. The aim of this work is to obtain the flavor decomposition of the up, down and strange quark distributions, by combining the total isoscalar and isovector contributions at each P 3 value. In Fig. 3 we show our final results at P 3 = 1.24 GeV for |x|∆q + (x) ≡ |x| (∆q + ∆q), for q = u, d, s, and compare with the JAM17 data [7]. We find that x∆d + (x) and x∆s + (x) nicely decay to zero at x = 1. While x∆u + (x) is also zero at x = 1 and in agreement with the JAM17 results for x 0.6, it decays slower than the JAM17-determined distribution. On the other hand, we find a remarkable agreement for x∆d + (x) for the whole x region. The lattice determination of the strange distribution x∆s + (x) is much more precise as compared to the one determined from the global analysis. Although small is non-zero for small values of x. This result provides a valuable input for phenomenological studies. Conclusions. Results for the up, down and strange quark helicity PDFs of the proton, within lattice QCD are presented for the first time using the quasi-PDFs approach. We compute matrix elements with nucleon states boosted to maximum momentum P 3 = 1.24 GeV. We verify that the ground state matrix elements are well-determined by using one-and two-state fits, confirming that t s = 1.13 fm is sufficiently large to suppress excited state contributions at this level of precision. The matrix elements are renormalized non-perturbatively, and matched to the light-cone PDFs using one-loop perturbation theory. For the flavor decomposition of the light quark PDFs we take into account, for the first time, both connected and disconnected diagrams and compute the totally disconnected strange PDF. The final results on |x|∆q + are shown in Fig. 3, and are compared with the global fits of the JAM Collaboration. We find a remarkable agreement for the case of ∆d + for all values of x and for case of ∆u + for x < 0.6. We also obtain ∆s + much more precise that the phenomenological determination and show that is clearly non-zero for small values of x. This work paves the way for a determination of these helicity PDFs using ensembles simulated with pion mass, which we plan to do in the near future.
In the near future, a number of sources of systematic uncertainties will be explored, using the particular ensemble, along the lines of the analysis of Ref. [34]. Other effects is the implementation of the mixing matching matrix between quark and gluon PDFs, that requires knowledge of the gluon matrix elements of non-local operators. Systematic uncertainties requiring more than one ensemble include discretization effects, volume effects, and pion mass dependence. We plan to assess a proper determination of all sources of systematic uncertainties for the individual flavor PDFs in the future. Once systematic uncertainties are addressed and quantified, lattice results can provide useful input in the global fits for the strange PDFs, as well as the individual light-quark PDFs. This calculation is a first step towards achieving this goal.
ropean Unions Horizon 2020 research and innovation programme under grant agreement No 765048. This research includes calculations carried out on the HPC resources of Temple University, supported in part by the National Science Foundation through major research instrumentation grant number 1625061 and by the US Army Research Laboratory under contract number W911NF-16-2-0189. Computations for this work were carried out in part on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725. The gauge configurations have been generated by the Extended Twisted Mass Collaboration on the KNL (A2) Partition of Marconi at CINECA, through the Prace project Pra13 3304 "SIMPHYS".