Calculation of transverse momentum dependent distributions beyond the leading power

We compute the contribution of twist-2 and twist-3 parton distribution functions to the small-$b$ expansion for transverse momentum dependent (TMD) distributions at all powers of $b$. The computation is done by the twist-decomposition method based on the spinor formalism for all eight quark TMD distributions. The newly computed terms are accompanied by the prefactor $(M^2b^2)^n$ and represent the target-mass corrections to the resummed cross-section. For the first time, a non-trivial expression for the pretzelosity distribution is derived.


Introduction
Transverse momentum dependent (TMD) distributions extend the parton model, including the transverse motion of hadron's constituents. Any TMD distribution is a function of two dynamical variables x and b. The variable x is the fraction of hadron's momentum carried by the parton. The variable b is the transverse distance that is Fourier conjugated to the transverse momentum of parton k T . In the limit b → 0, which corresponds to integrated or unobserved transverse momentum, a TMD distribution turns to the corresponding collinear parton distribution functions (PDFs) or fragmentation function (FFs). Technically, this relation, which is also known as the matching between TMD parton distribution functions (TMDPDFs) and PDFs (or TMD fragmentation functions (TMDFFs) and FFs), is received by the operator product expansion, and its leading power term is very well studied. In the present work, we extend this formalism beyond the leading power term and compute matching of TMDPDFs to PDFs of twist-2 and twist-3 at all powers of b 2 -expansion.
The matching of TMDPDFs to PDFs is an important part of the TMD factorization approach. The review of various aspects of TMD factorization can be found in [1][2][3]. On the theory side, the matching establishes the connection with the resummation formalism and allows interpolation between TMD factorized cross-section and fixed order computations, see f.i. [4][5][6]. On the phenomenological side, the matching essentially reduces the parametric freedom for TMD ansatzes. The modern phenomenology of TMD distributions is grounded on the matching relations and demonstrates perfect agreement with the large amount of experimental data [7,8].
So far, all studies of matching relations were restricted to the leading power term only. The leading power term is the most simple and numerically dominant contribution. Nonetheless, several aspects make the study of power corrections interesting. First of all, such a study carries a significant amount of methodological novelty. Indeed, the power corrections are generally considered as a complicated field, and their computation is an interesting theoretical task. In this work, we have computed the whole series of power correction with PDFs of twist-2 and twist-3, which is almost unprecedented. Methodologically, the closest example of similar computation is the computation of kinematic power corrections to Deeply Virtual Compton Scattering (DVCS) made by Braun and Manashov in ref. [9]. The second point of interest is the derivation of the matching relations for polarized TMD distributions. For many TMD distributions already, the leading power matching involves twist-3 functions and requires a non-trivial computation. The computations for different TMDPDFs have been made by different methods in refs. [10][11][12][13]. In ref. [14], all polarized TMDPDFs were systematically computed in a single scheme, and the agreement with previous computations had been shown. However, all these computations were unable to find the non-trivial matching of the pretzelosity distribution, which is derived for the first time in this work. The third point of interest is the comparison of the derived power corrections to the extracted ones. There are many examples, where the part of the power correction proportional to twist-2 PDFs (Wandzura-Wilczek approximation) is numerically dominant [15]. However, there are also known cases of opposite behavior. In this work, we demonstrate that TMD distributions belong to the latter case, and for them, the contribution of higher twist PDFs is essential. The forth point of interest is the targetmass dependence of TMD distributions. At higher powers of small-b series, the target mass is the only scale that compensates the dimension of b n for twist-2 and twist-3 distributions. Thus, the here derived corrections are target-mass corrections ∼ (M 2 b 2 ). Their knowledge is essential since much of the experimental data is measured on nuclei.
Formally, the matching is obtained by operator product expansion of the transverse momentum dependent operator (defined explicitly in (3.1)) at small values of b. The later has the schematic form where z is the distance between fields of the operator along the light-cone. The operators O n (z) are light-cone operators with the collinear twist n. For example, the leading power operator is O 2 (z) ∼q(zn)[zn, 0]q(0), where n is the light-cone vector, [a, b] is the straight gauge link, and q is the quark field. Each operator O n is the integral convolution of a coefficient function and an actual quantum-field operator. The coefficient functions for O 2 are all known at next-to-leading order (NLO) in α s -expansion [16,17], and NNLO [5,[18][19][20]. The leading power coefficient function for unpolarized distribution has been recently computed at N 3 LO [21,22]. Beyond the leading power the information is sparse. The tree order matching for O 3 has been derived in [14], see also [10][11][12][13] for particular cases. The only NLO computation for O 3 is made for the Sivers function [23]. In this work, we derive only the tree order matching, ignoring the α s -suppressed terms in the coefficient functions. For that reason, we do not specify the renormalization scales and omit corresponding arguments.
The operators with the collinear twist n can be presented as a sum of operators with different geometrical twists, where we omit indices µ 1 ...µ n for brevity. For shortness, we call this procedure as twist-decomposition. Generally, operators O t depend on many spatial points, which are parameterized by a set of variables {y}. They are mapped to the single variable z by the integral convolution ⊗. The geometric twist has a strict definition as the "dimension minus spin" of the operator. Operators with different geometrical twists have different transformation properties and thus represent independent physical observables. The matrix elements of operators with a given geometrical twist define a self-contained set of PDFs. Such a set of PDFs does not mix with other sets, and their evolution is autonomous. For the introduction to the twist decomposition see e.g. [24,25], the review of modern development can be found in [9]. Therefore, the central task is to derive the twist-decomposition for each operator on the right-hand-side (RHS) of (1.1). In turn, the small-b expansion of TMDPDFs in terms of collinear PDFs is received by evaluating the matrix element over the derived operator relation.
In the case of FF, the twist-decomposition operation is not well defined. The main reason is the absence of local expansion for fragmentation operators. In ref. [26] it has been shown that OPE for FF is defined up to terms that satisfy the Laplace equation (for the twist-2 part). Therefore, alternative methods such as differential equations [26], Feynman diagram correspondences [18,[27][28][29] and Lorentz invariant relations [13], should be used. For that reason, we do not consider the matching of TMDFFs to FFs in the present work. For an interested reader, we present some discussion on power corrections for TMDFF in appendix B.
In this work, we compute only quark TMD distributions, since they are of the prime practical interest. In total, there are eight TMDPDFs at the leading term of the factorization theorem [30]. They can be split into two classes with respect to the structure of matching relations. Four distributions, namely f 1 (unpolarized), g 1L (helicity), h 1 (transversity) and h ⊥ 1T (pretzelocity) have contributions of only even collinear twists where T t is a collinear distribution of twist-t and T 2 (x) = f (x) is the twist-2 PDF. Another four distributions, namely f ⊥ 1T (Sivers), g 1T (worm-gear T), h ⊥ 1 (Boer-Mulders) and h ⊥ 1L (worm-gear L), have contributions of only odd collinear twists The graphical representation of these sums is shown in fig.1 In the present work, we derive the first and the second terms of this sum for all eight TMD distributions. In fig.1, the shaded areas show the corresponding terms.
To perform this computation, we use the technique inherited from [31,32], where it was used for the analyses of twist-4 operators. The technique is based on the local equivalence of the Lorentz transformation group to SL(2,C) group (so-called spinor formalism) [33]. Within the spinor formalism, the twist-decomposition can be elegantly formulated as the action of a certain spinor-differential operator (see sec.3.2). In ref. [9] this method has been used to derivate kinematic power corrections (t/Q 2 ) in deeply-virtual Compton scattering. In contrast to ref. [9], TMD operators are essentially non-local. They contain the gauge link along a staple contour. To overcome this difficulty, we introduce a formal local expansion for the TMD operator. To our best knowledge, it is the first time when the series of local operators successfully describes the infinite staple contours. Almost all expressions presented in this article are novel. Only some of them, namely b 0 -and b 1 -terms, can be found in literature and agree with it. Additionally, we demonstrate that the series of corrections for the unpolarized TMDPDF f 1 can be derived differently using the results of ref. [34] and agrees with it.
The article is organized as follows. In sec.2 we articulate all definitions used in our work. This is an important part because we deal with functions that do not have a common definition such as polarized distributions and PDFs of twist-3. In sec.2.3 we specify the conventions for the spinor formalism used in our work. Sec.3 is devoted to the detailed description of the computation method. It is split into three subsections in accordance with three principal stages of the computation:(1.1), (1.2) and (1.3,1.4). So, in sec.3.1 the expansion (1.1) of the TMD operator in the series of (local) collinear operators is described. In sec.3.2 we explain the method of the twist-decomposition and derive it for TMD operators (1.2), with results for twist-3 part presented in the appendix A. The particularities of the computation of matrix elements for twist-decomposed operators are given in sec.3.3. The result of the computation and the discussion are given in sec.4.

Definitions and conventions
In this work, we operate with eight TMD distributions, three collinear twist-2 distributions and four collinear twist-3 distributions. They essentially depend on the conventions related to P-odd structures, such as Levi-Civita tensor, γ 5 , etc. An additional set of conventions is brought by the spinor formalism which is used for the twist-decomposition. There is no commonly accepted convention for all these subjects in the literature, e.g. compare conventions in refs. [12-14, 24, 33, 35-37]. Nicely, the final result is mostly independent on the conventions, because it is the relation between physical distributions. Nonetheless, these conventions play an important role in the intermediate steps. In order to structure the presentation we collect all used definitions and conventions in this section.

Definition of TMD distributions
The light-cone decomposition plays the central role. It is defined by two light-like vectors n µ and n µ (n 2 =n 2 = 0, (nn) = 1). We use the ordinary notation of vector decomposition where v + = (nv), v − = (nv), and v T is the transverse component (v T n) = (v Tn ) = 0. In what follows, the directionn µ is associated with the large-component of the hadron momentum p µ , where M is the mass of the hadron (p 2 = M 2 ). It is important, that the hadron's momentum does not have a transverse component, which gives the physical definition of the transverse plane. The spin of the hadron is parameterized by the spin-vector S µ , (S 2 = −1, (pS) = 0). Its light-cone decomposition is where the λ is the helicity of the hadron, and s µ T is the transverse component of the spin, s 2 T = λ 2 −1. The generic quark TMDPDF is defined as where the vector b µ is a transverse vector, (bp) = 0. The Wilson lines in the definition (2.4) are straight Wilson lines. Rigorously, one should add the T-and anti-T-ordering within the TMD operator. However, for the parton distributions (in contrast to fragmentation functions) it can be safely omitted (see e.g. discussion in [23]). TMDPDFs that appear in different processes have Wilson lines pointing into different direction, which is indicated by ∓∞n in (2.4). So, the TMD distributions which appear in semi-inclusive deep-inelastic scattering (SIDIS) have Wilson lines pointing to +∞n, while in Drell-Yan they point to −∞n. In the following, we distinguish these cases, and the upper sign refers to the Drell-Yan case, whereas the lower sign refers the SIDIS case. The open indices (ij) of the TMD operator in eq. (2.4) are to be contracted with different gamma-matrices, which we denote generically as Γ, There are only three Dirac structures that appear in the leading term of the TMD factorization theorem, these are Γ = {γ + , γ + γ 5 , iσ α+ γ 5 }. Here, the index α is transverse and with 0123 = − 0123 = 1. In the naive parton model interpretation, these gamma-structures are related to the observation of unpolarized (γ + ), longitudinally polarized (γ + γ 5 ) and transversely polarized (iσ α+ T γ 5 ) quarks inside the hadron. The standard parameterization of leading TMDPDFs in the position space reads The tensors g µν T and µν T are defined as such that g 11 T = g 22 T = −1, 12 T = − 21 T = 1, and the rest components are zero. The definition of the TMDPDFs coincides with the conventional one in [14,35,38]. In the following, we also compare to refs. [9,32,36], where the definitions of and s T have opposite sign, and ref. [14], where the definition of and T has opposite sign (so, component-wise the tensor µν T is the same). TMDPDFs defined in (2.7-2.9) are dimensionless functions which depend only on the modulus of b, but not on the direction. The conventional names for them are (see e.g. [35,38]): unpolarized (f 1 ), Sivers (f ⊥ 1T ), helicity (g 1L ), worm-gear T (g 1T ), transversity (h 1 ), worm-gear L (h ⊥ 1L ), Boer-Mulders (h ⊥ 1 ) and pretzelosity (h ⊥ 1T ) distributions. The position space representation of TMD distribution is advantageous, because the TMD evolution is multiplicative in the position space. For that reason, the phenomenological studies that incorporate the TMD evolution are made in the position space, for the most recent examples see [7,8]. TMD distributions in the momentum space are obtained by the Fourier transformation The transformation rules for particular TMD distributions can be found in refs. [14,38].
There is no conventional definition of the collinear distributions of twist-3. Here, we use the definition used in [14]. We define where we omit Wilson lines that connect the fields in operators for brevity. The definition of the distributions T and ∆T coincides with [39] and [36], taking into account the difference in conventions for the -tensor (explained after (2.10)). In ref. [39] one can also find comparison with other definitions. The integration measure [dx] is defined as The delta-function in this measure reflects the independence of the matrix element on the global position (z 1 +z 2 +z 3 ) of the field operator. Due to the delta-function in the measure, the distributions of twist-3 effectively depends on only two variables. Nonetheless, it is convenient to keep all three variables x 1,2,3 as independent. It reveals the symmetric properties [14] T Each range of x 1,2,3 ≶ 0 has a specific partonic interpretation [24].

Spinor formalism
The twist-decomposition for local operators consists in the decomposition of tensors with many indices into irreducible representations of the Lorentz group SO (3,1). This procedure is greatly simplified in the spinor formalism. The spinor formalism is used in many parts of quantum field theory, for a review see [33,37]. Here we remind only of the essential properties and introduce conventions that are necessary for the current work.
The spinor formalism is grounded on the local isomorphism of the Lorentz group to the group of complex unimodular matrices SL(2,C). The isomorphism is realized by the map of a four-vector to a hermitian matrix by the rule The scalar product of any two vectors is x µ y µ = x αα yα α /2. In the spinor formulation one must distinguish dotted and undotted indices because they are related to conjugated representations, (u α ) * =ūα. The scalar product of two spinors is defined as where 12 = − 12 = 1 as in ref. [9,31,32]. The light-like vectors n andn in the spinor formalism can be written as where λ and µ are independent spinors normalized as (λμ)(µλ) = 2. The spinors λ and µ form the basis, which can be used to decompose any tensor. In particular, the decomposition (2.1) for an arbitrary four-vector is The Dirac bi-spinors are written as composition of two-component spinors The decomposition of these spinors in the basis (2.23) is where ψ + = (λψ), ψ − = (µψ), etc. In the similar way we write down the decomposition of the gluon-strength tensor, where f αβ andfαβ are symmetric tensors,f = f † . In our computation we face only the gluon strength-tensor with one index transverse, and another contracted with the vector n. The decomposition of such tensor is The first term in (2.28) corresponds to F −+ components, whereas the last two terms describe F µ+ with µ being transverse index. Using the definitions (2.25) and (2.26) we write down the decompositions for the bi-spinor combinations. They arē where the order of the fields on LHS and RHS is preserved. Let us note that only "plus" components of quark fields appear in these decompositions. It is not accidental, but is the part of definition for the leading power TMD distributions. The components ψ + , χ + , etc, are known as "good" components of quark field in contrast to "bad" components (ψ − , χ − , etc) [24]. The operators build only from good components (including also good components of the gluon field f ++ andf ++ ) are called quasi-partonic operators. Their geometrical twist coincides with the collinear twist. All operators of twist-2 and twist-3 can be expressed as quasi-partonic operators with the help of EOMs [34]. The last ingredient needed for our computation are the equations of motion (EOMs) for the quark field Dq = 0 (for massless quarks). In the spinor notation EOMs are Dα α ψ α = 0, χ α ← − D αα = 0, similar for other spinors. Contracting these equations with basis spinors one receives EOMs for particular components. For our purposes we need the following EOMs where D ab = D αα a α bα.

Twist-decomposition for TMD operators
In this section, we present details of the twist-decomposition procedure. The base of the method is taken from refs. [31,32], to which we refer for extended details and the theory foundation. There are three principal steps of the computation: 1. The operator is presented as a series of local operators.
2. Local operators are sorted by irreducible representation of Lorentz group (twists), and simplified using EOMs.
3. The series of operators with the same twists are summed back into the non-local form.
This is a rather traditional approach for twist-decomposition. Examples of such consideration for collinear operators can be found in refs. [24,31,32,34,40,41]. For TMD operators each of these steps has a certain particularity. Let us list these particularities, and explain the methods that were used to resolve them: 1. The TMD operator has infinite Wilson lines. Therefore, it is not possible to present it as the series of local operators directly. We regularize the TMD operator by truncation of the length of Wilson lines by the parameter L. Then TMD operators can be presented as a limit of triple series of local operators (3.11).
2. The local operators for TMDPDFs have three sets of indices (s, n, t). They correspond to the number of light-cone derivatives acting onq(s), the number of transverse derivatives (n) and the number of light-cone derivatives acting on q(t). Such structure somewhat complicates the twist-decomposition algebra, in comparison to the case of local operators that describe (Mellin moments of) collinear distributions, where all indices are alike. To simplify the twistdecomposition procedure we use the method introduced in refs. [31,32], which is based on the spinor formalism. This method allows to perform the twist-decomposition for the most general operators using only differential operations (compare to refs. [14,42] where off-light-cone generalizations of operators and integral equations are used, ref. [34] where the differential equation are used, refs. [40,41] where an explicit procedure of index symmetrization is performed).
3. The summation of the series is made for the matrix-elements, i.e. for distributions. It simplifies certain steps of computation, and helps to resolve potential ambiguities in the limit L → ∞.
The following sections give details for each step in the order. The results of the computation are collected in the sec.4. We stress that such an approach is suitable only for TMDPDFs, but does not apply for TMDFF. The discussion for TMDFFs case is given in sec.B.

TMD operator as a series of local operators
Let us introduce the TMD operator in the form The upper(lower) sign corresponds to Drell-Yan(SIDIS) induced TMD. The transverse gauge link [∓∞n + b, ∓∞n] ensures the explicit gauge invariance of the operator. At the tree-order quantum fields can be treated as classical fields, and thus the small-b expansion is an ordinary Taylor expansion. It is convenient to write it in the form (1.1) The operators on RHS of (3.2) have collinear twist n + 2, which follows from their dimension. At the same time, these operators do not have a definite geometrical twist, and therefore, their matrix element is a complicated composition of collinear distributions with different properties. Our goal is to perform the twist-decomposition and express these operators in terms of operators with definite geometrical twist.
In contrast to collinear operators, the TMD operator (3.1) spans the infinite range along the light cone. It is the most famous feature of TMD operators, and it leads to many physical effects, such as rapidity divergences [43], double-scale nature of TMD evolution [39] and the sign-change of P-odd distributions [44]. In the b → 0 limit, the infinite Wilson lines are partially compensated, due to the transitivity property of Wilson lines. For example, at n = 0 the operator (3.3) simplifies to However, already at n = 1 the infinities enter the expressions, is the gluon-strength tensor. The operators on the RHS of these formulas are ordinary collinear operators. The infinite size of the TMD operator is reflected in the limit of integration in the last term of (3.6). In the form like (3.6), the twist-decomposition of operators is straightforward, although algebraically complicated [14]. The main reason of the complication is that each next order of expansion introduces new classes of operators, for instance at n = 2 the operators likeq(zn)F +µ (τ 1 n)F +ν (τ 2 n)q(0) andq(zn)[D ν F µ+ (τ n)]q(0) appear. Each of these new classes requires individual investigation, and therefore, this approach is ineffective.
To avoid these complications we operate directly with the operators (3.2) with the help of the following formal procedure. We introduce the regularized operator, This operator turns to the TMD operator in the limit L → ∓∞. Note, that the same regularization also regularizes rapidity divergences and can be used used to derive the non-perturbative definition of the Collins-Soper kernel [45]. In the regularized form the operator (3.3) is This operator can be written as the formal expansion, where In this way the TMD operator (3.1) is presented as a triple sum The RHS of this expression is suitable for the twist-decomposition procedure. The limit L → ∓∞ must be taken with caution. In particular, the summation over s and t must be done before the limiting operation, and the result of the summation should be presented in the form that does not allow any ambiguity. The significant simplification of the summation procedure comes from the possibility to neglect the total derivative terms. The matrix element of a total-derivative operator is proportional to the momentum transfer, (3.13) Consequently, the total derivative operators do not contribute to TMDPDFs, and we neglect such terms in the following. After elimination of the total derivative contribution the result of summation is independent on L in many cases. The simplest example is n = 0 case, where ∞ s,t=0 with z 1 − z 2 = z being independent on L. Only in the cases of Sivers and Boer-Mulders functions the limit L → ∓∞ is not so straightforward. Let us also note that the direction of the limit is the only difference between Drell-Yan (L → −∞) and SIDIS (L → +∞) cases in the representation (3.11).

Twist decomposition in the spinor formalism
The twist-decomposition for the local operators (3.12) is equivalent to the decomposition into irreducible representation of Lorentz group. The highest spin representation (completely symmetric and traceless in all vector indices) corresponds to the twist-2 term. The next representation (anti-symmetric for one pair of indices, and symmetric and traceless for all other vector indices) corresponds to twist-3 term. Despite the twist-decomposition is a straightforward operation, it is algebraically complicated, especially for the operators with many transverse indices (for which one has to subtract traces). Additional complication is caused by EOMs, which could relate operators. In refs. [31,32] it was observed that the decomposition of higher-indices tensors over irreducible representations is simpler in the spinor representation. The main simplification comes from the anti-symmetry of the scalar product (2.22). Due to it, any symmetric tensor (in the spinor space) is automatically traceless. The irreducible representations in SL(2,C) group are obtained by simple symmetrization or anti-symmetrization of spinor indices. This operation can be presented as a differential operator, what essentially simplifies the algebra. Here we present this methods in a slightly modified form, which makes its application more explicit. Let us consider an operator where all open indices are contracted with basis spinors: where Λ 1,2 are monomials of basis spinors, for example (3.27). The lowest geometrical twist contribution can be obtained by the symmetrization of dotted and undotted indices (independently). It can be done for the operator, or for the tensor Λ. Due to the irreducibility of the representation the symmetry properties of one are mapped to another in the convolution. The symmetrization of the tensor Λ can be made by the following operation where and n is the number of spinors µ in the tensor Λ. The logic behind the construction is the following. First, the action of derivatives (λ∂)/(∂µ) replaces all µ's by λ's. The resulting tensor is automatically symmetric in indices. Next, the derivatives (µ∂)/(∂λ) replace entries of λ by µ in the fully symmetric fashion. The higher twist representations are build by anti-symmetrizing pairs of indices. So, the antisymmetrization of a single pair can be done by the operator where (µ∂)/(∂µ) and (λ∂)/(∂λ) defined similarly to (3.17). The symmetrization of n pairs is done by The operators S and A commute, [A n , S k ] = 0, for any n and k. The complete decomposition of a tensor with n entries of spinor µ then reads where c n,k are numbers depending on the tensor Λ. Each term in this sum corresponds to an irreducible representation, and the operators S n−k A k are projectors onto this representation.
To find the coefficients c n,k of the decompositions (3.20) we need to normalize the operators A k S n , such that (A k S n ) 2 = A k S n . For example, for the symmetrization operator we compute The operators (λ∂)/(∂λ) and (µ∂)/(∂µ) count the number of λ's and µ's in the tensor. Therefore, S n is indeed the projector to the totally symmetric representation, and the normalization factor for S n is (N − n)!/(n!N !) where N is the total number of indices of the tensor. In the same way, one can check that A n S k are projectors, and find the corresponding normalization factors. For our computation we only need the first two terms of the expansion (3.20). They are 22) where N is the total number of indices, and n is the number of µ's in the tensor Λ. Using this construction we build the operators that extract a part with the certain geometrical twist from the operator (3.15) (3.23) The projectors to the twist-2 iŝ where n, k, N − n andN − k are the numbers of µ,μ, λ andλ in the operator. The operator that projects the twist-3 has two termsT is obtained by the interchange of barred and un-barred variables.
In this way, the twist-decomposition is reduced to purely algebraic manipulations with monomials. One should distinguish the chiral-even (given in (2.29)) and chiral-odd (given in (2.30)) compositions of spinors, because they have different number of barred and un-barred spinors. Apart of this, the evaluation of all cases is alike. In the following we demonstrate the results of action by T 2 and T 3 on O Γ s,n,t . For simplicity of presentation, we replace the indication of the Dirac structure in O Γ s,n,t , by the indication of the corresponding spinor combination. Additionally, we write operators using spinor notation only. For example where the prefactor is originated from the conventions (2.22-2.24). Also we eliminate all totalderivative terms without indication.
The derivation of the twist-3 part is slightly more cumbersome, due to the reduction procedure to the quasi-partonic form. We remind that the quasi-partonic operators consist only of "good" components of fields and D λλ . All twist-3 operators can be reduced to this form. The reduction procedure is done as follows. After the action of A 1 we obtain operators of the formψ + D N λλ D λμ D M λλ ψ + andψ − D N λλ ψ + . Next, we commute the derivative with the transverse index to the quark field, such that it can be replaced by appropriate EOM (2.31). After this procedure all "bad" components of quark field cancel, and the twist-3 operator has the quasi-partonic form. The expressions for the twist-3 parts of operators are simple but rather lengthy. For example,

T (μλ) 3
Oψ + ψ+ s,n,t = 2ig n k=1 µ∂ ∂λ Here, the summations over m are originated from the commutation procedure. Other spinor combinations and action ofT (µλ) 3 differ from this example by ±1 terms in the factorial factors and summation limits. For completeness these expressions are listed in the appendix A.

Assembling the final result
The last step of the computation is to sum the operators over s and t to the non-local form and perform the limit L → ∞. This procedure can be done directly in terms of operators. However, the resulting expression has an overcomplicated form, because the most part of the expression vanishes on the level of matrix element due to symmetry relations (2.19,2.20). To avoid these complications we first compute the matrix element and then sum the expression. Conveniently, the evaluation of the matrix elements can be done before we apply the derivatives (µ∂)/(∂λ) and (μ∂)/(∂λ), which then act on the parameterization of the matrix element. The process is identical for all distributions. Here we demonstrate the computation for the case Γ = γ + which contains unpolarized and Sivers TMD distributions. The peculiarities of the computation for other cases are discussed at the end of the section.
The matrix element for the twist-2 operator with Γ = γ + is expressed in the terms of unpolarized PDF (2.12) as (2.29) Substituting it into (3.28), we observe that the derivatives (µ∂)/(∂λ) and (μ∂)/(∂λ) now could act only on p λλ . At the same time, the vector p µ does not have a transverse part (by definition), i.e. p λμ = p µλ = 0, and thus only the terms with n = 2k are non-zero: Now, we can pass to the vector notation, and the expression for the matrix element takes the form where we used that bb = −b 2 /2, p λλ = 2p + and p µμ = 2p − = M 2 /p + . Presenting the factorial factor by the beta-function (for n > 0) we evaluate the sum over s and t. The summation produces the exponent exp(i(z 1 − z 2 )zp + ) that is independent on L due to (3.10). Therefore, in this case the limit L → ±∞ is trivial. The summed expression reads Finally, we make the inverse Fourier transformation and compare the result to the parameterization (2.7). We find that the twist-2 part of the operator contributes only to the unpolarized TMD distribution. The final expression is convenient to present in the form of Mellin convolution (4.12) tw4). (3.35) This is the complete part of the small-b expansion for f 1 (x, b) that is matches twist-2 PDFs. The corrections to this expression contain PDFs of higher twist, starting from twist-4 (the absence of twist-3 part is demonstrated later). Also each term in this sum receives the perturbative corrections. All these corrections are indicated by O(α s , tw4). The computation of the twist-3 part follows the same pattern, but involves more terms. We use the matrix elements (2.15,2.16), and present them in the form of double moments and similarly for other spinor combinations. Now, the derivatives (µ∂)/(∂λ) and (μ∂)/(∂λ) also act on the vector S, and thus the expression analogous to (3.32) contains two terms. One with n = 2k − 1 that is proportional to transverse part of S, and another one with n = 2k − 2 that is proportional to S µμ = −λM/p + . The parts proportional to λ add up to zero for Γ = γ + , however, in the case of Γ = γ + γ 5 they produce the twist-3 terms in the helicity TMD distribution g 1L . After these operations the summation over m reduces to geometric progressions. The resulting expression is rather lengthy where we have used that x 1 + x 2 + x 3 = 0 for simplifications. This expression is very representative and demonstrates many features of the computation. Notably, there are two types of contributions relative to the summation over s and t. The regular ones that are in the third and forth lines, and the irregular one that is in the last line.
The regular contributions have a general form of (see also (3.14,3.33)) ∞ s,t=0 (3.38) and are explicitly independent on L, since z 1 − z 2 = z (3.10). For regular contributions, the limit L → ∓∞ is trivial, and the resulting operators are spatially compact. The final expression has ordinary form with twist-3 distributions, see eqns. (4.3,4.4,4.5,4.7,4.8). In the current example with Γ = γ + , these term do not appear in the final expression, due to symmetry properties (2.19). Indeed, the change of variables {x 1 , x 2 , x 3 } → {−x 3 , −x 2 , −x 1 } flips the common sign of the prefactors in square brackets, but leaves T /x 2 2 and ∆T /x 2 unchanged. Thus, the contribution of the third and the fourth lines vanishes.
The evaluation of irregular terms requires special attention. The sums over s and t lead to the following expression To resolve this ambiguity, we use the symmetry of the integration measure [dx] and transfer the variable L to the limits of the integration It is important that the sign of this expression depends on the sign of the limiting direction. This is how the famous sign-flip for the Sivers function appears in this computation. To integrate over x 2 , we extract all entries of x 2 in (3.40) with the help of relation Clearly, only the term with m = 0 contributes. Performing the Fourier transformation and comparing it to (2.7) we conclude that it is the contribution to the Sivers function, which reads This is the matching of the Sivers function at all orders of b 2 -expansion to collinear distributions of twist-2 (which is null) and twist-3. Miraculously, only the Qui-Sterman function T (−x, 0, x) [46] contributes to this expression. As expected, the sign of the Sivers function depends on the direction of the gauge contour.
The computation of other Lorentz structures is done in the same way. The differences to the demonstrated computation are immaterial. The resulting expression are not as simple as for Γ = γ + case. In particular, generally the distributions have entries of both twist-2 and twist-3 PDFs. It happens due to non-vanishing contributions with derivatives of the vector S µ . For example, for Γ = γ + γ 5 the analog of (3.31) has S λλ which after the differentiation (3.32) produces the transverse components of S. These terms result into twist-2 PDFs within worm-gear function. In all P-even cases only regular terms contribute, whereas P-odd cases (i.e. Sivers and Boer-Mulders functions) are given by only irregular terms. The results of our computation together with additional comments are presented in the next section.

