Origin and Resummation of Threshold Logarithms in the Lattice QCD Calculations of PDFs

Xiang Gao, 2, ∗ Kyle Lee, 4, † Swagato Mukherjee, and Yong Zhao 5, ‡ Physics Department, Brookhaven National Laboratory, Bldg. 510A, Upton, NY 11973, USA Physics Department, Tsinghua University, Beijing 100084, China Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA Physics Department, University of California, Berkeley, CA 94720, USA Physics Division, Argonne National Laboratory, Lemont, IL 60439, USA


I. INTRODUCTION
Parton distribution functions (PDFs) are among the most basic quantities that describe the inner structure of the proton in high-energy experiments such as the deep inelastic scattering (DIS) and Drell-Yan (DY) processes. Being intrinsic to the hadrons, which is dominated by the strong interaction, the PDFs are nonperturbative functions and can only be determined from experiments or first-principles calculations in Quantum Chromodynacmis (QCD).
Determination of the PDFs from experimental measurements is based on QCD factorization theorems [1], where a cross section can be factorized as a convolution of a hard coefficient function, which can be calculated analytically in perturbation theory, and a PDF. For instance, the factorization of the DIS cross section of a lepton scattering off of a hadron of momentum P , with the exchange of a boson of momentum q, can be schematically written as where Q 2 = −q 2 , η = Q 2 /(2P · q) is the Bjorken variable and µ is the factorization scale. Here q(x, µ) is a PDF at the momentum fraction x, and C is the hard coefficient function. Using similar factorization formulas for different processes to fit to the enormous amount of cross section data accumulated from the worldwide highenergy physics programs over the past few decades, the * xgao@bnl.gov † kylelee@lbl.gov ‡ yzhao@bnl.gov PDFs have been extracted with remarkable precision [2][3][4][5][6][7]. Nevertheless, due to the limitations of experiments, there are still large uncertainties for the PDFs of the sea quarks, gluon, spin-dependent partons, and at small and large x regions. Among them, the extraction of the PDF at large x is complicated by the soft divergences in the hard coefficient functions as their partonic variables approach threshold, for instance, C(z, Q, µ) of Eq. (1) as z approaches one. Such divergence gives rise to threshold logarithms of (1 − z), which necessitates all-order resummation [8][9][10]. The global analysis of the PDF is then affected by whether one uses the fixed-order or the resummed version of the coefficient function C [11][12][13][14]. As can be seen from Eq. (1), since the PDF is in general convoluted with a hard coefficient function, the extraction of the PDF at even moderately large x may be affected by threshold resummation. For example, the NNPDF collaboration used DIS, DY, and inclusive top-quark pair production processes simultaneously to study the impact of threshold resummation in PDF extractions [13], and concluded that it suppresses the PDF in the large-x region around x 0.1, while at intermediate values of x, 0.01 x 0.1, the quark PDFs are somewhat enhanced 1 . The suppression of the extracted PDFs at large x can be intuitively understood as compensation for the enhancement of resummed coefficient functions near threshold.
Among the many PDFs extracted from experiments, those of the pions are of special interests as they are the pseudo-Nambu-Goldstone bosons of QCD and are the simplest hadrons for studying strong dynamics. For example, the behavior of the pion valence PDFs at large x, which is characterized by a power law behavior (1 − x) β , has been under debate among various approaches. As x → 1, the pion momentum is mostly carried by the active valence quark, so the spectator quark is at low energy and interacts strongly with other soft particles. In the asymptotic limit, the Brodsky-Farrar quark counting rules [15] predicted that β ∼ 2, but the analysis of the E-0615 Fermilab data suggested β ∼ 1 [16]. With the inclusion of threshold resummation, the Fermilab data was reanalyzed by the authors of Ref. [12] (ASV), and they found β ∼ 2. Moreover, the recent global Monte Carlo analysis from the JAM collaboration [17], as well as the analysis done using the xFitter program [18], suggest that β ∼ 1, although neither of them included threshold resummation.
Apart from the determination of the PDFs using experimental data, there have also been long-standing efforts to calculate the PDFs from the first-principles method of lattice QCD. Due to their explicit dependence on the real time, the PDFs cannot be directly simulated on the Euclidean lattice that uses imaginary time. Instead, their Mellin moments, which are matrix elements of local gauge-invariant operators, can be calculated on the lattice. However, these calculations are only viable for the lowest few moments due to power divergent operator mixings and worsening of the signal-to-noise ratio for the higher moments (for a review, see Ref. [19]).
In recent years, several approaches have been proposed to calculate either the higher Mellin moments or x-dependence of the PDF. They include the calculation of the hadronic tensor [20][21][22][23], restoration of rotational symmetry to access higher moments [24], operator product expansion (OPE) of a Compton amplitude in the unphysical region [25,26], the coordinate-space OPE or short-distance expansion of the current-current correlator [27,28], the quasi-PDF (qPDF) in large-momentum effective theory (LaMET) [29][30][31], and the pseudo-PDF (pPDF) approach [32,33]. Except for the first two approaches, all the other ones are based on the factorization of a Euclidean observable into the PDF (or its moments) and a perturbative matching coefficient, up to power-suppressed corrections. Among them, the LaMET approach has attracted wide interest and led to much progress towards systematic calculations of the xdependence of the PDF over the past few years. Within this approach, the PDF can be obtained from the qPDF, which is the Fourier transform of a spatial correlator in a boosted hadron state, through an effective theory expansion in powers of large hadron momentum [31]. The pPDF approach is closely related to LaMET as it uses an equivalent short-distance expansion of the same spatial correlator in the coordinate space [34][35][36]. The pPDF in this approach is a different Fourier transform of the same spatial corrrelator.
In both the large-momentum expansion for the qPDF and the short-distance expansion for the pPDF (or spa-tial correlator), one suffers from power corrections at the end-point regions as x → 0 and x → 1 [36][37][38][39]. Apart from the power corrections, the matching coefficient in the leading-power contribution includes large logarithms in the threshold region, which will also impact the determination of the large-x exponent for the PDF. Similar large logarithms are also observed in the matching coefficient for the short-distance expansion of the currentcurrent correlator [40], which, however, is not the focus of this work. The forms of such logarithms are very similar to the threshold logarithms in the DIS and DY cross section formulas, but their origins are not yet understood, especially for the pPDF or spatial correlator, since there is no clear correspondence to a soft-momentum limit in its variables. As will be shown below, both the qPDF and the pPDF are the special limits of a 3D momentum distribution defined with a straight gauge link. This 3D momentum distribution then has a dependence on the variable x with clear physical interpretation as momentum fraction, and includes threshold logarithms of (1−x) in the limit of x → 1. The relation between this 3D distribution and the qPDF or the pPDF then elucidates the origin of the large logarithms found in the matching coefficients mentioned above. Especially, the threhsold logarithm in the matching coefficient for pPDF matches the factorization scale to the ultraviolet (UV) fixed point in the x → 1 limit, which is opposite to the DIS and DY cases where the scale being matched to becomes infrared (IR) in this limit. After identifying the origin of these threshold logarithms, we will then resum them at next-to-leading logarithmic (NLL) accuracy by using the form motivated by the standard threshold resummation formulas in the literature [41], and also confirm our resummation formula with the next-to-next-to-leading order (NNLO) matching coefficients [42]. Unlike the DIS and DY cases, the resummed matching coefficient for the pPDF or spatial correlator is free from the Landau pole as it matches to the UV fixed point.
Since the large-x region of the PDF is sensitive to longrange correlations, which are usually beyond the reach of contemporary lattice data, the predictive power of lattice QCD is limited. Although two existing lattice calculations have tried to fit the large-x exponent of the pion valence PDF [40,43], neither of them has quantified the power corrections or threshold resummation effects, and their results were inconclusive. In this work, we investigate the effects of threshold resummation on the fitted Mellin moments and the large-x exponent of the pion valence PDF for the first time. Despite the fact that current lattice data are only sensitive to lower moments and PDF within a moderate range of x, this investigation is aimed to deepen our understanding of the systematics at large x and to provide the framework for future analyses of data that are sensitive to this region.
The rest of this work is presented as follows: In Sec. II, we show that threshold logarithms in the qPDF and pPDF (or spatial correlator) originate from the leading divergences in the 3D momentum distribution as x → 1.
As a byproduct, we find that the qPDF and pPDF (or spatial correlator) are special limits of this particular 3D distribution, revealing a new relation between the two in addition to Fourier transforms. Then in Sec. III, we resum the threshold logarithms at NLL accuracy with next-to-leading order (NLO) matching, and estimate its impact on the large-x exponent of the PDF. In Sec. IV, we use the resummed formula to fit the second and fourth moments, as well as to fit the pion valence PDF directly with the parametrization ∼ x α (1 − x) β , and show that the effect of threshold resummation on the moments or value of β is negligible compared to other uncertainties in the current lattice data. Finally, we conclude in Sec. V.
The MS matching coefficients C and C both have been derived up to two-loop order [34,36,42,47,48]. At one-loop order, the matching coefficients in the threshold limit behave as where C F = 4/3, and ξ can approach one from both above and below. The plus functions labelled by "+" are the usual ones defined within the region [0, 1], while the definition of those labelled by "⊕" is for arbitrary function g(x) with x ∈ (−∞, ∞). As one can see, both Eqs. (11) and (12) include threshold logarithm structures of the form ln(1 − x)/(1 − x). In the qPDF case, the matching coefficient C matches the factorization scale µ to (1−ξ) 1/2 p z , which becomes IR as ξ → 1, just like the DIS and DY cross sections. This can be understood as the probability for a collinear quark to emit a soft gluon with momentum fraction 1 − ξ, which is sensitive to IR physics. In contrast, the matching coefficient C for the spatial correlator and pPDF matches µ to (1 − ω) −1 z −1 0 with z 0 = |z|e γ E /2, which becomes the UV fixed point in the limit of ω → 1. Since the variable ω does not have a direct physical interpretation as momentum fraction, the physics of this UV fixed point is unclear at first sight. As we will see later, the resummation of such logarithms will be very different from the DIS and DY cases. Therefore, it is important to understand the physical origin of the threshold logarithms for qPDF and pPDF.
To understand their origins, let us consider a 3D momentum distribution [31,49], where b µ = (0, b ⊥ , z). We suppress the µ-dependence as the following results apply to both bare and renormalized quantities. Here,q(x, k ⊥ , P z ) is different from the standard transverse-momentum-dependent (TMD) distribution as it is defined with a straight gauge link, instead of the staple-shaped one in the semi-inclusive DIS or DY processes [50,51]. This new straight-gauge-link TMD can also be Fourier transformed to the transverse coordinate space as, Therefore, the 3D momentum distribution can be related to the qPDFq(x, P z ) through which trivially satisfies the definition in Eq. (5). Therefore, the variable x in the qPDF has the physical meaning as a momentum fraction, and the x → 1 limit corresponds to ξ → 1 in Eq. (11), which stands for the emission of a soft gluon from a collinear quark. Therefore, the physical origin of the threshold logarithms in the qPDF is identical to the usual case. As for the spatial correlatorh γ t (λ, −b 2 ) in Eq. (6), we first note that in the infinite momentum limit, where for any finite λ = zP z , the P z → ∞ limit leads to z → 0 so that b 2 → − b 2 ⊥ is completely decoupled from λ. Moreover, due to Lorentz invaraince,h γ t (λ, b 2 ⊥ ) can also be obtained with b µ = (0, b − , b ⊥ ) in light-cone coordinates, which is exactly the correlator that defines the primordial TMD F(x, k 2 ⊥ ) introduced in Ref. [32], where x is restricted to the physical region [−1, 1]. Moreover, the pPDF P(x, b 2 ⊥ ) is related to the primordial TMD through a Fourier transform in the transverse plane [32], Therefore, we can identifyq(x, k ⊥ , P z = ∞) as the primordial TMD andq(x, b ⊥ , P z = ∞) as the pPDF, thus the variable x always has a physical interpretation as the longitudinal momentum fraction. Finally, if we set b 2 ⊥ = z 2 , thenh γ t (λ, b 2 ⊥ ) coincides with the spatial correlator in Eq. (6). Therefore, the variable ω in the matching coefficient C(ω, z 2 µ 2 ) in Eqs. (8) and (10) indeed corresponds to the longitudinal momentum fraction, which lays the physical ground to explain the origin of the threshold logarithms. Now we use a one-loop calculation of the 3D momentum distribution to explicitly demonstrate the emergence of threshold logarithms of (1 − x) in the qPDF and pPDF through relations in Eqs. (16)(17)(18). For an on-shell massless external quark state with momentum p µ = (p z , 0, 0, p z ), the leading collinear and soft divergence of the one-loop 3D momentum distribution comes from the term To obtain the qPDF, we integrate over k ⊥ and find that where the collinear divergence is regulated by , and the limit x → 1 is approached from below as x ∈ (−∞, ∞) for the qPDF [44]. By expanding in , we can reproduce the leading threshold logarithm in Eq. (11). Here the factor (1 − x) −2 is crucial to reproduce the correct sign of leading threshold logarithm, which plays a similar role as the phase-space measure in DIS and DY cross sections [52].
In the hard sub-process in DIS and DY, consider a single emission from the outgoing particle, then the twoparticle phase-space measure with incoming momentum P and outgoing momenta Q and k in d dimensions is (1), and Q 2 = 0 corresponds to the emission of a massless particle, i.e., gluon; for DY-like emissions, , and P 2 is equal to the threshold value. The exponent − or −2 is responsible for the IR logarithms in the threshold limit.
To obtain the pPDF or primordial TMD, we take the limit p z → ∞ in Eq. (21), then the leading divergence in the x → 1 limit is where the exponent 2 has the opposite sign to the qPDF case. After Fourier transforming to b ⊥ space, we obtain lim x→1q (1) Thus the leading threshold logarithm in Eq. (12) is reproduced after expansion in . Since the factorization for the pPDF or spatial correlator works in the small b ⊥ limit, the physical scale in the threshold logarithm is proportional to (1 − x) −2 b −2 ⊥ , which approaches the UV fixed point in the x → 1 limit.
The different dynamical behaviors of the threshold logarithms in the qPDF and pPDF can be understood in this way: in the pPDF case, since it corresponds to the primordial TMD, the emitted gluon remains off-shell with virtual mass −k 2 ⊥ in the limit of x → 1. In coordinate space, small b ⊥ corresponds to large k ⊥ , so the gluon is in the UV region. Moreover, since the limit p z → ∞ has been taken first, only collinear and soft emissions with k ⊥ → 0 are allowed, and the emission of a gluon with finite k ⊥ is suppressed, which explains the factor (1 − x) 2 in analogy to the phase-space measure. However, in the qPDF case, since k ⊥ is integrated over, the limit x → 1 includes contributions from both hard and soft transverse momentum modes, with the latter being sensitive to IR physics. If we first take the soft transverse momentum limit of Eq. (21) at finite p z and x, then k ⊥ is bounded by both (1 − x)p z and xp z . Therefore, it is the integration over k ⊥ that changes the factor 2(1 − x) 2 to 1 + (1 − x) −2 in Eq. (22), which agrees with the fact that the probability of soft gluon emission is divergent. At the end, we note that there is a related comparison between k T -resummation of the standard TMD and threshold resummation of the PDF in Refs. [53,54].
To make more accurate comparisons, let us calculate the full one-loop expression ofq(x, k ⊥ , p z ). We first write down the spatial correlator derived in Ref. [36] where the leading-order (LO) splitting function is given as Then we use Lorentz invariance to generalize Eq. (26) to the case when b µ = (0, b ⊥ , z) in Eq. (14) to obtaiñ q(x, k ⊥ , P z ). By using the formula for the Fourier transform of the logarithmic term, As one can see, the threshold logarithms are evident at finite b ⊥ and p z in Eq. (29). The leading threshold logarithms come from the overlap of collinear and soft divergences of the 3D momentum distribution as shown in Eqs. (24) and (25). The convolution integral in the last line of Eq. (29) can affect the threshold logarithms depending on the limits in Eqs. (17) and (18).
First, we have Then, by explicitly carrying out the convolution integral, we obtain x < 0 plus its virtual part proportional to δ(1 − x). Plugging the above result into Eq. (29), we exactly reproduce the qPDF at one-loop with its threhsold logarithms [36].
On the other hand, we find that Therefore, as long as the above convolution integral vanishes as exponentially suppressed contributions. As a result, the 3D momentum distribution is restricted to the physical region 0 < x < 1, and we reproduce the threshold logarithms in the pPDF or spatial correlator in Eq. (6). According to Eq. (18), we can identify that the limits in Eq. (34) correspond to for the spatial correlatorh γ t (λ = zp z , −b 2 = z 2 ), which implies that the threshold logarithms are sensitive to the long-range correlations.
Last but not least, we note that the collinear divergence in the 3D momentum distribution in Eq. (29) is the same as that in the PDF. Therefore, if b ⊥ Λ −1 QCD , the 3D momentum distribution should also satisfy a factorization formula that goes as where the one-loop matching coefficient C can be obtained by subtracting the 1/ poles from Eq. (29).
Equivalently, the reproduction of the threshold logarithms in the appropriate limits in Eqs. (17) and (18) correspond to the reduction of the matching coefficient for the 3D momentum distribution to those for the qPDF and pPDF in the corresponding limits, Therefore, this 3D momentum distribution can be used as an alternative lattice observable to extract the PDF. Since both p and b can be of any spatial direction, the data will cover a wider range of (p · b, −b 2 ) space, and we can choose p to be along the diagonal directions on the lattice to reach larger boost momentum. Nevertheless, to increase the coverage of (p · b, −b 2 ) space, we also have to let b be off-axis, and thus cannot construct a smooth Wilson line along the b direction, which may complicate the renormalization and introduce discretization effects that break the spatial rotational symmetry. In Ref. [49], a step-like Wilson line geometry was considered, and the study there showed that the Wilson line operator is still multiplicatively renormalizable and the rotational symmetry is weakly broken with smeared gauge links.

