Renormalon cancellation and linear power correction to threshold-like asymptotics of space-like parton correlators

In this paper, we show that the common hard kernel of double-log-type or threshold-type factorization for certain space-like parton correlators that arise in the context of lattice parton distributions, the heavy-light Sudakov hard kernel, has linear infrared (IR) renormalon. We explicitly demonstrate how this IR renormalon correlates with ultraviolet (UV) renormalons of next-to-leading power operators in two explicit examples: threshold asymptotics of space-like quark-bilinear coefficient functions and transverse momentum dependent (TMD) factorization of quasi wave function amplitude. Theoretically, the pattern of renormalon cancellation complies with general expectations to marginal asymptotics in the UV limit. Practically, this linear renormalon explains the slow convergence of imaginary parts observed in lattice extraction of the Collins-Soper kernel and signals the relevance of next-to-leading power contributions. Fully factorized, fully controlled threshold asymptotic expansion for space-like quark-bilinear coefficient functions in coordinate and moment space has also been proposed.

3 Renormalon cancellation and power corrections to threshold asymptotics of space-like quark-bilinear coefficient functions 3.1 Collection of full results before threshold limit 3.2 Threshold limit in coordinate and moment space 3.3 Borel transform and conservation of "genuine" IR renormalons in [1] 3. 4 The NLP threshold soft factor and its UV renormalon 3.5 Exponentially small terms and resurgence of the full distribution from threshold asymptotics 3.6 A conjecture on threshold limit of hadronic quark-bilinear matrix element It is known that perturbative series induced by marginal perturbations to UV or IR conformal field theories (CFTs) are tricky and subtle.On one hand, fixed-order Feynmandiagram-based calculations suffer from uncontrolled large logarithms that must be resummed using renormalization group equation (RGE) or whatever structures to become controlled asymptotic expansions.On the other hand, even after resummation, normally the perturbative series can only be a divergent, Borel non-summable asymptotic series.Physically, this reflects the fact that in the presence of marginal perturbation, the decoupling between UV and IR scales is very slow, and the Borel non-summability is the price one pays for decoupling asymptotically.
In fact, Borel non-summable logarithmic asymptotic series is not unusual in nature.A famous example from mathematics is the prime number theorem (PNT), stating that the asymptotic distribution π(x) of prime numbers takes the form in the x → ∞ limit, which is similar in shape to a perturbative series with "one-loop running coupling constant" α(x) = 1 ln x .One can check that the Borel resummation of the above series has a pole at t = 1, and different prescriptions, such as the principle value or ±i0 prescriptions lead to O(1) ambiguity Can one read from this ambiguity any useful information regarding "power corrections" to the PNT?The answer is unfortunately no.In fact, one knows that the true "next-toleading-power" correction comes from non-trivial zeros of the Riemann zeta function and is of the order √ x (assuming Riemann hypothesis) In such a case, the Borel singularity at "leading power" carries no useful information regarding the power corrections at all.
Given this example, one may ask the following question: near UV fixed point of local quantum field theories (QFTs), is it also true that marginal perturbative series for IR safe quantities, even with their singularity structures in the Borel plane completely decoded, still can not tell where the non-trivial power corrections are?The answer to this question is fortunately no.In fact, although the perturbative series in QFT is usually much more complicated than Eq.(1.1) and in prior one can not even see where the Borel singularity is, due to the special context of local QFTs, important understandings can still be formed.In particular, the perturbative renomalon [2][3][4] is widely believed to be a crucial structure that can provide useful information regarding the high-order behavior of perturbative series as well as power corrections.
Based on our understanding, positive arguments and evidences in favor of perturbative renormalons can be summarized as below.
1. First of all, they are caused by a small amount of diagrams (such as the bubble chain diagram for QCD) which allow rigorous analysis as well as explicit calculations to show that they indeed lead to non-alternating factorial growths of the form ( β 0 2k ) n n! which would become Borel singularity at t = 2k β 0 if not miraculously canceled by other diagrams.
2. Second, at the level of Feynman integrals, the factorial growth is caused by running inside loop integrals which exactly drives the loop momenta towards the "dangerous direction", such as the small momentum region in QCD, suggesting natural correlation with power corrections or non-perturbative contributions.
3. Third, in several examples where the high power contributions are strongly believed to be caused by high-dimensional operators in the operator product expansion (OPE), a characteristic structure for CFTs that are very likely to be inherited in some form in their marginal deformations as well, the IR renormalon in the leading power can be shown to correlate [5][6][7][8][9][10][11][12][13] with UV renormalons caused by ambiguities in defining these high dimensional operators.
4. Fourth, in all the integrable none-supersymmetric asymptotically safe 2D massive QFTs that have been analyzed using the thermodynamical Bethe's ansatz (TBA) method so far [14][15][16][17][18][19][20][21], the leading singularities in the right-half Borel plane (for the perturbative series for ground state energy in the large "chemical potential" limit) and the place of NLP contribution all agree 1 with predictions based on leading perturbative IR renormalon.In particular, for the O(N ) non-linear sigma models, both the large-N expansion [6,7,11,23,24] at O( 1 N ) and the exact TBA method agrees [16,25,26] and works in favor of perturbative renormalon. 5. Finally, for the lattice perturbation series for the self-energy of a Wilson-line, the asymptotic behavior extracted from high-order numerical calculations agrees with the renormalon-based prediction [27,28].
To summarize, the perturbative renormalon, in particular, the bubble-chain-based renormalon analysis for QCD (sometimes called "large β 0 approximation" in literature) is a valuable method to determine leading Borel singularity of QCD perturbative series as well as to probe, or at least, to provide lower-bounds for high power corrections.
In this paper, based on this understanding, we investigate the Borel singularity using the bubble-chain diagrams, for a particular hard kernel that appears in several factorization formulas in the large momentum effective theory (LaMET) formulation [29][30][31][32] for lattice parton distributions.This hard kernel first appeared in the factorization of quasi-TMDPDFs [33][34][35][36][37][38], and then in TMD factorization of quasi light-front wave function (LFWF) amplitudes [39] and recently in threshold limit of space-like quark-bilinear coefficient functions [40].It is a double-logarithmic type or threshold-type hard kernel in the sense that it satisfies an RGE with a logarithmic term in its anomalous dimension of the form Γ cusp (α) ln ζz µ 2 controlled by the universal light-like cusp anomalous dimension Γ cusp (α) [41][42][43] .Moreover, its definition is very simple: it is just an on-shell form factor defined purely in PQCD with UV and IR divergences subtracted multiplicatively in the MS scheme, with an incoming on-shell light-quark with four-momenta p 2 = 0 and an outgoing on-shell heavy Wilson-line in space-like direction n 2 z = −1.The Lorentz invariant combination ζ z = 4(p • n z ) 2 ̸ = 0 serves as the only scale of the form factor.For this reason, it is called the heavy-light Sudakov form factor [40] or heavy-light Sudakov hard kernel due to its similarity with the standard light-light Sudakov form factor [44][45][46][47] with on-shell external light-partons.On the other hand, since the gauge link can be interpreted as the gauge transformation that transforms a generic gauge to the axial gauge , this form factor can also be regarded as the gauge-invariant definition of the quark field renormalization factor in the axial gauge, for which double-logarithmic asymptotics is known long ago [48].In all three factorization formulas, this form factor is caused by the Sudakov effect around the quark-link vertices and has recently been determined up to NNLO [40,49].
Naively looking, this hard kernel should contain only quadratic IR renormalon singularities which is the usual case for hard kernels.On the other hand, the presence of a heavy Wilson-line in its definition resembles quantities in heavy-quark effective theory, for which linear renormalon is common.It is the purpose of this paper to show using explicit calculations that the heavy-light Sudakov form factor contains linear IR renormalon in its phase angle.Moreover, we show that for both the threshold factorization of space-like quark correlator as well as the TMD factorization of quasi-LFWF amplitude, this linear renormalon cancels with corresponding UV renormalons for operators at next-to-leading power (NLP).As such, these examples can also serve as simple demonstration-of-principle for the general pattern of renormalon cancellation in the context of double-logarithmic-type asymptotic expansions.The case of quasi-LFWF also provides a simple example of the "OPE nonconvergence" phenomenon [50,51].We hope that these examples will be not only useful for the lattice extraction of partonic distributions and evolution kernels, but also interesting for their pedagogical value to demonstrate the renormalon cancellation/correlation pattern mentioned before.
We emphasize that the linear renormalon in the heavy-light Sudakov hard kernel requires "on-shell" space-like gauge-link and affects only the double-logarithmic type factorization of quasi-distributions.Before power expansions that lead to double-log evolution, the gauge links "cancel among themselves" at infinity between "real" and virtual diagrams, and the linear renormalons are hidden.After the TMD expansion in 1 b ⊥ p z or threshold expansion in 1  zp z , at leading power the cancellation is no-longer perfect, and the linear renormalon can only be cancelled between leading and higher powers.
The organization of the paper is as follows.
1.In Sec. 2, we calculate explicitly the bubble chain diagram for the hard kernel in the MS scheme following the method of [52][53][54][55].We show that the imaginary part suffers linear renormalon ambiguity, while the real part is free from linear renormalon.Furthermore, the strength for the linear renormalon ambiguity happens to be the same as that for the linear divergence of pole mass.As a by-product, we also obtained the corresponding MS scheme anomalous dimensions for the bubble chain diagram.
2. In Sec. 3, we investigate the role played by this linear renormalon in the threshold asymptotics of space-like "quasi" [29,30] or "pseudo" [56,57] quark-bilinear coefficient functions in coordinate space and moment space.We first review the factorization formalism, emphasizing the UV nature of the threshold limit in coordinate space and moment space at fixed z 2 .A fully-factorized controlled expansion in the N → ∞ limit has been written down explicitly.
We then use explicit analytical result of bubble chain diagram [1] to perform the threshold expansion both in coordinate and moment space.The threshold factorization formula at leading threshold power is verified at the bubble chain level.The linear renormalon cancels between the heavy-light hard kernel at leading threshold power and the NLP threshold soft factor.We identify the operator definition for the NLP contribution and demonstrate that its linear renormalon is UV.We show that the "genuine" IR renormalons related to zΛ QCD power corrections remain the same before and after threshold expansion.We further investigate the "resurgence" of full distribution from threshold asymptotics by noticing the correlation between Borel ambiguity of e −iλ terms (corresponding to α → 1 threshold asymptotics) and branch ambiguity of purely algebraic terms (corresponding to α → 0 asymptotics).
3. In Sec. 4, we perform the same analysis to the TMD factorization of quasi-LFWF amplitude.We show that the bubble chain diagram for the quasi-LWFW amplitude, before the TMD expansion, is free from the linear renormalon.After performing the power expansion, the linear renormalon in the heavy-light Sudakov form factor, playing the role of the hard kernel at leading power, cancels with the UV renormalon of NLP soft contribution.We determine the operator definition of the NLP soft factor and show that its UV renormalon is still correlated with linear UV divergence rather than the linear rapidity divergence.We show that the power expansion in this example fails to sum to the full result due to extra exponentially small contributions and comment on features of TMD expansion in b ⊥ and k ⊥ space.
4. In Sec. 5, we calculate the bubble chain diagram for the quark field wave function renormalization in Coulomb gauge and show that it suffers linear renormalon ambiguity as well.In Sec. 6, we discuss the practical implications of this linear renormalon to lattice extraction of Collins-Soper kernel and then conclude.
Useful technical details are collected in the Appendices.

