Isovector Parton Distribution Functions of Proton on a Superfine Lattice

We study isovector unpolarized and helicity parton distribution functions (PDF) of proton within the framework of Large Momentum Effective Theory. We use a gauge ensemble, generated by the MILC collaboration, with a superfine lattice spacing of $0.042$ fm and a pion mass of $310$ MeV, enabling us to simultaneously reach sub-fermi spatial separations and larger nucleon momenta. We compare the spatial dependence of quasi-PDF matrix elements in different renormalization schemes with the corresponding results of the global fits, obtained using 1-loop perturbative matching. We present determinations of first four moments of the unpolarized and helicity PDFs of proton from the Ioffe-time dependence of the isovector matrix elements, obtained employing a ratio-based renormalization scheme.

The status of this field is well summarized in recent review papers [39][40][41][42]. All these calculations for the nucleon, so far, have been carried out with lattice spacing a > 0.08 fm.
Having small lattice spacing plays a crucial role in calculation of PDF within the LaMET framework. To suppress the target mass and higher twist corrections the hadron momentum, P z , should be large. But to avoid large discretization effects one must ensure aP z 1. Furthermore, to obtain lightcone-PDF from qPDF one needs perturbative matching, which, presently, is known only up to 1-loop order. Applicability of 1-loop perturbative matching can be guaranteed only for spatial separations zΛ QCD 1 and, therefore, demands use of fine lattices. The main goal of the present work is to study systematic of the PDF calculations within the LaMET framework by going to the extreme limit with the use of a superfine lattice having a = 0.042 fm. The lattice spacing used in this study is at least twice smaller than that used in any previous lattice calculations of the nucleon PDF. The unpolarized and helicity PDFs of the nucleon are well constrained through global fits to experimental results. Thus, we study the systematic of our calculations by comparing P z -and z-dependence of renormalized qPDF matrix elements with the same reconstructed from the well-known phenomenological PDFs using the LaMET framework.
The rest of the paper is organized as follows. In section II we discuss the general features of LaMET and our lattice setup. In section III we discuss the nucleon 2-point functions for large values of P z and the determination of the energy levels of the fast moving nucleon. Section IV is dedicated to the analysis of the nucleon 3-point functions and the calculations of bare qPDF. Section V describes the non-perturbative RI-MOM renormalization. The comparison of the lattice results on qPDF with the results of global analysis of unpolarized and helicity PDF is discussed in sections VI and VII, respectively. Different from RI-MOM renormalization, we discuss the analysis of ratios of nucleon matrix elements in Section VIII. Finally, section IX contains our conclusions.

II. LATTICE SETUP AND LAMET
In this paper, we report the results of a lattice QCD calculation using clover valence fermions on an ensemble of N f = 2 + 1 + 1 gauge configurations with lattice spacing a = 0.042 fm, with space-time dimensions of 64 3 × 192 and pion mass M π ≈ 310 MeV in the continuum limit. The gauge configurations have been generated using Highly Improved Staggered Quarks (HISQ) [43] by the MILC Collaboration [44]. The gauge links entering the clover Wilson-Dirac operator have been smeared using hypercubic (HYP) smearing [45]. We used tree-level tadpole improved result for the coefficient of the clover term and the bare quark mass has been tuned to recover the lowest pion mass of the staggered quarks in the sea [46][47][48][49]. We use only one step of HYP smearing to improve the discretization effects, since it is possible that multiple applications of smearing could alter the ultraviolet results for the PDF. We use multigrid algorithm [50,51] in Chroma software package [52] to perform the inversion of the clover fermion matrix allowing us collect relatively high statistics sample. We collected a total of 3258 measurements using 6 sources per configuration and 543 gauge configurations. In the following, we elaborate on the steps of our computation.

A. Nucleon two-point correlators
The two crucial components of the lattice computation are the two-point function and the three-point function involving boosted nucleon and the qPDF operator. The two point function for the nucleon boosted to spatial momentum P is the standard operator C 2pt (t s ) = N s (P, t s )N † s (P, 0) ,N s (P, t) = x abc u (s) wherex = (x, t) and t s is the source-sink separation along the Euclidean time direction. The index 's' refers to the kind of quark smearing that is applied to improve the signal-to-noise of the boosted nucleon states. We either used point quark operators ψ(x) or we used the Gaussian momentum smeared [53] for the quark fields, ψ (s) (x) that enterŝ N s , where k is the momentum of the quark field, U j (x) are the gauge links in theĵ direction, and α is a tunable parameter as in traditional Gaussian smearing. The quark momentum should be chosen such that the signal-to-noise ratio is optimal for the given nucleon momentum. Naively one would expect that |k| should be one third of the nucleon momenta [53]. For this particular study, we use j = z and k z = 4π/L, and a large Gaussian-smearing parameter α = 10. Such a momentum source is designed to align the overlap with nucleons of the desired boost momentum, and we are able to reach higher boost momentum for the nucleon states with reasonable signals. In the nucleon two-point correlators, we can study multiple values of the nucleon momentum, P = {0, 0, P z } with without a significant increase in computational needs. These values of n z from 1 to 6 correspond to P z = 0.46, 0.92, 1.38, 1.84, 2.31 and 2.77 GeV in physical units respectively. We either used smeared fields for both the source and sink, which we refer to as SS, or smeared fields only for the source and point fields for the sink which we refer to as SP in the rest of the paper.

