A Complete Matching for Quasi-distribution Functions in Large Momentum Effective Theory

We complete the procedure of extracting parton distribution functions (PDFs) using large momentum effective theory (LaMET) at leading power accuracy in the hadron momentum. We derive a general factorization formula for the quasi PDFs in the presence of mixing, and give the corresponding hard matching kernel at $\mathcal O(\alpha_s)$, both for the unpolarized and for the polarized quark and gluon quasi-PDFs. Our calculation is performed in a regularization-independent momentum subtraction scheme. The results allow us to match the nonperturbatively renormalized quasi-PDFs to normal PDFs in the presence of mixing, and therefore can be used to extract flavor-singlet quark PDFs as well as gluon PDFs from lattice simulations.


I. INTRODUCTION
Understanding the internal structure of hadrons from quarks and gluons -the fundamental degrees of freedom of QCD Lagrangian -has been a key goal in hadron physics. However, this is profoundly difficult because it requires solving QCD at large distance scales and thus at strong coupling. In high energy collisions, the hadron and/or the probe moves nearly at the speed of light, the hadron structure greatly simplifies and can be characterized by certain parton observables such as the parton distribution functions (PDFs), distribution amplitudes (DAs) etc. The parton observables are defined as the expectation value of lightcone correlations in the hadron state and therefore can not be readily computed on a Euclidean lattice. Currently, the most widely used approach to determine them is to assume a smoothly parametrized form and fit the unknown parameters to a large variety of experimental data (for a recent review, see e.g. Ref. [1]). Lattice efforts on determining the parton observables have been mainly focused on the computation of their moments, which are matrix elements of local operators. The parton observables can be reconstructed in principle if all their moments are known. However, to date only the first few moments can be calculated in lattice QCD [2][3][4][5] due to power divergent mixing between different moments operators and increasing stochastic noise for high moments operators.
In the past few years, a breakthrough has been made to circumvent the above difficulty, which has now been formulated as large momentum effective theory (LaMET) [6,7]. According to LaMET, a parton observable, instead of its moments, can be directly accessed from lattice QCD using the following procedure: 1) Construct an appropriate static-operator matrix element (quasi-observable) that approaches the parton observable in the infinite momentum limit of the external hadron. The quasi-observable constructed in this way is usually hadron-momentum-dependent but time-independent, and thus can be readily computed on the lattice. 2) Calculate the quasi-observable on the lattice and renormalize it nonperturbatively in an appropriate scheme. 3) Match the renormalized quasi-observable to the parton observable through a factorization formula accurate up to power corrections that are suppressed by the hadron momentum. The existence of such a factorization is ensured by construction; for a proof in the case of isovector quark distribution, see Refs. [8][9][10].
So far the lattice calculations of PDFs have been focused on the isovector quark PDFs only, which do not involve mixing with gluon PDFs and therefore are the easiest to calculate. In the past few years, there has been increasing interest in calculating flavor-singlet quark PDFs and gluon PDFs from lattice QCD. Such calculations are possible only if the renormalization and mixing pattern of gluon quasi-PDFs are fully understood. In Ref. [53], we performed a systematic study of the renormalization property of the gluon quasi-PDF operator, and showed that with an appropriate choice it is multiplicatively renormalizable. We have identified four independent gluon quasi-PDF operators that have an easy implementation on the lattice. Also, a general factorization formula for the gluon as well as the quark quasi-PDF in the presence of mixing has been conjectured.
In this paper, we provide all necessary inputs for extracting both the flavor-singlet quark PDF and the gluon PDF from lattice QCD, thereby completing the procedure of calculating PDFs using LaMET at leading power accuracy in the hadron momentum. We explain how to nonperturbatively renormalize the quark and gluon quasi-PDFs, and derive a general factorization formula for the renormalized quasi-PDFs in the presence of mixing, following the operator product expansion (OPE) method in Refs. [9,10]. We then present the complete one-loop results for the hard matching kernels that appear in the factorization of quasi-PDFs.
The rest of the paper is organized as follows: In Sec. II, we briefly review the renormalization and factorization of quark and gluon quasi-PDFs. In Sec. III, we present our one-loop calculation of the hard matching kernel connecting the RI/MOM renormalized quasi-PDFs to the PDFs in MS scheme, with a particular focus on the unpolarized case. Sec. IV is devoted to the polarized case. We then conclude in Sec. V and give some computational details in the Appendix.

II. RENORMALIZATION AND FACTORIZATION OF QUARK AND GLUON QUASI-PDFS
In this section, we give a brief review of the renormalization and factorization of quark and gluon quasi-PDFs in LaMET.