A. NLL Resummation Formula
Albeit their differences from the DIS and DY cases, the previous section demonstrated that the threshold logarithms in qPDF and pPDF indeed arise from the leading divergences of the 3D momentum distribution in the x → 1 limit. Now, we turn our attention to the resummation of these potentially large logarithms. The threshold resummation is most conveniently performed in the Mellin-moment space. Let us work on the OPE of the MS spatial correlator at short distance, where we write b 2 = −z 2 for b µ = (0, 0, 0, z), where a N (µ) is the (N + 1)-th moment of the PDF, a N (µ) = 1 −1 dy y N q(y, µ) .
At NLO, the Wilson coefficient C N is given by The threshold limit corresponds to N → ∞, given as where N = N e γ E . As one can see, the leading and subleading threshold logarithms are, respectively, α s ln 2 N , α s ln N .
We propose to use the standard technique of threshold resummation to resum these large logarithms [41]. Usually, it requires resummation at NLL accuracy to guarantee the convergence when one performs the inverse Mellin transform to obtain the resummed matching coefficient C or C [41]. From now on let us denote a s = α s /(2π). The all-order resummed form of the Wilson coefficient C N inspired by the standard threshold resummation technique is where a = −2 as it corresponds to the logarithm ln((1 − w) 2 z 2 0 µ 2 ) in Eq. (12). As a consequence of the negative a, k 2 is evolved to the UV fixed point in the limit x → 1. In comparison, a = 2 and 1 correspond to threshold resummation for DIS and DY cases, respectively [55], which reaches the Landau pole. Additionally, Since the leading threshold logarithm arises from the overlap of collinear and soft divergences in the 3D momentum distribution, which is the common origin of Sudakov double logarithms [56,57], the coefficient A should be the universal cusp anomalous dimension [58,59]. Based on the one-loop Wilson coefficient in Eq. (42), For NLL resummation which neglects O(α 2 s ln N ) terms, we only need A (1) , the two-loop cusp anomalous dimension [58], where C A = 3, T F = 1/2, and n f is the number of active quark flavors. According to the NNLO matching coefficient C calculated in Ref. [42], which exactly verifies Eq. (48) and the expansion of Eq. (44) to NNLO, especially that a = −2.
Plugging the above results into Eq. (44), we obtain the resummed Wilson coefficient, where τ = β 0 a s ln N , L = ln(z 2 0 µ 2 ). The first term is of O(a s ), so it can be either exponentiated or kept in the NLO matching coefficient for NLL resummation. Here, we exponentiate it in accordance to the treatment in the ASV analysis [12]. The functions g 1 and g 2 are Since τ ∼ 1 at NLL accuracy and L is not resummed at the fixed order, Eqs. (51) and (52) correspond to the hierarchy L ln N . For practical applications to lattice QCD calculations, the range of available z is not necessarily small, so one might also need to include the DGLAP evolution [60][61][62] of the PDF or Mellin moments to resum L. To accomplish this, we start from the evolution equation for the Wilson coefficient C N , where γ N is the anomalous dimension which expands in a s , At one-loop [36], Therefore, the solution to Eq. (54) at leading logarithmic (LL) accuracy is where the superscript "evo" indicates the DGLAPevolved Wilson coefficient. Then, we can perform threshold resummation of the Wilson coefficient at µ = z −1 0 , C N α s (z −1 0 ), 1 , which is related to Eq. (50) by setting L = 0, where Using the inverse Mellin transform, we can eventually obtain the resummed matching coefficient with LL DGLAP evolution as which can be carried out numerically. The resummed matching coefficient C NLL+evo (ξ, µ/p z ) for the qPDF is related to C NLL+evo (w, z 2 µ 2 ) through a double Fourier transform [36], which becomes more involved and is beyond the scope of this work. Notably, both g 1 (τ, 0) and g 2 (τ, 0) are regular for all N . Therefore, unlike DIS and DY, there is no Landau pole along the path of integration in Eq. (61), which is due to the UV fixed point that we mentioned above.
In the previous work by some of the authors [43], the lattice data has been analyzed with the NLO Wilson coefficients C NLO N (α s (µ), z 2 0 µ 2 ) in Eq. (39). Now we can include the threshold resummation at NLL and DGLAP evolution by replacing C N in Eq. (39) with C evo N from Eq. (57) then further replacing C N at µ = 1/z 0 in Eq. (57) with C NLL N in Eq. (58). We can also combine the fixed-order and threshold resummation for all N by replacing C N with where C NLOexp , which is the same as the right-hand side of Eq. (42), so that at small N the latter reduces to C NLO N up to O(a 2 s ln 2 N , a 2 s ln N , a 2 s ) corrections. Note NLL threshold resummation is accompanied by the NLO Wilson coefficient.

