Transverse momentum dependent transversely polarized distributions at next-to-next-to-leading-order

We calculate the matching of the transversity and pretzelosity transverse momentum dependent distributions (TMD) on transversity collinear distribution at the next-to-next-to-leading order (NNLO). We find that the matching coefficient for pretzelosity distribution is zero, despite the matrix element for it is nontrivial. This result suggests that the pretzelosity matches a twist-4 distribution. The matching for transversity TMD distributions is provided for both parton distribution functions and fragmentation functions cases.


Introduction
The transverse momentum dependent distributions (TMD) offer the possibility to study the structure of hadrons in great detail [1]. The TMD parton distribution functions (TMDPDF) and fragmentation functions (TMDFF) enter factorization theorems for cross sections of such processes as Drell-Yan or semi-inclusive deep inelastic scattering (SIDIS) [2][3][4][5]. The usage of the highest available perturbative input is important for the successful description of the experimental data, and significantly increases the predictive power of the framework [6,7]. Recently, many efforts were made to evaluate elements of TMD factorization theorem at NNLO. As a consequence, the evolution of TMD distributions has been calculated at two and three loops [8][9][10] and the matching of unpolarized TMDPDFs and TMDFFs have been calculated up to two loops respectively in [11][12][13][14][15][16] and in [16,17]. The status of the polarized distributions is weaker, given also a large number of different distributions. In this work, we consider transversity and pretzelosity (also called quadrupole) TMD distribution, and we evaluate their twist-2 matching at two loop order.
Both these distributions have been recently subject of experimental, phenomenological and theoretical investigations. The SIDIS data relevant for this extractions come mainly from HERMES [18] and COMPASS [19,20]. Recently also RHIC has provided data in this direction [21] and we expect that transversity will be one of the central measurements in future EIC and LHCSpin. The transverse momentum dependent transversity has been extracted using SIDIS data by Anselmino et al. in [22][23][24] with Gaussian models without taking into account the TMD evolution. In refs. [25,26] the lowest order evolution has been considered. An issue of these extractions is the size of the theoretical error. The reduction of it essentially requires the inclusion of the higher order perturbative information. The detailed discussion on theoretical errors, their dependence on perturbative order and related issues has been recently produced in [7]. Let us also mention here the attempts to constrain the transversity distributions by Monte Carlo and lattice collaborations [27]. More work in this sense is expected in the future. The transversity TMDFF is also an interesting and practically important object, see [28] for a recent review. For what concern pretzelosity, we mention here the recent analysis made in [29,30]. According to this analysis, the pretzelosity distribution is very small and practically consistent with a null value.
The one loop results for twist-2 TMDPDFs matching have been obtained recently by our group in [31]. The transversity and pretzelosity belong to the set of TMD distributions which in the regime of large transverse momentum or, equivalently, the small transverse distance match on twist-2 functions. This fact allows us to perform the two-loop calculations in a manner similar to [16]. Namely, we start from the operator definition of transversely polarized TMD distribution and consider its quark matrix element perturbatively with successive matching on the collinear matrix element. The main parts of the calculation practically coincide with the unpolarized case, apart from the algebraic structures and several new master integrals. The calculation is done in the δ-regulator scheme introduced in [8]. The structure of ultraviolet and rapidity divergences is independent of polarization, as it is predicted by TMD factorization theorem [2,5], and confirmed by the present calculation. Therefore, the renormalization of these divergences is done using the universal TMD soft factor [8] and renormalization constants calculated in the unpolarized case [16].
We anticipate here that the matching coefficient obtained for pretzelosity distribution is zero both at one and two loops. This result is particularly surprising since in principle this observable is expected to have a twist-2 contribution [31]. As a counterexample, we recall that the linearly polarized gluon TMD, which has the same quadrupole tensor structure has a matching coefficient different from zero already at one loop order [31]. So, the theoretical result that we obtain is puzzling. Nonetheless, it confirms the phenomenological and experimental analyses done in [29,30].
For the case of transversity, we provide the matching results at next-to-next-to-leading order (NNLO) both for TMDPDF and TMDFF cases. In the case of TMDFF, we also report the NLO expression, which is missing in the literature to our knowledge. Thus, with the result of this work, the transversity distribution is evaluated at the same level of precision as the unpolarized distributions, and it is the first example of NNLO evaluation of polarized distribution in the TMD factorization formalism. Consequently, the phenomenology for related observables can be developed with a similar level of precision. We provide the essential notation in sec. 2 and more technical details of TMD distributions in sec. 3. The results for transversity and pretzelosity TMDPDFs are described in sec. 4 and 5 and the case of TMDFF is considered in sec. 6.