B. Nucleon three-point function
The three point function we compute is of the form C 3pt (t s , τ ) = P P N s (P, t s )O Γ (z; τ )N † s (P, 0) (4) where O Γ (z; τ ) is the u − d isovector qPDF operator W z is the straight Wilson line along the spatial z-direction, connecting lattice sitesx andx + z . The Dirac Γ used will determine the quantum numbers of the PDF -Γ = γ t for the unpolarized case and Γ = γ z γ 5 for the longitudinally polarized case. The projector operator, P P, is given by P P = 1+γt 2 for the unpolarized case and P P = iγ z γ 5 1+γt 2 for the longitudinally polarized case, respectively. We only use smeared quark sources for the computation of C 3pt . In order to reduce the computational cost, we only computed the C 3pt for two large values P z = 1.84 and 2.31 GeV, and for source-sink separations t s = 16a, 18a, 20a.
C. Extraction of nucleon matrix element and perturbative matching to PDF Using the three-point and two-point functions whose calculations are described above, we can extract the bare matrix element h(z, P z , Γ) = P z |O Γ (z)|P z , formally in the infinite source-sink separation t s limit of their ratio To obtain the matrix element h(z, P z ) from the above ratio, we calculate the nucleon three-point function with insertion of O Γ (z) operator at three nucleon three-point source-sink separations, approximately t s = 0.67, 0.76, 0.84 fm, and describe its t s -and τ -dependence through 2-and 3-state ansatz. In Sec. IV, we describe our extraction of bare matrix element from various extrapolations in detail. The next step of the computation is the renormalization of the bare matrix element h. One possible choice for O Γ is O γz . However, for this case of Γ = γ z there is a mixing with the quark bilinear operator containing the unit matrix, Γ = 1 if Wilson fermions are used [20,21,54]. This mixing is absent if we use Γ = γ t , and we will use this choice for the unpolarized PDF in this study. One way to perform the renormalization procedure on the lattice to use RI-MOM scheme [20,22], where in the renormalized matrix element is defined as where Z is the non-perturbatively determined RI-MOM renormalization factor at a renormalization scale whose norm is µ R and the z-component is p R z . The dependence on p R z arises because the z-component of the momentum now plays a special role. We will discuss the details of the RI-MOM renormalization in section V. We will also consider an alternate ratio scheme that has a well defined continuum limit in Sec. VIII. Here, the multiplicative renormalization factor Z ratio (z) can be taken as the hadron matrix element at a different fixed momentum P z i.e., Z ratio (z) = (h(z, P z , Γ)) −1 .
After the RI-MOM renormalization one obtains the renormalized matrix element h R (z, P z , µ R , p R z ), from which we can define the qPDF as a function of Bjorken-x From this formula it is clear that h R (z, P z , µ R , p R z ) can be considered as the coordinate space qPDF. For finite momentum P z ,q(x, P z , µ R , p R z ) has support in −∞ < x < ∞. Unlike the physical PDF, which is frame independent, the qPDF has a nontrivial dependence on the nucleon momentum P z . When the nucleon momentum P z {M, Λ QCD } with M being the nucleon mass, the qPDF in RI-MOM scheme can be matched to the PDF defined in MS-scheme, q(x, µ) through the factorization theorem [12,13,25], where r = (µ R /p R z ) 2 and C is the perturbative matching coefficient, O(M 2 /P 2 z ) is the target-mass correction due to the non-zero nucleon mass, and O(Λ 2 QCD /P 2 z ) stands for higher-twist contributions. The flavor indices of q,q, and C are implied. In what follows we will discuss the non-singlet case, and therefore, mixing with gluon and sea-quark PDFs is absent in the above formula. We use 1-loop expression of the kernel C. (The 1-loop matching including for the singlet case also has been worked out in Ref. [55,56].) The matching kernel C(x, r, P z /µ, P z /p R z ) for Γ = γ t was derived in Ref. [24] and depends on details of the RI-MOM scheme. It can be written in the following form The subscript '+' stands for the plus-prescription. Both the functions, f 1,Γ and f 2,Γ,P , depend on the choice of the Γ in the operator insertion [24]. On the other hand, f 1,Γ is independent of the projection operator (P) used in defining the RI-MOM renormalization condition, but f 2,Γ,P is different for different choices of the RI-MOM renormalization condition [24]. We also note that it is also possible to convert h R (z, P z , µ R , p R z ) to MS-scheme and define the corresponding qPDFq(x, P z , µ) that then can be directly matched to MS PDF [22].
To study the longitudinally polarized quark PDF one can use Γ = γ z γ 5 or Γ = γ t γ 5 . In the case Γ = γ z γ 5 there is no mixing with quark bilinear operators with Γ = 1 [22]. Therefore, we will use this choice to study the longitudinally polarized quark PDF and qPDF. The bare matrix element of O γzγ5 can be renormalized using RI-MOM scheme and then match to PDF in the same manner as this was done for unpolarized. The RI-MOM renormalization for the longitudinally polarized case will be discussed in section V, while details of the matching procedure, including the formulas for f 1 and f 2 functions will be give in section VII.