B. Estimating impacts of NLL resummation on the large-x exponent
Since the resummed Wilson coefficient is generally greater than the fixed-order one at large N , the fitted higher moments with threshold resummation should be smaller, which will suppress the PDF at large x or lead to a larger value of β.
For a more quantitative estimate of the effect of threshold resummation on β, we choose a PDF parametrized as the form where N is a normalization constant and G(x) is a smooth and well-behaved function within 0 < x < 1. It was shown [43] that at large N regardless of the function G(x). If x N is known, then β can be approximated as For given lattice matrix elements, if one solely changes the matching coefficient from C NLO N to C NLO+NLL N , then β will be changed by  For the purpose of this discussion we do not include DGLAP evolution, that is, we set µ = 1/z 0 so that ln(z 2 0 µ 2 ) = 0 in C NLO+NLL N and C NLO N . In Fig. 1, we plot δβ(N ) as a function of N for µ = 3.2 GeV and α s (µ) = 0.242, within the range of N ∈ [1,16]. We have to truncate at N = 16 because C NLO N would become negative at N > 16. At N = 4, threshold resummation increases the NLO fit of β by about 0.2, which is a weak effect. As N increases the impact of threshold resummation becomes stronger. The increment of β by threshold resummation was also observed in the ASV fit [12] of the DY data. Since C NLO N falls close to zero as N increases, it will significantly enhance the fitted higher moments and lead to smaller and even negative β, which could be compensated by the rapid growth of δβ shown in Fig. 1. This indicates that threshold resummation is necessary and may help stabilize the fitting of β.