Transversely polarized TMD distributions
The TMD distribution of the transversely polarized quark is where index α is transverse and n is a light-like vector. The gauge linksW T n (x) are rooted at the position x and continue to the infinity along the direction n. We use the standard notation for the light-cone components of vector v µ = n µ v − +n µ v + + g µν T v ν (with n 2 =n 2 = 0, n ·n = 1). The superscript T in Wilson lines indicates that they are incorporated by an additional transverse link at light-cone infinity, which ensures the gauge invariance of the operator.
The transversely polarized TMD distribution is parameterized in terms of four TMD parton distribution functions (TMDPDFs), which were originally introduced in momentum space [32,33]. For our purposes we need the equivalent parametrization in the position space. It reads where s T is the transverse part of the hadron spin, λ is helicity, M is the mass of hadron and b 2 = −b 2 > 0. The detailed relation between momentum and position space definitions can be found in [34,35]. In this work we consider the functions h 1 , which is known as TMD transversity distribution, and in the function h ⊥ 1T , which is known as the pretzelosity distribution. The small-b operator product expansion (OPE) allows the systematic expansion of the TMD operator in powers of b. The operators that are associated to the powers of b are then classified by twists. The evaluation of the matrix element of the small-b OPE results into an expression of the form where h are collinear distributions, C are coefficient functions and ⊗ is the integral convolution in the momentum fractions. The terms in eq. (2.3) incorporate all tensor structures of the TMD distribution parametrization in eq. (2.2). Extracting particular tensors, one can find the matching of individual TMDPDFs onto collinear functions. In particular, the tensor structure of h 1 and h ⊥

1T
appears in the twist-2 term, while the tensor structure for h ⊥ 1L and h ⊥ 1 can be produced only at the twist-3 [35]. In this work we concentrate on the twist-2 contributions.
The OPE in transversely polarized case has exceptionally simple structure, since there are no gluon operators. The only PDF that contribute to the matching of these twist-2 distributions is the collinear transversity PDF. Its expression reads [36] The PDF h 1 (x) can also be interpreted as the probability distribution to find a transversely polarized quark in a hadron. The coefficient functions of the OPE are dimensionless and the dependence on b enters only via logarithms, or via dimensionless tensors. Generally, the twist-2 coefficient functions can have structures ∼ g α,β T and ∼ b α b β /b 2 . It is natural to decompose it as here the parameter of dimension regularization (d = 4 − 2 ), and The pieces of this decomposition do not mix. In particular, the tensor in the second term of eq. (2.5) has zero trace (in d = 4 − 2 dimensions). Comparing the parametrization in eq. (2.5) with the parametrization for TMD distributions we find that the matching for individual TMDPDFs are where we explicitly express the Mellin convolution integral. In these formulas we also suppress the scale dependence of functions which is discussed in details in the next section. The sum over flavors runs only over non-singlet combinations, since there are no gluon operator with transverse polarized configuration that could mix with the quarks. The coefficient functions δC and δ ⊥ C can be evaluated perturbatively. Naturally, at tree order ∼ a 0 s , the coefficient function in eq. (2.5) is proportional to g αβ T , and thus only δC is non-zero. The terms proportional to ∼ b α b β /b 2 are generated by loop-diagrams and appear already at oneloop level [31]. In ref. [31] it has been found that the 1-loop contribution to δ ⊥ C is zero in four dimensions, but it can be different from zero in dimensional regularization at order O( ). This observation suggests that potentially this contribution does not vanish at two-loop level.

Evaluation of small-b OPE
In this section we discuss technical details of the matching procedure at the twist-2 level. In particular, we present the properties of the operator under renormalization and discuss the divergences that appear along the calculation. Since these details are universal for all TMDPDFs, we do not specify them, but instead we use the generic notation Φ, that in the context of this work mean h 1 or h ⊥ 1T .