III. ANALYSIS OF THE NUCLEON TWO-POINT FUNCTION
For the extraction of the qPDF matrix element of the nucleon at large momenta it is important to understand the contribution of different energy states to the nucleon two-point correlation function. We calculated nucleon twopoint function using smeared source and smeared sink (SS correlator), as well as smeared source and point sink (SP correlator), for seven values of the momenta aP z = 2π/L·n z , n z = 0, 1, 2, 3, 4, 5 and 6. From the two-point correlators, C i 2pt (t s , P z ), i = SS or SP, we define the effective mass Our results for the effective masses are shown in Fig. 1 for the SP and SS correlators. The effective mass should approach a constant corresponding to the ground state energy E 0 (P z ) at sufficiently large t s . The momentum dependence of the ground state energy is expected to be described by the dispersion relation E 0 (P z ) = P 2 z + M 2 , with M being the nucleon mass. Therefore, in Fig. 1, we show the expected ground state energy at different P z obtained from the dispersion relation as horizontal lines at the right for comparison. Along  with the expected asymptotic values at large t s , we also show the t s -dependence of the the effective mass based on an effective two-state fit to the two-point function, as we will explain shortly. Indeed, we see that the effective masses approach the corresponding values. The effective masses corresponding to the SP correlator reach a plateau at a slightly larger t s than the SS correlators. On the other hand, at small t s , the effective masses for the SP correlators are smaller than those for SS correlators. This implies that the contribution of the excited states is smaller for the SP correlator, for which a plausible reason could be that the different excited states contribute with different signs to the correlator. Thus, even though the ground and the excited state energies are the same in the SP and SS correlators, the two are affected differently by the higher excited states, which we can take advantage of to obtain the excited state spectrum reliably. In order to determine the energy levels, we fit the spectral decomposition of C 2pt (t s ), truncated at N state to the two-point function data over a range of values of t s between [t min , 32a]. Since the lattice extent in the time direction is 192, we did not find any effect of lattice periodicity in this range of t s to be important. We performed this fitting with one-state (N state = 1), two-state (N state = 2), and three-state (N state = 3) Ansätze. The ground state energies, E 0 from the fits of SS correlators for n z = 3 and 4 are shown in left panels of Fig. 2 as function of t min , where t min indicates that only C 2pt (t s > t min ) have been fitted. Similar results were obtained at the other values of the momenta. The horizontal lines in the figures correspond to the results from the dispersion relation for E 0 . The single exponential fits give a good description of the SS correlator for t min > 11a, while two exponential fits give stable results for the ground state energy already for t min > 5a. We found the determination of the excited state energies from the SS correlators to be more problematic than from SP correlators. The excited state energy for SS is not well-constrained by simple two exponential fits, and it is also not very stable with respect to the variation of t min . Since the SP and SS correlators receive different contributions from excited states, we performed a combined analysis of them to obtain more reliable results for the excited state energies. Since we were able to obtain the ground state energy E 0 reliably from one or two exponential fits to both the SS and SP correlators and they agree with the expectation from the dispersion relation well, we used E 0 as a prior to performed more stable two-exponential fits. The results from the two-state exponential fits, with E 0 as prior, for n z = 3 and n z = 4 are shown in middle and right panels Fig. 2 for the SP and SS correlators, respectively. For the SP correlators, the excited state energy E 1 seems to approach a plateau smoothly for t min > 13a. It is interesting to note that, empirically, we observe the values of the plateaus agree with the dispersion relation E 1 (P z ) = P 2 z + E 1 (P z = 0) 2 , which are shown as the horizontal lines. While being an interesting observation, such a stringent identification of this state is not important to our analysis and requires further studies to rigorously establish this. For the SS correlator, E 1 develops a pseudo-plateau for 5a < t min < 10a and it relaxes to the true plateau (i.e., as identified from the SP case) for t min > 11a. For n z = 4, it is actually difficult to identify the true plateau. To model the the excited state contributions to R(τ, t s ) in the range 0 < τ < t s /2, with t s = 16a, 18a, 20a, one might consider using the well-determined values of E 0 and E 1 from the SP correlator at large t s . However, as we will demonstrate now, such choices provide a less accurate description of the SS two-point function in the range 5a < t s < 10a. A better description of the excited-state contributions to the C 2pt (5a < t s < 10a) can be obtained by using the effective pseudo-plateau value of E 1 in the range of 5a < t min < 10a. Since, we observe E 1 to be well-described by a particle-like dispersion relation for sufficiently large t s , we perform three-state fits for both SP and SS correlators by imposing a prior on E 1 as well, using its best estimate from the SP correlator. The results are shown in middle and right panels in Fig. 2. We see that with the prior-based three-state exponential fits, we can obtain stable results for the first excited state energy, E 1 (P z ), already for relatively small t min which agrees with the dispersion relation value that we input via the prior. The value of the second excited state is also shown in Fig. 2 and it roughly agrees with the values of E 1 from the two-exponential fit (with prior only on E 0 ) at smaller t min . Since the value of E 2 is quite large, the third exponential probably corresponds to a combination of several excited states. In Fig. 1, we show the 1-σ bands for the effective mass corresponding to: (1) Two-state fit that uses values of E 0 and the true value of E 1 ; (2) two-state fit obtained by setting E 1 to be the effective value in the range from 5a < t min < 10a; (3) three-state fit that we described above. We find that the curves (2) and (3) agree quite well with each other in the range of 5a < t s < 10a and they extrapolate in the similar fashion to the asymptotic value E 0 . However, the curve (1) fails in capturing the data in the range 5a < t s < 10a. Since for our three-point calculations the source-sink separations were chosen to be t s = 16a, 18a, 20a, we must model the effective excited state contributions to the three-point functions in the range 0 < τ < t s /2. Thus, through this analysis on SP and SS correlators, we numerically demonstrate the usage of an effective value of E 1 in the range of 5a < t s < 10a that is higher than the true value of E 1 is justified, and is the best extrapolation one could perform for the extraction of bare matrix elements in the absence of enough data to perform a three-state fit.
Let us now, summarize the analysis of the nucleon two-point function. Using boosted Gaussian sources we were able to extract ground state energy levels up to momenta 2.7 GeV from SP and SS correlators. The ground state energy dependence on P z seem to follow the continuum dispersion relation. Using this fact, we performed prior-based fits using the energy from the dispersion relation as a prior and extracted the excited state energies as function of P z . For SP correlator the extracted value of E 1 agrees well with the one from the dispersion relation. We show this in the left panel of Fig. 3. Furthermore, we were able to extract an effective third energy level. These results are shown as blue points (E 1 ) and orange points (E 2 ) in the left panel of Fig. 3. Our results for the energy levels obtained from SS as function of P z are summarized in right panel of Fig. 3. Here, the effective values of E 1 from the pseudo-plateau and the true values are both showed. The main point of the elaborate analysis is that even though a third excited state contributes in the relatively shorter range of t s we use in the paper, it possible to describe the the SS correlator very well by a two state form with an effective value of E 1 , which is larger than the energy of the physical excited state. Further details on the analysis of the two-point functions are provided in Appendix A.