IV. APPLICATION TO LATTICE QCD CALCULATIONS
In this section, we study the impact of threshold resummation on the valence PDF of pion based on the recent lattice QCD results of Ref. [43]. Since our resum-mation has been derived in the Mellin-moment space, it is naturally implemented on the coordinate-space lattice data. Using the OPE formula in Eq. (39) with C NLO N , C NLO+NLL N with DGLAP evolution, and C NLOevo N , the coefficient with only DGLAP evolution included, we study their effects by fitting the data with a finite number of moments and also a parameterization model for the pion valence PDF.
The lattice results we use are from a previous publication by some of the authors [43]. These results are for a pion wtih 300 MeV mass, including two different lattices with spacing a = 0.04 fm and 0.06 fm, pion momenta P z = (2πn z )/L s with n z from 0 to 5 listed in Table I, and volumes (2.56) 3 fm 3 and (2.88) 3 fm 3 for each a respectively. Therefore, the largest boost momentum on the finer lattice is P z = 2.42 GeV.

A. Model-independent moments
In the OPE of the spatial correlator, the higher-twist contributions and higher moments contributions are suppressed by O(z 2 Λ 2 QCD ) and O(λ N /N !), respectively. The former restricts z to be short, which in turn limits the range of λ = zP z with finite P z . While the higher moments a N are further suppressed by the asymptotic behavior of PDF at large x, the currently available lattice results are only sensitive to the lowest ones. In Ref. [43] it was found that these lattice results are not sensitive to moments beyond N = 6. Furthermore, since the pion matrix elements for the valence PDF are real and symmetric in z, only even moments contribute. Thus, for this OPE-type analysis, reasonable signals are expected for only x 2 and x 4 .
The lattice data we fit to are renormalization-groupinvariant ratios constructed from bare matrix elements h γ t (λ = zP z , z 2 ) [43] with different pion momenta, which was first proposed in Ref. [33] with P z 0 = 0 and then adapted in Ref. [63] to use P z 0 Λ QCD to suppress higher-twist effects. For this analysis, we choose aP z 0 = 2πn 0 z /L s with n 0 z = 1, 2 for both the lattice spacings. The leading-twist approximation to the above ratio is where c N (z 2 µ 2 ) is the MS Wilson coefficient C N normalized by C 0 without changing the outcome. We consider c N for NLO without any DGLAP evolution (denoted as NLO), NLO with LL DGLAP evolution (denoted as NLOevo), and NLO+NLL with LL DGLAP evolution (denoted as NLOevo+NLL), respectively. We use one-loop QCD β function for the DGLAP evolution with n f = 3 to be consistent with the 2+1 flavor lattice ensemble. Here, x N P z is the Mellin moment with target mass correction included [46], The target-mass corrections have little effect for the pion because of the smallness of its mass m π = 300 MeV in our case. In the previous work [43], the data has been analyzed with only the NLO Wilson coefficients. In Fig. 2, we show c N (z 2 µ 2 ), with N = 2, 4, 16 and at a scale µ = 3.2 GeV, as a function of 1/z 0 for the leading order (LO) in α s , NLO, NLOevo and NLOevo+NLL. The DGLAP evolution is included here because for a large range of z, the logarithm L = ln(z 2 0 µ 2 ) can become large and reduce the predictive power of OPE with fixed-order Wilson coefficients. As one can see, the DGLAP evolution is indeed an important effect within the range of z we considered, and the LL resummation of L makes the Wilson coefficients change slower than the NLO ones as 1/z 0 is varied. However, one can also notice that as N increases the size of one-loop correction grows in the NLOevo Wilson coefficient, which even makes it turn negative at small scales of 1/z 0 . This is exactly caused by the large logarithmic term −α s ln 2 N in Eq. (42), especially with large α s (z −1 0 ), and it implies the necessity of threshold resummation for either large α s or large N appearing in the Wilson coefficients. As shown in Fig. 2, significant impact can be observed for large α s (small 1/z 0 ) or N when we compare the NLOevo+NLL coefficients to the NLOevo ones, though there is only mild difference for N = 2, 4 at 1/z 0 > 2 GeV.
As mentioned above, current lattice data is only sensitive to the first few lowest moments. In Fig. 3, we take the a = 0.04 fm lattice data as an example, and plot the fitted second moment x 2 and fourth moment x 4 at each single 1/z 0 . In the plots, we skip the data point at z = a (or 1/z 0 =5.5 GeV) to avoid the possible lattice discretization effects. In principle, the result x n (µ) should only depend on the MS scale µ, as the ln(z 2 )dependence of lattice data should be canceled by that of the Wilson coefficient. The cancellation is expected to work better at smaller z or at larger 1/z 0 since α s (1/z 0 ) is smaller. Focusing on the second moment x 2 , at LO one can clearly observe the z-dependence of x 2 LO (z), which grows faster at smaller scales in 1/z 0 . Beyond LO, one can find a plateau for x 2 which indicates that the coefficients c 2 (z 2 µ 2 ) can explain the evolution of the moments well in the plateau region. The fixed-order analysis can be found in Ref. [43]. After introducing the DGLAP evolution, one can observe the plateau for x 2 NLOevo when 1/z 0 0.8 GeV. Beyond 1/z 0 ∼ 0.8 GeV, the plateau remains stable up to around 1/z 0 ∼ 3 GeV.
For x 2 NLOevo+NLL one can see that threshold resummation slightly improves predictive power, as the plateau region is extended to 1/z 0 0.6 GeV within the errors. Nevertheless, the plateau value remains almost unchanged compared to x 2 NLOevo as we expect. The lower scale region is not showing a good plateau but is not of concern because α s (1/z 0 ) becomes large and is accompanied by large theoretical uncertainty. As for the fourth moment x 4 (z) extracted from each z, though noisy, one can observe the problematic behaivor of the x 4 NLOevo at small scale is caused by large −α s ln 2 N in the Wilson coefficients as discussed above, and the situation is improved for x 4 NLOevo+NLL with NLL resummation within the errors.
To stabilize the fit and extract higher moments, we then use NLOevo+NLL Wilson coefficients to perform a combined fit to all z within the range [z min , z max ]. We include the theoretical uncertainty in the fit by varying the factorization scale µ up and down by a factor of two as in Ref. [43]. The fitted second and fourth moments are shown in Fig. 4 as a function of z max with n 0 z = 1, 2 and z min = 2a, 3a for a = 0.04 fm (up) and 0.06 fm (down). We skipped z min = a to avoid possible discretization effects. We observed that the z min as well as the n 0 z dependence are mild within the errors. Therefore, we estimate the statistic and systematic errors as the darker and lighter bands using the same strategy as in Ref. [43]: for each bootstrap sample, we evaluate the average and standard deviation of observable A from a range of n 0 z , z min and z max as Mean(A) and SD(A). Then we can get the statistic errors of Mean(A), and take SD(A) as the systematic error. In this analysis, we use n 0 z = 1,2, z min = 2a, 3a and z max ∈ [0.36, 0.48] fm. One can observe a plateau for both x 2 and x 4 within the errors.
Since mild lattice spacings dependence of the moments can be observed in Fig. 4, a joint fit is performed for a continuum extrapolation by modifying the moments inserted in Eq. (68) as The continuum extrapolations are shown in Fig. 5. One observes that the second moment indeed shows some lattice spacing dependence, while the fourth moment shows little dependence within errors. For comparison, we also show the moments of our previous fixed-order estimate, and the xFitter [18], JAM [17] ASV [12] pion valence PDF, which are extracted from fits with data from experiments. Though the mean value of the second moment x 2 = 0.103(11)(3) is slightly higher than the NLO fit 0.099(7)(5) in our previous work [43], our result still shows better agreement with xFitter and JAM. The fourth moment x 4 shows overall consistency with the NLO fit [43] and the results extracted from experimental measurements because of the large errors. This mild change of the first few moments compared to the NLO analysis is mainly driven by the DGLAP evolution, while the threshold resummation has little impact but helps to stabilize the fit in the low scale region of 1/z 0 . To access higher moments, where the threshold resummation generally has a considerable impact, one has to increase the pion momentum while working at a small z. The fit results of the moments are summarized in Table II [43], xFitter [18], JAM [17] and ASV [12] valence PDF of pion are also shown as the lines.
of λ = zP z , so it is difficult for us to control the systematic uncertainties from model dependence like the global analyses. Nevertheless, we can still explore several models that are encoded with the power-law behavior (1−x) β , and see how β changes under threshold resummation. As we did previously in the NLO fit in Ref. [43], we use a simple ansatz for the pion valence PDF: where N (α, β) = Γ (2 + α + β) /(Γ(1 + α)Γ(1 + β)) is chosen such that 1 0 dxf π v (x) = 1. We construct moments from the models, which are determined by the parameters in the ansatz, and plug them into Eq. (68). Then we fit the lattice data by minimizing , where σ 2 stat (z, P z , P z 0 ) are the statistic errors and σ 2 sys (z, P z , P z 0 ) are the systematic errors defined by scale variation, More details of a similar procedure can be found in Ref. [43]. Though the models are used to fit the lattice data, the corresponding parameters are mainly determined by the first few moments. Therefore, we truncate the series in Eq. (68) at N = 16 for the fit, which is large enough for the result to stabilize. Similar to the moment fit, we use  6. Valence PDF of pion from NLOevo+NLL fits of lattice results [43] to Eq.refeq:model. Also shown are the same PDF from NLO fits [43], xFitter [18], JAM [17] and ASV [12].
results are shown in Table II. As one can see, the largex exponent β in this NLO+NLL analysis with DGLAP evolution (NLOevo+NLL) is around 1.0 with about 50% uncertainty. The central value of β is increased from our previous fix-order NLO analysis [43] by a small amount, which should be mainly driven by the DGLAP evolution, but they are still consistent within errors. This is due to the same reason that our current lattice data are only sensitive to the first few moments, which endures a small impact from threshold resummation. The situation can be improved if we manage to get more precise data and increase the pion momentum. One may notice that the results we get here have larger errors than our previous NLO fit [43], which is because we used more data with z max up to 0.6 fm in that analysis. Finally, in Fig. 6 we show the comparison of pion valence PDF by xFitter [18], JAM [17], ASV [12] as well as the fix-order estimate from the same lattice data [43].

