Factorization Theorem Relating Euclidean and Light-Cone Parton Distributions

In a large-momentum nucleon state, the matrix element of a gauge-invariant Euclidean Wilson line operator accessible from lattice QCD can be related to the standard light-cone parton distribution function through the large-momentum effective theory (LaMET) expansion. This relation is given by a factorization theorem with a non-trivial matching coefficient. Using the operator product expansion we prove the large-momentum factorization of the quasi-parton distribution function in LaMET, and show that the more recently discussed Ioffe-time distribution approach also obeys an equivalent factorization theorem. Explicit results for the coefficients are obtained and compared at one-loop. Our proof clearly demonstrates that the matching coefficients in the $\overline{\rm MS}$ scheme depend on the large partonic momentum rather than the nucleon momentum.


I. INTRODUCTION
Parton distribution functions (PDFs) are key quantities for gaining an understanding of hadron structure and for making predictions for the cross sections in highenergy scattering experiments. In QCD factorization theorems for hard scattering processes [1], the relevant PDFs are defined in terms of the nucleon matrix elements of light-cone correlation operators. For example, in dimensional regularization with d = 4−2 , the bare unpolarized quark PDF is where x is the momentum fraction, the nucleon momentum P µ = (P 0 , 0, 0, P z ), ξ ± = (t ± z)/ √ 2 are the lightcone coordinates, and the Wilson line is Most often the bare PDF is renormalized in the MS scheme to obtain q(x, µ), and this renormalized PDF is used to make predictions for experiment. The relation is q(x, ) = 1 x dy y Z MS x y , , µ q(y, µ) , where µ is the renormalization scale, and we have suppressed the flavor indices in the renormalization constant Z MS and PDFs. In light-cone quantization with A + = 0, the MS definition has an interpretation as a parton number density. So far our main knowledge of the PDF is obtained from global fits to deep inelastic scattering and jet data, see for example [2][3][4][5][6]. On the other hand, calculating the PDF from first principles with QCD has been an attractive subject, which can for example provide access to spin and momentum distributions that are hard to determine experimentally. Several different approaches to this have been considered using the lattice theory which is a nonperturbative method to solve QCD. Since the lattice theory is defined in a discretized Euclidean space with imaginary time, it is very difficult to calculate Minkowskian quantities with real-time dependence such as the PDF. The first and most well explored option is calculating the moments of the PDF [7][8][9][10][11] that are matrix elements of local gauge-invariant operators. However, since the lattice regularization breaks O(4) rotational symmetry, the consequent mixing between operators of different dimensions makes it difficult to compute higher moments, which in practice has limited the amount of information that can be extracted from this approach. A method to improve this situation by restoring the rotational symmetry has been proposed in Ref. [12]. Other proposals include extracting the PDF from the hadronic tensor [13][14][15][16] and the forward Compton amplitude [17], possibly with flavor changing currents [18], and the more general "lattice cross sections" [19,20]. Systematic lattice analyses of these approaches are under investigation, but challenges remain.
In Ref [21] Ji proposed that the x-dependence of the PDF can be extracted from a Euclidean distribution on the lattice, which can be understood in the language of the large momentum effective theory (LaMET) [22]. This Euclidean distribution is referred to as the quasi-PDF, whose bare matrix element is defined using a spatial cor-arXiv:1801.03917v2 [hep-ph] 15 Sep 2018 relation of quarks along the z direction, e ixP z z P |ψ(z)ΓU (z, 0)ψ(0)|P , where Γ = γ z , z µ = ze µ , e µ = (0, 0, 0, 1), and the Wilson line is For finite momentum P z ,q(x, P z , ) has support in −∞ < x < ∞. According to Ref. [23], there is a universality class of operators that can be considered. For example, for the quasi-PDF, one could also replace Γ = γ z by Γ = γ 0 in Eq. (4) as both definitions reduce to the PDF under an infinite Lorentz boost along the z direction. Unlike the PDF in Eq. (1) that is invariant under the Lorentz boost, the quasi-PDF depends dynamically on it through the nucleon momentum P z . When the nucleon momentum P z is much larger that the nucleon mass M and Λ QCD , which is an attainable window on the lattice, the quasi-PDF can be factorized into a matching coefficient and the PDF [21,22]. The factorization formula is where the renormalized quasi-PDFq(x, P z , µ R ) is defined in a particular scheme at renormalization scale µ R , and O(M 2 /P 2 z , Λ 2 QCD /P 2 z ) are power corrections suppressed by the nucleon momentum. In general the result for C will depend on the choice of Γ = γ z or γ 0 and renormalization schemes. For Γ = γ z the matching coefficient C has been computed for the iso-vector quark quasi-PDF at one-loop level, first with a transverse momentum cutoff in [24], confirmed in [19,25], and also recently determined in the regularization-invariant momentum subtraction (RI/MOM) scheme [26]. Matching for the gluon quasi-PDF is calculated in [27,28]. The matching coefficient C is independent of the choice of states used forq and q. 1 Since matching calculations are carried out with quark states of momentum p z , it can be tricky to know what the right choice to make is for C, and in some of the literature the choice of p z = P z has been suggested when utilizing C for the hadronic nucleon state. This is for example the case in the original quasi-PDF papers [21,22,24] and in the pioneering lattice calculations of the PDF from the quasi-PDF in [25,[29][30][31][32][33][34][35][36], which was summarized in Ref. [37]. In the quasi-generalized parton distribution analysis in [38] it was observed that one should take p z = |y|P z . Through our rigorous analysis of Eq. (6) we show that the correct result for this equation is indeed p z = |y|P z .
Recently, a different procedure [39] to extract PDFs from the same lattice QCD matrix element as in [21] has been proposed based on the Lorentz invariant variables of the spatial correlator (or pseudo-PDF), in place of the quasi-PDF. In this approach, one starts from the spatial correlatorQ γ µ defined for µ = 0 or µ = z by which depends on the two Lorentz invariants z 2 and ζ = −z · P = P z z; the latter is also called the Ioffe time. For an arbitrary Dirac matrix Γ the operatorÕ Γ is defined asÕ This is the same spatial correlator (calculable on lattice) used to define the quasi-PDF in Eq. (4), where P z is fixed and one Fourier transforms with respect to z. If instead z 2 is fixed, and we Fourier transform from the Ioffe time ζ-which is in principle integrating over P z -to the momentum fraction x, then one obtains the pseudo-PDF [39], For arbitrary finite z, the pseudo-PDF only has support in −1 ≤ x ≤ 1 [40,41], but has no parton model interpretation. (Again the pseudo-PDF can equally well be considered for Γ = γ z .) The spatial correlator or pseudo-PDF approach has been explored on the lattice [42,43], where the short distance behavior was explored. The PDF corresponds to the situation when z µ is light-like, in which case the space-time correlator is referred to as the Ioffe-time distribution [44], When Fourier transformed this correlation gives the PDF In short the quasi-PDF and pseudo-PDF are different representations of the Euclidean spatial correlator, as summarized in Table. I. It was pointed out in Ref. [45] that to obtain enough information to extract the PDF for the spatial correlator with small z 2 , one has to do lattice calculations with large momenta P z , which is the same requirement as for the quasi-PDF. Ref. [45] also proposed that the renormalized Distribution Fourier transform Arguments from spatial correlator Spatial correlation pseudo-PDF satisfies the following small z 2 factorization, which they verified at order O(α s ) for the unpolarized iso-vector case with Γ = γ 0 . (Again the coefficient C will depend on the choice of Γ = γ 0 or γ z .) In Ref. [19] a diagrammatic derivation of the factorization formula in Eq. (6) for the quasi-PDF was given. Here we derive this factorization formula for the quasi-PDF in an alternate manner, and also show that spatial correlator and pseudo-PDF are different representations of the same fundamental factorization. Our approach is based on the operator product expansion (OPE) for spacelike separated local operators [46]. For such operators the OPE has been proven for scalar field theory to all orders in perturbation theory [47][48][49], and is widely assumed to hold for any renormalizable quantum field theory including QCD. By introducing auxiliary fields in place of the Wilson line [50], the correlator in Eq. (8) is known to be equivalent to a product of local renormalizable operators of this type. Through our derivation we find the explicit form of the large P z and small z 2 factorization formulas in Eq.(6) and Eq. (12) respectively, as well as the relationship between the matching coefficients C and C. Since the requirement for large P z and small z 2 is the same for both the quasi-PDF and pseudo-PDF approaches, there is in principle no fundamental difference in applying either one to lattice calculations of the proton matrix element ofÕ Γ (z). It is interesting to compare both approaches utilizing the same lattice data, although they shall not yield different result in principle.
The rest of this paper is organized as follows: In Sec. II, we use an OPE ofÕ Γ (z) to derive the large P z factorization of LaMET in Eq. (6) and small z 2 factorization of the pseudo-PDF in Eq. (12). We prove that one must take p z = |y|P z in Eq. (6), so the corresponding argument in C is µ/(|y|P z ). (This OPE approach was used recently in Ref. [20] to prove the factorization theorem for the "lattice cross sections", and the OPE proof carried out here was done independently and first presented in Ref. [51].) In Sec. III, we derive the spatial correlator, pseudo-PDF, and quasi-PDF distributions and matching coefficients at one-loop in MS and analyze the Fourier-transform relation between the quasi-PDF and pseudo-PDF. Unlike earlier results for the quasi-PDF in MS, we also use dimensional regularization with minimal subtraction to renormalize divergences at x = ±∞. In Sec. IV, we discuss how renormalization schemes other than MS are easily incorporated into the factorization formulas. In Sec. V we carry out a numerical analysis of the matching coefficients, by computing the convolution in Eq. (6) numerically using the PDF determined by global fits [4]. We show that the difference between using p z = P z and p z = |y|P z in Eq. (6) is an important effect, and that our MS matching coefficients are insensitive to cutoffs in the convolution integral. In Sec. VI, we discuss the implications of our OPE analysis for the lattice calculation of the PDF in both the quasi-PDF and pseudo-PDF approaches. Finally, we conclude in Sec. VII.

II. FACTORIZATION FROM THE OPE
In this section we make use of the operator product expansion to derive the matching relation for the quasi-PDF, as well as the equivalent matching relations for the spatial correlator and pseudo-PDF. For simplicity these three equivalent cases are presented in separate subsections.

A. OPE and Factorization for the Spatial Correlator
The OPE is a technique to expand nonlocal operators with separation z µ in terms of local ones in the Euclidean limit of z 2 → 0. It can be applied to both bare regulated operators as well as renormalized operators, and our focus will be on the latter. For the gauge-invariant Wilson operatorÕ Γ (z), it was proven that it can be multiplicatively renormalized in coordinate space as [52,53] where δm captures the power divergence from the Wilson line self-energy, Z ψ,z only depends on the end points z, 0 and renormalizes the logarithmic divergences. This multiplicative renormalization was also discussed earlier in Refs. [54][55][56]. For simplicity, in this section we take Γ = γ z for Eq. (13).
In the MS scheme, the power divergence vanishes, and using the OPE the renormalizedÕ Γ (z, µ) can be expanded in terms of local gauge-invariant operators as where µ 0 = z, C n = 1 + O(α s ) and Here Z ij n+1 = Z ij n+1 (µ, ) are multiplicative MS renormalization factors and (µ 0 · · · µ n ) stands for the symmetrization of these Lorentz indices.
The above OPE is valid for the operator itself, where we implicitly constrain ourselves to the subspace of matrix elements for which the twist expansion is appropriate. In the iso-vector case, the mixing with the gluon operators is absent, which we will stick to for the rest of the paper. When O µ0µ1···µn 1 is evaluated in the nucleon state, where a n+1 (µ) is the (n + 1)-th moment of the PDF, and the explicit expression of the trace term have been derived in Ref. [30,57,58]. The inverse relation to Eq. (17) is that q(x, µ) has an expansion with terms proportional to the n'th derivative of the δ-function, as in δ (n) (x) a n (µ), without any nontrivial short distance Wilson coefficient. As pointed out in Ref. [45], to obtain enough information for the spatial correlator at |ζ| = |P z z| ∼ 1 at small z 2 , we have to choose P z large compared to the scale Λ QCD . When P 2 z {Λ 2 QCD , M 2 }, the trace terms in Eq. (16) are suppressed by powers of M 2 /P 2 z , while the contributions from higher-twist operators in Eq. (14) are suppressed by powers of Λ 2 QCD /P 2 z or z 2 Λ 2 QCD . Therefore, the twist-2 contribution is the leading approximation of the nucleon matrix element P |Õ γ z (z)|P at large momentum. From now on we will drop all the power corrections from our discussion.
The Wilson coefficients C n (µ 2 z 2 ) in the OPE ofÕ γ z (z) can be computed in perturbation theory for µ ∼ 1/|z| Λ QCD . In the MS scheme, the C n are log-singular near z 2 = 0, and so is P |Õ γ z (z, µ)|P . For this reason the x-moments of the quasi-PDFq(x, P z , µ) are proportional to C n | z=0 which is divergent, and the quasi-PDF will not simply become the PDF in the infinite P z limit. Instead, we need a factorization formula which matches the quasi-PDF to the PDF. In contrast, for the pseudo-PDF the moments do exist since we hold µ 2 z 2 fixed when taking the x-moment. However, we still need a factorization formula to match the pseudo-PDF to the PDF. We will comment further about this below. Based on Eqs. (14)(15)(16)(17), we can write down the leadingtwist approximation to the spatial correlator as It should be noted that the only approximation we have made so far is ignoring the higher-twist effects that are suppressed by small z 2 and the large momentum P z of the nucleon. In the limit of P 2 z M 2 , Λ 2 QCD , we have P 0 ∼ P z , so even if µ 0 = 0 in Eq. (14), the leading approximation ofÕ γ 0 (z) is still given by the twist-2 contributions in Eq. (18), just with modified coefficients C n .
Based on the OPE results in Eq. (18), we can derive a factorization formula for the Euclidean spatial correlator. First of all, let us define a function C(α, µ 2 z 2 ): From Eq. (18) and the renormalized analog of the Fourier-transform relation in Eq. (9) C(α, µ 2 z 2 ) corresponds to a pseudo-PDF in the special case where a n+1 (µ) = 1. The analysis of Refs. [40,41] implies that the support of C(α, we find from Eq. (18) that Finally, using the inverse transform of the renormalized analog of Eq. (11), we obtaiñ The result in Eq. (23) is the factorization formula for the lattice calculable spatial correlatorQ(ζ, µ 2 z 2 ) which expresses it in terms of the light-cone correlation Q(ζ, µ) that defines the PDF. It has the same structure as the factorization formula for the spatial correlator used for the calculation of the pion distribution amplitude in Ref. [59,60].

B. Factorization for the quasi-PDF
The renormalized quasi-PDF is defined as a Fourier transform of the renormalized spatial correlator, Note that we could use eitherQ γ z orQ γ 0 here. Using the result for the spatial correlator in Eq. (18) this gives Already, one can see that the matching kernel is a function of x/y and µ/(|y|P z ). We define the kernel as and then Eq. (25) can be rewritten as which is the MS factorization formula for the quasi-PDF. This result shows that the factorization formula in Eq. (6) must have p z = |y|P z for the quasi-PDF in the MS scheme. We will show that this remains true for any quasi-PDF renormalization scheme in Sec. IV. This differs from the choice p z = P z which had been conjectured and used in the early papers on the quasi-PDF [21,22,24]. Physically the correct result in Eq. (27) can be understood as the fact that the matching coefficient is only sensitive to the perturbative partonic dynamics, and hence it is the magnitude of the partonic momentum |y|P z which appears, rather than the hadronic momentum P z . Taking the moment of the quasi-PDF using Eq. (24) gives Since the C n coefficients have ln(ζ 2 ) dependence, the derivative for n = n will always have a logarithmic singularity as ζ → 0, and there will be even more singular terms for n < n. This explains why the short distance Wilson coefficient causes the moments not to exist for the quasi-PDF.

C. Factorization for the pseudo-PDF
The renormalized pseudo-PDF is the Fourier transform of the renormalized spatial correlator Since both the pseudo-PDF and spatial correlator are multiplicatively renormalized in a ζ-independent manner, this follows immediately from Eq. (9). If we take Eq. (23) and Fourier transform the spatial correlatorQ(ζ, µ 2 z 2 ) into the pseudo-PDF, and light-cone correlation Q(αζ, µ) into the PDF, then we immediately obtain the factorization formula for the pseudo-PDF, which is the small z 2 factorization formula in Eq. (12). The upper and lower limits of the integrals in Eq. (30) follow immediately from the support −1 ≤ α ≤ 1 of the matching coefficient C(α, z 2 µ 2 ), and we recall that we also have −1 ≤ x ≤ 1 for the pseudo-PDF on the LHS. Since the range of x is bounded for the pseudo-PDF the terms in the series expansion of the exponential in Eq. (29) exist, Comparing with Eq. (18) this implies that the moments of the pseudo-PDF are given by So far we have proven the large P z factorization of the quasi-PDF and small z 2 factorization of the spatial correlation and pseudo-PDFs. After deriving one factorization, it immediately leads to the others, since they are just different representations of the same spatial correlator. Indeed, we see that the quasi-PDF and pseudo-PDF are related at leading power by their definitions: where we have used z = ζ/P z . Based on Eq. (19) and Eq. (26), the Wilson coefficients in their factorization theorems also maintain the same relationship, For the relations in Eqs. (33) and (34) the same choice of Γ = γ 0 or γ z should be used in the quasi-and pseudo-PDFs, or their corresponding coefficients.
In summary, there is a unique factorization formula that matches the quasi-PDF, spatial correlator and pseudo-PDF to the PDF. Since their factorizations into the PDF all require small distances and have large nucleon momentum, the setup for their lattice calculations must also be the same. Therefore, the LaMET and pseudo distribution approaches are in principle equivalent to each other, and they differ perhaps only by effects related to their implementation on the lattice.
In Ref. [39] it was speculated that one can study a ratio functionQ on a lattice with spacing a, and the O(z 2 ) corrections may cancel approximately. This idea was tested in Ref. [42] in lattice QCD, and the results show that the ratio evolves slowly in z 2 at small values. It is then interesting to consider what type of non-perturbative information can be extracted from this ratio.
This question can be answered using the small z 2 factorization for the spatial correlator. According to where in the MS scheme to one-loop which was also derived recently in [61]. Then the ratio becomes Using Eq. (32) and our MS one-loop perturbative pseudo-PDF result in Eq. (54) below we find for Γ = γ 0 that where the Harmonic numbers are H n = n i=1 1/i and H which also modifies Eq. (37) for n = 0. At small z 2 where the perturbative expansion with µ 1/|z| is valid, the ratio in Eq. (38) has a weak logarithmic dependence on |z|, which is consistent with the lattice findings in Refs. [42,43]. The weak dependence on |z| can be quantitatively described by an evolution equation in ln z 2 [39,62]. According to our OPE analysis, herẽ Q(0, µ 2 z 2 ) only serves as an overall normalization factor which is contaminated by higher-twist corrections, and the ln z 2 evolution can be put in accurate terms with the factorization formula in Eq. (23) that enables us extract the PDF from the ratio function. The same point was demonstrated by work done very recently in [61], which appeared simultaneously with our paper.

III. EQUIVALENCE AT ONE-LOOP ORDER
As has been proven in Sec. II, the quasi-PDF and pseudo-PDF as well as their matching coefficients are related by a simple Fourier transform in Eq. (33) and Eq. (34). This relation is valid to all orders in perturbation theory. In this section we check the relations in Eq. (33) and Eq. (34) at one-loop order. We choose Γ = γ 0 for our main presentation, but also quote final results for the case Γ = γ z .
In the Feynman gauge, we calculate the quark matrix elements of the unpolarized iso-vector quasi-PDF, pseudo-PDF, and light-cone PDF at one-loop order in dimensional regularization with d = 4 − 2 . The external quark state is chosen to be on-shell and massless, and we regularize the UV and collinear divergences by 1/ UV ( UV > 0) and 1/ IR ( IR < 0) respectively. The one-loop order Feynman diagrams are shown in Fig. 1. One-loop Feynman diagrams for the quasi-PDF, spatial correlator and pseudo-PDFs. The first one is named "vertex", the second and third ones are named "sail", and the last one "tadpole". The standard quark self energy wavefunction is also included.
In an on-shell quark state with momentum p µ = (p 0 = p z , 0, 0, p z ), for Γ = γ 0 , each diagram gives where ι = e γ E /(4π) is included to implement µ in the MS scheme, C F = 4/3, and T a is the SU(3) color matrix in the fundamental representation. The second term in the brackets in the last line, which is proportional z, does not contribute to the loop integral as it is odd under the exchange of p z −k z → −(p z −k z ). The quark self-energy correction isQ After carrying out the loop integrals in Eq. (41) according to the method in Ref. [45], we obtaiñ For simplicity we have left out the tree level multiplicative spinor factorūγ 0 u when quoting one-loop results in Eq. (43), and will continue to do so for the spatial correlator, quasi-PDF, and pseudo-PDF results quoted below. Since α bare s = α s (µ)µ 2 Z 2 g = α s (µ)µ 2 + O(α 2 s ) is µ-independent we do not include µ as an argument in the bare functions. In the final result for each term we have also specified whether 1/ factors (that remain after expanding about → 0) are IR or UV divergences. Combined with the wavefunction corection, the bare spatial correlator Note that as → 0 the terms in the innermost curly brackets have no 1/ term. Also we can verify that in the local limit of the operator that the bare one-loop correction vanishes as expected by conservation of the vector current: where we note that it is important that the 1/ IR terms cancel since the assumption > 0 is only valid for the 1/ UV term. The corresponding bare pseudo-PDF is Since we will encounter plus functions over different domains below, we define a plus function at x = x 0 within a given domain D so that (See App. C for more details.) It is straightforward to confirm that the bare pseudo-PDF satisfies the local vector current conservation, lim z→0 dx P (1) (x, z 2 µ 2 , ) = 0, with the same cancellation as in Eq. (45). Now, according to the relations between the quasi-PDF and the spatial correlator or pseudo-PDFs in Eqs. (24,33), we can do a Fourier or double Fourier transform of the results in Eqs. (44,46) to get the quasi-PDF. Despite its straightforwardness, the Fourier transform is subtle and the details are provided in App. A. Here we simply quote the result for the bare quasi-PDF, where After some algebra one can confirm that the bare quasi-PDF satisfies local vector current conservation, with dxq (1) (x, p z , ) = 0. To verify this result one must carefully separate out 1/ UV factors arising from requiring > 0 to obtain convergence at x = ±∞, and 1/ IR factors that arise from requiring < 0 to obtain convergence at x = 1.
An alternate method of obtaining the quasi-PDF is to directly calculate it from the Feynman diagrams by first Fourier transforming z into xp z . As a result, the factors (e −ip z z −e −ik z z ) are transformed into [δ(p z −xp z )−δ(k z −xp z )], and all the loop integrals reduce to (d − 1)-dimensional ones. This is the procedure for the matching calculations of the quasi-PDF used in Refs. [24,26,27], and is distinct from doing the Fourier transformation after fully carrying out the integrals as in Eqs. (43)(44)(45)(46)(47)(48). As a cross-check we have confirmed in App. B that we obtain the exact same bare quasi-PDF in Eq. (48) from both procedures. Now we consider the expansion to obtain MS renormalized results for the spatial correlator, pseudo-PDF, and quasi-PDF. Expanding the spatial correlator in we obtaiñ with the MS counterterm and renormalized spatial correlator given by Here 3 F 3 is a hypergeometric function and the Fourier transform of (1 + For the position space PDF we have Next we expand the bare pseudo-PDF from Eq. (46) in to obtain the MS counterterm and renormalized pseudo-PDF as P (1) Note that the renormalized MS pseudo-PDF depends explicitly on µ 2 , and satisfies the relation to the renormalized MS spatial correlator given in Eq. (29). It is also interesting to note that having expanded in , local vector current conservation is no longer satisfied by the limit of the renormalized MS pseudo-PDF, since gives a divergent result. The same divergence is present in the one-loop MS renormalized spatial correlator in Eq. (51).
Although this is the case in MS, it does not need to be the case in other renormalization schemes. For the quasi-PDF there are two methods that we can consider for the renormalized calculation, either expanding the bare result in Eq. (48) and renormalizing in (x, p z ) space, or following our preferred definition in Eq. (24) and Fourier transforming the renormalized spatial correlator in Eq. (51). Although these two approaches will lead to the same final result for C for practical applications, there is a subtle difference that we will explain.
First consider the renormalization of the quasi-PDF done in (x, p z ) space. Expanding Eq. (48) in , and writing q (1) (x, p z , ) = δq (1) (x, µ/|p z |, UV ) +q (1) (x, µ/|p z |, IR ) + O( ) allows us to identify the MS counterterm and renormal-ized quasi-PDF as The details of working out the expansion of Eq. (48) are provided in App. C, including definitions of the plus functions and δ-functions at x 0 = ±∞ that appear in the result quoted here. The MS quasi-PDF obtained in Eq. (56) still satisfies vector current conservation This is obviously the case for the plus function terms which individually integrate to zero, and is also true for the combination of δ-functions which appears in Eq. (56). The renormalized MS quasi-PDF in Eq. (56) differs slightly from that obtained using our definition in Eq. (24). Using Eq. (24) and the renormalized spatial correlator in Eq. (51) we instead obtain To carry out this calculation we defined the Fourier transformation of the singular function ln(ζ 2 ) as where we have used the results in Eqs. (A1,C7) to derive the second and last equalities, and took the limit η → 0 + or η → 0 − when needed. Thisq (1) (x, µ/|p z |, IR ) does not satisfy vector-current conservation, and is different from q (1) (x, µ/|p z |, IR ) in Eq. (56) only by the δ-functions at x 0 = ±∞. Within the function domain −∞ < x < ∞, they are exactly the same. We will see below that both Eq. (56) and Eq. (58) eventually lead to the same result for the one-loop matching coefficient. The final ingredient we need for the matching calculations is the PDF, whose one-loop bare matrix element can be written as a sum of an MS counterterm and renormalized matrix element, q (1) (x, ) = δq (1) With the above results in hand we can now determine the matching coefficients up to one-loop order. Using Eq. (30) we find Therefore the matching coefficient relating the pseudo-PDF and PDF in the MS scheme with Γ = γ 0 is 2 This result is independent of the infrared regulator as it must be. We have also computed the matching coefficient for the Γ = γ z case, and it is C γ z (α, z 2 µ 2 ) = C(α, Due to the ln(z 2 µ 2 )δ(1 − α) term in Eq. (62), the matching coefficient for the MS pseudo-PDF again displays the fact that there is not a smooth local limit as z → 0. It is possible to define a scheme other than MS to ensure that this limit is smooth, reproducing a renormalization for z → 0 that agrees with the fact that the local operator corresponds with a conserved current. One such scheme would be to simply multiply all MS renormalization constants by C 0 (µ 2 z 2 ), which would lead to a spatial correlator renormalized in a different scheme, and a corresponding different matching coefficient in Eq. (62) with a smooth z → 0 limit. This is equivalent to studying the ratio of Eq. (38) from the start as advocated in Ref. [39,62]. We will give explicit results for this scheme choice below. This modified scheme should not be confused with the strict definition of the MS scheme. From Eq. (27) the corresponding relation for the matching coefficient for the quasi-PDF defined in Eq. (24) is Therefore using Eq. (58) the matching coefficient relating the quasi-PDF and PDF is 2 A one-loop analysis of the spatial correlator in the coordinate space also recently appeared in Refs. [62,63]. Our factorization result for the spatial correlator in Eq. (23) has a similar form to the hard part of the reduced spatial correlator found in Eq.(3.35) of Ref. [62] and Eq.(17) of Ref. [63]. It is therefore interesting to compare our C(α, z 2 µ 2 )/C 0 (µ 2 z 2 ) and this hard part. Again this result is independent of the IR regulator as it must be. Here the plus function terms g 1 (ξ) [1,∞] +(1) and g 2 (ξ) [−∞,0] + (1) have integrands that converge for ξ → ±∞, behaving as g i (ξ) ∼ 1/ξ 2 . Note that if we had instead used the renormalized MS quasi-PDF calculated in Eq. (56), we would obtain a different matching coefficient C with different δ-functions at ξ = ±∞. However, the δ-functions do not contribute to the convolution integral in Eq. (27) for any integrable PDFs. For example, to carry out the convolution with 1/ξ 2 δ + (1/ξ) we can use δ + 1 ξ = lim β→0 + δ 1 ξ − β , which when plugged into the factorization formula gives lim β→0 + dy |y| For the plus-function at ∞ using Eqs. (C9) and (C3) we have In the last line we dropped the θ(x/y − β) since at small y our PDF behaves as f u−d (y) ∼ y −1+a with 0 < a < 1. This also implies which means that the distribution contributions evaluated at ξ = ±∞ in the matching coefficient C give zero contribution. Therefore, the matching coefficients calculated from the quasi-PDFs in Eq. (56) and Eq. (58) are the same in effect, and we can simply drop all the δ-functions at ξ = ±∞ when plugging them into the factorization formula: The use of Eq. (68) in the factorization formula is valid for any PDF that behaves as lim y→0 f (y, µ) ∼ y −1+a with a > 0. We have also computed the matching coefficient for the Γ = γ z case, and it is given by C γ z (ξ, µ/(|y|P z )) = C(ξ, µ/(|y|P z )) + ∆C γ z (ξ, µ/(|y|P z )) with Note that our result for the quark matching coefficient in MS differs from that of Ref. [27] which is a pure plus function, but gives a convolution that does not converge, just as in the case of the quasi-PDF with a transverse momentum cutoff, see Ref. [26].
Besides, if one uses a scheme other than MS for the quasi-PDF, such as the scheme obtained by absorbing C 0 into the MS renormalization constant, then this will lead to a result for the matching coefficient that is a pure plus function and hence satisfies current conservation. Starting with Eq. (65) and using Eq. (37) together with Eq. (59) we obtain and for the Γ = γ z case, While retaining current conservation in the renormalized quasi-PDF, Eq. (70) can be used for example as input to the two-step matching procedure in the lattice calculation of PDF in Refs. [33]. For the matching step, an equivalent procedure is to study the ratio given in Eq. (38) in the MS scheme from the start, as advocated in Ref. [39,62], performing its matching onto the PDF, which will yield Eq. (70). This concludes our discussion of matching results and the equivalence between the quasi-PDF and pseudo-PDF at one-loop order.

IV. OTHER RENORMALIZATION SCHEMES
Although we derive the above matching formula assuming that the quasi-PDF is renormalized in the MS scheme, this is not a limitation to our result. Since the gauge-invariant Wilson line operatorÕ Γ (z) has been proven to be multiplicatively renormalizable in the coordinate space [52,53], one can convertQ Γ (z) from any other scheme to the MS scheme before using the above factorization formula. The renormalization of the quasi-PDF has been studied in many recent papers [26,28,33,34,[54][55][56][64][65][66]. We will discuss some of these results and show how they can be incorporated into the factorization formula in Eq. (27).
The MS scheme is convenient for our discussion of the OPE as it guarantees Lorentz and gauge invariances, but it is not practical for lattice renormalization. Since the lattice theory has a natural UV cut-off 1/a with a being the lattice spacing, the unrenormalized spatial correlator Q inherits the power divergence from the Wilson line selfenergy according to Eq. (13). For an arbitrary scheme X, the renormalized spatial correlator should be free of all the UV divergences and have a welldefined continuum limit as a → 0. This continuum limit, in particular, is independent of the UV regulator, so As a result, we can relateQ X (ζ, z 2 µ 2 R ) to the MS scheme by the conversioñ where the regulator dependence is completely canceled out between Z MS and Z X . The ratio Z X can be calculated perturbatively in QCD, which was done in [56] for several lattice schemes and the RI/MOM scheme. Thus the factorization formula we have proven in Sec. II still applies toQ X with a slight modification to the coefficient function, where the matching coefficient for the scheme X is related to that of MS by For the pseudo-PDF the modified result also involves this same coefficient Meanwhile, for the quasi-PDF we have, Here the modified coefficient for the X scheme is related to coefficient in the MS scheme by where hereZ X is defined by the Fourier transform Depending on the scheme X we note that slightly modified definitions ofZ X may be more appropriate.
One undesirable feature of the MS scheme for the renormalized spatial correlator is that it does not have a smooth z → 0 limit, and hence no simple connection with the fact that the local operator for z = 0 is a conserved current. To avoid this one can simply make use of a different scheme that has a simple relation to MS, such as by adding C 0 (µ 2 z 2 ) to the MS renormalization constant. This removes the offending ln(µ 2 z 2 ) terms and yields a scheme with a smooth connection to the conserved current.
Besides the MS scheme, the quasi-PDF has also been defined with a transverse momentum cut-off [19,24,25,29] and in the RI/MOM scheme [26,33,34,56,66]. The RI/MOM scheme has attracted strong interest recently as it can be implemented nonperturbatively on the lattice, so we consider it as an explicit example of the above relations. In this scheme, the renormalization constant Z OM is determined by imposing a condition on the spatial correlator in an off-shell quark state, where q denotes the quark state, p µ is the external momentum, and "(0)" in the superscript stands for the treelevel matrix element. As a result, and here we definē Then the matching coefficient in Eq. (78) becomes The choices of µ R and p z R are independent of µ and P z , and p z R = P z was used in Refs. [26,34]. It should be noted that on the lattice, due to the breaking of chiral symmetry, the vector-like quark Wilson line operatorÕ γ µ (z) can mix with the scalar operatorÕ 1 (z), as has been discussed in Refs. [33,34,56,66,67]. After considering the mixing effects, the same factorization formula can still be applied to the RI/MOM quasi-PDF from lattice QCD.

V. NUMERICAL RESULTS
In this section we numerically analyze the quasi-PDF, spatial correlator and pseudo-PDF by studying how the matching coefficients in Eqs. (27,30) change the PDF. The quasi-PDF has already been studied in this manner for the MS, transverse momentum cut-off, and RI/MOM schemes in Ref. [26]. Our new MS result for the matching is given in Eq. (68), and leads to stable convolution integrals. We also compare the differences between using hadron momentum p z = P z and the parton momentum p z = |y|P z for the matching coefficient in the MS scheme. We take Γ = γ 0 for the results here.
As an example we use for our analysis the unpolarized iso-vector parton distribution, where we include fū(−x, µ) = −fū(x, µ) and fd(−x, µ) = −fd(x, µ), the anti-parton distributions. For ease of comparison, we use the next-to-leading-order iso-vector PDF f u−d from MSTW 2008 [4] with the corresponding running coupling α s (µ).
To implement the plus functions in the numerical calculation, we impose a soft cutoff |y − x| < 10 −m and test the sensitivity of results to m. Since the limit of y → 0 corresponds to the asymptotic region |x/y| → ∞, we also impose a UV cutoff |y| > 10 −n to test the convergence of the convolution integral. We find that all the results presented below are insensitive to m and n. The fact that our result in Eq. (68) has terms outside the plus function at 1 in each of the ξ ∈ [1, ∞] and ξ ∈ [−∞, 0] intervals is important for ensuring that our MS result for C is insensitive to the |y| > 10 −n cutoff. This was not the case for the quasi-PDF that was defined with a transverse momentum cutoff [24]. The RI/MOM scheme result [26] also does not suffer from this issue.
In Fig. 2 we compare the PDF with the quasi-PDF in the MS scheme obtained from the convolution in Eq. using our one-loop result in Eq. (68). We observe that changing from p z = P z to the correct p z = |y|P z shifts the result in the physical region by a considerable amount.
The same type of comparison can be made for the pseudo-PDF in the MS scheme by applying the factorization formula in Eq. (30) and matching coefficient in Eq. (62). In Fig. 3a we compare the PDF and pseudo-PDF and their dependence on the factorization scale µ, while in Fig. 3b we include the dependence of the pseudo-PDF on the distance |z|. Since the matching coefficient in Eq.(62) is similar to the parton splitting function except for the nontrivial finite constants, matching the PDF to the pseudo-PDF is analogous to evolving the PDF from µ to the scale of 1/|z|. This evolution has been calculated in Refs. [39]. The variation of |z| has a similar effect to the PDF evolution, as is observed in the right panel of Fig. 3. When |z|µ = 1, the logarithm is zero, and the matching effect from C is determined by the nontrivial constants in Eq. (62), which shifts the PDF downward in the large-x region and upward in the small-x region.
Finally, we can make a similar comparison for the spatial correlator in the MS scheme obtained with Eq. (23) and Eq. (62). Its real and imaginary parts are even and odd functions of ζ respectively, and are shown in Fig. 4. Again we show the residual dependence on µ and |z| which are similar to that for the pseudo-PDF. The matching broadens the curves in the coordinate space. The spatial correlator renormalized in the MS scheme does not exhibit vector current (or particle number) conservation, which can be clearly seen from the fact that the real part of the distribution is not equal to 1 at ζ = 0 (except for the special case where |z|µ is tuned to cancel the constant terms in the one-loop C).

VI. IMPLICATIONS FOR LATTICE CALCULATIONS
Our proof in Sec. II makes clear the relationship between the renormalized quasi-PDF, spatial correlator, and pseduo-PDF distributions. As a practical matter there are a few different ways in which these equations can be used to convert a lattice calculation of the spatial correlatorQ into a PDF. Three examples are 1) first Fourier transform to the quasi-PDF with Eq. (24), and then use Eq. (27), 2) first Fourier transform to the pseudo-PDF with Eq. (29), and then use Eq. (30), and 3) first match to the Fourier transform the position space PDF Q(ζ, µ) using Eq. (23), and then transform it to the PDF with the inverse of Eq. (22). Since the numerical implementation of these steps may have slightly different systematics it is interesting to compare them, or to use more than one approach in order to reduce uncertainties.
According to the analysis in Sec. II, for the factorization formula of the Euclidean distributions to work, one must calculate the same spatial correlator with small distance z 2 and large momentum P z so that the dynamical and kinematic higher-twist effects are suppressed. For practical lattice calculations, this means that there is only a finite number of useful data points in (z, P z ) that we can use to extract the PDF.
To illustrate this, consider a 48 3 × 64 lattice with spacing a = 0.09 fm. The distance of the spatial correlation z is in units of a ∼ 1/2.2 GeV −1 , and the nucleon momentum P z is in units of 2π/(48a) ∼ 0.29 GeV. Let us take Λ QCD ∼ 0.3 GeV. In principal the target mass corrections can be subtracted. If we consider various values z = ma and P z = n * 2π/(48a) for integer m and n, then to satisfy zΛ QCD 1 and P z Λ QCD , we must have To control the higher-twist correction at 20%, i.e.
where the largest value for n is limited by what is practical in current lattice simulations. Six is the largest number of units attained in Ref. [42]. For quasi-PDF calculations, there are 4 × 2 + 1 = 9 useful data points for each fixed momentum |P z |; for pseudo-PDF calculation, there are only 4 × 2 = 8 useful data points for each fixed |z|. In either case, it is anticipated that a direct Fourier transform with respect to z or ζ = zP z will lead to oscillation in x-space and incorrect prediction for the small-x region due to the truncation in coordinate space. This has been observed in a recent lattice calculation of the quasi-PDF in Ref. [34]. Methods have been developed in recent works to eliminate the oscillation from the truncation effect [35,68] in the quasi-PDF, while the higher-twist contributions at large z still need to be systematically corrected. It should be noted that the above is a rough estimate of the higher-twist corrections since the prefactor of z 2 Λ 2 QCD could be smaller than 1. Their actual significance can only be quantitatively determined from lattice simulations.
To fully take advantage of all the useful data points, we can evolve them to either the same z 2 or P z according to the perturbative analysis, which has been studied in Refs. [42,43,62] for the spatial correlator. However, since the evolution equation of the spatial correlator in ln z 2 or ln P 2 z follows a nonlocal convolution in ζ = zP z or z, one has to know the full information in coordinate space to do the evolution. With limited number of data points, either large uncertainties or adopting a model-dependent assumption about the shape is inevitable.
To improve the precision of either approach, the only way forward is to have finer lattice spacing a so that we could have more data points which satisfy |z| Λ −1 QCD and larger nucleon momentum P z . With increasing P z , the valence distribution of the nucleon is contracted in the z direction, so the spatial correlation of valence quarks is shrinked into smaller distance in z. If P z is large enough, the spatial correlation will fall off quickly within |z| < Λ −1 QCD , then the truncation error from Fourier transform will be significantly reduced. On the other hand, if we interpret the spatial correlation as the spatial correlator, its shape will not change under a Lorentz boost because it is a scalar function of ζ = zP z and z 2 . Nevertheless, finer lattice spacing a allows for calculation with a wider range of P z , thus covering larger values of ζ = zP z to reduce the truncation error. Since the number of useful data points increases quadratically with 1/a, a more precise lattice calculation with controlled systematic errors will be available in the future.

VII. CONCLUSIONS
Starting with a Euclidean operator product expansion for products of gauge invariant operators in QCD, we have derived the factorization formulas for the quasi-PDF, spatial correlator and pseudo-PDF. The three Euclidean distribution functions are related observables, and all follow from the same fundamental factorization. For the spatial correlator this derivation implies that the ratio in Eq. (35) does scale in z 2 , but needs the small z 2 factorization formula in Eq. (23) to extract the PDF. Our derivation for the factorization formula applies when the renormalized spatial correlator is defined in any scheme. The OPE used here could also be used to systematically derive factorization formulas for power corrections to Eq. (6), which will involve matching onto higher-twist parton distributions. (The numerical relevance of these corrections is considered in Ref. [30].) Note that LaMET is not equivalent to the expansion from the OPE, as the former is more general and can be applied to the lattice calculations of other quantities, for example the TMD-PDF where a simple OPE does not exist.
Our derivation of the factorization formula for the quasi-PDF also verifies that the parton momentum p z in the matching coefficient in Eq. (6) has to be p z = yP z , which makes a considerable difference for the MS matching result when compared with p z = P z (see Fig. 2). The proper p z should therefore be used in lattice calculations of the PDF in the LaMET approach.
As a non-trivial test of relations between the various distributions and factorization formulas we have considered results at one-loop in the MS scheme. We have derived a corrected results for the coefficient C for the MS scheme, given in Eq. (68), which leads to convergent results in the convolution integral. We have also computed the one-loop MS result for the Wilson coefficient C appearing in the spatial correlator and pseudo-PDF factorizations. A numerical analysis of these one-loop corrections in MS has also been provided. The one-loop matching coefficient C has a smaller effect for the pseudo-PDF than C does for quasi-PDF, as can seen by comparing Figs. 2 and 3. Given systematic uncertainties in manipulating the lattice data, it is potentially interesting to consider using the same lattice data on the spatial correlator to extract the parton distribution function using both the quasi-PDF and pseudo-PDF approaches.
There are several different ways of implementing the factorization formula to calculate the PDF from lattice data for the spatial correlatorQ, which we have discussed in Sec. VI. One always has to work with short distance correlation and large nucleon momentum to reduce higher-twist corrections. This limits the number of useful data points from lattice calculations as described in Sec. VI. To achieve precision calculations without making model assumptions it will be highly desirable to move towards finer lattice spacing to increase the number of effective data points.
FG02-93ER-40762, de-sc0011090 and de-sc0012704, by the Laboratory Directed Research and Development (LDRD) funding of BNL under contract de-ec0012704, by a grant from the National Science Foundation of China (No. 11405104), and within the framework of the TMD Topical Collaboration. I.S. was also supported in part by the Simons Foundation through the Investigator grant 327942. T.I. was also supported by JSPS KAK-ENHI Grant Numbers JP26400261, JP17H02906, and also MEXT as "Priority Issue on Post-K computer" (Elucidation of the Fundamental Laws and Evolution of the Universe) and JICFuS.
while for the sail diagram q (1) sail (x, p z , ) and for the tadpole diagram q (1) tadpole (x, p z , ) After adding these expressions and carrying out the remaining integrations over y, we obtain the same result as in Eq. (48).

Appendix C: Expansion and Plus Functions
Since the support of the quasi-PDF ranges from −∞ to ∞, its asymptotic behavior as ∼ 1/|x| at |x| → ∞ implies a UV divergence which can be regularized by dimensional regularization. Therefore, the expansion of the quasi-PDF should account for this feature.
In general, we need to expand and it is well known that where < 0. Here we follow [69] and the plus functions L n (x) are defined as Note that 1 0 dx L n (x) = 0.
To define the expansion in the range x ∈ [1, ∞] we simply map this interval to [0, 1] via t = 1/x. Taking an arbitrary smooth test function g(x) we have Here > 0 and the superscript + on the δ + function indicates that its argument should be positive. Therefore δ + (1/x) has its support at x = +∞, not x = −∞. Since g is arbitrary we can identify Combining the above results and denoting which 1/ poles are UV or IR divergences we have