A. Quasi-PDFs in LaMET
In high-energy collisions, the PDFs are defined as the hadron matrix elements of quark and gluon nonlocal correlators along the lightcone. For example, the unpolarized quark distribution is defined as for a given flavor i, where x = k + /P + is the longitudinal momentum fraction carried by the quark of flavor i. µ is the renormalization scale in the MS scheme, P µ = (P 0 , 0, 0, P z ) is the hadron momentum, ξ ± = (t ± z)/ √ 2 are the lightcone coordinates, and is the Wilson line inserted to maintain the gauge invariance of the nonlocal correlator. The A + = A + a t a with t a being the generators in the fundamental representation of color SU (3) group.
Analogously, the unpolarized gluon distribution can be defined as [79] where F µν a = ∂ µ A ν a − ∂ ν A µ a − gf abc A µ b A ν c is the gluon field strength, and i runs over the transverse indices. The above Wilson line W takes a similar form as the quark case, but is defined in the adjoint representation.
The quark and gluon PDFs defined above can not be directly computed on the lattice due to their real-time dependence. However, according to LaMET, they can be extracted from lattice calculations of appropriately constructed quasi-PDFs via a factorization procedure. For the unpolarized quark PDF, a well-suited quasi-PDF candidate is given byf where z is a spatial direction, Γ = {γ z , γ t } is a Dirac matrix with the corresponding normalization factor N = {1, P z /P t }, respectively. As shown in Ref. [29], the renormalization of the quark quasi-PDF defined above is of a multiplicative form so that the matrix elements at different z do not mix with each other. In addition, the choice with Γ = γ t has the advantage of avoiding mixing with the scalar PDF when a non-chiral lattice fermion is used [19,33]. We will focus on this choice in the rest of the paper.
In comparison with the quark case, what is the most appropriate operator to define the gluon quasi-PDF is less obvious. In principle, one can use with µ, ν = {t, z} and α running either over all Lorentz indices or only over transverse indices. However, such a choice could in principle mix with other relevant operators under renormalization. Using the auxiliary field approach [80], we have explicitly shown [53] that different components of O µν indeed renormalize differently, which complicates the construction of appropriate gluon quasi-PDFs. A brief review of the formalism used in Refs. [53] and [29] will be given in the forthcoming subsections. Nevertheless, we have identified four gluon operators [53] that are multiplicatively renormalizable and therefore are suitable for defining the gluon quasi-PDF. These operators are where a summation over transverse (all) components is implied for i (µ). The corresponding gluon quasi-PDF is then defined asf The normalization factors are chosen by so that all partonic gluon PDFs at tree-level arẽ f (n,0) with the hadron state H being replaced by a gluon state. Note that in the above result (also in the sections below unless stated otherwise) we have ignored the contributions from the crossed diagrams, which correspond to interchanging the contraction between the two external gluons and gluon fields from the operators O (n) g . These crossed diagrams contribute to x < 0 and can be easily obtained fromf All above gluon quasi-PDF operators are defined in terms of an adjoint gauge link. Alternatively, these operators can also be parametrized using gauge links in the fundamental representation U (z 2 , z 1 ) [80][81][82][83][84][85]. Taking the operator O (3) g as an example, one could use Here F µν = F a µν t a and t a is the generator in the fundamental representation with tr[t a t b ] = 1/2δ ab . We wish to stress that the Eq. (5) facilitates the renormalization study of gluon quasi-PDFs, whereas Eq. (10) makes the implementation on the lattice simpler. In the following, we will mainly focus on the definition Eq. (5), as has been done in Ref. [53], but the results also apply to Eq. (10).
In the forthcoming subsections, we briefly review the renormalization of quasi-PDFs in the auxiliary field approach, following our earlier work in Refs. [29,53]. Other studies have been available using a similar formalism [20] or using the Feynman diagrammatic approach [30,54].

