Pion Valence Quark Distribution at Physical Pion mass of $N_f=2+1+1$ Lattice QCD

We present a state-of-the-art calculation of the unpolarized pion valence-quark distribution in the framework of large-momentum effective theory (LaMET) with improved handling of systematic errors as well as two-loop perturbative matching. We use lattice ensembles generated by the MILC collaboration at lattice spacing $a\approx 0.09$~fm, lattice volume $64^3\times 96$, $N_f=2+1+1$ flavors of highly-improved staggered quarks and a physical pion mass. The LaMET matrix elements are calculated with pions boosted to momentum $P_z\approx 1.72$~GeV with high-statistics of $O(10^6)$ measurements. We study the pion PDF in both hybrid-ratio and hybrid-regularization-independent momentum subtraction (hybrid-RI/MOM) schemes and also compare the systematic errors with and without the addition of leading-renormalon resummation (LRR) and renormalization-group resummation (RGR) in both the renormalization and lightcone matching. The final lightcone PDF results are presented in the modified minimal-subtraction scheme at renormalization scale $\mu=2.0$~GeV. We show that the $x$-dependent PDFs are compatible between the hybrid-ratio and hybrid-RI/MOM renormalization with the same improvements. We also show that systematics are greatly reduced by the simultaneous inclusion of RGR and LRR and that these methods are necessary if improved precision is to be reached with higher-order terms in renormalization and matching.