IV. NUCLEON THREE-POINT CORRELATORS
In order to obtain the nucleon qPDF matrix element we consider the ratio of the 3-point function to 2-point function, R(z, P z ; t s , τ ), at different source sink separations, t s , and operator insertion, τ . At fixed (z, P z ) we are interested in fitting the (t s , τ )-dependence as expected from the spectral decomposition of R. If only two states contribute to the correlation functions, the dependence of this ratio on τ and t s is given by the following form: Here B 0 is the desired matrix element h, and ∆E = E 1 − E 0 . Generically, B 1 and B 2 are independent fit parameters, except at z = 0, where B 1 = B 2 . If we assume that the terms proportional to A 1 are small, the denominator can be expanded to leading order to obtain a simpler form Finally, if the term proportional to B 3 is also small compared to other terms, we get an even simpler expression that depends only on three parameters, B 0 , B 1 and B 2 : For each (z, P z ), we fitted the (t s , τ )-dependence of R(z, P z ; t s , τ ) to Eqs. (14,15,16) and determined B 0 in each case. In all these fits we used a fixed value of ∆E(P z ) = E 1 (P z ) − E 0 (P z ), with the pseudo-plateau values of E 1 (P z ) and the ground state energies E 0 (P z ) determined from the 2-point SS correlation function, as shown in Fig. 3(right).
In the following, we discuss the ratio of the 3-point function to 2-point function, R(z, P z ; t s , τ ), and the corresponding fits for Γ = γ t and n z = 4. In Fig. 4, we show the lattice data on this ratio, together with the fit results for two representative values of z, namely z = 0 and z = 8a. The t s dependence of the lattice results is small compared to the statistical errors. In particular, the difference between t s = 12a and t s = 16a data is quite small. This means that the contribution of the excited states is not large even though the source-sink separation is below 1 fm. Given that the t s -dependence of the ratio is small, it is natural to set B 3 = 0 since the term is suppressed by e −ts∆E , and perform fits using R fit 3 (t s , τ ). We performed fits of the lattice results using the three different fit forms above and the value of ∆E The real part ratio of the three-point function to two-point function at z = 0 (left) and z = 8a (right) as a function of τ and for different ts. The red, orange and gray bands correspond to the bare matrix elements B0 extracted from fits to R fit 1 , R fit 2 and R fit 3 , respectively.   obtained from two-state fits of the 2-point functions with t min = 6a. We used operator insertion time τ τ min = 2a in the fits. The matrix elements, B 0 , were obtained using the three fit functions agree within errors. The real and imaginary part of the ratio R should be symmetric and anti-symmetric with respect to z = 0 at any fixed t s and τ . In general, our lattice data is compatible with this expectation, but in certain cases we see some tension between the values calculated for positive and negative z. This is probably due to limited statistics. Hence, we reduce this effect by symmetrizing and antisymmetrizing the real and imaginary part of the ratio with respect to z = 0, respectively. The z-dependence of the bare matrix elements is shown in Fig. 5 for all three types of fits. We see again that all three fits give the consistent results within errors. Since B 0 obtained from R fit 3 and R fit 1 are consistent with that obtained from R fit 1 , but with larger errors, in the following we will focus on the results obtained from R fit 1 . We also carried out additional checks for any systematic effects, as discussed below.
We performed several checks to understand the systematic effects in R fit 1 . First, we studied the dependence of the extracted matrix element on τ min and found no significant dependence on it. Second, we performed the fits using only a single source-sink separation t s and compared the corresponding results from the three values of t s . Interestingly, the matrix elements calculated for source sink separation t s = 16a, 18a and 20a agree within errors, though the t s = 20a results have very large errors. We also studied the variation of the extracting matrix element on ∆E by using E 1 obtained using different values of t min . We found no significant variation. Finally, we used the summation method to obtain the matrix element. This determination has very large statistical errors but it is still compatible with all other determinations. The above checks of systematic effects are discussed further in the Appendix B. We performed a similar analysis for the three-point functions corresponding to the helicity qPDF, i.e. for Γ = γ z γ 5 . Details of those analysis also are discussed in the Appendix B.