B. Auxiliary Field Approach
In the auxiliary field approach [80], one introduces an auxiliary "heavy quark" field into the QCD interaction such that the Wilson line can be reinterpreted as a two-point function of the auxiliary field. For the quark/gluon quasi-PDF, this auxiliary field is chosen to be in the fundamental/adjoint representation of color SU (3) group, respectively. Similar with the ordinary heavy quark, the auxiliary "heavy quark" has trivial spin degrees of freedom. An advantage of this approach is to convert the study of renormalization of nonlocal operators into the analysis of two local operators. In the following we present, as an example, the auxiliary Lagrangian that can be used to study quark quasi-PDFs, while for gluon quasi-PDFs the procedure is completely analogous.
The effective Lagrangian with an auxiliary fundamental "heavy quark" field (denoted as Q) can be written as where D µ = ∂ µ + igt a A a,µ is the covariant derivative in the fundamental representation. The unit vector n µ is chosen as n µ = (0, 0, 0, 1). As shown in Ref. [29], the two-point function of the auxiliary "heavy quark" is evaluated as The above equation holds up to a determinant det(in · D) which can be absorbed into the normalization of the generating functional [86]. The propagator S Q (x, y) in the above is the Green function of n · D operator with An solution can be derived with an appropriate boundary condition. One should notice that Eq. (14) is nothing but a spacelike Wilson line along the z direction. One can always restrict oneself to x z > y z , without loss of generality.
C. Renormalization of Quasi-PDFs in Auxiliary Field Approach

