Leading Power Accuracy in Lattice Calculations of Parton Distributions

In lattice-QCD calculations of parton distribution functions (PDFs) via large-momentum effective theory, the leading power (twist-three) correction appears as ${\cal O}(\Lambda_{\rm QCD}/P^z)$ due to the linear-divergent self-energy of Wilson line in quasi-PDF operators. For lattice data with hadron momentum $P^z$ of a few GeV, this correction is dominant in matching, as large as 30\% or more. We show how to eliminate this uncertainty through choosing the mass renormalization parameter consistently with the resummation scheme of the infrared-renormalon series in perturbative matching coefficients. An example on the lattice pion PDF data at $P^z = 1.9$ GeV shows an improvement of matching accuracy by a factor of more than $3\sim 5$ in the expansion region $x= 0.2\sim 0.5$.


Introduction
The parton distribution functions (PDFs) of the hadrons, describing the momentum distribution of quarks and gluons inside, serve as key inputs in theoretical predictions of the hadronic cross sections in high-energy scattering experiments [1].They have been determined to high precision by global fittings to experimental data [2][3][4][5][6][7].At present, a direct way to calculate the Bjorken-x-dependent PDFs non-perturbatively is through lattice QCD simulations of Euclidean correlation functions, interpreted, e.g., in the framework of the large momentum expansion or large-momentum effective theory (LaMET) [8,9].A recent summary of PDF lattice calculations with various other approaches can be found in Ref. [10].
Since the proposal of LaMET, much progress has been made within the framework in first-principle predictions of PDFs [10][11][12].The calculations in large-momentum expansion have been improved to a level that the precision control at various steps has become important [13].Apart from the usual requirements of lattice QCD, such as physical pion mass, large volume and small lattice spacing, one needs to properly renormalize the lattice data to achieve a reliable continuum limit and matching to physical PDFs with high precision.Since the LaMET expansion parameter is both in coupling α s (∼ P z ) and in powers of Λ QCD /P z , and since most of the calculations so far are limited to P z ∼ few GeV, it is critically important to understand the effects of power corrections, which have not been properly addressed in the literature so far.In a recent paper [14], it has been shown that the proper perturbative scale in LaMET matching is around 2xP z , which implies that the power correction of type Λ QCD /2xP z can reach ∼ 30% level at lower x limit of 0.2, and as large as 15-20% at x = 0.5.Therefore, achieving twist-three power accuracy is necessary for a precision prediction of PDFs in the lower x region.
Once done, the naive counting puts the twist-four power correction to order (Λ QCD /2xP z ) 2 which is below 10% in the LaMET expansion region.
In the matching formula for PDFs, there are two sources of the twist-three power corrections.First, the linearly-divergent self-energy of Wilson line has an intrinsic ambiguity O(Λ QCD ) due to long-range nonperturbative effects [15][16][17][18][19][20][21][22][23] and, corresponding to the freedom of a finite term for subtraction in renormalization.Different choices introduce a scheme dependence through a twist-three mass parameter m 0 in the form of a e m0z factor in quasi-PDFs, corresponding to a O(Λ QCD /xP z ) freedom in the momentum space quasi-PDF [24].Meanwhile, the perturbative matching kernel to the light-cone PDFs is an asymptotic QCD perturbation series in α s with the coefficients grow factorially at high orders [25], a phenomenon called infrared renormalons [26].The leading renormalon series can be summed/regularized in multiple ways that differ by a term linear in the spatial Wilson line length z, also introducing a twist-three scheme dependence.Leading power accuracy is achieved by choosing both renormalization and matching schemes consistently, thereby removing O(Λ QCD /xP z ) twist-three contribution entirely in the light-cone PDF extraction.
In this paper, we show how to achieve twist-three power accuracy in calculations of PDFs for the first time.We use perturbative matching coefficients accurate to next-to-next-to leading order (NNLO) plus the leading renormalon contribution estimated from the heavy quark pole mass calculations [27,28].The results of the latter at n f = 0 agree well with a precise lattice heavy-mass perturbative series up to 20th order in the QCD coupling [19].We adopt the renormalization condition that the continuum-limit lattice matrix elements have a shortdistance (z) expansion consistent with Wilson coefficients resummed through the principal value (PV) prescription.Using this, we extract the finite renormalization or twistthree mass parameter m 0 from the P z = 0 lattice matrix element in the perturbative z region [29].The scaleindependence of m 0 and fast convergence between NLO and NNLO perturbative calculations are important indi-cators that leading power contribution is under control.As a demonstrative example, we show that the uncertainty in x-dependent pion PDF is significantly reduced after the leading power correction.
Linear divergence and uncertainty in twist-three contribution Most LaMET calculations start from the coordinate space bare correlator [8] where W (z, 0) = P exp[ z 0 A z (z )dz ] with path ordering is a gauge link along the z-direction, making the operator gauge invariant.Due to the self-energy of the Wilson line, the above operator has a linear divergence which can be renormalized in the multiplicative factor of e −δm(a)z , equivalent to the quantum correction to the mass δm(a) of an auxiliary "infinitely-heavy quark" (whose propagator Q(z) Q(0) generates the Wilson line W (z, 0) [24,[30][31][32][33]).However, the quark carries a color charge, which leads to a ill-defined infrared contribution to the mass.Thus, one can choose the renormalization factor either with a "short-distance mass" defined by an infraredregulated "pole mass" [19], or with a non-perturbatively defined mass parameter in terms of a physical matrix element [13,29,34].Different mass subtractions (labelled by a parameter τ ) differ by the intrinsic non-perturbative scale, ∆ (δm(a, τ )) ∼ Λ QCD .Thus, besides the usual renormalization scale µ for logarithmic divergences, the quasi-PDF matrix element has the τ dependence from the renormalization of the linear divergence.
To achieve high precision, it has been suggested to eliminate the power-divergence from lattice data through a self-renormalization approach [29].The lattice ultraviolet (UV) regulator a dependence of the P z = 0 quasi-PDF matrix elements h B (z, a, P z = 0) has been parameterized by fitting the lattice data on multiple lattice spacings.Multiple sets of parameters were found corresponding to the same minimal χ 2 fit, can be regarded as different choices of τ .These different sets show an O(Λ QCD ) ambiguity in the mass renormalization δm(a, τ ), an example showing the τ -scheme dependence.
The τ -dependence entails that the short-distance (z) expansion [35] of the renormalized correlator must contain a twist-three non-perturbative parameter m 0 (τ ) of order Λ QCD contributing as a linear-z term (the sign is chosen so m 0 contributes positively as "residual mass"), where λ = zP z , a k are spin-k twist-2 matrix elements, and perturbative series C k are the associated Wilson coefficients.The twist-three contribution is universal in the sense that it multiplies the leading-twist term in the same manner independent of the spin of the local operators.On the other hand, m 0 (τ ) does depend on the external states in which the correlator matrix elements are taken as well as kinematic Lorentz invariants formed by external momenta.In the LaMET expansion, the above twist-three term translates to a linear one in inverse large momentum 1/P z .
To extract the twist-two PDFs up to leading-power accuracy, one needs to fix m 0 (τ ) from a non-perturbative matrix element.A natural choice is the P z = 0 case where only k = 0 term in Eq. ( 2) survives when λ = 0.With renormalized lattice data h R (z, P z = 0, µ, τ ), we can fit m 0 (τ ) with the expansion at linear-z accuracy using known a 1 = 1 and MS series C 0 (α s (µ), µ 2 z 2 ).
The fitted m 0 (τ ) can alternatively be absorbed into an additional finite renormalization to define a new τ -independent renormalizaiton factor Z R (z, a, µ) ∼ exp[δm(a, τ )z+m 0 (τ )z] and the renormalized quasi-PDF q(x, P z , µ) for P z > 0, the latter can be used in the standard LaMET matching up to leading power accuracy but without the explicit twist-three contribution.
The determination of m 0 (τ ) requires a rigorous relation between the bare lattice data h B (z, P z , a −1 ) and the renormalized matrix elements h R (z, P z , µ, τ ) in MS scheme, or equivalently, the rigorous definition of the We start from h R (z, P z = 0, µ = z −1 ) at small z in MS scheme without regularizing the twist-3 ambiguity, and solve the renormalization group (RG) equation to separate the dependence on z and the UV regulator a −1 (or UV renormalization scale µ in dimensional regularization), where the anomalous dimension γ(α(µ)) of h R (z, 0, µ) has been calculated up to 3-loop order for the quasi-PDF operator [36].The above equation can be solved when z is perturbative by evolving from the initial point µ = z −1 .The solution at scale µ is then where I(µ) is a analytic function of α(µ) because both γ(α) and β(α) are polynomials of α when truncated at a certain order in perturbation theory.Similarly, for the lattice scheme with explicit linear divergence, the matrix elements now have the following form: where h lat is the perturbation series evaluated at scale z −1 in the lattice scheme.In both schemes, the zdependence is now completely factorized.In the lattice matrix element, we have manifestly separated out the mass dependence with the linearly-divergent mass correction δm(a).The z-dependence is physical, thus should not depend on scheme choices.So we can identify h lat (z, 0, z −1 )e −I lat (z −1 ) in the lattice scheme and h(z, 0, z −1 ) exp −I(z −1 ) in the MS scheme.The above expressions allow us to express the lattice matrix element in terms of the MS perturbation series h B (z, 0, a −1 ) = h R (z, 0, z −1 )e −I(z −1 ) e I lat (a −1 ) e −δm(a)z , when combined with Eq. ( 5), indicating that a twist-3 accurate renormalization constant can be defined as which also works for P z > 0 lattice data.Now the mass correction δm(a, τ ) is fixed in the τ -scheme, and a corresponding non-perturbative mass parameter m 0 (τ ) in Eq. ( 2) is absorbed into the renormalization to guarantee the twist-3 accuracy.After a rearrangement, Eq. ( 7) becomes − ln h B (z, 0, a −1 )e −I lat (a −1 ) e δm(a,τ )z , that allows us to extract m 0 (τ ) as a slope of such a quantity.
In Ref. [29], the fixed (next-to-leading, NLO) order Wilson coefficient [35,37] C 0 (α s (µ), µ 2 z 2 ) = 1 + α s (µ) 3π 3 ln z 2 µ 2 e 2γ E 4 + 5 (10) has been used to fit m 0 (τ ).The use of this fixed-order expression in Eq. ( 9) will in principle introduce small deviations from unresummed α n s (µ) ln n (z 2 µ 2 ) terms thus slightly alters the z-dependence.The result will then depend on the choice of both µ and z.Here we first follow their strategy and examine the uncertainty in this fitting method.For a single lattice spacing, I 0 = e −I lat (a −1 ) is a constant thus is absorbed into the interception and not affecting the slope m 0 (τ ), when fitted to an approximate form of Eq. ( 9) at fixed order perturbation theory m 0 (τ )z − I 0 ≈ ln e −δm(a,τ )z h R (z, 0, µ, τ )/h B (z, 0, a −1 ) .(11) We fit m 0 (τ ) to Eq. ( 11) in a range [z − a, z + a] from the pion quasi-PDF lattice data from the ANL and BNL collaboration [13] with lattice spacing a = 0.04 fm.As shown by cyan and orange bands in the top panel in Fig. 1, the fixed-order (NLO or NNLO) C 0 introduces a large uncertainty from scale µ variation from 1 GeV to 4 GeV, corresponding to all possible relevant physics scales (2xP z ) in the problem.We can resum the large logarithmic terms α n s (µ) ln n (z 2 µ 2 ) in C 0 to reduce the µ dependence using the renormalization group resummed (RGR) renormalization constant defined in Eq. ( 8).We then fit m 0 to Eq. ( 9), shown as the hatched green band on the top panel, has a strong dependence on z, and becomes unreasonably large at larger z > 0.2 fm.This is a indication that there is a significant contamination effect from unaccounted higher-order terms in twist-two C 0 , which has logarithmic dependence in z.The z-dependence is altered by the truncation in α s (z −1 ).Note that the uncertainty shown is purely systematic error, because the data are so clean at P z = 0 and z < 0.2 fm that can be completely neglected in the extraction of m 0 .This large uncertainty in the twist-three mass parameter m 0 can be translated to that for light-cone PDFs in LaMET calculations.The extracted isovector light-cone PDF is shown in the bottom plot of Fig. 1.With the fixed order C k , the uncertainty in m 0 yields > 30% error near x ∼ 0.3.With renormalization group improvement, a significant uncertainty in the extracted PDF still exists, shown as the hatched green band with NLO+RGR label in the same plot.These large uncertainties indicate that improving calculations up to twist-three accuracy is crucially important for lattice data at P z ∼ 2 GeV.

Power accuracy in matching coefficients through infrared renormalons
It is well-known that perturbative Wilson coefficients C k (α s ) = n c kn α n s are asymptotic series because of infra-red contributions at large orders, a phenomenon called infrared renormalons (IRR) [26].For the quasi-PDF operator, the leading IRR comes from the long-distance contributions to the self-energy of the charged Wilson line, and resulting c kn growing factorially ∝ n!(β 0 /2π) n at large-n orders [25], where β 0 = 11 − 2n f /3 (n f is the number of active quark flavors) is the first term of the QCD beta function.This behavior is exactly the same as the perturbation series for the "pole" mass of a quark [19,27].
The perturbative runaway infrared contributions at large orders must be regularized, or equivalently, one has to specify a way to resum the perturbative series.Different methods of resummation/regularization yield results differing by the order of the minimal term in the series, at n ∼ 2π/β 0 α s [16].A simple estimation shows that this is twist-three contribution O(zΛ QCD ), a linear power in z (other renormalon poles corresponding to higher-powerz/twist terms).Therefore, the twist-two contributions themselves are ambiguous up to higher-twist contributions, and the twist-three parameter m 0 (τ ) depends on the resummation/regularization method for the leading renormalon series in C k (α s ).There have been a number of proposals for renormalon regularization in the literature [17].For example, one can define subtraction for infrared contributions at every order in perturbation series.One can also calculate the series in usual MS or lattice method and define an all-order sum through a Borel transformation.In this regard, it has been advocated to use the principal value prescription to regulate the renormalon poles in the Borel integral [22,23].We consider a prescription for C k (α s ) together with the UV regularization of the correlator as the complete τ -scheme.
Calculating C k to large orders is extremely challenging if not impossible.The only type of diagrams we can calculate to all orders analytically is the bubble-chain ones, which are leading in the large-n f limit in Abelian gauge theory.When restored to the full β function in non-Abelian theories, the leading asymptotic coefficient is n!(β 0 /2π) n , which is indeed a renormalon corresponding to the quark pole mass.On the other hand, the actual overall strength (represented by an overall factor N m ) of the renormalon can be quite different from the bubble-chain diagram result, N bubble chain m = 0.98 [38].Fortunately, lattice numerical calculations to perturbative large orders have become possible for very limited quantities [19,20].In particular, the pole mass on pure gauge ensembles (n f = 0) has been calculated to n ∼ 20, m = µ n r n α n+1 s , and confirmed the conjecture on the leading IR renormalon form at large n [19], where b = β 1 /2β 2 0 and c 1 = (β 2 1 − β 0 β 2 )/(4bβ 4 0 ) are from higher orders in the QCD beta function.Moreover, the mass renormalon strength has been determined to be N m (n f = 0) = 0.660(56) in MS scheme.Furthermore, using an analytical method in Ref. [27], one can estimate N m (n f = 0) = 0.622 consistent with the previous determination, and with dynamic fermions, we used N m (n f = 3) = 0.575 in this work.
Using the above state-of-art knowledge on the mass renormalon, we have the leading renormalon contribution for C k after a Borel transformation, where we have chosen a PV prescription for regulating poles on the Borel u-plane.We define a leading renormalon resummation (LRR) result by resumming the leading divergent contributions to all orders at µ = z −1 , while keeping the lower-order expansion of C LRR (α s ) the same as C 0 (α s (µ), z 2 µ 2 ), Now C LRR (α s ) contains not only the fixed-order results calculated explicitly, but also higher-order (twist-two) perturbative terms contributing to the leading factorial growth.
We compare the k = 0 Wilson coefficient C 0 (α s (µ), z 2 µ 2 ) for the fixed-order (NLO+NNLO), fixed order (NLO+NNLO) with RGR, and the LRR-improved formalism on the top panel in Fig. 2. The error bands in RGR are obtained by varying the resummation scale from 0.75z −1 to 1.5z −1 , corresponding to about 30% change in the strong coupling.Going beyond the lower bound in our data range, the perturbation theory breaks down.While there is a large difference from NLO to NNLO in fixed-order calculations with or without renormalization group improvement, the LRR results show much better convergence in the perturbative region z < 0.3 fm, and much smaller dependence on the scale variation, indicating that NNLO term is already dominated by the leading renormalon.
We show the NLO LRR-improved m 0 result as blue band together with fixed-order results on the lower panel in Fig. 2. By including the leading renormalon, there is now a clear window near z = 0.12 fm for a constant m 0 (τ ) = 0.161 +0.025 −0.002 GeV for NLO with much smaller uncertainty.Thus Eq. ( 2) achieves the linear-z accuracy when the leading renormalon series is resummed.We also show the NNLO renormalon-resummed results m 0 (τ ) = 0.164 +0.016 −0.003 GeV as the red band to demonstrate the good convergence with this method, consistent with the blue band and has smaller scale dependence at small z.The difference between the non-perturbative lattice result and the perturbation series is well described by the linear dependence in z in the perturbatively-reliable region.This gives us confidence that we have reached twist-three power accuracy for describing the P z = 0 matrix element.
PDF matching to leading power accuracy We commented after Eq. ( 2) that the leading-power correction term m 0 (τ ) multiplies the twist-two matrix elements in the same way independent of their spin.This observation is still valid when m 0 (τ ) plays the additional role to account for the scheme dependence in regularizing the leading renormalon divergence in the coefficient function C k (α s ).This is because all C k (α s ) has the exactly the same leading renormalon series as a quark "pole" mass.Moreover, this leading renormalon series exponentiates such that it matches exactly the mass renormalization of the Wilson line in the quasi-PDF operator.Therefore, if we renormalize the large-P z spatial correlators h B (z, P z , a) with the Z R (z, a, µ, τ ) and m 0 (τ ) from the previous section, the resulting h R (z, P z , µ, τ ) can be matched to light-cone PDFs with τ prescription, e.g., PV, for the leading renormalon in the matching coefficient without any explicit leading power corrections.
To extract the x-dependent PDF with linear accuracy in 1/P z , one needs either to match the Euclidean coorelation functions to the lightcone with linear-z accuracy in coordinate space then perform a Fourier transformation, or to first obtain the x-dependent quasi-PDF then match it to PDF with the linear-1/P z accuracy.The former approach faces the problems of the breaking down of twist expansion and the existence of the Landau pole when going beyond z ∼ Λ −1 QCD .So we take the second approach, which then requires the LRR correction to the momentum-space matching.In this approach, to avoid the Landau pole, we derive the regularized LRR correction with fixed renormalization scale µ, then apply the renormalization group resummation in momentum space, as shown in Ref. [14].Utilizing the fact that the LRR correction is universal to all C k , the correction to the momentum-space matching kernel C(x, y, µ, P z ) is just the Fourier transformation of the coordinate-space corrections: ∆C LRR (x, y, µ, P z , τ ) = yP z dz 2π e i(x−y)zPz Note that the integrand is linear in z, thus is integrated to a singular function which includes derivative of delta functions δ (x − y).To avoid the instability of its numerical implementation, it is helpful to introduce a small exponential decay to regularize the long-tail, where the smaller parameter m → 0 along with a finer discretization of the x ∈ [0, 1] region eventually recovers its continuum version.With finite but small m , the extra correction introduced by this regularization is of O(z 2 m Λ QCD ), thus is converted to O( m Λ QCD /(2xP z ) 2 ), not breaking the linear accuracy.In the hybrid scheme [24], where the integral is from z s to ∞, we obtain where we have used for simplicity, and a plus function to compensate a neglected δ(1−ξ) term in this correction.
Testing with some different m ∈ [20, 100] MeV, we find the results are consistent and stable.
We then apply the leading-renormalon resummed matching coefficients and the corresponding m 0 to the analysis of the pion PDF lattice data [13], with results shown in Fig. 3.The results from fixed-order perturbation theory from Fig. 1 are shown again for comparison.The m 0 (τ ) used in calculating the blue (red) band is from the bottom plot in Fig. 2. The error bands are obtained by varying the starting point of the RG evolution in both m 0 (τ ) extraction and perturbative matching.The results show much reduced error bands from LRR because of the much smaller uncertainty in m 0 (τ ).Interestingly, the NNLO+RGR+LRR result suggests a even smaller error after matching, because the scale variation in the RGR matching cancels most of the corresponding m 0 difference in coordinate space.Moreover, the consistency in x > 0.2 between NLO and NNLO suggests good convergence of the perturbation theory after LRR, the same as our observation in coordinate space.In conclusion, we made a first systematic study of the leading power accuracy in LaMET calculations of PDFs and suggested an approach to make practical progress.We show the importance to define the renormalization scheme for the linear divergence consistently in the renormalization and the perturbative matching.To achieve that, we resum the leading IRR in the the perturbative Wilson coefficients.The leading IRR is universal in all orders of Wilson coefficients.Thus after renormalization in the same scheme, the P z = 0 lattice correlators are also consistent with perturbative results up to linear accuracy.Then P z > 0 PDF extraction allows the same level of precision, provided a perturbative matching in the same scheme is applied.An application to the pion PDF shows that our approach significantly reduces the uncertainty from linear corrections.This is a necessary step towards the future high-precision calculation of PDFs on lattice.

FIG. 1 :
FIG. 1: Top: Uncertainty in m0(τ ) extracted from the pion Pz = 0 matrix element and fixed-order Wilson coefficients with and without RG resummation.The band width shows the renormalizaiton scale µ dependence.Bottom: The uncertainties in the pion light-cone PDF extracted from LaMET expansion with the above m0.The overlap region exhibits a darker color.

FIG. 3 :
FIG. 3:The effect of leading-renormalon resummation (the red and blue band) on the pion PDF, compared with fixedorder results in the background.