1 Introduction

One of the modern challenges of QCD is represented by the study of the effects of polarization in differential cross sections for semi-inclusive deep inelastic scattering (SIDIS) and Drell-Yan/vector/scalar-boson production. The cross sections can be factorized in terms of transverse momentum dependent (TMD) distributions [1,2,3,4] which describe the (transverse and collinear) momentum distribution of quarks and gluons inside a nucleon. The perturbative inputs of the factorization formula play an important role, and presumably should be included at the maximum allowed order to provide the best agreement/prediction of the theory with experiment. The importance of perturbative input both at high and low energy, was demonstrated for instance in [5] within the analysis of unpolarized TMD parton distribution functions (TMDPDFs).

Some perturbative parts of the TMD factorization are universal and independent of polarization. It concerns primary the hard coefficient function and evolution kernels, which nowadays are known up three-loop order [6,7,8]. The additional parts that require the perturbative input are the actual models for TMD distributions. The perturbative computation gives us some relevant information in the limit of large transverse momentum (that is, in the limit of small transverse distances or small-b). In this limit, the TMD distributions match collinear distributions providing a starting point for phenomenology and greatly increasing the agreement with high-energy data. Nowadays, the matching for the most part TMD distributions are known, although the information is often difficult to extract from the literature. The matching of TMD distributions to the twist-2 functions is known uniformly at the one-loop level [9], and some of them are known at two-loop level [10,11,12]. The matching onto the twist-3 functions is less elaborated. A preliminary study was provided in [13] prior to the modern formulation of factorization theorem. The main aim of this paper is to uniformly evaluate the matching of quark TMDPDFs at the twist-3 level, using solely the definition of the TMD distribution as it is provided by the TMD factorization theorem.

The desired matching expressions are given by the first terms of the operator product expansion (OPE) for the TMD operator at small-b, or in the vicinity of the light-cone. It can be systematically done starting from the field-theoretical definition of the TMD operator, using the algebra of fields and QCD equations of motion. We recall that the OPE is naturally formulated in position space, and can be performed without any explicit reference to a particular process. So, the method is universal and it allows one to calculate any TMD distribution at any order. Here we consider only the contributions of twist-2 and twist-3 functions, and it covers almost all (7 out of 8) quark TMDPDFs. The gluon distributions, as well as, TMD fragmentation functions (TMDFFs) can be considered in principle in the same fashion. It is important to point out that the OPE does not depend on hadronic states, therefore, many results of this work can be applied directly, or with a minimal effort, to closely related areas, such as studies of generalized transverse momentum distributions (GTMDs) and Wigner function [14, 15]. Additionally, we think that OPE approach is technically simpler and more systematic in comparison to a direct evaluation and factorization of cross-sections, which is often used see e.g. [16,17,18].

In order to realize our computations we start observing that the operators that define TMD distributions have a peculiar structure, which distinguishes them from more traditional parton distribution functions (PDFs) and fragmentation functions (FFs), distributions amplitudes and others. Namely, they include half-infinite light-like Wilson lines, and they are geometrically non-compact. Moreover, the direction of the Wilson lines depends on the underline process. This direction is future pointing for production mechanisms (such as fragmentation in semi-inclusive deep-inelastic scattering (SIDIS)), and past pointing for collision mechanisms (such as Drell-Yan process). It gives a superficial process dependence of the TMD distributions, in the form of sign-flip for P-odd distributions [19]. At the same time, the collinear distributions are perfectly independent on the process, therefore, any process dependence of TMD operator must reveal itself within OPE. We indeed observe this effect and demonstrate that it appears in the contributions specific for the TMD operators. Exactly these contributions give rise to the famous Efremov-Teryaev-Qiu-Sterman (ETQS) function [20,21,22,23]. Moreover, we have observed that these process-dependent terms of OPE, also contribute to the P-even, and hence process-independent, distributions. Altogether, to our best knowledge, the analysis presented here is the first study of OPE for TMD operators beyond the twist-2 accuracy. As a final result, we obtain all leading order matching expressions for the TMDPDFs in the quark sector.

The definition framework and all relevant TMD distributions and operators are given in Sect. 2. The light-cone OPE for TMD operators up to linear in transverse positions terms is given in Sect. 3. The parameterization of relevant collinear distributions is presented in Sect. 4. The assembling of OPE and its application to particular distributions is presented in Sect. 5. We also present by-product result for Mellin moments of worm-gear function in Sect. 6. The final results are shown and commented in Sect. 7.

2 Definition of TMD distributions

We outline our work for the transverse momentum dependent parton distribution functions (TMDPDFs) which, in the quark case, are defined by the matrix element [1, 24, 25]

$$\begin{aligned}&\Phi _{q\leftarrow h,ij}(x,{\varvec{b}})\nonumber \\&\quad =\int \frac{dz}{2\pi }e^{-ixz(pn)}\langle P,S|\bar{T}\{\bar{q}_j\left( zn+{\varvec{b}}\right) [\lambda n+{\varvec{b}},\pm \infty n]\}\nonumber \\&\qquad \times T\{ [\pm \infty n,0] q_i(0)\}|P,S\rangle , \end{aligned}$$
(2.1)

where \(\pm \infty n\) indicates different light-cone infinities. The TMD distributions which appear in SIDIS have Wilson lines pointing to \(+\infty n\), while in Drell-Yan they point to \(-\infty n\). The Wilson lines within the TMD operator are along a light-like direction n. Another light-like vector is associated with the large-component of the hadron momentum P,

$$\begin{aligned} p^\mu =(np)\bar{n}^\mu =P^\mu -\frac{n^\mu }{2}\frac{M^2}{(nP)}, \end{aligned}$$
(2.2)

where \((nP)=(np)\), and M is the mass of hadron (\(P^2=M^2\)). Together vectors n and \(\bar{n}\) define the scattering plane. The relative normalization of vectors is

$$\begin{aligned} n^2=\bar{n}^2=0,\quad (n\bar{n})=1. \end{aligned}$$
(2.3)

Thus, any four-vector can be decomposed into the components

$$\begin{aligned} v^\mu =v^+ \bar{n}^\mu +v^- n^\mu +v_T^\mu , \end{aligned}$$
(2.4)

where \(v^+=(nv)\), \(v^-=(\bar{n} v)\), and \(v_T\) is the transverse component orthogonal to the scattering plane \((v_Tn)=(v_T\bar{n})=0\). To specify the reference frame we state that \(v^\pm =(v^0\pm v^3)/\sqrt{2}\).

The transverse components play a special role in our consideration. The transverse subspace is projected out by the transverse part of the metric tensor

$$\begin{aligned} g_T^{\mu \nu }=g^{\mu \nu }-\frac{n^\mu p^\nu +p^\mu n^\nu }{(np)}. \end{aligned}$$
(2.5)

There are only two non-zero components, \(g^{11}_T=g_T^{22}=-1\). In the following, we also need the transverse part of the Levi-Civita tensor

$$\begin{aligned} \epsilon _T^{\mu \nu }=\frac{n_\alpha p_\beta }{(np)}\epsilon ^{\alpha \beta \mu \nu }, \end{aligned}$$
(2.6)

where \(\epsilon ^{\mu \nu \rho \sigma }\) is defined in the Bjorken convention (\(\epsilon _{0123}=-\epsilon ^{0123}=1\)). Consequently, we have \(\epsilon _T^{12}=-\epsilon _T^{21}=1\), which coincides with the definition [24, 26, 27], despite the opposite normalization of the four-dimension \(\epsilon \)-tensor. The tensor \(\epsilon _T^{\mu \nu }\) does not change sign when both indices are down, \(\epsilon _{T\mu \nu }=\epsilon _T^{\mu \nu }\), and \(\epsilon ^{\mu \nu }_T\epsilon _{T\mu \rho }=\delta ^{~~\nu }_{T~\rho }\). Since the transverse subspace is Euclidian, the scalar product transverse vectors is negative, \(v_T^2<0\). In the following, we use the bold font notation to designate the Euclidian scalar product of transverse vectors, i.e. \({\varvec{b}}^2=-b^2>0\), when it is convenient.

The spin of the hadron is parameterized by the spin-vector S,

$$\begin{aligned} S^2=-1,\quad (PS)=0. \end{aligned}$$
(2.7)

The light-cone decomposition of the spin vector is

$$\begin{aligned} S^\mu =\frac{\lambda }{M}p^\mu -\frac{\lambda }{2}\frac{M}{(np)}n^\mu +s_T^\mu , \end{aligned}$$
(2.8)

where the helicity \(\lambda \) of the hadron is

$$\begin{aligned} \frac{(nS)}{(np)}=\frac{\lambda }{M}. \end{aligned}$$
(2.9)

The vector \(s_T^\mu \) is the transverse component of the spin, \(s_T^2=\lambda ^2-1\). With the help of \(\epsilon _T\)-tensor we can introduce another useful (axial) vector

$$\begin{aligned} \tilde{s}^\mu _T=\epsilon ^{\mu \nu }_TS_\nu , \end{aligned}$$
(2.10)

and it implies \(\tilde{s}^2_T=s^2_T\).

The open spinor indices (ij) of the TMD operator in Eq. (2.1) are to be contracted with different gamma-matrices, which we denote generically as \(\Gamma \). The gamma-matrices that appear at the leading order of TMD factorization are

$$\begin{aligned} \Gamma =\{\gamma ^+,\gamma ^+\gamma _5,i\sigma ^{\alpha +}_T\gamma _5\}, \end{aligned}$$
(2.11)