Linear renormalon of the heavy-light Sudakov hard kernel
In this section, we calculate the bubble chain diagram for the heavy-light Sudakov hard kernel, or the common hard kernel of double-log type factorizations in the context of lattice parton distributions.
As a reminder, the heavy-light Sudakov hard kernel is a quantity defined purely in PQCD, with UV and IR divergences regulated in DR and subtracted multiplicatively in the MS scheme.It is similar to the standard light-light Sudakov hard kernel, but with one of the external light-parton states replaced by an "on-shell" gauge link in n 2 z = −1 direction.More precisely, one has where |p⟩ is an external collinear quark with the momentum p = (p z , 0, 0, p z ) and ψ(0) is a quark field.In this paper, we use the notation [y, x] to denote a Wilson line from x to y along the straight line between them.If y = ∞, then [±∞v, x] means a Wilson-line from x to infinity along the direction ±v.More precisely, one has In Eq. (2.1), p z = p z > 0 should be identified as quark-parton's momenta xP z in factorization formulas, and we use sign(z) to denote the direction of the n z .Due to Lorentz invariance as well as the freedom to rescale n z → ρn z , the natural scale of the form factor is z , which enters the argument of the form factor through the logarithm L z ≡ ln ζz µ 2 .Moreover, the H HL depends on the direction of n z , σ = sign(z) only through the combination The hard kernel H HL satisfies the standard form of double-log RGE [35,39,40], In the above RGE, besides the universal light-like cusp anomalous dimension Γ cusp , there is a single-log anomalous dimension γH , which can be expressed as [35,39,40] γH = 2γ F + γ V + 2γ HL − γ s . (2.5) In the above, γ F is the UV anomalous dimension for heavy-light current [58,59], γ V is the constant part for the anomalous dimension of the standard light-light Sudakov hard kernel [60], γ s is the soft anomalous dimension, which can be defined through UV anomalous dimension of a light-like Wilson-line cusp [61], and γ HL is the constant part of the UV anomalous dimension of a heavy-light Wilson-line cusp originally calculated in Ref. [62].They are known to NNNLO.
Given the RGE, it is possible to introduce the RGE re-summed form of the hard kernel.For this purpose, it is convenient to introduce ) In terms of the above, the RGE re-summation reads [40] where the various RGE resummation factors reads [60] (2.12) Notice that there are imaginary parts of the heavy-light Sudakov form factor, resulting in the iπsign(z) term in the anomalous dimension and the a Γ term in the RGE resummation.
As demonstrated in Ref. [39], such an imaginary part can be obtained by applying "cuttingrules" 1 k z ±i0 → ∓iπδ(k z ) on the eikonal propagators which make one of them on shell.On the other hand, the eikonal propagators are not pinched at such region.In fact, the ⊥ region in the loop integrals can be deformed away combining contour deformations k 2 .This implies that the heavy-light Sudakov hard kernel as well as the quasi-LFWF amplitudes to be introduced later are deeply inside the natural region of analyticity, in a way similar to the angle-dependent cusp anomalous dimension with a iπ 2 in the hyperbolic angle [41].Nevertheless, the presence of an imaginary part or a phase angle in the hard kernel still results in some inconvenience in practical applications [64,65].It is the purpose of this section to show that, the imaginary part of the heavy-light Sudakov hard kernel, or more precisely, the phase angle A(L z = 0, α( √ ζ z )), contains a linear renormalon ambiguity.Here we calculate the heavy-light bubble chain diagram in Fig. 1.The momentum space gluon propagator with n-bubbles inserted reads Here 3π , and u(p) is an on-shell Dirac spinor with p = (p z , 0, 0, p z ).Now, the diagram for sign(z) = σ = −1 reads (for σ = 1 one simply fips the imaginary part) where µ 2ϵ 0 = µ 2 ϵ e γ E 4π ϵ and n z = (0, 0, 0, 1).After introducing Schwinger parameters and after Wick rotation, one has (2.17) Now, the integrals all factorize and can be evaluated into a product of gamma functions.The result reads (p z = −p • n z ) .18)This implies that the original diagram V n reads This can be re-written as In particular, the bare matrix element has no factorial growth in the n → ∞ limit at all.

