The corresponding contribution to trace part of
$W(\alpha_q,\beta_q,x_\perp)$ takes the form

\begin{eqnarray}
g^{mn}W_{mn}^{\rm 6+7th}(q)~=~{s/2\over(2\pi)^4N_c}\int\!d^4x ~e^{-iqx} g^{mn}\cheW_{mn}^{\rm 6+7th}(x_\perp)
\label{w67trace}\\
=~-{2\over (N_c^2-1)\alpha_q\beta_qs}\!\int\! d^2k_\perp (k,q-k)_\perp[(j_2-\barj_2)(\alpha_q,k_\perp)(j^\star_2-\barj^\star_2)
(\beta_q,(q-k)_\perp)+{\rm c.c.}]
\nonumber
\end{eqnarray}

which agrees with Eq. (6.2) from Ref. [28] after
replacements $j_2=j_2^{\rm tw 3}-i\tilj_2^{\rm tw 3}$ and $\barj_2=j_1^{\rm tw 3}+i\tilj_1^{\rm tw 3}$.
It should be noted that the difference between $j_1$ and $j_2$ in traceless vs trace part is due to difference in formulas
(<ref>) and (<ref>).

Thus, the result is the sum of Eqs. (<ref>) and

\begin{eqnarray}
W_{\mu\nu}^{\rm 6+7th}(q)~=~{s/2\over(2\pi)^4N_c}\int\!d^4x ~e^{-iqx} \cheW_{mn}^{\rm 6+7th}(x_\perp)
\nonumber\\
=~{1\over (N_c^2-1)\alpha_q\beta_qs}\!\int\! d^2k_\perp\Big([(j_1-\barj_1)(\alpha_q,k_\perp)(j^\star_1-\barj^\star_1)
(\beta_q,(q-k)_\perp)+{\rm c.c.}]
\nonumber\\
\times~[k_\mu(q-k)_\nu+\mu\leftrightarrow\nu+g^\perp_{\mu\nu}(k,q-k)_\perp]
\nonumber\\
-~g^\perp_{\mu\nu} (k,q-k)_\perp[(j_2-\barj_2)(\alpha_q,k_\perp)(j^\star_2-\barj^\star_2)
(\beta_q,(q-k)_\perp)+{\rm c.c.}]
\label{w67otvet}
\end{eqnarray}

As we mentioned in the Introduction, in this paper we we will take into account only leading  and sub-leading terms in $N_c$
and leave the ${1\over N_c^2}$ corrections discussed above for future publications.

Finally, as proved in Appendix <ref>, we can neglect contributions proportional to the product of quark and gluon TMDs.

§ RESULTS AND ESTIMATES

§.§ Results

Assembling Eqs. (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) we get
the result for $W_{\mu\nu}(q)$ that consists of two parts:

\begin{equation}
\label{result}
\end{equation}

The first, gauge-invariant,  part is given by

\begin{eqnarray}
\nonumber\\
W^{1F}_{\mu\nu}(q)~=~\sum_f e_f^2W^{fF}_{\mu\nu}(q),~~~~
W^{fF}_{\mu\nu}(q)~=~{1\over N_c}\!\int\!d^2k_\perp F^f(q,k_\perp)\calW^F_{\mu\nu}(q,k_\perp),
\nonumber\\
W^{1H}_{\mu\nu}(q)~=~\sum_f e_f^2W^{fH}_{\mu\nu}(q),~~~~W^{fH}_{\mu\nu}(q)~=~ {1\over N_c}\!\int\!d^2k_\perp H^f(q,k_\perp)\calW^H_{\mu\nu}(q,k_\perp)
\label{resultginv}
\end{eqnarray}

where $F^f$ and $H^f$ are given by Eq. (<ref>) and

\begin{eqnarray}
\calW^{F}_{\mu\nu}(q,k_\perp)~=~
-g_{\mu\nu}^\perp +{1\over Q_\parallel^2}(q^\parallel_\mu q^\perp_\nu+q^\parallel_\nu q^\perp_\mu)
+{q_\perp^2\over Q_\parallel^4}q^\parallel_\mu q^\parallel_\nu+{\tilq_\mu\tilq_\nu\over Q_\parallel^2}[q_\perp^2-4(k,q-k)_\perp]
\nonumber\\
-~\Big[{\tilq_\mu\over  Q_\parallel^2}\Big(g^\perp_{\nu i}-{q^\parallel_\nu q_i\over Q_\parallel^2}\Big)(q-2k)_\perp^i
\label{WF}
\end{eqnarray}

\begin{eqnarray}
\label{WH}\\
+2{\tilq_\mu\tilq_\nu-q^\parallel_\mu q^\parallel_\nu \over Q_\parallel^4}k_\perp^2(q-k)_\perp^2
\nonumber\\
-~\Big({q^\parallel_\mu\over Q_\parallel^2}\big[k_\perp^2(q-k)^\perp_\nu+k^\perp_\nu(q-k)_\perp^2\big]
+~{\tilq_\mu\over Q_\parallel^2}\big[k_\perp^2(q-k)^\perp_\nu-k^\perp_\nu(q-k)_\perp^2\big]+\mu\leftrightarrow\nu\Big)
\nonumber\\
-~{\tilq_\mu\tilq_\nu+q^\parallel_\mu q^\parallel_\nu \over Q_\parallel^4}\big[q_\perp^2-2(k,q-k)_\perp\big](k,q-k)_\perp
-~{q^\parallel_\mu\tilq_\nu+ \tilq_\mu q^\parallel_\nu \over Q_\parallel^4}(2k-q,q)_\perp(k,q-k)_\perp
\nonumber
\end{eqnarray}

where $q^\parallel_\mu\equiv \alpha_qp_1+\beta_qp_2$ and $\tilq_\mu\equiv \alpha_qp_1-\beta_qp_2$. These are the same expressions as in Eq. (<ref>) if
one identifies $x_A$ with $\alpha_q$ and $x_B$ with $\beta_q$ and neglects $O\big({m^2\over s}\big)$ terms
in $p_A$ and $p_B$ and $O\big({q_\perp^4\over Q^4}\big)$ corrections due to difference between $Q^2$ and $Q_\parallel^2$.
It is easy to see that $q^\mu W^F_{\mu\nu}~=~0$ and  $q^\mu W^H_{\mu\nu}~=~0$. Note that
$q^\mu W^F_{\mu\nu}$ and  $q^\mu W^H_{\mu\nu}$ are exactly zero without any ${q_\perp^2\over Q^2}$ corrections.
This is similar to usual “forward” DIS, but different from off-forward DVCS where the cancellations of right-hand sides of Ward identities
involve infinite towers of twists [42, 43, 44]

The second part is

\begin{eqnarray}
=~{1\over N_c}\sum_f e_f^2{1\over Q^2}\!\int d^2k_\perp
\bigg[ {1\over m^2}\big\{[k^\perp_\mu(q-k)^\perp_\nu+\mu\leftrightarrow\nu](k,q-k)_\perp
\nonumber\\
-(q-k_\perp)^2k^\perp_\mu k^\perp_\nu+g_{\mu\nu}^\perp(k,q-k)_\perp^2
-g_{\mu\nu}^\perp k_\perp^2(q-k_\perp)^2\big]\big\}H_A^f(q,k_\perp)
\nonumber\\
+~{N_c\over N_c^2-1}\Big([k^\perp_\mu(q-k)^\perp_\nu+\mu\leftrightarrow\nu+g^\perp_{\mu\nu}(k,q-k)_\perp]J_1^f(q,k_\perp)
\nonumber\\
-~g^\perp_{\mu\nu} (k,q-k)_\perp J_2^f(q,k_\perp)\Big)~+~O\big({1\over N_c^2}\big)\bigg]~+~O\big({Q_\perp^4\over Q^4}\big)
\label{resultnoninv}
\end{eqnarray}

where $H_A$ is given by Eq. (<ref>) and

\begin{eqnarray}
J_1^f(q,k_\perp)~=~(j_1-\barj_1)(\alpha_q,k_\perp)(j^\star_1-\barj^\star_1)(\beta_q,(q-k)_\perp)+{\rm c.c.}
\nonumber\\
(\beta_q,(q-k)_\perp)+{\rm c.c.}
\label{W2}
\end{eqnarray}

These terms are not gauge invariant: $q^\mu W^2_{\mu\nu}(q)\neq 0$. The reason is that gauge invariance is restored after adding terms like ${m_\perp^2\over Q^2}\times{\rm Eq. (\ref{pc})}$  which we do not calculate in this paper. Indeed, for example,

\begin{eqnarray}
q^\mu W^2_{\mu\nu}(q)~\sim~{q^\perp_\nu q_\perp^2\over\alpha_q\beta_qs}~~~~
{\rm and}~~~q^\mu \times~{p_2^\mu q_\perp^\nu q_\perp^2\over \alpha_q^2\beta_qs^2}
~=~{q^\perp_\nu q_\perp^2\over\alpha_q\beta_qs}
\end{eqnarray}

They are of the same order so one should expect that gauge invariance is restored after calculation of the terms $\sim{p_2^\mu q_\perp^\nu q_\perp^2\over \alpha_q^2\beta_qs^2}$ which are beyond the scope of this paper.
For the same reason we see that all structures in Eq. (<ref>) except
${g_{\mu\nu}^\perp q_\perp^2\over\alpha_q\beta_qs}$
and ${q^\perp_\mu q^\perp_\nu \over\alpha_q\beta_qs}$ are determined by leading-twist TMDs $f_1$ and $h_1^\perp$.

Sometimes it is convenient to represent hadronic tensor in transverse coordinate space. Introducing

\begin{eqnarray}
\begin{array}{c}
\\
\barf(\alpha,b_\perp)
\end{array}\bigg\}~=~\int\!{d^2k_\perp\over 4\pi^2}~e^{i(k,b)_\perp}
\bigg\{\begin{array}{c}
\\
\barf(\alpha,k_\perp)
\end{array}
\end{eqnarray}

(and similarly for target TMDs) we get

\begin{eqnarray}
\label{w1fcoord}\\
=~4\pi^2\sum_f e_f^2
\Big\{-g_{\mu\nu}^\perp +{i\over Q_\parallel^2}(q^\parallel_\mu\partial^\perp_\nu+q^\parallel_\nu\partial^\perp_\mu)
-{q^\parallel_\mu q^\parallel_\nu+\tilq_\mu\tilq_\nu\over Q_\parallel^2}\partial_\perp^2\Big]f_1\barf_1
-~{4\tilq_\mu\tilq_\nu\over Q_\parallel^2}(\partial_if_1^f)(\partial^i \barf_1^f)
\nonumber\\
-~\Big[{\tilq_\mu\over  Q_\parallel^2}\Big(\delta_\nu^i-{q^\parallel_\nu \over Q_\parallel^2}i\partial^i\Big)
\nonumber
\end{eqnarray}

$f_1=f_1(\alpha_q,b_\perp)$, $\barf_1\equiv\barf_1(\beta_q,b_\perp)$ everywhere except
$f\leftrightarrow\barf$ terms where it is opposite (the question about rapidity cutoffs for TMDs will be addressed in Sect.

Similarly, we can write down $W^H$ contribution in coordinate space. For future use, however, it is convenient to define Fourier transform in a slightly different way. Introduce  $h_i(k_\perp)\equiv k_i h_1^\perp(k), ~\barh_i(k_\perp)\equiv k_i \barh_1^\perp(k)$ and

\begin{eqnarray}
\begin{array}{c}
\\
\barh_i(\alpha,b_\perp)
\end{array}\bigg\}~=~\int\!{d^2k_\perp\over 4\pi^2}~e^{i(k,b)_\perp}
\bigg\{\begin{array}{c}
\\
\barh_i(\alpha,k_\perp)
\end{array}
\label{defhi}
\end{eqnarray}

then $W^{1H}_{\mu\nu}$  can be represented as

\begin{eqnarray}
m^2W_{\mu\nu}^{1H}(\alpha_q,\beta_q,b_\perp)~=~4\pi^2\sum_{\rm f} e_f^2\Big(g_{\mu\nu}^\perp h^f_j\barh^{fj}-h^f_\mu\barh^f_\nu-h_\nu\barh^f_\mu
+~2{q^\parallel_\mu q^\parallel_\nu -\tilq_\mu\tilq_\nu\over Q_\parallel^4}\partial^i h_i^f\partial^j\barh_j^f
\nonumber\\
+~{q^\parallel_\mu\over Q_\parallel^2}\big[(i\partial^ih^f_i)\barh^f_\nu+h^f_\nu i\partial^i\barh^f_i\big]
+~{\tilq_\mu\over Q_\parallel^2}\big[(i\partial^ih^f_i)\barh^f_\nu-h^f_\nu i\partial^i\barh^f_i\big]~+~\mu\leftrightarrow\nu\Big)
-~{\tilq_\mu\tilq_\nu+q^\parallel_\mu q^\parallel_\nu \over Q_\parallel^4}
\nonumber\\
\times~\big[\partial_\perp^2(h^f_i\barh_f^i)+\partial^i h^f_j\partial^i\barh^f_j\big]
-~{q^\parallel_\mu\tilq_\nu+ \tilq_\mu q^\parallel_\nu \over Q_\parallel^4}\partial_i[(h^f_j\partial\barh_f^j-\barh^f_j\partial h_f^j)]~+~h\leftrightarrow\barh
\label{w1hcoord}
\end{eqnarray}

$h_i=h_i(\alpha_q,k_\perp)$ and $\barh_i\equiv\barh_i(\beta_q,(q-k)_\perp)$ everywhere except
$h\leftrightarrow\barh$ terms where it is opposite, cf. Eq. (<ref>).

§.§ Four Lorentz structures of hadronic tensor

The four Lorentz structures of hadronic tensor in Collins-Soper frame are given by Eq. (<ref>) where ($Q_\perp\equiv |q_\perp|$)

\begin{eqnarray}
Z~=~{\tilq\over Q_\parallel}~=~{1\over Q_\parallel}(\alpha_qp_1-\beta_qp_2),~~~~~~~
X~=~\Big[{Q_\perp\over Q_\parallel Q}(\alpha_qp_1+\beta_qp_2)+{Q_\parallel\over Q_\perp Q}q_\perp\Big]
\end{eqnarray}

such that $q\cdot X=q\cdot Z=X\cdot Z=0$ and $X^2=Z^2=-1$.

First, let us check the structure corresponding to the total cross section of DY pair production. From Eq. (<ref>) we get

\begin{eqnarray}
~=~-{2\over N_c}\sum_fe_f^2\!\int\! d^2k_\perp\Big\{
\big[1-2{(k,q-k)_\perp\over Q^2}\big]
\label{wmumu}\\
+~2{k_\perp^2(q-k)_\perp^2\over m^2_NQ^2}H^f(q,k_\perp)
+{N_c\over N_c^2-1}(k,q-k)_\perp J_2^f(q,k_\perp)\Big\}
\Big[1+O\big({1\over N_c^2}\big)\Big]+O\big({q_\perp^4\over Q^4}\big)
\nonumber
\end{eqnarray}

which agrees with Eq. (6.2) from Ref. [28]. This equation gives the sum of structures  $W_\mu^\mu~=~-(2W_T+W_L)$.

§.§.§ $W_L$

The easiest structure to get is $W_L$.  Multiplying Eq. (<ref>) by $Z_\mu Z_\nu$ and comparing to Eq. (<ref>) we get

\begin{eqnarray}
W_L(q)~=~Z^\mu Z^\nu W^1_{\mu\nu}(q)~
\sum_fe_f^2{1\over Q^2N_c}\!\int\! d k_\perp\Big\{ (q-2k)_\perp^2F^f(q,k_\perp)
\label{wl}\\
+~{1\over m^2}\big(2k_\perp^2(q-k)_\perp^2-[k_\perp^2+(q-k)_\perp^2](k,q-k)_\perp\big)
\Big[1+O\big({q_\perp^2\over Q^2}\big)+O\big({1\over N_c^2}\big)\Big]
\nonumber
\end{eqnarray}

Thus, one may say that  $W_L$ is known at LHC energies at $q_\perp^2\ll Q^2$ as far as $f_1$ and $h_1^\perp$ are known.

§.§.§ $W_\Delta$

Using formula
$q_\perp^\mu Z^\nu W_{\mu\nu}~=~(X\cdot q)_\perp W_\Delta~=~-{Q_\parallel Q_\perp\over Q}W_\Delta$
we get

\begin{eqnarray}
W_\Delta~=~{Q\over Q_\parallel^2Q_\perp N_c}\sum_fe_f^2
\!\int\!d^2k_\perp\Big\{(q,q-2k)_\perp F^f(q,k_\perp)
\nonumber\\
-~(q,q-2k){(k,q-k)_\perp\over m^2}H^f(q,k_\perp)\Big\}\Big[1+O\big({q_\perp^2\over Q^2}\big)+O\big({1\over N_c^2}\big)\Big]
\label{wd}
\end{eqnarray}

Again, we see that $W_\Delta$ is expressed via $f_1$ and $h_1^\perp$ with great accuracy.

§.§.§ $W_T$

Next, from Eqs. (<ref>), (<ref>) and $W_\mu^\mu~=~-(2W_T+W_L)$ one easily obtains

\begin{eqnarray}
W_T(q)~=~{1\over N_c}\sum_fe_f^2\!\int\! d^2k_\perp\Big\{
\big[1-{q_\perp^2\over 2Q_\parallel^2}\big]F^f(q,k_\perp)
\nonumber\\
+~{1\over 2m^2Q_\parallel^2}\big(2k_\perp^2(q-k)_\perp^2+[k_\perp^2+(q-k)_\perp^2](k,q-k)_\perp\big)H^f(q,k_\perp)
\nonumber\\
+~{N_c\over N_c^2-1}(k,q-k)_\perp J_2^f(q,k_\perp)\Big\}
\Big[1+O\big({1\over N_c^2}\big)\Big]+O\big({q_\perp^4\over Q^4}\big)
\label{wt}
\end{eqnarray}

§.§.§ $W_\dd$

Finally, the easiest way to pick out $W_\dd$ is to multiply $W_{\mu\nu}$ by $q^\perp_\mu q^\perp_\nu/q_\perp^2$. One obtains from Eq. (<ref>)
${q^\perp_\mu q^\perp_\nu\over q_\perp^2} W^{\mu\nu}(q)={Q_\parallel^2\over Q^2}(W_T-W_\dd)$. On the other hand, from Eqs. (<ref>) and (<ref>)
one gets

\begin{eqnarray}
{q_\perp^\mu q_\perp^\nu\over q_\perp^2} W^1_{\mu\nu}(q)~
\nonumber\\
&&\hspace{-1mm}=~{1\over N_c}\sum_f e_f^2\!\int\!d^2k_\perp \Big\{F^f(q,k_\perp)
+~\Big[{(k,q-k)_\perp\over m^2}-{2(q,k)_\perp(q,q-k)_\perp\over m^2q_\perp^2}\Big] H^f(q,k_\perp)\Big\}
\nonumber\\
\label{vkladotw1}
\end{eqnarray}

and from Eq. (<ref>)

\begin{eqnarray}
{q_\perp^\mu q_\perp^\nu\over q_\perp^2} W^2_{\mu\nu}(q)~
=~{1\over N_c}\sum_f e_f^2{1\over Q_\parallel^2q_\perp^2}\!\int d^2k_\perp\Big\{ {1\over m^2}
\big\{2(q,k)_\perp(q,q-k)_\perp(k,q-k)_\perp
\label{vkladotw2}\\
+q_\perp^2 k_\perp^2(q-k_\perp)^2\big]\big\}H_A^f(q,k_\perp)
\nonumber\\
+~{N_c\over N_c^2-1}\Big([2(q,k)_\perp(q,q-k)_\perp-q_\perp^2(k,q-k)_\perp]J_1^f(q,k_\perp)
+~q_\perp^2 (k,q-k)_\perp J_2^f(q,k_\perp)\Big)\Big\}.
\nonumber
\end{eqnarray}

Thus, we get

\begin{eqnarray}
W_\dd~=~W_T-{Q^2\over Q_\parallel^2}{q^\perp_\mu q_\perp^\nu\over q_\perp^2}W_{\mu\nu}~=~
{1\over N_c}\sum_f e_f^2\!\int\!d^2k_\perp \Big\{{q_\perp^2\over 2Q_\parallel^2}F^f(q,k_\perp)
\label{wdd}\\
+~\Big({2(q,k)_\perp(q,q-k)_\perp\over q_\perp^2}-(k,q-k)_\perp
+{1\over Q^2}\big\{k_\perp^2(q-k)_\perp^2
\nonumber\\
+~\half (k,q-k)_\perp[k_\perp^2+(q-k)_\perp^2]+q_\perp^2(k,q-k)_\perp-2(q,k)(q,q-k)_\perp\big\}
\Big){1\over m^2}H^f(q,k_\perp)
\nonumber\\
+~{1\over m^2Q_\parallel^2}\Big({k_\perp^2\over q_\perp^2}(q,q-k)_\perp^2+{(q-k_\perp)^2\over q_\perp^2}(q,k)_\perp^2
- k_\perp^2(q-k_\perp)^2+(k,q-k)_\perp^2
\nonumber\\
-~2(q,k)_\perp(q,q-k)_\perp{(k,q-k)_\perp\over q_\perp^2}\Big)H_A^f(q,k_\perp)
\nonumber\\
+~{N_c\over N_c^2-1}{1\over Q_\parallel^2}\Big[(k,q-k)_\perp-{2(q,k)_\perp(q,q-k)_\perp\over q_\perp^2}\Big]J_1^f(q,k_\perp)
~+~O\big({1\over N_c^2}\big)\Big\}+O\big({q_\perp^4\over Q^4}\big)
\nonumber
\end{eqnarray}

This is the only function which has a $O\big({1\over Q^2}\big)$, leading-$N_c$ contribution proportional to twist-three TMD $H_A$ not related to leading-twist TMDs
by equations of motion. The  functions $W_T,W_L$, and $W_\Delta$ do not have such contributions (although they have such contributions at the ${1\over N_c}$ level).

§.§ Estimates of $W_i(q)$ at $q_\perp^2\gg m^2$

§.§.§ Order-of-magnitude estimates

Following the analysis in Ref. [28], let us estimate the relative strength of Lorentz structures $W_i$ at $q_\perp^2\gg m^2$.
First, we assume that ${1\over N_c}$ is a good parameter and leave only terms leading in $N_c$.
Second,  at $q_\perp^2\gg m^2$ we probe the perturbative tails of TMD's $f_1\sim{1\over k_\perp^2}$ and
$h_1^\perp\sim{1\over k_\perp^4}$ [45].
So, as long as  $Q^2\gg q_\perp^2\gg m^2$ we can approximate

\begin{equation}
\hspace{-0mm}
f_1(\alpha_z,k_\perp^2)~\simeq~{f(\alpha_q)\over k_\perp^2},~~h_1^\perp(\alpha_q,k_\perp)~\simeq~{m^2_Nh(\alpha_q)\over k_\perp^4} ,
~~~~~~\barf_1\simeq{\barf(\alpha_q)\over k_\perp^2},~~\barh_1^\perp\simeq{m^2_N\barh(\alpha_q)\over k_\perp^4}
\label{fhestimate1}
\end{equation}

(up to logarithmic corrections). Similarly, for the target we can use the estimate

\begin{equation}
f_1(\beta_z,k_\perp^2)~\simeq~{f(\beta_z)\over k_\perp^2},~~h_1^\perp(\beta_z,k_\perp^2)~\simeq~{m^2_N h(\beta_z)\over k_\perp^4} ,
~~~~~~\barf_1\simeq{\barf(\beta_z)\over k_\perp^2},~~\barh_1^\perp\simeq{m^2_N\barh(\beta_z)\over k_\perp^4}
\label{fhestimate2}
\end{equation}

as long as $k_\perp^2\ll Q^2$.  Thus, we get an estimate

\begin{eqnarray}
F^f(q,k_\perp)~\simeq~{F^f(\alpha_q,\beta_q)\over k_\perp^2(q-k)_\perp^2},
\nonumber\\
H^f(q,k_\perp)~\simeq~m^4{H^f(\alpha_q,\beta_q)\over k_\perp^4(q-k)_\perp^4},
\label{fhestimate}
\end{eqnarray}

Note that due to the “positivity constraint” [46]

\begin{equation}
h_1^\perp(x,k_\perp^2)\leq {m\over |k_\perp|}f_1^\perp(x,k_\perp^2)~~~
\label{poscons}
\end{equation}

we can safely assume that the functions $f(x)$ and $h(x)$ defined in
Eqs. (<ref>) and (<ref>) are of the same order of magnitude.
Moreover,  both theoretical [47] and phenomenological [48, 49] analysis  indicate  that $h_1^\perp$ is several times smaller than $f_1$
so in numerical estimates we will disregard the contribution of $h_1^\perp$.

§.§.§ Power corrections for total DY cross section

Substituting the above approximations  to Eq. (<ref>) we get
the following estimate  of the strength of power corrections for total DY cross section [28]

\begin{eqnarray}
~=~-{2\over N_c}\sum e_f^2\!\int\! d^2k_\perp\Big\{
\Big[1-2{(k,q-k)_\perp\over Q^2}\Big]F^f(q,k_\perp)
+~2{k_\perp^2(q-k)_\perp^2\over m^2_NQ^2}H^f(q,k_\perp)  \Big\}
\nonumber\\
\simeq~-{2\over N_c}\sum e_f^2\!\int\! d^2k_\perp
\Big\{\Big[1-2{(k,q-k)_\perp\over Q^2}\Big]
{F^f(\alpha_q,\beta_q)\over k_\perp^2(q-k)_\perp^2}
+{2m^2\over Q^2}{H^f(\alpha_q,\beta_q)\over k_\perp^2(q-k)_\perp^2}\Big\}
\nonumber\\
\simeq~-{2\over N_c}\sum e_f^2\!\int\! d^2k_\perp
\Big[1-2{(k,q-k)_\perp\over Q^2}\Big]
{F^f(\alpha_q,\beta_q)\over k_\perp^2(q-k)_\perp^2}
\label{dyx1}
\end{eqnarray}

where we used estimates (<ref>) and the fact that $(k,q-k)_\perp\sim q_\perp^2\gg m^2$.
Thus, the relative weight of the leading term and power correction is determined by the factor
$1-2{(k,q-k)_\perp\over Q^2}$. Due to Eqs. (<ref>) and (<ref>), the integrals over $k_\perp$ are logarithmic and should be cut from below
by $m^2_N$ and from above by $Q^2$ so we get an estimate

\begin{eqnarray}
\!\int\! d^2k_\perp{1\over k_\perp^2(q-k)_\perp^2}~\simeq~{2\pi\over q_\perp^2}\ln{q_\perp^2\over m^2},~~~~~~
\!\int\! d^2k_\perp{(k,q-k)_\perp\over k_\perp^2(q-k)_\perp^2}~\simeq~-\pi\ln{Q^2\over q_\perp^2}
\label{loginteg1}
\end{eqnarray}

where we assumed that the first integral is determined by the logarithmical region $q_\perp^2\gg k_\perp^2\gg m^2_N$ and the second by $Q^2\gg k_\perp^2\gg q_\perp^2$.
Taking these integrals to Eq. (<ref>)  one obtains

\begin{equation}
\hspace{-0mm}
~=~-{4\pi \over N_c}\sum e_f^2\Big[{1\over q_\perp^2}\ln{q_\perp^2\over m^2_N}+{1\over Q^2}\ln{Q^2\over q_\perp^2}\Big]
\label{wmumuestimate}
\end{equation}

By this estimate, the power correction reaches the level of few percent at $q_\perp\sim {Q\over 4}$.

§.§.§ Power corrections for $W_T$

Let us now consider estimates described in  Sect. <ref>  for $W_T$ given by Eq. (<ref>).
At large $N_c$, we can omit the third line so

\begin{eqnarray}
W_T(q)~=~{1\over N_c}\sum_fe_f^2\!\int\! \dhd^2k_\perp\Big[1-{q_\perp^2\over 2Q^2}\Big]F^f(q,k_\perp)
\\
+~{1\over 2m^2Q^2}\big(2k_\perp^2(q-k)_\perp^2+[k_\perp^2+(q-k)_\perp^2](k,q-k)_\perp\big)H^f(q,k_\perp)
\nonumber\\
\simeq~{1\over N_c}\sum_fe_f^2\!\int\! \dhd^2k_\perp\Big\{\Big[1-{q_\perp^2\over 2Q^2}\Big]
{F^f(\alpha_q,\beta_q)\over k_\perp^2(q-k)_\perp^2}
\nonumber\\
+~{m^2\over Q^2}\Big(1+[k_\perp^2+(q-k)_\perp^2]{(k,q-k)_\perp\over 2k_\perp^2(q-k)_\perp^2}\Big)
{H^f(\alpha_q,\beta_q)\over k_\perp^2(q-k)_\perp^2}\Big\}
\label{wtestimat}
\end{eqnarray}

Again, due to $q_\perp^2\gg m^2$ the second term in braces can be neglected and we get

\begin{eqnarray}
W_T(q)~\simeq~{2\pi \over N_c}\Big[{1\over q_\perp^2}-{1\over 2Q^2}\Big]\ln{q_\perp^2\over m^2}\sum e_f^2F^f(\alpha_q,\beta_q)
\label{wtestimate}
\end{eqnarray}

Thus, for $W_T$ the power correction reaches 10% level at $q_\perp\sim {Q\over 2}$.

§.§.§ Estimate of $W_L$

Again, using estimates from Sect. <ref> one obtains

\begin{eqnarray}
\sum_fe_f^2{1\over Q^2N_c}\!\int\! {d k_\perp\over k_\perp^2(q-k)_\perp^2}
\Big\{ (q-2k)_\perp^2F^f(\alpha_q,\beta_q)
\nonumber\\
+~m^2\Big(2-[k_\perp^2+(q-k)_\perp^2]{(k,q-k)_\perp\over k_\perp^2(q-k)_\perp^2}\Big)
\label{wlestimat}
\end{eqnarray}

which gives approximately

\begin{eqnarray}
{2\pi\over Q^2N_c}\Big[\ln{q_\perp^2\over m^2_N}+2\ln{Q^2\over q_\perp^2}\Big]\sum_fe_f^2F^f(\alpha_q,\beta_q)
\label{wlestimate}
\end{eqnarray}

in agreement with Eqs. (<ref>) and (<ref>).
The estimate of the ratio of $W_L/W_T$ is

\begin{eqnarray}
{W_L(q)\over W_T(q)}~\simeq~
{q_\perp^2\over Q^2}\Big[1+2{\ln Q^2/q_\perp^2\over \ln q_\perp^2/m^2}\Big]
\label{wltestimate}
\end{eqnarray}

§.§.§ Magnitude of $W_\Delta$

It is easy to see that $W_\Delta$ vanishes if one uses the estimates (<ref>) and (<ref>). Indeed, with these formulas
$F(q,k_\perp)$ and $H(q,k_\perp)$ are symmetric under replacement $k_\perp\leftrightarrow (q-k)_\perp$ whereas $(q,q-2k)_\perp$ in
the integrand in Eq. (<ref>) is antisymmetric. Moreover, this vanishing of $W_\Delta$ will occur for any factorizable model of TMDs $f_1$ and $h_1^\perp$:
if $f_1(\alpha,k_\perp)~=~f(\alpha)\phi(k_\perp)$ and $h_1^\perp(\alpha,k_\perp)~=~h(\alpha)\psi(k_\perp)$ the integral (<ref>)
On the other hand, $W_\Delta$ is only $\sim {Q_\perp\over Q}W_T$ so without better knowledge of TMDs it is impossible
to tell whether $W_\Delta$ is smaller or bigger than, say,  $W_L$. Also, if the parameter $\alpha_q{Q\over Q_\perp}$ is not negligible, to compare
$W_\Delta$ and  $W_L$ one needs to take into account $O(\alpha_q)$
corrections to $W_\Delta$ defined by TMDs other than $f_1$ and $h_1^\perp$.

§.§.§ Estimate of $W_\dd$

Let us consider the relative weight of the terms in the r.h.s. of Eq. (<ref>). As we mentioned, we assume that ${1\over N_c}$ is a valid
small parameter so we can omit the last $J_1$ term. Also,  it is natural to assume that $H_A^f(q,k_\perp)$ is of the same order of magnitude as
$H^f(q,k_\perp)$ and,  since the term with $H_A$ is a power correction, it is not unreasonable to neglect this term in the first approximation.
Using now estimates (<ref>)  and the integrals

\begin{eqnarray}
\int\! d^2k {k_i(q-k)_j\over k_\perp^2(q-k)_\perp^4}\theta(k_\perp^2-m^2)\theta((q-k)_\perp^2-m^2)
~\simeq~{\pi\over 2q_\perp^2}\Big(g^\perp_{ik}+2{q_iq_k\over q_\perp^2}\Big)\ln{q_\perp^2\over m^2}
\nonumber\\
\int\! d^2k {k_i(q-k)_j\over k_\perp^4(q-k)_\perp^4}\theta(k_\perp^2-m^2)\theta((q-k)_\perp^2-m^2)
~\simeq~{\pi\over q_\perp^4}\Big(g^\perp_{ik}+4{q_iq_k\over q_\perp^2}\Big)\ln{q_\perp^2\over m^2}
\label{loginteg2}
\end{eqnarray}

one gets an estimate

\begin{eqnarray}
{1\over N_c}\sum_f e_f^2\!\int\!d^2k_\perp \Big\{{q_\perp^2\over 2Q_\parallel^2}F^f(q,k_\perp)
\nonumber\\
+~\Big({2(q,k)_\perp(q,q-k)_\perp\over q_\perp^2}-(k,q-k)_\perp
+{1\over Q^2}\big[k_\perp^2(q-k)_\perp^2
\nonumber\\
+~\half (k,q-k)_\perp[k_\perp^2+(q-k)_\perp^2]+q_\perp^2(k,q-k)_\perp-2(q,k)(q,q-k)_\perp\big]
\Big){1\over m^2}H^f(q,k_\perp)
\nonumber\\
\simeq~{\pi\over Q^2N_c}\ln{q_\perp^2\over m^2}\sum_f e_f^2\Big[F^f(\alpha_q,\beta_q)
+{4m^2Q^2\over q_\perp^4}\Big(1-{q_\perp^2\over 2Q^2}\Big)H^f(\alpha_q,\beta_q)\Big]
\label{wddest}
\end{eqnarray}

§.§.§ Lam-Tung relation

It is easy to see that if one neglects $H$ in Eq. (<ref>) the ratio of $W_L$ and $2W_\dd$ is approximately

\begin{eqnarray}
{W_L\over 2W_\dd}~\simeq~1+2{\ln Q^2/q_\perp^2\over\ln{q_\perp^2/m^2}}
\label{ltungestim}
\end{eqnarray}

It seems like the Lam-Tung relation works better if we move closer to the domain of collinear factorization  $Q^2\sim Q_\perp^2\gg m^2$.

§.§.§ Estimates of asymmetries

The differential cross section of DY process is parametrized as

\begin{equation}
\hspace{1mm}
\Big({d\sigma\over d^4q}\Big)^{-1}{d\sigma\over d\Omega d^4q}~=~{3\over 4\pi(\lambda+3)}\big(1+\lambda\cos^2\theta+\mu\sin 2\theta\cos\phi+
{\nu\over 2}\sin^2\theta\cos 2\phi\big)
\end{equation}

where  $\Omega$
is the solid angle of the lepton in terms of its polar and azimuthal angles in the center-of-mass system of the lepton pair.
The angular coefficients $\lambda$, $\mu$, and $\nu$ can be expressed in terms of the hadronic tensor:

\begin{equation}
\hspace{1mm}
\lambda~=~{W_T-W_L\over W_T+W_L},~~~\mu~=~{W_\Delta\over W_T+W_L},~~~\nu~=~{2W_\dd\over W_T+W_L}
\label{lamunu}
\end{equation}

For an estimate, let us take $s$=8 TeV and $Q$=90 GeV so that $x_A\sim x_B\sim 0.1$ in central region of rapidity.
Although we did not include the contribution of $Z$-boson, we can compare our order-of-magnitude estimates with experimental data
at this kinematics [50, 51].
Let us take $Q_\perp$= 20 GeV so
the power corrections $\sim{Q_\perp^2\over Q^2}$ are small but sizable, of order of few per cent.  At this kinematics, we obtain

\begin{equation}
\hspace{1mm}
1-\lambda~=~2{W_L\over W_T+W_L}~\simeq~ 2{1+2{\ln Q^2/q_\perp^2\over \ln q_\perp^2/m^2}\over
{Q^2\over q_\perp^2}-\half+2{\ln Q^2/q_\perp^2\over \ln q_\perp^2/m^2}} ~\simeq~0.19
\end{equation}

from Eq. (<ref>) which agrees with estimates in Ref. [52]. Next, in our kinematics the expression
in square brackets in the r.h.s. in Eq. (<ref>)  is approximately $F+0.17 H$. Since  the Boer-Mulders function
seem to be of order of few percent of $f_1$ (see the discussion in Sect. <ref>), the term with $H$ can be safely neglected and we get

\begin{equation}
\hspace{1mm}
\nu~=~{2W_\dd\over W_T+W_L}~\simeq~{1\over {Q^2\over q_\perp^2}-\half+2{\ln Q^2/q_\perp^2\over\ln{q_\perp^2/m^2}}}~\simeq~0.05
\end{equation}

As to $\mu$ coefficient, as we mentioned, we cannot estimate it since with factorization hypothesis
for TMDs it vanishes. Reversing the argument, if $\mu$ will be measured to be much smaller than $\nu$, it will be an argument in favor
of factorization hypothesis for TMDs $f_1$ and $h_1^\perp$. Actually, there are  experiments at much lower $q_\perp$ and $Q$ $\sim$ few GeV
which indicate that $\mu$ is very small [53].

Last but not least, let us estimate Lam-Tung relation. With our approximation in the above kinematics we get

\begin{eqnarray}
{W_L\over 2W_\dd}~\simeq~1+2{\ln Q^2/q_\perp^2\over\ln{q_\perp^2/m^2}}~\simeq~2.0
\label{ltungest}
\end{eqnarray}

so it seems to be violated at this kinematics. Again, these order-of-magnitude estimates do not include the contribution to DY cross section
mediated by the $Z$-boson.

§ COEFFICIENT FUNCTIONS AND MATCHING OF RAPIDITY CUTOFFS

The result (<ref>) is a tree-level formula and to fully understand Eq. (<ref>) we
should specify the rapidity cutoffs for $f_1$'s and $h_1^\perp$'s. As we discussed in section <ref>,
the rapidity cutoff for longitudinal momenta in $f_1(\alpha_q,k_\perp)$ is $\beta\leq \sigma_p$ and for $f_1(\beta_q,k_\perp)$ $\alpha\leq\sigma_t$,
where $\sigma_p$ and $\sigma_t$ are rapidity bounds for central fields. To avoid double counting, the region where both
$\alpha<\sigma_t$ and $\beta<\sigma_p$ should give only small power corrections. This is achieved if one takes
$\sigma_p,\sigma_t\sim{Q_\perp\over\sqrt{s}}$ so power corrections from double counting are ${Q_\perp\over Q}$.
In this case, the region $\alpha_q>\alpha>\sigma_t,~\beta_q>\beta>\sigma_p$ gives Sudakov double-log

\begin{eqnarray}
~\sim ~e^{-{\alpha_sc_F\over\pi}\ln{\alpha_q\over\sigma_t}\ln{\beta_q\over\sigma_p}}
\label{doublog}
\end{eqnarray}

where the coefficient ${\alpha_sc_F\over\pi}$ is two times $\gamma_{\rm cusp}$ for quarks.
A more precise formula can be obtained from the requirement
that the product of two TMDs and the coefficient function (<ref>) does not depend on the “rapidity divides”
$\sigma_p$ and $\sigma_t$. For simplicity, let us start with the leading-twist term $\sim g_{\mu\nu}^\perp F$.
Rapidity evolution of the function $f_1(\alpha_q,b_\perp;\vigma_p)$ was found in Ref. [54]
[As noted in Ref. [54], the factor $\sim \gamma_E$ depends on the exact way to cut integrals over $\alpha$
and $\beta$. Here the factor $-\gamma_E$ corresponds to “smooth” cutoffs $e^{-{\alpha\over\sigma_t}}$ and $e^{-{\beta\over\sigma_p}}$,
see the discussion in Ref. [54]

\begin{eqnarray}
{d\over d\ln\vigma_p}f_1(\alpha_q,b_\perp;\vigma_p)
~=~{\alpha_sC_F\over \pi}\big[-\ln\alpha_q\vigma_p-\half\ln {b_\perp^2s\over 4} -\gamma_E
\nonumber\\
{d\over d\ln\vigma_t}f_1(\beta_q,b_\perp;\vigma_t)
~=~{\alpha_sC_F\over \pi}\big[-\ln\beta_q\vigma_t-\half\ln {b_\perp^2s\over 4} -\gamma_E+O(\alpha_s)\big]f_1(\beta_q,b_\perp;\vigma_t)
\label{evoleqs}
\end{eqnarray}

where $\vigma_p=\sigma_pb_\perp\sqrt{s}$, $\vigma_t=\sigma_tb_\perp\sqrt{s}$ are $b_\perp$-dependent cutoffs providing conformal invariance of the
leading-order TMD rapidity evolution (in the coordinate space) and $\gamma_E$ is Euler's constant.
Similar equation holds true for $\barf_1$ since it is obtained from the evolution of the same operator.

Looking at  Eqs. (<ref>) and (<ref>) one can guess that
the coefficient function $\sim g_\perp^{\mu\nu}$ times two TMDs $f_1$ in the coordinate space has the form

\begin{eqnarray}
\big[f_1(\alpha_q,b_\perp;\vigma_p)\barf_1(\beta_q,b_\perp;\vigma_t)
+f_1\leftrightarrow \barf_1\big]
\label{match}
\end{eqnarray}

\begin{equation}
\hspace{-0mm}
\ln\big({\alpha_qb_\perp\sqrt{s}\over \vigma_t}e^{\gamma_E}\big)
\ln\big({\beta_qb_\perp\sqrt{s}\over \vigma_p}e^{\gamma_E}\big)
+{\alpha_sc_F\over 2\pi}\ln^2\vigma_p\vigma_t}\big[
\label{M}
\end{equation}

It is easy to check that with $M$ given by Eq. (<ref>)  we have ${d\over d\vigma_p}$(r.h.s. of Eq. (<ref>)) = 0
and ${d\over d\vigma_t}$(r.h.s. of Eq. (<ref>)) = 0 so our guess (<ref>) for the coefficient function
is correct up to $O(\alpha_s)$ terms.
[The Eq. (<ref>) is obtained in the leading order in $\alpha_s$ so
the argument of coupling constant is left undetermined.
One should expect Sudakov formula with running coupling constant [55, 11] at the NLO level.

To write precise matching for other parts of $W_{\mu\nu}$ is a more complicated task. Let us start with  $W^{1F}_{\mu\nu}$ terms considered in
the next Section.

§.§ Matching for $W^{1F}$ terms

We need to multiply Eq. (<ref>)  in coordinate space  by $M(\alpha_q,\beta_q,b_\perp;\vigma_p,\vigma_t)$.
First, recall that

\begin{equation}
\hspace{0mm}
[f_1(\alpha_q,b_\perp;\vigma_p)\barf_1(\beta_q,b_\perp;\vigma_t) +f_1\leftrightarrow \barf_1\big]
\end{equation}

does not actually depend on the “rapidity divides” $\vigma_p$ and $\vigma_t$. However, the differentiation
${\partial\over\partial b_i}$ affects evolution equations (<ref>).
In this case we  modify the derivative with respect to $b_i$ as follows

\begin{eqnarray}
\tartial_i f_1(\alpha_q,b_\perp;\vigma_p)~\equiv~\big(\partial_i-{\alpha_sc_F\over \pi}{b^i\over b_\perp^2}\ln\vigma_p\big)f_1(\alpha_q,b_\perp;\vigma_p),
\nonumber\\
\tartial_i f_1(\beta_q,b_\perp;\vigma_t)~\equiv~\big(\partial_i-{\alpha_sc_F\over \pi}{b^i\over b_\perp^2}\ln\vigma_t\big)f_1(\beta_q,b_\perp;\vigma_t)
\label{tartial}
\end{eqnarray}

(and similarly for $\barf$'s)
so that the l.h.s.' of these equations satisfy Eqs. (<ref>).
[Strictly speaking, the difference between $\tartial_i f$ and $\partial_i f$ is
$\sim O(\alpha_s)$  but since our matching is correct  at single-log limit (see Eq. (<ref>)) we keep $\tartial_i f$
to avoid corrections $\sim O(\alpha_s\ln\vigma)$ in our matching formulas (<ref>) (and (<ref>) below).
Note also that $\partial_i(Mf\barf)=M(f\tartial_i\barf+\barf\tartial_i f)$.
With this definitions, one can write $W^{1F}_{\mu\nu}$ in the double-log approximation in the form

\begin{eqnarray}
W^{1F}(\alpha_q,\beta_q,b_\perp)~=~\sum_{\rm flavors} e_f^2 W^{fF}(\alpha_q,\beta_q,b_\perp)~
\label{w1fcoordm}\\
\nonumber\\
=~ 4\pi^2
\Big\{-g_{\mu\nu}^\perp +{i\over Q_\parallel^2}(q^\parallel_\mu\partial^\perp_\nu+q^\parallel_\nu\partial^\perp_\mu)
-{q^\parallel_\mu q^\parallel_\nu+\tilq_\mu\tilq_\nu\over Q_\parallel^2}\partial_\perp^2\Big]Mf_1\barf_1
-~{4\tilq_\mu\tilq_\nu\over Q_\parallel^2}M(\tartial^if_1^f)(\tartial^i \barf_1^f)
\nonumber\\
-~\Big[{\tilq_\mu\over  Q_\parallel^2}\Big(\delta_\nu^i-{q^\parallel_\nu \over Q_\parallel^2}i\partial^i\Big)
\nonumber
\end{eqnarray}

where $M=M(\alpha_q,\beta_q, b_\perp;\vigma_p,\vigma_t)$ and
$f_1=f_1(\alpha_q,b_\perp,\vigma_p)$, $\barf_1\equiv\barf_1(\beta_q,b_\perp,\vigma_t)$ everywhere except
$f\leftrightarrow\barf$ terms where it is opposite.

It is easy to check gauge invariance:  $(\alpha_q p_1^\mu+\beta_q p_2^\mu+i\partial_\perp^\mu)W_{\mu\nu}^{1F}(\alpha_q,\beta_q,b_\perp)~=~0$.

§.§ Matching for $W^{1H}$ terms

First, with our definitions (<ref>)  the  Eq. (<ref>) reads

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(x_\bu,x_\perp)\sigma_{i \bu}\psi_f(0)|A\rangle
=~{1\over m_N} h_i^f(\alpha,k_\perp)
\nonumber\\
{1\over 8\pi^3s}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(0)\sigma_{i \bu}\psi_f(x_\bu,x_\perp)|A\rangle
=~-{1\over m_N}\barh_i^f(\alpha,k_\perp)
\label{hmaelnew}
\end{eqnarray}

and similarly for the target matrix elements. With such definition, the evolution equation for
$h_i(\alpha,b_\perp)\equiv {1\over 4\pi^2}\int\! d^2k_\perp e^{i(k,x)_\perp}h_i(\alpha,k_\perp)$
is the same as Eq. (<ref>)

\begin{eqnarray}
{d\over d\ln\vigma_p}h_i(\alpha_q,b_\perp;\vigma_p)
~=~{\alpha_sC_F\over \pi}\big[-\ln\alpha_q\vigma_p-\ln {|b_\perp|\sqrt{s}\over 2} -\gamma_E\big]h_i(\alpha_q,b_\perp;\vigma_p)
\nonumber\\
{d\over d\ln\vigma_t}h_i(\beta_q,b_\perp;\vigma_t)
~=~{\alpha_sC_F\over \pi}\big[-\ln\beta_q\vigma_t-\half\ln {b_\perp^2s\over 4} -\gamma_E\big]h_i(\beta_q,b_\perp;\vigma_t)
\label{hevoleqs}
\end{eqnarray}

and similarly for $\barh_i$. The reason is that
one-loop rapidity evolution for $\hbsi(x_\bu,x_\perp)\Gamma\hsi(0)$ in the
Sudakov region is the same for all matrices $\Gamma$  between $\hbsi(x_\bu,x_\perp)$
and $\hpsi_(0)$ due to the fact that the “handbag”  diagram in Fig. <ref>c is small
and in two other diagrams (as well as self-energy corrections) the  matrix $\Gamma$  between $\hbsi(x_\bu,x_\perp)$
and $\hpsi(0)$ just multiplies the result of calculation.

Typical diagrams for the rapidity evolution of quark TMD in the Sudakov regime.

Next, to write down  the product of $m^2W_{\mu\nu}^H$ and the coefficient function
we need modified derivatives of $h_i$'s of Eq. (<ref>) type:

\begin{eqnarray}
\tartial_i h_j(\alpha_q,b_\perp;\vigma_p)~\equiv~\big(\partial_i-{\alpha_sc_F\over \pi}{b^i\over b_\perp^2}\ln\vigma_p\big)h_j(\alpha_q,b_\perp;\vigma_p),
\nonumber\\
\tartial_i h_j(\beta_q,b_\perp;\vigma_t)~\equiv~\big(\partial_i-{\alpha_sc_F\over \pi}{b^i\over b_\perp^2}\ln\vigma_t\big)h_j(\beta_q,b_\perp;\vigma_t)
\label{htartial}
\end{eqnarray}

(and similarly for $\barh$'s) so that  $\tartial_i h_j$ will satisfy same evolution equations (<ref>)  as $h_j$. We get then

\begin{eqnarray}
W_{\mu\nu}^{1H}(\alpha_q,\beta_q,b_\perp)~=~\sum_{\rm flavors} e_f^2W_{\mu\nu}^{fH}(\alpha_q,\beta_q,b_\perp),
\nonumber\\
{m^2\over 4\pi^2}W_{\mu\nu}^{fH}(\alpha_q,\beta_q,b_\perp)~=~M\Big(g_{\mu\nu}^\perp h^f_j\barh^{fj}-h^f_\mu\barh^f_\nu-h_\nu\barh^f_\mu
+~2{q^\parallel_\mu q^\parallel_\nu -\tilq_\mu\tilq_\nu\over Q_\parallel^4}\tartial^i h_i^f\tartial^j\barh_j^f
\nonumber\\
+~{q^\parallel_\mu\over Q_\parallel^2}\big[(i\tartial^ih^f_i)\barh^f_\nu+h^f_\nu i\tartial^i\barh^f_i\big]
+~{\tilq_\mu\over Q_\parallel^2}\big[(i\tartial^ih^f_i)\barh^f_\nu-h^f_\nu i\tartial^i\barh^f_i\big]~+~\mu\leftrightarrow\nu\Big)
-~{\tilq_\mu\tilq_\nu+q^\parallel_\mu q^\parallel_\nu \over Q_\parallel^4}
\nonumber\\
\times~\big[\partial_\perp^2(Mh^f_i\barh_f^i)+M\tartial^i h^f_j\tartial^i\barh^f_j\big]
-~{q^\parallel_\mu\tilq_\nu+ \tilq_\mu q^\parallel_\nu \over Q_\parallel^4}\partial_i[M(h^f_j\tartial\barh_f^j-\barh^f_j\tartial h_f^j)]~+~h\leftrightarrow\barh
\label{w1hcoordm}
\end{eqnarray}

where $M=M(\alpha_q,\beta_q;\vigma_p,\vigma_t)$ and
$h_i=h_i(\alpha_q,k_\perp,\vigma_p)$, $\barh_i\equiv\barh_i(\beta_q,(q-k)_\perp,\vigma_t)$ everywhere except
$h\leftrightarrow\barh$ terms where it is opposite, cf. Eq. (<ref>).

Let us comment on  the choice of “rapidity divides” $\vigma_p$ and $\vigma_t$ in the product

$M(\alpha_q,\beta_q,b_\perp;\vigma_p,\vigma_t)f_1(\alpha_q,b_\perp;\vigma_p)\barf_1(\beta_q,b_\perp;\vigma_t)$ (and in similar $Mh\barh$ product).
As we mentioned in the beginning of this Section, in order to avoid double counting one should write down  factorization of the amplitude in projectile, target and central fields at $\vigma_p,\vigma_t\sim 1$. After that, as discussed in Ref. [54],
one can use the double-log Sudakov evolution (<ref>) until

\begin{equation}
\vigma_p\geq \chvigma_p,~~~~ \chvigma_p\equiv{1\over\alpha_qb_\perp\sqrt{s}},~~~~~~~~~\vigma_t\geq \chvigma_t,~~~~\chvigma_t\equiv{1\over\beta_qb_\perp\sqrt{s}}
\label{vigmas}
\end{equation}

At this point, the result of Sudakov evolution is

\begin{eqnarray}
\big[f_1(\alpha_q,b_\perp;\chvigma_p)\barf_1(\beta_q,b_\perp;\chvigma_t)+f_1\leftrightarrow \barf_1\big]
\nonumber\\
=~e^{-{\alpha_s c_F\over 2\pi}\ln^2{Q^2b_\perp^2}}\big[f_1(\alpha_q,b_\perp;\chvigma_p)\barf_1(\beta_q,b_\perp;\chvigma_t)+f_1\leftrightarrow \barf_1\big]
\end{eqnarray}

so in the final result (<ref>) one should take

\begin{equation}
Mf\barf ~\rightarrow~e^{-{\alpha_s N_c\over 2\pi}\ln^2{Q^2b_\perp^2}}\Big[1+O\big(\alpha_s\big)\Big]
\label{fsud}
\end{equation}

Similarly, for $W_{\mu\nu}^{1H}(\alpha_q,\beta_q,b_\perp)$ one should take

\begin{equation}
Mh_i\barh_j ~\rightarrow~e^{-{\alpha_s N_c\over 2\pi}\ln^2{Q^2b_\perp^2}}\Big[1+O\big(\alpha_s\big)\Big]
\label{hsud}
\end{equation}

at the end of Sudakov evolution (<ref>). It should be emphasized that since factor $M$ is universal for (<ref>) and (<ref>),
our estimates of asymmetries in Sect. <ref> are not affected by summation of Sudakov double logs.

§.§.§ Rapidity-only cutoff for TMD

As discussed in Refs. [29, 28], from the rapidity factorization (<ref>) we get
TMDs with rapidity-only cutoff $|\alpha|<\sigma_t$ or $|\beta|<\sigma_p$ (or with modifications (<ref>)). Such cutoff, relevant for small-$x$ physics, is
different from the combination of UV and rapidity cutoffs for TMDs used by moderate-$x$ community, see the analysis in two  [56, 57, 58]
and three [59] loops.
For the tree-level formulas of Sect. <ref>, this difference in cutoffs does not matter, but if one uses the formulas from Sect. <ref>
and integrates models for TMDs with Sudakov factor $M$ of Eq. (<ref>), one has to relate TMDs with rapidity-only cutoffs
to the TMD models with conventional cutoffs.
This requires calculations at the NLO level which are in progress.

§ CONCLUSIONS AND OUTLOOK

Main result of this paper is Eq. (<ref>) which gives the DY hadronic tensor for electromagnetic current at small $x$
with gauge invariance at the ${1\over Q^2}$ level.
The part (<ref>), determined by leading-twist TMDs $f_1$ and $h_1^\perp$, is manifestly gauge invariant.
The only non-gauge invariant term at the ${1\over Q^2}$ level is Eq. (<ref>) with  transverse $\mu$ and $\nu$
which is $\sim{q^\perp_\mu q^\perp_\nu\over Q^2}$ times twist-3 TMDs. Also, in the leading-$N_c$  approximation
the only structure affected by those terms is $W_\dd$, all other structures are calculated up to $O\big({q_\perp^4\over Q^4}\big)$ terms.
It is interesting to note that ${1\over Q^2}$ terms necessary for gauge invariance are calculated more than than two decades after
the calculation of ${1\over Q}$ corrections in Ref. [20].

It should be emphasized that,  as discussed above, our rapidity factorization is different
from the standard factorization scheme for particle production in hadron-hadron
scattering, namely splitting  the diagrams in collinear to projectile part, collinear to
target part, hard factor, and soft factor [9]. Here we factorize only in rapidity
and the $Q^2$ evolution arises from $k_\perp^2$ dependence of the rapidity evolution kernels, same
as in the BK (and NLO BK [60]) equations. Also, since matrix elements of TMD operators
with our rapidity cutoffs are UV-finite [61, 62], the only UV divergencies
in our approach are usual
UV divergencies absorbed in the QCD running coupling. For the tree-level result (<ref>) this does not matter,
but if one intends to use the result like (<ref>) with Sudakov logarithms for conventional TMDs with double
UV and rapidity cutoffs, one needs to relate our TMDs with rapidity-only cutoff to conventional TMDs.
Needless to say,
the gauge-invariant tree-level result (<ref>) should be correct for TMDs with any cutoffs.

An obvious outlook is to extend these results to the “real” DY process involving $Z$-boson contributions which are  relevant for
our kinematics.
Another outlook is the one-loop calculations in this rapidity-based factorization and comparison to resummations
of large $\ln x$ and $\ln{Q^2/Q_\perp^2}$ based on usual
collinear factorization, see e.g. Refs. [63, 34].  The study is in progress.

The author is grateful to  V. Braun,  A. Prokudin,  A. Radyushkin, J. Qiu, and A. Vladimirov for valuable discussions. This  work is
supported by Jefferson Science Associates, LLC under the U.S. DOE contract #DE-AC05-06OR23177
and by U.S. DOE grant #DE-FG02-97ER41028.

§ APPENDIX

§.§ Formulas with Dirac matrices

§.§.§ Fierz transformation

First, let us write down Fierz transformation for symmetric hadronic tensor

\begin{eqnarray}
\half[(\bsi\gamma_\mu\chi)(\bhi\gamma_\nu\psi)+\mu\leftrightarrow\nu]
\label{fierz}\\
-{1\over 4}\big(\delta_\mu^\alpha\delta_\nu^\beta+\delta_\nu^\alpha\delta_\mu^\beta-g_{\mu\nu}g^{\alpha\beta}\big)
\big[(\bsi\gamma_\alpha\psi)(\bhi\gamma_\beta\chi)
\nonumber\\
+~{1\over 4}\big(\delta_\mu^\alpha\delta_\nu^\beta+\delta_\nu^\alpha\delta_\mu^\beta-\half g_{\mu\nu}g^{\alpha\beta}\big)
(\bsi\sigma_{\alpha\xi}\psi)(\bhi\sigma_\beta^{~\xi}\chi)-{g_{\mu\nu}\over 4}(\bsi\psi)(\bhi\chi)+{g_{\mu\nu}\over 4}(\bsi\gamma_5\psi)(\bhi\gamma_5\chi)
\nonumber
\end{eqnarray}

§.§.§ Formulas with $\sigma$-matrices

It is convenient to define
[We use conventions from Bjorken & Drell where $\epsilon^{0123}=-1$ and
\gamma^\mu\gamma^\nu\gamma^\lambda=g^{\mu\nu}\gamma^\lambda +g^{\nu\lambda}\gamma^\mu-g^{\mu\lambda}\gamma^\nu
Also, with this convention $\tigma_{\mu\nu}\equiv\half \epsilon_{\mu\nu\lambda\rho}\sigma^{\lambda\rho}=i\sigma_{\mu\nu}\gamma_5$.

\begin{equation}
\epsilon_{ij}~\equiv~ {2\over s}\epsilon_{\ast\bu ij}~=~ {2\over s}p_2^\mu p_1^\nu\epsilon_{\mu\nu ij}
\label{eps2}
\end{equation}

such that $\epsilon_{12}~=~1$ and $\epsilon_{ij}\epsilon_{kl}~=~g_{ik}g_{jl}-g_{il}g_{jk}$.
The frequently used formula is

\begin{equation}
\hspace{-1mm}
\sigma_{\mu\nu}\sigma_{\alpha\beta}~=~(g_{\mu\alpha}g_{\nu\beta}-g_{\mu\beta}g_{\nu\alpha})-i\epsilon_{\mu\nu\alpha\beta}\gamma_5
\label{sigmasigma}
\end{equation}

with variations

\begin{eqnarray}
{2\over s}\sigma_{\bu i}\sigma_{\ast j}~=~g_{ij}-i\epsilon_{ij}\gamma_5
-i\sigma_{ij}-{2i\over s}g_{ij}\sigma_{\bu\ast},~~~
{2\over s}\sigma_{\ast i}\sigma_{\bu j}~=~g_{ij}+i\epsilon_{ij}\gamma_5
-i\sigma_{ij}+{2\over s}g_{ij}\sigma_{\ast\bu},
\nonumber\\
\sigma_{i j}\sigma_{\bu k}=-\sigma_{\bu k}\sigma_{i j}~=-ig_{ik}\sigma_{\bu j}+ig_{jk}\sigma_{\bu i},~~~
\sigma_{i j}\sigma_{\ast k}=-\sigma_{\ast k}\sigma_{i j}~=-ig_{ik}\sigma_{\ast j}+ig_{jk}\sigma_{\ast i}
\nonumber\\
\label{sigmasigmas}
\end{eqnarray}

We need also the following formulas with $\sigma$-matrices in different matrix elements

\begin{eqnarray}
\tigma_{\mu\nu}\otimes\tigma_{\alpha\beta}~
\nonumber\\
\label{tigmi}
\end{eqnarray}

\begin{eqnarray}
\tigma_{\mu\xi}\otimes\tigma_\nu^{~\xi}~=~-{g_{\mu\nu}\over 2}\sigma_{\xi\eta}\otimes\sigma^{\xi\eta}+\sigma_{\nu\xi}\otimes\sigma_\mu^{~\xi}
\label{formulaxz}\\
\sigma_{\mu\xi}\gamma_5\otimes\sigma_\nu^{~\xi}\gamma_5+\mu\leftrightarrow\nu-{g_{\mu\nu}\over 2}\sigma_{\xi\eta}\gamma_5\otimes\sigma^{\xi\eta}\gamma_5
~=~-[\sigma_{\mu\xi}\otimes\sigma_\nu^{~\xi}+\mu\leftrightarrow\nu-{g_{\mu\nu}\over 2}\sigma_{\xi\eta}\otimes\sigma^{\xi\eta}]
\nonumber
\end{eqnarray}

\begin{eqnarray}
\sigma_\ast^{~k}\otimes\gamma_i\sigma_{\bu k}\gamma_j~=~\hatp_2\gamma^k\otimes\notp_1\gamma_i\gamma_k\gamma_j
\nonumber\\
=~\hatp_2(g_{ik}\gamma_j+g_{jk}\gamma_i-g_{ij}\gamma_k)\otimes \notp_1\gamma^k~=~(\gamma_j\sigma_\ast^{~k}\gamma_i)\otimes\sigma_{\bu k}
\label{sisigaga}
\end{eqnarray}

We will need also

\begin{eqnarray}
\notp_2\otimes\gamma_i\notp_1\gamma_j+\notp_2\gamma_5\otimes\gamma_i\notp_1\gamma_j\gamma_5
\label{flagamma}
\end{eqnarray}

§.§.§ Formulas with $\gamma$-matrices and one gluon field

In the gauge $A_\bu=0$ the field $A_i$ can be represented as

\begin{equation}
A_i(x_\bu,x_\perp)~=~{2\over s}\!\int_{-\infty}^{x_\bu}\!dx'_\bu~A_{\ast i}(x'_\bu,x_\perp)
\label{ai}
\end{equation}

(see eq. (<ref>)). We define “dual” fields by

\begin{equation}
\tilde{A}_i(x_\bu,x_\perp)~=~{2\over s}\!\int_{-\infty}^{x_\bu}\!dx'_\bu~\tilde{A}_{\ast i}(x'_\bu,x_\perp),
~~~\tilde{B}_i(x_\ast,x_\perp)~=~{2\over s}\!\int_{-\infty}^{x_\ast}\!dx'_\ast~\tilde{B}_{\bu i}(x'_\ast,x_\perp),
\label{abefildz}
\end{equation}

where $\tilF_{\mu\nu}=\half\epsilon_{\mu\nu\lambda\rho}F^{\lambda\rho}$ as usual.
With this definition we have $\tilA_i=-\epsilon_{ij}A^j$  and $\tilB_i=\epsilon_{ij}B^j$ so

\begin{equation}
\notp_2\brA_i~=~-\notA\notp_2\gamma_i,~~~~\brA_i\notp_2~=~-\gamma_i\notp_2\notA,~~~~
\notp_1\breB_i~=~-\notB\notp_1\gamma_i,~~~~\breB_i\notp_1~=~-\gamma_i\notp_1\notB
\label{glavla}
\end{equation}

\begin{equation}
\brA_i\equiv A_i-i\tilA_i\gamma_5,~~~\breB_i\equiv B_i-i\tilB_i\gamma_5
\end{equation}

We also use

\begin{eqnarray}
~=~-\notp_2\brA_n\otimes \notp_1-\notp_2\brA_n\gamma_5\otimes \notp_1\gamma_5
\nonumber\\
~=~-\brA_n\notp_2\otimes \notp_1-\brA_n\notp_2\gamma_5\otimes \notp_1\gamma_5
\nonumber\\
\gamma_n\slashed{p}_2\gamma^i\otimes\slashed{p}_1 B_i
~+~ \gamma_n\slashed{p}_2\gamma^i\gamma_5 \otimes\slashed{p}_1\gamma_5 B_i
~=~-\slashed{p}_2\otimes \slashed{p}_1\breB_n -
\slashed{p}_2 \gamma_5\otimes \slashed{p}_1\breB_n\gamma_5
\nonumber\\
\gamma^i\slashed{p}_2\gamma_n \otimes\slashed{p}_1B_i
~+~\gamma^i \slashed{p}_2\gamma_n\gamma_5 \otimes\slashed{p}_1\gamma_5 B_i
~=~-\slashed{p}_2\otimes \breB_n\slashed{p}_1 -
\slashed{p}_2 \gamma_5\otimes \breB_n\slashed{p}_1\gamma_5
\label{gammas1fild}
\end{eqnarray}

\begin{eqnarray}
{2\over s}\big[\notp_1\notp_2\gamma_i\otimes B^i\gamma_n
+\notp_1\notp_2\gamma_i\gamma_5\otimes B^i\gamma_n\gamma_5\big]
\nonumber\\
{2\over s}\big[\gamma_i\notp_2\notp_1\otimes B^i\gamma_n
+\gamma_i\notp_2\notp_1\gamma_5\otimes B^i\gamma_n\gamma_5\big]
\nonumber
\end{eqnarray}

§.§.§ Formulas with $\gamma$-matrices and two gluon fields

With definition (<ref>), we have the following formulas

\begin{eqnarray}
A_i\otimes \tilB_j=g_{ij}\tilA_k\otimes B^k-\tilA_j\otimes B_i,~~~\tilde{A}_i\otimes B_j=g_{ij}A_k\otimes \tilB^k-A_j\otimes \tilB_i
\label{abeznaki}\\
\tilde{A}_i\otimes \tilB_j=-g_{ij}A_k\otimes B^k+A_j\otimes B_i, ~~~\Rightarrow~~~\tilde{A}_i\otimes\tilde{B}^i=-A_i\otimes B^i,~~\tilde{A}_i\otimes B^i= A_i\otimes\tilde{B}^i
\nonumber
\end{eqnarray}

Using these formulas, after some algebra one obtains

\begin{eqnarray}
\gamma_m\notp_2\gamma_jA^i\otimes\gamma_n\notp_1\gamma_iB^j
\nonumber\\
\gamma_j\notp_2\gamma_mA^i\otimes\gamma_n\notp_1\gamma_iB^j
\nonumber\\
\gamma_m\notp_2\gamma_jA^i\otimes\gamma_i\notp_1\gamma_nB^j
\nonumber\\
\gamma_j\notp_2\gamma_mA^i\otimes\gamma_i\notp_1\gamma_nB^j
\nonumber\\
\label{gammas7}
\end{eqnarray}

\begin{eqnarray}
\notp_2\brA_m\otimes\notp_1\breB_n+\notp_2\brA_n\gamma_5\otimes\notp_1\breB_m\gamma_5~=~g_{mn}\notp_2\brA_k\otimes\notp_1\breB^k
\nonumber\\
\notp_2\brA_m\otimes\breB_n\notp_1+\notp_2\brA_n\gamma_5\otimes\gamma_5\breB_m\notp_1~=~g_{mn}\notp_2\brA_k\otimes\breB^k\notp_1
\nonumber\\
\brA_m\notp_2\otimes\notp_1\breB_n+\gamma_5\brA_n\notp_2\otimes\notp_1\breB_m\gamma_5~=~g_{mn}\brA_k\notp_2\otimes\notp_1\breB^k
\nonumber\\
\brA_m\notp_2\otimes\breB_n\notp_1+\gamma_5\brA_n\notp_2\otimes\gamma_5\breB_m\notp_1~=~g_{mn}\brA_k\notp_2\otimes\breB^k\notp_1
\label{gammas6}
\end{eqnarray}

The corollary of Eq. (<ref>) is

\begin{eqnarray}
\notp_2\brA_k\gamma_5\otimes\notp_1\breB^k\gamma_5~=~\notp_2\brA_k\otimes\notp_1\breB^k,~~~~~~~~
\notp_2\brA_k\gamma_5\otimes\gamma_5\breB^k\notp_1~=~\notp_2\brA_k\otimes\breB^k\notp_1
\nonumber\\
\gamma_5\brA_k\notp_2\otimes\notp_1\breB^k\gamma_5~=~\brA_k\notp_2\otimes\notp_1\breB^k,~~~~~~~~
\gamma_5\brA_k\notp_2\otimes\gamma_5\breB^k\notp_1~=~\brA_k\notp_2\otimes\breB^k\notp_1
\label{gammas6a}
\end{eqnarray}

From Eqs. (<ref>) and (<ref>) one easily obtains

\begin{equation}
\gamma_m\notp_2\gamma_jA^i\otimes\gamma_n\notp_1\gamma_iB^j
~+~m\leftrightarrow n~=~2g_{mn}\notp_2\brA_k\otimes\notp_1\breB^k
\label{formula9}
\end{equation}

\begin{eqnarray}
\gamma_m\notp_2\gamma_jA^i\otimes\gamma_n\notp_1\gamma_iB^j
~-~m\leftrightarrow n~
\nonumber\\
=~2\notp_2\brA_n\otimes\notp_1\breB_m~-~m\leftrightarrow n,
\nonumber\\
\gamma_j\notp_2\gamma_mA^i\otimes\gamma_i\notp_1\gamma_nB^j
~-~m\leftrightarrow n~
\nonumber\\
=~2\brA_n\notp_2\otimes\breB_m\notp_1~-~m\leftrightarrow n
\label{formula9a}
\end{eqnarray}

We need also formulas

\begin{eqnarray}
{4\over s^2}A^i\notp_1\notp_2\gamma_j\otimes B^j\notp_1\notp_2\gamma_i
\nonumber\\
A^i\gamma_j\otimes B^j\gamma_i-iA^i\gamma_j\gamma_5\otimes \tilB^j\gamma_i
+i\tilA^i \gamma_j\otimes  B^j\gamma_i\gamma_5+\tilA^i\gamma_j\gamma_5\otimes \tilB^j\gamma_i\gamma_5,
\nonumber\\
{4\over s^2}\big(A^i\notp_1\notp_2\gamma_j\otimes B^j\notp_1\notp_2\gamma_i
+A^i\notp_1\notp_2\gamma_j\gamma_5\otimes B^j\notp_1\notp_2\gamma_i\gamma_5\big)
\nonumber\\
\nonumber\\
\gamma_i\brA_j\gamma_5\otimes\gamma_j\brA_i\gamma_5~=~\gamma_i\brA_j\otimes\gamma^i\breB^j-\gamma_i\brA^i\otimes\gamma_j\breB^j
\label{gammas10}
\end{eqnarray}

\begin{eqnarray}
A_k\gamma_i\slashed{p}_2\gamma^j\otimes B_j\gamma^i\slashed{p}_1\gamma^k~=~
\slashed{p}_2\brA_i\otimes \slashed{p}_1\breB^i
\nonumber\\
A_k\gamma^j\slashed{p}_2\gamma_i\otimes B_j\gamma^k\slashed{p}_1\gamma^i~=~
\brA_i\slashed{p}_2\otimes \breB_i\slashed{p}_1
\nonumber\\
A_k\gamma_i\slashed{p}_2\gamma_j\otimes B^j\gamma^k\slashed{p}_1\gamma^i
\nonumber\\
A_k\gamma_j\slashed{p}_2\gamma_i\otimes B^j\gamma^i\slashed{p}_1\gamma^k
\label{gammas11}
\end{eqnarray}

\begin{eqnarray}
A^k\gamma_m\slashed{p}_2\gamma_j\otimes B^j\gamma_n\slashed{p}_1\gamma_k +m\leftrightarrow n
-g_{mn}A^k\gamma_i\slashed{p}_2\gamma_j\otimes B^j\gamma^i\slashed{p}_1\gamma_k
\nonumber\\
=~\brA_m\notp_2\otimes\breB_n\notp_1+m\leftrightarrow n-g_{mn}\brA_k\notp_2\otimes\breB^k\notp_1,
\nonumber\\
A^k\gamma_j\slashed{p}_2\gamma_m\otimes B^j\gamma_k\slashed{p}_1\gamma_n+m\leftrightarrow n
-g_{mn}A^k\gamma_j\slashed{p}_2\gamma_i\otimes B^j\gamma_k\slashed{p}_1\gamma^i
\nonumber\\
=~\notp_2\brA_m\otimes\notp_1\breB_n+m\leftrightarrow n-g_{mn}\notp_2\brA_k\otimes\notp_1\breB^k,
\nonumber\\
A^k\gamma_m\slashed{p}_2\gamma_j\otimes B^j\gamma_k\slashed{p}_1\gamma_n+m\leftrightarrow n
-g_{mn}A^k\gamma_j\slashed{p}_2\gamma_i\otimes B^j\gamma_k\slashed{p}_1\gamma^i
\nonumber\\
=~\brA_m\notp_2\otimes\notp_1\breB_n+m\leftrightarrow n-g_{mn}\brA_k\notp_2\otimes\notp_1\breB^k,
\nonumber\\
A^k\gamma_j\slashed{p}_2\gamma_m\otimes B^j\gamma_n\slashed{p}_1\gamma_k+m\leftrightarrow n
-g_{mn}A^k\gamma_j\slashed{p}_2\gamma_i\otimes B^j\gamma_k\slashed{p}_1\gamma^i
\nonumber\\
=~\notp_2\brA_m\otimes\breB_n\notp_1+m\leftrightarrow n-g_{mn}\notp_2\brA_k\otimes\breB^k\notp_1,
\label{formulas67}
\end{eqnarray}

\begin{eqnarray}
{2\over s}\big[A_i\notp_1\notp_2\gamma^j\otimes B_j\gamma_n\notp_1\gamma^i
+A_i\notp_1\notp_2\gamma^j\gamma_5\otimes B_j\gamma_{\nu_\perp}\notp_1\gamma^i\gamma_5\big]
\label{gammas16}\\
=~-\gamma_i\brA_n\otimes \notp_1\breB^i - \gamma_i\brA_n\gamma_5\otimes \notp_1\breB^i\gamma_5
\nonumber\\
{2\over s}\big[A_i\gamma_n\notp_2\gamma^j\otimes B_j\notp_2\notp_1\gamma^i
+A_i\gamma_n\notp_2\gamma^j\gamma_5\otimes B_j\notp_2\notp_1\gamma^i\gamma_5\big]
\nonumber\\
\nonumber
\end{eqnarray}

§.§ Parametrization of leading-twist matrix elements

Let us  first consider matrix elements of operators without $\gamma_5$. The standard parametrization of quark TMDs reads
(see e.g. Ref. [64]))

\begin{eqnarray}
{1\over 16\pi^3}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(x_\bu,x_\perp)\gamma^\mu\psi_f(0)|A\rangle
\label{Amael}\\
=~p_1^\mu f_1^f(\alpha,k_\perp)
+k_\perp^\mu f_\perp^f(\alpha,k_\perp)+p_2^\mu{2m^2_N\over s}f_3^f(\alpha,k_\perp),
\nonumber\\
{1\over 16\pi^3}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(x_\bu,x_\perp)\psi_f(0)|A\rangle
\nonumber
\end{eqnarray}

for quark distributions in the projectile and

\begin{eqnarray}
{1\over 16\pi^3}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(0)\gamma^\mu\psi_f(x_\bu,x_\perp)|A\rangle
\label{baramael}\\
=~-p_1^\mu \barf_1^f(\alpha,k_\perp)
-k_\perp^\mu\barf _\perp^f(\alpha,k_\perp)-p_2^\mu{2m^2_N\over s}\barf_3^f(\alpha,k_\perp),
\nonumber\\
{1\over 16\pi^3}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(0)\psi_f(x_\bu,x_\perp)|A\rangle
\nonumber
\end{eqnarray}

for the antiquark distributions.
[In an arbitrary gauge, there are gauge links to $-\infty$ as displayed in  eq. (<ref>).]

The corresponding matrix elements for the target are obtained by trivial replacements $p_1\leftrightarrow p_2$, $x_\bu\leftrightarrow x_\ast$
and $\alpha\leftrightarrow\beta$:

\begin{eqnarray}
{1\over 16\pi^3}\!\int\!dx_\ast d^2x_\perp~e^{-i\beta x_\ast+i(k,x)_\perp}
~\langle B|\bsi_f(x_\ast,x_\perp)\gamma^\mu\psi_f(0)|B\rangle
\label{Bmael}\\
=~p_2^\mu f_1^f(\beta,k_\perp)+k_\perp^\mu f_\perp^f(\beta,k_\perp)
+p_1^\mu{2m^2_N\over s}f_3^f(\beta,k_\perp),
\nonumber\\
{1\over 16\pi^3}\!\int\!dx_\ast d^2x_\perp~e^{-i\beta x_\ast+i(k,x)_\perp}
~\langle B|\bsi_f(x_\ast,x_\perp)\psi_f(0)|B\rangle
\nonumber
\end{eqnarray}

\begin{eqnarray}
{1\over 16\pi^3}\!\int\!dx_\ast d^2x_\perp~e^{-i\beta x_\ast+i(k,x)_\perp}
~\langle B|\bsi_f(0)\gamma^\mu\psi_f(x_\ast,x_\perp)|B\rangle
\label{barbmael}\\
=~-p_2^\mu \barf_1^f(\beta,k_\perp)
-k_\perp^\mu\barf _\perp^f(\beta,k_\perp)-p_1^\mu{2m^2_N\over s}\barf_3^f(\beta,k_\perp),
\nonumber\\
{1\over 16\pi^3}\!\int\!dx_\ast d^2x_\perp~e^{-i\beta x_\ast+i(k,x)_\perp}
~\langle B|\bsi_f(0)\psi_f(x_\ast,x_\perp)|B\rangle
\nonumber
\end{eqnarray}

Matrix elements of operators with $\gamma_5$ are parametrized as follows:

\begin{eqnarray}
{1\over 16\pi^3}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(x_\bu,x_\perp)\gamma^\mu\gamma_5\psi_f(0)|A\rangle
~=~-i\epsilon_{\mu_\perp i}k^ig^\perp_f(\alpha,k_\perp),
\nonumber\\
{1\over 16\pi^3}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(0)\gamma^\mu\gamma_5\psi_f(x_\bu,x_\perp)|A\rangle
~=~-i\epsilon_{\mu_\perp i}k^i\barg^\perp_f(\alpha,k_\perp)
\nonumber\\
\label{mael5}
\end{eqnarray}

The corresponding matrix elements for the target are obtained by trivial replacements $p_1\leftrightarrow p_2$, $x_\bu\leftrightarrow x_\ast$
and $\alpha\leftrightarrow\beta$ similarly to eq. (<ref>).

The parametrization of time-odd Boer-Mulders TMDs are

\begin{eqnarray}
{1\over 16\pi^3}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(x_\bu,x_\perp)\sigma^{\mu \nu}\psi_f(0)|A\rangle
\nonumber\\
=~{1\over m_N}(k_\perp^\mu p_1^\nu -\mu\leftrightarrow\nu)h_{1f}^\perp(\alpha,k_\perp)
+{2m_N\over s}(p_1^\mu p_2^\nu-\mu\leftrightarrow\nu)h_{f}(\alpha,k_\perp)
\nonumber\\
+~{2m_N\over s}(k_\perp^\mu p_2^\nu -\mu\leftrightarrow\nu)h_{3f}^\perp(\alpha,k_\perp),
\nonumber\\
{1\over 16\pi^3}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(0)\sigma^{\mu \nu}\psi_f(x_\bu,x_\perp)|A\rangle
\nonumber\\
=~-{1\over m_N}(k_\perp^\mu p_1^\nu -\mu\leftrightarrow\nu)\barh_{1f}^\perp(\alpha,k_\perp)
-{2m_N\over s}(p_1^\mu p_2^\nu-\mu\leftrightarrow\nu)\barh_{f}(\alpha,k_\perp)
\nonumber\\
-~{2m_N\over s}(k_\perp^\mu p_2^\nu -\mu\leftrightarrow\nu)\barh_{3f}^\perp(\alpha,k_\perp)
\label{hmael}
\end{eqnarray}

and similarly for the target with usual replacements   $p_1\leftrightarrow p_2$, $x_\bu\leftrightarrow x_\ast$
and $\alpha\leftrightarrow\beta$.

Note that
the coefficients in front of $f_3$,  $g^\perp_f$, $h$ and $h_3^\perp$ in eqs. (<ref>),  (<ref>), (<ref>),  and  (<ref>)
contain an extra ${1\over s}$ since $p_2^\mu$ enters only through the direction
of gauge link so the result should not depend on rescaling $p_2\rightarrow\lambda p_2$. For this reason,  these functions do not contribute to $W(q)$ in our approximation.

Last but not least, an important point in our analysis is that any $f(x,k_\perp)$ may have only logarithmic dependence on Bjorken $x$ but
not the power dependence $\sim{1\over x}$.
Indeed, the low-$x$ behavior of TMDs is determined by pomeron exchange with the nucleon.  The interaction of TMD with BFKL pomeron
is specified by so-called impact factor and it is easy to check that the impact factors for all leading-twist TMDs are similar and do not give
extra ${1\over x}$ factors. The only ${1\over x}$ may had come from some unfortunate definition of TMD which includes factor $s$ artificially, but
from power counting (<ref>) we see that all definitions of leading-twist TMDs do not have such factors.

§.§ Matrix elements of quark-quark-gluon operators

In this section we will demonstrate that matrix elements of quark-antiquark-gluon operators
from section <ref> can be expressed in terms of leading-power
matrix elements from section <ref>.

First, let us note that operators ${1\over\alpha}$ and ${1\over\beta}$ in Eqs. (<ref>) are
replaced by $\pm{1\over\alpha_q}$ and $\pm{1\over\beta_q}$ in forward matrix elements. Indeed,

\begin{eqnarray}
\!\int\!dx_\bu ~e^{-i\alpha_q x_\bu}\langle\bar\Phi(x_\bu,x_\perp)\Gamma{1\over \alpha+\ie}\psi(0)\rangle_A
\label{maelqg1}\\
=~{1\over i}\!\int\!dx_\bu \!\int_{-\infty}^0 \!\!\!dx'_\bu~e^{-i\alpha_q x_\bu}\langle\bar\Phi(x_\bu,x_\perp)\Gamma\psi(x'_\bu,0_\perp)\rangle_A
=~{1\over\alpha_q}\!\int\!dx_\bu ~e^{-i\alpha x_\bu}\langle\bar\Phi(x_\bu,x_\perp)\Gamma\psi(0)\rangle_A
\nonumber
\end{eqnarray}

where $\bar\Phi(x_\bu,x_\perp)$ can be $\bsi(x_\bu,x_\perp)$ or  $\bsi(x_\bu,x_\perp)A_i(x_\bu,x_\perp)$ and $\Gamma$ can be any $\gamma$-matrix.

\begin{eqnarray}
\!\int\!dx_\bu ~e^{-i\alpha_q x_\bu}\langle\big(\bar\psi{1\over \alpha-\ie}\big)(x_\bu,x_\perp)\Gamma\Phi(0)\rangle_A
=~{1\over\alpha_q}\!\int\!dx_\bu ~e^{-i\alpha x_\bu}\langle\bar\psi(x_\bu,x_\perp)\Gamma\Phi(0)\rangle_A
\label{maelqg2}\\
\!\int\!dx_\bu ~e^{-i\alpha_q x_\bu}\langle\big(\bar\psi{1\over \alpha-\ie}\big)(x_\bu,x_\perp)\Gamma{1\over \alpha+\ie}\psi(0)\rangle_A
=~{1\over\alpha_q^2}\!\int\!dx_\bu ~e^{-i\alpha_q x_\bu}\langle\bar\psi(x_\bu,x_\perp)\Gamma\psi(0)\rangle_A
\nonumber
\end{eqnarray}

where $\Phi(x_\bu,x_\perp)$ can be $\psi(x_\bu,x_\perp)$ or  $A_i(x_\bu,x_\perp)\psi(x_\bu,x_\perp)$. We need also

\begin{eqnarray}
\!\int\!dx_\bu ~e^{-i\alpha_q x_\bu}\langle\big(\bar\psi{1\over \alpha-\ie}\big)(0)\Gamma\Phi(x_\bu,x_\perp)\rangle_A
=~-{1\over\alpha_q}\!\int\!dx_\bu ~e^{-i\alpha x_\bu}\langle\bar\psi(0)\Gamma\Phi(x_\bu,x_\perp)\rangle_A
\nonumber\\
\!\int\!dx_\bu ~e^{-i\alpha_q x_\bu}\langle\bar\Phi(0)\Gamma{1\over \alpha+\ie}\psi(x_\bu,x_\perp)\rangle_A
=~-{1\over\alpha_q}\!\int\!dx_\bu ~e^{-i\alpha x_\bu}\langle\bar\Phi(0)\Gamma\psi(x_\bu,x_\perp)\rangle_A
\label{maelqg3}
\end{eqnarray}

The corresponding formulas for target matrix elements are obtained by substitution $\alpha\leftrightarrow\beta$ (and $x_\bu\leftrightarrow x_\ast$).

Next, we will use QCD equation of motion to reduce quark-quark-gluon TMDs to leading-twist TMDs (see Ref. [20]).
Let us start with matrix element

\begin{eqnarray}
\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}\langle A|\bar\psi(x_\bu,x_\perp)\slashed{p}_2\brA_i(x_\bu,x_\perp)\psi(0)|A\rangle
\label{tw3mael1}\\
=~-\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}\langle A|\bar\psi(x_\bu,x_\perp)\notA(x_\bu,x_\perp)\slashed{p}_2\gamma_i\psi(0)|A\rangle
\nonumber\\
=~\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\nonumber\\
\times~\big[\langle A|\bsi(x_\bu,x_\perp)\slashed{k}_\perp\slashed{p}_2\gamma_i\psi(0)|A\rangle
\langle A|\bsi(x_\bu,x_\perp)\!\stackrel{\leftarrow}{\notD}_\perp\!\gamma_j\slashed{p}_2\gamma_i\psi(0)|A\rangle\big].
\nonumber
\end{eqnarray}

Using QCD equations of motion (<ref>) we can rewrite the r.h.s. of eq. (<ref>) as

\begin{eqnarray}
\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}\big[\langle A|\bsi(x_\bu,x_\perp)\slashed{k}_\perp\slashed{p}_2\gamma_i\psi(0)|A\rangle
+\alpha_q\langle A|\bsi(x_\bu,x_\perp)\slashed{p}_1\slashed{p}_2\gamma_i\psi(0)|A\rangle\big]
\nonumber\\
=~\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\Big[
-k_i\langle A|\bsi(x_\bu,x_\perp)\slashed{p}_2\psi(0)|A\rangle
+~\alpha_q{s\over 2} \langle A|\bsi(x_\bu,x_\perp)\gamma_i\psi(0)|A\rangle
\nonumber\\
-~i\epsilon_{ij}k^j\langle A|\bsi(x_\bu,x_\perp)\slashed{p}_2\gamma_5\psi(0)|A\rangle
+i{s\over 2}\alpha\epsilon_{ij}
\langle A|\bsi(x_\bu,x_\perp)\gamma^j\gamma_5\psi(0)|A\rangle\Big]
\nonumber\\
=~-k_i8\pi^3sf_1(\alpha_q,k_\perp)+8\pi^3s\alpha_q k_i\big[ f_\perp(\alpha_q,k_\perp)+g^\perp(\alpha_q,k_\perp)\big],
\label{tw3mael2}
\end{eqnarray}

where we used parametrizations (<ref>) and (<ref>) for the leading power matrix elements.

Now, the second term in eq. (<ref>) contains extra $\alpha_q$ with respect to the first term
[As discussed in the end of Sect. <ref>, all leading-twist TMDs can have only logarithmic dependence on
Bjorken $x$ (which is here either $\alpha_q$ for the projectile or $\beta_q$ for the target matrix elements).
, so
it should be neglected in our kinematical region $s\gg Q^2\gg q_\perp^2$  and we get

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}\langle A|\bar\psi^f(x_\bu,x_\perp)\slashed{p}_2\brA_i(x_\bu,x_\perp)\psi^f(0)|A\rangle
\label{11.42}\\
=~-{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}\langle A|\bar\psi^f\notA(x_\bu,x_\perp)\slashed{p}_2\gamma_i\psi^f(0)|A\rangle
\nonumber
\end{eqnarray}

By complex conjugation

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\perp dx_\bu~e^{-i\alpha_q x_\bu+i(k,x)_\perp}
\langle A|\bsi _f(x_\bu,x_\perp)\brA_i(0)\slashed{p}_2\psi_f(0)|A\rangle
\label{11.36}\\
=~-{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}\langle A|\bar\psi^f(x_\bu,x_\perp)\gamma_i\slashed{p}_2\notA\psi^f(0)|A\rangle
\nonumber
\end{eqnarray}

For the  corresponding antiquark distributions we get

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle A|\bsi_f(0)\brA_i(x_\bu,x_\perp)\slashed{p}_2\psi_f(x_\bu,x_\perp)|A\rangle
\nonumber\\
=~{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\Big[-\langle A|\bsi(0)\gamma_i\slashed{p}_2\slashed{k}_\perp\psi(x_\bu,x_\perp)|A\rangle
\nonumber\\
-~i\langle A|\bsi(0)\gamma_i\slashed{p}_2\notD_\perp\psi(x_\bu,x_\perp)|A\rangle\Big]
\label{9.25}
\end{eqnarray}

\begin{equation}
\hspace{0mm}
{1\over 8\pi^3s}\!\int\! dx_\perp dx_\bu~e^{-i\alpha_q x_\bu+i(k,x)_\perp}
\langle A|\bsi _f(0)\slashed{p}_2\brA_i(0)\psi_f(x_\bu,x_\perp)|A\rangle
-k_i\barf_{1f}(\alpha_q,k_\perp) .
\label{11.45}
\end{equation}

The corresponding target matrix elements are obtained by trivial replacements
$x_\ast\leftrightarrow x_\bu$, $\alpha_q\leftrightarrow\beta_q$ and

Next, let us consider

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi(x_\bu,x_\perp)\slashed{p}_2\notA(x_\bu,x_\perp)\psi(0)|A\rangle
\label{eqm1}\\
=~{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\nonumber\\
\times~\Big[\langle A|\bsi(x_\bu,x_\perp)\slashed{k}_\perp\slashed{p}_2\psi(0)|A\rangle
+i\langle A|\bsi(x_\bu,x_\perp)\stackrel{\leftarrow}{\notD}_\perp\slashed{p}_2\psi(0)|A\rangle\Big].
\nonumber
\end{eqnarray}

Using QCD equation of motion and parametrization (<ref>), one can rewrite the r.h.s. of this equation as

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\Big[\langle A|\bsi(x_\bu,x_\perp)\slashed{k}_\perp\slashed{p}_2\psi(0)|A\rangle
+\alpha_q\langle A|\bsi(x_\bu,x_\perp)\slashed{p}_1\slashed{p}_2\psi(0)|A\rangle\Big]
\nonumber\\
=~i{k_\perp^2\over m_N}h_{1}^\perp(\alpha_q,k_\perp)+\alpha_q m_N\big[e(\alpha,k_\perp)+i h(\alpha,k_\perp)\big].
\label{eqm2}
\end{eqnarray}

Again, only the first term contributes in our kinematical region so we finally get

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi^f(x_\bu,x_\perp)\slashed{p}_2\notA(x_\bu,x_\perp)\psi^f(0)|A\rangle
~=~i{k_\perp^2\over m}h_{1f}^\perp(\alpha_q,k_\perp).
\nonumber\\
\label{11.41}
\end{eqnarray}

By complex conjugation we obtain

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi^f(x_\bu,x_\perp)\slashed{p}_2\notA(0)\psi^f(0)|A\rangle
~=~i{k_\perp^2\over m}h_{1f}^\perp(\alpha_q,k_\perp).
\nonumber\\
\label{11.42}
\end{eqnarray}

For corresponding antiquark distributions one gets in a similar way

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi^f(0)\slashed{p}_2\notA(x_\bu,x_\perp)\psi^f(x_\bu,x_\perp)|A\rangle
i{k_\perp^2\over m}\barh_{1f}^\perp(\alpha_q,k_\perp),
\nonumber\\
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi^f(0)\slashed{p}_2\notA(0)\psi^f(x_\bu,x_\perp)|A\rangle
i{k_\perp^2\over m}\barh_{1f}^\perp(\alpha_q,k_\perp).
\nonumber\\
\label{11.43}
\end{eqnarray}

The target matrix elements are obtained by usual replacements
$x_\ast\leftrightarrow x_\bu$, $\alpha_q\leftrightarrow\beta_q$ and

Finally, we need

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi(x_\bu,x_\perp)\notA(x_\bu,x_\perp)\slashed{p}_2\notA(0)\psi(0)|A\rangle
\nonumber\\
=~{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi(x_\bu,x_\perp)\Big(\slashed{k}_\perp
\big(\slashed{k}_\perp
\nonumber\\
=~{k_\perp^2\over 16\pi^3}f_1(\alpha_q,k_\perp)~+~O(\alpha_q,\beta_q)
\label{AA1}
\end{eqnarray}

and similarly

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi(x_\bu,x_\perp)\notA(x_\bu,x_\perp)\sigma_{\ast i}\notA(0)\psi(0)|A\rangle
\nonumber\\
=~{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi(x_\bu,x_\perp)\Big(\slashed{k}_\perp
+i\stackrel{\leftarrow}{\notD}\big)\sigma_{\ast i}
\big(\slashed{k}_\perp
\nonumber\\
=~{1\over 16\pi^3}{k_ik_\perp^2\over m}h_1^\perp(\alpha_q,k_\perp)~+~O(\alpha_q,\beta_q)
\label{AA2}
\end{eqnarray}

For corresponding antiquark distributions we get

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi(0)\notA(0)\slashed{p}_2\notA(x_\bu,x_\perp)\psi(x_\bu,x_\perp)|A\rangle
\nonumber\\
=~-{k_\perp^2\over 16\pi^3}\barf_1(\alpha_q,k_\perp)~+~O(\alpha_q,\beta_q)
\nonumber\\
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle A|\bsi(x_\bu,x_\perp)\notA(x_\bu,x_\perp)\sigma_{\ast i}\notA(0)\psi(0)|A\rangle
\nonumber\\
=~-{1\over 16\pi^3}{k_ik_\perp^2\over m}h_1^\perp(\alpha_q,k_\perp)~+~O(\alpha_q,\beta_q)
\end{eqnarray}

Also, as we saw in Sect. <ref>, at the leading order in $N_c$ there is one quark-antiquark-gluon
operator that does not reduce to twist-2 distributions. It can be parametrized as follows (cf. Eq. (<ref>))

\begin{eqnarray}
{1\over 16\pi^3}{2\over s}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(x_\bu,x_\perp)\big[A_i(x)\sigma_{\ast j}-\half g_{ij}A^k\sigma_{\ast k}(x)\big]\psi_f(0)|A\rangle~
\nonumber\\
=~-(k_ik_j+\half g_{ij}k_\perp^2){1\over m}h_A^f(\alpha,k_\perp),
\nonumber\\
{1\over 16\pi^3}{2\over s}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle A|\bsi_f(0)[A_i(0)\sigma_{\bu j}-\half g_{ij}A^k\sigma_{\bu k}(0)]\psi_f(x_\bu,x_\perp)|A\rangle~
\nonumber\\
=~-(k_ik_j+\half g_{ij}k_\perp^2){1\over m}\barh_A^f(\alpha,k_\perp)
\label{maelsa}
\end{eqnarray}

and similarly for the target matrix elements.

§.§ Parametrization of TMDs from section <ref>

We  parametrize TMDs from section <ref> as follows

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^{x_\bu}\! dx'_\bu\langle A|\bsi (x_\bu,x_\perp){2\slashed{p}_2\over s}
\big[F_{\ast i}(0)+ i\gamma_5\tilF_{\ast i}(0)\big]\psi(x'_\bu,x_\perp)|A\rangle
\nonumber\\
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^{x_\bu}\! dx'_\bu\langle A|\bsi (x_\bu,x_\perp){2\slashed{p}_2\over s}
\big[F_{\ast i}(0)- i\gamma_5\tilF_{\ast i}(0)\big]\psi(x'_\bu,x_\perp)|A\rangle
\nonumber\\
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^{x_\bu}\! dx'_\bu\langle A|\bsi (x'_\bu,x_\perp){2\slashed{p}_2\over s}
\big[F_{\ast i}(0)- i\gamma_5\tilF_{\ast i}(0)\big]\psi(x_\bu,x_\perp)|A\rangle
\nonumber\\
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^{x_\bu}\! dx'_\bu\langle A|\bsi (x'_\bu,x_\perp){2\slashed{p}_2\over s}
\big[F_{\ast i}(0)+i\gamma_5\tilF_{\ast i}(0)\big]\psi(x_\bu,x_\perp)|A\rangle
\nonumber\\
\label{paramj}
\end{eqnarray}

By complex conjugation we get

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^0\! dx'_\bu\langle A|\bar\psi(x'_\bu,0_\perp){2\slashed{p}_2\over s}
[F_{\ast i}(x)-i\gamma_5\tilF_{\ast i}(x)]\psi(0)|A\rangle,
\nonumber\\
=~k_i j_1^\star(\alpha,k_\perp),
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^0\! dx'_\bu\langle A|\bar\psi(x'_\bu,0_\perp){2\slashed{p}_2\over s}
\big[F_{\ast i}(x)+i\gamma_5\tilF_{\ast i}(x)\big]\psi(0)|A\rangle
\nonumber\\
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^0\! dx'_\bu\langle A|\bar\psi(0){2\slashed{p}_2\over s}
[F_{\ast i}(x)+i\gamma_5\tilF_{\ast i}(x)]\psi(x'_\bu,0_\perp)|A\rangle
\nonumber\\
=~k_i \barj_1^\star(\alpha,k_\perp),
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^0\! dx'_\bu\langle A|\bar\psi(0){2\slashed{p}_2\over s}
\big[F_{\ast i}(x)-i\gamma_5\tilF_{\ast i}(x)\big]\psi(x'_\bu,0_\perp)|A\rangle
\nonumber\\
\label{paramjstar}
\end{eqnarray}

Note that unlike two-quark matrix elements,
quark-quark-gluon ones may have  imaginary parts.

Target matrix elements are obtained by usual substitutions
$\alpha\leftrightarrow\beta$, $\slashed{p}_2\leftrightarrow\slashed{p}_1$, $x_\bu\leftrightarrow x_\ast$, and $\hatF_{\ast i}\leftrightarrow \hatF_{\bu i}$.

For completeness let us present the explicit form of the gauge links in an arbitrary gauge:

\begin{eqnarray}
\bsi (x'_\bu,x_\perp)
F_{\ast i}(0)\psi(x_\bu,x_\perp)
~\rightarrow~\bsi (x'_\bu,x_\perp)[x'_\bu,-\infty_\bu]_x[x_\perp,0_\perp]_{-\infty_\bu}
\\
\times~[-\infty_\bu,0]_{0_\perp}F_{\ast i}(0)
\nonumber
\end{eqnarray}

§.§ Gluon power corrections from $J_{A}^\mu(x)J_{A\mu}(0)$ terms

There is one more type of contributions proportional to the product of quark and gluon TMDs

\begin{eqnarray}
\nonumber\\
=~\sum_{\rm flavors}\Big(\big[\Bxi_{1}(x)\gamma^\mu\psi_A(x)\big]\big[\bar\psi_A(0)\gamma^\nu\Xi_{1}(0)\big]
\nonumber\\
\label{5.15}
\end{eqnarray}

where we neglected terms which cannot contribute to $W$ due to the reason discussed after eq. (<ref>), i.e.
that one hadron (“A” or “B”) cannot produce the DY pair on its own.

Let us  consider the first term in the r.h.s of this equation

\begin{eqnarray}
\cheW_{\mu\nu}(x)~=~\langle A,B|\big[\Bxi_{1}(x)\gamma_\mu\psi_A(x)\big]\big[\bar\psi_A(0)\gamma_\nu\Xi_{1}(0)\big]|A,B\rangle
\label{5.16}\\
=~-{g^2\over s(N_c^2-1)}
\langle \big(\bsi{1\over\alpha}\big)(x)\gamma_i\slashed{p}_2\gamma_\mu\psi(x)
\bsi(0)\gamma_\nu\slashed{p}_2\gamma_j{1\over\alpha}\psi(0)\rangle_A
\langle A^{ai}(x)A^{aj}(0)\rangle_B
\nonumber
\end{eqnarray}

To estimate the magnitude of this contribution, first note that

\begin{eqnarray}
\!\int\! dx_\ast ~e^{-i\beta_qx_\ast}~\langle B|A^a_i(x)A^{a}_j(0)|B\rangle
\label{gluonvklad}\\
{4\over s^2}\int\! dx_\ast ~e^{-i\beta_qx_\ast}\!\int_{-\infty}^{x_\ast}\! dx'_\ast \!\int_{-\infty}^{0}\! dx''_\ast~
\langle B|F^a_{\bu i}(x'_\ast,x_\perp)F^a_{\bu j}(x''_\ast,0_\perp)|B\rangle
\nonumber\\
=~{4\over\beta_q^2s^2}\!\int\! dx_\ast ~e^{-i\beta_qx_\ast}\langle B|F^a_{\bu i}(x_\ast,x_\perp)F^a_{\bu j}(0)|B\rangle
\nonumber\\
=~-{1\over\beta_q}8\pi^2\alpha_s\Big[\cald_g(\beta_q,x_\perp)+{1\over m^2}(2\partial_i\partial_j+g_{ij}\partial_\perp^2\calh_g(\beta_q,x_\perp)\big]
\nonumber
\end{eqnarray}

where we used parametrization (3.26) from Ref. [29]. Since the gluon TMDs $\cald_g(x_B,x_\perp)$ and $\calh_g(\beta_q,x_\perp)$ behave
only logarithmically as $x_B\rightarrow 0$ [62], the contribution of eq. (<ref>) to $W(q)$ is of order of
${m_\perp^2\over\beta_q s}\ll{m_\perp^2\over Q^2}$. (As discussed in Ref. [29],  the projectile TMD in the r.h.s. of eq. (<ref>) does not give ${1\over\alpha_q}$ after Fourier transformation).  Also, this contribution is $\sim{1\over N_c}$ with respect to our leading terms.

Similarly, all other terms in eq. (<ref>) are either
${m_\perp^2\over\beta_q s}$ or ${m_\perp^2\over\alpha_q s}$ times ${1\over N_c}$ so they can be neglected.
It is worth mentioning that if the DY pair is produced in the region of rapidity close to the projectile,
the contribution (<ref>) may be the most important since gluon parton densities  at small $x_B$
are larger than the quark ones.]

[1]
S.D. Drell and Tung-Mow Yan.
Massive Lepton Pair Production in Hadron-Hadron Collisions at
Phys. Rev. Lett., 25:316–320, 1970.
[Erratum: Phys.Rev.Lett. 25, 902 (1970)].

[2]
Jen-Chieh Peng and Jian-Wei Qiu.
Novel phenomenology of parton distributions from the
Drell–Yan process.
Prog. Part. Nucl. Phys., 76:43–75, 2014.

[3]
A.V. Radyushkin.
Massive Lepton Pair Production and Asymptotically Free Gauge
Phys. Lett. B, 69:245, 1977.
[Erratum: Phys.Lett.B 77, 461 (1978)].

[4]
H.David Politzer.
Gluon Corrections to Drell-Yan Processes.
Nucl. Phys. B, 129:301–318, 1977.

[5]
K. Kajantie and Risto Raitio.
Gluon Effects in Muon Pair Production.
Nucl. Phys. B, 139:72–84, 1978.

[6]
Guido Altarelli, G. Parisi, and R. Petronzio.
Transverse Momentum in Drell-Yan Processes.
Phys. Lett. B, 76:351–355, 1978.

[7]
R.Keith Ellis, Howard Georgi, Marie Machacek, H.David Politzer, and Graham G.
Perturbation Theory and the Parton Model in QCD.
Nucl. Phys. B, 152:285–329, 1979.

[8]
A.V. Efremov and A.V. Radyushkin.
Field Theoretic Treatment of High Momentum Transfer Processes. 3.
Gauge Theories.
Theor. Math. Phys., 44:774, 1981.

[9]
John Collins.
Foundations of perturbative QCD.
Cambridge University Press, 2013.

[10]
John C. Collins and Davison E. Soper.
Parton Distribution and Decay Functions.
Nucl. Phys., B194:445–492, 1982.

[11]
John C. Collins, Davison E. Soper, and George F. Sterman.
Transverse Momentum Distribution in Drell-Yan Pair and W and Z Boson
Nucl. Phys., B250:199–224, 1985.

[12]
Xiang-dong Ji, Jian-ping Ma, and Feng Yuan.
QCD factorization for semi-inclusive deep-inelastic scattering at
low transverse momentum.
Phys. Rev., D71:034005, 2005.

[13]
Miguel G. Echevarria, Ahmad Idilbi, and Ignazio Scimemi.
Factorization Theorem For Drell-Yan At Low $q_T$ And Transverse
Momentum Distributions On-The-Light-Cone.
JHEP, 07:002, 2012.

[14]
Christian W. Bauer, Dan Pirjol, and Iain W. Stewart.
Soft collinear factorization in effective field theory.
Phys. Rev., D65:054022, 2002.

[15]
Christian W. Bauer, Sean Fleming, Dan Pirjol, Ira Z. Rothstein, and Iain W.
Hard scattering factorization from effective field theory.
Phys. Rev., D66:014017, 2002.

[16]
Thomas Becher, Matthias Neubert, and Gang Xu.
Dynamical Threshold Enhancement and Resummation in Drell-Yan
JHEP, 07:030, 2008.

[17]
Thomas Becher and Matthias Neubert.
Drell-Yan Production at Small $q_T$, Transverse Parton
Distributions and the Collinear Anomaly.
Eur. Phys. J., C71:1665, 2011.

[18]
John C. Collins and Davison E. Soper.
Angular Distribution of Dileptons in High-Energy Hadron Collisions.
Phys. Rev. D, 16:2219, 1977.

[19]
John Collins and Ted Rogers.
Understanding the large-distance behavior of
transverse-momentum-dependent parton densities and the Collins-Soper
evolution kernel.
Phys. Rev., D91(7):074020, 2015.

[20]
P. J. Mulders and R. D. Tangerman.
The Complete tree level result up to order 1/Q for polarized deep
inelastic leptoproduction.
Nucl. Phys., B461:197–237, 1996.
[Erratum: Nucl. Phys.B484,538(1997)].

[21]
Xiang-Dong Ji.
Off forward parton distributions.
J. Phys. G, 24:1181–1205, 1998.

[22]
Pierre A.M. Guichon and M. Vanderhaeghen.
Virtual Compton scattering off the nucleon.
Prog. Part. Nucl. Phys., 41:125–190, 1998.

[23]
I.V. Anikin, B. Pire, and O.V. Teryaev.
On the gauge invariance of the DVCS amplitude.
Phys. Rev. D, 62:071501, 2000.

[24]
M. Penttinen, Maxim V. Polyakov, A.G. Shuvaev, and M. Strikman.
DVCS amplitude in the parton model.
Phys. Lett. B, 491:96–100, 2000.

[25]
Andrei V. Belitsky and Dieter Mueller.
Twist- three effects in two photon processes.
Nucl. Phys. B, 589:611–630, 2000.

[26]
A.V. Radyushkin and C. Weiss.
DVCS amplitude with kinematical twist - three terms.
Phys. Lett. B, 493:332–340, 2000.

[27]
A.V. Radyushkin and C. Weiss.
DVCS amplitude at tree level: Transversality, twist - three, and
Phys. Rev. D, 63:114012, 2001.

[28]
I. Balitsky and A. Tarasov.
Power corrections to TMD factorization for Z-boson production.
JHEP, 05:150, 2018.

[29]
I. Balitsky and A. Tarasov.
Higher-twist corrections to gluon TMD factorization.
JHEP, 07:095, 2017.

[30]
Daniel Boer and P.J. Mulders.
Time reversal odd distribution functions in leptoproduction.
Phys. Rev. D, 57:5780–5786, 1998.

[31]
Maxim Nefedov and Vladimir Saleev.
Off-shell initial state effects, gauge invariance and angular
distributions in the Drell–Yan process.
Phys. Lett. B, 790:551–556, 2019.

[32]
Maxim A. Nefedov and Vladimir A. Saleev.
High-Energy Factorization for Drell-Yan process in $pp$ and $p{\bar
p}$ collisions with new Unintegrated PDFs.
9 2020.

[33]
M.A. Nefedov, N.N. Nikolaev, and V.A. Saleev.
Drell-Yan lepton pair production at high energies in the Parton
Reggeization Approach.
Phys. Rev. D, 87(1):014022, 2013.

[34]
Simone Marzani.
Combining $Q_T$ and small-$x$ resummations.
Phys. Rev. D, 93(5):054047, 2016.

[35]
I. I. Balitsky and Vladimir M. Braun.
Nonlocal Operator Expansion for Structure Functions of $e^+ e^-$
Phys. Lett., B222:123–131, 1989.

[36]
Ian I. Balitsky and Vladimir M. Braun.
The Nonlocal operator expansion for inclusive particle production in
e+ e- annihilation.
Nucl. Phys., B361:93–140, 1991.

[37]
Ian I. Balitsky and Vladimir M. Braun.
Valleys in Minkowski space and instanton induced cross-sections.
Nucl. Phys., B380:51–82, 1992.

[38]
Francois Gelis and Raju Venugopalan.
Large mass q anti-q production from the color glass condensate.
Phys. Rev., D69:014019, 2004.

[39]
Francois Gelis and Raju Venugopalan.
Particle production in field theories coupled to strong external
Nucl. Phys., A776:135–171, 2006.

[40]
Francois Gelis, Tuomas Lappi, and Raju Venugopalan.
High energy scattering in Quantum Chromodynamics.
Int. J. Mod. Phys., E16:2595–2637, 2007.

[41]
R. D. Tangerman and P. J. Mulders.
Intrinsic transverse momentum and the polarized Drell-Yan process.
Phys. Rev., D51:3357–3372, 1995.

[42]
V.M. Braun and A.N. Manashov.
Operator product expansion in QCD in off-forward kinematics:
Separation of kinematic and dynamical contributions.
JHEP, 01:085, 2012.

[43]
V.M. Braun and A.N. Manashov.
Kinematic power corrections in off-forward hard reactions.
Phys. Rev. Lett., 107:202001, 2011.

[44]
V.M. Braun, Yao Ji, and A.N. Manashov.
Two-photon processes in conformal QCD: Resummation of the
descendants of leading-twist operators.
11 2020.

[45]
Jian Zhou, Feng Yuan, and Zuo-Tang Liang.
Hyperon Polarization in Unpolarized Scattering Processes.
Phys. Rev., D78:114008, 2008.

[46]
A. Bacchetta, Mariaelena Boglione, A. Henneman, and P. J. Mulders.
Bounds on transverse momentum dependent distribution and
fragmentation functions.
Phys. Rev. Lett., 85:712–715, 2000.

[47]
Ignazio Scimemi and Alexey Vladimirov.
Matching of transverse momentum dependent distributions at twist-3.
The European Physical Journal C, 78(10), Oct 2018.

[48]
Vincenzo Barone, Stefano Melis, and Alexei Prokudin.
The Boer-Mulders effect in unpolarized SIDIS: An Analysis of the
COMPASS and HERMES data on the cos 2 phi asymmetry.
Phys. Rev. D, 81:114026, 2010.

[49]
Vincenzo Barone, Stefano Melis, and Alexei Prokudin.
Azimuthal asymmetries in unpolarized Drell-Yan processes and the
Boer-Mulders distributions of antiquarks.
Phys. Rev. D, 82:114025, 2010.

[50]
Albert M Sirunyan et al.
Measurement of the differential Drell-Yan cross section in
proton-proton collisions at $ \sqrt{\mathrm{s}} $ = 13 TeV.
JHEP, 12:059, 2019.

[51]
Georges Aad et al.
Measurement of the transverse momentum distribution of
Drell–Yan lepton pairs in proton–proton collisions at
$\sqrt{s}=13$ TeV with the ATLAS detector.
Eur. Phys. J. C, 80(7):616, 2020.

[52]
Martin Lambertsen and Werner Vogelsang.
Drell-Yan lepton angular distributions in perturbative QCD.
Phys. Rev. D, 93(11):114013, 2016.

[53]
L.Y. Zhu et al.
Measurement of Angular Distributions of Drell-Yan Dimuons in p + p
Interactions at 800-GeV/c.
Phys. Rev. Lett., 102:182001, 2009.

[54]
Ian Balitsky and Giovanni A. Chirilli.
Conformal invariance of transverse-momentum dependent parton
distributions rapidity evolution.
Phys. Rev. D, 100(5):051504, 2019.

[55]
G. Parisi and R. Petronzio.
Small Transverse Momentum Distributions in Hard Processes.
Nucl. Phys. B, 154:427–440, 1979.

[56]
Miguel G. Echevarria, Ignazio Scimemi, and Alexey Vladimirov.
Universal transverse momentum dependent soft function at NNLO.
Phys. Rev., D93(5):054004, 2016.

[57]
Ye Li, Duff Neill, and Hua Xing Zhu.
An exponential regulator for rapidity divergences.
Nucl. Phys. B, 960:115193, 2020.

[58]
Thomas Luebbert, Joel Oredsson, and Maximilian Stahlhofen.
Rapidity renormalized TMD soft and beam functions at two loops.
JHEP, 03:168, 2016.

[59]
Miguel G. Echevarria, Ignazio Scimemi, and Alexey Vladimirov.
Unpolarized Transverse Momentum Dependent Parton Distribution and
Fragmentation Functions at next-to-next-to-leading order.
JHEP, 09:004, 2016.

[60]
Ian Balitsky and Giovanni A. Chirilli.
Next-to-leading order evolution of color dipoles.
Phys. Rev., D77:014019, 2008.

[61]
I. Balitsky and A. Tarasov.
Rapidity evolution of gluon TMD from low to moderate x.
JHEP, 10:017, 2015.

[62]
I. Balitsky and A. Tarasov.
Gluon TMD in particle production from low to moderate x.
JHEP, 06:164, 2016.

[63]
Stefano Forte and Claudio Muselli.
High energy resummation of transverse momentum distributions: Higgs
in gluon fusion.
JHEP, 03:122, 2016.

[64]
S. Arnold, A. Metz, and M. Schlegel.
Dilepton production from polarized hadron hadron collisions.
Phys. Rev. D, 79:034005, 2009.

[28]

author = "Collins, John and Rogers, Ted",
title = "Understanding the large-distance behavior of transverse-momentum-dependent parton densities and the Collins-Soper evolution kernel",
eprint = "1412.3820",
archivePrefix = "arXiv",
primaryClass = "hep-ph",
reportNumber = "YITP-WSB-14-45, JLAB-THY-15-2001",
doi = "10.1103/PhysRevD.91.074020",
journal = "Phys. Rev. D",
volume = "91",
number = "7",
pages = "074020",
year = "2015"