where \(\sigma ^{\alpha +}_T=g^{\alpha \beta }_T \sigma _{\beta \gamma }n^\gamma \), and

$$\begin{aligned} \sigma ^{\mu \nu }= & {} \frac{i}{2}(\gamma ^\mu \gamma ^\nu -\gamma ^\nu \gamma ^\mu ),\nonumber \\ \gamma _5= & {} i\gamma ^0\gamma ^1\gamma ^2\gamma ^3=\frac{i}{4!} \epsilon _{\mu \nu \alpha \beta }\gamma ^\mu \gamma ^\nu \gamma ^\alpha \gamma ^\beta . \end{aligned}$$
(2.12)

In the naive parton model interpretation, these gamma-structures are related to the observation of unpolarized (\(\gamma ^+\)), longitudinally polarized (\(\gamma ^+\gamma ^5\)) and transversely polarized (\(i\sigma ^{\alpha +}_T\gamma ^5\)) quarks inside the hadron. Beyond the leading order factorization one expects that the power suppressed terms of TMD show also different gamma structures. However, currently, the TMD factorization theorem is not established beyond the leading order. Moreover, it is known that TMD distributions with a gamma-structure different from (2.11) contain rapidity divergences that are not renormalized by the standard TMD soft factor [9].

Historically, the TMD distributions have been introduced and parameterized in momentum space [13]. Denoting

$$\begin{aligned} \Phi ^{[\Gamma ]}_{q\leftarrow h}= & {} \frac{1}{2}\mathrm {Tr}\left( \Phi _{q\leftarrow h}\Gamma \right) , \end{aligned}$$
(2.13)

we have [27, 28]

$$\begin{aligned}&\Phi _{q\leftarrow h}^{[\gamma ^+]}(x,p_T)= f_1(x,p_T) -\frac{\epsilon _T^{\mu \nu }p_{T\mu }s_{T\nu }}{M}f_{1T}^\perp (x,p_T), \nonumber \\ \end{aligned}$$
(2.14)
$$\begin{aligned}&\Phi _{q\leftarrow h}^{[\gamma ^+\gamma _5]}(x,p_T)= \lambda \,g_{1L}(x,p_T)-\frac{p_{T\mu }s_T^\mu }{M}g_{1T}(x,p_T),\nonumber \\ \end{aligned}$$
(2.15)
$$\begin{aligned}&\Phi ^{[i\sigma ^{\alpha +}\gamma _5]}_{q\leftarrow h}(x,p_T)=s_T^\alpha h_1(x,p_T)+\lambda \frac{p_T^\alpha }{M}h_{1L}^\perp (x,p_T)\nonumber \\&\quad -\frac{\epsilon _T^{\alpha \mu }p_{T\mu }}{M}h_1^\perp (x,p_T)\nonumber \\&\quad +\frac{p_T^2}{M^2}\left( \frac{g_T^{\alpha \mu }}{2} -\frac{p_T^{\alpha }p_T^{\mu }}{p_T^2}\right) s_{T\mu }h_{1T}^\perp (x,p_T),\nonumber \\ \end{aligned}$$
(2.16)

where \(p_T^2=-{\varvec{p}}_T^2<0\). Note, that the functions \(f(x,p_T)\) depend only on the modulus of \(p_T\), but not on the direction. The functions presented here are traditionally called unpolarized (\(f_1\)), Sivers (\(f_{1T}^\perp \)), helicity (\(g_{1L}\)), worm-gear T (\(g_{1T}\)), transversity \((h_1)\), worm-gear L (\(h_{1L}^\perp \)), Boer-Mulders (\(h_1^\perp \)) and pretzelosity (\(h_{1T}^\perp \)) distributions.

For practical calculations it is convenient to write TMD distributions in momentum space as Fourier transform of distributions in position space in the usual manner

$$\begin{aligned} \Phi _{q\leftarrow h,ij}(x,p_T)=\int \frac{d^2 {\varvec{b}}}{(2\pi )^2}e^{+i({\varvec{b}}{\varvec{p}}_T)}\Phi _{q\leftarrow h,ij}(x,{\varvec{b}}), \end{aligned}$$
(2.17)

where the scalar product \(({\varvec{b}}{\varvec{p}}_T)\) is Euclidian. The decomposition in Eqs. (2.142.16) is then replaced by its analog it position space,

$$\begin{aligned}&\Phi _{q\leftarrow h}^{[\gamma ^+]}(x,{\varvec{b}})=f_1(x,{\varvec{b}})+i\epsilon _T^{\mu \nu } b_\mu s_{T\nu } M f_{1T}^\perp (x,{\varvec{b}}), \end{aligned}$$
(2.18)
$$\begin{aligned}&\Phi ^{[\gamma ^+\gamma _5]}_{q\leftarrow h}(x,{\varvec{b}})=\lambda g_{1L}(x,{\varvec{b}})+i b_{\mu }s^{\mu }_T M g_{1T}(x,{\varvec{b}}), \end{aligned}$$
(2.19)
$$\begin{aligned}&\Phi ^{[i\sigma ^{\alpha +}\gamma _5]}_{q\leftarrow h}(x,{\varvec{b}})=s_{T}^\alpha h_1(x,{\varvec{b}}) -i\lambda b^{\alpha }M h_{1L}^\perp (x,{\varvec{b}})\nonumber \\&\quad +i\epsilon _{T}^{\alpha \mu }b_\mu M h_1^\perp (x,{\varvec{b}})\nonumber \\&\quad +\frac{M^2{\varvec{b}}^2}{2}\left( \frac{ g_{T}^{\alpha \mu }}{2}+\frac{b^\alpha b^\mu }{{\varvec{b}}^2}\right) s_{T\mu } h_{1T}^\perp (x,{\varvec{b}}).\nonumber \\ \end{aligned}$$
(2.20)

This parameterization coincidesFootnote 1 with the parameterization given in [29]. The explicit transformation rules for all these functions can be found in Appendix A.

3 Light-cone expansion for TMD operator

The small-b matching of TMD distribution to the integrated distributions is obtained by the operator product expansion (OPE) at small-b. The OPE is independent from the hadronic states and for this reason it is universal. Let us introduce a separate notation for the TMD operators. The operator that produces TMD distributions in the Drell-Yan case is

$$\begin{aligned} \mathcal {U}_{\text {DY}}^{\Gamma }(z,{\varvec{b}})= & {} \bar{q}(z n+{\varvec{b}})[z n+{\varvec{b}},-\infty n+{\varvec{b}}]\nonumber \\&\times \Gamma [-\infty n-{\varvec{b}},-z n-{\varvec{b}}] q(-zn -{\varvec{b}}), \end{aligned}$$
(3.1)

where \(\Gamma \) represents the gamma-matrices of the leading set (2.11), and the half-infinite Wilson lines are defined as

$$\begin{aligned} {[}a_1 n+{\varvec{b}},a_2n+{\varvec{b}}]= & {} P\exp \left( ig\int _{a_2}^{a_1} d\sigma n^\mu A_\mu (\sigma n+{\varvec{b}})\right) .\nonumber \\ \end{aligned}$$
(3.2)

Here and in the following we also omit the T-ordering of the fields, since it is irrelevant in the tree order approximation. The operator that produces the TMD distributions for SIDIS reads

$$\begin{aligned} \mathcal {U}_{\text {DIS}}^{\Gamma }(x,{\varvec{b}})= & {} \bar{q}(z n+{\varvec{b}})[z n+{\varvec{b}},+\infty n+{\varvec{b}}]\nonumber \\&\times \Gamma [+\infty n-{\varvec{b}},-z n-{\varvec{b}}] q(-zn -{\varvec{b}}). \end{aligned}$$
(3.3)

Generally, the links which connect the end points of Wilson lines at a distant transverse plane must be added in both operators (for DY and for SIDIS). Here, we omit them for simplicity, assuming that some non-singular gauge (e.g. covariant gauge) is in use. In non-singular gauges the field nullifies at infinities, \(A_\mu (\pm \infty n)=0\), and the contribution of distant gauge links vanish.

The relation between the TMD distribution (2.1) and the TMD operator (3.1) is

$$\begin{aligned} \Phi ^{[\Gamma ]}_{q\leftarrow h}(x,{\varvec{b}})=\int \frac{dz}{2\pi }e^{-2ixzp^+}\langle P,S|\mathcal {U}^\Gamma \left( z,\frac{{\varvec{b}}}{2}\right) |P,S\rangle .\nonumber \\ \end{aligned}$$
(3.4)

The light-cone expansion of the TMD operators corresponds to the expansion in the variable \({\varvec{b}}\). The OPE has a generic schematic form

$$\begin{aligned}&\mathcal {U}(z,{\varvec{b}}) =\sum _i \left[ C_i*\mathcal {O}^{\text {tw2}}_i\right] (z)+{\varvec{b}}^\mu \sum _i \left[ \tilde{C}_i*\mathcal {O}^{\text {tw3}}_{\mu ,i}\right] (z) +O({\varvec{b}}^2),\nonumber \\ \end{aligned}$$
(3.5)

where C’s are Wilson coefficient functions which depend on \(\ln {\varvec{b}}^2\), \(\mathcal {O}\) ’s are light-cone operators, and the symbol \(*\) denotes some integral convolution between coefficient function and operators. Here, the superscripts tw2 and tw3 indicate the collinear twist, which in principle differs from the geometrical twist. We remind that the term collinear twist indicates the distributions which enter the same order of momentum expansion. It is not a well-defined quantum number, in contrast to the geometrical twist. The later is defined by “dimension-spin” value, and is a well-defined quantum number, in the sense that e.g. it conserves and does not mix under the scaling transformations. As we will see the operators \(\mathcal {O}^{\text {tw3}}\) are compositions of geometrical twist-2 and twist-3 operators. The coefficient functions are perturbatively calculable. In this work, we study the matching only at order \(\alpha _s^0\).

