Power corrections and renormalons in parton quasi-distributions

Perturbative expansions for short-distance quantities in QCD are factorially divergent and this deficiency can be turned into a useful tool to investigate nonperturbative corrections. In this work, we use this approach to study the structure of power corrections to parton quasi-distributions and pseudo-distributions which appear in lattice calculations of parton distribution functions. As the main result, we predict the functional dependence of the leading power corrections to quasi(pseudo)-distributions on the Bjorken $x$ variable. We also show that these corrections can be strongly affected by the normalization procedure.

this case, the analytic structure of higher-twist effects is well understood, so they could have been modelled by just one free parameter. The general situation may be more complicated. Having in mind the multitude of observables that can potentially be employed in lattice calculations, see e.g. [6], it is desirable to have a general method to estimate the corresponding higher-twist corrections and their expected Bjorken-x dependence. The purpose of this work is to point out that the problem of power-suppressed contributions for such observables can be addressed using the concept of renormalons [28,29].
In what follows we apply this technique to qPDFs and pPDFs.
The renormalon approach to the investigation of power corrections is founded on the fact that operators of different twist mix with each other under renormalization, due to the violation of QCD scale invariance through the running of the coupling constant. In cutoff schemes, this mixing is explicit, whereas in dimensional regularization, it manifests itself in factorial divergence of the perturbative series. Independence of a physical observable on the factorization scale implies intricate cancellations between different twists -the cancellation of renormalon ambiguities. In turn, the existence of these ambiguities in the leading-twist expressions can be used to estimate the size of power-suppressed corrections. Conceptually, it is similar to the estimation of the accuracy of fixedorder perturbative results by the logarithmic scale dependence. The renormalon approach was used before for the study of Bjorken-x dependence of higher-twist corrections in deep-inelastic scattering [30][31][32], fragmentation functions [33,34], pion LCDA [35] and transversemomentum-dependent parton distributions [36].
In order to explain how the concept of renormalons can be used to get insight in the structure of power corrections, let us consider the usual expression for the quasidistribution [22,37], where q(y, µ F ) is the quark PDF, p is the hadron mo-arXiv:1810.00048v1 [hep-ph] 28 Sep 2018 mentum, x refers to the momentum fraction. For brevity we do not show the dependence on the renormalization scale. The factorization scale µ F has to be taken of the order of x p to avoid large logarithms. The coefficient function C(x, p, µ F ) = δ(1 − x) + O(α s ) is given by the perturbative expansion. The correction of O(α s ) was first computed for the flavor-nonsinglet case in [38], see also [39].
To understand the role of renormalons, it is necessary to examine carefully the separation made in (1) between the leading term and the higher twist addenda. Let us assume for a moment that the factorization is done using a hard cutoff Λ QCD ≪ µ F ≪ p, i.e. the contributions with loop momenta k > µ F are included in the coefficient function, whereas the contributions with k < µ F are included in PDF. In this scheme, the coefficient function has the following expansion at p → ∞ where c k = c k (x, ln p 2 µ 2 F ) are the perturbative coefficients depending logarithmically on the scales and the D Q -term represents the leading power correction. Since the left-hand-side of Eq. (1) does not depend on µ F , any such dependence should cancel on the right-hand-side. In particular, the logarithmic dependence on the scale in c k (x, ln p 2 µ 2 F ) is canceled by the scale-dependence of PDF q(x, µ F ). The cancellation of the power dependence, on the other hand, must involve the twist-four contribution Q 4 (x, p). Thus, in this factorization scheme one expects that whereQ 4 depends on µ F at most logarithmically. Appearance of the term ∼ µ 2 F can be traced to quadratic ultraviolet divergence (in addition to the logarithmic ultraviolet divergence) of the twist-four operators that are responsible for the power correction, in the cutoff scheme. One can prove that the cutoff dependence ∼ µ 2 F of the higher-twist operators is indeed that of Eq. (2).
In practice, perturbative calculations are usually done using dimensional regularization. In this case, powerlike terms as in (2) do not appear. The price to pay is that the coefficients c k computed in a MS scheme grow factorially with the order k. The factorial growth implies that the sum of the perturbative series is only defined to a power accuracy and this ambiguity (renormalon ambiguity) must be compensated by adding a nonperturbative higher-twist correction. The detailed analysis shows [28,29] that the divergent large-k behavior (the renomalon) of the coefficients is in the one-to-one correspondence with the sensitivity to extreme (small or large) loop momenta. In particular, infrared renormalons in the leading-twist coefficient function are compensated by ultraviolet renormalons in the matrix elements of twist-four operators. In this way the same picture as in the cutoff scheme re-appears in dimension regularization.
Returning to (3), we observe that the quadratic term in µ F is spurious since its sole purpose is to cancel the similar contribution to the coefficient function. Therefore, it does not contribute to any physical observable. The idea of the renormalon model of the power corrections [30][31][32][33][34][35] is that, with a replacement of µ F by a suitable non-perturbative scale, this contribution reflects the order and the functional form of actual powersuppressed contribution. Assuming this "ultraviolet dominance" [28,29,40] one obtains the following model: with the dimensionless coefficient κ = O(1) which cannot be fixed within theory and remains a free parameter.
In this work we calculate the function D Q (x) for different versions of the qPDFs (pPDFs) in the so-called bubble-chain approximation [28], which is our main result. This calculation reveals that the power corrections to quasi-and pseudo-distributions have the following generic structure , respectively, and can be affected significantly by normalization. In particular, the normalization of the involved matrix elements to their value at zero momentum considerably reduces the power correction for the qPDFs at smaller-x at the cost of a strong enhancement at larger values. We emphasize that the leading power correction to qPDF at x → 0 is enhanced by two powers of the Bjorken x variable. For pPDFs the power corrections are suppresed at x → 1, which, unfortunately, does not hold after the normalization procedure of Ref. [24] (but can be upheld with a different choice). Additionally, as a byproduct of bubble-chain calculation, we have obtained the large-n f part, n k f α k+1 s , of the leading twist coefficient function to all orders in perturbation theory.
The presentation is organized as follows. In Sect. 2 we formulate our program in more precise terms. We define qPDFs and pPDFs as particular Fourier transforms of the position space (Ioffe-time) distributions, discuss briefly the light-ray operator product expansion (OPE) and the target mass corrections, and introduce the relevant techniques (Borel transform) and the systematic approximation (large-n f expansion) that will be used throughout the rest of the work. In Sect. 3 we present our result for the Borel transform of the leading-twist coefficient function and discuss the structure of its singularities. The leading power corrections to various versions of the quasi-distributions are obtained in Sect. 4. We collect there the relevant analytic expressions and also present the results of a numerical study using realistic parametrizations for the valence quark PDFs. The final Sect. 5 is reserved for the summary and conclusions.

