Lattice QCD Determination of the Bjorken-$x$ Dependence of Parton Distribution Functions at Next-to-next-to-leading Order

We report the first lattice QCD calculation of pion valence quark distribution with next-to-next-to-leading order perturbative matching correction, which is done using two fine lattices with spacings $a=0.04$ fm and $0.06$ fm and valence pion mass $m_\pi=300$ MeV, at boost momentum as large as $2.42$ GeV. As a crucial step to control the systematics, we renormalize the pion valence quasi distribution in the recently proposed hybrid scheme, which features a Wilson-line mass subtraction at large distances in coordinate space, and develop a procedure to match it to the $\overline{\rm MS}$ scheme. We demonstrate that the renormalization and the perturbative matching in Bjorken-$x$ space yield a reliable determination of the valence quark distribution for $0.03\lesssim x \lesssim 0.80$ with 5-20\% uncertainties.

Understanding the hadron inner structure remains one of the top fundamental questions in nuclear and particle physics. As the lightest hadrons in nature, pions are the Nambu-Goldstone bosons of quantum chromodynamics (QCD), and their quark and gluon structures can help understand the origins of hadron mass and dynamical chiral symmetry breaking. The parton distribution functions (PDFs), which describe 1D momentum densities of quarks and gluons in a hadron, are the simplest and most important quantities that have been extensively studied from global high-energy scattering experiments and will be probed at unprecedented precision at the future Electron-Ion Collider [1,2]. Besides the experimental efforts, the first-principles calculations of PDFs using lattice QCD are also expected to provide useful predictions.
Computation of the PDFs on a Euclidean lattice has been extremely difficult because they are defined from light-cone correlations with real-time dependence in Minkowski space. For a long time, only the lowest moments of the PDFs were calculable as they are matrix elements of local gauge-invariant operators. For reviews see Refs. [3,4]. Less than a decade ago, a breakthrough was made by large-momentum effective theory (LaMET) [5][6][7], which starts from a Euclidean "quasi-PDF" (qPDF) in a boosted hadron and obtains the PDF through a largemomentum expansion and perturbative matching of the qPDF in Bjorken-x (longitudinal momentum fraction) space. Over the years, LaMET has led to much progress in the calculation of PDFs and other parton physics [4,7], which reinvigorated the field as other proposals [8][9][10][11][12][13] are also being studied and implemented.
Despite substantial progress, lattice calculation of the PDF x-dependence has yet to achieve essential control * xgao@bnl.gov † yong.zhao@anl.gov of the systematic uncertainties [14]. In the LaMET approach, lattice renormalization is one of the most important sources of error. The nonlocal quark bilinear operator O Γ (z) ≡ψ(z)ΓW (z, 0)ψ(0), where Γ is a Dirac matrix and z µ = (0, 0, 0, z), which defines the qPDF, suffers from a linear power divergence in the Wilson line W (z, 0) that must be subtracted before taking the continuum limit. The most popular methods so far are the regularization independent momentum subtraction scheme [15][16][17][18] and other ratio schemes [19][20][21][22], which use the matrix element of O Γ (z) in an off-shell quark [15][16][17][18], a static/boosted hadron [19,22] or the vacuum state [20,21] as the renormalization factor. At small z the matrix elements in these schemes satisfy a factorization relation to the light-cone correlation [13,[23][24][25]. However, at large z they introduce nonperturbative effects [26] that propagate to the qPDF via Fourier transform (FT) of the matrix elements, which contaminates the LaMET matching in x-space. To overcome this limitation, the hybrid scheme [27] was proposed to subtract the linear divergence at large z and match the result to the MS scheme, thus preserving the LaMET matching after FT. To date, the hybrid scheme has not been used in calculating the PDFs, except for a recent work on meson distribution amplitudes [28]. Apart from renormalization, the accuracy of perturbative matching also controls the precision of the calculation. In all the existing lattice calculations, the matching was done at only next-to-leading order (NLO), and it is not until recently that the next-to-next-to-leading order (NNLO) matching was derived for the non-singlet quark qPDF in the MS scheme [21,29].
In this Letter we present a state-of-the-art calculation of pion valence quark PDF using high-statistics, superfine-spacing, and large-momentum lattice data [26], with an adapted hybrid-scheme renormalization and the first-time implementation of NNLO matching. The pion valence PDF has been extracted from global fits [30][31][32][33] and studied in lattice QCD [26,[34][35][36][37][38][39][40], with both at NLO accuracy. In this work, we subtract the linear divergence in O Γ (z) with sub-percent precision, and develop a procedure to match the lattice subtraction scheme to MS, a crucial step in the hybrid scheme to reduce the power corrections [27]. We derive the NNLO hybrid-scheme matching and apply it to the qPDF, showing good perturbative convergence and reduced scale-variation uncertainty compared to NLO matching. Finally, we demonstrate that our analysis yields a reliable determination of the PDF for 0.03 x 0.80 with 5-20% uncertainties.
Our lattice data was produced using gauge ensembles in 2+1 flavor QCD generated by the HotQCD collaboration [41] with Highly Improved Staggered Quarks [42], including two lattice spacings a = 0.04 and 0.06 fm, and volumes L 3 s ×L t = 64 4 and 48 3 ×64, respectively. We use tadpole-improved clover Wilson valence fermions on the hypercubic (HYP) smeared [43] gauge background, with a valence pion mass m π = 300 MeV. Furthermore, the Wilson line in O Γ (z) is constructed from HYP-smeared gauge links. We use pion momenta P z = (2πn z )/(L s a) with 0 ≤ n z ≤ 5, resulting in P z as large as 2.42 GeV.
The qPDFf v (x, P z , µ) is defined in a boosted pion state |P with four-momentum P µ = (P t , 0, 0, P z ): whereh(z, P z , µ) ≡ P |O γ t (z)|P /(2P t ), and µ is the MS scale. The operator O Γ (z) can be renormalized under lattice regularization as [44][45][46] where "B" and "R" denote bare and renormalized quantities. The factor Z O (a) includes all the logarithmic ultraviolet (UV) divergences which are independent of z, while the Wilson-line mass correction δm(a) includes the linear UV divergence ∝ 1/a and can be expressed as where m −1 (a) is a series in the strong coupling α s (1/a), and m 0 is an O(Λ QCD ) constant originating from the renormalon ambiguity in m −1 (a) [47]. The hybrid scheme is implemented as follows: For 0 ≤ z ≤ z S with a z S 1/Λ QCD , we form the ratiõ h(z, P z , a)/h(z, 0, a) to cancel the UV divergences and the cutoff effects from z ∼ a [19]; at z > z S we subtract δm(a) and determine Z O (a) by imposing a continuity condition of the renormalized matrix elements at z = z S . There are different ways to calculate δm(a) [27,46,[48][49][50][51]. We determine δm(a) from the combination of the static quark-antiquark potential, V lat (r) [41,52], and the free energy of a static quark at non-zero temperature [53][54][55], with the following normalization scheme, V lat (a, r = r 0 ) + 2δm(a) = 0.95/r 0 , where r 0 = 0.469 fm is the Sommer scale for 2+1 flavor QCD [41], and the constant 0.95 defines the scheme. The linear divergence m −1 (a)/a does not depend on the scheme, while m 0 does. The results are aδm = 0.1586 (8) and 0.1508 (12) for a = 0.06 and 0.04 fm, respectively. Since m 0 is scheme dependent, a factor of e m0|z| with m 0 ∼ O(Λ QCD ) is needed to match the lattice scheme to MS, otherwise the LaMET expansion of the qPDF will include a power correction ∝m 0 /P z [27], which slows down convergence to the PDF as P z grows. It was proposed that m 0 can be obtained by comparing the subtracted matrix elements of O Γ (z) [51] or W (z, 0) [49] with their MS operator product expansion (OPE), whose accuracy requires z 0.2 fm [27]. But due to discretization effects, the window of z that can be used is actually narrow.
Our new procedure for the hybrid scheme is distinct by the determination of m 0 . In order to use larger z, we construct the following ratio and compare it to a form motivated by the OPE ofh(z, 0, µ), where z, z 0 a, and the parameter Λ ∼ O(Λ 2 QCD ). The Wilson coefficient C 0 is known to NNLO [21,25,29], and m 0 and Λz 2 originate from the leading UV and infrared renormalons in C 0 [20]. According to Eq. (2), the l.h.s. of Eq. (5) must have a continuum limit if δm(a) includes all the linear divergences, which is renormalization group (RG) invariant. We choose z ≥ z 0 = 0.24 fm and find agreement between the a = 0.04 fm and a = 0.06 fm ratios at sub-percent level up to z ∼ 1 fm (see App. A). Then we extrapolate the lattice ratios to the continuum with a 2 -dependence [26], and fit the result to the r.h.s of Eq. (5). For z 0 ≤ z ≤ 0.4 fm, we obtain decent plateaus and χ 2 values for bothm 0 and Λ with the NNLO C 0 . By definitionm 0 cancels the lattice scheme dependence of δm(a), as changing the scheme only shifts δm(a) by a constant, butm 0 will inherit the O(Λ QCD ) ambiguity in the MS scheme. Since C 0 is at fixed order, bothm 0 and Λ depend on µ, which we vary to estimate the related uncertainty in the final result. At µ = 2.0 GeV,m 0 = 0.151(1) GeV and Λ = 0.041(6) GeV 2 , so the power correction is not negligible. Therefore, we modify the hybrid scheme by correcting the Λz 2 term inh(z, 0, µ) at short z as where δm = δm +m 0 , and N =h(0, 0, a)/h(0, P z , a) normalizesh(z, z S , P z , µ, a) to one at z = 0. Since C 0 is at fixed order,h(z, z S , P z , µ, a) depends on µ despite the fact that it should be RG invariant. Such a renormalization is performed through bootstrap loops so that the correlation between different P z and z is taken care of. The hybrid-scheme matrix elements are shown in Fig. 1. At small z,h(z, P z ) is dominated by the leadingtwist contribution. At large z, the spacelike correlator for pion valence quarks will exhibit an exponential decay ∝ e −m eff |z| where m eff is an effective mass related to the system [56]. When plotted as a function of λ = zP z , h(λ, P z ) should scale in P z at small λ, with slight violation due to QCD evolution. Its exponential decay will emerge at a larger λ with greater P z and with decay rate m eff /P z . In the P z → ∞ limit, the exponential decay vanishes at finite λ (z → 0), and only the leading-twist contribution remains, which almost scales in P z and features a power-law decay at large λ that corresponds to small-x PDF [27]. This picture is consistent with Fig. 1.
The next step is a FT. We truncate the matrix elements at z L or λ L = z L P z whereh(λ L ) ∼ 0, and extrapolate to ∞ to remove the unphysical oscillations from a truncated FT [27]. The extrapolation form is Ae −m eff |z| /|λ| d , where A, m eff and d are the parameters. Since m eff is independent of P z , by fitting to the P z = 0 matrix elements we find that it is around 0.1 GeV, which is not far from the phenomenological estimate of 0.2-0.5 GeV in HQET [57]. Therefore, we impose m eff > 0.1 GeV, as well as A > 0 and d > 0, to ensure a convergent FT on each bootstrap sample. Since the FT converges fast with the exponential decay, the extrapolation mainly affects the small-x region apart from removing the unphysical oscillations. To verify this we vary z L , which turns out to have little impact, and use different m eff bounds and extrapolation forms, which lead to consistent qPDFs down to x ∼ 0.05. (See App. B).
Then, we match the qPDFf v (x, λ S , P z , µ) to the MS PDF f v (x, µ) through LaMET [25,27,58,59]: where λ S = z S P z , z S = 0.24 fm, and the power corrections are controlled by the parton and spectator momenta xP z and (1 − x)P z [27]. Here C −1 is the inverse of the hybrid-scheme matching coefficient C, which we derive at NNLO [60] by conversion from the MS result [21,29]. Based on Eq. (7), we can directly calculate the PDF with P z -controlled power corrections for x ∈ [x min , x max ].
In Fig. 2 we show the results of perturbative matching. The matching drives the qPDF to smaller x and reduces the statistical errors at moderate x, because matching effectively relates the qPDF from finite P z to infinity, and the qPDF evolves to smaller x as P z increases. The NNLO correction is generally smaller than the NLO correction, which indicates good perturbative convergence, a crucial criterion for precision calculation. Besides, by varying µ and evolving the matched results to the same µ, we find that the scale-variation uncertainty is reduced at NNLO, which is further evidence of improved precision. The matching correction diverges as x → 0, implying that resummation of small-x logarithms is needed. A resummation is also necessary as x → 1 [40], but these resummations are not needed for moderate x.
We compare the PDFs obtained at different P z with NNLO matching in Fig. 3. At moderate x, the P zdependence is remarkably reduced, and the results appear to converge for P z ≥ 1.45 GeV, which strongly indicates the effectiveness of LaMET matching. At x 1, each PDF curve has a small non-vanishing tail due to the power corrections in Eq. (7), but they decrease with larger P z (see also App. C 3). To estimate the size of the power corrections, we fit the PDFs obtained at a = 0.04 fm, P z = {1.45, 1.94, 2.42} GeV and a = 0.06 fm, P z = {1.72, 2.15} GeV to the ansatz f v (x)+α(x)/P 2 z for each fixed x, where we ignore the a-dependence as it has been found that the matrix elements have O(a 2 P 2 z ) effects that are less than 1% [26]. Since this fit is mainly affected by the data sets at lower P z with smaller statistical errors, which have larger power corrections, we use the result at P z = 2.42 GeV instead of the fitted f v (x) as our final prediction. The power correction at P z = 2.42 GeV is estimated to be α(x)/[P 2 z f v (x)] < 0.10 for 0.01 < x < 0.80. It is surprising that the results are insensitive to P z for x as small as 0.01, nor do they show dependence on the extrapolation form in the FT as we have checked. This can be explained by that, under matching, the qPDF contributes to the PDF at larger x which has less dependence on P z or the extrapolation. Nevertheless, it must be pointed out that the smallness here is only relative, as α(x)/P 2 z still diverges as x → 0. Our final prediction for the pion valence quark PDF (BNL-ANL21) is shown in Fig. 4, which is obtained from the qPDF at a = 0.04 fm, z S = 0.24 fm, z L = 0.92 fm, µ = 2.0 GeV and P z = 2.42 GeV with exponential extrapolation and NNLO matching. The red band represents the statistical error, and the light purple band includes the error from scale variations, which is obtained by repeating the same analysis for µ = 1.4 GeV and 2.8 GeV and evolving the PDFs to µ = 2.0 GeV with the NLO DGLAP kernel. Since the hybrid-scheme parameterm 0 depends on µ, the small scale variation in the final result shows that the renormalization uncertainty is well under control. We require that the O(α 3 s ) matching correction at µ = 2.0 GeV be smaller than 5%, which propagates geometrically to < 37% at NLO and < 14% at NNLO, thus excluding x < 0.03 and x > 0.88. A list of the above uncertainties at selected x is shown in Table I. See also App. D. We neglect the FT uncertainty as it is extremely small. As for m π dependence, our associated calculation of the second PDF moment at m π = 140 MeV [61] shows consistency within 5% statistical uncertainty, which will be validated by a direct comparison in the future. Previous studies [62,63] also suggest that the finite volume correction is less than 1% for our lattice setup. At last, by limiting the estimated power corrections to be less than 10%, we determine the PDF at 0.03 x 0.80 with 5-20% uncertainties. Our result is in great agreement with the recent global fits by xFitter [31] and JAM21nlo [32] for 0.2 < x < 0.6, but deviates from the earlier GRVPI1 [30] and ASV [33] fits. When compared to a previous analysis of the same s ) and power corrections be smaller than 5% and 10%, respectively. lattice data (BNL20) [26], which used a short-distance factorization of the matrix elements at NLO, and a parameterization of the PDF, our new result has shifted central values and considerably reduced uncertainties at moderate x, but still agrees within errors. With finite P z and statistics, lattice QCD can only make predictions for x ∈ [x min , x max ]. The PDF parameterization correlates the information at all x ∈ [0, 1], so the larger uncertainties at moderate x in BNL20 could be propagated from the uncontrolled errors in the end-point regions. Besides, there is no practical estimate of the model uncertainty in the parameterization. Therefore, the LaMET calculation for x ∈ [x min , x max ] is more reliable as it does the power expansion and matching directly in x-space.
In summary, we have performed a state-of-the-art lattice QCD calculation of the x-dependence of pion valence quark PDF, where we developed a procedure to renormalize the qPDF in the hybrid scheme and match it to the MS PDF at NNLO. The final results show reduced perturbation theory uncertainty and converge at moderate x with pion momenta greater than 1.45 GeV, which allows us to reliably estimate the systematic errors. This calculation can be improved with physical pion mass, continuum extrapolation, and higher statistics for the matrix elements at long distances and at larger boost momenta.
Our renormalization procedure can also be incorporated into the lattice calculations of gluon PDFs, distribution amplitudes, generalized parton distributions and transverse momentum distributions. With the systematics under control, we can expect lattice QCD to provide reliable predictions for these quantities in the future. As has been described in the main text, the hybrid scheme renormalization includes two parts: • For z ≤ z S , we form the ratio of bare matrix elements [19],h which has a well-defined continuum limit and is renormalization group (RG) invariant.
• For z ≥ z S , the renormalized matrix element is which is equal to the ratio in Eq. (A1) at z = z S .
To determine δm(a) we use the additive renormalization constant, c Q (a) = δm(a), which is obtained in Ref. [54] from the analysis of the free energy of a static quark, F Q (T ), at non-zero temperature T with the normalization condition in Eq. (4). Recently F Q has been calculated using one step of HYP smearing [55], and it was found that HYP smearing does not affect the temperature dependence of F Q (T ), but only shifts it by an additive constant. Therefore, we have F B,1 with superscripts 0 and 1 referring to the number of HYP smearing steps in the bare free energy of the static quark. Using the lattice results for F B,0 Q (T ) and F B,1 Q (T ) obtained on N τ = 12 lattices and temperatures corresponding to a = 0.04 fm and a = 0.06 fm (where cutoff effects can be neglected), as well as the values of c Q from Table X of Ref. [54] for β = 7.825 (a = 0.04 fm) and β = 7.373 (a = 0.06 fm), we obtain δm(a). The results are aδm(a = 0.06 fm) = 0.1586(8) and aδm(a = 0.04 fm) = 0.1508 (12).
First of all, to test how well the subtraction of δm(a) can remove the linear divergences inh(z, P z , a), we construct the ratio in Eq. (5), where z 0 = 0.24 fm for both lattice spacings. According to Eq. (2), the renormalization factor Z O (a) cancels out in the ratio. Therefore, if δm(a) includes all the linear divergences, thenR(z, z 0 , a) should have a well-defined continuum limit. Our lattice results for the above ratio with z 0 = 0.24 fm is shown in Fig. 5. As one can see, the differences between the ratios at a = 0.04 fm and 0.06 fm are at sub-percent level, which clearly shows that the linear divergences have been sufficiently subtracted by δm(a). Therefore, the ratio in Eq. (A3) has a continuum limit which is RG invariant. Our next step is to match the lattice subtraction scheme to MS. When z, z 0 Λ −1 QCD , the MS matrix elementh MS (z, 0, µ) has an OPE that goes as where m MS 0 is the O(Λ QCD ) renormalon ambiguity from the Wilson line self-energy renormalization [20], O tw4 (µ) is a twist-four operator (for example,ψD 2 ψ or gψσ µν F µν ψ), C 0 and C 2 are perturbative coefficient functions, and ". . ." denotes contributions at higher twists. Since P z = 0, C 0 is the only Wilson coefficient that contributes at leading-twist. The leadingtwist contribution is proportional to P |ψγ t ψ|P /(2P t ) which is trivially one due to vector current conservation. Sinceh MS (z, 0, µ) is multiplicatively renormalizable, both C 0 (z 2 µ 2 ) and C 2 (z 2 µ 2 ) P |O tw4 (µ)|P must satisfy RG equations with the same anomalous dimension, which is known to next-to-next-to-next-to-leading order (N 3 LO) [64]. Due to the ambiguity in summing the perturbative series in C 0 (z 2 µ 2 ), there are O(Λ 2n QCD ) IR renormalons in the leading-twist contribution that should be cancelled by those from higher-twist condensates, along with the O(Λ QCD ) UV renormalon to be cancelled by m MS 0 [20,57]. Both the UV and IR renormalon contributions cannot be well defined unless one specifies how to sum the perturbative series in C 0 (z 2 µ 2 ) to all orders, which, however, is unknown as C 0 (z 2 µ 2 ) has been calculated to only NNLO so far [21].
Note that m MS 0 is analogous to the mass renormalization in heavy-quark effective theory (HQET) [57], which is of UV origin and cannot be attributed to any shortdistance condensate. Instead, it appears as a residual mass term in the HQET Lagrangian and exists iñ h MS (z, 0, µ) at all z, i.e., whereh MS 0 (z, 0, µ) at short distance reduces to the OPE series in the brackets of Eq. (A5).
In lattice perturbation theory, one has to compute the perturbative series to very high orders of α s in order to see the renormalon effects. Nevertheless, in the MS scheme, the OPE with Wilson coefficient at a few loop orders and the condensate term, turns out to be successful in describing the static potential at short distance up to ∼ 0.25 fm [69]. One explanation is that α s in the MS scheme is larger than that in lattice perturbation theory, so the renormalon effect which is of O(α n s ) with n ∼ (2π)/(β 0 α s ) becomes significant at lower orders. This situation is similar to the OPE in QCD sum rules [70][71][72][73][74], which works well in phenomenology. The reason behind such success is probably due to a proper choice of the renormalization scale µ so that α s (µ) is small enough for the perturbative series to converge, while the µ-dependent effects in the condensate remain insignificant as they should be of the same magnitude of highest order in the truncated perturbative series [73,74].
Therefore, we approximate Eq. (A5) as where "FO" stands for fixed order, Λ(µ) is a parameter of O(Λ 2 QCD ), and we ignore the higher power corrections by working at not too large z. The µ dependence of the parameters m MS 0 and Λ is understandable because this approximation is valid for a small window of µ, and they also depend on the perturbative orders in C FO 0 if the latter does not converge fast. Note that the although the model in Eq. (A7) is not guaranteed to satisfy the RG equation forh MS (z, 0, µ), we argue that within the range of µ where it can describe the physical results, the µ-dependence in the power correction term, which is already suppressed, is weak and can be ignored.
Based on the above approximation, we fit our lattice results of the ratio in Eq. (A3) to the following ansatz, where the mass shift Moreover, since the ansatz in Eq. (A8) can describe the short-distance matrix elements well, we can correct the which is equivalent to replacingh MS 0 (z, 0, µ) by the perturbative C 0 , as in Eq. (6). Eventually, the continuum limit of the matched matrix element in Eq. (6) is which is different from MS through a perturbative matching for all z as long as z S Λ −1 QCD . Therefore, the qPDF defined as FT ofh(z, z S , P z ) is still factorizable.
Note thatm 0 (µ) introduces the ambiguity m MS 0 (µ) to the matched matrix elements. Nevertheless, we argue that C FO 0 (µ 2 z 2 ) at NNLO is different from a particular summation prescription by O(α 3 s ) contributions, which cannot be smaller than the ambiguity in m MS 0 (µ) as the latter reflects the uncertainty in summing divergent perturbative series at sufficiently high orders. Therefore, we can attribute the renormalon ambiguity inm 0 (µ) to higher loop-order effects, and estimate the latter by varying µ by a factor of √ 2 and 1/ √ 2. The range of µ we vary from cannot be too large. If µ is too small, then α s (µ) becomes too large; if µ is too large, then we need to resum the large ln(z 2 µ 2 ) in C 0 (z 2 µ 2 ). In both cases the perturbative series converges slowly. In our analysis, we scan µ within [0.9, 2.0] GeV for C NLO 0 and [1.4, 3.2] GeV for C NNLO 0 to study the scale dependence and uncertainty from renormalon ambiguity.

Fitting ofm0 and Λ(µ)
Currently, the Wilson coefficient C 0 (µ 2 z 2 ) is known to NNLO [21,29] and its anomalous dimension has been calculated at three-loop order [64], where a s = α s /(2π), L = ln(µ 2 z 2 /b 2 0 ), and b 0 = 2e −γ E ,. The factor 400 in the last square bracket is a simple guess by assuming that the constant part of the perturbative correction grows as a geometric series in the order of a s .
We also consider the RG improved (RGI) Wilson coefficient [40] where γ O is the anomalous dimension of the operator O Γ (z, µ), and β(α s (µ)) = dα s (µ)/d ln µ 2 . In this way, we can first factor out the evolution factor in Eq. (A7) as it must be satisfied by the full matrix elementh MS (z, 0, µ), and therefore construct the ratioR(z, z 0 ) in an explicitly µ-independent way. We compare C 0 and C RGI 0 at NLO, NNLO and N 3 LO at µ = 2.0 GeV in Fig. 6. The strong coupling constants at each perturbative order are defined by the corresponding Λ MS QCD with one-, two-and three-loop β functions and n f = 3, which are fixed by matching to α s (µ = 2 GeV) = 0.293. The latter is obtained from Λ MS QCD = 332 MeV with five-loop β-function and n f = 3, as has been calculated using the same lattice ensembles [75]. As one can see, at z > 0.2 fm the RGI Wilson coefficients start to deviate significantly from the fixedorder ones, which is mainly due to the large value of α s as in RGI Wilson coefficients as we evolve from µ to 1/z z z. This indicates that at z > 0.2 fm, the scale uncertainty in the perturbative series is significant due to the enhancement of non-perturbative effects, and to use OPE we should work at very short distances (z < 0.2 fm). However, there will not be enough room for varying z to satisfy z a so that discretization effects are suppressed. Therefore, in our analysis we loosen our requirement for very small z by only using the ansatz in Eq. (A8) and not considering the RGI Wilson coefficients.
In Fig. 7a, we plot an effective massm eff 0 (z) which is defined as where µ = 2.0 GeV. If the twist-four condensate is negligible, then we should expect a plateau in z, but Fig. 7a shows that it has an almost constant nonzero slope at z from 0.24 fm up to 1.0 fm. In Fig. 7b we plot its slopē which is consistent with being constant for a wide range of z. This suggests that there is considerable quadratic z-dependence from the twist-four condensate, as inclued in the ansatz in Eq. (A8).
Our results form 0 and Λ fitted fromR(z, z 0 ) for z 0 < z < z max with z 0 = 0.24 fm are shown in Fig. 8. As one can see, the two parameters remain constant in z max up to around 0.5 fm within a small window of µ, which is different with the NLO and NNLO Wilson coefficients. At larger z, the higher-twist and α s ln(z 2 µ 2 ) effects become significant, which can no longer be described by the simple ansatz in Eq. (A8). In this work, we useR(z, z 0 ) at 0.24 fm < z < 0.4 fm to fit the parameters at all µ as input for the hybrid scheme renormalization and matching. To estimate the uncertainty from the choice of µ, we will match the qPDFs obtained at different µ to the corresponding PDFs, and then evolve the final results to µ = 2.0 GeV for comparison.

Appendix B: Fourier transform (FT)
The qPDF is defined as the FT ofh(z, z S , P z ) or h(λ, λ S , P z ), Sinceh(λ, λ S , P z ) is perturbatively matched from the MS scheme, the factorization formula should still be valid for the corresponding qPDFf (x, z S , P z ) [27]. Therefore, we should integrate over all z in the FT to obtain the xdependence of the qPDF. However, due to finite lattice size effects, worsening signal-to-noise ratio and other systematics at large z, we have to truncateh(z, z S , P z ) at z = z L and extrapolate to z → ∞ to complete the FT. As a result, the small-x (x 1/λ L ) region is the most sensitive to the extrapolation model, and the corresponding systematic uncertainty cannot be well controlled. On the other hand, the reliability of the x 1/λ L region depends on the premises that theh(z) is small at z = z L and exhibits an exponential decay when z L is large enough. The first condition is easy to understand as a truncated FT will lead to an unphysical oscillation in the x-space with amplitude proportional to |h(z L )|, while the exponential decay guarantees that the FT converges fast and the qPDF at x 1/λ L has very little dependence on the specific model used in the extrapolation.
In this section, we first derive that the equal-time correlator in a hadron state does exhibit an exponential decay at large distances, then we demonstrate that including this constraint in the extrapolation will lead to a reliable FT in the moderate-to-large x region. Finally, we perform the extrapolated FT on our lattice results.

Matrix elements at large z
To begin with, let us consider a current-current correlation in the vacuum, Ω|J 5 (x)J 5 (0)|Ω , where J 5 =qγ 5 q and x 2 < 0. If we ignore the existence of zero modes and only consider gapped vacuum excitations, then where Z n is the overlap between the operator J 5 (x) and intermediate sate |n . Here m n is the mass of the intermediate state particle, and K n is the modified Bessel The correlation function should, therefore, be dominated by the exponential decay of the lowest-lying state that overlaps with J 5 (x). When the external state is a static hadron, it has also been shown that the spacelike correlations exhibit an exponential decay at large distance [56].
We are interested in equal-time quark bilinear correlators in a boosted hadron state, which can be expressed in terms of the product of two "heavy-light" currents [44,46], where the "heavy quark" hx is an auxiliary field defined along thex direction, similar to that in HQET.
Let us choose the external state to be a pion. According to Lorentz covariance, we can decompose the correlation as π(p)|q(x)γ µ hx(x)hx(0)q(0)|π(p) where the scalar functions f p,x (p·x, x 2 ) are analytic functions of p · x and x 2 . We can select the index µ such that x µ = 0. For example, we can choose µ = z when x µ = (t, 0, 0, 0) or µ = t when x µ = (0, 0, 0, z). The HQET corresponds to the timelike case, as where h v is the effective heavy-quark field moving with velocity v µ and related to the QCD heavy quark Q by the projection The lowest intermediate state |H(v) is a heavy-light meson with mass m H = m Q +Λ and momentum k µ = m H v µ , where m Q is the heavy quark pole mass, andΛ can be interpreted as the mass of the constituent light quark or binding energy. BothΛ and m Q have O(Λ QCD ) renormalon ambiguities which cancel between each other. In the Λ QCD /m Q → 0 limit,Λ should be independent of the heavy quark mass, but can depend on the light quark mass. The matrix element π(p)|qΓh v |H(v) is given by the transition form factors [76], where the form factors f 1 and f 2 only depend on v · p in HQET, and the projetion operator M(v) depends on the spin of the heavy-light meson H(v), with µ being the polarization vector for vector mesons. Therefore, Then, the correlation function becomes constitutes a large phase, so the integrand is quickly oscillating and should be suppressed. To have a naive estimate, let us assume f 1 and f 2 are constant in v · p, and the remaining integral is simply where we first obtained the result for imaginary x 0 and then analytically continued back to the real axis. Then, using Lorentz invariance and analyticity, we can obtain the result for x 2 < 0, which corresponds to the equal-time correlator that we calculate in this work. At large separation, we have which also exhibits an exponential behavior with decay constantΛ. Moreover, the correlation also includes a phase e ip·x which becomes cos(p · x) in the case of the valence quark distribution. Another important takeaway is thatΛ is a Lorentz-invariant quantity and should be independent of the external momentum.
However, it must be pointed out that the conclusion in Eq. (B13) is based on a rather crude approximation that f 1 and f 2 are constant in v · p. In practice, the transition form factors could have a pole at the mass of a heavylight meson created by the currentqγ µ h v orh v q, which is different from m H for the intermediate state |H(v) . As a result, the binding energyΛ would also be different. If we take this into account in Eq. (B11), then the result will exhibit a more complicated asymptotic behavior at large distance, where the decay constantΛ should vary among the different binding energies for the heavy-light mesons, which is similar to the observation in Ref. [56], and g is a function that can have both oscillating and non-oscillating dependence on p · x. For large enough |x|, the exponential decay should suppress the correlation and make it or its extremes decrease monotonically in magnitude.
Note that after we match the hybrid scheme matrix elements to MS, the renormalon ambiguity in the Wilson line mass, m MS 0 , is subtracted out, so the matched result should exhibit an asymptotic behavior that goes as e −(Λ−m MS 0 )|z| at large z. Therefore, the sign of (Λ−m MS 0 ) becomes crucial in determining whether it is exponentially decaying or growing.
In Since the quarks have heavier-than-physical masses in our lattice calculation, one should expect a largerΛ, so it is very likely thatΛ−m MS 0 still remains positive. After all, this can be always put to test on the P z = 0 matrix elements sinceΛ − m MS 0 is a Lorentz-invariant quantity.

Extrapolation and FT
If z L is large enough for the correlationh(z) to reach the asymptotic region, then an extrapolation that encodes the exponential decay behavior we derived in App. B 1 should lead to reliable FT for moderate-to-large x.
To be more precise, there is a rigorous upper bound for the uncertainty of FT which decreases with x.
To prove the above statement, let us consider extrapolation based on the general model where g(λ L ) =h(λ L ), and c = m eff /P z with m eff being the effective mass for the exponential decay. Motivated by QCD sum rule results, we expect m eff ∼ 0.2 − 0.5 GeV, which can be larger since we have used heavierthan-physical quark masses. Therefore, for P z ∼ 2.0 GeV in the current work, we should have c ∼ 0.10 − 0.25 or higher. Now let us compare two extrapolations h 1 and h 2 with different g 1 and g 2 . The difference between the two extrapolations, should satisfy δh(λ L ) = 0 and δh(∞) = 0. The difference in the FT with extrapolation is therefore If we can approximate δh(λ) as a flat curve within one period of the oscillatory function cos(xλ), then the integral in that region vanishes. This condition can be satisfied if |δh (λ)| x, which should be reached very quickly due to the exponential suppression at large λ. For each x, there should be a minimal integer N x which satisfies |δh (λ L + N x 2π/x)| x, so that we can approximate δf (x) as Since δh(λ L ) = 0 and δh(∞) = 0, there must be at least one extremum of δh(λ) for λ L < λ < ∞, so we have the inequality According to our estimate ofΛ − m MS 0 , c 0.1 at P z ∼ 2 GeV, so and N x ∼ O(1) should be sufficient to satisfy |δh (λ L + N x 2π/x)| x with 0 < x < 1. Therefore, in Eq. (B19) we demonstrate that there is an upper bound for the model uncertainty in the FT with exponential extrapolation, which decreases in x. The error is also proportional to |δh(λ)| max which can be much smaller than |h(λ L )| that is already close to zero. If h(λ L ) = 0.1, |δh(λ)| max = 0.05, and N x = 1, then we have which is less than 0.15 at x = 0.5 and around 15% of the central value of the qPDF as we obtain below. It is worth pointing out that our estimate of the upper bound in Eq. (B19) can be highly overestimated, as δh(λ) has an oscillation from cos(λ) and sin(λ) which are out of pace with cos(xλ) for 0 < x < 1, and |δh(λ)| max could be much smaller than |h(λ L )| and at a sharp peak within λ L < λ < λ L + N x 2π/x. Therefore, the FT with exponential extrapolation is under control for moderate and large x. Whenh(λ L ) is small enough, the model uncertainty from the extrapolation can be controlled to be much smaller than the other systematic uncertainties which are about 10% − 20% in this work.
It is worth to compare with the extrapolation error when the correlation function decreases algebraically as 1/|λ| d , which corresponds to the generic model Suppose we truncate at λ L = 10, then The power d is related to the small-x behavior of the PDF. If we parameterize the PDF as ∼ x a (1 − x) b , then with LO matching one can derive that d = min{1 + a, 1 + b} [27], which is O(1) empirically. Therefore, it will take N x 1 for the factor in Eq. (B23) to decrease sufficiently to satisfy the condition |δh (λ L + N x 2π/x)| x. As a result, the uncertainty in the FT is of orders of magnitude larger than that of extrapolation with exponential decay.
To test our claim of controlled FT error with exponential decay, we choose a particular model Suppose that the extrapolation is done at λ L = 10 with h(λ L ) = 0.15, and the parameters c and d are fitted with errors δc and δd, then we analytically FT the extrapolated result to the x-space, and calculate its error using In Fig. 9, we plot the extrapolation error against x. We have chosen different central values of the parameters c and d and fairly large uncertainties in them. The parameter d cannot have a large negative value, otherwise it would makeh(λ) grow beyond λ L . In most of the scenarios considered, the error is 0.1 for x > 0.1. As we shall see below, the actual extrapolation error is much smaller than this estimate and thus negligible when compared to the other systematic errors.  In the following, we perform the extrapolation with four different models. The extrapolation is carried out on each bootstrap sample by a minimal-square fit. For each P z , we truncateh(z) at the largest z, z >0 , where the central value ofh(z) remains positive, and choose z max = {z >0 − 2a, z >0 − a, z >0 } to estimate the truncation error. The range of z used to fit the parameters is z min ≤ z ≤ z max where z min satisfiesh(z min ) < 0.2. The continuty condition between data and model was imposed in the middle point of the fit range, namely z L , which is listed in Table II. The extrapolation models are: Exponential decay model, or "model-exp". The model for extrapolation is We have tried to fit m eff from the same range of z for P z = 0 matrix elements with a similar form, Ae −m eff |z| /|z| d , and found that m eff is around 0.1 GeV, about the same scale as the phenomenological estimate. For the P z = 0 matrix elements, we do not fix m eff , but constrain it with a prior m eff ≥ m min . To test the dependence on this prior condition, we have set m min = {0, 0.1, 0.2} GeV. Besides, we also impose A > 0 and d > 0 to ensure that the extrapolated result is positive and decreases in λ.
Power-law decay model, or "model-pow". The model is defined by setting m eff = 0 in model-exp. As the P z → ∞ limit of model-exp, model-pow can be used to give a coarse estimate of the significance of higher-twist effects, although its FT error is not well under control as we discussed above. We impose the conditions A, d > 0 so that the fitted results decrease to zero as λ → ∞. Two-parameter model with exponential decay, or "model-2p-exp". As we can see from Fig. 1, the matrix elements at λ L ∼ 6 − 10 do not show a clear exponential decay, although they can be fitted by the latter with χ 2 /d.o.f < 1 due to the large errors. This may indicate that there is oscillation inh(λ). To incorporate such dependence, we ignore the higher-twist contributions and assume that the qPDF is parameterized as By doing an inverse FT into the λ-space, the asymptotic form ofh 2p (λ) at large λ reads, Then we multiplyh 2p (λ) with an exponential decay factor as our model for extrapolation, Two-parameter model, or "model-2p". Again, we ignore the exponential decay and useh 2p as the extrapolation model, which can help us estimate the significance of higher-twist effects.
In Fig. 10 we compare the FT with different z L for extrapolation with model-exp and condition m eff > 0.1 GeV. Except for very small x, the results are consistent, and those at smaller z L have smaller errors because the error of the matrix element grows with z. Therefore, for the rest of our analysis, we simply use the largest z L for each P z .
In Fig. 11 we show the extrapolations with different models, which have noticeable differences at λ > λ L . In Fig. 12 we compare the FT with different extrapolation models as well as with the discrete FT (DFT). As we can see, the DFT introduces unphysical oscillation in the qPDF which is due to the truncation ofh(λ) at λ L . In contrast, the extrapolations are free of such oscillation, and different models yield consistent qPDFs at moderate and large x, thought they differ significantly at small x. We notice that the qPDF from model-2p extrapolation still has slight oscillations despite its agreement with the others, because the extrapolated result decays too slowly at λ > λ L . As expected, the models with exponential decay lead to regular qPDFs at x = 0, whereas modelpow and model-2p give divergent qPDFs as x → 0.
Based on the above results, we use model-exp with m eff > 0.1 GeV for the FT in our following analysis. To have a coarse estimate of the uncertainties from extrapolation model and higher-twist contributions, we look into the difference between final PDFs matched from qPDFs with model-exp and model-pow extrapolations.
Recall that although the hybrid-scheme matrix elementsh(λ, λ S , P z ) should be RG invariant, they can still as the ansatz in Eq. (A8) appear to best describe the lattice matrix elements according to Fig. 8 at these scales. The results are almost identical to each other, which shows that the renormalon-inspired model with fixed-order Wilson coefficient can indeed describe the data within a specific window of µ. At NLO, smaller µ is favored as α s (µ) is larger so that the renormalon effects become important at lower orders. In Fig. 14 we show the µ-dependence of the qPDFs from NLO-and NNLO-matchedh(λ, λ S , P z , µ, a). As one can see, the results have mild dependence on µ which becomes more significant at lower scales. Therefore, the uncertainty from scale variation will also be larger in this region.

Appendix C: Perturbative matching
In this section we perform the perturbative matching to the qPDF. Recall that Eq. (7) relates the qPDF to the PDF, The matching kernel C can be expanded to O(α s ) as The inverse matching kernel C −1 can obtained by solving order by order in α s [60], and the result is It has been shown in Ref. [60] that the inverse matching coefficient satisfies the correct RG and P z -evolution equations.

Numerical implementation of matching
Since in the asymptotic regions, and is a plus function (with "r" denotes the x = y part) that regulates the singularity at y = x, the convolution integral in Eq. (C1) is convergent and insensitive to the cutoffs for y → 0, x, ∞, as long as the qPDF is integrable. Therefore, we are able to evaluate the integral numerically within a finite range of y with a target precision. The numerical integration in Eq. (7) is time consuming, especially when we have to perform the matching for the qPDF on each bootstrap sample. Therefore, to speed up the matching procedure, we discretize the integral in Eq. (7) and reexpress it as multiplication of a matching matrix and the qPDF vector. In our implementation, our integration domain is −2.0 < y < 2.0 discretized with a step size δy = 0.001. Since the qPDF falls very close to zero at |y| = 2.0, the corresponding uncertainty is negligible as we have varied the truncation point. Note that the matching coefficient is a plus function, the step size δy also serves as a soft cutoff for the singularity near |x/y| = 1 in the plus functions. To test how well the matrix multiplication can reproduce the exact numerical intergration, we compare the NLO corrections to the qPDF from one bootstrap sample using the two methods in Fig. 15. With our current step size, the results are almost indentical for x as small as 0.01.
Moreover, to test the reliability of our inverse matching coefficient, which is obtained through expansion in α s , we compare it to direct matrix inversion. To be specific, we construct a square matching matrix C in x and y with x, y ∈ [−2, 2], which is asymmetric but has dominant diagonal elements, and then invert it to obtain the inverse matching matrix C −1 . At small α s , the matrix C can be schematically expressed as where I is an identity matrix, whereas E is O(α s ), so that its inverse can be expanded as In Fig. 16a we first test the convergence of the solution in Eq. (C8) for the NLO matching matrix. By expanding the solution to order n, we calculate the NLO matching correction to a qPDF sample, and then compare it to the result from direct matrix inversion. Since our main purpose is to compare the two inversion methods, we increase the step size to δy = 0.01 to reduce the computing time regardless the accuracy of numerical integration. We find that by increasing n, the expansion method gradually coverges to direct inversion, as expected. Of course, in perturbation theory, we should calculate the matching coefficient to n-loop accuracy for consistency, for α s is the actual power-counting parameter.
In Fig. 16b we compare the NLO and NNLO matching corrections to a qPDF sample using direct matrix inversion and the α s -expansion methods. The results are basically consistent with each other for almost the entire range of x ∈ (0, 1), except for small deviations. This is because direct matrix inversion includes all-order terms in α s , and the deviations reflect the size of higher-order effects, whose smallness shows that the perturbation series is convergent. With our current two-loop accuracy, we adopt the α s -expansion method.

Perturbative convergence
In Fig. 17 we show the matched results for the PDF from the qPDF obtained from model-exp extrapolation (with m eff > 0.1 GeV) of the NNLO-matched h(λ, λ S , P z , µ, a). As one can see, the NNLO correction is generally smaller than the NLO correction for moderate x, which indicates good perturbative convergence. Near the end-point regions, the NLO and NNLO corrections become larger than 50%, which suggests that higher-order corrections or resummation effects become important.
To see whether the NNLO matching reduces the uncertainty from scale variation, we match qPDFs at different µ to the corresponding PDFs, and then use DGLAP equation to evolve the results to µ = 2.0 GeV. We use NLO matching coefficient and LO DGLAP evolution kernel for the qPDF obtained from the NLO-matched h(λ, λ s , P z , µ, a), and NNLO matching coefficient and NLO DGLAP evolution kernel for the qPDF obtained from the NNLO-matchedh(λ, λ s , P z , µ, a). The NLO DGLAP evolution formula takes the following form, where t = ln(µ 2 /µ 2 0 ), β 0 = (11C A − 2n f )/6, P qq is the LO splitting kernel, and P V (1) qq is the NLO splitting kernel [77] for the valence quark PDF.
Since there are only a few common µ values for the NLO-and NNLO-matchedh(λ, λ s , P z , µ, a), we choose µ = 1.4 and 2.0 GeV for our comparison. In Fig. 18 we show the scale variation of the PDFs from NLO and NNLO matching, where only the central values are plotted for our purpose. As one can see, the NNLO matching correction significantly reduces the uncertainty for x 0.4 at NLO, while for x 0.4 the NNLO uncertainty band is still about a factor of one half of the NLO case. Therefore, the NNLO matching does indeed improve the perturbation theory uncertainty.
Finally, for the NNLO matching we vary µ = 2.0 GeV by a factor of √ 2 and 1/ √ 2, and then use NLO DGLAP equation to evolve the matched results to µ = 2.0 GeV, whose central vavlues are shown in Fig. 19. As one can see, there is virtually no difference between choosing  In Fig. 20 we show the P z -dependence of the PDF with NNLO matching correction. We find that despite the considerable differences between the qPDFs at P z ≤ 1.45 GeV and those at P z ≥ 1.94 GeV, the matching corrections bring the final results closer, which shows the effectiveness of LaMET. Note that the matching drives the qPDF closer to the smaller x region, so the error -FIG. 20. The PDFs from NNLO matching of the qPDFs at different P z , which is obtained from model-exp extrapolation of the NNLO-matchedh(λ, λS, P z , µ, a). bands of the PDFs also shrink after matching as they are contributed from the larger x region. Moreover, we find that the PDFs start to converge at P z ≥ 1.29 GeV, which corresponds to a boost factor of ∼ 4. As P z increases, the results becomes smaller as x → 1, which agrees with our expectation that large momentum suppresses the highertwist contributions. It is worth mentioning that both the P z -dependence and matching correction appear to be small for x as low as 0.05, which hints that the power correction and resummation effects are less severe than our naive estimate through power counting.
In Fig. 21 we compare the PDFs matched from the qPDFs with model-exp (with m eff > 0.1 GeV) and model-pow extrapolations. For a = 0.04 fm and P z = 1.94 GeV, we also added comparison to the model-2p-exp and model-2p extrapolations. Despite the differences between the qPDFs at small x, the matched results are almost identical even at the smallest x shown in the plot. Again, this is the outcome of the PDF receiving contributions from the qPDF at larger x through matching, which suggests that the extrapolation error can still be under control for x as small as ∼ 0.01. Note that the result from model-2p also shows agreement, but it includes slight oscillations in the x-space, because the extrapolatedh(λ) decays too slowly in the coordinate space. Therefore, in the region where other systematic errors are under control, the difference between model-exp and other extrapolations is negligible, and we will use the model-exp extrapolation to obtain the final results. 0.01 < x < 0.70. According to Fig. 21, the qPDF from power-law extrapolation leads to almost identical PDF after the matching correction for x as small as 0.01. Our explanation is that the matching correction drives the qPDF to smaller x, so the PDF at a given x receives contributions from the larger-x region of the qPDF which has less P z dependence. Although there are logarithms of µ/(xP z ) in the matching coefficient which become large at small x, they are always multiplied by the DGLAP splitting function, which when convoluted with the qPDF always drives the result to smaller x, thus the perturbative correction remains small even at x = 0.03.
In Fig. 23 we show uncertainty of the PDF, δf v (x)/f v (x), where δf v (x) includes both statistical and scale-variation errors. The uncertainty is ≤ 20% for 0.01 ≤ x ≤ 0.93, as x = 0.01 is the smallest x that we show in the plot, and ≤ 10% for 0.08 ≤ x ≤ 0.45. Therefore, by combining the estimates of power correction, higher-order perturbative correction, statistical and scale-variation errors, we determine the PDF at 0.03 x 0.80 with ≤ 20% uncertainty and at 0.08 x 0.45 with ≤ 10% uncertainty, which is shown in Fig. 4.