V. CONCLUSION
In this work, we studied the threshold resummation in the qPDF, the pPDF and the corresponding spatial correlator, which have been used for the lattice calculation of PDFs in the LaMET and pPDF approaches. We find that the threshold logarithms in them originate from the interplay between the collinear and soft divergences in a 3D momentum distribution, and they are sensitive to long-range correlations of the spatial correlator in coordinate space. Moreover, the 3D momentum distribution, when Fourier transformed from k ⊥ to b ⊥ space, reduces to the qPDF and the pPDF in the b ⊥ → 0 and infinite momentum limits, respectively, which offers us new insights on the relation between the latter two quantities.
Motivated by the standard procedure, we propose to resum these threshold logarithms in the OPE of the spa-  (52)(15) TABLE II. Summary of results for valence PDF of pion at the scale µ = 3.2 GeV, obtained using NLO+NLL matching with DGLAP evolution (denoted as NLOevo+NLL) on the lattice QCD results of Ref. [43]. We show results for two different lattice spacings a as well as the continuum estimate denoted by a → 0. For comparison, the 2-parameter fit results with NLO matching and without DGLAP evolution are α = -0.55(15)(08) and β = 0.66 (34) (22) in Ref. [43]. tial correlator at NLL accuracy, and confirm our results by comparing them with the NNLO Wilson coefficients. Unlike DIS and DY, the resummed Wilson coefficients in this case match the factorization scale to the UV fixed point in the threshold limit, and therefore are free from the influence of Landau poles. In addition, we included DGLAP evolution for applications to more realistic lattice calculations where the range of available z is not small. We showed that, apart from the DGLAP evolution, our proposed resummation of threshold logarithms leads to the suppression of PDFs at large x, similar to the DIS and DY cases as observed in the ASV analysis [12].
We then reanalyzed the pion valence PDF using recent lattice data from [43] by applying our proposed resummation formula at NLO+NLL accuracy with LL DGLAP evolution. We found that threshold resummation makes little impact on the second and fourth Mellin moments, as well as on the large-x exponent of the PDF from model fits, compared to the previous fixed-order analysis [43] with NLO matching. This is because the currently available lattice data are sensitive only to the lowest moments or moderate values of x, and they are not precise enough. On the other hand, the effects of DGLAP evolution in the NLO matching constitutes an important systematic uncertainty that should be included in the analysis. Finally, with increased pion momentum and data precision in the future, we expect threshold resummation to become a necessary systematic correction in the calculation of higher moments and the large-x PDF.