A. Parton quasi-distributions
Let us start with the following nucleon matrix element where p µ is the nucleon momentum, p 2 = m 2 , v µ is a given four-vector and z a real number. All notations correspond to Minkowski space. In the following, we keep the normalization of the four-vector v µ arbitrary, keeping in mind that v 2 < 0 in the lattice calculation. We suppress flavor indices and tacitly assume considering the flavornon-singlet combination of the matrix elements in what follows. We also neglect quark masses. The operator product in Eq. (6) suffers from ultraviolet (UV) divergences and has to be renormalized. The argument µ of matrix elements I refers to the renormalization scale. The renormalization-scale dependence can be studied by going over to an effective theory [41][42][43] such that the Wilson line is substituted by the propagator of an auxiliary field. For time-like separations the resulting theory is the heavy-quark effective theory (HQET) and the renormalization factor of interest is the renormalization factor (squared) for the heavy-light quark current, see e.g. [44]. We are not aware of a calculation for such a renormalization factor at space-like separations beyond one-loop order, however, to this accuracy there is no difference from the time-like case.
The "longitudinal" and "transverse" invariant functions in Eq. (6) (with respect to v µ ) correspond to particular projections of the matrix element that are employed in lattice calculations: where ( · v) = 0. Eq. (8) is the starting point for the construction of quasi-distributions. The invariant functions I ∥ , I ⊥ (6) coincide at the tree level. Assuming the power counting they can be written in terms of the position-space quark PDF [45][46][47]: where µ F is the factorization scale, and Here q(x, µ F ) for x > 0 is the quark PDF, and for x < 0 is the antiquark PDF. The position-space PDFs I are known as the Ioffe-time distribution (ITD). To distinguish ITD I from the functions I, we refer to the functions I ∥ , I ⊥ as the Ioffe-time quasi -distributions (qITDs), following the terminology introduced in [5]. The qPDFs Q ∥ , Q ⊥ [5] and the pPDF P [23] are defined in terms of the qITDs by the Fourier transform In what follows we will tacitly assume that (pv) > 0 in the discussion of qPDFs and z > 0 for pPDFs and drop the absolute value sign, for brevity. In renormalization schemes with an explicit regularization scale, the Wilson line in Eq. (6) suffers from an additional linear UV divergence that has to be removed. In dimensional regularization, this UV linear divergence reveals itself as a factorial growth of high orders of perturbative series [48]. The linear UV divergence can be removed by a mass renormalization associated with the Wilson line [41,42,49] or by a regularization-independent renormalization [14,39,50]. Given the multiplicative renormalizability of the quasi-PDF operator [43,51,52], it can also be removed by considering a suitable ratio of matrix elements involving the same operator, e.g. by normalizing to the same matrix element at zero proton momentum [24], or, alternatively, to the vacuum expectation valuê where Eqs. (13) and (14) will be our main focus in the present paper. By forming the ratio, the scale dependence cancels out (including the usual logarithmic renormalization) and one can define the scale-independent qPDF/pPDF and similarly forQ(x, p) andP(x, z). The difference between these two options in the present context is that the normalization to the vacuum correlator does not affect the leading O(v 2 z 2 ) power corrections that are subject of this work (since there is no gauge-invariant operator), whereas the normalization to the value at zero momentum, as we will see, has a substantial effect.

