TMD PDFs in the Laguerre polynomial basis

We suggest the modified matching procedure for TMD PDF to the integrated PDF aimed to increase the amount of perturbative information in the TMD PDF expression. The procedure consists in the selection and usage of the non-minimal operator basis, which restricts the expansion to desired general behavior. The implication of OPE allows to systematic account of the higher order corrections. In the case of TMD PDF we assume the Gaussian behavior, which suggests Laguerre polynomial basis as the best for the convergence of OPE. We present the leading and next-to-leading expression of TMD PDF in this basis. The obtained perturbative expression for the TMD PDF is valid in the wide region of bT (we estimate this region as bT ≲ 2 − 3 GeV−1 depending on x).


Introduction
Transverse momentum dependent (TMD) parton distribution functions (PDFs) and fragmentation functions (FFs) (we will refer them collectively as TMDs) contain a mixture of both perturbative and non-perturbative contribution. There are no strick rules and prescriptions for the separation of these contribution. Depending on the details on the definition one can add or subtract some parts of perturbative expansion to a non-perturbative contributions. In this paper we suggest new point of view on the perturbative contribution to TMDs. The main goal of our reformulation is to increase the role of the perturbative input in the description of the TMDs.
The TMD factorization theorems give the relation between the TMDs and transversely differential cross-sections in the limit Q 2 ≫ b −2 T , Λ 2 , where Q 2 is the large external momentum, b T is the Fourier conjugated variable to the relative transverse momentum of hadrons, and Λ is a typical intrinsic hadronic scale. As a result of factorization procedure, TMDs depend on x (the longitudinal part of the parton momentum), b T , and scales of the factorization procedure µ and ζ. While the dependence on the scales µ and ζ is known, via the corresponding evolution equations, the x-and b T -dependence are the subject of fitting and modeling. For a theoretical introduction to the TMD factorization theorems and TMDs see [1,2], while for the recent phenomenological review see e.g. [3][4][5] and references therein.
In order to get in touch with the integrated parton distribution function, as well as to increase the perturbative QCD input the additional factorization procedure is applied [6].