I. INTRODUCTION
Parton distribution functions (PDFs) were first introduced by Feynman in 1969 [1] and describe the distribution of longitudinal momentum among a hadron's constituent quarks and gluons.They are of particular interest due to their use as inputs for the computation of scattering cross sections in high-energy hadron collisions, since the sum of convolutions of PDFs with the cross section of the corresponding parton produces a leadingorder approximation of the hadronic cross section in the collinear framework [2].In addition, they provide insight into the internal structure of the corresponding hadron.A great deal has been learned about nucleon PDFs from analysis of hard-scattering experiments since the 1960s (for reviews and latest results, see Refs.[3][4][5][6][7]), and these measurements have provided a standard against which theoretical calculations can be judged.Still greater experimental precision is anticipated at the future Electron-Ion Collider [8][9][10].The pion is the lightest hadron and is of special interest due to its interpretation as a pseudo-Nambu-Goldstone boson from the spontaneous breaking of approximate chiral symmetry; hence, an understanding of the pion PDF is of high value.
PDFs can be determined numerically (as well as from experimental data) in lattice quantum chromodynamics (QCD) with the current-current correlator method [11,12], the pseudo-PDF method [13][14][15] or large-momentum effective theory (LaMET) [16][17][18], the third of which we use in this work.LaMET examines the behavior of equaltime, spatially separated correlators on a Euclidean lat-tice and recovers the lightcone physics through perturbative matching.A great deal has been achieved applying the LaMET method to calculations of PDFs .
The x-dependent pion valence-quark PDF has been studied on the lattice in recent years [42,[46][47][48][49][50][51], and the field has matured to the point where it is important to control sources of systematic errors.Improvements have been made in the development of the hybrid renormalization scheme [52] (compared to the pure ratio-or pure regularization-independent momentumsubtraction (pure RI/MOM)-schemes), lightcone matching up to two-loop order and the inclusion of Wilsonline extrapolation in the renormalized matrix elements which reduces the presence of unphysical oscillations in the PDF.Two recent developments which we examine in this paper are the presence of large logarithms and the renormalon ambiguity both in the renormalization and matching processes.The large logarithms are handled with renormalization group resummation (RGR) [53] and the renormalon ambiguity is handled with of leading renormalon resummation (LRR) [54].The RGR method accounts for the fact that lightcone matching depends on the longitudinal momentum of the parton as well as the final renormalization scale, and the difference between them results in large logarithms, which require resummation.The RGR technique, as applied to the matching process, chooses the renormalization scale such that the logarithms vanish, and the result is evolved to the final desired scale with the DGLAP (Dokshitzer-Gribov-Lipatov-Altarelli-Parisi) equation.The LRR method is applied to the renormalization of the bare matrix elements by demanding that the short-distance behavior agrees with the theoretical predictions of the operator product expansion (OPE).These are known as Wilson coefficients and are computed as a perturbation series in the strong coupling; however, the series is not con-vergent due to the presence of an infrared renormalon (IRR).The LRR method resums the terms due to the IRR and improves the convergence of the perturbative calculation.The two methods of RGR and LRR in both the renormalization and matching processes performed up to next-to-next-to-leading-order (NNLO) result in greatly reduced systematic uncertainties in the computed PDF.The method of RGR in the matching process has been applied to the pion PDF [53] as well as the nucleon transversity PDF [55].The RGR and LRR methods were applied simultaneously to the pion PDF in Ref. [54] for the renormalization and lightcone matching processes; only the systematic errors were presented.
In this paper we apply both RGR and LRR methods on top of different renormalization schemes to examine their effects on the pion valence-quark PDFs and their systematic errors.We use clover fermions at a lattice spacing of a ≈ 0.09 fm and box length L = 64a ≈ 5.76 fm with N f = 2 + 1 + 1 flavors of highly-improved staggered quarks.The parameters are tuned so as to produce a pion of mass m π ≈ 130 MeV [56].The lattice configurations were generated by the MILC collaboration [57][58][59].The pion matrix elements is calculated with a boosted momentum of P z = 8 × 2π L ≈ 1.72 GeV with the number of measurements for each source-sink separation up to O(10 6 ).However, Ref. [56] only reports pion valence-quark PDF using the hybrid-ratio renormalization scheme with NNLO matching process; the large logarithms and the renormalon ambiguity were not included in the systemics error estimation.
This paper is organized as follows: in Sec.II we outline the methods of RGR and LRR applied both to the renormalization method and the lightcone matching.In Sec.III we display and discuss our results of the renormalized matrix elements and the x-dependent PDF, and compare with previous results in the literature.We conclude our paper in Sec.IV.

II. METHODOLOGY
In this section we outline the method of hybrid renormalization [52] with the additions of RGR and LRR described in Ref. [54] and with particular emphasis on the linear divergence and renormalon ambiguity.We then summarize the lightcone matching with both RGR and LRR improvements [53,54] which we will be used in Sec.III.

A. Improved Hybrid Renormalization Scheme
The matrix elements used in LaMET require renormalization to remove both ultraviolet (UV) and infrared (IR) divergences.The bare matrix elements used in LaMET are spatially separated correlators where ψ, ψ are the bare quark field, |π(P z )⟩ is the pion state with boost momentum P z in the z direction, γ t is a Dirac matrix and W (z, 0) = P exp −ig z 0 dz ′ A z (z ′ ) is the path-ordered Wilson line with gluon field A µ (z) connecting the two spatially separated coordinates (0, 0, 0, 0) and (0, 0, 0, z).
Early renormalization methods for LaMET matrix elements used in lattice parton calculations were the regularization-independent momentum-subtraction (RI/MOM) scheme [60] in Refs.[21-23, 25-28, 35-45, 61-63] and the ratio scheme similar to those used in the pseudo-PDF methods [13,15].The ratio scheme renormalizes the non-zero momentum matrix elements by dividing by them by the equivalent zero-momentum matrix element.The RI/MOM scheme divides the non-zero momentum matrix elements by the tree-level matrix element at a given momentum and energy scale.Later, improved renormalized schemes, such as hybrid-ratio and hybrid-RI/MOM schemes [52], were proposed to use the ratio and RI/MOM methods, respectively, up to a distance z s ≈ 0.3 fm [19,20,52,55], and for large distance, z > z s , the bare matrix elements are instead multiplied by an exponential term designed to remove both the linear divergence and the renormalon ambiguity.In such a scheme, the hybrid renormalized matrix elements h R π (z, P z ) are given by where Z(z) can be the bare matrix element at zero momentum h B π (z, P z = 0) for the hybrid-ratio or RI/MOM factor or Z RI/MOM (z, µ RI , p z R = 0) for the hybrid-RI/MOM scheme; δm and m 0 are the linear divergence and renormalon ambiguity, respectively; the normalization N = Z(0)/h B π (z = 0, P z ) sets the matrix element h R π (0, P z ) = 1, satisfying conservation of charge.The constants Z(z s ) and e −(δm+m0)zs enforce continuity at z = z s .
The linear divergence arises from the self-energy of the Wilson line [52,64] and can be determined by fitting Z(z) = Be −δm z , whether it is in the hybrid-ratio or hybrid-RI/MOM scheme.The determination of the linear divergence would appear to be a source of systematic errors; however, the renormalon ambiguity conspires to cancel this error such that the sum δm + m 0 that appears in Eq. ( 2) is fixed.To determine the renormalon ambiguity, we demand that the renormalized matrix elements agree with the OPE at short distances, z ≲ 0.2 fm.These are functions of the Wilson coefficients C 0 (z, µ), which can be computed as a perturbation series in the strong coupling.However, such series are, in general, not convergent to all orders.The asymptotic nature of the perturbation series results in an uncertainty known as the "renormalon ambiguity" [65].
Having determined the linear divergence, we determine the renormalon ambiguity by fitting the bare matrix ele-ments to at short distances as in Ref. [54].The unpolarized Wilson coefficient in Eq. ( 3) has been computed up to NNLO [66,67] and can be improved with one or both of RGR and LRR [54].The different schemes result in different central values and uncertainties for the renormalon ambiguity.
The unpolarized Wilson coefficients with Wilson length z calculated in the modified minimal-subtraction (MS) scheme at scale µ, are at NLO [67] and at NNLO [66], where l(z, µ) = ln z 2 µ 2 e 2γ E /4 , γ E is the Euler-Mascheroni constant, α s (µ) is the strong coupling at energy scale µ, C F is the quadratic Casimir for the fundamental representation of SU(3) and n f is the number of fermion flavors.The Wilson coefficients depend on the renormalization scale µ, as well as the intrinsic physical scale, and the difference between them results in logarithmic terms that need to be resummed.We perform the resummation using the renormalization-group equation (RGE) where γ(µ) is the anomalous dimension, which has been calculated up to three loops [68].We can set the energy scale such that the logarithms vanish and then evolve the Wilson coefficient to the desired energy scale by solving Eq. ( 6): where β(α s ) is the QCD β function, and γ must be computed to the same order as C 0 (z, z −1 ).The Wilson coefficients themselves are determined by a perturbation series in the strong coupling which, in general, is not convergent to all orders resulting in the renormalon ambiguity.We can account for this by improving the Wilson coefficients with LRR.Reference [54] (motivated by Refs.[69,70]) suggests modifying the Wilson coefficient according to where k = 1 for NLO, k = 2 for NNLO, and The subscript "PV" denotes the Cauchy principal value to regulate the poles in the integrand, and the various parameters are , where β n is the n th coefficient of the QCD beta function, and Γ(n) is the Euler Gamma function.The expression in Eq. ( 10) is applicable to all spin states at twist-two.We can combine the two methods of RGR and LRR to make the final high-quality Wilson coefficient Having improved the Wilson coefficients with RGR and LRR, we can determine the renormalon ambiguity m 0 using Eq. ( 3) and fully renormalize the matrix elements in the hybrid scheme using Eq. ( 2).
The lightcone PDF q v π (x, µ) is obtained by applying the perturbative matching to the quasi-PDF, qv π (x, P z ) which is related to the coordinate space matrix elements via a Fourier transform: The variable x is Bjorken x, the fraction of the hadron's longitudinal momentum carried by the valence quark.
In order to prevent unphysical oscillations in the quasi-PDF, we first extrapolate the renormalized matrix element h R (z, P z ) to infinite distance before Fourier transforming using the model proposed in Refs.[19,20]: where A, m and d are fitting parameters.This extrapolation model is inspired by the PDF having the functional form q v π (x, µ) ∼ x d−1 at small x values, which corresponds to the anticipated large-distance behavior in the renormalized matrix element in the above equation.This was also applied to the nucleon in Ref. [20].With the renormalized coordinate-space matrix element computed and extrapolated to infinite length of Wilson-line displacement, we can then compute the quasi-PDF with Eq. ( 12) and follow by performing the lightcone matching.

B. Lightcone matching
The matching process aligns the UV behavior of the quasi-PDF with that of the lightcone PDF q v π (x, µ).The two quantities are related through the matching formula where K is the matching kernel.The corrections to the matching process arise from the fact that the quasi-PDF is computed at finite momentum, whereas the lightcone PDF is defined at infinite momentum [19,71].The full matching kernel is depending on whether we renormalize in the hybrid-ratio or hybrid-RI/MOM scheme.The matching kernel has been computed up to NNLO for unpolarized quasi-PDFs renormalized in the hybrid-ratio scheme with LRR [53,66,72] as well as to NLO in the RI/MOM scheme in Refs.[63,73]. 1  The LRR modification [54] in the matching kernel can be written down as with k = 1 for NLO and k = 2 for NNLO.Since the matching process can be numerically expensive, we convert the matching kernel K of Eq. ( 15) into a matrix in x and y, K xy , then multiply a vector of quasi-PDF values by the inverse K −1 xy .The method of RGR can also be applied to the matching process; the algorithm has been derived in Ref. [53].This time, the intrinsic physical scale is that of the parton (∼ 2xP z ).We perform the lightcone matching at the scale µ = 2xP z so the logarithms vanish, and we then evolve to the desired energy scale using the DGLAP equation: where P(z) is the DGLAP kernel, which has been calculated up to three loops [74].It should be noted that the DGLAP evolution formula begins to break down at x ≈ 0.2 where α s (µ = 2xP z ) becomes nonperturbative.

A. Renormalization Matrix Elements
The first stage in the calculation is to renormalize the bare matrix elements according to Eq. ( 2).We determine the linear divergence by fitting Z(z) of Eq. ( 2) to an exponential decay Be −δm z .Our calculation from examining the exponential decay of Z(z) yields δm = 0.713 (13) and δm = 0.668 (10) GeV with Z(z) inputs from the RI/MOM renormalization factors and the zero boostmomentum bare pion matrix elements respectively.Our δm parameter are of the same order of magnitude as ANL/BNL's [19] which was δm = 0.7439(59) GeV at N f = 2 + 1 a = 0.04 fm on 310-MeV ensemble.We then proceed to determine the renormalon ambiguity from the bare boosted-momentum pion matrix elements h B π (z, P z ) and Wilson coefficients as described in Eq. ( 3).The renormalon ambiguity m 0 is obtained by fitting to a linear function in z according to We show the values of the renormalon ambiguity in Fig. 1 with the Wilson coefficient determined to different orders and in different fitting ranges.The error bars are derived from "scale variation".When we use the RGE, we set the initial scale to µ = z −1 so as to eliminate the logarithms and then evolve to the final desired energy scale of µ = 2.0 GeV.Scale variation involves setting the initial scale to c ′ × z −1 for c ′ = 0.75 and c ′ = 1.5 as was used in Ref. [54].This corresponds to a variation of approximately 15% either side of α s (µ = 2.0 GeV) in the strong coupling.The central value corresponds to c ′ = 1.0.The numerical results for the renormalon ambiguity with errors derived from scale variation are shown in Tables I and II for the hybrid-ratio and hybrid-RI/MOM renormalization scheme.The linear divergence is included in each table for the corresponding scheme.It can be seen from both the plots in Fig. 1 as well as Tables I and II that the calculations with both the LRR and RGR improvements are the most reliable due to their nice plateau behavior in the region z c = [0.12,0.2] fm as well as their greatly reduced errors compared to the other schemes.The improvement of RGR without LRR causes the renormalon ambiguity to become quite unstable and the error bars are very large.This is an evidence (also present in Ref. [54]) that the resummation of large logarithms on its own results in the Wilson coefficient being dominated by the renormalon ambiguity.The renormalon divergence begins to emerge at order k; k ∼ 1/α s (µ).At an energy scale of µ = 2.0 GeV, 1/α s ≈ 3 and so the renormalon ambiguity will not begin to emerge until N 3 LO.However, when RGR is used, the energy scale is set to an initial value of µ = 2e −γ E /z.At short distances, this is a large energy, hence a small coupling and the renormalon will not emerge at order NLO or NNLO.However, the reverse is the case at large z for which this is a small energy and the strong coupling is large.The renormalon divergence can be significant even at NLO in this region.This is the reason that the Wilson coefficients with RGR but not LRR have very large systematic errors at z ≳ 0.2 fm.The renormalon divergence must be included when RGR is used for a reliable measurement.
The lattice data renormalized in the hybrid-ratio and hybrid-RI/MOM schemes are shown in Fig. 2. The plot shows data with the renormalon ambiguity omitted and determined to (NLO+LRR)×RGR both in the hybridratio and hybrid-RI/MOM schemes.We show both statistical errors and, in the cases of (NLO+LRR)×RGR, combined statistical and systematic errors.It is clear from Fig. 2 that the central values between the hybridratio and hybrid-RI/MOM methods to an otherwise fixed order are compatible within statistical errors.The cen-tral values of the δm matrix elements can differ by up to 20% between the two renormalization schemes but have very large statistical errors.By contrast, the (NLO+LRR)×RGR results differ by no more than 7% between the two schemes and have much smaller statistical errors.The systematic errors are also very small in the (NLO+LRR)×RGR case.For this reason, the subsequent steps in the calculation of the pion PDF use the hybrid-ratio scheme since the matching kernel has been computed up to NNLO.
In Fig. 3, we show the pion matrix elements renormalized in the hybrid-ratio scheme with δm only (red) and with the renormalon ambiguity determined to NLO (blue), NLO×RGR (green) and (NLO+LRR)×RGR (purple).We can see that the systematic errors are minimized when both RGR and LRR improvements are applied simultaneously.In addition, the statistical and systematic errors are very large when RGR is applied on its own due to the RGR process enhancing the presence of the renormalon ambiguity.Between the results at NLO and those at (NLO+LRR)×RGR, the central values are compatible within statistical errors but the latter have much smaller systematic errors compared to the former.The overall z-behavior changes among the four schemes due to m 0 in the exponential term that renormalizes the matrix elements for z > z s in Eq. (2).Only NLO×RGR has a positive m 0 value; hence, the corresponding renormalized matrix elements have the largest central values.By contrast, m 0 at (NLO+LRR)×RGR is negative and NLO is even smaller.Hence, the (NLO+LRR)×RGR renormalized matrix elements are smaller than the δm terms and the NLO are smaller still.A greater suppression in z results in less pronounced oscillations in the Fourier transform (i.e. the quasi-PDF), so the NLO and (NLO+LRR)×RGR x-dependent PDFs are more steady.
We then compare the renormalized matrix elements at orders NLO, (NLO+LRR)×RGR, NNLO and (NNLO+LRR)×RGR in Fig. 4.While the lower systematic errors are insignificant with relative to the statistical errors, the upper systematic errors increase from 0.13 to 0.21 at z = 0.36 fm and from 1.5 to 3.4 at z = 0.99 fm from NLO to NNLO.However, the same errors decrease to 0.08 at z = 0.36 fm and 0.8 at z = 0.99 fm from NNLO to (NNLO+LRR)×RGR.In addition, the relative systematic errors decrease from (NLO+LRR)×RGR to (NNLO+LRR)×RGR.Not only does this reaffirm the benefits of including RGR and LRR in our calculations, but it also shows that if one wishes to renormalize to order NNLO or higher, it is necessary to account for both large logarithms and the renormalon ambiguity.In other words, handling of the systematic errors must improve in parallel with higher order terms in the renormalization process.If systematics are not accounted for in the calculations, the results will deteriorate at higher orders.This demonstrates the need to account for the large logarithms and the renormalon divergence in the renormalization of LaMET matrix elements.The errors in m0 are derived from scale variation described in Sec.III A. Each m0 is added to the linear divergence δm = 0.713 (13) GeV in the hybrid-RI/MOM scheme for z > zs.
Note the greatly reduced errors in m0 when RGR and LRR are applied simultaneously.

B. Parton Distribution Functions
To obtain the pion quasi-PDF, we first need to extrapolate the renormalized matrix elements at large Wilson-

FIG. 2:
δm only renormalized matrix elements in the hybrid-ratio (red) and hybrid-RI/MOM (blue) schemes; (NLO+LRR)×RGR renormalized matrix elements in the hybrid-ratio (green) and the hybrid-RI/MOM (purple) schemes.In the two cases of (NLO+LRR)×RGR, the solid error bars are statistical and the dashed error bars are combined statistical and systematic, the latter arising from the scale variation.Except for the δm-only calculation in the hybrid-ratio scheme, the plotted data have been offset from their exact z value to allow for readability.
line displacement to infinity.The Wilson-line displacement must be large enough (≳ 0.5 fm) that we can realistically use the small-x model corresponding to the large distance behavior; in this work, we choose the fitting range [8a, 14a] ≈ [0.72, 1.26] fm.Several extrapolation models were tested by the ANL/BNL collaboration in Ref. [19] for the pion matrix elements and the final one used was that of Eq. ( 13) with the constraints A > 0, d > 0 and m > 0.1 GeV; the last two constraints ensure that the large distance behavior decays sufficiently quickly to obtain a convergent Fourier transform.Once we obtain quasi-PDF, the final step is to apply the match- FIG.4: NLO (blue), (NLO+LRR)×RGR (purple), NNLO (magenta) and (NNLO+LRR)×RGR (orange) matrix elements renormalized in the hybrid ratio scheme.The solid error bars are statistical and the dashed error bars are combined statistical and systematic errors, the latter arising from the scale variation.All plotted data except for NLO have been offset from their exact z values to allow for readability.
Before we study the full statistical and systematic errors of pion valence-quark PDF, we first examine their x-dependent NLO, NNLO, NLO×RGR, (NLO+LRR)×RGR and (NNLO+LRR)×RGR systematic errors in Fig. 5.We see that the lower systematic errors are larger for the two orders NNLO and NLO×RGR than they are at NLO, which was also found by Ref. [54].In addition, the systematic errors decrease significantly from NLO, to (NLO+LRR)×RGR, and further reduction from (NLO+LRR)×RGR to (NNLO+LRR)×RGR.We also observe the central values for (NLO+LRR)×RGR and (NNLO+LRR)×RGR are much closer to each other than those for NLO to NNLO showing better convergence going to higher order with LRR and RGR improvements.These qualitative behaviors are also consistent with what was found in the earlier pion-PDF study [54].
We compare the pion PDFs from hybrid-RI/MOM and hybrid-ratio in Fig. 6; the renormalized matrix elements can be found in Fig. 2. The RGR matching process begins to break down at small-x, where the strong coupling becomes nonperturbative at energy scale µ = 2xP z .For this reason, we do not plot the PDF for x ≲ 0.22.We noted that the two (NLO+LRR)×RGR renormalization schemes PDFs produce very similar central values, while the PDFs from the two δm-only renormalization schemes differ noticeably.This is understandable, since the PDFs from the δm-only schemes have different linear-divergence contributions, as discussed in Sec.III A. Furthermore, the two corresponding (NLO+LRR)×RGR renormalized matrix elements have very similar central values, as shown in Fig. 2. In addition, the lightcone matching is the same for the two schemes when the RI/MOM matrix elements are evaluated at p z R = 0, as was demonstrated in Ref. [73] at NLO, so the matching process does not cause a further differences in the xdependent results for hybrid-ratio and hybrid-RI/MOM schemes.Overall, we found the pion valence-quark PDF to be consistent between the hybrid-RI/MOM and hybrid-ratio renormalization schemes.For the rest of this work, we will focus on the pion valence-quark PDF results with hybrid-ratio renormalization scheme.
We compare x-dependent PDFs in the hybridratio scheme with δm-only, NLO, NLO×RGR and (NLO+LRR)×RGR improvement in Fig. 7. Going from NLO to NLO×RGR, the relative systematic errors increase by as much as 40%, showing that the RGR process on its own can enhance the presence of the infrared renormalon.However, the relative systematic errors are reduced by as much as a factor of two when going from NLO×RGR to (NLO+LRR)×RGR, demonstrating that FIG.6: The lightcone PDFs renormalized in hybrid-ratio and hybrid-RI/MOM schemes with δm only (dotted cyan and dashed orange) and at (NLO+LRR)×RGR (hatched blue and hatched red).In each PDF, the inner darker bands show statistical error, while the outer lighter bands show combined statistical and systematic errors (due to scale variation).The dark-gray region is where the LaMET calculation begins to break down, and the light-gray region is where the RGR matching begins to break down.
it is necessary to include the LRR method when RGR is used.We also see that the central values of the two RGR plots (NLO×RGR and (NLO+LRR)×RGR) are lower than those of NLO in the region x ≲ 0.6.This is to be expected, since the logarithmic terms that are resummed in the RGR matching process are ln µ 2 /4x 2 P 2 z , which become larger as x decreases.Thus, their effect is more significant at small x, and the effect of resumming them is to lower the central value.
The pion valence-quark PDFs for one-and two-loop treatments are compared in Fig. 8: NLO (dotted cyan), (NLO+LRR)×RGR (solid blue), NNLO (dashed orange) and (NNLO+LRR)×RGR (solid red) improvements.We would normally expect a calculation to yield more precise results as we go to higher order (e.g.NLO to NNLO).However, this is not the case in Fig. 8; the systematic errors increase from NLO to NNLO.It is necessary to account for the effects of large logarithms and the renormalon divergence with the methods of RGR and LRR, respectively.We can see that the systematic errors decrease from (NLO+LRR)×RGR to (NNLO+LRR)×RGR, showing that, in this case, the calculation becomes more precise.Figure 8 is a demonstration that it is necessary to control the sources of systematic errors if higher-order terms are to be included in the lightcone matching.Since the systematic errors decrease from (NLO+LRR)×RGR to (NNLO+LRR)×RGR by 10% to 15%, we see the benefits of including higher-order terms in the matching and renormalization.This demonstrates that it is necessary to include both RGR and LRR if we wish to renormalize and match to two-loop level.
We compare our numerical results with those of global fits performed by the JAM [75] and the xFitter [76] col-FIG.7: The lightcone pion valence-quark PDFs as a function of x renormalized in hybrid-ratio schemes with δm only (dotted cyan), NLO (dashed orange), NLO×RGR (hatched blue) and (NLO+LRR)×RGR (hatched red) improvement.In each PDF, the inner darker bands show statistical error, while the outer lighter bands show combined statistical and systematic errors (due to scale variation).The dark-gray region is where the LaMET calculation begins to break down, and the lightgray region is where the RGR matching begins to break down.laborations in Fig. 9.We scale the results of the xFitter collaboration to match our convention of valence quark number conservation: This is also the same convention that the JAM collaboration uses.We focus on the JAM results with nextto-leading logs (NLL) since they offer a better systematic improvement of the valence quark PDF.Our results include both statistical errors and combined statistical and systematic errors.We show our results at NNLO and (NNLO+LRR)×RGR to show the difference between two-loop matching with and without RGR and LRR improvements.Our (NNLO+LRR)×RGR results agree within two sigma with the JAM results in the interval x ≈ [0.[76] ("xFitter'20") at NLO (dotted orange) and [75] ("JAM'21") at NLO with next-to-leading logarithms (dashed blue).The xFitter'20 results have been scaled so as to satisfy Eq. (19).

IV. CONCLUSION
In this paper we have computed the pion valencequark PDF with physical quark mass at boost momentum 1.72 GeV with the improvements of RGR and LRR.We use a physical pion mass ensemble generated by the MILC collaboration [57][58][59] with lattice spacing a ≈ 0.09 fm, and 2+1+1 flavors of highly improved staggered quarks in the sea and clover fermion for the valence sector.We compute and compare the LaMET matrix elements renormalized in the hybrid-ratio and hybrid-RI/MOM schemes as well as the corresponding xdependent PDFs with matching performed at both oneand two-loops levels.
We report the impacts of different levels of improvement in the renormalization and matching from δm-only to implementation of RGR and LRR.We found that the renormalized matrix elements in both the hybrid-RI/MOM and hybrid-ratio scheme are consistent with each other within the statistical errors, but the former have slightly higher central value across all the Wilsonline displacements we studied.The systematic errors from scale variation in the renormalized matrix elements in the hybrid-RI/MOM and hybrid-ratio schemes are greatly reduced by the simultaneous application of RGR and LRR at one-and two-loop level respectively.However, the application of RGR on its own at either level increases the systematic errors, due to its enhancement of the renormalon ambiguity.We found our pion valencequark PDF in hybrid-RI/MOM and hybrid-ratio scheme to be consistent with each other within one sigma, both in the case of δm only and with RGR and LRR improvement.Unfortunately, there are no results in the literature with which we can compare directly, since the only previously calculated NNLO pion valence-quark PDF was given by Ref. [56], renormalized in the hybridratio scheme, but these LaMET systematics were not included.In studying the x-dependent PDFs, we also show that there is an increase in systematic errors between NLO and NNLO but a decrease in systematic errors from (NLO+LRR)×RGR to (NNLO+LRR)×RGR and better convergence of the central values.This demonstrates that the inclusion of higher orders in both the renormalization and matching processes must be supplemented with an improved handling of the systematic errors if the results are to be made more precise.Our (NNLO+LRR)×RGR result is the most reliable, since it includes the highest-order lightcone matching and accounts for both the large logarithms and the renormalon ambiguity.We also compare our results to global fits performed by the JAM and xFitter collaborations.Overall, we have reasonable agreement up to differences due to improvements from RGR and LRR.
) at the three z values {z c − 0.02 fm, z c , z c + 0.02 fm} for different central values z c within the range of validity of the OPE (z c ≲ 0.2 fm).The terms m 0 and c are fitting parameters.The Wilson coefficient C 0 (z, µ) in the formula can be determined to NLO or NNLO and may have the improvements of LRR and/or RGR as detailed in Sec.II A. The coefficient of z in the linear fit is the renormalon ambiguity m 0 and the term c is the y-intercept, which is irrelevant for this calculation.

FIG. 1 :
FIG. 1: Renormalon ambiguity determined in the hybrid-ratio scheme in multiple z ranges.The value m0(z) is determined by fitting Eq. (18) to the z values {z −0.02 fm, z, z +0.02 fm}.The vertical width of each band corresponds to the systematic error derived from scale variation described in Sec.III A. We show m0(z) determined at NLO, NLO×RGR, (NLO+LRR)×RGR (left figure), NNLO, NNLO×RGR and (NNLO+LRR)×RGR (right figure) in blue, green, red, orange, cyan and purple, respectively.Note that the errors are minimized in the two cases in which RGR and LRR are applied simultaneously.

FIG. 3 :
FIG.3: δm-only (red), NLO (blue), NLO×RGR (green) and (NLO+LRR)×RGR (purple) matrix elements renormalized in the hybrid-ratio scheme.In the cases of NLO, NLO×RGR and (NLO+LRR)×RGR, the solid error bars are statistical and the dashed error bars are combined statistical and systematic errors, the latter arising from the scale variation.The same three sets of plotted data have been offset from their exact z values to allow for readability.

FIG. 8 :
FIG.8:The lightcone pion valence-quark PDFs as a function of x renormalized in hybrid-ratio scheme with NLO (dotted cyan), (NLO+LRR)×RGR (hatched blue), NNLO (dashed orange) and (NNLO+LRR)×RGR (hatched red) improvement.In each PDF, the inner darker band shows statistical error, while the outer lighter band shows combined statistical and systematic errors (due to scale variation).The dark-gray region is where the LaMET calculation begins to break down, and the light-gray region is where the RGR matching begins to break down.

TABLE I :
(10)rmalon ambiguity m0 for hybrid-ratio scheme determined in the fitting range z ∈ [0.14, 0.18] fm.The errors in m0 are derived from scale variation described in Sec.III A. Each m0 is added to the linear divergence δm = 0.668(10)GeV in the hybrid-ratio scheme for z > zs.Note the greatly reduced errors in m0 when RGR and LRR are applied simultaneously.
[77]0.83] and with the xFitter results in the interval x ≈ [0.35, 0.95].The difference we see in the mid-x region is likely due to the different levels of improvement.Note that the JAM results[75]computed with and without NLL also shows 1-2 sigmas difference within the same analysis frame; the latter one has better agreement with the xFitter results.Likely we are seeing difference due to the different systematic improvements from NLL, RGR and LRR.In the large-x region, our results are larger than those of JAM and xFitter who use the parametrization form of x a (1 − x) b to ensure the PDF goes to 0 at x = 1, while ours does not.This same behavior has occurred in other LaMET framework calculations of the valence-quark pion PDF done by ANL/BNL group[77]at NNLO level.Overall, there is reasonable agreement between our numerical results, previous LaMET calculations and the global fits of JAM and xFitter.