Renormalization of TMDPDF
The TMD operator has two types of divergences. The ultraviolet (UV) and rapidity divergences. Both these divergences are renormalized by appropriate renormalization constants [5]. Consequently, the renormalized, and hence physical, TMDPDF depends on two scales. Traditionally, UV renormalization scale is denoted by µ and the rapidity renormalization scale is denoted by ζ. The renormalized TMDPDF has the form where we explicitly show the dependence on regularization parameters. In particular, is the parameter of dimensional regularization (d = 4 − 2 ) and regularizes UV divergences, which are renormalized by the factor Z. The δ is the parameter of δ-regularization [8,16], and R is the rapidity renormalization factor. The singularities in and δ cancels in the product of eq. (3.1). The renormalization factors are independent of the Lorentz structure but they change according the parton color representation. The renormalization factors Z and R are intertwined and it is important to specify in which order the divergences are subtracted. Here we work in the scheme where renormalization of rapidity divergences is made prior to the renormalization of UV divergences. The final result for Φ ren is of course independent of the subtraction order. The renormalization factors are scheme dependent. For the UV renormalization we use the MSscheme. To specify the renormalization scheme for rapidity divergences we recall that within the TMD factorization theorem the rapidity divergences are compensated by the soft factor [2,5,37]. The soft factor is defined as following where S n andSn stand for soft Wilson lines along n andn (for the precise definition of W T and S T see e.g. [16]). The factors R introduced in eq. (3.1) also renormalize the soft factor, such that the whole factorization expression is finite, see the proof and detailed derivation in [5]. It can be shown that within a properly defined scheme the renormalized soft factor is trivial, i.e.
where δ ± regularize rapidity divergences in the corresponding direction. In this scheme the rapidity renormalization factor has an exceptionally simple form [3,5,8] Such a scheme is very natural since it does not leave any remnant of soft factor in the factorization theorem, and therefore, coincides with other popular schemes of rapidity renormalization, e.g. with the one suggested in [2]. Although here we adopt the δ-regularization notation, these expressions could be translated to other regularization schemes, e.g. the translation dictionary of δ-regularization to the regularization by tilted Wilson lines is given in [38].

Scaling properties
The dependence on renormalization scales is given by the pair of evolution equations The anomalous dimensions are defined via the corresponding renormalization constants. A detailed study of this system has been recently presented in ref. [7]. The values of anomalous dimensions are known up to three-loop order inclusively [9,10,39,40]. The transversity PDF evolves with the DGLAP kernel evolution equation where the kernel δP is known up to two loop order [41,42]. Combining together eq. (3.5, 3.6) with eq. (3.7) we obtain the evolution properties of the matching kernels. They are where l ζ = ln(µ 2 /ζ). The analogous expression holds for δ ⊥ C coefficient function. In perturbation theory, the expression for the coefficient function can be presented as where a s = g 2 /(4π) 2 . The coefficients δC (n;k,l) with k + l > 0 are fixed order-by-order with the help of the renormalization group equations in eq. (3.8, 3.9). The expressions for these coefficients in their generic form up to two-loop are given, e.g. in ref. [17] (see appendix D.1) or in ref. [6] (see appendix B.2). The explicit expressions for the transversity function are presented in the supplementary file. Thus, the only non-trivial part to evaluate is δC (n;0,0) .

Evaluation of matching coefficient at NNLO
In order to find the expression for the matching coefficients, we evaluate the matrix elements with free quark states. The spinor indices of the in/out-going quarks can be contracted with iσ −β γ 5 , since it is the only non-vanishing spinor structure at twist-2. The resulting tensor diagram can be projected on tensor structures of transversity or pretzelosity. The obtained result is then compared with the matched expressions in eq. (2.7, 2.8). Schematically, we deal with the equation where h 1;f ←q (x) are the PDF evaluated on free-quark states. This equation can be solved recursively starting from the first non-zero contribution.
Evaluating Feynman diagrams we keep the momentum of quark collinear, p µ = p +nµ . This choice of kinematics significantly simplifies the calculation. In particular, it implies that b 2 is the only (Lorentz invariant) scale that is presented in the diagrams. Since the scaleless loop-integrals are zero in the dimensional regularization, many diagrams vanish. For example, this is the case of all pure virtual correction diagrams. Also all loop-integrals contributing to PDF are zero. Thus the only non-zero part of h 1 (x) is the renormalization constant. That is, the needed expression for PDF evaluated on quarks is a pure singularity, which can be found using renormalization group equations. It reads f ←r ⊗ δP where δP [n] are perturbative coefficients of transversity DGLAP kernel at a n s -order [41,42], and β 0 is the QCD β-function. Note that δP [1] q←q (x) = 0. The evaluation of TMD is made using the δ-regularization, which is described in details in [8,16]. It is a very convenient form of regularization, in particular, because it allows a clear separation of divergences. The outcome of each diagram at NNLO has a generic form where the functions f i are regulated in the limit x → 1 with +-distributions. The second and the third terms here represent the IR divergence. Therefore, the functions f 2 and f 3 exactly cancel in the sum of all the diagrams (and this fact can be also traced in the sum of sub-classes of diagrams). The last two terms represent the rapidity diverging pieces and thus the functions f 4 and f 5 are canceled by the rapidity renormalization factor. Altogether these cancellations serve as a good intermediate check of the computation. Summing together the diagrams we obtain the un-subtracted expression for TMDPDF on freequark states. Let us denote it as and it is UV and rapidity divergent. Starting from eq. (3.1) the renormalization procedure reads (here we omit the suffix ren for simplicity) where superscript in square brackets indicates the perturbative order of a quantity. In this expression we have Z = Z −1 2 Z q with Z 2 the quark-field renormalization constant, and Z q the TMD renormalization constant. The expression for them can be found e.g. in [16]. The soft function up to NNLO in δ-regularization is calculated in [8].