V. NONPERTURBATIVE RENORMALIZATION
We calculated the nonperturbative renormalization of the qPDF operator in the RI-MOM scheme using off-shell quark states in the Landau gauge [20,22]. The matrix element of O Γ (z) in an off-shell quark state |p is where n µ = (0, 0, 0, 1) is the unit vector along the z direction, and the summation is over all lattice sites w. The quark propagators are defined as and γ 5 is inserted on both sides of S † (p, w + zn) in Eq. 17 to get the necessary propagator y e −ipy ψ (y)ψ(w + zn) . For the unpolarized qPDF, we use the RI-MOM renormalization constant defined via where Λ tree (p, z, γ t ) = γ t e −izpz is the tree level matrix element in the momentum space. Furthermore, P = γ t − (p t /p x )γ x is the projection operator corresponding to the so-called minimal projection, where only the term with the Dirac structure proportional to γ t is kept [23,24]. Hence, we use the subscript 'mp' for the renormalization constant. The renormalization constant Z mp (z, p R z , a −1 , µ R ) depends on the lattice spacing a, as well as the other two scales p R z and µ R . We followed a similar procedure for the longitudinally polarized case, where the RI-MOM renormalization constant is defined as with Λ tree (z, p z , γ z γ 5 ) = γ z γ 5 e −ipzz . The projection operator P in this case was chosen to be P = γ 5 γ z /4. We calculated the non-perturbative RI-MOM renormalization constants in Landau gauge. The calculations were performed using 14 gauge configurations. The relative uncertainties of the renormalization constants for z = 0, 16, 32 are 0.02%, 1% and 10%, respectively. Such precision is much better than that of our nucleon matrix elements with the same z, so it is enough at the present stage. We used the following values of the momenta for the off-shell quark state: ap = 2π L (5, 5, 5, 0), 2π L (6, 2, 1, 17/3) and 2π L (7, 4, 3, 1/3), L = 64 being the spatial size the of the lattice. These momenta correspond to µ R = |p| = 3.99 GeV, 3.94 GeV and 3.97 GeV, i.e. to µ R ∼ 4 GeV within 1.5%. Since all the spatial directions are equivalent, each of them could be considered as the z-direction and, therefore, with the above choice of the three momenta we have p R z = 0.46 × {0, 1, 2, . . . , 7} GeV.  The renormalization constant at z = 16a ≈ 0.67 fm is plotted in Fig. 6 as a function of p R z . Due to the linear divergence, the renormalization constant can be far from 1 at a large z ≈ 0.67 fm, making the nonperturbative renormalization unavoidable. Fig. 6 also shows that the renormalization constant will be sensitive to the value of p R z , while such a dependence should be canceled by the matching in the continuum. We will consider the residual p R z dependence in the final PDF prediction as a systematic uncertainty.
Having determined the renormalization constants Z mp and Z γzγ5 we obtained the renormalized matrix elements, i.e., coordinate space qPDF. For the unpolarized case, and for longitudinally polarized case, In the above equations, Z q is the quark wavefunction renormalization factor. In Fig. 7 we show the renormalized matrix elements, modulo the factor Z q , in the RI-MOM scheme at p R z = 0, µ R = 4 GeV. We find that the errors are large. We can achieve substantial error reductions at z = 0, by redefining the renormalized matrix elements as The errors of the matrix elements for z = 0 are reduced due to the strong correlations between z = 0 and z = 0 matrix elements for each gauge configurations. The effectiveness of this procedure in can be seen from Figs. 12 and 15. As one can see the error reduction due to this division is very significant. In fact, with this method, the errors are reduced enough that the z-dependence of the matrix element is well constrained also for n z = 5. Since for the extraction of the qPDF we are only interested in the z-dependence of the matrix element, and we know that the unpolarized isovector nucleon matrix element at z = 0 is the isospin of the nucleon, which is unity (in our convention, c.f. Eq. 5 ) after renormalization, we can consider the above improved ratio of renormalized matrix elements. However, the effect of taking this ratio is not trivial in the case of the matrix element of the helicity qPDF-the value of the renormalized matrix element at z = 0 should be g A ≈ 1.3; this procedure is equivalent to studying a helicity PDF with the first moment normalized to unity, i.e., in a normalization where g A = 1.