At leading order in \(\alpha _s\) the quantum fields can be considered as classical fields, that satisfy QCD equations of motion. In this approximation, the small-b OPE is just the Taylor expansion at \(b=0\). Expanding \(\mathcal {U}\) in powers of \({\varvec{b}}\) up to the linear order we obtain

$$\begin{aligned} \mathcal {U}^\Gamma (z,{\varvec{b}})=\mathcal {U}^\Gamma (z,{\varvec{0}})+b^\mu \frac{\partial }{\partial b^\mu } \mathcal {U}^\Gamma (z,{\varvec{b}})\Big |_{{\varvec{b}}=0}+O({\varvec{b}}^2). \end{aligned}$$
(3.6)

The leading term reads

$$\begin{aligned} \mathcal {U}_{\text {DY}}^{\Gamma }(z,{\varvec{0}})= & {} \bar{q}(z n)[z n,-\infty n]\Gamma [-\infty n,-z n] q(-zn) \nonumber \\= & {} \bar{q}(z n)[z n,-z n]\Gamma q(-z n), \end{aligned}$$
(3.7)

where the half-infinite segments of Wilson line compensate each other due to the unitarity of a Wilson line. The same holds for the SIDIS operator

$$\begin{aligned} \mathcal {U}_{\text {DIS}}^{\Gamma }(z,{\varvec{0}})= & {} \bar{q}(z n)[z n,+\infty n]\Gamma [+\infty n,-z n] q(-zn)\nonumber \\= & {} \bar{q}(z n)[z n,-z n]\Gamma q(-zn). \end{aligned}$$
(3.8)

Therefore, we obtain that \(\mathcal {U}_{\text {DY}}^{\Gamma }(z,{\varvec{0}}) =\mathcal {U}_{\text {DIS}}^{\Gamma }(z,{\varvec{0}})\), which is well known.

The term linear in \({\varvec{b}}^\mu \) is given by the derivative of the operator. Explicitly, it reads

$$\begin{aligned} \frac{\partial }{\partial b^\mu } \mathcal {U}^{\text {DY}}_\Gamma (z,{\varvec{b}})\Big |_{{\varvec{b}}=0}= & {} \bar{q}(z n)[z n,-\infty n](\overleftarrow{\partial _{T\mu }}-\overrightarrow{\partial _{T\mu }})\nonumber \\&\times \Gamma [-\infty n,-z \lambda ] q(-z), \end{aligned}$$
(3.9)

where \(\partial _{T\mu }\) is the derivative with respect to the transverse components only. This expression can be written as

$$\begin{aligned}&\frac{\partial }{\partial b^\mu } \mathcal {U}_{\text {DY}}^\Gamma (z,{\varvec{b}})\Big |_{{\varvec{b}}=0}\nonumber \\&\quad = \bar{q}(zn)\overleftarrow{D_\mu }[zn,-\infty n]\Gamma [-\infty n,-zn]q(-zn)\nonumber \\&\qquad +ig\int _{-\infty }^z \bar{q}(zn)[zn,\tau n]F_{\mu +}(\tau n)[\tau n,-\infty n]\nonumber \\&\qquad \times \Gamma [-\infty n,-zn]q(-zn)\nonumber \\&\qquad -ig\int ^{-\infty }_{-z} \bar{q}(zn)[zn,-\infty n]\nonumber \\&\qquad \times \Gamma [-\infty n,\tau n]F_{\mu +}(\tau n)[\tau n,-zn]q(-zn) \nonumber \\&\qquad -\bar{q}(zn)[zn,-\infty n]\Gamma [-\infty n,-zn]\overrightarrow{D_\mu }q(-zn), \end{aligned}$$
(3.10)

where the covariant derivative and the field-strength tensor are defined as usual

$$\begin{aligned} \overrightarrow{D}_\mu= & {} \overrightarrow{\partial }_\mu -igA_\mu ,\quad \overleftarrow{D}_\mu =\overleftarrow{\partial }_\mu +igA_\mu ,\nonumber \\ F_{\mu \nu }= & {} \partial _\mu A_\nu -\partial _\nu A_\mu -ig[A_\mu ,A_\nu ]. \end{aligned}$$
(3.11)

To obtain the expression (3.10) we have used the assumption thatFootnote 2 \(A(-\infty n)=0\), and the explicit expression for the total derivative of a Wilson line,

$$\begin{aligned} \partial _\mu \{[z_1 n,z_2 n]\}= & {} \frac{d}{dy^\mu }[z_1n+y,z_2n+y]\Big |_{y=0}\nonumber \\= & {} ig\left( A_\mu (z_1n)[z_1n,z_2n]-[z_1n,z_2n]A_\mu (z_2n)\right. \nonumber \\&\left. +\int _{z_2}^{z_1}d\tau [z_1n,\tau n]F_{\mu +}(\tau n)[\tau n,z_2 n]\right) , \end{aligned}$$
(3.12)

where the vector n can be arbitrary.

The segments of Wilson line between \(-\infty \) and \(\tau \) cancel and we obtain

$$\begin{aligned}&\frac{\partial }{\partial b^\mu } \mathcal {U}_{\text {DY}}^\Gamma (z,{\varvec{b}})\Big |_{{\varvec{b}}=0}\nonumber \\&\quad = \bar{q}(zn)\left( \overleftarrow{D_\mu }[zn,-zn] -[zn,-zn]\overrightarrow{D_\mu }\right) \Gamma q(-zn) \nonumber \\&\qquad +ig\left( \int _{-\infty }^z +\int _{-\infty }^{-z}\right) d\tau ~ \bar{q}(zn)[zn,\tau n]\Gamma F_{\mu +}(\tau n) \nonumber \\&\qquad \quad [\tau n,-zn]q(-zn). \end{aligned}$$
(3.13)

In the case of SIDIS kinematics the Wilson lines point the future light-like infinity, and therefore, the same derivation gives

$$\begin{aligned}&\frac{\partial }{\partial b^\mu } \mathcal {U}_{\text {DIS}}^\Gamma (z,{\varvec{b}})\Big |_{{\varvec{b}}=0}\nonumber \\&\quad =\bar{q}(zn) \left( \overleftarrow{D_\mu }[zn,-zn]-[zn,-zn]\overrightarrow{D_\mu }\right) \Gamma q(-zn)\nonumber \\&\qquad -ig\left( \int ^{\infty }_z +\int ^{\infty }_{-z}\right) d\tau ~ \bar{q}(zn)[zn,\tau n]\Gamma F_{\mu +}(\tau n)\nonumber \\&\qquad \quad [\tau n,-zn]q(-zn). \end{aligned}$$
(3.14)

Comparing the results for DY in Eq. (3.13) and SIDIS in Eq. (3.14) kinematics we observe that the first term is the same, while the last terms differ because of the limits of integration and a common sign. Therefore, already at this stage it is clear that the operator in the first term does not contribute to P-odd distributions (i.e. Sivers and Boer-Mulders functions) which is known to change sign in different kinematics.

As expected, the non-compact (in the sense that it spans an infinite range in position space) TMD operator is expressed via a set of compact light-cone operators. The light-cone operators in Eqs. (3.133.14) are not very well defined, in the sense, that they are of indefinite geometrical twist (more specifically, this issue concerns the first terms of Eqs. (3.133.14)). At the next stage of the OPE we need to classify the contributions with respect to twist and decompose over independent components. These components are parameterized via parton distributions functions, which are universal and can be measured in different experiment.

As the key point here is the twist-expansion we provide some additional discussion. The standard approach to twist-decomposition of operators is to consider their local expansion. In the local expansion the contributions of different twists can be separated by the permutation algebra, and summed back to a non-local representation, see e.g. the detailed decomposition of similar operators in [34]. However, a much simpler approach consists in taking the operator directly in a non-local form [35, 36]. In this approach, one starts with operators off the light-cone, and makes the twist-decomposition, and then performs the limit to the light-cone.

In principle, the procedure of twist-decomposition can be made at the level of operators, see e.g. [35]. However, practically it is involved, especially for tensor gamma-structure. The evaluation is significantly simpler in terms of distributions, e.g. as it is done in Ref. [36]. Here we are going to follow this second approach. In fact, the derivation presented in the next sections closely follows the procedure described in details in [36] for the case of meson distribution amplitudes. The difference in kinematics does not allow us to use the powerful method of conformal basis, but there is no principle difference in other aspects.

Prior to the parameterization and twist-decomposition let us prepare the operator for this procedure, and make its off-light-cone generalization. At our order of accuracy (twist-3) we do not need the most general form of the three-point operators, since they are already of geometrical twist-3 and do not contain admixture with twist-2 operators. Therefore, the generalization should be done only for the two-point operators, and it can be simply achieved by the replacement \(z n^\mu \rightarrow y^\mu \) with \(y^2\ne 0\). The result is conveniently re-written in the following form

$$\begin{aligned}&\bar{q}(y)\left( \overleftarrow{D_\mu }[y,-y]-[y,-y] \overrightarrow{D_\mu }\right) \Gamma q(-y)\nonumber \\&\quad = \frac{\partial }{\partial y^\mu } \bar{q}(y)[y,-y]\Gamma q(-y) \nonumber \\&\qquad -ig \int _{-1}^1 dv \, v y^\nu \bar{q}(y)[y,vy]\Gamma F_{\mu \nu }(v y)[vy,-y]q(-y), \end{aligned}$$
(3.15)