Matching of transversity TMD distribution at NNLO
The evaluation of the transversity matching coefficient is very similar to the evaluation of the unpolarized matching coefficient made in [16]. The main difference, which only simplifies the evaluation, is the absence of the mixing with the gluon operator.
where the (..) + -distribution is defined as usual This result agrees 1 with the ones obtained in refs. [31,38,43]. It is easy to see that logarithmic part of eq. (4.5) satisfies renormalization group equation eq. (3.8, 3.9). The finite part is Note, that in order to evaluate NNLO matching coefficient one needs the terms suppressed by , since they interfere with the singularities of the PDF, and produce a non-zero finite and 1/ contribution to eq. (4.4). The complete expression at all orders of can be found in [31]. The expression for δC [2] is lengthy. Since the only non-trivial part is the finite part δC (2;0,0) we restrict to it here. The logarithmic part can be restored using the renormalization group, and in the explicit form it is given in the supplementary Mathematica file. At NNLO we have the mixing with anti-quark operator, therefore, the matching is split into two channels Here, is the regular part of LO DGLAP kernel. In this expression the singularity at x → 1 is treated as the (..) + -distribution. These expressions are one of the main result of this work. It is intriguing to observe that the parts of eq. (4.9, 4.10) enclosed by the square brackets literally coincide with corresponding parts for the unpolarized matching coefficient, see [16] eqns. (7.3) and (7.8). In other words, the matching coefficient has the form where P [1] (x) is LO DGLAP kernel, for the corresponding PDF. Then we observe that the function F 1 (x) and the constant F 3 are the same for unpolarized and transversity kernels (for both flavor channels). Such behavior is expected since the contributions proportional to 1/(1 − x), as well as δ(x) contributions, that primary form the LO DGLAP kernel, comes from the diagrams where the quarks interact with the Wilson lines. Such diagrams are insensitive to the polarization structure of the operator. The rest of diagrams are not singular in the limit x → 1, and thus form a regular contribution. For more detailed discussion on the internal structure of transversity kernel see [42].