VI. UNPOLARIZED PDF: PERTURBATIVE MATCHING AND COMPARISONS WITH h R (z, Pz)
In this section, we will discuss how the renormalized coordinate-space qPDF, h R (z, P z , µ R , p R z ), can be related and compared with phenomenological unpolarized nucleon PDF, such as the CT18 [57] and NNPDF3.1 [58], extracted from the global analysis of experimental data. The unpolarized quark PDF in the valence region is well constrained through global analysis. Therefore, it is natural to start from these phenomenological PDFs as a function of Bjorken-x, use the perturbative matching to reconstruct the corresponding coordinate-space qPDF as a function of z for different P z values, and compare with our results for h R (z, P z , µ R , p R z ). The reason for comparing in the z-space, rather than constructing the x-dependent PDF from our h R (z, P z , µ R , p R z ) and then comparing with the phenomenological PDFs, is the following: As can seen from Figs. 12 and 15, h R (z, P z , µ R , p R z ) is quite noisy for z 0.5 fm. Thus, the Fourier transformation which is needed to calculate the qPDF in x -space is difficult to perform. Similar approach also had been used for pion PDF [35].
Even at the leading α 0 s order the qPDF and the PDF differ due to the trace term in the small z-expansion [12,25].   This difference was explicitly calculated in Ref. [59]. In the context of DIS, such corrections have been studied long ago [60], and are known as target-mass corrections. Following Ref. [59], we introduce the target-mass corrected PDF is the usual PDF that corresponds to P z → ∞. In our analysis we use two sets of NNLO PDF for the u and d quark and anti-quark distributions, the CT18 [57], and NNPDF3.1 [58], evaluated at scale µ = 3.2 GeV. If the matching was known to all orders of perturbation theory, the prediction for real space qPDF should have been independent of the value of µ at which the PDF was evaluated. Since the matching only known to 1-loop order, we chose a scale µ = 3.2 GeV that is of the same order of the other momentum scales used in our computations and, thereby, avoided corrections due to large logarithms. The lightcone quark PDF for u quark is calculated as q u (x) = u(x), x > 0 and q u (x) = −ū(−x), x < 0. The isovector nucleon PDF, q u (x) − q d (x) is shown in Fig. 8.
In Fig. 8, we also show the target-mass corrected isovector nucleon PDF for the two momenta used in our study, namely 1.84 GeV and 2.31 GeV. We see from the figure that target-mass correction is small for the values of P z used in this study. The CT18 and NNPDF3.1 results for isovector nucleon PDF agree well for positive x, while there is a significant difference between the two for −0.1 < x < 0. Using the target-mass corrected NNPDF3.1 isovector nucleon PDF obtained from Eq. (24) and the 1-loop matching to RI-MOM we obtained the corresponding qPDF for P z = 1.84 GeV and P z = 2.31 GeV, µ R = 4 GeV, and p R z = 0, 0.47, 0.93 GeV. The functions f 1,γt and f 2,γt,mp in Eq. (11) for the 1-loop matching to RI-MOM scheme with minimal projection were taken from Eq. (28) Eq. (31) of Ref. [24]. Fig. 9 shows comparisons of the NNPDF3.1 with the corresponding qPDFs. In these comparisons α s was evaluated at scale µ = 3.2 GeV, which resulted in a value α s = 0.25. We see significant differences between the PDF and qPDF. For large positive x, the qPDF is larger than PDF, while for negative x qPDF can turn negative. The qPDF strongly depends on the choice of the RI-MOM scales. It is possible to choose the RI-MOM scale such that the qPDF is negative for x < −0.2, even though the PDF is positive.
By Fourier transforming the CT18 and NNPDF3.1 target-mass corrected qPDFs with respect to x we obtained the corresponding distributions as a function of the so-called Ioffe-time, zP z , i.e. the corresponding ITDs [61]. Since the matching is only up to 1-loop order, the scale entering α s is not fixed. We considered three choices of the scale for α s , namely µ/2, µ, 2µ. The corresponding variations in the ITDs can be considered as estimates of the perturbative uncertainties, and are shown as bands in Fig. 10. In the same figure, also we compare with the lattice results for the ITDs in RI-MOM renormalization, at the renormalization scales of µ R = 4 GeV and p R z = 0 GeV. Albeit large errors, for both values of P z the real parts of the ITDs compare well at least up to zP z 5. However, lattice results for the imaginary parts of ITDs undershoot the phenomenological ITDs even for zP z 2.
Albeit the significant difference between CT18 and NNPDF3.1 PDFs in the negative-x region (c.f. Fig. 8), Fig. 10 do not show any visible difference in their corresponding ITDs. To understand this better, Fig. 11 we explore these ITDs in an extended range of Ioffe-time. The difference between the PDFs in the negative-x region is only reflected in < 10% difference in the imaginary part of the ITDs for zP z > 25, with essentially showing no difference in the real part ITDs even up to zP z = 50.
To explore the dependence of the lattice results on the choice of RI-MOM scale p R z and the range of validity of the 1-loop matching, in Fig. 12 we show comparisons between the qPDFs as a function of z obtained in the lattice calculations and from the global analysis of PDF for two different choices of renormalization scale, namely p R z = 0, 0.93 GeV. Very little dependence on the p R z was observed. While the real part of the the qPDF obtained from the global analysis agrees with the lattice results up to z ∼ 1 fm within relative large errors, for the agreement is limited only for z 0.2 fm. For P z = 2.31 GeV the agreements seem to extend to larger values of z, partly because of larger errors. However, it is encouraging that the central value seems to shift towards the global analysis results as P z is increased from 1.84 GeV to 2.31 GeV. In any case, at large z, we see clear tension between the imaginary part of the lattice qPDF lattice and the results of global analysis. This suggests that the range of applicability of 1-loop matching is perhaps limited to z 0.2 fm in the case of the nucleon. It remains to be seen if this agreement gets better with the addition of higher-loop corrections, or this observed discrepancy arises because of contamination of higher-twist effects at larger z. This observation has an important implication for our ability to described the x-dependence of PDF within the LaMET framework. For example, if the 1-loop perturbative matching works only for z 0.2 fm, reliable calculations of nucleon PDF down to x 0.1 will need P z 10 GeV.