where we have used the formula for the stretch derivative of the Wilson line

$$\begin{aligned} \frac{\partial }{\partial y^\mu } {[}y,-y]= & {} ig\left( A_\mu (y)[y,-y]+[y,-y]A_\mu (-y)\right. \nonumber \\&\left. +\int _{-1}^1dv v y^\nu [y,v y]F_{\mu \nu }(v y)[v y,-y]\right) .\nonumber \\ \end{aligned}$$
(3.16)

Note, that this expression is the same for DY and SIDIS operators. The last term of (3.15) is again pure twist-3 operator, and thus one can set it directly on the light-cone.

Let us conclude this section with an intermediate summary of our main results. For convenience we introduce the generic notation for two- and three-point operators

$$\begin{aligned}&\mathcal {O}_{\Gamma }(z)=\bar{q}(z n)[z n,-z n]\Gamma q(-zn), \end{aligned}$$
(3.17)
$$\begin{aligned}&\mathcal {T}_\Gamma ^\mu (z_1,z_2,z_3)=g\bar{q}(z_1 n)[z_1 n, z_2 n]\Gamma F^{\mu +}(z_2 n)\nonumber \\&\quad \times [z_2 n,z_3 n]q(z_3 n). \end{aligned}$$
(3.18)

The expression for the first terms of small-b expansion for TMD operator reads (at leading order in \(\alpha _s\) )

$$\begin{aligned}&\mathcal {U}_{\text {DY}}^\Gamma (z,{\varvec{b}})\nonumber \\&\quad =\mathcal {O}_\Gamma (z)+b_\mu \left\{ \lim _{y\rightarrow zn}\frac{\partial }{\partial y_\mu }\mathcal {O}_\Gamma (y)-i\int _{-1}^1 dv \,v z \,\mathcal {T}_\Gamma ^\mu (z,vz,-z) \right. \nonumber \\&\qquad \left. +i\left( \int _{-\infty }^z+\int _{-\infty }^{-z}\right) d\tau \mathcal {T}_\Gamma ^\mu (z,\tau ,-z)\right\} +O({\varvec{b}}^2), \end{aligned}$$
(3.19)
$$\begin{aligned}&\mathcal {U}_{\text {DIS}}^\Gamma (z,{\varvec{b}})\nonumber \\&\quad =\mathcal {O}_\Gamma (z)+ b_\mu \left\{ \lim _{y\rightarrow zn}\frac{\partial }{\partial y_\mu }\mathcal {O}_\Gamma (y)-i\int _{-1}^1 dv \,vz\,\mathcal {T}_\Gamma ^\mu (z,vz,-z) \right. \nonumber \\&\left. \qquad -i\left( \int ^{\infty }_z+\int ^{\infty }_{-z}\right) d\tau \mathcal {T}_\Gamma ^\mu (z,\tau ,-z)\right\} +O({\varvec{b}}^2). \end{aligned}$$
(3.20)

The limit \(y\rightarrow z n\) implies \(y^2\rightarrow 0\) such that the light-like separation between fields is z.

4 Collinear distributions

Evaluating the matrix elements of Eqs. (3.193.20), and hence the matching of TMD distributions, we meet with a number of collinear parton distributions. In this section we present the parameterization of two and three point parton distributions that appear in the final result. In fact, the functions that we find represent a complete set of geometrical twist 2 and 3 quark distributions. For the two-point functions we use the standard parameterization by [37]. For the three point functions there is not a commonly accepted parameterization, therefore, we introduce a parameterization inspiring in [38].

4.1 Parameterization of quark–quark correlators

The standard parameterization of light-cone quark–quark correlators is given [37] and reads

$$\begin{aligned}&\langle P,S|O^{\gamma ^\mu }(z)|P,S\rangle \nonumber \\&\quad =2\int dx e^{2ix z p^+} \Big \{p^\mu f_1(x)+\frac{n^\mu }{(np)}M^2 f_4(x)\Big \}, \end{aligned}$$
(4.1)
$$\begin{aligned}&\langle P,S|O^{\gamma ^\mu \gamma ^5}(z)|P,S\rangle \nonumber \\&\quad =2\int dx e^{2ixzp^+}\Big \{\lambda p^\mu g_1(x)+s_T^\mu M g_T(x)\nonumber \\&\qquad +\lambda M^2 \frac{n^\mu }{(np)}g_3(x)\Big \}, \end{aligned}$$
(4.2)
$$\begin{aligned}&\langle P,S|O^{i\sigma ^{\mu \nu }\gamma ^5}(z)|P,S\rangle \nonumber \\&\quad = 2\int dx e^{2ixzp^+}\Big \{(s_T^\mu p^\nu -p^\mu s_T^\nu ) h_1(x)\nonumber \\&\qquad +\lambda \frac{M}{(np)}(p^\mu n^\nu -n^\mu p^\nu )h_L(x)\nonumber \\&\qquad +(s_T^\mu n^\nu -n^\mu s_T^\nu )\frac{M^2}{(np)} h_3(x)\Big \}, \end{aligned}$$
(4.3)

where the operators \(\mathcal {O}^\Gamma \) are defined in Eq. (3.17). The twist-2 PDFs \(f_1\), \(g_1\) and \(h_1\) are known as unpolarized, helicity and transversity PDFs. The PDFs \(g_T\) and \(h_L\) are of collinear twist-3. The PDFs \(f_4\), \(g_3\) and \(h_3\) are of collinear twist-4, and do not appear in the current final results. The collinear twist-3 PDFs are not independent as they are combinations of PDFs of twist-2 and three-point PDFs. The derivation of this relation can be done with the help of QCD equations of motion and is presented in the Appendix C.

The PDF defined by Eqs. (4.1, 4.2, 4.3) are non-zero for \(-1<x<1\) and zero for \(|x|>1\) [39]. They can be represented by

$$\begin{aligned} f_1(x)=\theta (x)q(x)-\theta (-x)\bar{q}(x), \end{aligned}$$
(4.4)

where q(x) and \(\bar{q}(x)\) are the usual quark and anti-quark parton densities in the infinite momentum frame. A similar interpretation holds for \(g_1\) and \(h_1\).

At \(z\rightarrow 0\) the operators turn to local operators. The matrix elements of local operator can be parameterized in terms of the corresponding charges. This implies the existence of exact relations relations among the first moments of PDFs. In the present case the important relations are

$$\begin{aligned} \int _{-1}^1 dx g_1(x)= & {} \int _{-1}^1 dx g_T(x),\nonumber \\ \int _{-1}^1 dx h_1(x)= & {} \int _{-1}^1 dx h_L(x), \end{aligned}$$
(4.5)

and they are another form of the Burkhardt-Cottingham sum rule [40].

In order to proceed with the matching, we need also a parameterization of off light-cone collinear functions. In general, the parameterization of matrix elements off light-cone does not coincide with the parameterization of light-cone matrix elements, which is given in Eqs. (4.14.24.3). However, on and off light-cone parameterizations can be related to each other order by order in the expansion over \(y^2\) (where y is the distance between quark fields), see e.g. discussion in [41]. Such relations up to linear terms in y are presented in Appendix B. Using the off-light-cone parameterization of Eqs. (B1B2B3) we derive the matrix elements of the first terms in the small-b OPE in EqS. (3.193.20). We find

$$\begin{aligned}&n_\alpha g_T^{\mu \nu }\lim _{y\rightarrow zn}\frac{\partial }{\partial y^\nu }\langle P,S|O^{\gamma ^\alpha }(y)|P,S\rangle =0, \end{aligned}$$
(4.6)
$$\begin{aligned}&n_\alpha g_T^{\mu \nu }\lim _{y\rightarrow zn}\frac{\partial }{\partial x^\nu }\langle P,S|O^{\gamma ^\alpha \gamma ^5}(y)|P,S\rangle \nonumber \\&\quad = 2s_T^\mu M\int du e^{2ix zp^+} \frac{g_1(x)-g_T(x)}{z}, \end{aligned}$$
(4.7)
$$\begin{aligned}&n_\gamma g_T^{\alpha \beta } g_T^{\mu \nu }\lim _{y\rightarrow zn}\frac{\partial }{\partial y^\nu }\langle P,S|O^{i\sigma ^{\beta \gamma }\gamma ^5}(y)|P,S\rangle \nonumber \\&\quad = 2\lambda M g_T^{\mu \alpha }\int dx e^{2ix zp^+} \frac{h_1(x)-h_L(x)}{z}. \end{aligned}$$
(4.8)

Moreover these expressions depend on the particular off-light-cone parameterization that is used. In any case, the functions \(g_T\) and \(h_L\) are not independent, and must be expressed in terms of distributions with definite geometrical twist. Such a re-expression is also dependent on the parameterization. In the final result all (intermediate and off-light-cone) parameterization dependence cancels, and the result is uniquely defined using definite twist distributions.

4.2 Parameterization of quark-gluon-quark correlators

The parameterization of matrix elements of a three-point operator has the following general structure

$$\begin{aligned}&\langle P,S|\mathcal {T}^\mu _\Gamma (z_1,z_2,z_3)|P,S\rangle \nonumber \\&\quad =\sum _{i}t_\Gamma ^{i;\mu ...}(P,S,n,g_T,\epsilon _T)\nonumber \\&\qquad \times \int {[}dx]e^{-ip^+(x_1z_1+x_2z_2+x_3z_3)}F(x_1,x_2,x_3), \end{aligned}$$
(4.9)

where the integration measure is [39]