Matching of pretzelosity distribution at NNLO
The calculation for the matching of the pretzelosity TMDPDF over the transversity integrated PDF is in principle similar to the one of the transversity TMDPDF. In this case one has a different projector, see eq. (2.5), to be compared to g µν T used in the transversity calculation. Moreover the relation allows a simplification of many diagrams. In particular, the diagrams with a non-interacting quark line are exactly zero, according to the expression (5.2). This feature reduces the number of diagrams that we have to calculate for the pretzelosity TMD distribution. The pretzelosity projector is built as a sum of two terms. The first one is g µν T and it is the same as in the transversity calculation. As the topology of the diagrams is the same in both cases the integrals that appear in the calculation of the diagrams are also the same. The second term of the s b µ b ν /b 2 and this implies new types of master integrals, that has scalar products (b · q) 2 in the numerator (here, q µ is a loop-momentum). Such structures appears due to the convolution of a generic diagram with open indices with the projector (5.1).
The small-b expression fo the matching of the pretzelosity distribution is written in a form equivalent to the transversity case, Note, that the collinear function in eq. (5.3) is the transversity PDF. As in the transversity case, at NLO we have only the quark-to-quark channel, and at NNLO we have quark-to-quark and quark-to-antiquark channels. Due to eq. (5.2), the un-subtracted pretzelosity distribution is zero at LO, i.e. δ ⊥ Φ [0] (x) = 0. Consequently, the LO matching coefficient is also zero, i.e. δ ⊥ C q←q (x) = 0. This fact induces a simplification in the renormalization of the pretzelosity TMDPDF at NLO, demanding the absence of any divergences at this order. Moreover, due to the absence of the tree order collinear counterpart for the matching procedure the pretzelosity is suppressed by a s . As a result, the expression for the matching coefficient is given solely by the one-loop TMD matrix element We see that the obtained matching coefficient is -suppressed, so, in the limit → 0 it is zero, i.e. δ ⊥ C [1] q←q (x, b) = 0. This result is given in [31]. According to eq. (5.4), the pretzelosity distribution is suppressed numerically however this result does not ensure that a non-trivial coefficient can be obtained at higher orders.
The nullity of the LO pretzelosity distribution and the -suppressed behavior of the NLO contribution yields in a simple expression for the renormalized pretzelosity TMDPDF at NNLO In this expression it is important to keep all -terms of δ ⊥ Φ [1] f ←f , since they are multiplied by factors with leading 1/ 2 -behavior (the TMD renormalization factor, and soft factor). Thus, these terms produce 1/ singularities, despite the suppressed behavior of δ ⊥ Φ [1] f ←f . Naturally, these terms cancel the corresponding ultraviolet singularities of the un-subtracted TMD matrix element. We have also checked the exact cancellation of infrared divergences, and rapidity divergences. It is interesting to trace the distribution of the contributions between diagrams with different color factor. There are four types of contribution The contribution A N is -suppressed in a similar manner as the one-loop expression. The contribution A F is canceled by the renormalization factor entirely up terms suppressed in . Thus, the only non-zero contribution to the TMD matrix element comes from A F A and A A , which we find to be equal up to higher powers of . So, concluding we have found

7)
Therefore, the contribution proportional to C A disappears from the final expression, despite only these diagram are non-trivial. The resulting TMD matrix element in the pretzelosity channel is proportional to C 2 F only, and it reads h 1T,q←q = 0. (5.10) Expanding eq. (5.3) up to order a 2 s we obtain the following expressions for matching coefficient The convolution term that appears in eq. (5.11) is different from zero because the NLO matching coefficient is -suppressed but the NLO transversity integrated PDF is -divergent. So, the result for the convolution term is finite, which is the same that we get in eq. (5.9). Using eq. (5.11) we obtain a null value for the NNLO pretzelosity to transversity matching coefficient. So, q←f (x, b) = 0 + O( ), (5.14) where f = q,q.
We stress once more that the cancellation that lead to the zero result has a non-trivial structure. Because the topologies of diagrams that contribute to convolution term in (5.11) and to the TMD matrix element (5.9) are completely different. In the first case, these are ladder diagrams, while in the second case we have diagrams with tree-gluon vertex and non-planar diagrams. All this indicates the presence of a not yet understood concept behind these cancellations, and it suggests that such cancellations take a place at higher orders as well. Therefore we conjecture that at all orders of perturbation theory. To support this conjecture we also performed the calculation of the matching coefficient in the large-N f approximation using the same approach as in ref. [44]. We obtained that the resummed expression is also -suppressed. It gives additional confirmation of the conjecture in eq. (5.15).