VII. HELICITY PDF: PERTURBATIVE MATCHING AND COMPARISONS WITH ∆h R (z, Pz)
Our analysis of helicity qPDF closely follows the analysis performed in the unpolarized case, namely we start from the helicity PDF obtained in global analyses, reconstruct the corresponding target-mass corrected qPDF, and then compare with the lattice results. The helicity PDF have been extracted from the global analysis by NNPDF collaboration using DIS, inclusive W ± and jet production data from RHIC, as well as the open charm data from COMPAS resulting in NNPDFpol1.1 [62]. The JAM collaboration used the DIS and SIDIS data in their global analysis, combined with e + e − data to constrain the fragmentation functions at NLO [63]. The resulting PDF parameterization is called JAM17. In Fig. 13, we show the isovector helicity PDF ∆q u − ∆q d . The positive-x region corresponds to quark contribution, while the negative-x region corresponds to anti-quark region. The target-mass corrected helicity PDF, ∆q (x, P z ), was obtained from helicity PDF, ∆q(x), following Ref. [59]: where c = M 2 /4P 2 z , f ± = √ 1 + 4c ± 1, and for the integration limits +∞ (-∞) correspond to x > 0 (x < 0). Although the matching for helicity qPDF has not been explicitly presented in the literature before, it was straightforwardly deduced from the results presented in Ref. [24]. The key observation here was the fact that, owing to the chiral symmetry, for a mass-less quarks in 1-loop perturbation theory Tr[γ 5 γ z Λ(p, z, γ z γ 5 )] = Tr[γ z Λ(p, z, γ z )]. Thus, the 1-loop matching of the helicity qPDF in the RI-MOM scheme with minimal projection is same as that for the unpolarized qPDF with Γ = γ z (instead of the Γ = γ t used before), and with the RI-MOM renormalization condition corresponding to the projection operator P = γ z (instead of the minimal projection). The 1-loop matching for the Γ = γ z operator is known for two different RI-MOM projections, the minimal projection and the / p projection, corresponding to P = γ z − (p z /p x )γ x and P = / p/(4p z ), respectively [24]. The function that depends on the RI-MOM projection operator, i.e. f 2,γz,γz , entering the matching coefficient in Eq. 11 was simply deduced from these known results. The Lorentz structure of Λ(p, z, γ α ) for a general γ α , α = x, y, z, t is given by and f 2,γz,mp =f t +f z and f 2,γz,/ p =f t +f z +f p [24]. Here, the subscript '+' refers to the standard plus-prescription and ρ = −p 2 /p 2 z . The functionsf t ,f z andf p have been calculated in Ref. [24], and we use the same notations here. Therefore, for the case of P = γ z the RI-MOM projection-dependent function is given by Thus, for the helicity qPDF the 1-loop matching RI-MOM function in the minimal projection scheme is the same as in Eq. 11, but with f 2,γz,mp given by Eq. 27, and f 1,γz , f 2,γz,mp and f 2,γz,/ p are given by Eqs. (A6-A8) of Ref. [24]. Using the matching discussed above, we can obtain the isovector helicity qPDF from the target-mass corrected NNPDF1.1pol and JAM17. As before, the 1-loop matching we used α s evaluated at scale µ = 3.0 GeV, and the scale was varied between µ/2 to 2µ to estimate the the scale uncertainty. We found noticeable difference between the isovector helicity PDFs and the corresponding qPDFs: for positive x the qPDFs are larger than the PDFs, while for x < 0 it is mostly negative for the values of p R z and µ R considered by us. By Fourier transforming the qPDFs we obtained the isovector helicity ITDs and compared it with our lattice results in Fig. 14. Since we normalized our lattice results by the value of matrix element at z = 0, we normalized the phenomenological ITDs by dividing with g A = 1.25. Within the large statistical errors, we did not find a significant P z dependence of the lattice results. While the real parts of the lattice results agree with that obtained from the phenomenological PDFs up to zP z 3, the imaginary parts do not agree quantitatively, but found to be qualitatively similar. We also explored the dependence of our result on the choice of RI-MOM scales. In Fig. 15, we compare the qPDFs for µ R = 4 GeV, and p R z = 0 GeV and p R z = 0.93 GeV. From the figure, we see that the comparison between the results of lattice calculation, as well as the global analyses are not sensitive to the choice of the renormalization scales. For both values of P z , the agreement between the lattice and the global analyses extends to values of |z| of about 0.3 fm for the real parts, but not for the imaginary parts. In the next section we will discuss how these disagreements show up in the moments of the PDFs.

VIII. MOMENTS OF PDF FROM RATIO OF IOFFE-TIME DISTRIBUTIONS
In the previous sections, we analyzed the boosted nucleon matrix matrix elements renormalized in the RI-MOM scheme and matched it to the PDFs in the MS scheme. Due to the multiplicative renormalizability of h(z, P z , γ t ) and h(z, P z , γ z γ 5 ), we can form well-defined renormalized quantities by taking the ratios of matrix elements at two different momenta P z and P z as The second factor on the right hand side of the above definition normalizes the z = 0 matrix element to unity, as we did in the case of the RI-MOM scheme. The choice P z = 0 in the ratio is usually referred to as the reduced Ioffe-time distributions [16], and one should think of P z = 0 as a generalization of this choice. Here, we take P z = 2.31 GeV and P z = 1.84 GeV, respectively. Since both P z , P z > Λ QCD and the nucleon mass, we expect this ratio to be simply described by the leading twist expression [25], Following Ref. [59], the target-mass corrected unpolarized PDF moments x n can be obtained by relation: and for the helicity case, where C i n is the binomial function, c = M 2 /4P 2 z . In Eq. 29, c n (µz) is the 1-loop order Wilson coefficients in the MS scheme. The Wilson coefficients describes the z dependence of the twist-2 local operator associated with the n th moment of the PDF, x n (µ), in the MS scheme and at a factorization scale µ. As in our RI-MOM analysis, we will use µ = 3.2 GeV for the unpolarized case and µ = 3 GeV for the helicity case in the following analysis. Now, we can perform an independent analysis that avoids the usage of RI-MOM procedure completely and compare the outcome to the prediction for M(z, P z , P z , Γ) from the knowledge of NNPDF and CTEQ PDF moments. We perform such a comparison in Fig. 16. For this, we used the values of x n (µ) up to an order n = n max for NNPDF31 in Eq. 29, , and the complete result for CT18, to obtain the phenomenological expectation for the ratio M(z, P z , P z , γ t ). The results obtained by using the truncation order n max = 2, 3, 4, 20 using the NNPDF31 values for x n are shown as different colored bands in Fig. 16. It is clear that inclusion of up to n max = 20 moments is sufficient for convergence to the correct PDF within z 0.5 fm. For z < 0.3 fm, which is where the lattice data has a good signal to noise ratio, we find that N = 4 is sufficient to describe the lattice results. This gives us an idea of which moments are being probed by our lattice data at different z. We observe some discernible differences between the phenomenological expectations and our lattice M(z, P z , P z , γ t ) for z > 0.2 fm, as we also observed in the case of RI-MOM scheme in Fig. 10. To understand this, we estimate the values of the moments x n that best describe our lattice data. To avoid overfitting the data, we truncate the expansion in Eq. 29 at most by n = 4. In order to avoid lattice artifacts that might be present for z of the order of lattice spacing, we fit the data only from z = 2a to a value z max . The variation of the best fit values of x n with z max is a source of systematic error. In Fig. 17, we show the z max dependence of our estimates for x 1 , x 2 , x 3 and x 4 . From Fig. 16, we note the noisy determination of the imaginary part of M. As a consequence, we find our estimates of x 1 and x 3 to be noisy as well. On the contrary, we were able to determine x 2 and x 4 reasonably well. In addition to z max dependence, we also studied whether our determination of the moments is affected by the order of truncation used in Eq. 29. We observe no significant variations with truncation. For comparison, the NNPDF and CT18 values of these moments are shown by the horizontal lines. Further, when we fix the values of x 1 and x 3 from NNPDF to reduce the number of fit parameters, we find the estimates for x 2 to be slightly elevated in value and in the direction away from NNPDF,CT18 value. To a small extent, this is seen in x 4 as well. Thus, the observed difference between our lattice result and the NNPDF, CT18 results could be attributed to this tendency for our lattice values of x 2 , x 4 to be slightly higher than the corresponding phenomenological values.
We repeated similar analysis for the helicity matrix element, Γ = γ 5 γ z . In this case, the Wilson coefficients c n (µz) are the same as in the case of unpolarized case with Γ = γ z . Since we are setting the value of the matrix elements at z = 0 to be 1 through the ratio, we only obtain the values of x n / x 0 in the expansion Eq. 29, with x 0 = g A . In Fig. 18, we compare the results corresponding to the NNPDF11pol and JAM17 with the lattice result for the ratio. As in the case of the unpolarized matrix element, we also test the dependence of this comparison on the truncation order n max . The sensitivity to higher moments is a bit more than that for the unpolarized case, and we find convergence at only n max = 6 at z < 0.3 fm. Surprisingly, the global fit expectation agrees quite well with our lattice result even though there is a little tension in the imaginary parts. As explained above, we also obtained the best fit values of x 1 /g A , x 2 /g A , x 3 /g A and x 4 /g A that describes our lattice data via Eq. 29 truncated at most by 4 th order. In Fig. 19, we show the results as a function of the largest z used in the fits, z max . Like the unpolarized PDF case, x 1 /g A is noisy, but seems agree with the global fit results. The more precisely determined value of x 2 /g A is quite robust to various ways of fitting the data and agrees nicely with the global fit values.