$$\begin{aligned} {[}dx]= & {} dx_1dx_2dx_3\delta (x_1+x_2+x_3),\nonumber \\&-1<x_{1},x_2,x_3<1. \end{aligned}$$
(4.10)

In the rest of the paper we use the tilde notation for Fourier images of the functions

$$\begin{aligned} \tilde{F}(z_1,z_2,z_3)= & {} \int [dx]e^{-ip^+(x_1z_1+x_2z_2+x_3z_3)}F(x_1,x_2,x_3)\ . \nonumber \\ \end{aligned}$$
(4.11)

In Eq. (4.9) we have introduced a tensor t built of \(P^\mu \), (single entry of)\(S^\mu \), \(n^\mu \), \(g_T^{\mu \nu }\) and \(\epsilon _T^{\mu \nu }\), and their scalar products, such that it preserves the permutation symmetry of indices on left-hand side, and it is invariant under rescaling \(z\rightarrow \alpha z\). Such a tensor contains significant number of terms, which can be restricted by discrete symmetries, such as parity, time-reversal and charge-conjugation (which can be replaced by hermiticity due to CPT theorem). The parity invariance results into a relation among the terms of Eq. (4.9)

$$\begin{aligned}&t_\Gamma ^{i;\mu ...}(P,S,n,g_T,\epsilon _T)F(x_1,x_2,x_3)\nonumber \\&\quad = \eta _{P}^{\Gamma ;\mu ...} t_\Gamma ^{i;\mu ...}(\bar{P},s_T,\bar{n},g_T,-\epsilon _T)F(x_1,x_2,x_3), \end{aligned}$$
(4.12)

where the bar denotes the parity transformation of a vector \(\bar{v}^\mu =v_\mu \), and \(\eta _P^{\Gamma ,\mu ...}\) is the sign factor that appears in the parity transformation of the operator \(\mathcal {P}\mathcal {T}_\Gamma \mathcal {P}^\dagger =\eta _P^\Gamma \mathcal {T}_\Gamma \). The time reversal transformation results into

$$\begin{aligned}&t_\Gamma ^{i;\mu ...}(P,S,n,g_T,\epsilon _T)F(x_1,x_2,x_3)\nonumber \\&\quad = \eta _{T}^{\Gamma ;\mu ...} t_\Gamma ^{i;\mu ...}(\bar{P},-s_T,-\bar{n},-g_T,-\epsilon _T)F(-x_3,-x_2,-x_1),\nonumber \\ \end{aligned}$$
(4.13)

where \(\eta _T^{\Gamma ,\mu ...}\) is the sign factor that appears in the time-reversal transformation of the operator \(T\mathcal {T}_\Gamma T^\dagger =\eta _T^\Gamma \mathcal {T}_\Gamma \). In contrast to the two-point functions the time-reversal symmetry does not restrict the number of tensor structures \(t_i\), because the functions on left- and right-hand sides of Eq. (4.13) are of different arguments. Additionally one has the hermiticity relation which gives

$$\begin{aligned} \eta ^\Gamma _H F^*(-x_1,-x_2,-x_3)=F(x_3,x_2,x_1), \end{aligned}$$
(4.14)

where \(\eta _H\) is sign of hermitian conjugation of the operator \((\mathcal {T}_\Gamma )^\dagger =\eta _H^\Gamma \mathcal {T}_\Gamma \) (here we expect that the tensors t are real). Together the time-reversal (4.13) and hermiticity (4.14) relations dictates the complex and symmetry properties of the functions F.

In general the number of tensors t is very large. However, for the current work we need only the tensors which are non-zero if open indices are transverse, and the rest of indices are contracted with \(n^\mu \). In other words, we require the tensor structure of collinear twist-3. We find four such functions

$$\begin{aligned}&\langle P,S|\mathcal {T}^\mu _{\gamma ^+}|P,S\rangle \nonumber \\&\quad =2(p^+)^2 \tilde{s}^\mu _T M \int {[}dx]e^{-ip^+(x_1z_1+x_2z_2+x_3z_3)}T(x_1,x_2,x_3),\nonumber \\ \end{aligned}$$
(4.15)
$$\begin{aligned}&\langle P,S|\mathcal {T}^\mu _{\gamma ^+\gamma ^5}|P,S\rangle \nonumber \\&\quad =2i(p^+)^2 s^\mu _T M \int [dx]e^{-ip^+(x_1z_1+x_2z_2+x_3z_3)}\Delta T(x_1,x_2,x_3),\nonumber \\ \end{aligned}$$
(4.16)
$$\begin{aligned}&\langle P,S|\mathcal {T}^\mu _{i\sigma ^{\alpha +} \gamma ^5}|P,S\rangle \nonumber \\&\quad =2(p^+)^2 \epsilon _T^{\mu \alpha } M \int {[}dx]e^{-ip^+(x_1z_1+x_2z_2+x_3z_3)}\delta T_\epsilon (x_1,x_2,x_3)\nonumber \\&\qquad +2i(p^+)^2 \lambda g_T^{\mu \alpha } M \int {[}dx]e^{-ip^+(x_1z_1+x_2z_2+x_3z_3)}\delta T_g(x_1,x_2,x_3).\nonumber \\ \end{aligned}$$
(4.17)

Here, the factors M are set to have dimensionless three-point PDFs T. The definition of distributions T and \(\Delta T\) coincidesFootnote 3 with the definition used in [38], up to a factor M. The comparison to ETQSFootnote 4 functions (here we compare to definitions in Eq. (12) and Eq. (21) of [42]) gives

$$\begin{aligned} \widetilde{\mathcal {T}}_{q,F}(x,x+x_2)= & {} M T(-x-x_2,x_2,x),\nonumber \\ \widetilde{\mathcal {T}}_{\Delta q,F}(x,x+x_2)= & {} M \Delta T(-x-x_2,x_2,x). \end{aligned}$$
(4.18)

The distribution T are real dimensionless functions. According to Eq. (4.13) they obey the following symmetry properties

$$\begin{aligned} T(x_1,x_2,x_3)= & {} T(-x_3,-x_2,-x_1), \end{aligned}$$
(4.19)
$$\begin{aligned} \Delta T(x_1,x_2,x_3)= & {} -\Delta T(-x_3,-x_2,-x_1), \end{aligned}$$
(4.20)
$$\begin{aligned} \delta T_\epsilon (x_1,x_2,x_3)= & {} \delta T_\epsilon (-x_3,-x_2,-x_1), \end{aligned}$$
(4.21)
$$\begin{aligned} \delta T_g(x_1,x_2,x_3)= & {} -\delta T_g(-x_3,-x_2,-x_1). \end{aligned}$$
(4.22)

The Fourier transform of these distributions obey the same symmetry properties. These four functions are the only genuine twist-3 distributions in the quark sector.

It appears very convenient to introduce the following integral combinations,

$$\begin{aligned}&T^{(n)}(x)\nonumber \\&\quad =\int \frac{[dx]}{x_2^n}\left( \delta (x-x_3)+(-1)^n \delta (x+x_1)\right) T(x_1,x_2,x_3),\nonumber \\ \end{aligned}$$
(4.23)
$$\begin{aligned}&\Delta T^{(n)}(x)\nonumber \\&\quad =\int \frac{[dx]}{x_2^n}\left( \delta (x-x_3)-(-1)^n \delta (x+x_1)\right) \Delta T(x_1,x_2,x_3), \nonumber \\\end{aligned}$$
(4.24)
$$\begin{aligned}&\delta T_\epsilon ^{(n)}(x)\nonumber \\&\quad =\int \frac{[dx]}{x_2^n}\left( \delta (x-x_3)+(-1)^n \delta (x+x_1)\right) \delta T_\epsilon (x_1,x_2,x_3),\nonumber \\ \end{aligned}$$
(4.25)
$$\begin{aligned}&\delta T_g^{(n)}(x)\nonumber \\&\quad =\int \frac{[dx]}{x_2^n}\left( \delta (x-x_3)-(-1)^n \delta (x+x_1)\right) \delta T_g(x_1,x_2,x_3). \end{aligned}$$
(4.26)

The one-variable functions \(T^{(n)}\), \(\Delta T^{(n)}\) and \(\delta T^{(n)}\) are in some aspects similar to the usual PDFs. For example, they have zero boundary conditions,

$$\begin{aligned} T^{(n)}(\pm 1)= & {} 0, \quad \Delta T^{(n)}(\pm 1)=0,\nonumber \\ \delta T_\epsilon ^{(n)}(\pm 1)= & {} 0,\quad \delta T_g^{(n)}(\pm 1)=0. \end{aligned}$$
(4.27)

In the following, we intensively use the functions in Eqs. (4.234.26), since they naturally arise and describe the worm-gear functions and allow a simplification of formulas.

5 Leading matching of TMD distributions

In this section we assemble the result for the leading matching of TMD distributions up to terms linear in \({\varvec{b}}\). For this purpose we need to evaluate the matrix element of the operators in Eqs. (3.193.20) using the parameterizations introduced in the previous section. Here we should take into account the decomposition of collinear twist-3 distributions over the distributions with definite geometrical twist. In the following subsections we consider each gamma-structure individually, and discuss the features of its evaluation. For convenience we also collect the final results in Sect. 7.

5.1 Vector operator

We start with the study of the vector operator, i.e. with \(\Gamma =\gamma ^+\), in the DY kinematics. Taking the forward matrix element of the operator relation in Eq. (3.19) we obtain