Renormalization, RGE and anomalous dimensions
Given the above, the partially renormalized bubble chain diagram reads To proceed, one needs the formulas [52][53][54] n k=0 ) Here, k is the nth harmonic number.They are all special examples of the second Stirling numbers n!S(j + 1, n + 1) and we collect useful formulas for them in Appendix A. Given these, the partially renormalized bubble chain reads Now, using the regular Taylor expansion of V (ϵ, s) at ϵ = 0 where V j (s) are regular at s = 0, the fully renormalized result in MS scheme reads Now, since V (ϵ, 0) and d ds V (ϵ, s)| s=0 are analytic functions in the neighborhood of ϵ = 0, the first two terms are Borel summable.On the other hand, the last term contains renormalon ambiguity due to the pole of V 0 (s) at s = 1 2 .One can see that it is the sub-divergence subtraction that introduces the factorial growth, explaining the name for "renormalon".
Given the full results for the bubble chain, as a by-product, we can check that they indeed satisfy the RGE in the large β 0 limit and extract all the anomalous dimensions.More precisely, we consider the exponentiation of the bubble chain diagrams (in the large β 0 or large n f expansion, exponentiation or not makes no difference up to O( 1 n f )) From the above, the scale µ dependence reads Thus one can obtain the anomalous dimension in the large β 0 limit through calculating dH HL d ln µ .For this purpose one needs the Taylor expansion of V (ϵ, s) at s = 0 where g j (ϵ) are analytic at ϵ = 0.One also needs the scale µ dependencies for the objects Based on the above identities, after certain simplifications we obtain where the identity H n = 1 0 dt 1−t n 1−t is used in the above process.The expressions for g 0 and g reads , . (2.36)This implies that the exponentiation of the bubble chain diagram e ∞ n=0 RVn is closed under RGE.From the above, we can extract the cusp anomalous dimension for the bubble chain diagram [54] Γ whose small α expansion (one also replaces The result is also consistent with that from [54,60] in large n f limit.The imaginary part of the single log anomalous dimension is which is consistent with imaginary part of the known result.On the other hand, the single-log anomalous dimension γH for the bubble chain diagrams reads The small α expansion of γH (one also replaces which reproduces the single log anomalous dimension in the large n f limit up to NNNLO.
In addition, we can also check the constant terms in H HL up to NNLO (setting which are consistent with the n f part for constant terms c a and c H [40].

Borel transform and linear renormalon
We now move to the key part of this section.Given the fully renormalized hard kernel in Eq. (2.29), the Borel transform reads (u = β 0 t 2 ) where R(u) is an entire function of u.Now, one has The Borel transform of the imaginary part reads now This implies that the Borel transform has a pole at u = 1 2 , with the residue Clearly, the linear renormalon appears only in the imaginary part.As such, the asymptotic behavior of the phase angle A(L z , α) can be predicted as where the factor 2 here comes from the factor 2 in the definition Eq. (2.7).The coefficient is In the MS scheme, through comparing Eqs.(2.47) (2.49), one has It is interesting to notice that the above amplitude of the renormalon ambiguity is exactly the same as that for the linear divergence or the heavy-quark pole mass in the large β 0 limit [66][67][68].Eq. (2.47),(2.49)and (2.50) are the main results of this section.
Through requiring that the renormalon ambiguity is renormalization scale independent, one can generalize Eq. (2.49) beyond large β 0 limit [66,67], ) Here we propose a conjecture that the amplitude N m for the phase angle A is the same as that for the linear divergence or the heavy-quark pole mass even beyond large β 0 limit, which means where those values are obtained based on the method in Ref. [67] using the perturbation series between the MS and the pole quark masses [69].Here the N m = 0.622 for n f = 0 is consistent with that obtained from the linear divergence series N m = 0.660(56) [70].There are several evidences for this conjecture.First, they share the same N m in the large β 0 limit.Second, as we will discuss later, the renormalon series for the heavy light Sudakov hard kernel matches that of NLP TMDWF, whose Borel transformation is very similar to that of the static potential, see Eq. (4.26).As shown in [71], the renormalon series of the static potential cancels that of the self energy in a Wilson loop.

Renormalon cancellation and power corrections to threshold asymptotics of space-like quark-bilinear coefficient functions
In this section, we investigate the role played by the linear renormalon of H HL in an important place where it appears, the threshold limit of the hard function or coefficient function that relates the space-like quark-bilinear matrix elements ("qPDF" or "pPDF") to twisttwo light-front correlators.Throughout this section, |z|Λ QCD ≪ 1 is always understood.The main object is the perturbative quark correlator in coordinate space where |p⟩ is the on-shell external quark state with p z > 0 and λ = zp z .All the IR divergences are subtracted in the M S scheme through convolution and the tree level is normalized to e −iλ .In this paper we only consider quark non-singlet un-polarized contributions with ⃗ p ⊥ ≡ 0. Clearly, this correlator is the leading-twist OPE coefficient function for the space-like quark-bilinear transformed into coordinate space.Our convention for the moment space read As usual, the support in α is in [−1, 1] (this implies that the coordinate space version is analytic in the whole complex λ plane) and for the purpose of extracting the α → 1 threshold asymptotics through the N → ∞ limit, in principle one should integrate from 0 to 1.However, due to lack of contributions at O(1 + α) −1 as α → −1 for quark non-singlet un-polarized coefficient functions, to leading power in 1 N , the [−1, 0] part can be added as well.Also notice that at the level of bubble chain diagrams, supports are always in [0, 1].
As shown in [40] to NNLO, the coordinate space coefficient function in the threshold limit λ → ∞3 , has the following factorization formula at leading threshold power in Here H HL is nothing but the HL Sudakov hard kernel, and the LP space-like threshold soft factor J or "jet function" is a specific Wilson-loop average involving space-like and light-like gauge-links defined in [40] J where l z ≡ ln e γ E |z|µ 2 is the natural logarithm in z.Here we emphasize that in the threshold λ → ∞ limit, the threshold logarithms are due to the hard scale |λ| |z| ≫ 1 |z| ≫ Λ QCD and are exactly resummed using the heavy-light Sudakov hard kernel.Notice that the threshold power is finer the standard "twist-power" counted in zΛ QCD , it is counted in 1 λ or 1 N .A single twist in the threshold limit can split into infinitely many threshold powers and there is no conflict between these two expansions due to the fact that the threshold limit is UV in our case.
Although the initial work is in coordinate space, we can convert all the results to the moment space.More precisely, to leading threshold power in 1  N , one has for N -th moment for the coefficient function which can be converted into the less convenient "A-B-C" form.Again, the threshold scale

2N
|z| in moment space is hard.To obtain the fully resumed form, one needs the RGEs In the above, γJ is the single-log anomalous dimension for the space-like threshold soft factor and equals to 2γ HL − 2γ s .The RGEs of the hard kernel and the threshold soft factor allow the systematic resummation of threshold-logarithms in a way similar to Eq. (5.5) of [54] without using the "A-B-C" form.Moreover, after partial integration, all the α(µ) dependencies can be factorized out completely and the controlled asymptotic expansion in the threshold limit in moment space takes a simple factorized form in terms of α(z as we show now.The controlled expansion in the coordinate space is similar.
To write out the controlled expansion in the N → ∞ limit for the H N , it is convenient to define the following regular functions (β(α Γ α 2 + .... , (3.9) Γ α 2 + ... and The definition of f γ is similar for any anomalous dimension γ.They are all renormalon-free provided that the corresponding β function and anomalous dimensions are renormalon-free.
In terms of the above, one has the fully-resummed form in the Mellin space to leading power in 1 N : (3.12) In the above, the µ independent re-summation factors ÔH and ÔJ are naturally written in terms of α(z Notice that the α(z N ) factorizes completely from α(z), explaining the name of threshold factorization.Moreover, in the N → ∞ limit, the re-summation factors ÔH is not only a controlled asymptotic expansion, but also free from any renormalon.Nevertheless, the initial condition H 2 HL (0, α(z N )) of the heavy-light hard kernel, although still a controlled asymptotic expansion, has linear renormalon ambiguity of the order This reflects the intrinsic order of ambiguity for the separation between the leading and the next-to-leading threshold powers or equivalently, the separation between two hard scales |λ| |z| and 1 |z| .Notice that unlike the "standard" case of the DY or DIS threshold limits, the threshold expansion here do not compete with the standard twist expansion due to the fact that N is in the denominator instead of numerator.Indeed, higher-order renormalons in the HL hard kernel correspond to , k ≥ 2 ambiguities and are always controlled in the large N limit, see Eq. (3.23).Clearly, this is due to the fact that the threshold scale z N for the space-like parton distribution is hard and the N → ∞ limit is a UV limit.To demonstrate the above statements concerning the structure of renormalon, in this section we investigate the bubble chain diagram for the coefficient function H in the threshold limit.The purpose is to show 1.At the level of bubble chain diagram, the threshold factorization formula is correct to leading power in 1 λ or 1 N (but to all logarithmic orders).More precisely, we perform the threshold expansion of the bubble chain diagram to leading power both in coordinate and moment space.We show that the resulting hard kernel exactly reproduces that calculated in Sec. 2, and the resulting "jet function" exactly reproduces the bubble chain diagrams for the space-like threshold soft factor.
2. Given the above, the threshold factorization formula is established at the level of bubble chain diagram.We further investigate the role played by the linear renormalon in the heavy-light hard kernel H HL .We show that the linear renormalon cancels between the leading and next-to-leading threshold powers by expanding the full bubble chain results to NLP.We further identify the operator definition of the obtained NLP contribution (at O( 1 |λ| )) and show that the s = 1 2 renormalon for the NLP soft factor is UV in nature and is due to implicit ambiguity of subtracting UV divergences.Notice this has nothing to do with the linear UV divergence for space-like gauge-link self-interactions.
To summarize, the factorization formula Eq. (3.12) can be cast into a controlled asymptotic expansion in the threshold limit using the RGEs.Nevertheless, the separation between leading and next-to-leading threshold power contributions (more precisely, separation between |λ| |z| and 1 |z| ) requires a scheme, reflected in the (IR) renormalon of the leading threshold power hard contribution of order and (UV) renormalon of the next-to-leading threshold power "soft" contribution of the same order ( 1 N × Λ QCD |z|).The scheme dependency cancels between LP and NLP, or equivalently, the renormalon cancels between LP and NLP.
It is also helpful to make the following comparisons to the DY threshold limit investigated in [54]: 1.In the DIS or DY case 4 , in the z → 1 limit the threshold scales Therefore, the resummation factors in the "A-B-C" form are sensitive to the Landau pole and finite-order truncation of the anomalous dimensions can create renormalons.In our case, in the λ → ∞ or N → limit the threshold scales λ |z| or N |z| are UV and the resummation factors are always far away from Landau pole.Moreover, we use the RGEs for the hard kernel and the soft factor to perform the resummation in a way similar to Eq. (5.5) of [54], instead of truncating the "A-B-C".As such, the linear renormalon observed here has absolutely no relation to resummation.
2. Due to the above, in [54] the authors focus on soft and collinear contributions in the threshold limit.They show that the "soft-collinear" contribution which dominates the DLA cancels with "wide-angle" soft emissions in terms of linear renormalon.
In particular, they show that the bubble chain diagram for the DY threshold soft factor is free from linear renormalon and this is consistent with the power correction extracted using finite gluon mass.In our case, the similar object that captures "IR" contributions is the space-like threshold soft factor J(l z , α(µ)) (similar to the DY case, it is also a Wilson-loop) at scale 1 |z| and it is free from linear renormalon as well (except the overall UV renormalon for the space-like gauge-link inherited from the definition of pPDF and have nothing to do with threshold expansion).But this is not the object that is responsible for the threshold logarithms in our case.
3. In our case, the object that is responsible for the threshold logarithms is in fact the heavy-light hard kernel.It is purely virtual instead of real, and it is not purely eikonal.Unlike the standard light-light Sudakov hard kernel (namely, the H(α(Q)) in Eq. (1.4) of [54] ), in our hard kernel, the heavy gauge-link and the light parton are not symmetric, as a result, when a finite gluon mass is added, it develops linear power correction, consistent with the observed linear renormalon.This also means that the cancellation of the linear renormalon must involve sub-eikonal effects, therefore must be between leading and next to leading threshold powers, consistent with our observation.
Here we also comment in terms of the relation to the original calculation in Ref. [1].In fact, the only question we need to answer is that: why in the original work there is almost no linear renormalon (except for the trivial one localized into the gauge-link self energy diagrams), but in this work there is.The reason is again due to the projection to the leading threshold power.Before the threshold expansion, both the virtual and the real diagrams are of the same twist in |z|Λ QCD and the linear renormalon in the virtual diagram is completely cancelled by the real diagram.After projecting to the leading threshold power, the virtual diagrams survive, but the real diagrams eikonalize and become part of the space-like threshold soft factor.The contributions in the real diagrams that cancel the linear renormalon are sub-eikonal and move to the next-to-leading threshold power 1  N .The linear renormalon then simply measures the intrinsic ambiguity of the threshold-power separation (but not z, Λ QCD separation).
Let's emphasize here that the order of our linear renormalon ambiguity is and is much smaller in the N → ∞ limit than the "genuine" high-twist contribution O(|z| 2 Λ 2 QCD ).In particular, the threshold expansion is still performed within a given twist, including the order of ambiguities.Also notice that although we call the linear renormalon in the HL hard kernel "IR", it is only IR with respect to the scale |λ| |z| for the hard kernel, but not to 1 |z| .In fact, we will see that the linear "IR" renormalon cancels with UV renormalon of NLP "soft" contributions at scale 1 |z| .
Due to the above, the linear renormalon should not be regarded as a prediction of presence of genuine non-perturbative effects enhanced in the threshold limit.Instead, it should be interpreted as predicting the existence of next-to-leading threshold power contributions to the leading twist coefficient function that has linear UV renormalon.Moreover, all the "threshold renormalons" generated after threshold expansion always cancel between different threshold powers, and the "twist-expansion renormalons" separating the leading and high powers in the zΛ QCD expansion remain the same as given by Eq. ( 48) and Eq. ( 49) of the original work [1].In the threshold region, they start contributing only from at 1 N and (zΛ QCD ) 4 at N 0 .Finally, both the "threshold renormalons" and the "genuine" renormalons in [1] never increases as N → ∞.Genuine non-perturbative effects are uniformly bounded by (zΛ QCD ) 2 itself in an N -independent way.

Collection of full results before threshold limit
In this section we consider the bubble chain diagrams for Eq.(3.1) between the incoming quark and the space-like gauge-link corresponding to Fig. 1a in [1] (another diagram Fig. 1b that couples the outgoing quark and the gauge link contributes equally).Only these two diagrams contribute to non-trivial threshold logarithms in the threshold limit.The gaugelink self-interaction diagram Fig. 1d depends only on z 2 and is part of the LP threshold soft factor.The one gluon exchange diagrams Fig. 1c between quark and anti-quark start contributing only at 1  N 2 and contain no linear renormalons either before or after threshold expansion so will not be considered.Notice that all results are normalized in a way such that the tree-level diagram reads e −iλ and only un-polarized component is considered.In the notation of [1], our results always corresponds to H ∥ , although for the diagrams we considered there is no difference between H ∥ and H ⊥ .
Since the calculation for all the diagrams has already been performed in [1], we will not present calculation details.Instead, we start from the following expression for unrenormalized n-bubble chain contribution (corresponding to Fig. 1a in [1]) one has After integrating over x and t, the above leads to ( 2 F2 is the regularized hypergeometric function and s = (n + 1)ϵ) (2, 2s + 1; −ϵ + s + 3, 2s + 2; −λ) . (3.17) Notice that the above is equivalent to Eq. (A6) in [1] and we will come back to this in Sec.3.5.Converting into the moment space, one has for the moments (the notation for h n,N (z 2 ) follows the same convention as h n (z 2 , λ)) Here we adopted the notation The above serve as the starting point of the following discussions.Notice that the above has no pole at s = 1 2 , ϵ = 0 , thus after renormalization in standard manner, no linear renormalon can be generated.This is consistent with the original result in [1].

Threshold limit in coordinate and moment space
Now we study the threshold limit of the full results above.We perform large λ and large N expansions up to NLP.We show that the threshold logs at leading power are due to hard scale λ 2 z 2 and the corresponding hard contributions factorize exactly into the heavy-light hard kernel calculated in Sec. 2. The soft contribution at LP corresponds exactly to the bubble chain diagram between n + and n z for the space-like threshold soft factor Eq. (3.5).This verifies the threshold factorization formulas Eq. (3.4) and Eq.(3.6) at the level of bubble chain diagram (previously it has been verified to NNLO).The pole at s = 1  2 cancel between LP hard and NLP soft contributions.
We first consider the threshold limit in the coordinate space.Given Eq. (3.17), its asymptotic expansion in large λ limit reads It is clear now that the first term, after renormalization, simply generates all the threshold logarithms at leading power.From the combination |z| 2λ 2s it is clear that threshold logs are due to the hard scale 4λ 2 z 2 and is UV in nature.Moreover, this term is nothing but the n-bubble diagram for the heavy-light hard kernel given in Eq. (2.18).On the other hand, it is not hard to show that the first term in the second line simply reproduces the n-bubble diagram between n + and n z for the LP space-like threshold soft factor in Eq. (3.5) and is free from singularity at s = 1 2 .Therefore, the leading power threshold factorization formula Eq. (3.4) is verified to all orders at the level of bubble chain diagrams.Moreover, the s = 1 2 pole for the HL hard kernel is simply cancelled by the NLP contribution We will show that the s = 1 2 pole for the NLP contribution is in fact UV in nature, in subsection.3.4.
Similarly, starting from Eq. (3.18) one can perform the threshold limit N → ∞ in the moment space as well.The result reads The first line again corresponds to the heavy-light hard kernel, while the second line corresponds to soft factors at LP and NLP.Therefore, the threshold factorization formula Eq. (3.6) is established at LP in the moment space as well.Notice that the term − s(2s+1) N in the first line is due to "kinematic" 1 N contributions and is absent in the coordinate space version Eq. (3.20).They can be generated naturally when converting the coordinate space version to moment space.The full contribution to the moment from the hard kernel simply reads which reproduces the first line in Eq. (3.22) after expanding in N .Notice the above is similar to Eq. (5.19) of [54], in fact in our notation Ref.
[54] has 1 The difference in sign of the 2s in the denominator reflects the deep fact that the threshold scale in [54] 5 .Also notice that to Fourier transform the half-plane lim In this way one maintains the correct support property and one can show that the result Eq. (3.20) after converting to the moment space (regarding the issue of support in α space, see Sec.To summarize, the threshold factorization formulas Eq. (3.4) and Eq.(3.6) are verified exactly at the level of un-renormalized n-bubble chain level.After renormalization one obtains the corresponding renormalized version.Notice that the 1 s 2 poles at LP simply cancel between the hard and soft contributions, and the total anomalous dimension remains single-log in nature and simply compensates the N → ∞ limit of the γ qq N up to 2γ F 1 0 All these anomalous dimensions are free from renormalons at the level of bubble chain diagrams.

Borel transform and conservation of "genuine" IR renormalons in [1]
To check consistency, one can also perform the threshold expansion at the level of the Borel-transforms of the renormalized results and verify the threshold factorization.Here we only present the results in coordinate space.Up to terms which are entire functions, by performing the threshold expansion on the M S scheme Borel transform, in terms of the natural Borel argument u = β 0 t 2 , the first two orders in the large λ expansion reads In this equation, the first term is nothing but the Borel transform of the hard kernel in Eq. (2.44), (2.45).The second term is the Borel transform of the (bubble chain diagrams between n + and n z ) space-like threshold soft factor.The u = 1 2 singularity only appears in the hard contribution.The first two terms represent the leading power contribution in the threshold expansion.On the other hand, the term is the NLP soft contribution.They are clearly in one-to-one correspondence with the three terms in Eq. (3.20) and can also be directly obtained from Eq. (3.20) using standard methods [52][53][54][55].Again, the linear renormalon at u = 1 2 cancels between LP hard kernel and NLP soft factor.
It is important to notice that although the u = 1 2 renormalon for the hard kernel generated during the threshold expansion cancel between leading and next to leading threshold powers, the renormalons for the LP threshold soft factor is not cancelled by anything and exactly corresponds to the "genuine" renormalons in [1] projected to the leading threshold power.They can only be cancelled by UV renormalons for high-twist operators (true "non-perturbative" effect) in the zΛ QCD expansion.To show this we need to add to twice of Eq. (3.28) the contribution from the gauge-link self interactions (namely, Fig. 1d of [1]) Notice the striking similarity of the Borel transform for the self-energy contribution and the NLP soft factor.Combining the above, one has the Borel transform for all the three bubble chain diagrams for the LP threshold soft factor Eq. (3.5) (3.31) Clearly, the above contains IR renormalons at u = k for k ≥ 2, with residues read which exactly corresponds to the δ(1 − α) term in Eq. ( 49) of [1] (the factor 4π difference is due to the fact that in [1] the Borel transform is performed with respect to αs 4π ).To conclude, the LP threshold soft factor Eq. (3.5) correctly captures all the genuine IR renormalons for the pPDF OPE coefficients projected to the leading threshold power.There is essentially no enhancement of "true" non-perturbative effect due to scales softer than 1 |z| in the threshold limit.

The NLP threshold soft factor and its UV renormalon
In this subsection we address the nature of the NLP contribution Eq. (3.21).We show that it can be generated from a NLP threshold soft factor at the level of bubble chain diagrams.Moreover, we show that its s = 1 2 pole is UV.In fact, the NLP soft factor one needs simply reads (the NLP contribution for Fig. 1b of [1] is identical.Fig. 1c starts contributing only at NNLP.) Here, n + = 1 √ 2 (1, 0, 0, 1) is the light-front "plus" direction vector, and [y + , x + ] denotes a light-like Wilson line in the n + direction from x + n + to y + n + .This can be explained as the sub-eikonal insertion at the collinear quark line (explaining the subscript q in the definition) One can show that the bubble chain for e −iλ p + J q (z, α(µ)) reproduces the NLP threshold power contribution Eq. (3.21).More precisely, for z > 0 one has with Now, after integration one has which is exactly the NLP contribution in Eq. (3.21) obtained from direct expansion.For z < 0 one can see it reproduces Eq. (3.21) as well.
Here we show that the linear s = 1 2 pole for the J q above is UV in nature.For this purpose one can change the ∞ 0 λ 2 dλ 2 into ∞ a λ 2 dλ 2 with a > 0 playing the role of a UV cutoff.This only changes the short distance singularities of the integral.Then, the full result reads As expected, the singularity at s = 1 2 simply vanishes.Now, after power-expansion in a, one has Clearly, the second term simply corresponds to the finite result evaluated at a = 0 and the first term corresponds to the "power UV-divergence" due to the fact that the operator dimension for J q equals to 1.The pole at s = 1 2 cancel between the "UV divergent" term in the first line and the "finite term" in the second line, implying that the s = 1  2 renormalon of the NLP soft factor is actually a UV renormalon.More precisely, this corresponds to the ambiguity of separating leading (UV divergence) and next to leading power (finite renormalized term) in az.Thus, the renormalon cancellation pattern in the coordinate space threshold expansion is fully in agreement with the general pattern of renormalon cancellation between LP and NLP for controlled asymptotic expansions in UV limit.

Exponentially small terms and resurgence of the full distribution from threshold asymptotics
In fact, the full asymptotic expansion for Eq.(3.17) in the threshold limit reads .40)where in the last line we included all the exponentially small contributions related to α → 0 + asymptotics.Clearly, the threshold expansion in coordinate space projected only to algebraic terms is not convergent, but we will show that the "Borel ambiguity" for the 1 λ expansion of these terms cancels the branch ambiguities of (−λ) ϵ−s for the exponentially small term.Notice the "threshold scale" for the exponentially small term is Here we should stress that the overall coefficient (−λ) ϵ−s in the last line is in fact at the branch cut for λ > 0 and requires a scheme to be regarded as an analytic function in the Re(λ) > 0 plane.This is due to the fact that only for Re(λ) < 0, this term is the leading large λ asymptotics and in the Re(λ) > 0 case it can only be determined when the leading terms (threshold terms) are determined.Indeed, the infinite sum in the second line is not Borel summable and requires a scheme to be well-defined.The scheme choice for the Borel sum then transmutes to the branch choice for the (−λ) ϵ−s term to make sure the total result is analytic in λ in the whole complex plane, or equivalently, to make sure the α < 0 part is completely cancelled between the threshold and small α contributions.
On the other hand, since in this case the Fourier transforms for the exponentially small contributions are only supported in α < 0, one can simply take the threshold contributions and truncate the Fourier transforms to 0 < α < 1 to obtain the full distribution in the α space.Notice in this region the Fourier transform for the infinite sum in the second line is free from ambiguities since Fourier transform introduces additional 1  (k−1)!and the convergence radius for the Fourier transform is equal to the convergence radius for the Borel sum for the 1 λ expansion.In this way only the threshold asymptotics in λ space is needed to recover the full distribution due to the correlation between the algebraic and exponentially small terms, demonstrating the principle of "resurgence".
To check the above statements, it is convenient to Fourier transform back to the α space.Then the threshold terms read (they are supported in α ≤ 1) Now, the contribution for the exponentially small contributions reads (it is supported only in α < 0) To show that it compensates the threshold contributions for α < 0, it is convenient to use the connecting formula for hypergeometric functions that changes 1 − α to α where a = 1−2s, b = ϵ−s, c = 2−2s.Now, for α < 0 the original 2 F 1 at the left-hand side of the equation has a branch cut, reflected in the α c−a−b term.One can see from the above that if the branch for this α c−a−b is chosen in agreement with the (−1) s−ϵ in Eq. (3.42), then the α < 0 part exactly cancels between Eq. (3.41) and Eq.(3.42).As a result, the full α space distribution is given by Eq. (3.41) in 0 < α ≤ 1.Notice that in the region 0 < α < 1, after applying the connecting formula, Eq. (3.41) is directly equal to Eq. (A6) of [1].To show this, first converting the 2 F 1 to the "plus-distribution" Then combing all δ(ᾱ) terms one has in coordinate space Now, using the connecting formula for the hypergeometric function, the ᾱ−1+2s term is cancelled, left finally with which is identical to Eq. (A6) of [1].In fact, the above procedure, when reverted, can also be used to derive the threshold limit.Nevertheless, the fact that we started from Eq. (3.17) and recovered the full distribution from threshold asymptotics provides a strong consistency check of the whole procedure as well.
3.6 A conjecture on threshold limit of hadronic quark-bilinear matrix element Notice that although up to now, we are discussing the threshold-limit of leading twist coefficient function for quark-bilinears, the results allow extension to threshold limit of the hadronic matrix element [1] ⟨N For simplicity we only consider the flavor non-singlet valence-quark induced hadronic matrix element.Clearly, for z > 0 the above matrix element allows analytic continuation in v with imaginary part in the negative forward light-cone, in particular, with v = −i(1, 0, 0, 0) in the Euclidean time direction.Compared with v = n z , this does not change −v 2 = 1, but changes λ → iλ, which is exactly the positive imaginary axis along which threshold limit can be extracted.Furthermore, using one can see that in the large λ → +∞ limit, at leading twist, one has for the where Q PDF (iλ, α(µ)) is the threshold limit of the quark PDF (without the e λ factor) and z λ = z λ = 1 p 0 is the hard scale.Resummation factors are defined in Eq. (3.12), Eq. (3.13).In fact, due to the fact that the µ dependency cancels between the PDF and the α(µ) part of the resummation factor, one can combine them to form the µ-independent version of quark PDF, in term of which the leading twist threshold limit takes the form It is clear now that all the remaining µ dependency is due to the UV anomalous dimension of the quark-bilinear itself.
Now, one may question if the leading twist always dominate the large λ limit.In our opinion, we believe that the leading twist dominance holds to λ → +∞ in the sense of Eq. (3.50) to be introduced.The reason is that, in the large λ limit, due to the exponential suppression caused by the Euclidean time evolution, the large energy λ z of the incoming and out going hadron must flow out and in at the quark-line vertices, leaving the intermediate state (namely, the "cut") light.In this way, the large λ limit is dominated by threshold contributions, for which the two vertices are always in Sudakov-like hard kernels with at least one incoming collinear parton and one external soft gauge-link.Clearly, increasing the number of incoming collinear partons (except longitudinal gluons which factorize to eikonal lines in the PDF) and external soft partons to the hard kernel always introduce additional power suppressions.The Sudakov hard kernel with one incoming collinear valence quark and one external soft gauge-link is the hard kernel at leading power, on top of which collinear contributions factorize into the threshold limit of twist-two PDF (this explains the name of "twist-two dominance") and soft exchanges factorize into the threshold soft factor.Moreover, this picture for the dominant contribution in the large λ limit should hold for moderate z as well, due to the fact that the threshold soft factor, defined as a Wilson-loop average, allows non-perturbative generalization.
As such, we conjecture the following λ → +∞ factorization for the hadronic quarkbilinear matrix element from small to moderate z: Here J(zΛ QCD ) is the appropriate non-perturbative version of the threshold-soft factor Eq. (3.5) with µ dependency remove by evolution factors.The small z expansion of the threshold soft factor is exactly Clearly, the fully factorized form of the threshold limit provides considerable simplification compared to the leading twist factorization formula at generic λ.Provided that Q PDF (iλ) does not decay very fast in the λ → +∞ limit, in principle it allows the lattice verification of the celebrated "marginally-relevant" Sudakov-like asymptotics Eq. (3.13) which dates back to late 1970s [48,73].
4 Renormalon cancellation and power correction to TMD limit of quasi-LFWF amplitude In this section, we consider the role played by the linear renormalon of the hard kernel in the TMD expansion of quasi-LFWF amplitude.We show that this linear renormalon cancels with UV renormalon of NLP soft contribution.We also show that in this case the TMD expansion fails to converge to the full result due to exponentially small contributions.As a reminder, the quasi-LFWF is defined as the Fourier transform of the following correlator [32,39] where λ = zP z is the natural scaleless Lorentz invariant length conjugating to the quark's momentum fraction x.The staple-shaped gauge-link W z is defined as On the other hand, Z E is the square root of the vacuum expectation value of a flat rectangular Euclidean Wilson-loop that subtracts out the gauge-link self-interactions in the large L limit.The Z E itself has nothing to do with the factorization.As shown in Ref. [39], the quasi-LFWF amplitude has the following leading-power factorization in the double-logarithmic evolution limit P 2 z → ∞ with fixed x and b ⊥ : where the hard kernel H ± is nothing but a product of two heavy-light Sudakov form factors in Eq. (2.1) The natural scales ζ z = 4x 2 P 2 z and ζz = 4(1 − x) 2 P 2 z are provided by the momenta p = xP of the collinear quark and p = (1 − x)P of anti-quark.
In this section, we only consider the quasi-LFWF amplitude in PQCD for the external states with a pair of on-shell quark anti-quark (with momenta p 2 = p2 = 0), which serves as the small b OPE coefficients that match the mesonic quasi-LFWF to collinear distribution amplitudes (DA).We calculate the bubble chain diagrams for the quasi-LFWF amplitude that couples the external quark and anti-quark to the space-like gauge links, since at bubble chain level only these two diagrams know about the presence of multiples scales at LP.To avoid complications caused by the λ (or momentum fractions), we simply calculate the λ = 0 version of the quasi-LFWF amplitude.It is not hard to see that at the bubble-chain level, the factorization for λ = 0 simply reduces to additions: Here λ = 0 is omitted for simplicity.Below we calculate the bubble chain diagram for the λ = 0 quasi-LFWF Ψ+ (the minus version simply differs in sign for the imaginary part) and demonstrate the renormalon cancellation between leading and next to leading power.

Bubble chain diagram for quasi-LFWF amplitude
Here we consider the bubble chain diagram for the quasi-LFWF amplitude at λ = 0.As explained before, we only consider the two bubble chain diagrams which couple the external quark/anti-quark line with the space-like gauge link.For quasi-LFWF amplitude, two diagrams contribute equally with replacements ζ z → ζz , so we only present the calculation details for one of them, the one which connects the collinear quark to the Wilson line.The diagram contains two scales, ζ z = 4p 2 z and b ⊥ .Using the same notation as Eq.(2.15), introducing one has the integral representation Now we consider the imaginary part (as one can check, only the imaginary part contributes to NLP at 1 ), for which the λ integral can be evaluated exactly, leading to where one have used s = (n+1)ϵ.Now, the integral for ρ can be further performed, leading to It is easy to see that before performing the power-expansion, the above function is analytic at s = 1 2 and ϵ = 0, therefore not contributing to linear renormalon.It still contributes to quadratic renormalon at s = 1.This implies that in the deeply UV region , the OPE coefficient that matches the quasi-LFWF to the collinear DA has no linear renormalon.The linear renormalon is purely a result of the ζ z ≫ 1 b ⊥ expansion in the double-logarithmic limit.

Power expansion and renormalon cancellation
One expects that the linear renormalon in the HL Sudakov cancels with linear renormalon in "high-twist" TMDWF amplitudes or soft factors.Indeed, if one perform the b ⊥ ≪ 1 Λ QCD expansion, then the renormalon ambiguity of twist-three distribution will be of the structure 1 b ⊥ pz × b ⊥ Λ QCD , where 1 b ⊥ pz is the natural power for all the perturbative twist-three contributions, and b ⊥ Λ QCD is the linear renomalon in the twsit-three TMD distributions.This cancellation can actually be demonstrated.Indeed, for this purpose, it is convenient to perform Mellin transform on q = ζ z b 2 ⊥ in Eq. (4.9).Using the Mellin transform reads with the strip of convergence and the condition ϵ > − 1 n .It is easy to see now, that the heavy-light Sudakov contribution at leading power Eq.(2.18) exactly corresponds to the pole at t = 0, while the IR contribution at leading power corresponds to t = −s : Imψ The NL power correction corresponds to the pole at t = −s + 1 2 , with residue Given this, contributions at various powers can be obtained now, by shifting the vertical line for Mellin inverse transform to the right of poles and pick up the residue at the pole.More precisely, the leading power contribution reads The first line corresponds to the hard contribution while the second line corresponds to the soft contribution.The NLP contribution reads It has the desired linear renormalon at s = 1 2 and one can show that it precisely cancels the linear renormalon in the leading power hard contribution.Indeed, at leading power one has for the quasi-LFWF At NLP one has Furthermore, one can check that the real part has no contribution at power 1 pzb ⊥ .This implies that the Borel transform of the LFWF at LP reads (u = β 0 t 2 , L b = ln with residue at u = 1 2 reads This is the same as Eq.(2.47) for z < 0. On the other hand, the Borel transform at NLP reads with the residue at u = It cancels the linear renormalon at leading power.Moreover, the Borel transform is almost identical to the bubble chain diagram for the heavy-quark potential, except for the 2u + 1 term.
Note that the above analytical results may improve the accuracy in extracting the CS kernel, especially in suppressing the imaginary part.See Appendix C for more discussions.

NLP soft factor. Linear UV vs rapidity divergences
Now the question arises, can one attribute the above NLP contribution to twist-three TMD distributions or soft factors?In fact, the NLP contribution at the bubble chain level 6 can be explained by the soft factor which is similar to Eq. (3.33) for pPDF in the threshold limit.Given the above, it is not difficult to show that the bubble-chain diagram for 1 p + S q exactly reproduces the NLP contribution Eq. (4.26).Another bubble chain diagram between the antiquark and spacelike Wilson-line has a similar S q soft factor which contributes the same as S q in case of quasi-LFWF but cancels with S q for quasi-TMDPDF.The only scale for S q is b ⊥ , explaining its name as a "soft factor".More precisely, one has with In this above, ν is the rapidity regulator which should be sent to 0 later, η is a mass-scale similar to the δ + in the delta regulator.Now, after Wick rotation and performing the integrals, one has From this, it is clear that the ν → 0 limit is regular and will not generate any logarithmic rapidity divergence.Sending ν to zero, one has which simply reproduces the results Eq. (4.19) before.A special feature of the above NLP soft factor is that in its definition one needs a "high-dimensional Wilson-line cusp" with dimension 1, therefore in principle can contribute to power-like UV divergences.Moreover, the correlator Eq. (4.28) suffers from linear rapidity divergence which requires special treatment.Although in regulators without introducing new scales into Feynman integrals such as the analytic regulator [75] adopted above, the linear rapidity divergence in Eq. (4.28) is simply removed when removing the regulator ν → 0 like the linear divergence in DR when ϵ → 0, the existence of such divergence strongly suggests that the object in definition still suffers rapidity regularization scheme dependency.So which one among these two linear divergences is related to the linear UV renormalon?The answer is, it is not the linear rapidity divergence, but the linear UV divergence that correlates with the linear renormalon.Indeed, one can introduce a displacement a in z direction in a way similar to the exponential regulator [76,77] to regulate the rapidity divergence as well.It is not hard to show that the bubble chain diagram after introducing this regulator still contains s = 1 2 pole at a 0 that can not be canceled by "linear rapidity divergences" at a −1 .Only after introducing the cutoff regulator for the "high-dimensional Wilson-line cusp" in a way similar to that in Sec.3.4 the s = 1 2 pole can be removed, demonstrating the UV nature of the renormalon.Thus, the linear rapidity divergence is not sufficient to generate linear renormalons, at least in the current case.

Exponentially small contributions and divergence of TMD expansion
In this example, one can make the interesting observation: for any s, the power expansion based on Mellin's transform is not convergent.In fact, the power-expansion corresponds to poles at with n ≥ 0. The non-convergence of the OPE can be read from the negative residue which diverges factorialy in the n → ∞ limit.The above implies that the quasi-LFWF amplitude at power (p z b ⊥ ) −n , n ≥ 1 has the Borel transform which implies that the power expansion diverges for a generic u > 0. For example, the derivative of the above at u = 0 reads for n ≥ 2 which growth factorially.The Borel transform of the above reads which has a singularly at t = b ⊥ p z , suggesting naruall correlation with exponential small contributions at e −b ⊥ p z .
In fact, the OPE non-convergence as well as existence of such an exponential small contribution can be shown even at one-loop level for s = ϵ.In this case, the imaginary part of the LFWF is UV finite and has the representation Now, after a few simplifications, the integral reads Clearly, the first line corresponds to the standard power expansion at leading power and NLP, while the second line is exponential-small contributions that can not be cast into the form of a power expansion.They are caused by "tunneling" of hard momenta between the two spatially separated hard kernels at 0 and ⃗ b ⊥ .On the other hand, when Fourier transformed into k ⊥ space, such exponentially small contributions will result in algebraically power-suppressed terms starting from power 1 p 2 z .This implies that the power expansion in b ⊥ space and k ⊥ space is not likely to be commuting starting at quadratic power.Moreover, compared to the case of 1D threshold-limit in Sec. 3, reconstruction of exponentially small terms from algebraic asymptotics for TMD expansion, due to the 2D nature, is less trivial even in the presence of support property (such as the Wightman-like current-current correlator for DY process, for which one has |q ⊥ | ≤ q 0 due to spectral condition).As such, the b ⊥ space TMD expansion, due to the clean separation between "universal" algebraic terms and "non-universal" exponentially small terms, has the advantage in terms of distilling the massless scaling limit.On the other hand, the k ⊥ space TMD expansion condenses more information into power-log terms at small k ⊥ and has better phenomenological performance.

Linear renormalon of quark wave function renormalization in Coulomb
gauge.
Here we study the hard kernel for the recently proposed Coulomb-gauge approach to the lattice-TMD [78,79].In this case, the gauge link is replaced by the following gauge transformation (called "dressing factor" in [79] and references therein) which transforms a generic gauge to Coulomb gauge Notice that in QED the O(g 2 ) corrections vanishes.The Coulomb gauge hard kernel is defined as the following object where all the divergences are regulated in DR and subtracted multiplicatively in the MS scheme.As for the heavy-light Sudakov hard kernel, the above Coulomb gauge hard kernel can be interpreted as a gauge-invariant definition of the Coulomb gauge quark field renormalization.
We now calculate the bubble chain diagram of the above hard kernel.The diagram with n bubble-chain insertions reads where one has Now, for p = (p 0 , 0, 0, p 0 ), using rotational symmetry and equation of motion, one can write Introducing the parameters for | ⃗ k| 2 , after Wick rotation one has the following representation Now, in terms of s = (n + 1)ϵ, the result reads . (5.7) One can see that the above simply reproduces the hard kernel obtained in the recent paper [79] for s = ϵ.Moreover, using methods in previous sections, one can see from the factor Γ(−2s) that the Borel transform for the renormalized hard kernel has a linear renormalon at u = 1 2 as well.This renormalon is expected to be canceled by the UV renormalon for linear power corrections O Here, we comment on the strength of the renormalon.Given the above equations, the Borel transform can be extracted and reads (5.9) The residue at u = 1 2 reads 4(u − 1 2 ) . (5.10) Compared with the case for heavy-light Sudakov hard kernel (Eq.(2.47)), the linear renormalon moves to the real part, but in the absence of a factor of π and enhanced by a factor of π 2 .Currently, we do not have a physical understanding of this linear renormalon, but the different transcendentality implies that the physical origin for the linear renormalons in Coulomb and axial gauges might be different.

Discussion and conclusion
Before ending the paper, we would like to make two comments.First, the phase angle in the axial gauge can be heuristically explained by collinear quark exposed in an " imaginary potential" gA z (t, ⃗ created by the space-like gauge link, which is nothing but the analytical continuation of the standard linear potential into the time-like region t > ⃗ x ⊥ .Since the linear potential is known to suffer from the mass renormalon, δV ∼ δm, one expects that correction to the "quark field renormalization" in the external field gA z ∼ iδm takes the form which corresponds to the old-fashioned diagram where an anti-quark is created at t = 0 and annihilates with the incoming quark at a large time t > |⃗ x ⊥ | > 0. Second, we should mention that the observed slow convergence for the imaginary part in lattice extraction of Collins Soper kernel [64,65] can be explained by this renormalon.We hope that by using some form of Borel-resummed hard kernel taking into account the renormalon asymptotics, numerical uncertainties related to imaginary parts can be reduced and the convergence speed of perturbative matching can be improved.In our preliminary numerical test, reduction in the imaginary part has indeed been observed, see Appendix C.
To summarize, we found that the common hard kernel for double-logarithmic factorization formulas in the context of lattice parton distributions, the heavy-light Sudakov hard kernel or axial gauge quark field renormalization, has linear renormalon ambiguity in its imaginary part.In the bubble chain or large β 0 approximation, the strength of the renormalon ambiguity coincides with that for heavy-quark pole mass or heavy Wilson-line linear divergence.For both the threshold expansion of quasi-PDF matching kernel or the TMD expansion of quasi-LFWF amplitude, the linear renormalon anticipates the presence of linear power correction and correlates with UV renormalon at NLP.The quark wave function renormalization in the Coulomb gauge also suffers linear renormalon.
Here H (i) n = n k=1 1 k i are the higher harmonic numbers.The last two formulas are required in the presence of 1 ϵ 3 and 1 ϵ 4 poles.
B Absence of linear renormalon in the light-light Sudakov hard kernel In this appendix, for the convenience of the readers, we show that the standard light-light Sudakov hard kernel for j µ = ψγ µ ψ is free from linear renormalon ambiguity.
We calculate the bubble chain diagram for the vertex correction of the light light Sudakov factor.The diagram dressed up with n bubbles is where the calculation is performed under d = 4 − 2ϵ dimension.u(p) and ū(p ′ ) are the Dirac spinors with p 2 = p ′2 = 0. We define q = p ′ − p and Q 2 = −q 2 .D where G j (ϵ) is analytic at ϵ = 0. Now, using the standard method, the partially renormalized result is The above expression is analytical at u = 1 2 thus the light light Sudakov form factor doesn't suffer from linear renormalon.
The light light Sudakov hard kernel in the large β 0 limit is (B.10) Following the same strategy as the heavy light Sudakov form factor, one can calculate the anomalous dimensions in the large β 0 limit as follows which agrees with Eq. (2.37).The small α expansion (one also replaces β 0 → − reproduces the cusp anomalous dimension in the large n f limit up to NNNLO which is consistent with that from [60] in large n f limit.The single log anomalous dimension is whose small α expansion reproduces the single log anomalous dimension in the large n f limit up to NNNLO which is consistent with that from [60] in large n f limit.Setting β 0 → − n f 3π and µ = Q in Eq. (B.7), one can also obtain the constants term in the large n f limit, which is consistent with that from [60] in large n f limit.

C "Imaginary part" of Collins-Soper kernel
The Collins-Soper (CS) kernel describes the rapidity dependence of the TMDPDF or TMDWF [80,81].It was proposed that the CS kernel can be extracted on lattice through quasi-TMDPDF [82] or quasi-TMDWF [39] with perturbative matching.The CS kernel is purely real by definition.However, a non-vanishing imaginary part was observed during the extraction of the CS kernel through quasi-TMDWF [64,65,83].According to these papers, the ratio of the quasi-TMDWFs is nearly real and the matching correction is a dominant source of the imaginary part of the CS kernel.
To study the imaginary part of the matching correction, we define a R factor following the similar strategy in Appendix C in Ref. [65], R[µ, x, P 1 , P 2 ] = 1 2 ln P 1 /P 2 A ln 4x where A is the phase angle defined in Eq. (2.7).We perform the renormalization group resummation (RGR) for the phase angle A based on Eq. (2.9).The R factors with the RG resummed phase angles are shown as the dashed curves in Fig. 2. The dashed curves deviate from zero, which causes the non-vanishing imaginary part for the CS kernel.Similar observations have been found in Ref. [65].We make use of the asymptotic form Eq. (2.51) to perform the leading renormalon resummation (LRR) [84] for the phase angle A. The R factors with the RGR+LRR phase angles are shown as the solid curves in Fig. 2.After the leading renormalon resummation, the imaginary part of the matching correction is suppressed during the moderate x range.Thus the imaginary part of the CS kernel is expected to be suppressed with leading renormalon resummed matching coefficient.
In the above LRR approach, we regulate the renormalon divergence in Borel plane with the principle value prescription [85].One can choose different contours in Borel plane, corresponding to different results with the renormalon ambiguity at O Λ QCD x(1−x)P 1 − Λ QCD x(1−x)P 2 .According to Eqs. (4.26)(4.27),this renormalon ambiguity is mixed with the NLP TMDWF.Thus, to achieve the leading power accuracy for the extraction of CS kernel through TMDWF, one needs to consider NLP TMDWF.One approach is to parametrize the NLP TMDWF and perform some phenomenological fits with lattice data, which requires further investigation.

D The U †
C (⃗ x) U C (0) correlator in Coulomb gauge The Coulomb gauge quasi-PDF was proposed in Ref. [78], which is expected to be free from the linear divergence and the associated renormalon.The numerical tests in Ref. [78] indicate that there is no linear divergence in the lattice matrix elements.To check the linear renormalon, we calculate the bubble chain diagram for the U † C (⃗ x) U C (0) correlator, where U C (⃗ x) is defined in Eq. (5.1).
The U † C (⃗ x) U C (0) correlator is similar to the self-energy diagram for quasi-PDF in standard gauge and will also be encountered in the TMD objects playing the role of the "gauge-link self-interactions".The bubble chain diagram with n bubbles is where f (ϵ) is defined in Eq. (2.14).The relevant Feynman integral f n (⃗ x) is As before, s = (n + 1)ϵ.Based on the standard procedure, it is free from linear renormalon.
of the heavy-light Sudakov hard kernel 2.1 Bubble-chain diagram for the heavy-light hard kernel 2.2 Renormalization, RGE and anomalous dimensions 2.3 Borel transform and linear renormalon

4
Renormalon cancellation and power correction to TMD limit of quasi-LFWF amplitude 4.1 Bubble chain diagram for quasi-LFWF amplitude 4.2 Power expansion and renormalon cancellation 4.3 NLP soft factor.Linear UV vs rapidity divergences 4.4 Exponentially small contributions and divergence of TMD expansion 5 Linear renormalon of quark wave function renormalization in Coulomb gauge.6 Discussion and conclusion A Second Stirling numbers B Absence of linear renormalon in the light-light Sudakov hard kernel C "Imaginary part" of Collins-Soper kernel D The U † C (⃗ x) U C (0) correlator in Coulomb gauge 1 Introduction

Figure 1 :
Figure 1: Bubble chain diagram for the heavy light Sudakov hard kernel.

Figure 2 :
Figure 2: The imaginary part of the matching correction during extracting the CS kernel through quasi-LFWF amplitude.The dashed curves are the results only with RGE resummation.The solid curves are the results with the RGE and leading renormalon resummation (LRR).We choose the typical parameters under current lattice techniques, µ = 2 GeV, P 1 = 1.72 GeV, P 2 = 2.15 GeV and n f = 3.