IX. SUMMARY AND CONCLUSIONS
In this paper we studied isovector unpolarized and helicity PDFs of proton using the LaMET approach. The lattice calculations have been performed for an unphysically large pion mass of 310 MeV. On the other hand, our lattice study was carried out using lattice spacing a = 0.042 fm, which is the smallest lattice spacing used in such studies. We argued that such small lattice spacing is essential for the validity of 1-loop perturbative matching between PDF and qPDF, which is a key ingredient of LaMET.
Extracting the nucleon matrix elements for such large momenta and small lattice spacing is challenging because of poor signal to noise ratio. To deal with this problem we performed a detailed study of the nucleon two-point function with momentum smeared source and sink, as well as with momentum smeared source and point sink to better control the excited state contributions. We demonstrated that the ground state can be reliably isolated up to the highest momenta used in this study. Furthermore, for the Euclidean time separations used that are relevant for our lattice analysis the two-point function is very well described by the ground state and and an 'effective' excited state contribution, with the energy that is larger than the true excited state energy. Therefore, we argued that the two-state Ansätze is sufficient to describe the dependence of the 3-point function on the source-sink separation and on the operator insertion time. We showed that the qPDF matrix elements can be extracted in this way, and the results do no depend on the choices of the fit interval used in our study, demonstrating the robustness of our analysis procedure.
After non-perturbative RI-MOM renormalizations we compared the lattice calculations of the spatial, z, dependence of qPDFs with that from the phenomenological PDFs, obtained from the global pQCD-based analyses of pertinent experimental data performed by different collaborations. Working in z-space allowed us to test the LaMET approach. The comparisons showed that there is a rough agreement between the lattice results and the results of global analysis, but only at quite small distances. Even for the very small lattice spacing used in this study, there was not enough data points to constrain the x-dependence of the PDFs. Instead, to translate our z-space comparisons to x-dependence, we introduced a new ratio-based renormalization scheme for the Ioffe-time distributions. Using our lattice calculations for Ioffe-time distributions, renormalized via this new ratio-based scheme, we determined the first moments of the isovector unpolarized and helicity PDFs of proton, and compared these moments with that from the corresponding phenomenological PDFs.         results for R f it 1 . As one can see from the figure R f it 1 can describe the data well for all values of t. In Fig. 26 we show the same analysis but for n z = 5.
As discussed in the main text we performed R f it 1 using single value of source sink separation. The results are shown in Fig. 27 for the real part of the matrix element. As one can see from the figure the results obtained from this fit for t = 16, 18 and 20 agree within errors. We performed fits using the form f it 1 with τ > τ min and taking the value of E 1 from the 2-point function fit with t > t min . The results are shown in Fig. 28. We see no significant dependence on τ min and t min .
Another way to obtain the matrix element is to use the summation method. The summation method is illustrated in Fig. 29 for n z = 4. The results obtained from the summation method agree with those from R f it 1 but have much larger errors. The statistical errors of the n z = 5 data are too large to use the summation method. Furthermore, we could also reduce the error in the summation method by dividing by the matrix element at z = 0 as can be seen in Fig. 30 Similar analysis of the ratio if the three point function to two point function was carries out for longitudinally polarized qPDF operator. The results are summarized in Figs. 31, 32 33. To take the advatage of correlation between different z and cancel the field renormalization factor, we divided the bare matrix elements by the matrix     Re(h(z, P z , 5 z ))