$$\begin{aligned}&\langle P,S|\mathcal {U}_{\text {DY}}^{\gamma ^+}\bigg (z,\frac{{\varvec{b}}}{2}\bigg )|P,S\rangle \nonumber \\&\quad =2p^+ \int dx e^{2ix zp^+}f_1(x)+2(p^+)^2 M \tilde{s}^\mu \frac{b_\mu }{2} \nonumber \\&\qquad \times \bigg [-i\int _{-1}^1 dv vz \tilde{T}(z,vz,-z)\nonumber \\&\qquad +i\left( \int _{-\infty }^z+\int _{-\infty }^{-z}\right) d\tau \tilde{T}(z,\tau ,-z)\bigg ]+O({\varvec{b}}^2), \end{aligned}$$
(5.1)

where the contribution of the two-point correlator vanishes in accordance to Eq. (4.6).

The function \(T(z,vz -z)\) is symmetric in v due to the symmetry relation in Eq. (4.19). Therefore, the anti-symmetric integral, which is the first in the square brackets of Eq. (5.1), vanishes,

$$\begin{aligned} \int _{-1}^{1} dv\,vz\,\tilde{T}(z,vz,-z)= & {} 0. \end{aligned}$$
(5.2)

In this way, the contributions linear in b are represented by a single entry, namely, by the last term of Eq. (5.1). Using the reflection of coordinates in Eq. (4.19) we present it as

$$\begin{aligned} \left( \int _{-\infty }^z+\int _{-\infty }^{-z}\right) d\tau \tilde{T}(z,\tau ,-z) = \int ^{\infty }_{-\infty } d\tau \tilde{T}(z,\tau ,-z).\nonumber \\ \end{aligned}$$
(5.3)

Taking into account these simplifications we find

$$\begin{aligned}&\langle P,S|\mathcal {U}_{\text {DY}}^{\gamma ^+} \left( z,\frac{{\varvec{b}}}{2}\right) |P,S\rangle \nonumber \\&\quad = 2p^+ \int dx e^{2ix zp^+}f_1(x)+i(p^+)^2 M \tilde{s}^\mu b_\mu \nonumber \\&\qquad \times \int _{-\infty }^{\infty }d\tau \tilde{T}(z,\tau ,-z)+O({\varvec{b}}^2). \end{aligned}$$
(5.4)

In the case of SIDIS kinematic the operators are given by Eq. (3.20). Applying the same procedure we find

$$\begin{aligned}&\langle P,S|\mathcal {U}_{\text {DIS}}^{\gamma ^+} \left( z,\frac{{\varvec{b}}}{2}\right) |P,S\rangle \nonumber \\&\quad =2p^+ \int dx e^{2ix zp^+}f_1(x)-i(p^+)^2 M \tilde{s}^\mu b_\mu \nonumber \\&\qquad \times \int _{-\infty }^{\infty }d\tau \tilde{T}(z,\tau ,-z)+O({\varvec{b}}^2), \end{aligned}$$
(5.5)

where we have used

$$\begin{aligned} \left( \int ^{\infty }_z+\int ^{\infty }_{-z}\right) d\tau \tilde{T}(z,\tau ,-z)= & {} \int ^{\infty }_{-\infty } d\tau \tilde{T}(z,\tau ,-z).\nonumber \\ \end{aligned}$$
(5.6)

The only difference between Drell-Yan, Eq. (5.4) and SIDIS, Eq. (5.5) cases is the sign of the linear term. It corresponds to the famous process dependence of the Sivers function [19].

The TMD distribution is obtained by Fourier transformation over the light-cone distance, Eq. (3.4). Performing it we obtain

$$\begin{aligned}&\text {(DY)}\qquad \Phi ^{[\gamma ^+]}_{q\leftarrow h}(x,{\varvec{b}})\nonumber \\&\quad = f_1(x)+i b_\mu \tilde{s}_T^\mu M\, \pi T(-x,0,x)+O({\varvec{b}}^2), \end{aligned}$$
(5.7)
$$\begin{aligned}&\text {(SIDIS)}\qquad \Phi ^{[\gamma ^+]}_{q\leftarrow h}(x,{\varvec{b}})\nonumber \\&\quad =f_1(x)-i b_\mu \tilde{s}_T^\mu M\, \pi T(-x,0,x)+O({\varvec{b}}^2). \end{aligned}$$
(5.8)

Here we have used,

$$\begin{aligned} \int _{-\infty }^\infty \frac{dz}{2\pi }\int _{-\infty }^\infty d\tau e^{-2ixzp^+}\tilde{T}(z,\tau ,-z)=\frac{\pi }{(p^+)^2}T(-x,0,x).\nonumber \\ \end{aligned}$$
(5.9)

These expressions represent the leading matching of vector TMD distribution. Comparing it to the parameterization in Eq. (2.18) we find the matching of individual functions. Naturally, the unpolarized TMDPDF matches the unpolarized PDF, \(f_1(x,{\varvec{b}})=f_1(x)+O({\varvec{b}}^2)\). The Sivers function matching is process dependent and it reads

$$\begin{aligned} \text {(DY)}\qquad f_{1T}^\perp (x,{\varvec{b}})= & {} \pi T(-x,0,x)+O({\varvec{b}}^2), \end{aligned}$$
(5.10)
$$\begin{aligned} \text {(SIDIS)}\qquad f_{1T}^\perp (x,{\varvec{b}})= & {} -\pi T(-x,0,x)+O({\varvec{b}}^2). \end{aligned}$$
(5.11)

Note, that the correction term is proportional to \({\varvec{b}}^2\), and therefore, generically, contains twist-5 functions (and twist-4 functions for unpolarized distribution).

These expression, albeit in the different form, are well-known. In the two-point notation for ETQS function (4.18), the central value of three-point function \(T(-x,0,x)\) corresponds to the diagonal value \(\widetilde{\mathcal {T}}_{q,F}(x,x)\). Therefore, we can compare (5.105.11) to the expressions given in literature, where certain momentum space moments are calculated. Using the transformation rules presented in Appendix A, one can check that

$$\begin{aligned} \int d^2{\varvec{p}}_T \frac{{\varvec{p}}_T^2}{M^2}f_{1T}^\perp (x,p_T)= & {} 2\pi T(-x,0,x),\nonumber \\ \int d^2{\varvec{p}}_T e^{-i ({\varvec{b}} {\varvec{p}}_T)} \frac{p_T^\alpha }{M} f_{1T}^\perp (x,p_T)= & {} i\pi b^\alpha T(-x,0,x).\nonumber \\ \end{aligned}$$
(5.12)

Here the sign is given for the DY case, and should be changed for the SIDIS case. To our best understandingFootnote 5 these expression coincide with ones presented in [16, 17, 43, 44].

5.2 Axial operator

Taking the forward matrix element of the operator in Eq. (3.19) with \(\Gamma =\gamma ^+\gamma ^5\), we obtain

$$\begin{aligned}&\langle P,S|\mathcal {U}_{\text {DY}}^{\gamma ^+\gamma ^5} \left( z,\frac{{\varvec{b}}}{2}\right) |P,S\rangle \nonumber \\&\quad =2\lambda p^+ \int dx e^{2ix zp^+}g_1(x)\nonumber \\&\qquad +2M s^\mu _T \frac{b_\mu }{2} \bigg [\int du e^{2iuzp^+}\frac{g_1(u)-g_T(u)}{z}\nonumber \\&\qquad +(p^+)^2\int _{-1}^1 dv vz \Delta \tilde{T}(z,vz,-z) \nonumber \\&\qquad -(p^+)^2\left( \int _{-\infty }^z+\int _{-\infty }^{-z}\right) d\tau \Delta \tilde{T}(z,\tau ,-z)\bigg ]+O({\varvec{b}}^2),\nonumber \\ \end{aligned}$$
(5.13)

where we have used the parameterizations in Eqs. (4.24.16) and the relation of Eq. (4.7).

To proceed further we take the inverse Fourier transform. We have observed that these integrals naturally enter into the moments of the three-point functions, which are defined in Eqs. (4.234.26). Moreover, it is convenient to present them as a Mellin convolution. Using these tricks we find

$$\begin{aligned}&\int \frac{dz}{2\pi }e^{-2ixzp^+} \int _{-1}^1 dv vz \Delta \tilde{T}(z,vz,-z)\nonumber \\&\quad =\frac{i}{(p^+)^2}\left[ \frac{\Delta T^{(1)}(x)}{2}+\int _{-1}^1 du \int _0^1 dy \, u\Delta T^{(2)}(u) \delta (x-yu)\right] , \nonumber \\ \end{aligned}$$
(5.14)
$$\begin{aligned}&\int \frac{dz}{2\pi }e^{-2ix zp^+}\int du e^{2iu zp^+} \frac{g_1(u)-g_T(u)}{z}\nonumber \\&\quad =i\int _{-1}^{1} du \int _0^1 dy\, u(g_1(u)-g_T(u))\delta (x-uy). \end{aligned}$$
(5.15)

The last integral in Eq. (5.13) over the process-dependent term does not vanish,

$$\begin{aligned}&\int \frac{dz}{2\pi }e^{-2ixzp^+} \left( \int _{-\infty }^z+\int _{-\infty }^{-z}\right) d\tau \Delta \tilde{T}(z,\tau ,-z)\nonumber \\&\quad =\frac{i}{(p^+)^2}\frac{\Delta T^{(1)}(x)}{2}, \end{aligned}$$
(5.16)
$$\begin{aligned}&\int \frac{dz}{2\pi }e^{-2ixzp^+} \left( \int ^{\infty }_z+\int ^{\infty }_{-z}\right) d\tau \Delta \tilde{T}(z,\tau ,-z)\nonumber \\&\quad =\frac{-i}{(p^+)^2}\frac{\Delta T^{(1)}(x)}{2}, \end{aligned}$$
(5.17)