Results and discussion
Following the method derived in the previous section we routinely computed all (quark) TMD distributions. Here is the list of expressions. The chiral-even distributions are The chiral-odd distributions are In these expression u = 1 − u, and where Φ is the Lerch function. The twist-3 PDFs are conveniently gathered into the following combinations where we use the shorthand notation ( The point u = 1 is regular for all integrals (4.1-4.8). The parameter of the expansion x 2 M 2 b 2 /4 is negative, due to b 2 < 0. In the following we discuss particularities of each TMD.
Unpolarized TMD f 1 . The unpolarized TMD f 1 (x, b) is given in (4.1). The leading term of this expression is well-known. Moreover, the perturbative corrections for the leading term are know up to α 3 s -order [21,22]. Peculiarly, the contributions proportional to the twist-3 are absent at all powers. So, the next contribution to f 1 contains twist-4 PDFs.
The series (4.1) can be summed to the Bessel function with the result where |b| = √ −b 2 . However, the value of the expansion parameter (xM |b|/2) 2 is so small that the numerical difference between the summed series and the first term is vanishing (for non-extreme values of b). Moreover, since the cross-section for TMD observables is dominated by small values of b (b < 1GeV −1 ), the derived corrections are negligible for many experiments (e.g. for the LHC, where typically x ∼ 10 −2 ). The comparison of (4.13) truncated series (4.1) is shown in fig.2. Fig.2 also shows the unpolarized TMDPDF f 1 (x, b) extracted in ref. [7]. The difference between theoretical and extracted curves demonstrates that the target-mass corrections give negligibly small portion of power corrections. The larger portion of corrections is given by the twist-4 (and higher) distributions. In ref. [47] it has been demonstrated that infrared renormalon poles for unpolarized TMD f 1 are proportional to (xM 2 b 2 ) n , and therefore, higher-twist contributions must be enhanced by at least a power of x −n .
The unpolarized distribution is the only case where the expression for target-mass corrections can be compared to the literature. To make the comparison, we mention that the twist-2 operators (at the tree order) are not affected by the direction of the Wilson lines. Therefore, to obtain the In all cases, the PDF is taken in the convolution with NNLO perturbative coefficient function. The curve labeled as "sum" is the summed expression for target-mass corrections (4.13). The curve with a band is the extraction of non-perturbative TMD distribution f1 made in ref. [7].
matching to the twist-2 operators (at the tree order), one could ignore the staple contour of the TMD operator and connect the quark-fields by the straight gauge link. This operator is the usual QCD string-operator, and the small-b expansion is the ordinary x 2 -expansion used in DIS. The twist-2 part of such operator at all powers of x 2 can be found in eqn.(5.1) of ref. [34], where it has been derived by means of differential equations. The expression in [34] is given in the coordinate space, and after the Fourier transformation coincides with (4.1). It gives a non-trivial check of the computation method and results. Sivers f ⊥ 1T and Boer-Mulders h ⊥ 1 functions. Sivers and Boer-Mulders functions are given in (4.2) and (4.6) correspondingly. Both functions are P-odd, and thus they have a different sign for Drell-Yan and SIDIS configurations [44]. In our computation, the sign-flip arises due to the order of integration limits in the irregular terms (3.40). These are the only terms where the infinite size of the Wilson lines (the parameter L) plays a role. The leading power matching for the Sivers function is known for a long time [11,12], whereas the leading power matching Boer-Mulders function was derived in [14]. Our computation agrees with these references.
In the limit L → ∓∞ the coefficient function reduces to δ(x 2 ) and thus Sivers and Boer-Mulders functions are expressed through T (−x, 0, x) and δT (−x, 0, x) at all powers of b 2 . The function T (−x, 0, x) is known as a Qiu-Sterman (QS) function [46], and δT (−x, 0, x) is the analog of QS function for the chiral-odd operator. Both TMDPDFs are expressed via QS functions at all orders of b-expansion. This observation is a non-trivial fact because the QS functions are not autonomous, and mix with the bulk of twist-3 distribution during the evolution [36]. The NLO expression for the leading power coefficient function [23] contains the non-QS but only in the logarithmic terms (i.e., the terms responsible for the evolution effects), whereas the finite part is proportional to the QS function. Based on this information, we make a conjecture that Sivers and Boer-Mulders functions can be expressed entirely through the QS function, up to logarithmic terms.
Helicity g 1L and transversity h 1 TMDs. The helicity and tranversity TMDPDFs are given in (4.3) and (4.5) correspondingly. These distributions have leading power matching to the twist-2 distributions. The perturbative coefficients for the leading power matching are known up to NLO [16,17] (for helicity) and NNLO [19] (for transversity). In contrast, to the unpolarized TMD f 1 (4.1) distribution g 1L and h 1 have non-zero contribution of twist-3 PDFs. The expressions given in this work do not account the twist-4 PDFs, and thus already the n = 1 term is incomplete and contains additional twist-4 contribution (see also fig.1).
Worm-gear TMDs g 1T and h ⊥ 1L . The worm-gear distributions g 1T and h ⊥ 1L are given in (4.4) and (4.7), correspondingly. They have the generic structure of distributions with the leading power matching at collinear twist-3 operator (see fig.1). Both worm-gear TMDPDFs have the factor x, which suppresses them at small-x. It is interesting to mention that despite both distribution have similar form, TMD g 1T is expected to be much larger according to a recent measurement [48]. The leading power expressions for both distributions were derived in [13] (using Lorentz invariance relations) and [14] (using the off-light-cone parametrization), and agree with our computation.
Pretzelosity TMD h ⊥ 1T . The expression for the pretzelosity TMD is given in (4.8). The expression (4.8) demonstrates a non-trivial value for pretzelosity distribution for the first time. The leading power term for the pretzelosity reads This expression is incomplete, because at the same level of accuracy PDFs of twist-4 could contribute. Thus eqn.(4.14) is only the Wandzura-Wilchek approximation to the full expression. Note, that unless T h (x) has some anomalously strong behavior, (that could happen only at x → 0, since T h (1) = 0) the pretzelosity distribution is a very small function. At x → 1 the convolution is suppressed by the factor (1 − u 2 ), whereas at x → 0 there is a common factor x 2 . This observation is in a general agreement with the experimental data [48,49], which is compatible with zero for the pretzelosity-involving observables.
There are no contribution of twist-2 PDFs. The absence of matching to twist-2 PDFs is true for all orders of the perturbation theory that were checked in ref. [19] up to α 2 s -order, and discussed in ref. [50]. This fact can be proven comparing the following combinations of chiral-odd TMDPDFs The non-zero pretzelosity appears only if the matrix element Φ [iσ α+ γ5] has terms with different parity (so the expressions in (4.15,4.16) in brackets are different). At the twist-2 level there is only one PDF, and thus the pretzelosity is null. At the twist-3 level there are pairs of operators with different parity: e.g.ψ + f ++χ+ andψ +f++χ+ which produce asymmetry in (4.15,4.16), because p, S|ψ +f++χ+ |p, S = 0. Since QCD Lagrangian conserves parity, these statements are generally preserved at any order of perturbative expansion in a properly defined regularization scheme. In the dimension regularization the Levi-Civita tensor is not uniquely defined and thus the symmetry between relations (4.15,4.16) could be violated. In this case, one observes the non-zero pretzelosity [16,19] in the -suppressed terms. However, it is only an artifact of the dimensional regularization, and it must vanish at → 0. The same statement holds for the quasi-TMDs for which the non-trivial but -suppressed matching to pretzelocity has been observed in ref. [51] at α s -order.

Conclusion
TMD distributions are related to the collinear distributions in the limit of small-b. In the present work, we have studied this relation at all powers of b-expansion and derived the contributions with twist-2 and twist-3 PDFs. Our computation includes all eight polarized TMDs. From the perspective of the resummation approach, the computed corrections are target-mass corrections. The summary of the here-derived and known results is presented in table 1. It is the first study of the matching between TMDPDF and collinear PDFs beyond the leading power, to our best knowledge. Thus, many results and conclusions made in this paper are novel.
The list of expression for TMDPDFs is presented in sec.4. All expressions have the following common structure where F is a TMD distribution, f is a (combination of) collinear distribution, and ⊗ is the Mellin convolution. In addition to the power-suppression each term is accompanied by the factor x 2n . Due to its numerical value of computed correction is almost negligible. The appearance of targetmass corrections in the combination (x 2 M 2 ) justifies the usage of proton TMDPDF for nuclei. For nuclear observables the target-mass is Z-times larger, but the measured x is Z-times smaller (with Z being the atomic number). Together these factors compensate each other, and thus the nuclear TMDPDF is roughly a nucleon TMDPDF. The knowledge of target-mass corrections is also essential for lattice computations of polarized quasi-TMDPDFs [51,52]. One of the central results of this work is the elaboration of the twist-decomposition-method. The method is taken from [32,36] and uses the simplifications of the spinor formalism. In the present study, it is applied at the tree-order, but similarly it can be used together with loop-calculation. Let us note that the computation is made in the position space, which is the only straightforward way to compute such power corrections. It is clear that each term with n > 0 of the expansion (5.1) is power-divergent in the the momentum space. Moreover, the resummed series of power corrections (4.13) also has a divergent Fourier transform. It demonstrates the well-known fact that at certain b, the perturbative series (in any form) fails to describe non-local objects, and truly non-perturbative effects must be accounted.
For the first time, we derive the non-zero matching of pretzelosity distribution to the collinear distributions. Its leading power expression contains twist-3 (given in (4.14)) and twist-4 PDFs (not computed). Previously there were unsuccessful attempts to find a non-zero matching of pretzelosity to twist-2 PDFs [16,19]. In sec.4 we provide the argumentation of why the pretzelosity does not have a contribution of twist-2 PDFs at all orders of power and α s corrections.

A Twist-3 part of the operators O s,n,t
In this appendix, we collect the expressions for twist-3 parts of the operators O s,n,t . The definition of operatorT 3 is given in (3.25,3.26). The definition of the operator O s,n,t is given in (3.12,3.27). The route of derivation is explained in sec.3.2.

B Power corrections for fragmentation functions
In the case of TMDFF the matching to the collinear FF is not entirely defined. The reason is the absence of local-operator expansion for FF-type operators. Only indirect methods of twistdecomposition for FF operators are possible, such as differential equations [26], Feynman diagram correspondences [18,[27][28][29] and Lorentz invariant relations [13]. However, all these methods are ambiguous, and allow an addition boundary contributions. In this appendix we demonstrate the result which appears if we partially ignores these problems. For detailed discussion we refer to [26].
The quark TMDFF is defined as [1] ∆ where we use x for collinear fraction of momentum instead of traditional z to avoid the clash of notation. The Wilson lines are connected to −∞ for SIDIS-like process, and to +∞ for e + e −annihilation-like processes. For spinless particles there are only two TMDFFs that contribute to the leading term of the factorization theorem. They read These functions are called unpolarized (D 1 ) and Collins (H ⊥ 1 ) functions. In contrast to the TMD-PDF operator, the TMDFF operator cannot be presented as a single T-ordered operator. This is the central issue because it prevents application of OPE. The definition of collinear FF of twist-2, Note, that collinear FF is defined with relative factor x 2 in comparison to TMDFF [1].
In the case of FF-type operators the procedure described in sec.3 should be modified. The reason is the operator in (B.1) has two parts which could not be joined under a single T-order sign. So, one must distinguish causal (indicated by (+)-sign) and anti-causal (indicated by (−)-sign) fields which then interact by a cut-propagator. Altogether it is known as Keldysh technique. In Keldysh technique the TMDFF-operator can be written analogously to ( ν ]. These terms cannot be written in a "quasi-partonic" form. It is an unsolved issue. The next principal difficulty appears when we combine operators to the non-local form. We have not found a way to perform this procedure on the operator level such that the limit L → ±∞ can be taken. Alternatively, we can use the analog of (3.31) Tr color N c X 0|γ + − → D N + q|h(p, s) + X h(p, s) + X|q Let us note that this is not a very well defined expression, because it assumes that LHS is dependent on N + M only, which generally does not hold. Nonetheless, operating in this way we receive the all-order expression for unpolarized TMDFF (compare to (3.31)) The most notable part here is the "improper" range of the integration over u. The ambiguity in the definition of OPE for FF, allows us to add any function that satisfies Laplace equation. For the detailed discussion we refer to sec.3 and 4 of ref. [26], where it is shown that for the unpolarized case a convenient addendum is expression (B.7) with integration range for u extended to (0, ∞). Subtracting this expression we get