Unpolarized isovector quark distribution function from lattice QCD: A systematic analysis of renormalization and matching

We present a detailed lattice QCD study of the unpolarized isovector quark parton distribution function (PDF) using a large-momentum effective theory framework. We choose a quasi-PDF defined by a spatial correlator which is free from mixing with other operators of the same dimension. In the lattice simulation, we use a Gaussian-momentum-smeared source at M π ¼ 356 MeV and P z ∈ f 1 . 8 ; 2 . 3 g GeV. To control the systematics associated with the excited states, we explore five different source-sink separations. The nonperturbative renormalization is conducted in a regularization-independent momentum subtraction scheme, and the matching between the renormalized quasi-PDF and MS PDF is calculated based on perturbative QCD up to one-loop order. Systematic errors due to renormalization and perturbative matching are also analyzed in detail. Our results for light-cone PDF are in reasonable agreement with the latest phenomenological analysis.


I. INTRODUCTION
Parton distribution functions (PDFs) of nucleons are not only important quantities characterizing the internal hadron structures but are also key ingredients to make predictions for high-energy scattering processes [1][2][3]. Thus calculating PDFs from first principles has been a holy grail in nuclear and particle physics. Since PDFs are embedded with the low-energy quark and gluon degrees of freedom in the hadron, they involve infrared (IR) dynamics of strong interactions and can only be determined by nonperturbative methods such as lattice QCD.
Within QCD factorization [4], the quark PDF is defined as where jPi denotes the nucleon state with momentum P μ ¼ ðP t ; 0; 0; P z Þ. x is the quark momentum fraction, μ is the renormalization scale in the MS scheme. ξ AE ¼ ðt AE zÞ= ffiffi ffi 2 p are the light-cone coordinates. The lightlike Wilson line is introduced to maintain the gauge invariance: PDFs are defined with light-cone coordinates, but the lattice simulation can only be conducted in Euclidean space with no proper treatment for light-cone quantities, which involves real time. Thus simulating PDFs on a Euclidean lattice is an extremely difficult task. Early studies based on operator product expansion were only able to derive the lowest few moments of the PDFs [5][6][7][8].
Recently, a novel approach that allows to directly access the x dependence of PDFs from lattice QCD was proposed in Ref. [9], now formulated as large-momentum effective theory (LaMET) [10]. Within this framework, one can extract PDFs-as well as other light-cone quantities-from the correlations of certain static operators in a nucleon state. On the one hand, the static correlations, often referred to as quasiobservables, can be directly calculated on a Euclidean lattice and depend dynamically on the nucleon momentum. On the other hand, at large momentum, the quasiobservables can be factorized into the parton observable and a perturbative matching coefficient, up to corrections suppressed by powers of the large nucleon momentum. Equating the results from the two sides provides a straightforward way to determine the light-cone PDFs.
To calculate the quark PDF in LaMET, one starts with a "quasi-PDF" which is defined as an equal-time correlation of quarks along the z direction [9]: In the above, O Γ ðzÞ ¼ψðzÞΓUðz; 0Þψð0Þ with Γ ¼ γ z or Γ ¼ γ t , and the spacelike Wilson line is For finite but large momentum P z ,qðx; P z Þ has support in −∞ < x < ∞. Unlike the light-cone PDF that is boost invariant, the quasi-PDF has a nontrivial dependence on the nucleon momentum P z . After renormalizing the quasi-PDF in a scheme such as the regularization-independent momentum subtraction (RI/MOM) scheme, one can match the renormalized quasi-PDF to the MS PDF through the factorization theorem [9][10][11][12][13][14]: qðx; P z ; p R z ; μ R Þ ¼ where p R z and μ R are introduced in the RI/MOM scheme. p R z is the momentum of the involved parton and μ R is renormalization scale. r ¼ μ 2 R =ðp R z Þ 2 , C is the perturbative matching coefficient, and OðM 2 =P 2 z ; Λ 2 QCD =x 2 P 2 z Þ denotes nucleon mass and higher-twist contributions suppressed by powers of the large nucleon momentum. The flavor indices in q;q, and C are implied. When −1 < y < 0, the distributions refer to the antiquark distributions.
Since the proposal of LaMET, remarkable progress has been made in both theoretical aspect and lattice calculations. It should be pointed out that these developments are achieved in an interactive way. The LaMET was first used to calculate the proton isovector quark distribution f u−d [15][16][17][18][19][20], including the unpolarized, polarized, and transversity cases, and subsequently to the meson distribution amplitudes [21,22]. The first lattice studies used the matching coefficients at one-loop order in a transversemomentum cutoff scheme [23][24][25]. However, as was found in Ref. [23], the original quasi-PDF suffers from an ultraviolet (UV) linear divergence which might pose a severe problem for the renormalization of its lattice matrix elements [26][27][28]. Then much attention has been paid to the renormalization property [29][30][31][32][33][34][35][36][37], and finally the multiplicative renormalizability of quasi-PDF in coordinate space in the continuum was proven to all orders in the strong coupling constant α s [34,35]. This finding has further motivated the lattice analysis of nonperturbative renormalization (NPR) of the quasi-PDF [36,38,39] in the RI/ MOM scheme [40], and the calculation of the matching coefficients between the RI/MOM quasi-PDFs and MS PDFs [11]. Besides the renormalization, the finite nucleon mass corrections were also worked out to all orders of M 2 =P 2 z [17], and higher-twist OðΛ 2 QCD =x 2 P 2 z Þ effects were numerically removed by extrapolating the results at several P z values to infinite momentum [15,17]. Based on these studies, calculations of the isovector quark PDF at physical pion mass have become available [20,[41][42][43]. Potential operator mixing in the lattice renormalization of the quasi-PDF has also been investigated [33,36,38,39], and the mixing pattern classified in Ref. [44]. Ways to reduce the systematic uncertainties from Fourier transforming the spatial correlation at long distance were proposed in Refs. [41,45]. The LaMET was also attempted to study transverse-momentum-dependent distributions [46][47][48][49][50][51][52][53], as well as the gluon PDF [54][55][56][57][58][59].
In addition to LaMET, other interesting approaches have been proposed in recent years to calculate the PDFs from lattice QCD. For example, one can extract the PDFs from a class of "lattice cross sections" [13,14], while a smeared quasi-PDF in the gradient flow method was proposed to sweep the power divergence in the lattice calculation [60,61]. One can also study a pseudodistribution [62], which shows interesting renormalization features [12,[63][64][65]. Moreover there are proposals using current-current correlators to compute the hadronic tensor [66,67], or the higher moments of the PDF, light-cone distribution amplitudes, etc. [67][68][69][70][71][72]. These different approaches are subject to their own systematics, but they can be compared to each other.
It was argued that the power divergent mixing between local moment operators may spoil the renormalization of quasi-PDFs [27,28], however, such a problem dissolves in LaMET since one first needs to take the continuum limit of the quasi-PDF after renormalization on the lattice, and then match it to obtain the x dependence of the PDF. The factorization has been derived rigorously [12,13] in the continuum, and one only needs to focus on the renormalization of the nonlocal spatial correlator only. Thus the renormalization of local moment operators is irrelevant to quasi-PDF. Besides, there are also confusions on the LaMET matching between Minkowskian and Euclidean matrix elements of the quasi-PDF [73], which have been clarified in Refs. [74,75].
Most of the available lattice calculations have used Γ ¼ γ z (except [42,43]) for the unpolarized quasi-PDF, which is now known to mix with the scalar quasi-PDF operator O I at Oða 0 Þ [33,36,44]. This operator mixing introduces an additional systematic uncertainty in nonperturbative renormalization [36,38,39,41], thus limiting the accuracy of the extracted PDF. On the contrary, the Γ ¼ γ t case is free from operator mixing with O I at Oða 0 Þ [33,36,44]. Therefore, it is highly desirable to start from the quasi-PDF with Γ ¼ γ t . This is one main motif of this study.
In this work, we will carry out a lattice calculation of the unpolarized isovector quark distribution from the quasi-PDF with Γ ¼ γ t with the same nonperturbative renormalization procedure as for the Γ ¼ γ z case in Ref. [39]. The calculation is performed using clover fermions on a CLS ensemble of gauge configurations with N f ¼ 2 þ 1 (degenerate up/down, and strange) flavors under an open boundary condition [76] with pion mass M π ¼ 356 MeV and lattice spacing a ¼ 0.086 fm [77]. We will examine the dependence on the nucleon momentum P z and the RI/MOM scales p R z , μ R , as well as on choices of the projection operator for the amputated Green's function in RI/MOM renormalization. Due to large uncertainties, it is hard to see the sea quark asymmetry observed in early studies which were performed without lattice renormalization [15][16][17][18]. In the future we plan to analyze CLS ensembles with a better accuracy and down to physical masses and a < 0.04 fm, both for (m s þ m u þ m d ) fixed to its physical value and for physical m s using flavor SU(3) and SU(2) extrapolations.
The rest of this paper is organized as follows. In Sec. II, we briefly review the procedure of nonperturbative renormalization and matching of the quasi-PDF in the RI/MOM scheme, in particular the explicit one-loop matching coefficient for the Γ ¼ γ t case. In Sec. III, we describe the details of lattice simulation of the hadronic matrix elements as well as its nonperturbative renormalization. Systematic errors in the calculation are also discussed in this section. In Sec. IV, we present our results on the x dependence of the unpolarized isovector quark PDF with the statistical and systematic uncertainties, and the last section contains the summary of our work.

II. NONPERTURBATIVE RENORMALIZATION AND MATCHING
To recover the continuum limit of a quasi-PDF matrix element, nonperturbative renormalization on the lattice is required to deal with linear and logarithmic UV divergences. In this work, we follow the RI/MOM scheme elaborated in Refs. [11,39], and match the result to the MS PDF with the one-loop matching coefficient [11].

A. RI/MOM renormalization on the lattice
The spatial correlator O Γ ðzÞ has been proven to be multiplicatively renormalizable in coordinate space in the continuum [34,35], which enables the renormalization in RI/MOM scheme [40].
For each value of z, the RI/MOM renormalization factor Z is obtained by requiring loop corrections for the matrix element of a quasi-PDF operator vanish in an off-shell quark state at a given momentum: The bare matrix element P s hp; sjO γ t ðzÞjp; si will be calculated on the lattice from the amputated Green's function Λ γ t ðp; zÞ of O γ t ðzÞ, with a projection operator P for the Dirac matrix: Due to the breaking of Lorentz covariance in O Γ ðzÞ, the RI/MOM subtraction depends on two scales μ R and p R z . As a result, the renormalization factor Zðz; p R z ; a −1 ; μ R Þ depends on the lattice spacing as well as on the two RI/MOM scales μ R and p R z .
Based on the symmetry of O Γ ðzÞ on the lattice, the amputated Green's function Λ γ t ðp; zÞ is not only proportional to the tree-level result γ t , but also includes two other independent Lorentz structures: In the aboveF i 's are independent form factors that are invariant under the hypercubic group Hð4Þ. According to Eq. (8) the RI/MOM renormalization factor Z will also depend on the projection operator P. One can choose to single outF t only [11], which we call the minimal projection. This projection has the simplest form but captures all the UV divergence in Λ γ t ðp; zÞ. Optionally, one can choose P ¼ p=ð4p t Þ, which we call the p projection. The renormalization factors Z with the minimal and p projections are defined as The bare nucleon matrix element from a lattice calculation in coordinate spacẽ Hereh R ðz; P z ; p R z ; μ R Þ is the continuum limit of the renormalized matrix element. Consequently, the quasi-PDFq R ðx; P z ; p R z ; μ R Þ in the RI/MOM scheme is obtained through the Fourier transformation ofh R ðz; P z ; p R z ; μ R Þ: In RI/MOM,h R ðz; P z ; p R z ; μ R Þ andq R ðx; P z ; p R z ; μ R Þ are independent of the UV regulator, and the one-step matching between the quasi-PDF and MS PDF can be carried out in the continuum theory with dimensional regularization [11].
The quasi-PDFs will eventually be matched to the same MS PDF, and the two projections with Z mp and Z p should generate the same result. However, the matching coefficient can only be calculated at a fixed loop order, hence remanent dependence on the projection operator is inevitable.
Using the same logic one reaches the conclusion that the quasi-PDF's dependence on the RI/MOM scales μ R and p R z should also be fully canceled by the matching coefficient. Any fixed-order matching calculation will inevitably lead to a residual μ R , p R z , and P z dependence of the final result for the PDF. These dependencies should be carefully studied and included in the systematic uncertainties.

B. One-loop matching for quasi-PDF and PDF
To obtain the matching coefficient between the quasi-PDFq R ðx; P z ; p R z ; μ R Þ and light-cone PDF qðx; μÞ, one can calculate the off-shell quark matrix elements in perturbation theory. In the following, the calculation will be conducted in Landau gauge for both minimal and p projections. See Sec. 1 of the Appendix for the results in a general covariant gauge with a general Lorentz structure.
The lowest-order quark quasi-PDF is At one-loop order, it is Thef i s aref with iϵ giving the prescription to analytically extrapolate ρ from ρ < 1 (Minkowski) to ρ > 1 (Euclidean). Notice that the vector current conservation guarantees that vertex corrections and wave function contributions can be combined into generalized plus functions [11]. These functions are defined with two arbitrary functions hðxÞ and gðxÞ: For the light-cone PDF with the same off-shell IR regulation in Landau gauge, the tree-level contribution is and the one-loop correction in the MS scheme is q ð1Þ ðx;p;μÞ ¼ Tr Here To match the quasi-PDF to the light-cone PDF, one needs to take the on-shell limit (p 2 → 0 or ρ → 0) and the large-momentum limit (p t → p z ) for the bare quasi-PDF q ð1Þ B ðx; ρÞ ¼q ð1Þ ðx; ðp t → p z ; ⃗p ⊥ ; p z Þ; ρ → 0Þ: One can observe that both terms proportional to γ t and γ z in Eq. (14) approach light-cone operators in the large-momentum limit and the combination of them captures the correct collinear behavior. Therefore the bare quasi-PDF in minimal projection is defined to pick up the coefficient of γ t and γ z in Eq. (14): For the light-cone PDF, the coefficient of γ þ in Eq. (21) is used for minimal projection: The bare matching coefficient is then derived as where In RI/MOM, the quasi-PDF is renormalized with an additional counterterm. We find that in the jxj → ∞ limit, onlyf t ðx; ρÞ behaves as 1=jxj. When integrating over x, this term recovers UV divergence in the local limit z ¼ 0. Therefore, it is a natural choice to pick up the γ t term in Eq. (14) as a counterterm:q Here r ¼ μ 2 R =ðp R z Þ 2 , and Finally, the matching coefficient C in the factorization formula given in Eq. (5) is derived as Here the coupling α s ðμÞ is in the standard MS scheme. Note that the antiquark distribution is mapped into the region −1 < y < 0 by setting qðyÞ ¼ −qð−yÞ.
For p projection, one has the bare quasi-PDF q ð1Þ B ðx; ρÞj p ¼ ½f t ðx; ρÞ þf z ðx; ρÞ þf p ðx; ρÞ þ : ð32Þ The light-cone PDF with a similar projection is Under this projection, the matching coefficient for the bare quasi-PDF coincides with Eq. (27): The counterterm can be obtained by calculating with P ¼ p=ð4p t Þ: The corresponding RI/MOM matching coefficient is obtained by replacing mp with p in Eq. (31), and the difference between f 2;p and f 2;mp vanishes in the p R z ¼ 0 limit. The matching coefficient with Γ ¼ γ z is also given in Sec. 2 of the Appendix.

A. Lattice matrix elements
In this section, we give the results of a lattice-QCD calculation using clover valence fermions on the CLS 32 3 × 96 2 þ 1 flavor clover fermion ensemble H102 with lattice spacing a ¼ 0.086 fm, pion mass M π ¼ 356 MeV, and box size L ≈ 2.7 fm (M π L ≈ 4.9) [77]. We use κ l ¼ 0.136865 and C SW ¼ 1.98625 for the valence clover fermion. We apply APE smearing [78] with size ¼ 2.5a twice in the source/sink smearing and also in the quasi-PDF operator O Γ , but not in the fermion propagators.
First of all, we will explore the nonperturbative renormalization in the RI/MOM scheme. Following Ref. [39], we use Landau gauge fixed wall sources (while limiting the source in the time slice range t ∈ ½32; 64 to avoid the boundary effect from the open boundary condition at t ¼ 0), and generate the propagators with the momentum modes listed in Table I. The four digits in brackets correspond to three spatial components and the energy component in units of 2π=L. The z direction can be selected by setting the Wilson link along any of the three spatial directions. Thus these results approximately cover three values of μ R ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi −ðp R Þ 2 p (2.4, 3.2, and 3.9 GeV, corresponding to a 2 μ 2 R ¼ 1.1, 2.0, and 2.8), and p R z ¼ f0; 1; 2; …; g2π=L up to the upper limit p R z < μ R . Note that in deriving these momentum modes we have required the spatial components of a given momentum mode to differ from each other, and adjusted p R t to ensure μ R invariant (within 2%). These choices allow us to explore the dependence on each component of p, but one should be cautious that the results may suffer from sizable discretization errors since the normal constraint P μ a 4 p 4 μ ð P μ a 2 p 2 μ Þ 2 < 0.3 is not respected. These discretization errors will be investigated in the future.
As shown in Eq. (31), the one-loop matching formula primarily depends on the combination r ¼ μ 2 R =ðp R z Þ 2 but is independent of p R t . In Fig. 1, the Z mp and Z p − Z mp at fixed z ∼ 0.5 fm are plotted as a function of p R z =μ R ¼ 1= ffiffi ffi r p . From this figure, one can see that the NPR factors, both real and imaginary parts, only show the dependence on r regardless of the values of p R z or μ R , with p R z =μ R < 0.4. In Fig. 2, we show 1=Z mp and 1=Z p − 1=Z mp as functions of the Wilson link length z, with the same μ R ¼ 3.2 GeV and p R z ¼ 1.4 GeV and two different values of p R t ¼ 0.9 and 2.7 GeV. As shown in this figure, the 1=Z mp and 1=Z p − 1=Z mp with the two different p R t 's are close to each other for all z (the curves with the same color). This is consistent with the one-loop matching formula. At z < 0.3 fm, the real part of 1=Z p − 1=Z mp with the two different p R t 's can be slightly nonzero, but it is still smaller than 1=Z p by 2 orders of magnitude.
In the following we will take p R z to be zero and estimate the systematic uncertainty from the p R z dependence by varying the p R z .
In the calculation of the nucleon matrix element, we use Gaussian momentum smearing [79] for the quark field ψðxÞ → S mom ψðxÞ , with various μ R . At p R z ¼ 0, Z is real and δZ=Z is less than 5%. where k is the desired momentum, U APE j ðxÞ are the APE smeared gauge links in the j direction, and α is a tunable parameter as in traditional Gaussian smearing.
Such a momentum source is designed to increase the overlap with nucleons of the desired boost momentum and we are able to reach higher-boosted momentum for the nucleon states than in the previous work [39]. Although in the exploratory study we varied the Gaussian smearing radius to better overlap with the largest momentum used in the calculation, the field smearing is still centered around zero momentum in momentum space. When we switch to the momentum smearing, the smearing center will be shifted to momentum OðkÞ, which will immediately allow us to reach higher boost momenta with better signal-tonoise ratios in the matrix elements. In this work, we use two values of nucleon boost momenta, P z ¼ n 2π L , with n ∈ f4; 5g, which corresponds to 1.8 and 2.3 GeV.
On the lattice, we calculate the time-independent and nonlocal in the z direction correlators of a nucleon with finite-P z boost Here the state j0; ⃗ Pi represents the ground (nucleon) state with momentum ⃗ P ¼ f0; 0; P z g. Γ ¼ γ t is used for the unpolarized parton distribution.
As the nucleon boost momentum increases, one anticipates that excited-state contributions are more severe; therefore, a careful study of the excited-state contamination is necessary. To do so, we calculate the nucleon matrix elementh lat at five source-sink separations t seq ∈ f7; 8; 9; 10; 11g × 0.086 fm, with f4; 4; 8; 8; 16g measurements on each of 2005 gauge configurations respectively in FIG. 3. The real (left panels) and imaginary (right panels) parts of the isovector nucleon matrix elements for unpolarized PDFs as functions of z at different momenta, with P z ¼ 8π L ¼ 1.8 GeV (top panels) and 10π L ¼ 2.3 GeV (bottom panels) respectively. The RI/MOM renormalization factors with fμ; p R z g ¼ f3.2; 0g GeV and the normalization at z ¼ 0 are applied on the bare matrix elements to improve the visibility at large z. At a given positive z value, the data are slightly offset to show the ground-state matrix element from the fits using different ranges; from left to right they are: t seq ∈ ½7; 9, [7,10], [7,11], [8,11], and [9,11]. Different analyses are consistent within statistical errors, while the fits with separation 7 and 8 have smaller uncertainties compared to other cases. the P z ¼ 1.8 GeV case, and of 2000 configurations in the P z ¼ 2.3 GeV case. We use a multigrid algorithm [80,81] with the CHROMA software package [82] to speed up the inversion of the quark propagator. Following Ref. [83], each three-point (3pt) correlator C ð3ptÞ Γ ðt; t seq Þ can be decomposed as (assuming the source is at t ¼ 0) where jni with n > 0 represents the excited states. The operator is inserted at time t, and the nucleon state is annihilated at the sink time t seq (which is also the sourcesink separation). The spectrum weights A 0;1 and energies E 0;1 in Eq. (39) can be obtained from the two-point (2pt) correlator: Eventually we apply the joint fit with the 3pt functions at several t seq and 2pt function using the following form [83]: with ΔE ¼ E 1 − E 0 . C 0;1;2;3 and E 0;1 are free parameters. We limit the range of t as t ∈ ½1; t seq − 1 for the 3pt/2pt ratio and t ∈ ½7; 11 for the 2pt to make the χ 2 =d:o:f: of the fit to be Oð1Þ. Using the ratio of 3pt/2pt instead of the 3pt function itself can improve the stability of the fit, especially when C 2pt ðtÞ with t < 7 is included in the fit.
FIG. 4. The real (left panels) and imaginary (right panels) parts of the renormalized isovector nucleon matrix elements for unpolarized PDFs with P z ; z ¼ 8π=L; 12 (top panels) and 10π=L; 10 (bottom panels), which correspond to zP z ∼ 9.5. The data points and the band predicted by the fit using t seq ∈ ½7; 11 agree with each other well.
In Fig. 3, we show the ground-state nucleon matrix elementsh lat ðz; P z ; γ t Þ obtained from five fits: using the separations t seq ∈ ½7; 9, [7, 10], [7,11], [8,11], and [9,11]. (The data points correspond to the same z but are shifted horizontally to enhance the visibility.) The data are further normalized by multiplying the renormalization factor with fμ R ; p R z g ¼ f3.2; 0g GeV and the real part normalized to 1 at z ¼ 0. From this figure, one can see that there is no clear signal for excited-state contributions in any of these analyses. If the data with smallest two separations are dropped, uncertainties are getting much larger. In the fit, we keep the C 3 term to make a moderate estimate of the uncertainty even when this term is not statistically significant.
For a comparison between data and the fit, we show our results at large z like ðP z ;zÞ ¼ ð8π=L;12aÞ and ð10π=L; 10aÞ with t seq ∈ ½7; 11 in Fig. 4. In these spatial separations, the real part of the matrix element seems to be negative. The ground-state contribution obtained from the fit is shown as the black band. As one can see, most data can be well described in the fit and thereby we use the two-state fits and the interval t seq ∈ ½7; 11 to obtain the results in the rest of this paper.
The renormalized quasi-PDF matrix elements with two values of P z are plotted in Fig. 5, as a function of zP z for p R z ¼ 0 and μ R ¼ 3.2 GeV. The results with different P z are consistent with each other within statistical uncertainties. This indicates that power corrections due to highertwist effects might not be sizable.

B. Systematic uncertainties
In this section, we will consider four systematic uncertainties, from Fourier transformation (FT), unphysical scales p R z and μ R , projection used in the RI/MOM scheme, and inversion of the matching coefficient.
In the following, we explain the details to include these systematic uncertainties.
(1) Fourier transformation. As shown in Fig. 5, theh R ðzÞ with P z ¼ 2.3 GeV is consistent with zero when z > 12a. Thus in the standard matching from quasi-PDF to PDF, it is reasonable to truncate the results at z ¼ 12a. With this spirit, the quasi-PDF and matched PDF using the standard FT are shown in the upper panel of Fig. 6, from which one can see the matched PDF shows an oscillatory behavior. A "derivative" method was proposed in Ref. [41] to cure this oscillatory behavior. To be concrete, one takes the derivative of the renormalized nucleon matrix elements ∂ zhR ðzÞ, whose Fourier transform differs from the original matrix element in a known way: From the above equation, one can see that whenh R ðzÞ goes to zero as jzj → ∞, performing an integration by part and neglecting the boundary term will lead to the ordinary Fourier transformation. In implementing the derivative method in Eq. (41), one can perform the integration over z in a discretized way, or interpolate theh R first. We have checked that the two methods give consistent results. With the same truncation as in the ordinary Fourier transformation, the result is shown in the lower panel of Fig. 6 and apparently the oscillatory behavior is less severe. Besides, results obtained using the derivative method are consistent with the standard FT method in most of the kinematics region except at small x. This is anticipated as the two methods only differ in the large-z region where we have made the truncation. We show the difference as the dotdashed-blue line in Fig. 7, together with the error from varying the truncation from z ¼ 10a to 14a (dottedgreen line).
Extending the lattice calculations restricted to a finite number of determinations of matrix elements to recapture all information on the PDF is a sophisticated problem. Thus it is necessary to point out that the derivative method we adopted in Eq. (41) can be viewed as a model-dependent treatment and might introduce some discretization uncertainties. This can possibly be overcome by using reconstruction methods; see Ref. [84] for detailed discussions.
FIG. 5. The renormalized quasi-PDF matrix elements with P z ¼ 1.8 and 2.3 GeV, using the minimal projection with p R z ¼ 0 and μ R ¼ 3.2 GeV, as function of zP z .
(2) Unphysical scales p R z and μ R . There are two unphysical scales p R z and μ R introduced in RI/MOM. In principle, when matching the quasi-PDF matrix element onto a lightcone PDF, the dependence on these two scales in the matrix element should exactly cancel with that in the matching kernel. However, since the quasi-PDF matrix element is nonperturbatively renormalized on the lattice, while the matching coefficient is calculated at one-loop order in perturbation theory, there will be residual dependence on these two scales after the perturbative matching. To estimate the residual p R z and μ R dependence, we choose p R z ¼ 0 GeV and μ R ¼ 3.2 GeV as the central value, and vary p R z from −1.4 to 1.4 GeV (dashed-red line) and μ R from 2.4 to 3.9 GeV (thick-solid-orange line). The difference between these matched PDFs is treated as the systematics of the residual dependence on unphysical scales, in Fig. 7. As shown in the figure, the systematic uncertainty due to the μ R dependence is small compared to the other sources, but the residual p R z dependence could be sizable.
(3) Dependencies on the projection. There are two projections discussed in this work: the minimum and "p" projections. With p R z ¼ 0, the projection dependence in both the NPR factor are less than 5% for all the z, and vanishes in the one-loop perturbative matching. Thus one can expect the difference due to the projections to also be small, as depicted in the long-dashed-magenta lines in Fig. 7.
(4) Inversion of matching. To extract the PDF from the quasi-PDF, one needs to invert the factorization formula Eq. (5). This could be done by changing the sign of α s in C, and convoluting the new matching coefficient withq. More explicitly, we have where C 0 ¼ Cðα s → −α s Þ. We estimate the error due to inverting the factorization formula by starting from the PDF from a global analysis [85], applying Eq. (5)  the systematic error coming from the inversion and higher-order corrections. This is shown in Fig. 8, from which one can find that the error only becomes sizable when jxj is small. This is expected because the relevant momentum scale is xP z such that higher-order corrections become large at small x.
There are more sophisticated methods to invert the factorization, such as using a recursion procedure. However, as we can see in Fig. 7, the systematic error caused by the matching procedure (thick-dashed-cyan line) is also smaller than those from the first two sources in most regions.
As shown in Fig. 7, one can find that the dominant uncertainties arise from the p R z dependence and from FT. With p R z ¼ 0, the uncertainties from different projections, and inversion, and the matching are typically less than 10% except in the region with very small jxj. The uncertainty from the μ R dependence is even smaller. With p R z ¼ 1.4 GeV, uncertainties from these three sources are getting larger in magnitude, but still smaller than those from the two major sources.

IV. FINAL RESULTS FOR PDF
With the derivative method of FT and the matching using p R z ¼ 0 GeV and μ R ¼ 3.2 GeV, we show the dependence on the nucleon boosted momentum in Fig. 9 with the statistical uncertainties. They are consistent with each other as we can expect from the consistency of the quasi-PDF matrix element results in Fig. 5.
Finally, we show our results for the PDF and a comparison with global analysis [85][86][87] in Fig. 10. As can be FIG. 8. Effects of the inversion matching formula using minimal projection. The solid-black, dotted-red, dotted-blue, and dot-dashed-green lines represent the CT14nnlo PDF, applying inverse matching from the CT14nnlo PDF [85] to the quasi-PDF, applying matching again to get back to the PDF, the difference between the PDF with iterative matching and the original CT14nnlo PDF. The upper (lower) figure corresponds to p R z ¼ 0 GeV (1.4 GeV). These plots show that the method we used to invert the matching formula is less reliable for small jxj. The difference shown by the dot-dashed-green curve is taken into account in our systematic error.
FIG. 9. Nucleon boost momentum dependence of the matched unpolarized isovector PDFs: the dotted-green and solid-blue lines correspond to the nucleon momentum P z being 1.8 and 2.3 GeV, respectively. The cutoff of Fourier transformation is chosen to be zP z ∼ 12 (z ¼ 15a for P z ¼ 1.8 GeV and z ¼ 12a for P z ¼ 2.3 GeV).
FIG. 10. Results for PDF at μ ¼ 2 GeV calculated from the RI/MOM quasi-PDF at nucleon momentum P z ¼ 2.3 GeV compared with CT14nnlo (90CL) [85], NNPDF3.1 (68CL) [86], and MMHT2014 (68CL) [87]. Our results agree with the global analysis within uncertainties. seen from the plot, our results show a reasonable agreement in the large-x region, but in the small-x region there exists a notable difference mainly due to the systematic uncertainties from the FT truncation method and also the p R z dependence.

V. SUMMARY
In this paper, we have studied the quasi-PDF defined with γ t which is free from mixing at Oða 0 Þ. We have used M π ¼ 356 MeV lattice data to demonstrate the matching procedure and show that the excited-state contamination is well under control. The one-loop matching coefficient has been calculated and we have discussed the sources of systematic errors as well as the choice of the projection in detail.
We have found that the systematic uncertainties from the FT truncation method and also the p R z dependence are sizable. But those uncertainties from μ R , inversion of matching, and choice of projection are relatively minor with p R z ¼ 0. At the same time, the significant change from a quasi-PDF to a matched PDF suggests that higher-loop corrections are needed, as exhibited in Fig. 6.
Controlling systematic uncertainty from the excited state is very challenging since the relative uncertainty grows very fast when either source-sink separation t seq or nucleon momentum P z becomes large. The two-state fit with smaller separation provides a possibility to obtain a precise result in small t seq < 1 fm region, while for an accurate measurement at large separation using very high statistics, estimating the systematic uncertainty of such a fit is still needed.
Besides the uncertainties that we have studied, in the future we plan to investigate other systematics such as lattice discretization and finite volume effects [88] as well as higher-twist contributions that affect the small-x result. The latter can be improved with larger nucleon momentum and estimated by extrapolating to infinite nucleon momentum.
Our final result for light-cone PDF agrees with the global analysis in the large-x region, which gives an encouraging signal that LaMET may allow us to precisely access parton physics in the future.