Power corrections and renormalons in Transverse Momentum Distributions

We study the power corrections to Transverse Momentum Distributions (TMDs) by analyzing renormalon divergences of the perturbative series. The renomalon divergences arise independently in two constituents of TMDs: the rapidity evolution kernel and the small-b matching coefficient. The renormalon contributions have a non-trivial dependence on the Bjorken variable and the transverse distance. We discuss the consistency requirements for power corrections for TMDs and suggest inputs for the TMD phenomenology in accordance with this study. Both unpolarized quark TMD parton distribution function and fragmentation function are considered.


Introduction
The transverse momentum dependent (TMD) distributions are fundamental non-perturbative objects that appear in many relevant processes at LHC, EIC and e + e − colliders, like Vector Boson Production, Higgs production, Semi-Inclusive Deep Inelastic Scattering, e + e − → 2 hadrons. The factorization theorems which establish definitions of TMD distributions in QCD and/or in effective field theory have been formulated recently in [1][2][3][4], using different regularization schemes.
The perturbative properties of unpolarized TMDs, such as evolution and operator product expansion (OPE) in the regime of small transverse separation, have been deduced by several groups using different frameworks (see e.g. [1,2,[4][5][6][7][8][9]). The explicit direct calculation of the TMD evolution function D at NNLO has been provided in [10,11] and recently it was obtained at N 3 LO [12]. Therefore, nowadays the perturbative knowledge of the unpolarized TMDs parton distribution functions (PDFs) and fragmentation functions (FFs) is comprehensive, thanks to the results obtained by various groups [11,[13][14][15][16][17][18][19] On the contrary, the study of the non-perturbative properties of TMDs has been based mainly on phenomenological arguments which combine the perturbative information on TMDs with their perturbatevely incalculable part [4,8,[20][21][22][23][24][25]. These works have lead to different forms of implementation of TMDs which in general are not easy to compare. For instance, on one hand, the well-known phenomenological considerations of Drell-Yan by [26] and [27] (the so-called BLNY model) implement an ansatz within the standard CSS approach with b * -prescription in the impact parameter space (or b-space). They introduce a set of non-perturbative parameters g 1,2,3 and all these parameters (including the definition of b * prescription) are fundamental for these fits. The same model is also the core of the RESBOS program package [28] which is widely used in applications. On another hand, the implementation of TMDPDFs by [29] does not use b * -prescription. They have found that part of the non-perturbative corrections (essentially to the TMD evolution kernel) are negligible. They were able to describe the same data with a different shape of non-perturbative input parametrized by two parameters λ 1,2 . Fits by other groups that limited themselves to the analysis of Vector Boson Production and Higgs production, are less sensitive to the non-perturbative input (although it is still necessary) [30,31]. Additional problems arise in the consideration of TMDFFs which are known to have very different and/or incomparable (in comparison to TMDPDFs) non-perturbative input.
This work is devoted to the study of the leading power corrections to TMD distributions. With this aim we perform an analysis of the leading renormalon structure of TMD distributions. A renormalon analysis of the perturbative series gives an important check of theoretical consistency for any phenomenological ansatz, although it cannot give too stringent restrictions on the fitting parameters. The study of renormalon poles allows to understand the asymptotic behavior of the perturbative series and to deduce the form of the leading non-perturbative corrections [32][33][34][35].
An explicit analysis of the renormalon structure for TMDs has never been done to our best knowledge, although assumptions on its structure were used even before the actual field-theoretical definition of TMDs. We refer here, for instance, to the seminal work of [33] about the Sudakov factor in differential cross-section which is usually referred to justify a Gaussian behavior for the non-perturbative part of the TMD evolution kernel [25]. In order to describe this effect in the modern TMD framework, we recall that the definition of TMDs requires the combination of the Soft Function matrix element with the transverse momentum dependent collinear function. As we show in this work, the renormalon divergences arise in the perturbative consideration of both of these functions. These renormalon contributions have different physical meaning and should be treated independently. Firstly, the renormalon divergence of the soft factor results to a power corrections within the TMD evolution kernel, which are strictly universal for any TMD due to the universality of the soft factor itself. The leading power correction that we derive here is quadratic, that actually suggests the Gaussian shape of the evolution kernel. Secondly, the renormalon divergences naturally arise within the coefficients of the small-b OPE. A study of those contributions gives access to the next-twist corrections of small-b matching, and specifies the shape and the general scaling of TMD.
The paper is built as the following. We provide the necessary concepts and definitions in Sec. 2. In Sec. 3 we perform the calculation of various TMD constituents (such as anomalous dimension and coefficient functions) within the large-β 0 approximation. In the end of this section we provide a collection of the main lessons that follows from our results. The impact of the renormalon divergences on the perturbative series and renormalon subtracted series are studied in Sec. 4. One of the main outcome of the study, namely a consistent ansatz for TMDs is presented in Sec. (4.2).

Notation and Basic Concepts
Throughout the paper we follow the notation for TMDs and corresponding functions introduced in [14]. The quark TMDPDFs and TMDFFs are given by the following matrix elements where R q and Z q are rapidity and ultraviolet renormalization constants, q are quark fields and W T are Wilson lines, and ξ = {0 + , ξ − , b}. The TMDs depend on the Bjorken variables (x for TMDPDFs and z for TMDFFs), the impact parameter b and the factorization scales ζ and µ. The considerations of the TMDPDF and TMDFF are similar in many aspects. Therefore, in order to keep the description transparent we mostly concentrate on the case of the TMDPDFs, while the results for TMDFFs are presented without derivation. The dependence on the factorization scales is given by the evolution equations, which are the same for TMDPDF and TMDFF, namely Through the article we consider only the quark TMDs, therefore in the following we suppress the subscript q on the anomalous dimensions. The values for both anomalous dimensions can be deduced from the renormalization constants [14]. Also γ and D are related to each other by the cross-derivatives where Γ cusp is the honored cusp anomalous dimension. The solution of the evolution equations Eq.
where R is the evolution kernel, .
The final values of scaling parameters is dictated by the kinematic of the TMD cross-section. The variable ζ f ∼ Q 2 (with Q being a typical hard scale) is the scale of the rapidity factorization, and the variable µ f is the scale of hard subprocess factorization. The intriguing point is that the evolution kernel R is not entirely perturbative, but contains a non-perturbative part. An estimate of the non-perturbative contribution to R is necessary in order to obtain the cross section in the momentum space where it is usually measured. The non-perturbative part of the evolution kernel is encoded in the D-function which can be obtained from the rapidity renormalization constant R q . The definition of the rapidity renormalization constant differs from scheme to scheme. In this work we use the δ-regularization scheme defined in [10,14]. In this scheme, the δ-regularization is used to regularize the rapidity divergences, and the dimensional regularization regularizes the rest of divergences. Such a configuration appears to be very effective for the TMD calculus. In particularly, the rapidity renormalization factor R q is expressed via the soft factor S as R q = S −1/2 [14]. In the coordinate space the soft factor is given by the following matrix element where we explicitly denote the ordering of operators and S T are Wilson lines, as defined in [10]. Considering the relation between renormalization constants one can show [10], that where l δ = ln µ 2 /|δ δ δ| . Eq. (2.8) can be used as the formal definition of the TMD evolution function D. In this way, a non-perturbative calculation of the SF gives access to the non-perturbative structure of D. The soft function is perturbatively universal for both Semi Inclusive Deep Inelastic Scattering and Drell-Yan type processes. Therefore, the perturbative part of the anomalous dimension D is universal for TMDPDF and TMDFF. One can also expect its universality in the non-perturbative regime. The TMDs are entirely non-perturbative functions. They cannot be evaluated in perturbative QCD, due to the non-perturbative origin of hadron states. The main subject of the paper is the dependence of TMDs on the parameter b which is generically unrestricted, since it is a variable of Fourier transformation. However it is interesting and numerically important to specially consider the range of small b (here and later b = √ b 2 ). In this range the TMDs can be matched onto corresponding integrated parton distributions. At the operator level the small-b matching is given by the leading term of the small-b OPE. The small-b OPE is a formal operator relation, that relates operators with both light-like and space-like field separation to operators with only light-like field separation. It reads where C n are Wilson coefficient functions, the µ b is the scale of small-b singularities factorization or the OPE matching scale (for simplicity we omit in Eq. (2.9) other matching scales included in the definitions of each component of this equation). Generally, the operators O n are all possible operators with proper quantum numbers and can be organized for instance according to a power expansion, i.e. twists. In this case the matching coefficients behave as where f is some function. The value of the parameter B is unknown, and its origin is entirely non-perturbative. In other words, the unknown scale B represents some characteristic transverse size of interactions inside a hadron B O(1GeV). In practice it is reasonable to consider only the leading term (n = 0) of Eq. (2.9) for b B . In this case f is an integrated parton distribution (or fragmentation function), and coefficient function is called the matching coefficient. So far, the power suppressed terms in Eq. (2.9) has been not considered, to our best knowledge.
For completeness we recall here the renormalization group properties of the TMD Wilson coefficients that we use in the following sections. The evolution equations for the matching coefficients (at µ b = µ) with respect to ζ is where f = q, g species, C f ←f are the matching coefficients on PDFs. It is practically convenient to extract the ζ-dependence from the matching coefficient. We introduce the notation Here and further we use the following notation for logarithms The ζ-free coefficient functionĈ satisfies the following renormalization group equation where the kernel K is and P (x) is the splitting function (DGLAP kernel). The matching coefficient for TMDFF C f →f satisfies the same set of evolution equation with only substitution of PDF splitting function P (x) by the FF ones, P(z)/z 2 [14]. Using these equations one can find the expression for the logarithmic part of the matching coefficients at any given order, in terms of the anomalous dimensions and the finite part of the coefficient functions. The expressions for the anomalous dimensions, the recursive solution of the RGEs and the explicit expressions for the coefficients C and C can be found, e.g. in [14].

TMD in large-β 0 approximation and renormalon divergences
The leading non-perturbative contribution to the perturbative series is commonly associated with renormalons. The renormalon contributions were intensively studied for various matrix elements and in different regimes, for review see [36,37]. A typical signature of renormalons is the factorial divergence of the perturbative series. These divergences are often discussed in terms of the corresponding singularities in the Borel plane. The best representative, and the only stable way to study the renormalon divergence within perturbative QCD, is the large-β 0 approximation. The large-β 0 expression can be obtained from the large-N f expression through the procedure of "naive Abelianization" [38,39]. In this section we present the calculation of large-β 0 correction to TMDs. Since the technique of large-N f calculus is well-known, we skip the detailed evaluation (redirecting the reader to the related literature) and present only intermediate expressions.

The soft function in the large-β 0 approximation
The soft function matrix elements is a key structure for the TMD construction and as such it is a good starting point for the renormalon analysis. The large-β 0 calculation of the soft factor runs in parallel to the calculation of the integrated soft factor for Drell-Yan, which is presented in [32] (see Sec.5.3). Here we present our results of the evaluation.
To begin with, we evaluate the large-N f contribution to the soft factor, which is given by the "bubble" resummed diagram, shown in Fig.1.A. The expression for the (renormalized) diagram with n-bubble insertion is the being parameters of rapidity regularization for Wilson lines pointing in n(n)-direction [14]. The function G is a standard function that appears in large-β 0 calculation [32,36,38,39], and is given by the expression Here, the Euler-Mascheroni constant is a result of the MS scheme. For n = 1, 2 this expression agrees with the direct calculation of the soft factor in δ-regularization [10]. We also introduce an additional function for the double-pole part The functions G andG have the following Taylor series These expressions define the coefficients g k and G j . Note, that g k and g [1] k =g [1] k . The procedure of "naive Abelianization" consists in the replacement of N f by the corresponding β 0 expression [38], i.e.
In this way, we obtain the large-β 0 expression for the soft factor where we have introduced the large-β 0 coupling constant Note, that in Eq. (3.7) the terms suppressed in are dropped. Eq. (3.7) gives access to the anomalous dimension D, which we study in Sec. 3.3, and to the rapidity renormalization factor R q . The factor R q (we recall that it is equal to R q = S −1/2 in the δ-regularization [14]) from the perspective of the large-β 0 approximation has the same perturbative combinatorics as the one-loop-truncated pertrubation series. It is given by at δ − = ζ/p + and SF given in Eq. (3.7). This expression is used in the next section to extract the large-β 0 expression of the Wilson coefficients of small-b OPE.

The TMD in the large-β 0 approximation
To obtain the TMD matching coefficient one should evaluate the diagrams B and C, which are shown in Fig. 1. The result for the sum of these diagrams and their Hermitian conjugations is where we have used the same notation as in Eq. (3.1) andx = 1 − x. The last term in square brackets represents the rapidity divergence which appears in the diagram C. For n = 0, 1 this expression reproduces the result of explicit calculation made in [13]. Using Eq. (3.8) and Eq. (3.9) we can complete the result for the large-N f expression of the TMDPDF, Here, we observe the cancellation of the rapidity divergences that leaves the residual l ζ dependence. In order to extract the matching coefficient of the TMDPDF onto the PDF one has to proceed to the renormalization of Eq. (3.10). This is greatly simplified in the δ-regularization scheme, where all virtual graphs and integrated graphs are zero. The only non-zero contribution is the UV counterterm which is a pure -singularity. The accounting of this part eliminates terms singular in , leaving the finite part unchanged. The latter provides the coefficient function. Performing the "naive Abelianization" as in Eq. (3.6) we obtain the large-β 0 result The additional variable in the square brackets for the functions g indicates the modified value of B B B µ to be substituted. The calculation of TMDFFs matching coefficient proceeds in the same way as for TMDPDFs.
The result of the calculation is One can see that the expression for TMDPDF Eq.  .12) is specific for TMDFF and it is an effect of the expansion of the normalization factor z −2 .
One can check that at n = 0, 1 the expressions (3.11) and (3.12) coincide with the one calculated in [14].

The TMD anomalous dimensions at large-β 0 and renormalon singularities of D
In the articles [10,14] it was shown that in the δ-regularization scheme the anomalous dimension D can be obtained from the rapidity singular part of the soft factor as in (2.8). Considering the Eq. (3.7) we obtain the anomalous dimension D in the large-β 0 approximation The first term in the brackets of Eq. (3.13) behaves ∼ n! at large n, and represents the renormalon singularity. At this point it is convenient to consider the Borel transformation of the result. We define the Borel transformation of a perturbative series in the usual way (3.14) A perturbative series is Borel summable if an integral exists. Performing the Borel transformation on the D function and applying Eq. (3.15), we find The first term is analytical and reproduces the cusp-anomalous dimension at large-β 0 [32] The function which appears in the second term  G(0, −u)).
There are multiple possibilities to define the sum Eq. (3.13), e.g. one can slightly shift the integration contour for Eq. (3.16) into the complex plane. The difference between integrals passing from the lower and upper sides of poles is called infrared (IR)-unambiguity and is given by a (−π) times the residue at the pole. For the anomalous dimension D it reads where c = πC F 2β 0 e 5 3

1.2.
(3.20) The IR-unambiguity represents the typical scale of the error for perturbative series.
The anomalous dimension γ V can be extracted from the coefficient function Eq. (3.11). We consider the derivative of coefficient function at l ζ = L µ where we have dropped the mixing among flavors. The DGLAP kernel at large-β 0 is given by the expression Considering the derivative of Eq. (3.11) and comparing right and left hand sides of Eq. (3.21) we obtain This expression contains no singularity, and hence it is renormalon-free, as it is usually expected for an ultraviolet anomalous dimension.

TMD matching coefficient at large-β 0
Before the evaluation of the sums in Eq. (3.11-3.12) we extract the part related to the anomalous dimension D to obtain the coefficientsĈ defined in Eq. (2.12). This procedure is important since the function D contains its own renormalon singularities, as described in Eq. (3.19). The contribution of D is easily recognized in the third lines of (3.11-3.12) (compare with Eq. (3.13)).
The result of the Borel transform for the coefficientĈ, Eq. (3.11) iŝ where by bold font we denote the Borel transformed functions, The terms in Eq. (3.24) are collected such that every bracket is finite at u → 0. The expression for TMDFF coefficient functionĈ can be obtained using the crossing transformation (x → z −1 ) and the addition of the normalization contribution (the last line in Eq. (3.12)). In the last two lines of Eq. (3.24) we have the infrared renormalon poles at u = 1, 2, ... One can see that the third line contains only first order poles, while the last line contains second order poles at G(0, −u)ψ(−u). Considering the infrared unambiguity at u = 1 we obtain where constant c is given in Eq. (3.20). The x−dependence of this expression exactly reproduces the x−dependence of the leading terms of the next power correction in small-b OPE, see detailed description in [40]. The consideration of unambiguites of higher renormalon poles gives access to the higher-power corrections. We obtain However, these expressions can be modified by the infrared renormalon contributions of the highertwist terms. The most important information of the higher-power corrections is that the renormalons scale as xb 2 , but not as b 2 which is a naive assumption. The consequences of this fact are discussed in the next sessions. The corresponding calculation for TMDFF gives which is the same as Eq. (3.26) with the crossing change x → 1/z. One can see that the difference in normalization which spoils the crossing between TMDPDFs and TMDFFs, disappears in the renormalon contribution. The higher poles unambiguites are provided using the crossing relation x → 1/z in Eq. (3.27).

Lessons from large-β 0
The expressions in Eq. (3.19, 3.26, 3.28) are one of the main results of this work. These expressions represent the leading power correction to the small-b regime, where all perturbative properties of TMDs are derived. These expressions give access to a general structure of the next-to-small-b regime. The practical implementation of results Eq. (3.19, 3.26, 3.28) is given in the next section, while here we collect the most important observation that follows from the large-β 0 calculation and which should be taken into account for TMD phenomenology. The first, and the most obvious, observation is that the leading renormalon poles in the Borel plane appear at u = 1, and not u = 1/2 as it could possibly be. That implies that there is not a correction ∼ √ b 2 at small b. In turn, it implies that an exponential decay of the TMDs that is sometimes suggested in phenomenological studies [41,42] can in no way affect the small-b region. Such type of corrections cannot overlap with the small-b region although may occur in the large-b regime, which is unaccessible by perturbative considerations.
Second, one can see that the renormalon corrections to TMDPDFs matching coefficient scales like xb 2 , and not as simply b 2 (as it is usually assumed), nor as x 2 b 2 (as suggested by Laguerre polynomial decomposition [9]). Therefore, the contributions of higher twist terms in small-b OPE for TMDPDF are largely functions of xb 2 . Correspondingly, TMDFFs matching coefficients are function of b 2 /z. This is important in respect of phenomenological implementation of the TMDs. For instance the b * -prescription which is often adopted does not respect this scaling and so, in this sense, it is not fully consistent with the estimated higher twist effects.
Third, the renormalon contributions to the anomalous dimension D and to matching coefficients have different physical origin and do not mix with each other. In fact, the anomalous dimension D is an universal object that is the same for all regimes of b and for TMDs of different quantum numbers [25]. Thus, the renormalon contribution to D represents a generic universal non-perturbative contribution, alike in the case of heavy quark masses. On the other hand, the (infrared) renormalon divergences within the matching coefficients are to be canceled by the corresponding (ultraviolet) renormalon contributions of higher twists. Therefore, while Eq. (3.19) represents a size of a universal non-perturbative contribution, the Eq. (3.26, 3.28) give the form of the twist-four contribution to small-b OPE. In other words, Eq. (3.26, 3.28) very accurately estimates the x-behavior of subleading correction to small-b OPE.
The consideration of the anomalous dimension D for gluon distributions is identical to those of quarks (apart of trivial replacing of common factor C F by C A ). Contrary, the calculation of the renormalon contribution for gluon and quark-gluon matching coefficient is much more complicated than the one presented here, and is beyond the scope of this paper. In general we can expect a nontrivial dependence of the renormalon contribution on the Bjorken variables. At present, we cannot find arguments which suggest a location for the renormalon poles and an xb 2 scaling different from that of quarks.

Renormalon substraction and power corrections
Our analysis is limited to the quark TMDs only, nonetheless we can advance some considerations on possible inputs, which are consistent with our findings and evaluate their impact on the nonperturbative structure of TMDs. The suggested ansatz for TMDs does not pretend to be unique and moreover is inspired by other popular models. We postpone to a future publication a more dedicated study on the subject.
We recall here that the form of the TMDPDFs which emerges at small-b is where the evolution kernel R is given in Eq. (2.6). The argument ζ b of R is collected from the combination of two exponents: the original factor R , Eq. (2.6) and the exponential prefactor ofĈ, Eq. (2.12), and it takes the value The analogue equation for TMDFFs is obtained replacing consistently the PDF f j←N by the fragmentation function d j→N and the coefficient function C q←j by C q→j , while the evolution kernel remains the same. This expression is usually taken as an initial ansatz for TMD phenomenology.
As we pointed earlier there are two places where the non-perturbative effects arise. The first one is the evolution kernel D which is a part of the evolution prefactor R, and it is common for all TMDs (TMDPDFs and TMDFFs of various polarizations). The second one is the higher twist corrections to the small-b OPE. These non-perturbative contributions are of essentially different origin and should not be mixed. In particular it is important to realize that the non-perturbative contribution of D enters Eq. (4.1) as a prefactor, while the higher order terms of OPE are added to the convolution integral. Therefore, the structure of non-perturbative corrections to TMD that we keep in mind is the following Here, D N P is the non-perturbative addition to the anomalous dimension D, and f N P is the cumulative effect of the higher twist corrections to the small-b OPE. At small (perturbative) b, the non-perturbative parts should turn to zero, such that Eq. (4.2) reproduces Eq. (4.1). In the following subsections we construct a minimal non-contradicting anzatz for TMD distributions that respects the study of large-β 0 approximation.

Non-perturbative corrections to the anomalous dimension D
The non-perturbative part of the anomalous dimension D is one of the most studied in the literature and the one for which a general consensus is achieved. It is assumed to have the Gaussian form at least in the moderate-b region. As we show in Eq. (3.19) the Gaussian form is also suggested by the large-β 0 approximation. A more subtle issue concerns the amount of non-perturbative correction to D, which can be very different depending on the implementation of the TMDs. A check of the renormalon contribution, as provided in this section, gives an estimate of such correction and it is so useful for practical implementations. Let us present the perturbative series for D in the form where d n ∼ n!g can be obtained from Eq. (3.13) and δ n is the large-β 0 suppressed part. The numerical comparison of the large-β 0 expression Eq. (3.16) and the exact expression for D is given in the Tab. 1. One can see that generally the large-β 0 expression overestimates the exact numbers, which is typical for this approximation.
In order to study the properties of the large-β 0 series we introduce a function for its partial sum  We investigate the convergence of the partial sums of M N to its Borel resummed value M ∞ , in order to find the scale at which the non-perturbative corrections associated with renormalons become important. In Tab. 2 the numerical values of partial sums at µ = 10 GeV and at several values of b are presented. The graphical representation for these values is shown in Fig. 2. The convergence of the series is perfect (in the sense that it converges at M 7 that is far beyond the scope of modern perturbative calculations) for the range of b 2 GeV −1 , it becomes weaker at b ∼ 3 GeV −1 , and it is completely lost at b 4 GeV −1 . These are the characteristic scales for switching the perturbative and non-perturbative regimes in D. In other words, the perturbative series can be trustful for b 2 GeV −1 , but it completely looses its prediction power for b 4 GeV −1 . The number N at which convergence is lost depends on the value of µ, however the interval of convergence in b is µ-independent, e.g. at µ = 50 GeV the series converges to M 8 in the region b 2 GeV −1 , but again looses stability at ∼ 4 GeV −1 .
In order to proceed to an estimate of the non-perturbative part of D we write it in the form where D P T is given by the perturbative expression at µ 0 scale, D N P encodes the non-perturbative part. The parameter µ 0 depends on b and should be selected such that a s (µ 0 ) is a reasonably small number. The non-perturbative part D N P is independent on µ (since the evolution part of D is renormalon-free) but depends on the choice of µ 0 . In principle, the best value of the parameter µ 0 can be extracted from the large-β 0 calculation. Indeed, the resummation of bubble-diagrams modifies the coupling in the interaction vertex, such that a loop integral appears to be naturally regularized in the infrared region. Practically, the effect of such resummation can be presented as a freezing of the coupling constant at large b. Particularly popular is the b * prescription [43] defined as At large b the parameter µ 0 approaches C 0 /b max , which should be chosen much less then Λ, i.e. b max C 0 /Λ ∼ 4 GeV −1 . For large-b (say b 3 GeV) the non-perturbative part of D dominates the perturbative one. The large-β 0 calculation allows to estimate the leading contribution (from the side of small-b's) to D N P from the infrared unambiguity Eq. (3.19), Since the large-β 0 approximation overestimates the exact values this number can be considered as an upper bound for non-perturbative input.
To estimate the parameters of the D more accurate, we consider a kind of renormalon subtraction scheme for the anomalous dimension D. We construct a renormalon subtracted expression D(µ, b) = D RS (µ, b) by explicitly summing the large-β 0 contribution in Eq. (4.3) The scale µ here should be chosen such that the logarithm L µ is reasonably small, otherwise the large-β 0 expansion is significantly violated. Using the model Eq. (4.10) we fit the parameters of Eq. (4.6) at µ = 10 GeV in the range b < 3 GeV, with g D = constant≡ g K , at all known perturbative orders. It appears that the result is very stable with respect to b max whose best value we find to be Concerning the non-perturbative part, it appears to be lower then the rude estimation Eq. (4.9) and actually consistent with 0, This value is generally smaller then the typical values presented in the literature, e.g. Ref. [8] quotes g K 0.17GeV 2 , Ref. [44] quotes g K 0.045 ± 0.005GeV 2 . But Ref. [29] finds g K consistent with 0, which agrees with the present findings. However, one should take into account that contrary to standard fits, the present considerations are purely theoretical. Moreover in fits with experimental data, one should consider the extra non-perturbative part of the TMD distribution itself (which is discussed in the next section).
Finally we comment on the possibility of a more sophisticated renormalon subtraction scheme as in the MSR scheme of [45]. In this scheme one provides a subtraction of the renormalon from a perturbative series which depend on an additional scale µ R . The new renormalon subtraction scale can result into large logarithms which, in turn should be resummed. Such a consideration can result to a more accurate restrictions on parameters.

Renormalon consistent ansatz for TMDs
The non-perturbative corrections to the matching coefficients are necessary for all analysis which include low energy data. These corrections have not been deeply studied in QCD and up to now only a phenomenological treatment has been provided. The renormalon calculation of the previous sections allows to formulate a consistent ansatz for the coefficient, separating the effects of the evolution kernel. The leading twist-four contribution to the small-b expansion has a form where the LO coefficient function of the renormalon contribution was calculated in Sec. 3.4 and reads The constants g in,out q are of order cΛ 2 within the large-β 0 approximation, however the actual value should be estimated from data. The non-perturbative scale Λ is the same as in the case of the evolution kernel. The renormalon correction very accurately restores the x-dependence for the twist-four corrections (see detailed explanation e.g. in [37,40,47]).
At larger values of b Eq. (4.14) is corrected by the higher orders of the OPE, and at a particular scale B (which defines the convergence radius of small-b OPE Eq. (2.10)) it is replaced by a single entirely non-perturbative function. Particularly, it is popular to assume a Gaussian behaviour at large-b for this function. A viable model, which takes into account the Gaussian suppression of the TMDs at large-b and consistent with small-b regime (4.14) can have the form Here, we took into account that according to our analysis the natural scale of power corrections to coefficient function is xb 2 (or b 2 /z), which is reflected in the form of exponents. In the figures 3 -5 we illustrate several features of the renormalon consistent model that we propose. In all the plots we fix the µ scale at the value µ = µ * = C 0 /b * . In Fig. 3-left we show that the change ofF with respect to the perturbative order of matching coefficient. On the right hand side of Fig. 3 we show the dependence on the choice of the scale b max , which we find very mild for 1 GeV −1 b max 2 GeV −1 .
The shape of the TMDs can strongly depend of the values of the non-perturbative constants g b,q for b ≥ 2 GeV −1 as shown in Fig. 4-5. Therefore we expect that the values of these constants can be fitted easily. The values that we have assigned in the plots have been inspired by the fit in [29], however they can also change in a real fit with the present model.
Finally we study the x-dependence ofF . For b ≤ 1 GeV −1 the non-perturbative model does not really affect the behavior of the TMD. In Fig. 5 we show instead that for instance at b ∼ 1.5 GeV −1 the model parameter can start to have their impact.
To conclude this section, we observe that in the literature we have not found any non-perturbative input for TMDs fully consistent with Eq. (4.13). For instance the standard b * prescription which is used in many phenomenological analysis [8,25,28] is inconsistent with Eq. (4.13). In fact, using the standard the b * -prescription the higher-twist corrections are simulated by replacing b → b * , and including an additional non-perturbative factor aŝ one does not reproduce Eq. (4.13). The main difference comes from the renormalon scaling, xb 2 vs. b 2 . On top of this the non-perturbative exponent that modulates the TMD at large-b is positioned outside of convolution integral (although it is commonly taken to be x−independent).

Conclusion
In this work we have studied the non-perturbative properties associated with renormalons for the soft function and unintegrated matrix elements. With this aim we have evaluated all constituents of TMD distributions (soft factor, matching coefficient and anomalous dimensions) within the large-β 0 approximation. The (factorial) divergences of the large-β 0 series are associated with the renormalon contribution and allow to estimate the leading non-perturbative contributions. We have found two independent renormalon structures in the perturbative description of TMD: the soft function and small-b matching coefficients. The consideration of the soft function allows to fix the power behavior of the evolution kernel of TMDs. We show the evidence of infrared renormalons at u = 1, 2, .. (u being the Borel parameter). Our results agree with the popular assumption about a Gaussian behavior of the evolution kernel. However the impact of the non-perturbative corrections is estimated to be not very significative for experiments where TMDs are evaluated at scales higher then a few GeV.
The nature of the renormalon contribution to the evolution kernel is peculiar, in the sense that it is generated by the non-perturbative part of a matrix element. In some aspects this is very similar to the renormalon contribution to heavy quark masses. We have discussed also an ansatz which implements a consistent renormalon subtraction for the TMD evolution kernel, which can be useful for phenomenology.
The most promising conclusion of the paper comes from the analyisis of the renormalon contribution to the small-b expansion of TMDs. The discussion of these results can be found in Sec. 3.5. We demonstrate that the power corrections to small-b behave as a function of xb 2 for TMDPDFs and as b 2 /z for TMDFFs. This observation should have a significant impact on the joined TMDPDF -TMDFF phenomenology. Additionally, the large-β 0 computation unveils the form of x−dependence for the leading power correction to the small-b matching. This behavior should be incorporated in realistic and consistent models for TMDs.
We have discussed a possible way to incorporate all collected information into a viable phenomenological ansatz, which in some aspects is inspired and similar to others in the literature. However, we find that generally in the literature the assumptions on the non-perturbative part of TMDs are inconsistent with our conclusions, mainly, due to the naive assumption that these corrections are largely functions of b 2 (contrary to xb 2 ). We postpone to a future work the fit of available data using the present results.