where we have used the assumption that the integrand goes to zero at infinity. The sign difference between these integrals, is compensated by the common sign difference in the operators for DY, Eq. (3.19) and SIDIS, Eq.  (3.20) kinematics. Therefore, the contribution of seemingly process-dependent terms is the same for both operators. It is exactly compensated by the contribution of Eq. (5.14), and thus the function \(\Delta T^{(1)}\) drops out of calculation.

Combining all together we obtain the same result for DY and SIDIS kinematics, which is

$$\begin{aligned}&\Phi ^{[\gamma ^+\gamma ^5]}(x,{\varvec{b}}) \nonumber \\&\quad =\lambda g_1(x)+i b_\mu s_T^\mu M\int _{-1}^{1} du \nonumber \\&\qquad \times \int _0^1 dy \delta (x-uy)u \left( g_1(u)-g_T(u)+\Delta T^{(2)}(u)\right) .\nonumber \\ \end{aligned}$$
(5.18)

Comparing to parameterizations in Eq. (2.19) we find that the matching for the helicity TMD distribution \(g_{1L}(x,{\varvec{b}})=g_1(x)+O({\varvec{b}}^2)\), and for the worm-gear-T distribution is

$$\begin{aligned} g_{1T}(x,{\varvec{b}})= & {} \int _{-1}^{1} du \int _0^1 dy \delta (x-uy)\nonumber \\&\times u\left( g_1(u)-g_T(u)+\Delta T^{(2)}(u)\right) +O({\varvec{b}}^2).\nonumber \\ \end{aligned}$$
(5.19)

The expression (5.19) is not the final one, because the function \(g_T\) can be rewritten via functions of definite twist,

$$\begin{aligned} g_T(x)= & {} \int _0^1 dy \int _{-1}^1 du \delta (x-yu)\nonumber \\&\times \left[ g_1(u) +\frac{T^{(1)}(u)-\Delta T^{(1)}(u)-\varepsilon _+ h_1(u)}{2u}(1-\delta (\bar{y}))\right. \nonumber \\&\left. +\Delta T^{(2)}(u)\right] , \end{aligned}$$
(5.20)

where \(\varepsilon _+=2m/M\) with m being the mass of a quark and \(\bar{y}=1-y\). The derivation of this decomposition is given in the Appendix C.1. It is straightforward to check that it obeys the Burkhardt-Cottingham sum rule Eq. (4.5). Inserting the function \(g_T\) into Eq. (5.19) and using the associativity of Mellin transformation (see also Eq. (C6)) we obtain

$$\begin{aligned} g_{1T}(x,{\varvec{b}})= & {} x\int _{-1}^{1} du \int _0^1 dy \delta (x-uy)\nonumber \\&\times \left( g_1(u)+\Delta T^{(2)}(u)\right. \nonumber \\&\left. +\frac{T^{(1)}(u)-\Delta T^{(1)}(u)-\varepsilon _+ h_1(u)}{2u}\right) +O({\varvec{b}}^2).\nonumber \\ \end{aligned}$$
(5.21)

This is the final form of the matching of the worm-gear function to the twist-2 and twist-3 functions. The Mellin convolution, which is presented in Eq. (5.21) by \(\delta \)-function, can be explicitly integrated. It gives the following representation

$$\begin{aligned} g_{1T}(x,{\varvec{b}})= & {} x \int _x^1 \frac{du}{u}\nonumber \\&\times \left( g_1(u)+\Delta T^{(2)}(u)\right. \nonumber \\&\left. +\frac{T^{(1)}(u) -\Delta T^{(1)}(u) -\varepsilon _+ h_1(u)}{2u}\right) \nonumber \\&+O({\varvec{b}}^2),\qquad \text {for }x>0, \end{aligned}$$
(5.22)
$$\begin{aligned} g_{1T}(x,{\varvec{b}})= & {} x \int ^{x}_{-1} \frac{du}{|u|}\nonumber \\&\times \left( g_1(u)+\Delta T^{(2)}(u)\right. \nonumber \\&\left. +\frac{T^{(1)}(u)-\Delta T^{(1)}(u)-\varepsilon _+ h_1(u)}{2u}\right) \nonumber \\&+O({\varvec{b}}^2),\qquad \text {for }x<0, \end{aligned}$$
(5.23)

The obtained result can be compared to the first transverse momentum moments of TMD distribution derived in Ref. [45] (see Eq. (47)), and agrees with it.

5.3 Tensor operator

The matrix element of the tensor operator in Eq. (3.19) (i.e. with \(\Gamma =i\sigma ^{\alpha +}_T\gamma ^5\)) has a more complicated form

$$\begin{aligned}&\langle P,S|\mathcal {U}_{\text {DY}}^{i\sigma ^{\alpha +}_T\gamma ^5} \left( z,\frac{{\varvec{b}}}{2}\right) |P,S\rangle \nonumber \\&\quad =2s_T^\mu p^+ \int dx e^{2ix zp^+}h_1(x) +2M \frac{b_\mu }{2} \nonumber \\&\qquad \times \bigg [\lambda g_T^{\mu \alpha }\int du e^{2iuzp^+}\frac{h_1(u)-h_L(u)}{z}\nonumber \\&\qquad +(p^+)^2\int _{-1}^1 dv vz \left( \lambda g_T^{\mu \alpha }\delta \tilde{T}_g(z,vz,-z)\right. \nonumber \\&\qquad \left. -i\epsilon _T^{\mu \alpha }\delta \tilde{T}_\epsilon (z,vz,-z)\right) \nonumber \\&\qquad -(p^+)^2\left( \int _{-\infty }^z+\int _{-\infty }^{-z}\right) d\tau \left( \lambda g_T^{\mu \alpha }\delta \tilde{T}_g(z,\tau ,-z)\right. \nonumber \\&\qquad \left. -i\epsilon _T^{\mu \alpha }\delta \tilde{T}_\epsilon (z,\tau ,-z)\right) \bigg ]+O({\varvec{b}}^2), \end{aligned}$$
(5.24)

where we have used the parameterizations of Eqs. (4.34.17) and the relation (4.8). Its structure repeats the structure discussed during evaluations of the vector operator (for terms proportional to \(\delta T_\epsilon \)) and axial operator (for terms proportional to \(\delta T_g\)). Therefore, we skip the discussion on the Fourier integrals and write down the final expression for the matching of transversally polarized TMD distribution. We obtain (compare to Eqs. (5.85.7) and (5.18))

$$\begin{aligned}&\Phi ^{[i\sigma ^{\alpha +}\gamma ^5]}(x,{\varvec{b}})\nonumber \\&\quad =s_T^\alpha h_1(x)\pm ib_\mu \epsilon ^{\mu \alpha }_T \pi \delta T_\epsilon (-x,0,x) \nonumber \\&\qquad +i\lambda b^\alpha M \int _{-1}^1 du \int _0^1 dy \delta (x-uy)\nonumber \\&\qquad \times u\left( h_1(u)-h_L(u)+\delta T_g^{(2)}(u)\right) +O({\varvec{b}}^2), \end{aligned}$$
(5.25)

where the upper sign should be taken for the DY kinematics, and lower for the SIDIS kinematics.

Comparing Eq. (5.25) with the parameterization of Eq. (2.20) we obtain the matching of individual TMD distributions. The transversity distribution is \(h_1(x,{\varvec{b}})=h_1(x)+O({\varvec{b}}^2)\). The Boer-Mulders functions depends on the underling process and reads

$$\begin{aligned} \text {(DY)}\qquad h_1^\perp (x,{\varvec{b}})= & {} -\pi \delta T_\epsilon (-x,0,x)+O({\varvec{b}}^2), \end{aligned}$$
(5.26)
$$\begin{aligned} \text {(SIDIS)}\qquad h_1^\perp (x,{\varvec{b}})= & {} \pi \delta T_\epsilon (-x,0,x)+O({\varvec{b}}^2). \end{aligned}$$
(5.27)

The worm-gear L function is independent on the process and has the expression

$$\begin{aligned} h_{1L}^\perp (x,{\varvec{b}})= & {} -\int _{-1}^1 du \int _0^1 dy \delta (x-uy)\nonumber \\&\times u\left( h_1(u)-h_L(u)+\delta T_g^{(2)}(u)\right) +O({\varvec{b}}^2).\nonumber \\ \end{aligned}$$
(5.28)

The pretzelosity distribution has no matching at this level of accuracy, despite the fact that the matrix element over free quarks is non-zero at \({\varvec{b}}^2\rightarrow 0\) [12]. It is expected that the first non-zero contribution to the pretzelosity is of twist-4.

As in the case of the worm-gear T function, the expression for the worm-gear L function should be rewritten via a definite twist function. The derivation of function \(h_L\) is given Appendix C.2. It reads

$$\begin{aligned} h_L(x)= & {} -\int _0^1 dy \int _{-1}^1 du \delta (x-y u)\nonumber \\&\times \left[ 2y(h_1(u)+\delta T^{(2)}_g(u))\right. \nonumber \\&\left. -\left( \frac{\delta T^{(1)}_g(u)}{u} +\frac{\varepsilon _+g_1(u)}{2u}\right) (2y-\delta (\bar{y}))\right] . \end{aligned}$$
(5.29)

Consequently, the worm-gear L function is