Quark Quasi-PDFs
From the discussions above, one can see that the Wilson line W (z 2 , z 1 ) appearing in the quark quasi-PDFs can be replaced by the product of two auxiliary "heavy quark" fields Q(z 2 )Q(z 1 ). The quark bilocal operator then reduces to the product of two local composite operators Since the "heavy quark" has trivial spin degrees of freedom, one can also move the Dirac matrix Γ into j(z 1 ). In dimensional regularization (DR), the local operatorsj(z 2 ), j(z 1 ) are "heavy-to-light" like and are multiplicatively renormalized:j with (D = 4 − 2 ) When the auxiliary field is integrated out, the nonlocal operator renormalizes as [29] O qi,R (z 2 , z 1 ) = Z −1 j Z −1 jq i (z 2 )ΓW (z 2 , z 1 )q i (z 1 ).
In lattice regularization, when going beyond leading-order perturbation theory, the self-energy of the heavy quark generates a linear divergence that does not show up in DR. Such a linear divergence can be absorbed into an effective mass counterterm, where δm ∼ O(1/a) with a being the lattice spacing [87]. As shown in Ref. [29], apart from the structures given in the Lagrangian Eq. (11), this is the only possible renormalizable counterterm allowed by the symmetry of the theory. Moreover, Becchi-Rouet-Stora-Tyutin (BRST) invariance requires a dependnce of δm on the signature of n in Eq. (11) [80], For a spacelike n µ , δm = iδm is imaginary. Including the effective mass term Eq. (21) into the Lagrangian and integrating out the auxiliary "heavy quark", we obtain the following renormalization for the nonlocal quark bilinear operator [29] O qi,R (z 2 ,

Gluon Quasi-PDFs
For the nonlocal gluon quasi-PDF operators, the desired auxiliary Lagrangian has exactly the same form as that for the quark, except that now the auxiliary "heavy quark" and the covariant derivative are defined in the adjoint representation. To distinguish from the fundamental auxiliary field used in the previous subsection, we denote the adjoint field as Q below.
With the auxiliary Q, one can decompose the nonlocal gluon operator in Eq. (6) into the product of two local composite operators. For example, the O An explicit one-loop calculation in Ref. [81] has indeed verified the above expectation. Since J µz 2 is not independent, it can be ignored in the studies of operator renormalization. In addition, Eqs. (26) and (27) indicate that J zµ 1 and J ti 1 (i = 1, 2) renormalize independently. As a result, the renormalization pattern can be simplified to J zµ The reason that (J ti 1 , J ij 1 ) and J zµ 1 have different renormalizations is due to the Lorentz symmetry breaking in the presence of a four-vector n µ along the z direction.
To extract the UV divergences, in particular the genuine power divergences inherited from the operator J µν 1 , one should introduce a proper UV regulator in a gauge-invariantly manner. In Ref. [53], we worked in DR and kept track of the linear divergences by expanding the results around d = 3, as the linear divergences appear as poles around d = 3 at one-loop.
The one-loop diagrams that give rise to linearly divergent contributions to the operator J µν 1 are shown in Fig. 1, and other diagrams are neglected. We have performed an detailed calculation in coordinate space in Ref. [53], and the result reads F ρσ a n ν n σ − F νσ a n ρ n σ Q a /n 2 + 1 2 (A ρ a n ν − A ν a n ρ )n · ∂Q a /n 2 where µ is the regularization scale and reg. denotes regular terms at both d = 4 and d = 3. Combining the results in Eq. (29), we find the linear divergences cancel. Our results show an identical mixing pattern as in Ref. [81] (note the difference in the normalization of the direction vector). Based on the renormalization analysis above, one can derive useful building blocks for the construction of appropriate gluon quasi-PDFs. To this end, we may use one of the indices in J µν 1 with z or t and let the other indices run either over all Lorentz components or over the transverse components only. It is necessary to point out that the operator J µν 3 only yields contact terms when integrating out the "heavy quark" field, since the equation of motion operator acting on the "heavy quark" propagator yields a δ-function. The nonvanishing contact terms at z 2 = z 1 indicate that an extra renormalization is required when the distance between two local composite operators shrinks to zero. When z 1 = z 2 , the operator J µν 3 is irrelevant and can be ignored. In a cutoff scheme like the lattice regularization, the mass term of the Q could appear beyond leading order in perturbation theory even if it does not exist at leading order. This is indeed what happens here. In perturbation theory, m = δm starts from O(α s ). Such a mass term serves the purpose to absorb power divergences arising from the Wilson line self-energy. Apart from this, there is no other power divergence in the theory. Therefore in a gauge-invariant cutoff scheme, the operator renormalization remains the same as Eq. (28).
With J zi 1,R , J ti 1,R , J zµ 1,R , and their conjugate as the building blocks, four multiplicatively renormalizable unpolarized gluon quasi-PDF operators have been constructed [53], and their explicit form has been given in Sec. II A. To illustrate how the gluon quasi-PDF operators renormalize, let us take as an example. When the auxiliary "heavy quark" field is integrated out, the O g,R (z 2 , z 1 ) operator renormalizes multiplicatively as (δm = iδm) The renormalization of other operators is analogous with different renormalization factors [53]. Actually, the operators O (i) g,R (i = 1, 2, 3, 4) belong to the same universality class [89] and differ only by power corrections in the large momentum limit. Bearing in mind the different renormalizations of O (i) g,R , one may use any combination of them to study the gluon quasi-PDF. A notable example is This operator (minus the trace term) has been used in a recent simulation [71]. Since the renormalizations for O (1) g,R (z 2 , z 1 ) is not multiplicatively renormalizable.

D. Renormalization in RI/MOM Scheme and Implementation on Lattice
From the discussions above, it is clear that the nonlocal operators at different z do not mix under renormalization. This allows us to carry out a nonperturbative renormalization of the quasi-PDF in the following manner: 1) Calculate the endpoint renormalization factors (e.g. Z {11,22} in Eq. (31)) and the Wilson line mass counterterm (δm in Eq. (31)) nonperturbatively. The calculation of the former is rather straightforward, while the latter can be determined by using the static-quark potential for the renormalization of Wilson loops [90]. This has been used in early studies of nucleon PDFs and meson DAs [26,28,65]. 2) Calculate the renormalization factors as a whole for each z. This is analogous to the renormalization of local composite operators, which is usually carried out in the RI/MOM scheme [91] on the lattice. In the RI/MOM scheme, the renormalization of local composite operators is done by demanding that the counterterm cancels all loop contributions to their matrix element between off-shell external states at specific momenta [18,31] (for the application to quark and glue momentum fractions see Ref. [92].) For multiplicatively renormalizable nonlocal correlators such as the quasi-PDFs given above, the renormalization is similar but now one requires calculating the renormalization factors at each z.
The quark and gluon quasi-PDFs can, in general, mix with each other under renormalization. In Ref. [53], we have argued that inserting the gluon quasi-PDF operator into a quark state only yields finite mixing as long as all subdivergences have been renormalized (note the difference from the quark and gluon lightcone PDF operators which mix with each other under renormalization [93,94]). The mixing effect can, in principle, be deferred to be considered at the factorization stage. Here we find that taking into account the mixing at the renormalization stage will help improve the convergence in the implementation of the matching in the RI/MOM scheme. To this end, it suffices to consider the following mixing of quasi-PDFs where O s q (z 1 , ] is the C-even combination of quark operators, Z ij (z) are dimensionless factors, and z compensates for the different mass dimension between the quark and gluon quasi-PDF operators. In the limit z → 0 (taken after combining the entries of the mixing matrix and the operators), the above mixing pattern reduces to the mixing pattern of local operators.
The renormalization factors in the above mixing matrix can be determined using the following renormalization conditions in an offshell gluon and quark state, respectively. P and P ab ij are projection operators that are associated with the quark and gluon matrix elements and define the RI/MOM renormalization factors. µ R and p R z are unphysical scales introduced in the RI/MOM scheme to specify the subtraction point. b, c are color indices and i, j Lorentz indices. In the non-singlet quark PDF case with Γ = γ t [49], the amputated Green's function has the following structure and P was chosen there in such a way that it projects out the coefficient of γ t only, which captures all terms in Λ γ t (p, z) that lead to UV divergences in the local limit. However, in general both the coefficient of γ t and γ z can lead to UV divergences in the local limit. This is the case e.g. in the mixing diagram to be considered below. We will need both coefficients to define the RI/MOM counterterm. As for P ab ij , a simple choice is P ab ij = δ ab g ⊥,ij /(2 − D), where g ⊥,ij denotes the transverse metric tensor and D is the spacetime dimension.
Defining the inverse of the renormalization matrix in Eq. (33) as we then have from Eqs. (33), (34) and (36) Denoting the hadron matrix element of O(z, 0) as h(z, P z , 1/a), i.e., h i (z, P z , 1/a) = P |O i (z, 0)|P , i = q, g, the renormalized hadron matrix elements then read The renormalized quasi-PDF in the RI/MOM scheme can be obtained from the above renormalized matrix elements by a Fourier transform given in Eqs. (4) and (7), respectively. Note that we can take the continuum limit a → 0 in h R since all terms singular in a have been removed by the renormalization procedure. This means that the factorization of the renormalized matrix element can be studied in the continuum, as will be done in the next subsection.