B. Light-ray OPE and target mass corrections
The general approach to collinear factorization of QCD amplitudes in the position space is provided by the lightray OPE [53][54][55][56][57][58]. For illustration consider the "longitudinal" projection. Specializing to the present case (forward matrix elements) we writē where we include the renormalization factor in the coefficient function H ∥ = δ(1−α)+O(α s ) and set the renormalization scale µ to be equal to the factorization scale µ F .
is the leading-twist projection operator defined below and ellipses stand for the higher-twist contributions.
The leading-twist projection of the nonlocal quarkantiquark operator is defined as the generating function of renormalized local leading-twist operators (traceless and symmetrized over all indices) where Here we indicate trace subtraction and symmetrization by enclosing the involved Lorentz indices in parentheses, The light-ray OPE differs from the usual short-distance Wilson expansion in local operators by imposing a different power counting. In the short-distance expansion one assumes that the distance between the quarks is small, z ∼ ηΛ −1 QCD with η → 0, and the operator matrix elements are of order unity in this limit, ⟨O n µ1...µn ⟩ ∼ Λ n QCD . In this case only a finite number of terms contribute to the r.h.s. of Eq. (18), whereas the rest as well as the higher-twist operators must be added to higher orders of OPE, starting from O(η 2 ). In other words, the relevant expansion parameter is the mass-dimension of an operator. The light-ray OPE assumes instead that the leading-twist operators scale as ⟨O n µ1...µn ⟩ ∼ η −n Λ n QCD so that z n v µ1 . . . v µn ⟨O n µ1...µn ⟩ = O(1). In this case, the series in (18) must be resummed to all orders, revealing its non-local "light-ray" structure. For a generic hadronic matrix element of leading twist operators one has where the reduced matrix element ⟨⟨O n ⟩⟩ = O(1). Therefore, the light-ray OPE is an adequate approximation if the hadron has large momentum, pv = O(η −1 ) and hence z pv = O(1). Higher-twist operators of the same dimension have smaller spin (by definition) and as a consequence, their matrix elements are power-suppressed [59]. Note that the above power counting is applicable both in Minkowski and Euclidean space. In Minkowski space, one can go over to a different reference frame where all components of the momentum are of order Λ QCD and simultaneously the separation between the quarks is al- On the calculation level, the light-ray OPE provides one with a convenient framework to operate with the leading-twist projected operators (18) avoiding the local expansion. Light-ray operators can be viewed as analytic operator functions of the separation between the quarks (all short-distance and light-cone singularities are subtracted). They satisfy the Laplace equation [ with the boundary condition Π µ F l.t. → 1 at v 2 → 0. Explicit expressions for the projection operator Π µ F l.t. can be found in [55,58,60,61]. The light-ray OPE combined with the background field method is the standard technique, e.g., in light-cone sum rules [62], where it is used for the calculation of higher-twist contributions, and for the derivation of the evolution equations for off-forward parton distributions [63,64]. The renormalization group kernel for the evolution of flavor-singlet light-ray operators in the general off-forward kinematics is known to three-loop accuracy [65].
The nucleon matrix element of the leading-twist projected operator (18) defines the leading-twist quark PDF, where [57] The second term in the square brackets is the leading nucleon mass correction, which is the position-space counterpart of the Nachtmann target mass correction [66] in the deep-inelastic scattering (DIS). The all-order expression in powers of the nucleon mass for the leading-twist projected exponential function can be found in [57]. Taking the nucleon matrix element of the operator relation in Eq. (17) we obtain the factorization theorem for the "longitudinal" qITD in the form where I(αzpv, µ F ) is ITD (11) and ellipses stand for the higher-power target mass and "genuine" higher-twist corrections.
Making the Fourier transformation (12) one obtains, to the leading twist accuracy, the factorization theorem for the "longitudinal" qPDF, where the coefficient function is given by In particular, at the leading order where q ′ (x) = (d dx)q(x). Note that the derivative applied to the quark PDF lowers the power q(x) x→1 ∼ (1−x) p by one unit so that the mass correction is effectively enhanced by the factor 1 (1 − x) as compared to the quark distribution itself at large Bjorken x. The similar enhancement of the target mass correction at x → 1 is familiar from DIS [66].
The target mass correction for the "transverse" qPDF can be calculated in a similar way, starting from the nonlocal operator with an open Lorentz index (6) and using the operator identity [55] where the Wilson line between the quarks is implied. One obtains, at the tree level, and the target mass correction to the pPDF Interestingly, it is suppressed as O(1 − x) at x → 1 and not enhanced O (1 (1 − x)) in contrast to the target mass correction to the qPDFs (27) and the structure functions in DIS.

C. Borel transform and renormalons
The coefficient function H ∥ in Eq. (17) and the similar coefficient function H ⊥ in the MS scheme have the perturbative expansion with factorially growing coefficients h k ∼ k!.
A convenient way to handle such a series is to consider the Borel transform where powers of β 0 = 11 3N C − 2 3n f are inserted for the later convenience. The Borel image can be used as a generating function for the fixed-order coefficients Moreover, the sum of the series can be obtained as the integral over positive values of the Borel parameter w As it stands, the integral is not defined because the Borel transform generally has singularities on the integration path, known as (infrared) renormalons. One can adopt a definition of the integral deforming the contour above or below the real axis, or as the principle value. These definitions are arbitrary, and their difference, which is exponentially small in the coupling, must be viewed as an intrinsic uncertainty of perturbation theory that has to be removed by adding power-suppressed nonperturbative corrections. Another potential problem concerns the convergence of the Borel integral at w → ∞. Since the quantity of interest depends on the single hard scale 1 zv , dimension counting requires that the Borel transform can be written as (µ 2 z 2 v 2 ) w times a function F (w) of the Borel parameter and dimensionless kinematic variables. Combining (µ 2 z 2 v 2 ) w with e −w (β0as(µ)) one sees that the (principal value) integral is convergent, provided the distance between the quarks zv is sufficiently small compared to 1 Λ QCD and F (w) does not increase exponentially at w → ∞, which is the case of all known examples. Naturally, a full all-order calculation cannot be performed. Instead, we employ the approximation [28,29] restricting ourselves to the perturbative series generated by the running-coupling effects in the one-loop diagrams, i.e. using QCD coupling at the scale of the gluon virtuality. Such contributions can be traced by computing the diagrams with the insertion of k fermion loops in the oneloop diagram and replacing − 2 3 n f ↦ β 0 = 11 3 N c − 2 3 n f , see Appendix A. Another, equivalent technique [67] is based on the calculation of one-loop diagrams with an effective gluon mass.
The singularity structure of the Borel transform can be extracted separately without explicit evaluation of the bubble-chain. It can be done by replacing the running coupling constant in the loop diagrams by its effective form, where Λ ≡ Λ MS QCD and the factor e At one-loop level, the pole structure in w reproduces the pole structure for the Borel plane. In our work, we have used both methods of calculation for the cross-check of the result.

III. LARGE-β0 COEFFICIENT FUNCTIONS
The leading contributions to the renormalon singularities in the coefficient functions for qITDs are shown in Fig. 1, see Appendix A for technical details. We obtain and the functionR(w) is defined as the series expansion in terms of another function such that The Taylor expansion of the Borel transform at w = 0 gives the perturbative expansion for the coefficient functions in terms of the coupling constant. The O(α s ) term is and the n f -part of the O(α 2 s ) correction reads where We have checked that the expressions in (42) reproduce after the Fourier transform (12) the one-loop correction to the qPDFs calculated in [38,39].

A. Singularities of the Borel transform
The structure of singularities of the Borel transform of the coefficient functions is illustrated in Fig. 2. There are an UV renormalon singularity at w = 1 2 and a series of IR renormalons at positive integer w = 1, 2, . . .. To the accuracy of our calculation (single bubble chain) all singularities are simple poles. These singularities obstruct the Borel integral in (34) (Borel-non-summable renormalons) and must be matched by the nonperturbative corrections. Note, that theR(w) term in (37) is an analytic function of w. Thus, it is irrelevant for the discussion of singularities.

Ultraviolet renormalon at w = 1 2
The singularity at w = 1 2 is generated by the contribution of large momenta in the self-energy insertions in the Wilson line and is part of the renormalization factor This singularity is well-known [48] and is in the one-toone correspondence to the linear UV divergence in the Wilson line's self energy (see also discusion prior to (13)). It can be removed by considering normalized qPDFs, (13) or (14), and will not be considered further in this work.

IV. POWER CORRECTIONS
A singularity on the integration path in Eq. (34) means that the perturbation theory is incomplete and the sum of the series is ill-defined. It is customary [28] to estimate the corresponding ambiguity as where w 0 is the position of the singularity and Res w=w0 B[H](w) is the corresponding residue. Note that e −w0 (β0as) = (Λ 2 µ 2 ) w0 . Following the standard logic [28,29] we assume that this ambiguity must be canceled by adding a non-perturbative correction of the same order of magnitude.

A. Ioffe time quasi-distributions
Considering δH (1), we obtain the leading power correction to the qITDs as functions of the "Ioffe-time" where κ is a real number of order one. The renormalon ambiguity (48) corresponds to but this number is only indicative. Alternatively, one can put κ = 1 and think of Λ as a certain nonperturbative parameter of the order of Λ QCD that determines an overall normalization of the power correction and cannot be fixed in this approach. Note that also the sign of the correction is not determined. For the normalized qITDs defined in Eq. (13) we obtain instead where the "plus" distribution is defined as usual, The leading power corrections to the qITDsÎ ∥(⊥) (14) that are normalized to the vacuum correlator, are the same as for the un-normalized distributions (52). The expressions in (52) and (54) present the starting point for the following analysis. In order to visualize the functional dependence of the power correction on the "Ioffe time" τ relative to the leading-twist result, we write where For the normalized qITDs one obtains and, obviously,R so that we do not consider them separately. In general, the qITDs I(τ ) and the higher-twist coefficients R(τ ) and R I (τ ) are complex functions, but their imaginary parts appear to be small. The real parts, Re R I (τ ) and Re R I (τ ), for a simple model of the valence quark distribution are plotted in Fig. 3. The power correction to the "transverse" qITD turns out to be roughly factor four larger than for the "longitudinal" qITD. In both cases the Rfunctions flatten out at large Ioffe times τ ≥ 10 that are, however, hardly accessible in present day lattice calculations. The normalization to zero proton momentum corresponds to the subtraction of the value at τ = 0, so that for z → 0 there is no power correction by construction.

B. Parton quasi-distributions
Making the Fourier transformation of the above results for the qITDs we obtain the qPDFs and Assuming for what follows x > 0, these expressions can be rewritten in the form with Note that we have extracted the prefactor 1 (x 2x ) for the power correction anticipating that it is enhanced as 1 x 2 and 1 (1 − x) in the regions of small x → 0 and large x → 1 Bjorken variable, respectively. The normalized QPDFs Q ∥(⊥) (x, p) are obtained by replacing the kernels in the first lines in Eqs. (61) and (62) by the plus distributions, cf. (54). Writing the result in the form whereq(x) is the quark PDF normalized to the unit integral, one obtains Note that the additional terms are proportional to the second derivative of the quark PDF and thus enhanced as 1 (1 − x) 2 at x → 1. As already mentioned, the normalization to the vacuum correlator does not affect the 1 p 2 power corrections so thatR Q (τ ) and we do not need to consider this case separately.
For a numerical study we have used the MSTW NLO valence u-and d-quark distributions [68] at the scale 2 GeV and the simple model in Eq. (60). It turns out that the power corrections for (unnormalized) "longitudinal" and "transverse" qPDFs are similar in size so that we show R ∥ Q (x) on the left panel in Fig. 4, and the difference R ⊥ Q (x) − R ∥ Q (x) on the right panel. The x-dependence of R Q (x) is similar for all quark PDF models: The power correction is small for x → 0 (but non-zero, R ∥ Q (0) ∼ 0.2, R ⊥ Q (0) ∼ 0.4) and increases steeply with x (almost linearly). For the uquark, the result is very similar to the simple model in Eq. (60) whereas for the d-quark the power correction is, roughly, factor two larger. Note that the difference depends on the PDF model only very weakly.
Constructing the qPDFs from the qITDs normalized to the value at zero momentum has a large effect. This is illustrated in Fig. 5 where in the upper panels we show the results for the "longitudinal" case and in the lower panels the difference between the "longitudinal" and "transverse" distributions. The three panels (from left to right) correspond to the MSTW valence u-quarks, d-quarks, and the simple model in Eq. (60), respectively. The red solid lines stand for R Q (x) (i.e. the same as in Fig. 4), and the result for the normalized qPDFs, R Q (x), is shown by the black dashed curves.
We see that the normalization to the qITD at zero momentum significantly reduces the power correction at moderate values of x ≲ 0.5 at the cost of dramatical increase at higher x values. This normalization procedure is, therefore, not suitable to access the large-x behavior of the PDFs, but apparently minimizes power corrections in the not-so-large x region.

C. Parton pseudo-distributions
Power corrections for the pPDF can be obtained easily from the corresponding expression for the "transverse" qITDs (52), (54). Writing the result as and similarly for the pPDF normalized to zero momentum, P(x, z), we obtain The numerical results are shown in Fig. 6. Note that R P (x) is very similar for all considered models for the valence quark PDFs and decreases at x → 1. Indeed, it is easy to see that R P (x) = O(1 − x) in this limit, similar to the target mass correction, Eq. (30). This suppression is removed once pPDF is normalized to zero momentum [24] (which adds a negative constant), but can be upheld if normalized to the vacuum expectation value of the same operator.

V. CONCLUSIONS
We have presented an analysis of power-suppressed (higher-twist) contributions to qPDFs and pPDFs based on the study of factorial divergences (renormalons) in the corresponding coefficient functions within the bubblechain approximation. Factorial asymptotic implies that the sum of the series is only defined to a power accuracy and therefore, the QCD perturbation theory must be corrected by nonperturbative power-suppressed contributions to produce unambiguous predictions. Our results have to be considered as a "minimal model" for the higher-twist corrections that captures effects that are necessary for the selfconsistency of the theory, but possibly misses other nonperturbative corrections that are, e.g, protected by symmetries and not "seen" through perturbative expansions. Our main conclusions are as follows: • Position space PDFs (qITDs) have flat power corrections at large Ioffe times. Generally, power corrections are much larger for the "transverse" projection as compared to the "longitudinal" projection.
• Power corrections for qPDFs have a generic behavior Note that the corresponding target mass corrections, (27) and (29), do not show up the 1 x 2 enhancement. This behavior is commensurate to the suppression of Nachtmann's target mass corrections ∼ x 2 m 2 Q 2 at small x for the DIS structure functions. The normalization of the underlying qITDs to unity at the zero momentum considerably reduces the size the power correction to the qPDFs at x ≲ 0.5 at the cost of a strong enhancement at x ≳ 0.5 − 0.6. Thus, such a normalization procedure is not suitable to access the large-x behavior of the PDFs, but apparently minimizes power corrections in the intermediate x-region. Note that the above discussion is based on a normalization procedure different from the nonperturbative renormalization used in ref. [17], thus the power corrections for the latter might behave differently from what has been shown here.
• Power corrections for pPDFs have a generic behavior The suppression at x → 1 is lifted by the zeromomentum normalization factor. However, it can be upheld by the normalization to the vacuum matrix element of the same operator. We conclude that pPDFs can offer an interesting alternative to qPDFs for the study of large-x region.
As a byproduct of this study, we provide the results for the coefficient functions in the large-n f approximation, terms ∼ n k f α k+1 s , to all orders in perturbation theory. These results can be useful to estimate the effects of uncalculated higher orders and scale-setting using the BLM-type procedure. The corresponding analysis goes beyond the scope of this paper and will be presented in a future publication. with 1 ⩽ m ⩽ n−1. This equation reflects the fact that the composition of two bubble-chains is again a bubble-chain. The -poles in ∆ n correspond to UV poles of fermion loops, and to be renormalized later.
The leading large-n f diagrams structurally reproduce the one-loop diagrams calculated in the Landau gauge. Due to the gauge invariance, the transverse part of the gluon propagator cancels in the sum of diagrams, and therefore, the complete result can be obtained by the use of only longitudinal part, which significantly simplifies the calculation. The reduced propagator reads ∆ µν n = −1 4π d 2 Γ(1 − (n + 1) ) Γ(1 + n ) The expressions for diagrams evaluated with this propagator are (notation corresponds to Fig.1) .
The renormalization of these diagrams is straightforward and discussed in details e.g. in ref. [30]. Next, we present the diagram-by-diagram renormalized expression for the coefficient function. One obtains r n+1 (ᾱ) + r n (ᾱ) n + 1 + (−1) n n!g 0 (γz 2 ) O µ v (ᾱz, βz) + 2(−1) n n!g where z 2 = −z 2 v 2 µ 2 4, γ = 1 − α − β. In these expressions we have promoted the factors (−2n f 3) to the β 0 coefficient of QCD, β 0 = 11 3C A − 2n f 3. The functions g Note, that operators on the r.h.s. of (A10-A13) are renor-malized operators. Taking the hadronic matrix element of these diagrams and performing the Borel summation we obtain the final result given in Eq. (37). The factorial divergences ∼ β n 0 n! that correspond to Borel non-summable renormalons can be inferred directly from the expressions in (A10-A13) by inspection. These terms give rize to the renormalon contribution in (37). The remaining terms are regular and contribute to the functionR(w) in Eq. (37). Finally, specializing to the particular n values one obtains the ∼ a n+1 s n n f contributions to the coefficient function. At n = 0 one obtains the NLO coefficient function, quoted in (42). The n = 1 result corresponds to the n f part of the NNLO coefficient function, Eq. (44).