Matching of transversity TMD fragmentation function at NLO and NNLO
The transversely polarized TMDFFs can be treated similarly to the case of transversely polarized TMDPDF.Despite the different origin and interpretation of these distributions, their perturbative treatment is analogous. Therefore, in this section we collect only the necessary definitions and results, and we do not provide the intermediate details.
The unsubtracted transversity TMDFFs are defined with the following hadronic matrix elements, which can be parameterized by the form factor decomposition where again s T is the transverse part of the hadron spin, λ is helicity, M is the mass of hadron and b 2 = −b 2 > 0. We recall that the parametrization presented here is valid for produced hadrons with spin-1/2. For the scalar or pseudo-scalar produced hadrons, the functions, are absent. The transversity TMDFF is represented by the function H 1 . The matching onto the fragmentation function is done as The factor z 2−2 is added to meet the common normalization of collinear FF function, that is defined as (nλ)|P, S; X iσ α+ γ 5 P, S; X|T q iW T n a (0)|0 .
The evolution kernels for the collinear FFs are known at two loops [41,42]. The main difference in evaluation of TMDFFs in comparison to TMDPDFs is the origin of parton momentum in diagrams, which is incoming in the PDF case, and outgoing in the FF case. Therefore, the expressions for TMDFFs could be obtained by the application of the crossing symmetry x → z −1 at the diagram level. This, however, should be done with caution since there is a brunch cut for x > 1, which should be transformed into a branch cut for z > 1. Additionally, one should take into account the z factors that are present in definitions of FFs, and that mix in the -expansions with various contributions. For the detailed discussion on relation between PDFs and FFs see [45] and references within. To avoid these complications, we re-evaluate the PDF master integrals with x → z −1 and reassemble the final result. The LO matching of transversity TMDFF is elementary f →f = δ f f δ(z).
(6.4) Therefore, the renormalization properties of the FF case are similar to the case of PDFs and the matching procedure follows the same pattern as in the unpolarized case [16]. For this reason we skip the details of evaluation and present the final result. For the presentation of the NLO and NNLO coefficient functions we introduce again logarithmic decomposition (see eq. (3.10)) δC f →f (z, L µ , l ζ ) = ∞ n=0 a n s n+1 k=0 n l=0 L k µ l l ζ δC (n;k,l) f →f (z). (6.5) The terms accompanied by logarithms, i.e. with k+l > 0 can be restored with renormalization group equation. For completeness, we present these lengthy expressions in the attached Mathematica file. The finite part of the coefficients at NLO is given by where is the LO DGLAP kernel for the transversity FF. At NNLO we have the same mixing with antiquark operator. The matching is split into two channels where, The singularity at z → 1 is understood as a (..) + -distribution (4.6). Similarly to the PDF case, the expressions enclosed by square brackets in eq. (6.9, 6.10) literally coincide with the ones of the unpolarized fragmenting function case, (see eq. (7.11) and (7.17) in [16]). In other words, it can be written the form (4.12), and the functions F 1 (z) and F 3 coincide for polarized and unpolarized cases.
We have not evaluated the pretzelosity TMDFF, because its calculation is rather involved. However, we have no doubts that it matching coefficient is zero at NNLO alike the matching of pretzelosity TMDPDF.

Conclusions
In this work, we have studied the twist-2 matching of transversity and pretzelosity (or quadrupole) TMD distributions. We have derived the matching coefficients for these distributions at next-tonext-to-leading order (NNLO) in the strong coupling. We have checked that the renormalization of rapidity divergences works exactly in the same way as for unpolarized distributions, as it is predicted by the transverse momentum dependent factorization theorem. The present calculation has a structure similar to the one of the NNLO matching of unpolarized TMDs made in [16]. In the article, we present only the finite part of the coefficient functions, while the logarithmic part can be restored with the renormalization group equations. The full result, including the logarithmic part, is also reported in the attached Mathematica file.
In the case of transversity, we have considered both the TMDPDF and the TMDFF cases. We have found several analogies and identities between the matching coefficients of transversity and unpolarized distributions, that can serve as a cross-check of both results. The matching coefficients for transversity (given in eqs. (4.8, 4.9, 4.10) and eq. (6.8, 6.9, 6.10)) can be readily used in phenomenological applications. A recent review of the phenomenology of transversity in fragmentation can be found in [28]. Our result is the first calculation of the NNLO matching for transversely polarized TMD operator. To our knowledge the NLO matching coefficient for transversity TMDFF (6.6)) is also calculated here for the first time. We stress that this is also the first NNLO evaluation of the matching for a polarized TMD distribution. Therefore, given the result of this work, the transversity TMD distribution is known to the same perturbative order as unpolarized distributions. This fact is important to establish phenomenologically the universality of TMD evolution.
For the pretzelosity, we have found that the expected two-loop matching coefficient is actually zero, despite the fact that the matrix element over free quark for pretzelosity distribution is non-zero. It is an unexpected result, since the analogous quadrapole distribution in the gluon sector (namely linearly polarized TMD gluon distribution) has a non-zero matching already at one-loop level. We have also checked that the LO of large-N f expansion (given by diagrams with an arbitrary number of fermion bubble insertions, for details see [6]) is also null. Although these facts do not demonstrate completely that the twist-2 part of pretzelosity is zero at all orders in perturbation theory, certainly they are an evidence of this statement. We conjecture that the pretzelosity distribution does not match the twist-2 distribution, and thus has the leading matching only at the twist-4 level. At the moment we have not found an argument to justify this fact beyond the present calculation. Nevertheless, it agrees with the phenomenological and experimental results that suggest a highly suppressed pretzelosity distribution [29,30]. We have performed the calculation of pretzelosity only for the case of TMDPDF, nonetheless, we expect the same result for the TMDFF.