JHEP08(2014)089
This procedure is based on the operator product expansion (OPE) at short transverse distances (we will refer it simply as OPE), and it is well-founded for b −2 T ≫ Λ 2 . At extremely small-b T the expansion is saturated by the first term of OPE, which is proportional to the integrated parton distribution. For larger b T the higher terms of the OPE should be taken into account. These terms result to some unknown parton distributions.
Usually, the contributions of higher order terms of OPE are replaced by a single unknown function, which at b T → 0 reduces to unity, see e.g. [7]. This function is called non-perturbative factor. For the details of the non-perturbative factor introduction and its modern status see [4,[8][9][10]. Within the perturbative QCD the non-perturbative factor cannot be obtained analytically, but only extracted from the comparison with data. Frequently, the non-perturbative factor is taken in the form of a Gaussian exponent see e.g. [4,5,8,10], so we may conclude that the contribution of the higher OPE terms is significant. At the same time, the details of the fine structure and the intrinsic scale dependence of non-perturbative factor are unknown and vary between studies, see e.g. discussion in [4], and references therein.
In order to proceed further we should clarify the notions of perturbative and nonperturbative contributions. Generally speaking, within QCD one can perform the small-b T OPE of the TMD operator up to any given α s -and operator-order. Hence, this part of information, namely, the coefficient functions of operators, is entirely perturbative. The truly non-perturbative objects are the hadronic states. Therefore, taking the matrix element of the small-b T OPE we obtain products of perturbative coefficient functions with non-perturbative parton distributions. In such a way, the non-perturbative factor contains a mixture of perturbative and non-perturbative parts. The important question is how to maximize the perturbative input with the minimum number of non-perturbative functions. In the following we suggest a reorganization of the small-b T expansion, which may increase the amount of perturbative information in TMDs.
The standard small-b T OPE is ordered by the powers of b T . Symbolically, it can be written as where O(x, b T ) is a TMD operator, O n is a transversally local operator (operator with only light-cone nonlocality), and ⊗ represents the Mellin convolution with respect to x. The coefficient function G n is proportional to b n T . In the absence of interaction the right-handside of (1.1) represents the Taylor series of the operator O(x, b T ) at b T = 0 (therefore, we mark the operators and coefficient function in (1.1) with superscript T). In this context, the operators O n are proportional to n'th power of transverse derivative, O n ∼ ∂ n .
The Taylor-like OPE is well-founded in the presence of an extreme parameter, e.g. at small b T , or large k T . In these cases the contribution of the higher terms of OPE is under control. In the absence of the extreme parameter, the contribution of the higher OPE terms is uncontrolled.
In the case of middle b T the Taylor-like OPE cannot guaranty the smallness of higher terms. Moreover, it is well-known that at the middle b T TMDs have rapid behavior, which JHEP08(2014)089 cannot be described by several first terms of a Taylor-like expansion. Altogether it shows that the Taylor-like OPE is inefficient for the description of TMDs. And probably, the OPE performed in some other basis of operators, would show better convergence at middle-b T .
The operator basis for the small-b T expansion should satisfy several general assumptions. First of all, the operator basis should be transversally local. Second, the operator basis should be orthogonal, at least in the free theory. This demand is necessary for the universal definitions of the parton distributions. Third, the operators should be defined on the two dimensional plane for the general case, and on the ray from the zero to infinity for the unpolarized case. Additionally, one can impose symmetry or other constraints, which follow from the auxiliary guidelines.
The first two assumptions imply that the basis for the OPE should be chosen in the class of orthogonal polynomials. The third assumption gives some constraints on these set of operators, but do not fix the basis unambiguously. In the absence of the external assumptions we limit ourself to the classical orthogonal polynomials. Within classical orthogonal polynomials there are only two types of polynomials which satisfy the demands [11].
The first one is the Hermite polynomials H n (x). They are orthogonal on the range x ∈ (−∞, +∞), and therefore, suites for the expansion of some general TMD operator with subscript 1 and 2 denoting the components of the transverse vector. Such choice of the operators is not very convenient.
The second one is the Laguerre polynomials L n (x). They are orthogonal on the range x ∈ (0, ∞). Therefore, they suite for the OPE of the direction independent TMD operators n ∼ L n (∂ 2 ). In particular, due to the simple relations between Laguerre polynomials and Hermite polynomials there is a natural relation between expansion (1.2) and (1.3). Namely, averaging (1.2) over angles one obtain (1.3). In this article we consider only the unpolarized TMD PDF, and therefore, we are concentrated on the Laguerre polynomials only. For brevity, we will call the OPE with the operators in the form the Laguerre polynomials, as Laguerre-based OPE.
We want to stress that there is no definite choice of OPE basis. Our choice of the Laguerre and Hermite polynomials is based only on their simplicity and common knowledge of the these polynomials. Nonetheless, it is not a bad choice. In particular, the usage of Laguerre or Hermite polynomials for OPE guaranties the dominating Gaussian behavior of coefficient functions at middle and large b T . This is because the generation function for the Laguerre and Hermite polynomials have Gaussian behavior. The Gaussian anzatz is often used for fitting TMDs and it describes b T -dependence well enough. Therefore, we can expect that the OPE over operators in form of Laguerre polynomials would saturate experimental data by less terms comparing to the Taylor-like expansion.

JHEP08(2014)089
In the both cases, Taylor-like expansion (1.1) and Laguerre-based expansion (1.3), one faces many operators of higher orders, which matrix elements are unknown. Only the leading term in both expansions results to a known distribution, namely, integrated PDF. Therefore, for the practical application one uses only the first term of OPE. Thus, the Laguerre-based version of OPE does not eliminate the necessity to introduce the nonperturbative factor. However, probably, the new version of OPE reduces the significance of the non-perturbative factor, and increases the amount of calculable input.
The most important part of OPE consideration is the loop-corrections to the coefficient function. The corrections give the deviation of the functional form from the free-theory limit. In the Taylor-like OPE the corrections can contain only the logarithms of b T , due to the renormalizability of the theory. In the non-Taylor-like expansion, the other type of corrections are possible, e.g. power corrections. These corrections are of special interest, because they produce the perturbative deviation from the Gaussianity. The observation of such a deviation in the experimental data can justify the usage of the Laguerre-based approach.
In the paper we derive the coefficient function for leading term of the Laguerre-based OPE of the unpolarized TMD PDF operator, and discuss its properties. The paper is structured as follow. In order to introduce the necessary algebraic relation and to discuss the technical details of the Laguerre-based OPE, in section I we consider TMD operator in the toy-example of φ 3 -theory in six dimensions. Then in section II we present calculation of the coefficient function for the Laguerre-based OPE in QCD. The discussion of the behavior of TMD PDF in the Laguerre-based OPE is presented in section II.2.

Small-b T expansion in φ 3 theory
In order to demonstrate the details of the Laguerre-based OPE without excessive technicalities we start from the consideration of OPE in the φ 3 theory in six dimensions. For shortness we refer this model as φ 3 D=6 -model. It has similar to QCD structure of diagrams and with a slight effort of imagination it can be transformed to QCD almost in everything, except gauge links. The action of the model reads An analog of QCD TMD PDF operator is where ξ denotes an arbitrary light-front vector with components {0, ξ − , b T }. The coefficient Z φ denotes the renormalization constant for the field φ, and φ r denotes the renormalized field φ. At small b T the OPE of (2.2) has the form JHEP08(2014)089 and the dots stay for the terms with operators of the higher dimensions. The operators of the higher dimensions include the operators with larger number of derivatives, as well as operators with larger number of fields. The coefficient functions C n are dimensionless, they are functions of x, g(µ) and the logarithms of b T κ. The scale κ in the OPE (2.3) is the scale of the subtractions of the operator singularities. This scale is independent on the renormalization scale of the theory, which is denoted as µ. However, usually it is simpler to set κ = µ. In particular, in the operators (2.4) we suppose the renormalization of the field φ at the scale κ. The difference between the normalization points is placed into the coefficient function. Such a composition of normalization points can be always achieved with the help of corresponding evolution equations, if the difference between κ and µ is small enough.
Our aim is to redefine the operators of OPE (2. 3) in such a way, that the dominating behavior of coefficient functions would be Gaussian. This implies the reorganization of the whole OPE. However, since the matrix elements of the higher terms of OPE are unknown (in contrast to the matrix element of operator O 0 , which is the integrated PDF), our main aim is the leading term of reorganized OPE. In particular, it means that the terms with four and higher number of fields are uninteresting to us. At the same time, we have to consider operators with arbitrary number of derivatives in the sector of two-field operators, because the desired reorganization involves coefficient functions of all such operators.

OPE in the free theory
To start with, let us consider the OPE in the free theory. In the free theory the fields can be interpreted as a classical fields and the OPE at small b T is just a Taylor expansion at b T = 0. Thus, the expression (2.3) in the free-theory reads (2.5) Here we omit the field-renormalization constants and the field labels. The subscript T on the vectors and scalar products denotes the transverse component of the vectors and Euclidian scalar product.
In the following we are interesting only in the unpolarized TMDs. Therefore, we neglect the operators, which do not contribute to the unpolarized case. It can be archived by averaging over the directions of b T . Averaging both sides of (2.5) we obtain

JHEP08(2014)089
where d ⊥ = 4 is the dimension of the transverse space. We will keep d ⊥ arbitrary for the φ 3 D=6 model in order to have clear connection with the QCD case considered in the section 3.
In order to rewrite the Taylor-like OPE via the Laguerre polynomial basis we use the following relation n (x) are the associated Laguerre polynomials of order n and degree d ⊥ −2 2 . The scale B T is some constant scale which needed for the compensation of the dimension in the argument of the Laguerre polynomial. In fact, the right-hand-side of (2.7) is independent on B T .
The expression (2.7) is an algebraic expression, and therefore, it is valid for the interacting fields, as well as for the free theory. With the help of this expression one can rewrite the two-field sector of any Teylor-like OPE via the Laguerre-based operators. Inserting the expansion (2.7) into (2.6) we obtain One can see that, indeed, the dominating behavior of the coefficient functions in Laguerre-based OPE is Gaussian. Therefore, we can expect that such an expansion would saturate faster, comparing with the expansion (2.6). Such an expectation cannot be proved by any perturbative calculation due to the lack of extreme parameters. The only reason for our hope is the experimental observation that the dominating behavior of the TMDs is Gaussian.
The expression (2.8) is independent on the parameter B T . This independence is of the algebraic origin and, therefore, one can derive the relation between the operators with different B T . However, such relation involves operators of different orders, and therefore, it is useless in the absence of information about the higher order terms.
The original OPE can be obtained by reexpanding (2.8) at b T = 0 or by taking the limit B T → +∞. Lowering the parameter B T one takes away the parts from the lower terms of OPE and distributes them between higher terms. In this way, one may think of the variation over B T , as about redistribution of the coefficient functions between the operators.

Laguerre-based OPE at one loop
The consideration of the OPE in the free theory is almost an algebraic exercise. The more interesting subject is consideration of the perturbative corrections to the coefficient functions. We will follow the strategy of the previous paragraph: firstly, we derive the coefficient functions for the Taylor-like expansion at arbitrary order; secondly, using the relation (2.7) we reorganize the series via the operators in the form of Laguerre polynomials. In principal, the perturbative corrections should not be of the same b T -behavior as the leading term. This is because the coefficient functions of the non-Taylor-like OPE can contain non-logarithmical dependence on the parameter of expansion. In this way, the theory dictates the form of b T -behavior, although at the leading order the b T -behavior is almost of our choice.
OPE in the theory of interacting fields contains operators with more then two fields. In order to define the coefficient functions self-consistently, the consideration of operators with the same dimension should be performed simultaneously. Since we are interested in the operators of arbitrary high dimension, we should consider all possible operators. However, at a given order of the perturbative expansion the number of operators is limited. In particular, for the definition of the coefficient functions of operators with two fields at order O(g 2 ) one needs to consider only two-field operators. While, for the definition of the same operators at O(g 4 ) one needs to consider additionally four-field operators (and, hence, four-point Green function of the non-local operator). The details of this analysis are presented in the appendix A.
We limit our-self to the definition of the coefficient functions of the order O(g 2 ). As it is shown in the appendix A, we need to consider only the two-point Green function of the operator O(x, b T ). Since we are interested in the higher-derivative terms and our operator breaks Lorentz invariance, the external fields of Green function should have an arbitrary momentum. Let us introduce the notation for the two-point Green function and its perturbative expansion where α g (µ) = g 2 (µ)/(4π) 3 . Correspondingly, the notation for the two-point Green function of the operators O n of OPE reads Finally, the perturbative expansion of the coefficient function reads where the leading term can be found from the free-theory expression (2.5) is not a part of the coefficient function. Therefore, the coefficient function (2.12) is a dimensionless scalar function. In the appendix A it shown that the coefficient C [2] can be obtain by solving the relation (A.7), which in the notations (2.10)-(2.12) reads The expression for the G [2] is given by a single diagram shown in figure 1, it reads where d = 6 − 2ǫ the parameter of dimensional regularization. The straightforward evaluation of the integral (2.15) gives us wherex = 1−x, and K is the modified Bessel function of the second kind. We bring special attention to the signs of i0's in arguments of the function (2.16), since in the neighborhood of the interesting to us point b 2 T = 0 Bessel functions tend to change its kind. In our case b 2 T > 0 one has where H (1) is the Hankel function of the first order, and the M S scheme is applied.
The expression for G [2] n can be found by independent calculation of the loop-correction to the operator O n or, alternatively, can be derived from the general expression (2.17). In the later case one needs to extract the term with the corresponding (integer) power of p and b T from (2.17) and then take the limit b T → 0.
The limit b T → 0 should be taken with special attention, because it does not commute with the limit ǫ → 0. One can reveal the non-commutation of the limits by considering the series representation for the Hankel function. For the expression (2.17) the expansion JHEP08(2014)089 in the power of p 2 reads [11] G [2] Indeed, taking limit b T → 0, while keeping ǫ > 0, one eliminates all terms of the first sum in (2.18). The second sum contains the poles of ǫ which should be subtracted in the operator of OPE. Taking the limits in the opposite order, i.e. ǫ → 0 and then b T → 0, one does not obtain the ǫ-poles, but obtains the logarithms of b T . The difference between the order of limits results into the untrivial coefficient functions of OPE. Indeed, the solution of (2.14) can be schematically written as where G n denotes the projection of the Green function onto the expression with the same Lorentz composition of b T and p T as in the operator O n . Considering the expression (2.18) one can see that our set of two-field operators (2.5) is incomplete. The set of two-field operators presented in (2.5) contains only operators with derivatives in the transverse directions. Whereas, the expression (2.18) contains the powers of p 2 as well. Therefore, the operators with different combinations of ∂ 2 operators should be added to the OPE. However, there are several observations which allows us to skip the considerations of such operators.
First of all, the operators containing (∂ 2 ) n do not mix with operators which contains (∂ 2 ) k with k < n. This is a consequence of the Lorentz invariance, e.g. an operator with ∂ µ T ∂ ν T can mix with an operator g µν T ∂ 2 , but not the opposite. Therefore, such operators do not influence on the definition of the coefficient function for operators with only transverse derivatives.
Second, the matrix elements of operator with ∂ 2 more suppressed in compare with matrix elements of only transverse derivatives. Indeed, we can implement the counting rules ∂ T ∼ k T , while ∂ 2 ∼ Λ 2 . In the TMD regime k 2 T ≫ Λ 2 , the matrix elements of the additional operators negligible.
Therefore, we omit operators with power of ∂ 2 in OPE, although, in general, their contribution to the TMDs is non-zero. We note that one should keep the non-zero p 2 during the calculation of the Green functions, in order to prevent artificial infrared singularities.
Applying the procedure (2.19) we find the expressions for the coefficient functions where

JHEP08(2014)089
We note that the factor e 2γ E in the logarithm appear due to the finite part of the gammafunctions in the first sum (2.18). Thus, the expression for the Taylor-like OPE has the following form and the dots stay for the two-field operators with ∂ 2 and operators with more fields.
The µ-dependence of the operator (2.23) is given by the equation where at n = 0 the kernel is the DGLAP kernel in φ 3 D=6 , see e.g [12]. The first term in the kernel comes from the subtraction of the ultraviolet-singular term to the operator during the matching procedure (2.19). It coincides (with the opposite sign) with the coefficient infront of ln µ 2 in the coefficient function (2.20). So, in total, the µ-dependence of these terms cancels. The second term comes from the µ-dependence of the field renormalization constant. And therefore, it resembles the renormalization group kernel for the TMD operator.
Averaging over the directions and reorganizing the expansion with the help of (2.7) we obtain: and the dots stay for the two-field operators with ∂ 2 and operators with more fields.
The µ-dependence of the operators O n can be obtained via µ-dependence of the operators O n and the algebraic relation between them. The operators O n satisfy the following equation

JHEP08(2014)089
Taking the matrix element of the operator in the hadronic brackets we obtain the modified small-b T expansion for TMD PDF where f (x, µ) is the integrated PDF and the dots denote terms involving the higher order operators.
The expression (2.28) is singular at b T = 0. The singularity appears because in the regime b T → 0 the TMD factorization should be replace by the collinear factorization. Different factorizations sums up different parts of the original OPE of the hadronic tensor. Therefore, the one should perform additional matching procedure in order to join the regimes of collinear factorization (small-b T ) and TMD factorization (middle-b T ).
The matching procedure can be done by performing OPE on the scale different from µ, e.g. on the scale κ like in the expression (2.3). The scale κ should be kept in the range from Λ 2 and b −2 T in order to prevent the logarithms of (b T κ) be large. Since, in this regime the logarithms of (κ/µ) in the coefficient functions became large and they should resummed with the help of DGLAP equation.
Another way of matching the regimes is to keep the logarithms L T small at small-b T by setting the renormalization scale µ ∼ b −1 T at b T → 0. The most popular choice of the scale µ b was suggested in [7] and reads where b max should be chosen such that α(µ b ) b T →∞ and L T have reasonably small values. Then performing evolution of the TMD PDF down to the scale µ b we obtain This expression is finite at b T → 0 for finite x.
The correction to the leading result is larger at smaller x. One can see that the dominating behavior in the form of the Gaussian factor with x 2 in the argument, would appear in all orders of the perturbation expansion. This indicates that there is a natural bound of convergence for Laguerre-based OPE. The bound is dependent on x and B T and follows from the demand that the correction should be smaller in comparison with the leading term. In the case of φ 3 D=6 this requirement leads to 3 Laguerre-based expansion for unpolarized TMD PDF

TMD operator definition
The gauge invariance of QCD leads to significant complication of the TMD factorization procedure. The soft gluon exchange between hadron states results to additional logarithms JHEP08(2014)089 to be factorized. On the operator level, the gauge invariance is achieved by Wilson lines connecting the quark fields. The geometry of the Wilson lines is dependent on the details of process kinematics. The factorization of the parton distributions related to different hadron states leads to additional divergences, called rapidity divergences, in every TMD operator. The cancelation of the rapidity divergences takes a place inside the product of TMDs and the soft factor. The soft factor is given by the vacuum expectation of the Wilson lines combination and also depends on the kinematics of the TMD process, see e.g. [1,13]. In order to have the divergence-free definition of TMDs, the soft factor can be included into the definition of TMDs operator in the form of the renormalization multiplier [14][15][16].
The definition of the TMD PDF is dependent on the kinematics of the process. In the following we consider the Drell-Yan kinematics and TMD PDF for the hadron propagating along the direction n. Then, the quark TMD PDF reads where µ is the renormalization and hard factorization scale, and ζ n is the scale of rapidity divergences subtraction. The quark TMD operator reads The factor Z q is the field renormalization constant, while the factor S is the rapidity divergences subtraction factor. The factor S in (3.2) can be defined in different ways. Two major definitions are socalled, Collins' definition [1,14], and so-called, Echevarria-Idilbi-Scimemi definition [17]. These approaches are equivalent [16,18], although they are quite different in the form. In the following we are going to use the Echevarria-Idilbi-Scimemi definition, in which the factor S is defined as [16] S where α = 2p + p − ζ n and ∆ is the parameter of the δ-regularization [17,19], see also appendix B. The δ-regularization is used for the regularization of the rapidity divergences and the dependence on the parameter ∆ is canceled in the product (3.2). The expression for theS is given by the following vacuum matrix element The parameters δ ± regularize the rapidity divergences of W (n) and W (n), respectively.

Coefficient functions in the Laguerre basis
In order to obtain OPE in the Laguerre polynomial basis, we use the same strategy as for the φ 3 D=6 model. First, we obtain the Taylor-like expression for the OPE. Second, we reorganize it in the Laguerre polynomials.
In order to obtain the coefficient functions for the operators with two fields at the order O(g 2 ) we need to calculate the two-point Green function where a s = g 2 (4π) 2 = αs 4π . The calculation of the function G [2] q/q is presented in the appendix B. In QCD case OPE contains operators with various Lorentz structures. There are operators with γ − and γ µ T , and various components of derivatives. Since we consider only unpolarized TMDs, we can neglect many of such operators, because they vanish after averaging over directions. The final expression contains three types of operators. They are the operators with γ + (∂ 2 T ) n , ∂ T (∂ 2 T ) n , and γ − (∂ 2 T ) n (there are also operators with ∂ 2 , but they can be omitted due to the same reasons as in φ 3 D=6 , see discussion after (2.19)). The operators belonging to these three types do not mix with each other due to different Lorentz transformation properties. Also the matrix elements of operators with γ T and γ − are suppressed due to the equation of motions (by Λ/p + and (Λ/p + ) 2 respectively). Therefore, we omit these operators at our level of accuracy.
There are also pure gluon operators which contribute to the OPE. The coefficient function of gluon operators can be found by evaluation of the Green function with external gluon fields The amount of various Lorentz structures for the gluon operator is larger. However, the most of operators do not contribute to the unpolarized TMD, and do not mix with the leading operators. Therefore, at our level of accuracy we can average over the polarizations of gluons. The calculation of the function G g/q is presented in appendix B. Thus, the OPE for the TMD operator (3.2) reads where the operators are The dots in (3.7) denote the operators with more then two fields (without counting the Wilson lines), and the operators with other Lorentz structures. The coefficient functions for the individual terms of OPE can be obtained with the help of (2.19). These operators are gauge-invariant in any non-singular gauge, because the parts of operators on the left-and the right-sides of the transverse derivatives are gauge-invariant independently. In the singular gauges, such as the light-cone gauge, transverse derivatives ∂ µ T should be extended to a covariant derivative at light-cone infinity: ∂ µ T − igA µ T (−n∞). These additional transverse gluons are produced by the expansion of the transverse links at light-cone infinity, which make the definition of the TMD PDF gauge-invariant [20].
The coefficient functions of (3.7) can be calculated in the same way as in the φ 3

D=6
theory. The intermidiate expressions for the Green function are presented in appendix B.

JHEP08(2014)089
Averaging over the directions of the b T , and applying the operator redefinition (2.7) we obtain the OPE in the Laguerre basis where the coefficient functions are The evolution equations for the operators O n can be easily derived, and they have the DGLAP structure.

Modified expression for TMD PDF
With the help of the OPE (3.13) we can obtain the modified expression for the TMD PDF, which is presumably valid in some wide region of b T . It reads where the symbol O 1 denotes the order of eliminated contribution. The estimation of the size for O 1 is impossible within the perturbative QCD, but can be obtained from the comparison with experimental data or other theoretical approaches. In the following we suppose that O 1 is negligible in comparison with the first term of (3.16). The function f j/H in (3.16) is integrated PDF. The coefficient function C can be obtained from equations (3.14)-(3.15) at n = 0, At b T → 0 this expressions reveal the standard expressions for the matching coefficients of TMD PDF to integrated PDF. For simplicity we take the size of the Gaussian distribution B T the same for the quark and gluon contributions. However, due to algebraic independence of these parts of OPE, the parameters B T can be different for these distributions. The derived expression can be compared to the standard expression used in the phenomenology (see e.g. [4,[8][9][10]) where g's are unknown functions, and C q/j is the coefficient function for Taylor-like expansion at n = 0 (3.8)-(3.9). The typical expression for the function g is g = −b 2 T P 2 T /4. One can see that expansion over the Laguerre basis at O(a 2 s ) implies the following expression for functions g which should be considered as a part of coefficient function (i.e. in the convolution with integrated PDF). The main difference of this expression from the frequently used Gaussian is the common multiplier x 2 . This multiplier spreads the coefficient function at smaller x. The logarithm factor in (3.20) produces small relative effect in the region b T < B T , where our calculation is applicable. The TMD PDF is dependent on the scale parameters µ and ζ. This dependence is given by the Collins-Soper equation and renormalization group equation. Together they result to the equations On the right-hand side of the expression (3.17) such behavior is not transparent. The µ-dependence of the coefficient function is collected from the general µ-dependence of the TMD operator (3.2) and the µ-dependence induced by OPE. The later cancels between JHEP08(2014)089 the operators and the coefficient functions in the expression (3.13), so the resulting µdependence of the operator (3.2) and its OPE coincide.
In the Taylor-like OPE every subset of operators with the same canonical dimension satisfies the same renormalization group equations as the operator. This is not the case of the Laguerre-based OPE, because the later are constructed from the operators of different dimensions. Therefore, the µ-dependence and ζ-dependence of every individual term of Laguerre-based OPE should not reproduce these dependencies of TMD operator, as it takes a place in Taylor-like OPE. Nonetheless, one can see from (3.13), that the ζ-dependence is universal for all terms of OPE, as well as b T -independent part of µ-dependence. These terms of OPE (they are collected in the first and second line of (3.17), and the expression (3.18) in total) behave in the usual way, i.e. every term reproduces the renormalization group behavior of the initial operator. The µ-dependence of the rest terms (the third line of (3.17)) cancels only in the sum of all operators. Therefore, during the evolution with respect to µ some portion of the truncated expression (3.16) "flows away" to the higher order terms. However, this "lost" part is of order O 1 so it can be neglected in the range of application of expression (3.16).
The OPE (3.13) is performed at such a scale that keeps the logarithms L T small. The common choice of the factorization scale µ is µ b , introduced in (2.29). Then the value of L T is weakly dependent on b T , and small in the region b T b max . The value of b max is typically taken of order 1 GeV −1 , and characterizes the boundary of the perturbative region. Simultaneously, the parameter ζ should be taken such that ln ζ/µ 2 b is also small. The best choice is ζ = µ b .
In practice, the values of parameters µ and ζ are dictated by the TMD factorization theorem, which results to µ 2 ∼ ζ ∼ Q 2 . Therefore, in order to apply the OPE we first need to evolve the TMD PDF from µ 2 ∼ ζ ∼ Q 2 down to µ 2 = ζ = µ 2 b . The evolution is given by equations (3.21), (3.22). We obtain The TMD evolution is independent on x. Indeed, since the kernels of equations (3.21)-(3.22) are independent on x, the evolution exponent in (3.23) is also x-independent. Therefore, we have two free parameters, namely, B T and b max .
The value of the b max can be fixed from two assumptions. First assumption is that a s (µ b ) should not touch the Landau singularity at any b T (this assumption restricts b max from above). Second assumption is that L b T should be reasonably small in some wide region of b T (this assumption restricts b max from below). These assumptions lead to the typical range of b max from 0.5 GeV −1 to 2 GeV −1 .
We have an additional restriction on b max . The calculated correction should be reasonably small in the comparison with the leading term. Since the coefficient function is the generalized function for the comparison we convolute it with some test function. In  The character feature of the result is that the obtained correction is negative. This is the general result, which can be deduced without calculation. Indeed, the main contribution to the coefficient function comes from the term proportional to L T , which, in its own turn, is proportional to DGLAP kernel with the negative sign (we remind that the negative sign is the result of the general renormalization point independence of OPE). Since the DGLAP kernel is strictly positive, the main correction to coefficient function is negative. At some very large b T the power corrections became significant and they can result into the positive expression. In figure 3 we compare our expression with the expression (3.19) at b max = JHEP08(2014)089 The deviation from the standard approach is small for the large-x, even for very large b T . While, for the smaller-x the difference between the expressions (3.16) and (3.19) is more significant. At b T = 2 GeV −1 (which is an optimistic limit of applicability for the expression (3.16)) the deviation is of order 8% for x = 0.1, whereas at x = 10 −2 it is already 16%.
We expect that the next order correction does not bring any significant deviation from the general form of the function g (3.20). One can show that the leading large- This implies only some additional terms to the argument of the logarithm in (3.20). Such corrections would not significantly influence on the region of the applicability of the calculation, because it is mostly defined by the exponent factor.

Conclusion
The TMDs are involved objects which cannot be obtained within the perturbative QCD. For the description of the b T -dependence one should use some non-perturbative information.
In the standard description (see e.g. [1]) the b T -dependence of TMDs is given by unknown function, called non-perturbative factor, which is matched with the integrated PDFs with the help of the first few terms of the small-b T OPE.
We suggest the modification of the standard approach, which consists in the consideration of the small-b T OPE in different operator basis. While the standard approach uses the power (Taylor-like) expansion, we suggest to consider the expansion which leading term reproduces the desired form of TMD. Such approach does not spoil the standard properties of TMDs, such as evolution properties. The main prediction of the approach is the "fine" structure upon the general non-perturbative input.
In this article we assume that at middle-and large-b T the dominating behavior is Gaussian. This is often used in phenomenology conjecture, (see e.g. [3,4,8,9]). In order to reproduces this behavior within the leading term of OPE we uses the operator basis based on Laguerre polynomials (2.9). The Laguerre polynomials are also singled out by the fact that the they are the only classical polynomials with the support b T ∈ (0, ∞), which is needed in the case of unpolarized TMD PDFs. We want to emphasize that there are no theoretical restrictions on the b T behavior and symmetry properties of TMDs, therefore, basis of another (may be non-classical) polynomials might describe the experimental data better. The general structure of consideration of TMDs in another basis is the same as presented in the paper.
The coefficient functions of Laguerre-based expansion have leading behavior C n ∼ (3.14), (3.15)). Neglecting the n > 0 terms of expansion we obtain the perfect Gaussian behavior at the leading term of perturbative expansion, which is modified by the quantum corrections. In the contrast to the Taylor-like expansion where the corrections take the form of logarithms, in the Laguerre-based expansion the corrections change the initial exponential behavior and also contain powers and logarithms of b T . We interpret these corrections as the perturbative modification of the non-perturbative (Gaussian) input. The corrections reproduce the physical effects which were not incorporated in JHEP08(2014)089 the initial simple model. For example, the initial Gaussian behavior is x-independent, but the corrections essentially depend on x.
In the paper we have presented the detailed calculation of OPE in basis of Laguerre polynomials. The main technical details are given by the example of the OPE for TMD operator in φ 3 D=6 theory. We have calculated the coefficient functions at tree and oneloop accuracy for the unpolarized quark TMD PDF, the results are presented in equations (3.13)-(3.15).
The obtained expression for TMD PDF (3.16) has only two free parameters b max and B T . The parameter b max describes the matching of small-b T regime with large-b T regime and can be restricted by the assumption that the size of the perturbative correction should be small in comparison to the leading order. The parameter B T is an unknown parameter of the OPE. In general, the OPE is independent on B T : it cancels in the sum of all terms. However, restricting our consideration by the first term we introduce artificial dependence on this scale. Therefore, the parameter B T is to be found from the comparison with experimental data.
The obtained perturbative results show the general tendency to spread the initial Gaussian input. Moreover the spreading is stronger for smaller x, which is in agreement with the parton picture of hadron. Indeed, the smaller x partons should be distributed in the wider spatial region in comparison to larger-x partons. At one-loop accuracy, the large-b T dominating term of coefficient function is proportional to exp(−x 2 b 2 T )b 4 T . One can show that at two loop order the large-b T dominating term is proportional to exp(−x 2 b 2 T )b 8 T . So, the higher order contributions can violate the Gaussian behavior suggested by leading term of Laguerre-based OPE.
The influence of the non-perturbative corrections grows with the growing of b T . For large b T the perturbative contribution to TMD is negligible, and can only be parameterized by some function. Therefore, the Laguerre-based OPE can be viewed as alternative form of the standard matching procedure, which allows to reduce the significance of the nonperturbative factor at small b T , but does not neglect it.

JHEP08(2014)089
The small-b T OPE is an exact relation of the form where O n,k is the transversely local operator which contains k transverse derivatives and n fields φ. The functions C are dimensionsionless, but they contain logarithms of b T . The expansion (A.1) is schematic, it counts only the orders of dimensions. In order to obtain the coefficient functions C one should consider all possible Green functions of (A.1) and match both sides of the relation. The matching can be done orderby-order in powers of b T and g. With this aim we introduce the notations for the terms of perturbative expansion. The N -field Green function G N (b T ) of the operator O(x, b T ) has the following pertubative expansion where summation runs with step 2. The N -field Green functions R N,n,k of operators O n,k have the following perturbative expansions Matching the perturbative expansion on the both sides of (A.1) we obtain the expressions for coefficient functions. For example, at the zeroth order of b T we have: where a N − 2, and * represents the Mellin convolution. This system is overcomplete, but the existence of solution is guarantied by the existence of OPE itself. Therefore, one can choose any complete subset of equations and use it for the evaluation of the coefficients. Let us find at the relations responsible for the coefficient functions of n = 2 operators with arbitrary number of derivatives. The coefficient functions at the leading orders can be obtained from the equations at (N, a) = (2, 0), they read 2,k * R [1] 3,2,k + C [1] 3,k * R However, the diagrams responsible for G [1] 3,k do not contain loops, therefore, this relation is fulfilled by the coefficient function C 2,k = C [2] 2,k * R 2,k * R [2] 2,2,k . (A.7) We conclude that one needs to consider only the 2-point Green function at one-loop in order to obtain C 2,k . The two-loops expression for C 2,k is already effected by operators with four fields (k > 1): G [4] 2,k = C [4] 2,k * R 2,2,k + C [2] 2,k * R [2] 2,2,k + C 2,k * R [4] 2,2,k + C [3] 3,k * R [1] 2,3,k + C (2) 4,k * R This exercise shows us that the coefficient function for the operators with two-field cannot be defined without simultaneous definition of all other operators. Nonetheless, the situation is not hopeless since the mixture between the operators increases order-by-order. It means that the coefficients C [2] can be obtained by considering the two-fields operators only. While the coefficients C [4] require the consideration of two-and four-fields operators.

B Expression for the diagrams in QCD
In this section we collect the expressions needed for the calculation of coefficient functions. The calculations are performed in the Feynman gauge within the dimension regularization (for the regularization of ultraviolet divergences) and δ-regularization (for the regularization of the rapidity divergences). Since the calculation is done for the general external momentum we do not need any additional infrared regularization.
The δ-regularization is introduced in the fashion which resembles the regularization used in the SCET-calculation [17,19]. Namely, the radiation of a gluon with momentum p by a half-infinite Wilson line (say, pointing in the direction n µ ) results to the integrals of the form d 4 z 1 − e ±ik + τ k + ←→ lim The superscript "plus"("minus") marks the δ-regulator for the eikonal propagator of the Wilson line in the n-(n)direction. Keeping the explicit δ-dependence we regularize the rapidity divergences. The contribution of the soft factor is represented by diagrams in figure 5. Their expression reads

JHEP08(2014)089
The obtained results in the limit b 2 T → 0 can be compared with the results presented in [17]. They are in agreement up to a finite constant, which is result of different infrared regularization. Also we note the difference by factor two (in comparison with corresponded expression in [17]) for the argument of logarithms of δ + δ − in (B.10). This is a result of different normalization of light-cone vectors (we use normalization (nn) = 1, whereas in [17] (nn) = 2 is used).
The expression for the second order of the Green function reads G [2] q/q (x, b T , p) = G [2] A + G [2] B + G [2] C + Z [2] q G [0] q/q S [2] A + S The parameters δ ± should be taken according to prescription (3.3). In the limit p 2 → 0 this expression reads G