E. Factorization
In Ref. [53], we have given a general factorization formula for the quark and gluon quasi-PDFs in the presence of mixing. In this subsection, we give a detailed derivation of it using the operator product expansion (OPE), along the same line as that used for the isovector quark quasi-PDF [10]. For illustration purposes, we choose Γ = γ t for the quark quasi-PDF and O (4) g for the gluon quasi-PDF. The derivation for other operators follows straightforwardly from what is presented below.
The renormalized quark and gluon nonlocal operator matrix elements can be expanded in terms of gauge-invariant local operator matrix elements to the leading-twist approximation as where we have introduced extra normalization factors so that the two matrix elements have the same mass dimension. For simplicity, we have also denoted all renormalization scales with µ. n t,ρ = (1, 0, 0, 0) and n ρ = (0, 0, 0, 1), C where {· · · } denotes a symmetrization of the enclosed indices. Their matrix elements are related to the moments of quark and gluon PDF, respectively with Owing to the symmetry of the gluon PDF, a g,n does not vanish only for even n.
Let us first considerh qi,R (z, P z , µ). Ignoring all trace terms, we can writẽ where we have introduced the Ioffe-time ν = −z · P = zP z and ν = −ξ · P = −P + ξ − , h qi/g,R denote the coordinate space matrix elements used to define the quark and gluon PDFs at lightlike separation ξ 2 = 0. Defining with u being in the range (−1, 1) [38,95], we then havẽ This is the general factorization of the coordinate space matrix element in the presence of mixing. To convert it to the factorization of quasi-PDFs, we need a Fourier transform of the above relatioñ where we have defined Now let us turn toh qi,R (z, P z , µ). By ignoring all trace terms, one can write as beforẽ Defining we then have the following factorization in coordinate spacẽ The factorization in mometum space reads where we have defined Restoring all renormalization scales, the general factorization of the quark and gluon quasi-PDFs reads where a summation of j over all quark flavors is implied. The factorization for the polarized quasi-PDFs has the same form as Eq. (53) with all unpolarized distributions being replaced by the polarized ones and also different hard coefficients. It is worthwhile to point out that the higher-twist contributions shall behave like 1/[x 2 (1 − x)(P z ) 2 ] instead of 1/(P z ) 2 , as demonstrated in Ref. [55].

III. ONE-LOOP MATCHING FOR UNPOLARIZED QUASI-PDFS IN RI/MOM SCHEME
As shown in the previous section, when the hadron momentum P z is much larger than the hadronic scale, the highet-twist contributions get suppressed (except for very small/large x), the quasi-PDFs can be factorized into the lightcone PDFs with perturbatively calculable hard matching coefficients. In this section, we present the one-loop calculation of the hard matching coefficients for unpolarized quark and gluon quasi-PDFs in the presence of mixing. The polarized case will be discussed in the next section. Our result is obtained in the RI/MOM scheme, which can be used to connect the RI/MOM renormalized quasi-PDFs to the PDFs in MS scheme. Since the matching depends on UV physics only and not on the external state, we can calculate it in quark or gluon external states |q(p) , |g(p) . The infrared (IR) divergences can be regularized using their offshellness.