$$\begin{aligned}&h_{1L}^\perp (x,{\varvec{b}})\nonumber \\&\quad =-x \int _{-1}^1 du \int _0^1 dy \delta (x-uy)\nonumber \\&\qquad \times y\left( h_1(u)+\delta T_g^{(2)}(u)-\frac{\delta T_g^{(1)}(u)}{u}-\frac{\varepsilon _+g_1(u)}{2u}\right) +O({\varvec{b}}^2).\nonumber \\ \end{aligned}$$
(5.30)

Let us point that the expressions for worm-gear L and worm-gear T functions have very similar structure, compare Eq. (5.30) to Eq. (5.21). The main difference is the factor y that appears in Eq. (5.21). The integral over \(\delta \)-function could be evaluated with the result

$$\begin{aligned} h_{1L}^\perp (x,{\varvec{b}})= & {} -x^2 \int _x^1 \frac{du}{u^2}\left( h_1(u)+\delta T_g^{(2)}(u)-\frac{\delta T_g^{(1)}(u)}{u}+\frac{\varepsilon _+g_1(u)}{2u}\right) \nonumber \\&+O({\varvec{b}}^2),\qquad \text {for }x>0, \end{aligned}$$
(5.31)
$$\begin{aligned} h_{1L}^\perp (x,{\varvec{b}})= & {} x^2 \int _{-1}^x \frac{du}{u^2}\left( h_1(u)+\delta T_g^{(2)}(u)-\frac{\delta T_g^{(1)}(u)}{u}+\frac{\varepsilon _+g_1(u)}{2u}\right) \nonumber \\&+O({\varvec{b}}^2),\qquad \text {for }x<0. \end{aligned}$$
(5.32)

This expression can be compared to the first transverse momentum moment of TMD distribution derived in Ref. [45] (see Eq.(49) in this reference). We find that our expression agrees with the result of Ref. [45].

Table 1 The summary of the quark TMD distributions and their leading matching at small-b

6 Mellin moments of worm-gear functions

The final expressions for worm-gear function, as well as, all intermediate expressions are naturally expressed via Mellin convolutions. This fact suggests a simple form for the Mellin moments of the worm-gear functions, which we present in this section.

First of all, let us point that functions \(T^{(n)}\), \(\Delta T^{(n)}\) and \(\delta T^{(n)}\) defined in Eq. (4.234.26) obey certain relations which simplify in the algebra of Mellin moments. The Mellin moment of \(T^{(n)}\) is

$$\begin{aligned} \int _{-1}^1 dx x^k T^{(n)}(x)=\left\{ \begin{array}{ll} 0,&{} \quad k=0,\quad n \text { odd},\\ 0,&{} \quad k\text { odd},\\ T^{(n,k)}, &{} \quad \text {otherwise}, \end{array}\right. \end{aligned}$$
(6.1)

which follows from the symmetry properties Eq. (4.19). The same relations hold for the function \(\delta T_\epsilon \). For antisymmetric functions \(\Delta T\) and \(\delta T_g\) we have

$$\begin{aligned} \int _{-1}^1 dx x^k \Delta T^{(n)}(x)=\left\{ \begin{array}{ll} 0,&{} \quad k=0,\quad n \text { even},\\ 0,&{} \quad k\text { even},\\ \Delta T^{(n,k)}, &{}\quad \text {otherwise}, \end{array}\right. \end{aligned}$$
(6.2)

which follows from symmetry property Eq. (4.20).

Evaluating the Mellin moments of the worm-gear functions in Eqs. (5.20) and (5.29) we find simple expression

$$\begin{aligned} g^{(k)}_{1T}({\varvec{b}})= & {} \int _{-1}^1 dx x^k g_{1T}(x,{\varvec{b}})\nonumber \\= & {} \frac{1}{k+2}\left\{ \begin{array}{ll} g_1^{(k+1)}-\frac{1}{2}\Delta T^{(1,k)},&{} k\text { odd},\\ g_1^{(k+1)}+\Delta T^{(2,k+1)}+\frac{1}{2}T^{(1,k)},&{} k\text { even}, \end{array}\right. \nonumber \\ \end{aligned}$$
(6.3)
$$\begin{aligned} h^{\perp (k)}_{1L}({\varvec{b}})= & {} \int _{-1}^1 dx x^k h_{1L}(x,{\varvec{b}})\nonumber \\= & {} \frac{-1}{k+3}\left\{ \begin{array}{ll} h_1^{(k+1)}-\delta T_g^{(1,k)},&{} k\text { odd}, \\ h_1^{(k+1)}+\delta T_g^{(2,k+1)},&{} k\text { even}, \end{array}\right. \end{aligned}$$
(6.4)

where we omit the quark mass correction. The peculiar feature of this expressions is that functions with odd and even index n enter different moments independently.

Such relations can be important for lattice studies of TMD distributions, where only Mellin moments of distributions can be evaluated. For example, in Ref. [46] the lattice calculation of the first moment of \(g_{1T}\) is performed. It has been found that

$$\begin{aligned} \frac{g_{1T}^{(0)}({\varvec{b}}\simeq 0.34)}{f_1^{(0)}({\varvec{b}}\simeq 0.34)}\Bigg |_{[42]}\approx 0.2. \end{aligned}$$
(6.5)

The calculation has been done for the isovector combination of operators \(q=u-d\). Here, the scales of TMD distributions are not clear since the translation rules between lattice scales and TMD evolution scales are not elaborated so far. Nonetheless, the evolution factors for both distributions are the same, and up to the first order of approximation the scale dependence of Eq. (6.5) can be omitted. In Ref. [46] it is shown that b-dependence of the ratio in Eq. (6.5) is very weak. In particular, the value at \(b\simeq 0.46\) practically coincides with Eq. (6.5). This suggests that the small-b expansion is a good approximation.

We can estimate the ratio in Eq. (6.5) from our calculation. Using Eq. (6.3) we find

$$\begin{aligned} \frac{g_{1T}^{(0)}({\varvec{b}})}{f_1^{(0)}({\varvec{b}})}=\frac{g_1^{(1)}+\Delta T^{(2,1)}}{2 f^{(0)}_1}+O(\alpha _s)+O({\varvec{b}}^2). \end{aligned}$$
(6.6)

Assuming that the contribution of \(\Delta T^{(2,1)}\) is small, i.e. in the Wandzura-Wilczek approximation, we find this ratio is \(\sim 0.13\), which at this level of comparison is a good agreement.

7 Conclusion

In this work we have evaluated the operator product expansion for the quark TMD operators up to linear in \({\varvec{b}}\) terms. This order of expansion includes the majority of the polarized distributions. The summary of the matching relations is presented in Table 1. The main result of this study is the leading matching of Sivers, Boer-Mulders and worm gear function. We resume all the matching here for simplicity

$$\begin{aligned} f_{1T}^\perp (x,{\varvec{b}})= & {} \pm \pi T(-x,0,x)+O({\varvec{b}}^2), \end{aligned}$$
(7.1)
$$\begin{aligned} g_{1T}(x,{\varvec{b}})= & {} x\int _{-1}^{1} du \int _0^1 dy \delta (x-uy)\nonumber \\&\times \left( g_1(u)+\Delta T^{(2)}(u)+\frac{T^{(1)}(u)-\Delta T^{(1)}(u)}{2u}\right) +O({\varvec{b}}^2), \nonumber \\ \end{aligned}$$
(7.2)
$$\begin{aligned} h_1^\perp (x,{\varvec{b}})= & {} \mp \pi \delta T_\epsilon (-x,0,x)+O({\varvec{b}}^2), \end{aligned}$$
(7.3)
$$\begin{aligned} h_{1L}^\perp (x,{\varvec{b}})= & {} -x \int _{-1}^1 du \int _0^1 dy \delta (x-uy)\nonumber \\&\times y\left( h_1(u)+\delta T_g^{(2)}(u)-\frac{\delta T_g^{(1)}(u)}{u}\right) +O({\varvec{b}}^2), \end{aligned}$$
(7.4)

where the upper sign corresponds to the Drell-Yan operator, and lower sign corresponds to the SIDIS operator. The functions \(g_1\) and \(h_1\) are helicity and transversity PDFs. The functions T are collinear distributions of twist-3. Their definition is given in (4.154.17, 4.234.26).

The expressions presented here are only the leading order perturbative QCD terms. The sub-leading terms include the power corrections in \({\varvec{b}}^2\) and perturbative corrections. The perturbative corrections can be accumulated into the coefficient functions. For the distributions that match solely to twist-2 PDF these coefficient functions are already known at higher perturbative orders; next-to-next-to-leading order (NNLO) for unpolarized [10, 11, 18] and transversity [12] distributions) and at NLO for helicity distribution [9, 47]. The polarized TMD distributions such as Sivers, Boer-Mulders, Collins and worm-gear functions, matches the twist-3 and twist-2 distributions. For these distributions, the coefficient functions are known only at LO and presented here altogether.

We find that the results obtained by us agree with expressions that we have found in the literature (as far as we can trace necessary definitions of various components). However, there are several essential differences since all known expressions are given in momentum space and they are presented in terms of certain integrals of TMD distribution over \(p_T\). This fact complicates the comparison since such integrals are not well-defined within perturbation theory, and require some regularization procedure. In contrast, our calculation is done directly in position space, and in this aspect, it represents a complete novelty. Another important distinctive fact of our work is that our calculation is based solely on the definition of TMD operators, whereas the majority of higher-twist calculations are based on the evaluation of particular cross-sections with successive re-interpretation in terms of TMDs. The only known example that we have found of a direct calculation is the one presented in Ref. [45], where the leading matching for worm-gear functions is calculated. Finally, we consider all TMD distributions on the same foot and in the same framework which provides a consistent relative normalization of all distributions improving their comprehension.