A. Gluon in Gluon
Let us start with the gluon matrix element of the gluon quasi-PDF operator, which is the most complicated among all calculations. At tree-level one finds: with ρ = −p 2 /p 2 z . As before, we have ignored the crossed terms which can be obtained from {f , f }(x) = −{f , f }(−x). Ignoring such terms has no impact on the extraction of the matching coefficient. The above results lead to the following tree-level matching coefficient: At one-loop level, the partonic quasi-PDF can be written as follows: with n = 1, 2, 3, 4, and the "+" subscript denotes the usual plus-prescription Integrating Eq. (56) over the momentum fraction, one arrives at which corresponds to the matrix element of local operators Before we proceed, a few general remarks on the calculation below are in order.
• The above equations apply to bare operator matrix elements. One can write down similar equations for the renormalized ones. In our calculation of the matching coefficients, the PDF is renormalized in MS scheme while the quasi-PDF is renormalized in the RI/MOM scheme. The renormalized local operator matrix elements in these two schemes differ from each other in general.
• The offshell gluon matrix elements of gauge invariant operators can mix with those of gauge variant operators.
To illustrate this point, it is worthwhile to consider the UV divergence from the offshell gluon matrix element of the local gluon operator F µα (0)F νβ (0): × 1 12 p 2 9g αν g βσ g µρ − 9g αβ g µρ g νσ − g αν g βρ g µσ + g αβ g µσ g νρ + g ασ g βρ g µν − g βµ g νρ +g αρ 9g βµ g νσ − 9g βσ g µν − 2g αν g βµ g ρσ + 2g αβ g µν g ρσ + 1 6 p µ p ν 4g ασ g βρ + 10g αρ g βσ − 7g αβ g ρσ − 1 6 p β p µ (4g ασ g νρ + 10g αρ g νσ − 7g αν g ρσ ) − 1 6 p α p ν 10g βσ g µρ + 4g βρ g µσ − 7g βµ g ρσ + 1 6 p α p β (4g µσ g νρ + 10g µρ g νσ − 7g µν g ρσ ) − 3 4 p µ p ρ g αν g βσ − g αβ g νσ − 3 4 p ν p σ g αρ g βµ − g αβ g µρ + 3 4 p α p ρ g βσ g µν − g βµ g νσ with the crossed contributions being neglected. This leads to the following contributions to the UV divergences inc (n) :c (1,g) = α s C A 12π if a physical projection P ab ij = δ ab g ⊥,ij /(2 − D) is employed. As can be seen from the above equations, the UV divergences might depend on the offshellness of external gluons, which is a sign of the potential mixing with gauge variant operators. It is interesting to note that the UV divergence ofc (3,g) is independent of p 2 . This is because it corresponds to the tz component of the gluon energy momentum tensor for which all gauge variant operators to mix turn out to vanish [96,97]. As we will see below, such a behavior is consistent with the asymptotic behavior at large x of the quasi-PDF defined with O g,R (z, 0), which does not depend on p 2 either. This feature turns out to help achieve a better convergence in the implementation of the matching. Thus, in the following we will focus on O g,R (z, 0), and present the one-loop matching calculation for the gluon quasi-PDF defined with this operator. For completeness and comparison purposes, the results for other definitions are also collected in Appendix A.
• In pure Yang-Mills theory, O g (0, 0) does not renormalize, as shown by the results in Eq. (61) 1 . In QCD, quarks can enter the gluon diagrams relevant for the above calculation, but only through gluon wave function renormalization at one-loop level, and lead to the following contribution toc (3,g) and c (3,g) (the counterpart of c (3,g) for the gluon PDF) after renormalizatioñ c (3,g) This will be needed in the calculation of the matching coefficient below.
Now we present the one-loop results for the partonic quasi-PDF and PDF. The calculation is carried out in Landau gauge, and the steps are similar to those presented in Refs. [22,23]. Given Eqs. (56) and (62), we only present the distribution part, i.e. the first term in Eq. (56). To this end, we need to calculate the one-loop matrix element of O (3) g (z, 0) in an offshell gluon state. The relevant Feynman diagrams are shown in Fig. 2, and the result reads As in the quark case [18,49], the bare quasi-PDF result is obtained by taking the onshell limit ρ → 0 of the above expression except where it has to be kept as an IR regulator xf (3,1) The renormalized lightcone PDF can be calculated analogously, and gives where the result in the first square bracket is the same as the Feynman gauge result. The one-loop matching coefficient is given by the difference in the renormalized quasi-PDF and lightcone PDF where the ln(−p 2 ) dependence in each individual term cancels out in the combination on the r.h.s., and the counterterm in the RI/MOM scheme can be determined from the renormalization condition above as (xf with

B. Quark in Quark
This case has already been considered at one-loop in Ref. [49]. For completeness, we also quote the results here and briefly explain how it was obtained. As we will see below, our definition of the counterterm differs from that defined in Ref. [49] by a finite piece. The relevant Feynman diagrams are shown in Fig. 3.
Owing to the offshellness of the external quark, the one-loop quark quasi-PDF contains two more Dirac structures apart from the tree-level one γ t , and is given by the following projection [49] Tr f (1) q/q,t (x, ρ) where the coefficients of γ t and γ z read in Landau gaugẽ In Ref. [49], a so-called minimal projector for P has been used, which determines the bare quark quasi-PDF as with the following explicit form Note that there is no extra local term likec (3,g) RI/MOM above due to vector current conservation. The renormalized lightcone quark PDF has the following expression The matching coefficient can then be extracted as where again the ln(−p 2 ) dependence cancels out in the combination on the r.h.s., and the counterterm in the RI/MOM scheme is given by Note that the counterterm defined here differs from that given in Ref. [49] by a finite piece. In Ref. [49] the projector used to define the counterterm was chosen differently from that used to define the bare quasi-PDF, and projected out the coefficient of γ t only since only this coefficient contributes to 1/|x| in the asymptotic limit x → ∞. In the present paper, we use the same projector to determine the counterterms in the quark matrix elements of quark and gluon quasi-PDF operators. As can be seen from the next subsection, for the latter we need P to project out both the coefficients of γ t and γ z . Therefore the same projection shall apply to the former. In fact, projecting out both coefficients is more natural, since in the infinite momentum limit both γ t and γ z approach γ + , therefore both may contribute to UV divergences. From a different point of view, we can always rewrite γ z in terms of γ t and / p if the external quark has no transverse momentum. This also implies that taking both the coefficients of γ t and γ z to define the counterterm is more natural.

C. Gluon in Quark
Now we turn to the mixing contributions. Let us first consider the quark matrix element of the gluon quasi-PDF operator, whose one-loop diagram is given in Fig. 4.
To illustrate the kinematic dependence of the mixing terms, it is useful to begin with the one-loop quark matrix element of the local operator F µα (0)F νβ (0) From the above result we obtain As the tz component of the gluon energy momentum tensor, O g,R (0, 0) in general mixes with the same component of the quark contribution where (· · · ) denotes an antisymmetrization of the enclosed indices. The above operator has the same momentum dependence as Eq. (77) when sandwiched in a quark state. This indicates that the mixing matrix element in Fig. 4 has the same momentum dependence as the tree-level quark contribution, which is indeed needed to define an appropriate RI/MOM counterterm. The renormalized mixing contribution from the lightcone gluon PDF has the following form For the quasi-PDF, we follow the decomposition as in the quark case: Tr xf and choose the projector P such that it projects out the coefficients of both γ t and γ z . We therefore have xf (3,1) g/q = xf (3,1) g/q,t + xf leading to xf (3,1) In the limit ρ → 0, we have for the bare quasi-PDF In the limit x → ∞, the above expression behaves asymptotically as If one integrates over the momentum fraction with DR, it is straightforward to see that the above behavior is consistent with the local result in Eq. (77). As before, the matching coefficient can be extracted as where the counterterm in the RI/MOM scheme is determined as (xf

D. Quark in Gluon
Now let us consider the gluon matrix element of the quark quasi-PDF operator, and the one-loop diagram is shown in Fig. 5. We again start with the local matrix element If µ = t and physical polarizations are used for the external gluons, one has the result: which also has the same momentum dependence as the gluon matrix element of O g,R (0, 0).
For the lightcone PDF, the result of the mixing diagram in Fig. 5 reads while for the quasi-PDF one has: Taking ρ → 0 gives the bare quasi-PDF result The matching coefficient is then given by

A. Gluon in Gluon
Now we turn to the polarized case. The calculation can be done in complete analogy with that presented in the previous section. As demonstrated in Ref. [53], to study the polarized gluon PDF we may use the following three operators to define the corresponding quasi-PDF where ⊥,ij is the two-dimensional antisymmetric tensor: with the convention 0123 = 1, and n µ t = (1, 0, 0, 0). The projection operator for the polarized gluon quasi-PDF is chosen as: As before, we decompose the polarized quasi-PDF as Integrating over the momentum fraction one obtains the matrix element of the corresponding local operators: with The local matrix elements have the following UV divergence structure where only the UV divergence of ∆c (3) does not depend on the external momentum. For the same reason as the unpolarized case, we choose ∆O to define the polarized gluon quasi-PDF and present the corresponding one-loop matching kernel below.
The light-cone PDF yields the following real contribution whereas the quasi-PDF gives x∆f (3,1) In the ρ → 0 limit, the above result gets simplified x∆f (3,1) The virtual contribution is the same as the unpolarized case, whereas the real contribution differs in the asymptotic limit x → ∞ as x∆f (3,1) Integrating over x in DR, this gives the UV divergence in Eq. (106) as expected.
The matching kernel can be written using the matching kernel for the unpolarized gluon quasi-PDF as where again the ln(−p 2 ) dependence in each individual term cancels out in the combination on the r.h.s., and the counterterm in the RI/MOM scheme is determined as:

B. Quark in Quark
For completeness, we also give the result for the polarized quark quasi-PDF and PDF defined as following The result for the polarized quark PDF is the same as that for the unpolarized one: For the quasi-PDF, the one-loop result can be decomposed as Tr ∆f where we define P to project out both the coefficients of γ t γ 5 and γ z γ 5 . For the Dirac matrix in Eq. (113), ∆f (1) q/q,t

vanishes, and ∆f
(1) q/q,z is given as: In the limit ρ → 0, it reduces to: The matching coefficient can then be extracted as where the counterterm in the RI/MOM scheme is determined as:

C. Gluon in Quark
The matrix element of the local gluon operator between the polarized quark state reads: where we have used the following identity Projecting onto ∆O (3) g,R (0, 0) gives: which has the same momentum dependence as the quark matrix element of the quark quasi-PDF operator in Eq. (113). The light-cone result for the polarized PDF is given as The corresponding quasi-PDF reads In the limit ρ → 0, we have x∆f (3,1) The matching coefficient can be extracted as x∆C (3,1) g/q (x, r, where the counterterm in the RI/MOM scheme is determined as:

D. Quark in Gluon
In this case, the light-cone result is whereas the quasi PDF result reads In the limit ρ → 0, we have The matching coefficient is then given by with (∆f V. CONCLUSION In this paper, we have studied how to extract the flavor-singlet quark PDF and the gluon PDF from LaMET, both in the unpolarized and in the polarized case. After briefly reviewing the auxiliary "heavy quark" formalism used in our earlier work to prove the multiplicative renormalizability of quark and gluon quasi-PDF operators, we explained how a nonperturbative RI/MOM renormalization can be carried out for the quark and gluon quasi-PDFs on the lattice in the presence of mixing. Using OPE, we also derived the factorization formulas that connect them to the usual quark and gluon PDFs in MS scheme. We then performed a one-loop calculation of the hard matching kernel appearing in the factorization. We found that certain gluon quasi-PDF operators are more favorable than others in the sense that the mixing with gauge variant operators can be avoided. We then focused on these operators and presented the corresponding one-loop matching kernel. Our results can be used to extract the flavor-singlet quark PDFs as well as the gluon PDFs from lattice simulations of the corresponding quasi-PDFs. We therefore completed the procedure of extracting quark and gluon PDFs from LaMET at leading power accuracy in the hadron momentum. It is interesting to note that the matrix elements of those non-favorable gluon quasi-PDF operators have nontrivial momentum dependence in their asymptotic behavior at large x, which is also exhibited in the UV divergences of their local limit. This is a sign of the potential mixing with gauge variant operators. For these operators, it shall also be possible to work out an appropriate RI/MOM renormalization and matching, but one needs to take into account the gauge variant operators that are allowed to mix with the original operators. This makes the situation much more complicated and is beyond the scope of the present paper. We leave it to future work. of the gluon quasi-PDF operators, the distribution part reads xf (1,1) xf (4,1) x∆f (1,1) The quark matrix elements of the gluon quasi-PDF operators are given as follows: xf (1,1) xf (1,1) g/q,t (x, ρ) = (1−ρ) 5/2 + ρ(2ρ−5)+8(ρ+2)x 3 −4(ρ+8)x 2 +2(ρ+8)x 2(ρ−1) 2 (ρ+4x 2 −4x) , x < 0, xf (4,1) g/q,t (x, ρ) = 0, (A17) xf (4,1) xf (4,1) x∆f (1,1) g/q,z (x, ρ) = g/q,p (x, ρ) = x∆f (2,1) g/q,p (x, ρ) = x∆f (3,1) g/q,p (x, ρ) = x < 0.