\cheW^{\rm I1S}_{\mu\nu}(x)
~=~{N_c\over 2s}\sum_{f,f'}
\langle A,B|
\nonumber\\
\nonumber\\
\label{wi2s}
\end{eqnarray}

where $\psi=\Psi_1+\Psi_2$ in our approximation.

§.§.§ Annihilation-type terms

Let us start from the “annihilation” part

\begin{eqnarray}
\cheW^{\rm I1San}_{\mu\nu}(x)
~=~{N_c\over 2s}\sum_f e_fc_f
\langle A,B|2a_f[\Bsi_1^f\gamma_\mu\Psi_2^f(x)][\Bsi_2^f\gamma_\nu\Psi_1^f(0)]
\label{wi1sa}\\
\nonumber
\end{eqnarray}

The first term can be copied from photon case (<ref>) so we need to consider the last two terms. Using Fierz transformation (<ref>) we get

\begin{eqnarray}
{N_c\over 2s}\langle A,B|[\Bsi_1\gamma_\mu\Psi_2(x)][\Bsi_2\gamma_\nu\gamma_5\Psi(0)]
\nonumber\\
\big(g_{\mu\nu}g^{\alpha\beta}-\delta_\mu^\alpha\delta_\nu^\beta-\delta_\nu^\alpha\delta_\mu^\beta\big)
{N_c\over 2s}\langle A,B|[\Bsi_1(x)\gamma_\alpha\Psi_1(0)][\Bsi_2(0)\gamma_\nu\gamma_5\Psi_2(x)]
\nonumber\\
|A,B\rangle+x\leftrightarrow0~=~ \cheW_{\mu\nu}^{\rm 5}(x)
\nonumber
\end{eqnarray}

Since we proved in Sect. <ref> that $\cheW_{\mu\nu}^{\rm 5}$ is negligible (see Eq. (<ref>)) we get

\begin{eqnarray}
\cheW^{\rm I1San}_{\mu\nu}(x)
~=~\sum_f e_fc_f a_f\chew^f(x)
\nonumber\\
\Rightarrow~W^{\rm I1San}_{\mu\nu}(q)~=~\sum_f e_fc_f a_f
\big[W^{\rmF f}_{\mu\nu}(q)+W^{\rmH f}_{\mu\nu}(q)+W^{\rmH 2f}_{\mu\nu}(q)\big]
\end{eqnarray}

where $W^{\rmF f}_{\mu\nu}$, $W^{\rmH f}_{\mu\nu}$, and $W^{\rmH 2f}_{\mu\nu}$ are given
by Eqs.  (<ref>),  (<ref>) , and  (<ref>), respectively.

§.§.§ Exchange-type power corrections

Let us consider now the exchange-type power corrections. From Eq. (<ref>) we get

\begin{eqnarray}
\cheW^{\rm I1Sex}_{\mu\nu}(x)
~=~{N_c\over 2s}\sum_f
\langle A,B|(e_{f'}c_fa_f+e_fc_{f'}a_{f'})[\Bsi_1^f\gamma_\mu\Psi_1^f(x)][\Bsi_2^{f'}\gamma_\nu\Psi_2^{f'}(0)]
\label{wi1sa}\\
\nonumber\\
\nonumber\\
=~\sum_{f,f'} \Big(\half(e_{f'}c_fa_f+e_fc_{f'}a_{f'})\chew^{ff'{\rm ex}}_{\mu\nu}(x)
-{e_{f'}c_f\over 2}(\chew^{{\rm S}ff'}_{5\rma})^{\rm ex}_{\mu\nu}(x)
-{e_fc_{f'}\over 2}(\chew^{{\rm S}ff'}_{5\rmb})^{\rm ex}_{\mu\nu}(x)\Big)
\nonumber
\end{eqnarray}

where $\chew^{ff'{\rm ex}}_{\mu\nu}(x)$ is defined by  Eq. (<ref>)
while $(\chew^{{\rm S}ff'}_{5\rm a,b})^{\rm ex}_{\mu\nu}(x)$ are defined in Eq. (<ref>).

\begin{eqnarray}
W^{\rm I1Sex}_{\mu\nu}(q)=\sum_{f,f'} \Big[{e_{f'}c_fa_f+e_fc_{f'}a_{f'}\over 2}W^{ff'{\rm ex}}_{\mu\nu}(q)
-{e_{f'}c_f\over 2}(W^{\rmS ff'}_{5\rma})^{\rm ex}_{\mu\nu}(q)
-{e_fc_{f'}\over 2}(W^{\rmS ff'}_{5\rmb})^{\rm ex}_{\mu\nu}(q)\Big]
\nonumber\\
=~{N_c\over 2(N_c^2-1)Q_\parallel^2}\!\int\! d^2k_\perp\sum_{f,f'} \Big[(e_{f'}c_fa_f+e_fc_{f'}a_{f'})
\Big([k^\perp_\mu(q-k)^\perp_\nu+\mu\leftrightarrow\nu+g^\perp_{\mu\nu}(k,q-k)_\perp]
\nonumber\\
\times~J^{1ff'}_{--}(q,k_\perp)
-~g^\perp_{\mu\nu} (k,q-k)_\perp J^{2ff'}_{--}(q,k_\perp)\Big)
\label{w1sexresult}\\
-~e_{f'}c_f[\epsilon_{\mu m}k^m(q-k)_\nu+\mu\leftrightarrow\nu]
[k_\mu\epsilon_{\nu n}(q-k)^n+\mu\leftrightarrow\nu]
\nonumber
\end{eqnarray}

where $W^{ff'{\rm ex}}_{\mu\nu}(q)$, $(W^{\rmS ff'}_{5\rma})^{\rm ex}_{\mu\nu}(q)$, and
$(W^{\rmS ff'}_{5\rmb})^{\rm ex}_{\mu\nu}(q)$
are given in Eqs. (<ref>), (<ref>), and (<ref>), respectively.

§.§.§ The result for the symmetric part of $W^{\rm I1}_{\mu\nu}(x)$

The result for $W^{\rm I1S}_{\mu\nu}(q)$ can be represented as a sum of “annihilation” and “exchange” parts

\begin{eqnarray}
W^{\rm I1S}_{\mu\nu}(q)~=~W^{\rm I1San}_{\mu\nu}(q)~+~W^{\rm I1Sex}_{\mu\nu}(q)
\nonumber\\
W^{\rm I1San}_{\mu\nu}(q)~
=~\sum_f e_fc_f a_f
\big[W^{\rm F f}_{\mu\nu}(q)+W^{\rm H f}_{\mu\nu}(q)+W^{\rm H 2f}_{\mu\nu}(q)\big]
\nonumber\\
W^{\rm I1Sex}_{\mu\nu}(q)~=~{\rm r.h.s.~of~Eq.~(\ref{w1sexresult})}
\label{wi1s}
\end{eqnarray}

where $W^{\rm F f}_{\mu\nu}(q)$ is given by Eq. (<ref>),
$W^{\rm H f}_{\mu\nu}(q)$ by Eq. (<ref>), and $W^{\rm H2 f}_{\mu\nu}(q)$ by Eq. (<ref>).

§.§ Antiymmetric part of interference hadronic tensor $W^{\rm I1}$

From definitions (<ref>) and (<ref>) we get

\begin{eqnarray}
\cheW^{\rm I1A}_{\mu\nu}(x)
~=~{N_c\over 2s}\sum_{f,f'}
\langle A,B|
\nonumber\\
\nonumber\\
\label{wi2s}
\end{eqnarray}

where $\psi=\Psi_1+\Psi_2$ in our approximation.

Let us start from the annihilation part. After Fierz transformations (<ref>) and (<ref>) we get

\begin{eqnarray}
\cheW^{\rm I1Aan}_{\mu\nu}(x)
~=~{N_c\over 2s}\sum_f e_fc_f
\langle A,B|2a_f[\Bsi_1^f\gamma_\mu\Psi_2^f(x)][\Bsi_2^f\gamma_\nu\Psi_1^f(0)]
\label{wi1sa}\\
\nonumber\\
=~\sum_f e_fc_f\big[ia_f\epsilon_{\mu\nu}^{~~\alpha\beta}\chepizw_{\alpha\beta}^{5f}(x)
+a_f\cheW^{{\rm as}f}_{\mu\nu}(x)
-{i\over 2}\epsilon_{\mu\nu}^{~~\alpha\beta}\chepizw_{\alpha\beta}^{\calf f}(x)\big]
\nonumber
\end{eqnarray}

where we used definitions (<ref>). As demonstrated in Sect. <ref>,
$\chepizw_{\mu\nu}^{5f}=\cheW^{{\rm as}f}_{\mu\nu}=0$ with our accuracy
and $\pizw_{\alpha\beta}^{\calf f}(q)$ is given by Eq. (<ref>) so we obtain

\begin{eqnarray}
W^{\rm I1Aan}_{\mu\nu}(q)
~=~-i\epsilon_{\mu\nu\alpha\beta}q^\alpha\!\int\! d^2k_\perp~(\tilq+q_\perp-2k_\perp)^\beta
\sum_f e_fc_f\pizf^f(q,k_\perp)
\end{eqnarray}

The exchange-type power corrections are

\begin{eqnarray}
\cheW^{\rm I1Aex}_{\mu\nu}(x)~=~{N_c\over 2s}\sum_{f,f'}
\langle A,B|
\nonumber\\
\\
-\mu\leftrightarrow\nu-x\leftrightarrow 0
\nonumber\\
=~\sum_{f,f'}\Big[{e_{f'}c_fa_f+e_fc_{f'}a_{f'}\over 2}(\chew^{\rmA ff'})^{\rm ex}_{\mu\nu}(x)
-~{e_{f'}c_f\over 2}(\chew^{{\rm A}ff'}_{5\rma})^{\rm ex}_{\mu\nu}(x)
-{e_fc_{f'}\over 2}(\chew^{{\rm A}ff'}_{5\rmb})^{\rm ex}_{\mu\nu}(x)\Big]
\nonumber
\end{eqnarray}

where we used Eq. (<ref>). In the momentum space this gives

\begin{eqnarray}
W^{\rm I1Aex}_{\mu\nu}(q)~=~
{iN_c\over 2(N_c^2-1)Q_\parallel^2}\sum_{f,f'}\!\int\! d^2k_\perp
\Big\{c_fe_{f'}
[\epsilon_{\mu m}k^m(q-k)_\nu-\mu\leftrightarrow\nu]
\nonumber\\
-~c_{f'}e_f[k_\mu\epsilon_{\nu n}(q-k)^n
-\mu\leftrightarrow\nu]  J^{2ff'}_{-+}\Big\}
\label{wi1aex}
\end{eqnarray}

where we used Eqs.  (<ref>) and (<ref>).

§.§.§ The result for antisymmetric part of interference hadronic tensor $W^{\rm I1}$

As usual, the result consists of “annihilation” and “exchange” parts

\begin{eqnarray}
W^{\rm I1A}_{\mu\nu}(q)~=~W^{\rm I1Aan}_{\mu\nu}(q)+W^{\rm I1Aex}_{\mu\nu}(q),
\nonumber\\
W^{\rm I1Aan}_{\mu\nu}(q)~=~-i\epsilon_{\mu\nu\alpha\beta}q^\alpha\!\int\! d^2k_\perp~(\tilq+q_\perp-2k_\perp)^\beta\sum_f e_fc_f\pizf^f(q,k_\perp),
\nonumber\\
W^{\rm I1Aex}_{\mu\nu}(q)~=~{\rm r.h.s. ~of~Eq.~(\ref{wi1aex})}
\label{wi1asy}
\end{eqnarray}

where $\pizf^f(q,k_\perp)$ is given by Eq. (<ref>). Obviously,
$q^\mu W^{\rm I1Aan}_{\mu\nu}(q)=0$.

§.§ Symmetric part of interference hadronic tensor $W^{\rm I2}$

From definitions (<ref>) and (<ref>) we get

\begin{eqnarray}
\cheW^{\rm I2S}_{\mu\nu}(x)
~=~{N_c\over 2s}\sum_{f,f'}\Big((e_{f'}c_fa_f-e_fc_{f'}a_{f'})
\langle A,B|[\bsi(x)\gamma_\mu\psi(x)]^f[\bsi(0)\gamma_\nu\psi(0)]^{f'}
\nonumber\\
+~e_fc_{f'}\langle A,B|[\bsi(x)\gamma_\mu\psi(x)]^f[\bsi(0)\gamma_\nu\gamma_5\psi(0)]^{f'}
\nonumber\\
-~e_{f'}c_f\langle A,B|[\bsi(x)\gamma_\mu\gamma_5\psi(x)]^f[\bsi(0)\gamma_\nu\psi(0)]^{f'}|A,B\rangle\Big)
\label{wi2s}
\end{eqnarray}

where $\psi=\Psi_1+\Psi_2$ in our approximation.

§.§.§ Annihilation-type power corrections

Let us start with annihilation-type power corrections. Since in this case $f=f'$, the first term in the r.h.s.
of Eq. (<ref>) vanishes and the second turns to

\begin{eqnarray}
\chew^{\rm I2San}_{\mu\nu}(x)~=~\sum_{f}e_fc_f\chew^{\rm If}_{\mu\nu}(x)
\nonumber\\
\chew^{\rm If}_{\mu\nu}(x)
~=~{N_c\over 2s}\langle A,B|[\Bsi_1(x)\gamma_\mu\Bsi_2(x)][\Bsi_2(0)\gamma_\nu\gamma_5\Bsi_1(0)]
\nonumber\\
+\mu\leftrightarrow\nu|A,B\rangle-x\leftrightarrow 0
\label{wisa}
\end{eqnarray}

After Fierz transformation (<ref>) it turns to

\begin{eqnarray}
\cheW^{\rm I}_{\mu\nu}(x)~=~
-{N_c\over 4s} \langle p_A,p_B|[\Bsi_1^m(x)\sigma_{\mu\xi}\Psi^n_1(0)][\Bsi_2^n(0)\sigma_\nu^{~\xi}\gamma_5\Psi_2^m(x)]
\label{WI1}\\
\nonumber\\
{N_cg_{\mu\nu}\over 2s} \langle p_A,p_B|[\Bsi_1^m(x)\Psi^n_1(0)][\Bsi_2^n(0)\gamma_5\Psi_2^m(x)]
\nonumber\\
\nonumber
\end{eqnarray}

where we suppressed flavor label.

Let us start from the second term. Obviously, the leading-twist contribution vanishes with our accuracy.
Next, consider one-gluon terms and start from

\begin{eqnarray}
{N_cg_{\mu\nu}\over 2s} \langle p_A,p_B|\big([\bsi_A^m(x)\Xi^n_1(0)]+[\Bxi_1^m(x)\Psi^n_1(0)\big)[\bsi_B^n(0)\gamma_5\psi_B^m(x)]
\\
\nonumber\\
=~{ig_{\mu\nu}\over 2s^2} \big\{
\langle \bsi(x)\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A\langle\bsi B^i(0)\gamma_5\psi(x)\rangle_B
-\langle \big(\bsi{1\over\alpha}\big)(x)\sigma_{\star i}\psi(0)\rangle_A\langle\bsi(0) \gamma_5B^i\psi(x)\rangle_B
\nonumber\\
-~\langle \bsi(x)\sigma_{\star i}\gamma_5{1\over\alpha}\psi(0)\rangle_A\langle\bsi B^i(0)\psi(x)\rangle_B
+\langle \big(\bsi{1\over\alpha}\big)(x)\sigma_{\star i}\gamma_5\psi(0)\rangle_A\langle\bsi(0) B^i\psi(x)\rangle_B\big\}
\nonumber
\end{eqnarray}

The projectile matrix element can bring one factor of $s$ while the target one cannot so the contribution to
$\cheW^{\rm I}_{\mu\nu}(x)$ 	is $\sim{g_{\mu\nu}\over \alpha_qs}$ which is $O(\beta_q)\times{q_{\mu\nu}\over Q_\parallel^2}$.
Similarly, the contributions coming from $\Xi_2$ and $\Bxi_1$ are  $\sim{g_{\mu\nu}\over \beta_qs}$ and can be neglected.

Let us consider now two-gluon term coming from $\Xi_1$ and $\Xi_2$

\begin{eqnarray}
{N_cg_{\mu\nu}\over 2s} \langle p_A,p_B|[\bsi_A^m(x)\Xi^n_1(0)][\bsi_B^n(0)\gamma_5\Xi_2^m(x)]
\nonumber\\
\nonumber\\
=~{g_{\mu\nu}\over 2s^3} \big\{
\langle \bsi A^j(x)\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A\langle\bsi B^i(0)\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
-~\langle \bsi A^j(x)\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A\langle\bsi B^i(0)\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B\big\}
\nonumber\\
=~{g_{\mu\nu}\over 2s^3} \big\{
\langle \bsi(x) (A^j(x)\sigma_{\star i}-A^i(x)\sigma_{\bu j}){1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle \bsi A^i(x)\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A\langle\bsi B^j(0)\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B\big\}
\nonumber\\
=~-{g_{\mu\nu}\over 2s^3} \langle \bsi(x)\notA(x)\notp_2{1\over\alpha}\psi(0)\rangle_A
\langle\bsi(0) \notB(0)\notp_1\gamma_5{1\over\beta}\psi(x)\rangle_B
\label{WI4}
\end{eqnarray}

where again we used Eq. (<ref>) without two last terms and Eq. (<ref>).
Now, from equations of motion (<ref>) and (<ref>) we see that the target matrix element vanishes:

\begin{eqnarray}
\int\! dx~e^{-i\beta_q x_\star+i(k,x)_\perp}\langle\bsi(0) \notB(0)\notp_1\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
=~{k^\perp_i\over\beta_q}\!\int\! dx~e^{-i\beta x_\star+i(k,x)_\perp}\langle\bsi(0)\gamma^i\notp_1\gamma_5\psi(x)\rangle_B
\label{b5vanishes}
\end{eqnarray}

It is easy to see that the  term coming from $\Bxi_1$ and $\Bxi_2$ vanishes for the same reason.

\begin{eqnarray}
{N_cg_{\mu\nu}\over 2s} \langle p_A,p_B|[\Bxi_1^m(x)\Psi^n_1(0)][\Bxi_2^n(0)\gamma_5\psi_B^m(x)]
\label{WI6}\\
\nonumber\\
=~-{g_{\mu\nu}\over 2s^3} \langle (\bsi{1\over\alpha}\big)(x)\notp_2\notA(0)\psi(0)\rangle_A
\langle\big(\bsi{1\over\beta}\big)(0)
\notp_1\notB(x)\gamma_5\psi(x)\rangle_B
\nonumber
\end{eqnarray}

Next, two-gluon contribution from $\Bxi_1$ and $\Xi_1$ vanishes since $\Bxi_1 \Xi_1=0$, and  similarly
$\Bxi_2\Xi_2=0$. Finally, the contribution coming from $\Bxi_1$ and $\Xi_2$ is $\sim {1\over N_c^2}$ as demonstrated in Sect.
<ref> (see Eq. (<ref>)), so the second term in Eq. (<ref>) vanishes and we get

\begin{eqnarray}
\cheW^{\rm I}_{\mu\nu}(x)~=~
-{N_c\over 4s} \langle p_A,p_B|[\Bsi_1^m(x)\sigma_{\mu\xi}\Psi^n_1(0)][\Bsi_2^n(0)\sigma_\nu^{~\xi}\gamma_5\Psi_2^m(x)]
\label{WI7}\\
\nonumber
\end{eqnarray}

This contribution is similar to Eq. (<ref>) up to extra $\gamma_5$ and relative signs.
As we discussed above, extra
$\gamma_5$ cannot change the power of our small parameters ${q_\perp^2\over Q^2}$ and $\alpha_q,\beta_q$ so we can consider only terms
which gave leading contributions to $V^{\rm H}_{\mu\nu}$.

The leading-twist contribution

\begin{eqnarray}
\cheW^{\rm I,l.t.}_{\mu\nu}(x)~=~
-{1\over 4s} \langle p_A,p_B|[\bsi_A(x)\sigma_{\mu\xi}\psi_A(0)][\bsi_B(0)\sigma_\nu^{~\xi}\gamma_5\psi_B(x)]
\label{WI1l}\\
\nonumber
\end{eqnarray}

is easily obtained from parametrizations (<ref>)

\begin{eqnarray}
W^{\rm I,lt}_{\mu\nu}(x)(q)
~=~{1\over 16\pi^4}\!\int\! dx_\bu dx_\star d^2x_\perp~e^{-i\alpha_qx_\bu-i\beta_qx_\star+i(q,x)_\perp}\cheW^{\rm I,lt}_{\mu\nu}(x)
\nonumber\\
~=~-{i\over 2m^2}\!\int\! d^2k_\perp
\big(\epsilon_{\mu j}[k^j(q-k)^\perp_\nu+(q-k)_j k^\perp_\nu]+\mu\leftrightarrow\nu\big)\pizh^f(q,k_\perp)
\label{WIlt}
\end{eqnarray}

\begin{eqnarray}
\pizh^f(q,k_\perp)~=~h^\perp_{1}(\alpha_q,k_\perp)\barh^\perp_{1}(\beta_q,(q-k)_\perp) ~-~h_{1f}^\perp\leftrightarrow\barh_{1f}^\perp
\label{pizh}
\end{eqnarray}

As usually, the term with $h_{1f}^\perp\leftrightarrow\barh_{1f}^\perp$ comes from $x\leftrightarrow0$ contribution.

Next, we need to consider terms in Eq. (<ref>) with one or two gluon operators and generalize the calculations from Sect.  <ref>
to our case.

§.§.§ One-gluon terms

Let us first consider term coming from  $\Xi_1(0)={i\over s}\sigma_{\star i}B^i{1\over \alpha}\psi_A$. Separating color-singlet
matrix elements in Eq. (<ref>), we get

\begin{eqnarray}
\cheW^{\rm 1I(1)}_{1\mu\nu}(x)~=~
-{i\over 4s^2}\big\{\langle\bsi(x)\sigma_{\mu\xi}\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_\nu^{~\xi}\gamma_5\psi(x)]\rangle_B
\label{WI9}\\
\nonumber
\end{eqnarray}

As we mentioned, this contribution is similar to the one considered in Sec. <ref> so we
need to take only the case
of longitudinal $\mu$ and transverse $\nu$, or vice versa - all other cases are power-suppressed. We get

\begin{eqnarray}
\cheW^{\rm 1I(1)}_{1\mu\nu}(x)~=~
-{ip_{2\mu}\over 2s^3}\big\{
\langle\bsi(x)\sigma_{\bu\xi}\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_{\nu_\perp}^{~\xi}\gamma_5\psi(x)]\rangle_B
\nonumber\\
+~\langle\bsi(x)\sigma_{\nu_\perp\xi}\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_\bu^{~\xi}\gamma_5\psi(x)]\rangle_B
\nonumber\\
\nonumber\\
-{ip_{2\mu}\over 2s^3}\big\{
\langle\bsi(x)\sigma_{\nu_\perp j}\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_\bu^{~j}\gamma_5\psi(x)]\rangle_B
\nonumber\\\
\langle\bsi(x)\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_{\nu_\perp\bu}\gamma_5\psi(x)]\rangle_B
\nonumber\\\
\langle\bsi(x)\sigma_{\bu j}\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_{\nu_\perp}^{~j}\gamma_5\psi(x)]\rangle_B
\nonumber\\
+~{2\over s}\langle\bsi(x)\sigma_{\nu_\perp\bu}\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_{\bu\star}\gamma_5\psi(x)]\rangle_B
\nonumber\\
\label{WI10}
\end{eqnarray}

where we used power-counting results from
Ref. [16]  to eliminate terms proportional to $p_{1\mu}$, cf Eq. (<ref>).
Moreover, similarly to Eq. (<ref>) case,  Eq. (<ref>) shows that two last terms in the r.h.s. of Eq. (<ref>)
are small and therefore

\begin{eqnarray}
\cheW^{\rm 1I(1)}_{1\mu\nu}(x)~=~
-{p_{2\mu}\over 2s^3}\big\{
\langle\bsi(x)\big[g_{i\nu_\perp}\sigma_{\star j}-g_{ij}\sigma_{\star \nu_\perp}\big]{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_\bu^{~j}\gamma_5\psi(x)\rangle_B
\nonumber\\\
\langle\bsi(x)\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_{\bu\nu_\perp}\gamma_5\psi(x)]\rangle_B
\nonumber\\\
\nonumber\\\
=~-{p_{2\mu}\over 2s^3}\big\{\langle\bsi(x)\sigma_{\star j}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi(0)[B_\nu(0)\sigma_\bu^{~j}-B^j(0)\sigma_{\bu\nu_\perp}]\gamma_5\psi(x)\rangle_B
\nonumber\\\
\langle\bsi(x)\sigma_{\star \nu_\perp}\gamma_5{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_{\bu i}\psi(x)\rangle_B
\big\}
\label{WI11}
\end{eqnarray}

where we used
Eqs. (<ref>), (<ref>), and (<ref>).
Next, the target matrix element in the first term in the r.h.s. can be rewritten as

\begin{eqnarray}
\langle\bsi(0)[B_\nu(0)\sigma_\bu^{~j}-B^j(0)\sigma_{\bu\nu_\perp}]\gamma_5\psi(x)]\rangle_B
~=~\epsilon_{\nu j}\langle\bsi(0)\notB(0)\notp_1\psi(x)\rangle_B
\end{eqnarray}

so we get

\begin{eqnarray}
\cheW^{\rm 1I(1)}_{1\mu\nu}(x)~
=~-{p_{2\mu}\over s^3}\epsilon_{\nu j}\langle\bsi(x)\sigma_{\star j}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi(0)\notB(0)\notp_1\psi(x)\rangle_B
\label{WI13}
\end{eqnarray}

where we used $\sigma_{\star \nu_\perp}\gamma_5=i\epsilon_{\nu j}\sigma_\star^{~j}$.
Now, from equation of motion (<ref>) and parametrization (<ref>) we get
the corresponding contribution to
$\cheW^{\rm 1I}_{\mu\nu}(x)$ in the form

\begin{eqnarray}
\cheW^{\rm 1I(1)}_{1\mu\nu}(x)~
=~-i\epsilon_{\nu j}{p_{2\mu}\over \alpha_qsm^2}\!\int\! d^2k_\perp~k^j_\perp(q-k)_\perp^2\pizh(q,k_\perp)
\label{WI14}
\end{eqnarray}

The term in Eq. (<ref>) coming from $\Bxi_1$ is similar:

\begin{eqnarray}
\cheW^{\rm 1I(1)}_{2\mu\nu}(x)~=~
{i\over 4s^2}\big\{\langle\big(\bsi{1\over\alpha}\big)(x)\sigma_{\mu\xi}\sigma_{\star i}\psi(0)\rangle_A
\langle\bsi(0)\sigma_\nu^{~\xi}\gamma_5B^i(x)\psi(x)]\rangle_B
\nonumber\\
\nonumber\\
{ip_{2\mu}\over 2s^3}\big\{
\langle\big(\bsi{1\over\alpha}\big)(x)\sigma_{\star i}\sigma_{\bu\xi}\psi(0)\rangle_A
\langle\bsi B^i(x)\sigma_{\nu_\perp}^{~\xi}\gamma_5\psi(x)]\rangle_B
\nonumber\\
+~\langle\big(\bsi{1\over\alpha}\big)(x)\sigma_{\star i}\sigma_{\nu_\perp\xi}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(x)\sigma_\bu^{~\xi}\gamma_5\psi(x)]\rangle_B
\nonumber\\
\nonumber\\
{ip_{2\mu}\over 2s^3}\big\{
\langle\big(\bsi{1\over\alpha}\big)(x)\sigma_{\star i}\sigma_{\nu_\perp j}\psi(0)\rangle_A
\langle\bsi B^i(x)\sigma_\bu^{~j}\gamma_5\psi(x)]\rangle_B
\nonumber\\\
\langle\big(\bsi{1\over\alpha}\big)(x)\sigma_{\star i}{1\over\alpha}\psi(0)\rangle_A
\langle\bsi B^i(x)\sigma_{\bu\nu_\perp}\gamma_5\psi(x)]\rangle_B
\nonumber\\\
\nonumber\\
-{p_{2\mu}\over 2s^3}\big\{
\langle\big(\bsi{1\over\alpha}\big)(x)\big[g_{i\nu_\perp}\sigma_{\star j}-g_{ij}\sigma_{\star \nu_\perp}\big]\psi(0)\rangle_A
\langle\bsi B^i(x)\sigma_\bu^{~j}\gamma_5\psi(x)\rangle_B
\nonumber\\\
\langle\bsi{1\over\alpha}\big)(x)\sigma_{\star i}\psi(0)\rangle_A
\langle\bsi B^i(x)\sigma_{\bu\nu_\perp}\gamma_5\psi(x)]\rangle_B
\nonumber\\\
\nonumber\\\
=~-{p_{2\mu}\over 2s^3}\big\{\langle\bsi{1\over\alpha}\big)(x)\sigma_{\star j}\psi(0)\rangle_A
\langle\bsi(0)[B_\nu(x)\sigma_\bu^{~j}-B^j(x)\sigma_{\bu\nu_\perp}]\gamma_5\psi(x)\rangle_B
\nonumber\\\
\langle\big(\bsi{1\over\alpha}\big)(x)\sigma_{\star \nu_\perp}\gamma_5\psi(0)\rangle_A
\langle\bsi B^i(x)\sigma_{\bu i}\psi(x)\rangle_B
\big\}
\nonumber\\\
=~{p_{2\mu}\over s^3}\epsilon_{\nu j}\langle\big(\bsi{1\over\alpha}\big)(x)\sigma_{\star j}\psi(0)\rangle_A
\langle\bsi(0)\notp_1\notB(x)\psi(x)\rangle_B
\label{WI15}
\end{eqnarray}

where we left only the terms similar to the leading terms in Eq. (<ref>)
and made the same transformations.

Now, from equation of motion (<ref>) and parametrization (<ref>) we get
the corresponding contribution to
$\cheW^{\rm 1I}_{\mu\nu}(x)$ in the form

\begin{eqnarray}
\cheW^{\rm 1I(1)}_{2\mu\nu}(x)~
=~-i\epsilon_{\nu j}{p_{2\mu}\over \alpha_qsm^2}\!\int\! d^2k_\perp~k^j_\perp(q-k)_\perp^2\pizh(q,k_\perp)
\label{WI15a}
\end{eqnarray}

so it doubles the contribution (<ref>).

Let us now consider term in Eq. (<ref>) coming form $\Xi_2$.
For longitudinal $\mu$ and transverse $\nu$ we get

\begin{eqnarray}
\cheW^{\rm 1I(2)}_{1\mu\nu}(x)~=~
-{i\over 4s^2}\big\{\langle\bsi A^i(x)\sigma_{\mu\xi}\psi(0)\rangle_A
\langle\bsi(0)\sigma_\nu^{~\xi}\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\label{WI16}\\
\nonumber\\
=~-{ip_{1\mu}\over 2s^3}\big\{\langle\bsi A^i(x)\sigma_{\star\xi}\psi(0)\rangle_A
\langle\bsi(0)\sigma_{\nu_\perp}^{~\xi}\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle\bsi A^i(x)\sigma_{\nu_\perp\xi}\psi(0)\rangle_A
\langle\bsi(0)\sigma_{\star}^{~\xi}\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
-~\psi(0)\otimes\gamma_5\psi(x)\leftrightarrow\gamma_5\psi(0)\otimes\psi(x)\big\}  ~-~x\leftrightarrow 0
\nonumber\\
=~-{ip_{1\mu}\over 2s^3}\big\{\langle\bsi A^i(x)\sigma_{\star j}\psi(0)\rangle_A
\langle\bsi(0)\sigma_{\nu_\perp}^{~j}\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
-~i\langle\bsi A^i(x)\sigma_{\nu_\perp\star}\psi(0)\rangle_A
\langle\bsi(0)\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle\bsi A^i(x)\sigma_{\nu_\perp j}\psi(0)\rangle_A
\langle\bsi(0)\sigma_{\star}^{~j}\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
+~{2\over s}\langle\bsi A^i(x)\sigma_{\star\bu}\psi(0)\rangle_A
\nonumber\\
\times~
\langle\bsi(0)\sigma_{\nu_\perp\star}\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
-~\psi(0)\otimes\gamma_5\psi(x)\leftrightarrow\gamma_5\psi(0)\otimes\psi(x)\big\}   ~-~x\leftrightarrow 0
\nonumber
\end{eqnarray}

where we neglected contribution $\sim p_{2\mu}$ since it is $\sim~p_{2\mu}q^\perp_\nu{m_\perp^2\over\beta_qs^2}$, see power counting (without $\gamma_5$) in Ref. [16] .
Next, similarly to Eq. (<ref>),  the last two terms in the r.h.s. of Eq. (<ref>) are
$\sim O\big({m_\perp^2\over s}\big)$ in comparison to the first two terms so we get

\begin{eqnarray}
\cheW^{\rm 1I(2)}_{1\mu\nu}(x)~
=~-{p_{1\mu}\over 2s^3}\big\{\langle\bsi A^i(x)\sigma_{\star j}\psi(0)\rangle_A
\langle\bsi(0)[g_{i\nu_\perp}\sigma_{\bu j}-g_{ij}\sigma_{\bu \nu_\perp}]\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
-~\langle\bsi A^i(x)\sigma_{\star\nu_\perp}\psi(0)\rangle_A
\langle\bsi(0)\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
-~\psi(0)\otimes\gamma_5\psi(x)\leftrightarrow\gamma_5\psi(0)\otimes\psi(x)\big\}   ~-~x\leftrightarrow 0
\nonumber\\
=~{p_{1\mu}\over s^3}\big\{\langle\bsi(x)\big(A_\nu(x)\sigma_{\star j}-A_j(x)\sigma_{\star \nu_\perp}\big)\gamma_5\psi(0)\rangle_A
\langle\bsi(0)\sigma_\bu^{~ j}{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle\bsi A^i(x)\sigma_{\star i}\psi(0)\rangle_A
\langle\bsi(0)\sigma_{\bu\nu_\perp}\gamma_5{1\over\beta}\psi(x)\rangle_B\big\}  ~-~x\leftrightarrow 0
\nonumber\\
=~-{ip_{1\mu}\over s^3}\epsilon_{\nu_\perp j}\langle\bsi A^i(x)\sigma_{\star i}\psi(0)\rangle_A\langle\bsi(0)\sigma_\bu^{~ j}{1\over\beta}\psi(x)\rangle_B
\label{WI18}
\end{eqnarray}

From Eq. (<ref>) we get the corresponding contribution to $W^{\rm I}_{\mu\nu}(q)$ in the form

\begin{eqnarray}
\cheW^{\rm 1I(2)}_{1\mu\nu}(q)~=~-i\epsilon_{\nu j}{p_{1\mu}\over \beta_qsm^2}\!\int\! d^2k_\perp~(q-k)^j_\perp k_\perp^2\pizh(q,k_\perp)
\label{WI19}
\end{eqnarray}

Similarly to Eq. (<ref>) one can demonstrate that the term coming from  $\Bxi_2$ doubles the result (<ref>) of $\Xi_2$
so we finally get

\begin{equation}
\hspace{-0mm}
W^{\rm 1I}_{\mu\nu}(q)
=~-i{2\epsilon_{\nu j}\over Q_\parallel^2m^2}\!\int\! d^2k_\perp~[\beta p_{2\mu}k^j_\perp(q-k)_\perp^2+\alpha p_{1\mu}(q-k)^jk_\perp^2]\pizh(q,k_\perp)
\label{WI20}
\end{equation}

where we have added the contribution of transverse $\mu$ and longitudinal $\nu$.

§.§.§ Two-gluon terms

Let us start from the contribution to $\cheW^{\rm I}_{\mu\nu}(x)$ of Eq. (<ref>) coming from $\Xi_{A}$ and $\Xi_{B}$.
After separation of color-singlet matrix elements, it takes the form

\begin{eqnarray}
\cheW^{\rm 2I(1)}_{\mu\nu}(x)~=~
{1\over 4s^3}\big\{\langle\bsi A^i(x)\sigma_{\mu\xi}\sigma_{\star j}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^j(0)\sigma_\nu^{~\xi}\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\label{WI21}\\
\nonumber
\end{eqnarray}

This equation resembles $\cheV_{1\mu\nu}^{\rm (2a)H}(x)$ of
Eq. (<ref>) calculated in Sect. <ref> so we will
use power counting from that Section and
consider only two transverse or two longitudinal indices.

First, let us consider  transverse $\mu$ and $\nu$.

\begin{eqnarray}
\cheW^{\rm 2I(1)}_{\mu\nu}(x)~=~
{1\over 4s^3}\big\{\langle\bsi A^i(x)\sigma_{\mu_\perp k}\sigma_{\star j}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^j(0)\sigma_{\nu_\perp}^{~k}\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\label{WI22}\\
+~{2\over s}\langle\bsi A^i(x)\sigma_{\mu_\perp\bu}\sigma_{\star j}\psi(0)\rangle_A
\langle\bsi B^j(0)\sigma_{\nu_\perp\star}\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
\nonumber
\end{eqnarray}

Similarly to Eq. (<ref>),  from Eqs.  (<ref>) and  (<ref>),
we see that the second term in the r.h.s.
can be neglected in comparison to the first one and after some algebra  we get

\begin{eqnarray}
\cheW^{\rm 2I(1)}_{\mu_\perp\nu_\perp}(x)~
\nonumber\\
{-1\over 4s^3}\big\{\langle\bsi A^i(x)[g_{\mu j}\sigma_{\star k}-g_{jk}\sigma_{\star\mu_\perp}]{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^j(0)[g_{\nu i}\sigma_\bu^{~k}-\delta_i^k\sigma_{\bu\nu_\perp}]\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
\nonumber\\
=~{1\over 2s^3}\Big[\langle\bsi(x)[A_\nu(x)\sigma_{\star i}-A_i(x)\sigma_{\star \nu}]{1\over \alpha}\psi(0)\rangle_A
\langle\bsi (B_\mu(0)\sigma_\bu^{~i}-B^i(0)\sigma_{\bu \mu})\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle\bsi A_\nu(x)\sigma_{\star k}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B(\delta_i^k\delta_\mu^j-\delta_i^j\delta_\mu^k)
\nonumber\\
+~\langle\bsi A^i(x)\sigma_{\star j}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B_\mu(0)\sigma_{\bu k}\gamma_5{1\over\beta}\psi(x)\rangle_B(\delta_i^k\delta_\nu^j-\delta_\nu^k\delta_i^j)
\nonumber\\
~+~\mu\leftrightarrow\nu~-~x\leftrightarrow 0
\label{WI23}
\end{eqnarray}

Moreover, it is easy to see that the first term in the r.h.s. can be omitted:
either projectile or target matrix element vanishes due to Eq. (<ref>). For the next two terms in the r.h.s. of Eq. (<ref>) we use

\begin{eqnarray}
{1\over 2s^3}\big\{\langle\bsi A_\nu(x)\sigma_{\star k}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^i(0)\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B(\delta_i^k\delta_\mu^j-\delta_i^j\delta_\mu^k)
\nonumber\\
\nonumber\\
=~{1\over 2s^3}\big\{\langle\bsi(x)A_\nu(x)\sigma_\star^{~ i}\gamma_5{1\over\alpha}\psi(0)\rangle_A
\langle\bsi [B_\mu(0)\sigma_{\bu i}-B_i(0)\sigma_{\bu\mu_\perp}]{1\over\beta}\psi(x)\rangle_B
\nonumber\\
-~\langle\bsi A_\nu(x)\sigma_{\star\mu_\perp}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^j(0)\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~{s\over 4}\langle\bsi A_\nu(x)\sigma_{mn}\gamma_5{1\over \alpha}\psi(0)\rangle_A \langle\bsi B_\mu(x)\sigma^{mn}{1\over\beta}\psi(x)\rangle_B
\nonumber\\
-~{s\over 2}\langle\bsi(x)A_\nu(x)\sigma_{\mu_\perp}^{~ j}\gamma_5{1\over\alpha}\psi(0)\rangle_A\big\}
\langle\bsi B_i(x)\sigma^{ij}{1\over\beta}\psi(x)\rangle_B~=~0
\label{WI24}
\end{eqnarray}

since the target matrix element in the first term vanishes due to Ee. (<ref>), the one in the second term due to Eq. (<ref>),
and the last two terms are small by power counting, $\sim {q_\perp^2\over s}\big({ q^\perp_\mu q^\perp_\nu\over Q_\parallel^2}\big)$.

\begin{eqnarray}
{2\over s^3}\big\{\langle\bsi A^i(x)\sigma_{\star j}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B_\mu(0)\sigma_{\bu k}\gamma_5{1\over\beta}\psi(x)\rangle_B(\delta_i^k\delta_\nu^j-\delta_\nu^k\delta_i^j)
\nonumber\\
\nonumber\\
=~\langle\bsi(x)[A_i(x)\sigma_{\star \nu}-A_\nu\sigma_{\star i}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B_\mu(0)\sigma_\bu^{~i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle\bsi A^i(x)\sigma_{\star i}\gamma_5{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B_\mu(0)\sigma_{\bu\nu_\perp}{1\over\beta}\psi(x)\rangle_B
\nonumber\\
-~{s\over 4}\langle\bsi(x)A_\nu\sigma_{mn}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B_\mu(0)\sigma^{mn}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~{s\over 2}\langle\bsi(x)A^i\sigma_{ij}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B_\mu(0)\sigma_{\nu_\perp}^{~j}\gamma_5{1\over\beta}\psi(x)\rangle_B~=~0
\end{eqnarray}

for the same reason: the projectile matrix element in the first term in the r.h.s. vanishes due to Eq. (<ref>),
the one in the second term due to Eq. (<ref>),
and the last two terms are small by power counting. Thus, we get the result that there is no contribution to
$W^{\rm 2I(1)}_{\mu_\perp\nu_\perp}(q)$ with our accuracy.

If both $\mu$ and $\nu$ are longitudinal,
using formula

\begin{eqnarray}
\sigma_{\mu_\parallel\alpha}\sigma_{\star j}\otimes\sigma_{\nu_\parallel}^{~\alpha}\sigma_{\bu j}+\mu\leftrightarrow\nu
~=~-2s g^\parallel_{\mu\nu}\big(\sigma_{\star j}\otimes\sigma_{\bu j}-{1\over s}\sigma_{\bu k}\sigma_{\star j}\otimes \sigma_\star^{~k}\sigma_{\bu j}\big)
\label{WI26}
\end{eqnarray}

we get

\begin{eqnarray}
\cheW^{\rm 2I(1)}_{\mu\nu}(x)~=~
-{g^\parallel_{\mu\nu}\over 2s^2}\big\{\langle\bsi A^i(x)\sigma_{\star j}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^j(0)\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\label{WI27}\\
\nonumber
\end{eqnarray}

where we omitted contribution from the second term in the r.h.s. of Eq. (<ref>) due
to power counting coming from Eq,  (<ref>). Now, using Eq. (<ref>) one obtains

\begin{eqnarray}
\cheW^{\rm 2I(1)}_{\mu\nu}(x)~=~
-{g^\parallel_{\mu\nu}\over 2s^2}\big\{\langle\bsi(x) [A_i(x)\sigma_{\star j}-A_j(x)\sigma_{\star i}]{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^j(0)\sigma_{\bu i}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle\bsi A^i(x)\sigma_{\star i}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^j(0)\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B  \big\}
\label{WI28}
\end{eqnarray}

due to Eqs. (<ref>) and  (<ref>).  Thus, we get

\begin{eqnarray}
\cheW^{\rm 2I(1)}_{\mu\nu}(x)~=~0
\label{WI29}
\end{eqnarray}

Similarly, one can demonstrate that  the contribution to $\cheW^{\rm I}_{\mu\nu}(x)$ of Eq. (<ref>) coming from $\Bxi_{A}$ and $\Bxi_{B}$

Let us now consider term coming from  $\Bxi_2$ and  $\Xi_2$. After separating color-singlet contributions, it takes the form

\begin{eqnarray}
\cheW^{\rm 3I}_{\mu\nu}(x)~=~
-{1\over 4s^3} \big\{\langle \bsi A^j(x)\sigma_{\mu\xi}A^i\psi(0)\rangle_A
\langle\big(\bar\psi_B{1\over \beta}\big)\sigma_{\bu i}\sigma_\nu^{~\xi}\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle \bsi A^j(x)\sigma_{\mu\xi}\gamma_5A^i\psi(0)\rangle_A
\langle\big(\bar\psi_B{1\over \beta}\big)\sigma_{\bu i}\sigma_\nu^{~\xi}\sigma_{\bu j}{1\over\beta}\psi(x)\rangle_B\big\}
\label{WI30}
\end{eqnarray}

The power counting for similar Eq. (<ref>) in Sect. <ref> shows that we need
to take $\mu$ and $\nu$ both longitudinal:

\begin{eqnarray}
\cheW^{\rm 3I}_{\mu\nu}(x)~=~
-{p_{1\mu}p_{1\nu}\over s^5} \big\{\langle \bsi A^j(x)\sigma_{\star k}A^i\psi(0)\rangle_A
\langle\big(\bar\psi_B{1\over \beta}\big)\sigma_{\bu i}\sigma_\star^{~k}\sigma_{\bu j}\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle \bsi A^j(x)\sigma_{\star k}\gamma_5A^i\psi(0)\rangle_A
\langle\big(\bar\psi_B{1\over \beta}\big)\sigma_{\bu i}\sigma_\star^{~k}\sigma_{\bu j}{1\over\beta}\psi(x)\rangle_B\big\}
\nonumber\\
=~{ip_{1\mu}p_{1\nu}\over s^4} \big\{\langle \bsi A^j(x)\sigma_{\star k}A^i\psi(0)\rangle_A
\langle\big(\bar\psi_B{1\over \beta}\big)\gamma_i\sigma_\bu^{~k}\gamma_j\gamma_5{1\over\beta}\psi(x)\rangle_B
\nonumber\\
\nonumber\\
=~{ip_{1\mu}p_{1\nu}\over s^4} \big\{\langle \bsi(x)\notA(x)\sigma_{\star k}\notA(0)\psi(0)\rangle_A
\langle\big(\bar\psi_B{1\over \beta}\big)\sigma_\bu^{~k}\gamma_5{1\over\beta}\psi(x)\rangle_B~=~0
\label{WI31}
\end{eqnarray}

where we used formula $\sigma_{\star k}\otimes\sigma_\bu^{~k}\gamma_5-\sigma_{\star k}\gamma_5\otimes\sigma_\bu^{~k}=0$ following from Eq.  (<ref>). Also, one can demonstrate that the contribution to $\cheW^{\rm I}_{\mu\nu}$ coming from $\Bxi_{A}$ and  $\Xi_{A}$ vanishes.
Finally, similarly to Eq. (<ref>), we can neglect terms $\sim{1\over N_c^2}$
coming from $\Bxi_1,\Xi_{B}$ and  $\Bxi_2,\Xi_{A}$, see the discussion in Sect.
<ref>. Thus, we get the result that the contribution to  $\cheW^{\rm I}_{\mu\nu}$  coming from two quark-quark-gluon TMDs vanishes with
our accuracy:

\begin{equation}
\cheW^{\rm 2I}_{\mu\nu}~=~0
\label{WI32}
\end{equation}

§.§.§ Exchange-type power corrections

The exchange-type power corrections to eq. (<ref>) are

\begin{eqnarray}
\cheW^{\rm I2Sex}_{\mu\nu}(x)
~=~{N_c\over 2s}\sum_{f,f'}\langle A,B|\Big((e_{f'}c_fa_f-e_fc_{f'}a_{f'})
\nonumber\\
\nonumber\\
+\mu\leftrightarrow\nu-x\leftrightarrow 0
\label{wi2sex}
\end{eqnarray}

The terms in the r.h.s.  differ from those in
Eqs. (<ref>),  (<ref>), and  (<ref>) by
replacement of “$+x\leftrightarrow 0$” by
“$-x\leftrightarrow 0$” which leads to change sign of “$\pm$  c.c.” terms in
those equations so we get instead of Eq. (<ref>)

\begin{eqnarray}
W^{\rm I2Sex}_{\mu\nu}(q)
~=~{iN_c\over 2(N_c^2-1)Q_\parallel^2}\sum_{f,f'}\!\int\! d^2k_\perp\Big\{
\nonumber\\
+~g^\perp_{\mu\nu}(k,q-k)_\perp]I^{1ff'}_{--}(q,k_\perp)-g^\perp_{\mu\nu} (k,q-k)_\perp I^{2ff'}_{--}(q,k_\perp)
\label{wexi2s}\\
+~c_fe_{f'}[\epsilon_{\mu m}k^m(q-k)_\nu+\mu\leftrightarrow\nu]J^{1ff'}_{+-}(q,k_\perp)
[k_\mu\epsilon_{\nu n}(q-k)^n+\mu\leftrightarrow\nu]  J^{1ff'}_{-+}(q,k_\perp)\Big\}
\nonumber
\end{eqnarray}

where $J^i_{\pm\pm}$ and $I^i_{--}$ are defined in Eqs. (<ref>) and (<ref>).

§.§.§ Result for symmetric interference term $W^{\rm I2San}$

As usual, we represent the result for hadronic tensor $W^{\rm I2S}_{\mu\nu}(q)$ as a sum
of the “annihilation” and “exchange” parts:

\begin{eqnarray}
W^{\rm I2S}_{\mu\nu}(q)~=~W^{\rm I2San}_{\mu\nu}(q)~+~W^{\rm I2Sex}_{\mu\nu}(q)
\nonumber\\
\end{eqnarray}

where the exchange-type corrections are presented in Eq. (<ref>) above while
$W^{\rm I2San}_{\mu\nu}(q)$ is given by
the sum of  Eqs. (<ref>) and (<ref>)

\begin{eqnarray}
W^{\rm I2San}_{\mu\nu}(q)~=~\sum_{f}e_fc_fW^{\rm If}_{\mu\nu}(q),
\nonumber\\
W^{\rm If}_{\mu\nu}(q)
~=~-{i\epsilon_{\mu j}\over 2m^2}\!\int\! d^2k_\perp
\Big([k^j(q-k)^\perp_\nu+(q-k)_\perp^j k^\perp_\nu]
\nonumber\\
+~{4\over Q_\parallel^2}[\beta p_{2\nu}k^j(q-k)_\perp^2+\alpha p_{1\nu}(q-k)^jk_\perp^2]\Big)\pizh^f(q,k_\perp)
\label{wi2sresult}
\end{eqnarray}

Let us check gauge invariance of annihilation part of interference hadronic tensor

$q^\mu W^{\rm I2San}_{\mu\nu}(q)~=~0$. First, let us rewrite it in as follows:

\begin{eqnarray}
W^{\rm If}_{\mu\nu}(q)
~=~-{i\epsilon_{\mu j}\over 2m^2}\Big(\delta_\nu^i-2{q^\parallel_\nu q^i\over Q_\parallel^2}\Big)
\!\int\! d^2k_\perp
[k^j(q-k)_i+(q-k)^j k_i+\delta_i^j(k,q-k)_\perp]\pizh^f(q,k_\perp)
\nonumber\\
+~{i\epsilon_{\mu j}\over m^2}{\tilq_\nu\over Q_\parallel^2}\!\int\! d^2k_\perp
\label{wi2sresult}
\end{eqnarray}

The integral in the first term is proportional to $(2q_i q^j+\delta_i^q q_\perp^2)$ and it is easy to see

\begin{eqnarray}
q^\mu\Big[\epsilon_{\mu j}\Big(\delta_\nu^i-2{q^\parallel_\nu q^i\over Q_\parallel^2}\Big)
(2q_iq^j+\delta_i^j) +\mu\leftrightarrow\nu\Big]~=~0
\end{eqnarray}

The integral in the second term is proportional to $q^j$ so

\begin{eqnarray}
q^\mu\Big[\epsilon_{\mu j}q^j
\end{eqnarray}

and therefore $q^\mu W^{\rm If}_{\mu\nu}(q)=0$.

§.§ Antisymmetric interference term of tensor $W^{\rm I2}$

Similarly to the symmetric case (<ref>),
from definitions (<ref>) and (<ref>) we get

\begin{eqnarray}
\cheW^{\rm I2A}_{\mu\nu}(x)
~=~{N_c\over 2s}\sum_{f,f'}\Big((e_{f'}c_fa_f-e_fc_{f'}a_{f'})
\langle A,B|[\bsi(x)\gamma_\mu\psi(x)]^f[\bsi(0)\gamma_\nu\psi(0)]^{f'}
\nonumber\\
+~e_fc_{f'}\langle A,B|[\bsi(x)\gamma_\mu\psi(x)]^f[\bsi(0)\gamma_\nu\gamma_5\psi(0)]^{f'}
\nonumber\\
-~e_{f'}c_f\langle A,B|[\bsi(x)\gamma_\mu\gamma_5\psi(x)]^f[\bsi(0)\gamma_\nu\psi(0)]^{f'}|A,B\rangle\Big)
\label{wi2a}
\end{eqnarray}

§.§.§ Annihilation-type power corrections

Again, let us start with annihilation-type power corrections. Since in this case $f=f'$, the first term in the r.h.s.
of Eq. (<ref>) vanishes and the second can be written as

\begin{eqnarray}
\chew^{\rm I2Aa}_{\mu\nu}(x)~=~\sum_{f}e_fc_f\chew^{\rm IAf}_{\mu\nu}(x)
\nonumber\\
\chew^{\rm IAf}_{\mu\nu}(x)
~=~{N_c\over 2s}\langle A,B|[\Bsi_1(x)\gamma_\mu\Bsi_2(x)][\Bsi_2(0)\gamma_\nu\gamma_5\Bsi_1(0)]
\nonumber\\
-\mu\leftrightarrow\nu|A,B\rangle+x\leftrightarrow 0
\label{wiaa}
\end{eqnarray}

where we made the usual replacement $\psi\rightarrow\Psi_1+\Psi_2$.
After Fierz transformation (<ref>) the r.h.s. of the above equation turns to

\begin{eqnarray}
\cheW^{\rm IA}_{\mu\nu}(x)~
\label{wi2aa}\\
i{N_c\over 2s} \langle p_A,p_B|
\nonumber\\
\nonumber
\end{eqnarray}

where we suppressed flavor label.

Let's us demonstrate that $\cheW^{\rm IA}_{\mu\nu}(x)$  is small in our approximation.
After  Fierz transformation (<ref>)
one obtains

\begin{eqnarray}
\cheW^{\rm IA}_{\mu\nu}(x)~
\nonumber\\
&&\hspace{-1mm}=~{iN_c\over 2s}
\langle A,B|-[\Bsi_1^m(x)\gamma_5\Psi_1^n(0)][\Bsi_2^n(0)\sigma_{\mu\nu}\Psi_2^m(x)]
\nonumber\\
|A,B\rangle+x\leftrightarrow 0
\nonumber\\
\label{chew0e}
\end{eqnarray}

It is convenient to convolute $W^{\rm IA}_{\mu\nu}(x)$ with $\epsilon_{\mu\nu\alpha\beta}$ and consider

\begin{equation}
\pizW^{\rm IA}_{\mu\nu}(x)~=~{i\over 2}\epsilon_{\mu\nu}^{~~\alpha\beta}W^{\rm IA}_{\alpha\beta}(x)
\label{pizwei}
\end{equation}

\begin{eqnarray}
\chepizW^{\rm Ia}_{\mu\nu}(x)~
\label{pizwei1}\\
=~{iN_c\over 2s}\langle A,B|-[\Bsi_1^m(x)\gamma_5\Psi_1^n(0)][\Bsi_2^n(0)\sigma_{\mu\nu}\gamma_5\Psi_2^m(x)]
\nonumber\\
|A,B\rangle+x\leftrightarrow 0
\nonumber
\end{eqnarray}

This is similar to $\chew^{\rm as}_{\mu\nu}(x)$ of Eq. (<ref>) studied in previous Section. The only difference is
the relative sign between the first and the second term in the r.h.s. of these equations (and replacement of
“$-x\leftrightarrow 0$” by
“$+x\leftrightarrow 0$” which does not change power counting).
Let us qiuckly check that this relative sign does not change the result that the contribution is small.
First, let us consider the one-gluon contribution similar to Eq. (<ref>)

\begin{eqnarray}
{i\over s^2}\big\{\langle\bsi(x)\notp_2\gamma_i{1\over\alpha}\psi(0)\rangle_A\langle\bsi(0)B^i(0)\sigma_{\mu\nu}\psi(x)\rangle_B
\label{pizwei2}\\
-\psi(0)\otimes \psi(x)\leftrightarrow \gamma_5\psi(0)\otimes \gamma_5\psi(x)\big\}-x\leftrightarrow 0
\nonumber
\end{eqnarray}

As discussed in previous Section after Eq. (<ref>), the two  terms in the r.h.s. of Eq. (<ref>) vanish separately for transverse $\mu$ and $\nu$,
are both small for one longitudinal and one transverse index, and the term which changed sign is neglected in Eq. (<ref>) so in all cases
the relative sign does not matter. Similarly, one can check this for other one-gluon terms.

Let us now consider two-gluon term coming from $\Xi_1$ and $\Xi_2$

\begin{eqnarray}
{i\over 2s^3}\big\{\langle\bsi A^j(x)\sigma_{\star i}{1\over \alpha}\psi(0)\rangle_A\langle\bsi B^i(0)\sigma_{\mu\nu}\sigma_{\bu j}
{1\over\beta}\psi(x)\rangle_B~+~ \langle\bsi A^j(x)\sigma_{\mu\nu}\sigma_{\star i}{1\over \alpha}\psi(0)\rangle_A
\nonumber\\
\times~\langle\bsi B^i(0)\sigma_{\bu j}{1\over\beta}\psi(x)\rangle_B
-\psi(0)\otimes \psi(x)\leftrightarrow \gamma_5\psi(0)\otimes \gamma_5\psi(x)\big\}+x\leftrightarrow 0
\label{pizwei3}
\end{eqnarray}

Again, this differs from Eq. (<ref>) by relative sign between two terms. Looking at the derivation of Eq. (<ref>) we see that
both terms in the r.h.s. vanish separately due to Eq. (<ref>) so

\begin{eqnarray}
{\rm Eq. ~(\ref{pizwei3})}
~=~{1\over 2s^3}\big\{\langle\bsi [A_\mu(x)\sigma_{\star \nu_\perp}-\mu\leftrightarrow\nu]{1\over \alpha}\psi(0)\rangle_A\langle\bsi B^i(0)\sigma_{\bu i}{1\over\beta}\psi(x)\rangle_B
\nonumber\\
+~\langle\bsi  A^i(x)\sigma_{\star i}{1\over \alpha}\psi(0)\rangle_A\langle\bsi [B_\mu(0)\sigma_{\bu \nu_\perp}-\mu\leftrightarrow \nu]{1\over\beta}\psi(x)\rangle_B\big\}
+x\leftrightarrow 0~=~0
\label{pizwei4}
\end{eqnarray}

If $\mu$ and $\nu$ are longitudinal, we get

\begin{eqnarray}
{\rm Eq. ~(\ref{pizwei3})}~=~{2i\over s^5}(p_{1\mu}p_{2\nu}-\mu\leftrightarrow\nu)
\big\{\langle\bsi A^j(x)\sigma_{\star i}{1\over \alpha}\psi(0)\rangle_A\langle\bsi B^i(0)\sigma_{\star\bu}\sigma_{\bu j}
\nonumber\\
+~ \langle\bsi A^j(x)\sigma_{\star\bu}\sigma_{\star i}{1\over \alpha}\psi(0)\rangle_A\langle\bsi B^i(0)\sigma_{\bu j}
\nonumber\\
-~\psi(0)\otimes \psi(x)\leftrightarrow \gamma_5\psi(0)\otimes \gamma_5\psi(x)\big\}-x\leftrightarrow 0
\label{pizwei5}
\end{eqnarray}

If now one of indices is longitudinal and the other transverse, the contribution is small due to power counting as discussed
after Eq. (<ref>) so the  two-gluon term coming from $\Xi_1$ and $\Xi_2$ vanishes with our accuracy.
The corresponding contribution coming from $\Bxi_1$ and $\Bxi_2$
vanishes for the same reason.

Finally, as shown in Eq. (<ref>),  terms coming from $\Bxi_1$, $\Xi_1$ and from $\Bxi_{2B}$, $\Xi_{2B}$ are small due to power counting.
Also, as usually we neglect the terms coming from  $\Bxi_1$, $\Xi_{2B}$ and from $\Bxi_{2B}$, $\Xi_1$ are $\sim{1\over N_c^2}$ so

\begin{equation}
\chepizW^{\rm IA}_{\mu\nu}(q)~=~\cheW^{\rm IA}_{\mu\nu}(q)~=~0
\label{wiazero}
\end{equation}

with our accuracy.

§.§.§ Exchange-type power corrections and the result for $W^{\rm I2}$

The exchange-type power corrections to eq. (<ref>) are

\begin{eqnarray}
(\cheW^{\rm ex})^{\rm I2A}_{\mu\nu}(x)
~=~{N_c\over 2s}\sum_{f,f'}\langle A,B|\Big((e_{f'}c_fa_f-e_fc_{f'}a_{f'})
\nonumber\\
-\mu\leftrightarrow\nu+x\leftrightarrow 0
\label{wi2aex}
\end{eqnarray}

Similarly to the symmetric case, one can use formulas from Sect. <ref>
with change of signs of complex conjugations.
Indeed, the terms in the r.h.s. of Eq. (<ref>) differ from those in
Eqs. (<ref>) and  (<ref>) by
replacement of “$+x\leftrightarrow 0$” by
“$-x\leftrightarrow 0$” which leads to change sign of “$\pm$  c.c.” terms in
those equations.  We get

\begin{eqnarray}
(W^{\rm I2A})^{\rm ex}_{\mu\nu}(q)~=~-{N_c\over (N_c^2-1)Q_\parallel^2}\sum_{f,f'}c_fc_{f'}
\!\int\! d^2k_\perp
\Big\{a_f[\epsilon_{\mu m}k^m(q-k)_\nu-\mu\leftrightarrow\nu]
\nonumber\\
a_{f'}[k_\mu\epsilon_{\nu n}(q-k)^n-\mu\leftrightarrow\nu]  I^{2ff'}_{-+}(q,k_\perp)
\Big\}
\label{wi2aexotvet}
\end{eqnarray}

cf. Eq. (<ref>). As usually, exchange power corrections exist only for transverse $\mu$ and $\nu$.

Finally, since we proved in previous Section that annihilation-type power corrections $W^{\rm I2A}_{\mu\nu}$ vanish, the total result for $W^{\rm I2A}_{\mu\nu}$ is equal to the “exchange” part:

\begin{equation}
W^{\rm I2A}_{\mu\nu}(q)~=~(W^{\rm I2A})^{\rm ex}_{\mu\nu}(q)~=~{\rm r.h.s~of~Eq.~(\ref{wi2aexotvet}) }
\end{equation}

§ RESULTS

In this section we will take into account only gauge-invariant terms coming from
annihilation-type terms proportional to TMDs $f_1$ and $h_1^\perp$.  The reason to neglect
annihilation-type $\sim W^{\rm 2H}(q)$ is that the twist-three matrix elements (<ref>) are virtually unknown, and exchange-type power corrections can presumably be neglected due to extra ${1\over N_c}$.
Anyway, taking into account leading-twist contributions and their “gauge-invariance-restoring”
counterparts appears to a good start for estimations of DY hadronic tensors.

§.§ Hadronic tensors in Collins-Soper frame

In Collins-Soper frame the hadronic tensors are parametrized in terms of $q$ and three unit
vectors $X,Y,Z$ orthogonal to $q$ and to each other. In terms of Sudakov variables they are

\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}q+{Q\over Q_\perp Q_\parallel}q_\perp\Big)
\label{xyz}\\
Y_\mu~=~-{1\over Q}\epsilon_{\mu\nu\lambda\rho}X^\nu Z^\lambda q^\rho
~=~\epsilon_{\mu i}{q_\perp^i\over Q_\perp}
\end{eqnarray}

where $Q_\perp\equiv |q_\perp|$.

§.§.§ Hadronic tensor for photon-mediated DY process

We  parametrize photon-mediated hadronic tensor in a standard way (up to extra $N_c$ and flavor factors), separately for $W^{\rmF f}_{\mu\nu}$ and $W^{\rmH f}_{\mu\nu}$ defined in Eq. (<ref>)

\begin{eqnarray}
W^{\rmF f}_{\mu\nu}(q)
~=~-\big(g_{\mu\nu}-{q_\mu q_\nu\over q^2}\big)(W^{\rmF f}_T+W^{\rmF f}_{\dd})
-2X_\mu X_\nu W^{\rmF f}_\dd
\nonumber\\
+~Z_\mu Z_\nu(W^{\rmF f}_L-W^{\rmF f}_T-W^{\rmF f}_\dd)
-(X_\mu Z_\nu+X_\nu Z_\mu)W^{\rmF f}_\Delta
\label{WFdec}
\end{eqnarray}

and similarly

\begin{eqnarray}
W^{\rmH f}_{\mu\nu}(q)
~=~-\big(g_{\mu\nu}-{q_\mu q_\nu\over q^2}\big)(W^{\rmH f}_T+W^{\rmH f}_{\dd})
-2X_\mu X_\nu W^{\rmH f}_\dd
\nonumber\\
+~Z_\mu Z_\nu(W^{\rmH f}_L-W^{\rmH f}_T-W^{\rmH f}_\dd)
-(X_\mu Z_\nu+X_\nu Z_\mu)W^{\rmH f}_\Delta
\label{WHdec}
\end{eqnarray}

The expressions for $W_i(q)$ can be easily obtained from Eqs. (<ref>) and (<ref>):

\begin{eqnarray}
W^{\rmF f}_T(q)~=~\!\int\! d^2k_\perp
\big[1-{q_\perp^2\over 2Q_\parallel^2}\big]F^f(q,k_\perp),~~~~~
W^{\rmF f}_L(q)~=~\!\int\! d k_\perp {(q-2k)_\perp^2\over Q_\parallel^2}F^f(q,k_\perp)
\label{WFs}\\
W^{\rmF f}_\dd(q)~=~\!\int\!d^2k_\perp {q_\perp^2\over 2Q_\parallel^2}F^f(q,k_\perp),~~~~~
W^{\rmF f}_\Delta(q)~=~{Q\over Q_\parallel^2Q_\perp}
\!\int\!d^2k_\perp (q,q-2k)_\perp F^f(q,k_\perp)
\nonumber
\end{eqnarray}

\begin{eqnarray}
W^{\rmH f}_T(q)
~=~{1\over 2m^2Q^2}\!\int\!d^2k_\perp \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\\
W_L^{\rmH f}(q)~=~
{1\over m^2Q^2}\!\int\! d k_\perp
\big(2k_\perp^2(q-k)_\perp^2-[k_\perp^2+(q-k)_\perp^2](k,q-k)_\perp\big)
\nonumber\\
W_\dd^{\rmH f}(q)~=~{1\over m^2}
\!\int\!d^2k_\perp
\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)H^f(q,k_\perp)
\nonumber\\
W_\Delta^{\rmH f}(q)~=~-{Q\over Q_\parallel^2Q_\perp m^2}
\!\int\!d^2k_\perp (q,q-2k)(k,q-k)_\perp H^f(q,k_\perp)
\end{eqnarray}

These expressions were obtained in Ref. [16] .

§.§.§ Hadronic tensor for $Z$-mediated DY process

The symmetric part of $W^{\rm ZS}_{\mu\nu}$ is given by Eq. (<ref>)

\begin{equation}
\hspace{-0mm}
W^{\rm ZS}_{\mu\nu}~=~e^2\sum_fc_f^2\big[(a_f^2+1)W_{\mu\nu}^{f\rm F}(q)+(a_f^2-1)W_{\mu\nu}^{f\rm H}(q)\big]
\label{wzsymae}
\end{equation}

where $W_{\mu\nu}^{\rmF f}(q)$ and  $W_{\mu\nu}^{\rm H f}(q)$ are
expressed in terms of $X$ and $Z$ vectors in Eqs. (<ref>) and  (<ref>) above.

The antisymmetric part (<ref>) can be parametrized as

\begin{eqnarray}
W^{\rm ZA}_{\mu\nu}(q)~=~-2i\epsilon_{\mu\nu\lambda\rho}q^\lambda\sum_f a_fc_f^2
\bigg[{Z^\rho\over Q_\parallel}\pizW_4^{\rmF f}(q)
+X^\rho{Q_\perp\over Q_\parallel Q}\pizW_3^{\rmF f}(q)\bigg]
\label{wasym}
\end{eqnarray}

\begin{equation}
\hspace{0mm}
\pizW_4^{\rmF f}(q)~=\int\! d^2k_\perp\pizf
\pizW_3^{\rmF f}(q)~=\int\! d^2k_\perp\Big(1-{2(q,k)_\perp\over q_\perp^2}\Big)\pizf^f(q,k_\perp)
\label{pizwan}
\end{equation}

with $\pizf^f(q,k_\perp)$ given by Eq. (<ref>).

§.§.§ Interference tensors

The symmetric part of the interference tensor $W^{\rm I1}_{\mu\nu}(q)$ is given by Eq. (<ref>)

\begin{eqnarray}
W^{\rm I1S}_{\mu\nu}(q)~
=~\sum_f e_fc_f a_f
\big[W^{\rm F f}_{\mu\nu}(q)+W^{\rm H f}_{\mu\nu}(q)\big]
\label{wi1sym}
\end{eqnarray}

where $W_{\mu\nu}^{\rmF f}(q)$ and  $W_{\mu\nu}^{\rm H f}(q)$ are
given by Eqs. (<ref>) and  (<ref>) above.

The antisymmetric part of $W^{\rm I1}_{\mu\nu}(q)$ is given by Eq. (<ref>) which can be represented
similarly to Eq. (<ref>)

\begin{eqnarray}
W^{\rm I1A}_{\mu\nu}(q)~
=~-i\epsilon_{\mu\nu\lambda\rho}q^\lambda\sum_f e_fc_f
\bigg[{Z^\rho\over Q_\parallel}\pizW_A^{\rmF f}(q)
+X^\rho{Q_\perp\over Q_\parallel Q}\pizW_B^{\rmF f}(q)\bigg]
\label{wi1asym}
\end{eqnarray}

where $\pizW_{A,B}^{\rmF f}(q)$ are given by Eq. (<ref>).

Next, the symmetric part of  $W^{\rm I2}_{\mu\nu}(q)$, given by Eq. (<ref>), can be parametrized as

\begin{eqnarray}
W^{\rm If}_{\mu\nu}(q)
~=~i\sum_f e_fc_f\Big({Q\over Q_\parallel}Y_\mu X_\nu W_1^{\rm If}(q)
+Y_\mu Z_\nu{Q_\perp\over Q_\parallel}W_2^{\rm If}(q)\Big)
\label{wi2sresulte}
\end{eqnarray}

\begin{eqnarray}
W_1^{\rm If}(q)~=~{1\over m^2}\!\int\! d^2k_\perp\Big[{2(q,k)_\perp(q,q-k)_\perp\over q_\perp^2}-(k,q-k)_\perp\Big]\pizh^f(q,k_\perp)
\nonumber\\
W_2^{\rm If}(q)~=~{1\over m^2q_\perp^2}\!\int\! d^2k_\perp\big[(q,k)_\perp(q-k)_\perp^2-(q,q-k)_\perp k_\perp^2\big]\pizh^f(q,k_\perp)
\end{eqnarray}

Finally, the antisymmetric part of  $W^{\rm I2}_{\mu\nu}(q)$ vanishes, see Eq. (<ref>).

\begin{eqnarray}
W^{\rm ZA}_{\mu\nu}(q)~=~-2i\epsilon_{\mu\nu\lambda\rho}q^\lambda\sum_f a_fc_f^2
\bigg[{Z^\rho\over Q_\parallel}\pizW_1^{\rmF f}(q)
+X^\rho{Q_\perp\over Q_\parallel Q}\pizW_2^{\rmF f}(q)\bigg]
\label{wasym}
\end{eqnarray}

§.§.§ Angular coefficients

Rewriting the differential cross section (<ref>) in terms of hadronic tensors listed in previous Section,
we obtain

\begin{eqnarray}
d\sigma~=~{d^4q\over q^4}{e^4\over N_cs\pi^2}\!\int\!d^3ld^3l'\delta(q-l-l')~\matW(q,l,l')
~=~{e^4\over 16\pi^2sN_c}{dQ^2\over Q^2}dYd^2q_\perp d\Omega_l~\matW(q,l,l')
\nonumber\\
\label{desigmae}
\end{eqnarray}

\begin{eqnarray}
\matW(q,l,l')~=~{2\over Q^2}L^{\mu\nu}\sum_f e_f^2[W^{\rmF f}_{\mu\nu}(q)+W^{\rmH f}_{\mu\nu}(q)]
\nonumber\\
+~{2c_e^2\over Q^2}\phi_1(Q^2)\Big\{(a_e^2+1)L^{\mu\nu}\sum_f
c_f^2\big[(a_f^2+1)W_{\mu\nu}^{f\rm F}(q)+(a_f^2-1)W_{\mu\nu}^{f\rm H}(q)\big]
\nonumber\\
-~4a_ea_fc_f^2\epsilon^{\mu\nu\alpha\beta}l_\alpha {l'}_\beta\epsilon_{\mu\nu\lambda\rho}q^\lambda
\Big[{Z^\rho\over Q_\parallel}\pizW_1^{\rmF f}(q)
+X^\rho{Q_\perp\over Q_\parallel Q}\pizW_2^{\rmF f}(q)\Big]\Big\}
\nonumber\\
+~{4c_e\over Q^2} \phi_2(Q^2)
\sum_f c_fe_f\Big\{a_e\big[W^{\rm F f}_{\mu\nu}(q)+W^{\rm H f}_{\mu\nu}(q)\big]
\nonumber\\
-~\epsilon^{\mu\nu\alpha\beta}l_\alpha {l'}_\beta\epsilon_{\mu\nu\lambda\rho}q^\lambda
\Big[{Z^\rho\over Q_\parallel}\pizW_1^{\rmF f}(q)
+X^\rho{Q_\perp\over Q_\parallel Q}\pizW_2^{\rmF f}(q)\Big]\Big\}
\nonumber\\
-~{4c_ea_e\over Q^2}\phi_3(Q^2)
\sum_f c_fe_f\Big[{Q\over Q_\parallel}Y_\mu X_\nu W_1^{\rm If}(q)
+Y_\mu Z_\nu{Q_\perp\over Q_\parallel}W_2^{\rm If}(q)\Big]L^{\mu\nu}
\label{matwe}
\end{eqnarray}

\begin{equation}
\phi_1(Q^2)\equiv{Q^4\over  |m_Z^2-Q^2|^2+\Gamma_Z^2m_Z^2} ,~~~~
\phi_2(Q^2)={(Q^2-m_z^2)\over Q^2} f_1(Q^2),~~~~\phi_3(Q^2)={\Gamma_z m_z\over Q^2} f_1(Q^2)
\end{equation}

The angular dependence in the CS frame can be displayed using
the convolutions of leptonic and hadronic tensors are presented in Sect. <ref>.
One obtains

\begin{eqnarray}
\matW(q,l,l')~=~\sum_f\Big\{\big[e_f^2+c_e^2(a_e^2+1)c_f^2(a_f^2+1)\phi_1
\nonumber\\
\times~\big[\big(W_T^{\rm Ff}+{W_L^{\rm Ff}\over 2}\big)(1+\cos^2\theta)
+{W_L^{\rm Ff}\over 2}(1-3\cos^2\theta)
+W_\Delta^{\rm Ff}\sin2\theta\cos\phi
+W_\dd^{\rm Ff}\sin^2\theta\cos2\phi \big]
\nonumber\\
\nonumber\\
\times~\big[\big(W_T^{\rm Hf}+{W_L^{\rm Hf}\over 2}\big)(1+\cos^2\theta)
+{W_L^{\rm Hf}\over 2}(1-3\cos^2\theta)
+W_\Delta^{\rm Hf}\sin2\theta\sin\phi
+W_\dd^{\rm Hf}\sin^2\theta\cos2\phi \big]
\nonumber\\
\Big[\pizW_1^{\rmF f}\cos\theta+{Q_\perp\over Q}\pizW_2^{\rmF f}\sin\theta\cos\phi\Big]
\nonumber\\
+~a_ea_fc_ec_f\phi_3\Big[W_1^{\rm If}\sin^2\theta\sin 2\phi
+{Q_\perp\over Q}W_2^{\rm If}\sin 2\theta\sin\phi\Big]
\Big\}
\label{wotvet}
\end{eqnarray}

where $\phi_i=\phi_i(Q^2)$ and $W^i=W^i(\alpha_q,\beta_q,Q^2)$.

The above formula (<ref>) is the main result of the paper. It should be compared to
standard representation of angular distribution of DY leptons (<ref>)

§.§ Comparison with LHC measurements

The LHC measurements are integrated over the region of invariant mass of DY pair between 80 and 100 GeV.
In this kinematic region the most important contribution comes from terms in Eq.  (<ref>) multiplied
by $\phi_1$. Indeed, for an estimate we can consider $W^i(\alpha_q,\beta_q,Q^2)\simeq W^i(\alpha_q,\beta_q,m_z^2)$ and get

\begin{equation}
\hspace{-1mm}
\int_{80}^{100}{dQ^2\over Q^2}\simeq 0.2, ~~
\int_{80}^{100}{dQ^2\over Q^2}\phi_1(Q^2)\simeq 95, ~~
\int_{80}^{100}{dQ^2\over Q^2}\phi_2(Q^2)\simeq 0,~~
\int_{80}^{100}{dQ^2\over Q^2}\phi_3(Q^2)\simeq 2.6
\end{equation}

Since the accuracy of our small-x approximation is $\alpha_q=x_A,~\beta_q=x_B\sim 0.1$ we can neglect
all contributions coming form terms not multiplied by $\phi_1$.
Moreover,  both theoretical [28] and phenomenological [29, 30] analysis  indicate  that $h_1^\perp$ is is of order of few percent of $f_1$
and hence in numerical estimates we will disregard the contribution of $h_1^\perp$.
Introducing notations

\begin{eqnarray}
\pizW^{\rmF f}(q)~=~\!\int\! d^2k_\perp F^f(q,k_\perp),~~~~~
\pizW^{\rmF f}_L(q)~=~\!\int\! d k_\perp {(q-2k)_\perp^2\over q_\perp^2}F^f(q,k_\perp)
\label{WFs}\\
\pizW^{\rmF f}_1(q)~=~
\!\int\!d^2k_\perp {(q,q-2k)_\perp\over q_\perp^2} F^f(q,k_\perp)
\nonumber
\end{eqnarray}

we get

\begin{eqnarray}
\matW(q,l,l')~=~c_e^2c_f^2\phi_1(Q^2)\sum_f\Big\{(a_e^2+1)(a_f^2+1)\Big(
\big[\pizw^{\rm Ff}-{Q_\perp^2\over 2Q^2}(\pizw^{\rm Ff}-\pizw_L^{\rm Ff})\big]
\nonumber\\
+~{Q_\perp^2\over 2Q^2}\pizW_L^{\rm Ff}(1-3\cos^2\theta)
+{Q_\perp\over Q}\pizW_1^{\rm Ff}\sin2\theta\cos\phi
+{Q_\perp^2\over 2Q^2}\pizW^{\rm Ff}\sin^2\theta\cos2\phi \big]\Big)
\nonumber\\
+~8a_ea_f\Big[{Q_\perp\over Q}\pizW_3^{\rmF f}\sin\theta\cos\phi+\pizW_4^{\rmF f}\cos\theta\Big]
\Big\}
\label{wotvete}
\end{eqnarray}

where  $\pizW^{\rmF f}_3$ and $\pizW^{\rmF f}_4$ are defined in  Eq. (<ref>).
Since we neglected exchange-type power corrections we should expect the accuracy of order of ${1\over N_c}\sim 30$%.
Let us discuss now some qualitative and semi-quantitative predictions of this equation.
First, let us evaluate $\pizw_i$ at $Q_\perp^2\gg m^2$ following  Ref. [14].

§.§.§ Logarithmical estimates of $\pizW_i$ at $Q_\perp^2\gg m^2$

At $q_\perp^2\gg m^2$ we probe the perturbative tail of TMD $f_1$ which is $\sim{1\over k_\perp^2}$
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},
~~~~~~\barf_1\simeq{\barf(\alpha_q)\over k_\perp^2}
\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},
~~~~~~\barf_1\simeq{\barf(\beta_z)\over k_\perp^2}
\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\\
\pizf^f(q,k_\perp)~\simeq~{\pizf^f(\alpha_q,\beta_q)\over k_\perp^2(q-k)_\perp^2},
\label{festimate}
\end{eqnarray}

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{eqnarray}
\pizW^{\rmF f}(q)~\simeq~{2\pi\over Q_\perp^2}\ln{q_\perp^2\over m^2} F^f(\alpha_q,\beta_q),~~~~~
\pizW_L^{\rmF f}(q)~\simeq~{2\pi\over Q_\perp^2}\Big[\ln{Q_\perp^2\over m^2}+2\ln{Q^2\over Q_\perp^2}\Big] F^f(\alpha_q,\beta_q)
\nonumber\\
\pizW_4^{\rmF f}(q)~\simeq~{2\pi\over Q_\perp^2}\ln{q_\perp^2\over m^2} \pizf^f(\alpha_q,\beta_q),~~~~~~\pizW^{\rmF f}_1(q)~\simeq~\pizw_3^{\rmF f}(q)~\simeq~0.
\label{WFests}
\end{eqnarray}

For numerical estimates of $A_0$ and $A_2$ we will need the ratio

\begin{equation}
r(Q_\perp^2)~\equiv~{\pizw_L^{\rm Ff}(q)\over\pizw^{\rm Ff}(q)}~\simeq~1+2{\ln Q^2/Q_\perp^2\over\ln Q_\perp^2/m^2}~~~~
\label{pizratio}
\end{equation}

which does not depend on quark flavors. Substituting these formulas to Eq. (<ref>), we get
we get

\begin{eqnarray}
\matW(q,l,l')~=~c_e^2(a_e^2+1)c_f^2\phi_1(Q^2)\sum_f(a_f^2+1)\pizw^{\rm Ff}\Big\{\Big(
\Big[1-{Q_\perp^2\over 2Q^2}(1-r)\Big]
\nonumber\\
+~{Q_\perp^2\over 2Q^2}r(1-3\cos^2\theta)
+{Q_\perp^2\over 2Q^2}\sin^2\theta\cos2\phi \big]\Big)
+~{8a_ea_f\over (a_e^2+1)(a_f^2+1)}{\pizf^f\over F^f}\cos\theta
\Big\}
\label{westimate}
\end{eqnarray}

§.§.§ Estimates of $A_0$ and $A_2$

Since our estimate (<ref>) does not depend on flavor, from Eq. (<ref>) we get

\begin{equation}
A_0~=~{Q_\perp^2\over m_z^2}{1+2{\ln m_z^2/Q_\perp^2\over\ln Q_\perp^2/m^2}\over 1+{Q_\perp^2\over m_z^2}{\ln m_z^2/Q_\perp^2\over\ln Q_\perp^2/m^2}}
\label{anot}
\end{equation}

where we replaced $Q^2$ by $m_Z^2$ since this function varies slowly between 80 and 100 GeV. This formula is compared to LHC measurements
in Fig. <ref>.

Comparison of prediction (<ref>)  with  lines depicting angular coefficient $A_0$ in bins  of $Q_\perp$ and $Y<1$ from
Ref. [18]. (long bins) and Ref. [17] (short bins).

The accuracy of our estimate is about 20% which is reasonable since we neglected corrections $\sim{1\over N_c}$.

Next, in our logarithmic approximation

\begin{equation}
A_2~=~{Q_\perp^2/m_z^2\over 1+{Q_\perp^2\over m_z^2}{\ln m_z^2/Q_\perp^2\over\ln Q_\perp^2/m^2}}
\label{adva}
\nonumber
\end{equation}

The comparison with LHC data is shown in Fig. <ref>.

Comparison of prediction (<ref>) with  lines depicting angular coefficient $A_2$ in bins  of $Q_\perp$ and $Y<1$ from
Ref. [18]. (long bins) and [17] (short bins).

The accuracy here is about 30%.

§.§.§ Qualitative checks for other angular coefficients

First, note that between 10 and 30 GeV the coefficient $A_1$ from Table 15 of Ref. [17] is an order of magnitude smaller than $A_0$ from table 14
(or $A_2$ from Table 16) in accordance with our estimate $A_1\simeq 0$. This is another argument in favor of factorization hypothesis  $f_1(x_B,k_\perp)\simeq f_1(x_B)g(k_\perp)$
which is frequently used in current TMD literature.

Second, from Eq. (<ref>) we see that our estimate for $A_4$ involves non-trivial flavor structure, but one thing immediately seen
is that $A_4$ does not depend on $Q_\perp$ in our region 10-30 GeV. As one can see from Table 18 of Ref. [17] , the experimental numbers for  $A_4$
almost do not change in this kinematical region.  Moreover, from Eq. (<ref>)  $\pizf^f(\alpha_q,\beta_q)=0$ at $Y=0$ (i.e., when $\alpha_q=\beta_q$) so we should expect that
$A_4$ in rapidity bin $2>Y>1$ is greater than in bin $Y<1$, and Table 18 confirms that. On the contrary,  the results for $A_0$ and $A_2$ do not change much between those
rapidity beans in accordance with our prediction that $F^f(\alpha_q,\beta_q)$ does not change radically.

Third, our estimate $A_3=0$ is in accordance to the fact that experimental numbers for $A_3$ from Table 17 of Ref. [17] are order
of magnitude smaller than numbers $A_4$ from Table 18.

Finally, the coefficients $A_5$ to $A_7$ do not appear in our estimate (<ref>), again in accordance with Tables 19-21 of Ref. [17]
where the numbers for $A_5,A_6$, and $A_7$ are actually two orders of magnitude smaller than the numbers for $A_0$ or $A_2$ in our kinematical region.

Summarizing, it looks like our  $f_1$-based estimates  point in the right direction, but of course more phenomenological work
is required. Note also that in Ref. [16]  it was demonstrated that TMD evolution  in the double-log approximation does not affect the predictions for angular coefficients,
but the single-log corrections to TMD evolution may change the results for asymmetries. The study is in progress.

§ CONCLUSIONS AND OUTLOOK

To my knowledge,
this is the first analysis of DY angular coefficients in the framework of TMD factorization, and, as we saw, the first signs are encouraging.
There are two questions about the TMD analysis of DY process: can we predict/explain angular dependence of DY cross
section from TMD factorization at LHC kinematical range and reciprocally, can we learn something about proton structure from this?

The answer to first question is probably yes.
We see that when we have sufficient knowledge about TMDs responsible for a particular angular coefficient, that
coefficient comes out with reasonable accuracy even with back-of-the-envelope estimates.
For example, the very naive estimate of the coefficients $A_0$ and $A_2$
agrees with experiment at 30% accuracy. The error  may be due to our approximation of TMDs by perturbative tails or
maybe it is due to ${1\over N_c}$ corrections proportional to higher-twist operators. I hope that careful analysis
involving established models of $f_1(x,q_\perp)$ will distinguish between these two possibilities.
As to the rest of angular coefficients, more details about quark TMDs and at least some guesses about quark-quark-gluon
TMDs are necessary to make quantitative predictions. Also, one needs to take into account fiducial power corrections, see
recent paper [19] for a review.

This bring us to the second question, namely how much we can learn about proton structure from $Z$-boson experiments.
First, as was demonstrated
above, the coefficients $A_0$ and $A_2$ are determined by  TMD $f_1$ with ${1\over N_c}\sim 30$% accuracy.
We used the naive logarithmic estimates  and it is natural to assume  that  the realistic models for $f_1$ will move the curves in Figs.
<ref> and <ref> closer to the experimental points. Also, it should be mentioned, tat the  “factorization hypothesis”
for the LT TMDs like $f(x_B,k_\perp)~\simeq~f(x_B)g(k_\perp)$ seems to be confirmed by experimental data in  a sense that the
angular coefficients which vanish in this approximation are smaller than the non-vanishing ones.
The coefficients $A_1,A_3$, and $A_4$ are determined by non-factorized  part of $f_1$ or higher-twist exchange contributions.

Next, the coefficients $A_5$-$A_7$ seem to be determined by  exchange-type corrections proportional
to higher-twist quark-quark-gluon TMDs. Unfortunately, comprehensive analysis of exchange-type power
corrections requires calculation of higher-twist contributions restoring EM gauge invariance of these corrections,
see the discussion in Sect. <ref>.

Let us also discuss the perturbative corrections to asymmetries. As demonstrated in Ref.  [16]
our estimates of asymmetries  are not affected by summation of Sudakov double logs, but single logs
may bring some changes to tree-level results. It should be also emphasized that,
as discussed in Refs. [15, 14], from the rapidity factorization (<ref>) we get
TMDs with rapidity-only cutoff $|\alpha|<\sigma_t$ or $|\beta|<\sigma_p$. 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  [31, 32, 33]
and three [34] loops.
This difference in cutoffs does not matter for the tree-level formulas of Sect. <ref>, but if one goes beyond the tree level,
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.

The author is grateful to  A. Prokudin  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 transformations

First, let us write down Fierz transformations for symmetric and antisymmetric combinations

\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}

\begin{eqnarray}
\half[(\bsi\gamma_\mu\chi)(\bhi\gamma_\nu\psi)-\mu\leftrightarrow\nu]
{i\over 4}\epsilon_{\mu\nu\alpha\beta}[(\bsi\gamma_\alpha\gamma_5\psi)(\bhi\gamma_\beta\chi)
\nonumber\\
-~{i\over 4}(\bsi\psi)(\bhi\sigma_{\mu\nu}\chi)+{i\over 4}(\bsi\gamma_5\psi)(\bhi\sigma_{\mu\nu}\gamma_5\chi)
+{i\over 4}(\bsi\sigma_{\mu\nu}\psi)(\bhi\chi)-{i\over 4}(\bsi\sigma_{\mu\nu}\gamma_5\psi)(\bhi\gamma_5\chi)
\label{fierzasy}
\end{eqnarray}

Using Eq. (<ref>) one can obtain the following formula

\begin{eqnarray}
{1\over 4}[(\bsi\gamma_\mu\chi)(\bhi\gamma_\nu\gamma_5\psi)+(\bsi\gamma_\mu\gamma_5\chi)(\bhi\gamma_\nu\psi)+\mu\leftrightarrow\nu]
\nonumber\\
-{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\gamma_5\psi)(\bhi\gamma_\beta\chi)
\label{fierz5sym}
\end{eqnarray}

We need also

\begin{eqnarray}
{1\over 4}\big[(\bsi\gamma_\mu\chi)(\bhi\gamma_\nu\gamma_5\psi)-(\bsi\gamma_\mu\gamma_5\chi)(\bhi\gamma_\nu\psi)+\mu\leftrightarrow\nu\big]
~=~{g_{\mu\nu}\over 4}(\bsi\psi)(\bhi\gamma_5\chi)
\label{fierz4}\\
-~{g_{\mu\nu}\over 4}(\bsi\gamma_5\psi)(\bhi\chi)
-{1\over 8}
\big[(\bsi\sigma_{\mu\xi}\psi)(\bhi\sigma_\nu^{~\xi}\gamma_5\chi)-(\bsi\sigma_{\mu\xi}\gamma_5\psi)(\bhi\sigma_\nu^{~\xi}\chi)
\nonumber
\end{eqnarray}

for symmetric tensors and

\begin{eqnarray}
{1\over 4}\big[(\bsi\gamma_\mu\chi)(\bhi\gamma_\nu\gamma_5\psi)+(\bsi\gamma_\mu\gamma_5\chi)(\bhi\gamma_\nu\psi)-\mu\leftrightarrow\nu\big]
\label{fierz5}\\
{i\over 4}\epsilon_{\mu\nu\alpha\beta}\big[(\bsi\gamma_\alpha\psi)(\bhi\gamma_\beta\chi)
\end{eqnarray}

\begin{eqnarray}
{1\over 4}\big[(\bsi\gamma_\mu\chi)(\bhi\gamma_\nu\gamma_5\psi)-(\bsi\gamma_\mu\gamma_5\chi)(\bhi\gamma_\nu\psi)-\mu\leftrightarrow\nu\big]
\label{fierz6}\\
-{i\over 4}(\bsi\gamma_5\psi)(\bhi\sigma_{\mu\nu}\chi)+{i\over 4}(\bsi\psi)(\bhi\sigma_{\mu\nu}\gamma_5\chi)
+{i\over 4}(\bsi\sigma_{\mu\nu}\gamma_5\psi)(\bhi\chi)-{i\over 4}(\bsi\sigma_{\mu\nu}\psi)(\bhi\gamma_5\chi)
\nonumber
\end{eqnarray}

for the antisymmetric ones.

§.§.§ 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_{\star\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_{\star j}~=~g_{ij}-i\epsilon_{ij}\gamma_5
-i\sigma_{ij}-{2i\over s}g_{ij}\sigma_{\bu\star},~~~
{2\over s}\sigma_{\star i}\sigma_{\bu j}~=~g_{ij}+i\epsilon_{ij}\gamma_5
-i\sigma_{ij}+{2\over s}g_{ij}\sigma_{\star\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_{\star k}=-\sigma_{\star k}\sigma_{i j}~=-ig_{ik}\sigma_{\star j}+ig_{jk}\sigma_{\star 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}
\tigma_{\star j}\otimes\tigma_{\bu k}~=~-{s\over 4}g_{jk}\sigma_{\xi\eta}\otimes\sigma^{\xi\eta}
+{s\over 2}\sigma_{k\xi}\otimes\sigma_j^{~\xi}+g_{jk}\sigma_{\bu\xi}\otimes\sigma_\star^{~\xi}-\sigma_{\bu k}\otimes\sigma_{\star j}
\nonumber\\
=~-{s\over 4}g_{jk}\sigma_{mn}\otimes\sigma^{mn}-g_{jk}\sigma_{\star l}\otimes\sigma_\bu^{~l}
+\sigma_{\star k}\otimes\sigma_{\bu j}
+{s\over 2}\sigma_{k l}\otimes\sigma_j^{~l}
\label{gammas13}
\end{eqnarray}

\begin{eqnarray}
\sigma_{\star i}\otimes\sigma_{\bu j}-\sigma_{\star i}\gamma_5\otimes\sigma_{\bu j}\gamma_5
\label{gammas13a}\\
=~-g_{ij}^\perp\sigma_{\star l}\otimes\sigma_\bu^{~l}+\sigma_{\star j}\otimes\sigma_{\bu i}+\sigma_{\star i}\otimes\sigma_{\bu j}
-{s\over 4}g_{ij}^\perp\sigma_{mn}\otimes\sigma^{mn}+{s\over 2}\sigma_{j l}\otimes\sigma_i^{~l}
\nonumber
\end{eqnarray}

\begin{eqnarray}
\sigma_\star^{~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_\star^{~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 “projectile” field $A_i$ can be represented as

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

and similarly for the “target” field

\begin{equation}
B_i(x_\bu,x_\perp)~=~{2\over s}\!\int_{-\infty}^{x_\star}\!dx'_\star~F^{(B)}_{\star i}(x'_\star,x_\perp)
\label{bi}
\end{equation}

in the $B_\bu=0$ gauge. It is convenient to define

\begin{equation}
\tilde{A}_i(x_\bu,x_\perp)~=~{2\over s}\!\int_{-\infty}^{x_\bu}\!dx'_\bu~\tilde{F}^{(A)}_{\star i}(x'_\bu,x_\perp),
~~~\tilde{B}_i(x_\star,x_\perp)~=~{2\over s}\!\int_{-\infty}^{x_\star}\!dx'_\star~\tilde{F}^{(B)}_{\bu i}(x'_\star,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}

We also used

\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\\
{2\over s}(\hatp_2\gamma^i\hatp_1\otimes \hatp_1 B_i+\hatp_2\gamma^i\hatp_1\gamma_5\otimes \hatp_1\gamma_5 B_i)
= -\gamma^i\otimes \breB_i\hatp_1-\gamma^i\gamma_5\otimes \breB_i\hatp_1\gamma_5
\nonumber\\
\gamma_k\gamma^i\hatp_2\otimes B_i\gamma^k
+ \gamma_k\gamma^i\hatp_2\gamma_5\otimes B_i\gamma^k\gamma_5
~=~\hatp_2\otimes \gamma^i\breB_i
\nonumber\\
\slashed{p}_2 \gamma^j\gamma^i\otimes\slashed{p}_1 B_j
~+~ \slashed{p}_2\gamma_5 \gamma^i\gamma^j\otimes\slashed{p}_1\gamma_5 B_j
=\slashed{p}_2\otimes \slashed{p}_1\breB_i +
\slashed{p}_2\gamma_5 \otimes \slashed{p}_1\gamma_5\breB^i,
\nonumber\\
\slashed{p}_2 \gamma^i\gamma^j\otimes\slashed{p}_1 B_j
~+~ \slashed{p}_2\gamma_5 \gamma^i\gamma^j\otimes\slashed{p}_1\gamma_5 B_j
=\slashed{p}_2\otimes \slashed{p}_1\breB^i + \slashed{p}_2\gamma_5\otimes \slashed{p}_1\gamma_5\breB^i.
\label{gammas17}
\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}

In addition, it is convenient to define

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

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 the formula

\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
\label{fla1010}
\end{eqnarray}

\begin{eqnarray}
{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\\
\label{gammas10}
\end{eqnarray}

\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

\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}

We used also

\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_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
\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,
\label{formulas67}
\end{eqnarray}

\begin{eqnarray}
\gamma^j\slashed{p}_2\gamma_m A^k \otimes \gamma^k\slashed{p}_1\gamma_n B^j-m\leftrightarrow n
~=~-\brA_m\slashed{p}_2\otimes\breB_n\slashed{p}_1-m\leftrightarrow n
\nonumber\\
\gamma_m\slashed{p}_2\gamma^j A^k \otimes \gamma^k\slashed{p}_1\gamma_n B^j-m\leftrightarrow n
~=~-\slashed{p}_2\brA_m\otimes\breB_n\slashed{p}_1-m\leftrightarrow n
\nonumber\\
\gamma^j\slashed{p}_2\gamma_mA^k \otimes \gamma_n\slashed{p}_1\gamma^k B^j-m\leftrightarrow n
~=~-\brA_m\slashed{p}_2\otimes\slashed{p}_1\breB_n-m\leftrightarrow n
\nonumber\\
\gamma_m\slashed{p}_2\gamma^j A^k \otimes \gamma_n\slashed{p}_1\gamma^k B^j-m\leftrightarrow n
~=~-\slashed{p}_2\brA_m\otimes\slashed{p}_1\breB_n-m\leftrightarrow n
\label{A1.33}
\end{eqnarray}

Next, using formula (<ref>) we get

\begin{eqnarray}
A^i\sigma_{\star j}\otimes B^j\sigma_{\bu i}-A^i\sigma_{\star j}\gamma_5\otimes B^j\sigma_{\bu i}\gamma_5
\label{gammas14}\\
=~(\barA^j\notp_2\gamma_k-\barA^k\notp_2\gamma_j)\otimes \barB^j\notp_1\gamma_k
-\bar\notA\notp_2\otimes \bar\notB\notp_1-{s\over 2}\bar\notA\gamma_i\otimes \bar\notB\gamma^i
-{s\over 4}\barA^i\sigma_{jk}\otimes \barB_i\sigma^{jk}
\nonumber
\end{eqnarray}

We frequently use this formula in matrix elements like
$\langle \bsi A^i(x)\sigma_{\star j}{1\over \alpha}\psi(0)\rangle_A
\langle\bsi B^j(0)\sigma_{\bu i}{1\over\beta}\psi(x)\rangle_B$. It is easy to see that in such matrix elements
the contribution of two last terms in Eq. (<ref>) is $O\big({q_\perp^2\over s}\big)$ in comparison to the first two ones
so for our purposes

\begin{equation}
\hspace{-0mm}
A^i\sigma_{\star j}\otimes B^j\sigma_{\bu i}-A^i\sigma_{\star j}\gamma_5\otimes B^j\sigma_{\bu i}\gamma_5
\simeq~(\barA^j\notp_2\gamma_k-\barA^k\notp_2\gamma_j)\otimes \barB^j\notp_1\gamma_k
-\bar\notA\notp_2\otimes \bar\notB\notp_1
\label{gammas14a}
\end{equation}

We also need

\begin{eqnarray}
{4\over s^2}\big[A_i\notp_1\notp_2\gamma_j\otimes B_j\notp_2\notp_1\gamma_i
+A_i\notp_1\notp_2\gamma_j\gamma_5\otimes B_j\notp_2\notp_1\gamma_i\gamma_5\big]~
\nonumber\\
=~\brA_i\gamma_j\otimes \breB_j\gamma_i+\brA_i\gamma_j\gamma_5\otimes \breB_j\gamma_i\gamma_5
\label{gammas15}
\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}

The last formula which we need is

\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{gammas17}
\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

\begin{eqnarray}
{1\over 16\pi^3}\!\int\!dx_\bu d^2x_\perp~e^{-i\alpha x_\bu+i(k,x)_\perp}
~\langle \hbsi_f(x_\bu,x_\perp)\gamma^\mu\hsi_f(0)\rangle_A
\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 \hbsi_f(x_\bu,x_\perp)\hsi_f(0)\rangle_A
\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 \hbsi_f(0)\gamma^\mu\hsi_f(x_\bu,x_\perp)\rangle_A
\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 \hbsi_f(0)\hsi_f(x_\bu,x_\perp)\rangle_A
\nonumber
\end{eqnarray}

for the antiquark distributions.
In a general gauge for
projectile and target fields these matrix elements read

\begin{eqnarray}
\langle\hsi_{f}(x)\gamma_\mu \hsi_{f}(0)\rangle_A
\nonumber\\
\langle\hsi_{f}(x)\gamma_\mu \hsi_{f}(0)\rangle_B
\label{gaugelinks}
\end{eqnarray}

and similarly for other operators.

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

\begin{eqnarray}
{1\over 16\pi^3}\!\int\!dx_\star d^2x_\perp~e^{-i\beta x_\star+i(k,x)_\perp}
\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_\star d^2x_\perp~e^{-i\beta x_\star+i(k,x)_\perp}
\nonumber
\end{eqnarray}

\begin{eqnarray}
{1\over 16\pi^3}\!\int\!dx_\star d^2x_\perp~e^{-i\beta x_\star+i(k,x)_\perp}
\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_\star d^2x_\perp~e^{-i\beta x_\star+i(k,x)_\perp}
\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 \hbsi_f(x_\bu,x_\perp)\gamma^\mu\gamma_5\hsi_f(0)\rangle_A
~=~-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 \hbsi_f(0)\gamma^\mu\gamma_5\hsi_f(x_\bu,x_\perp)\rangle_A
~=~-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_\star$
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 \hbsi_f(x_\bu,x_\perp)\sigma^{\mu \nu}\hsi_f(0)\rangle_A
\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 \hbsi_f(0)\sigma^{\mu \nu}\hsi_f(x_\bu,x_\perp)\rangle_A
\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_\star$
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, at small $x$ the cutoff of corresponding longitudinal integrals comes from
the rapidity cutoff $\sigma_a$, see footnote <ref> and corresponding discussion in Ref. [16]. Thus, at small $x$ one can safely put $x=0$ and the corresponding
logarithmic contributions would be proportional to powers of $\alpha_s\ln\sigma_a$ (or, in some cases, $\alpha_s\ln^2\sigma_a$, see e.g.
ref. [35]).

§.§ 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}\\
=~-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(x_\bu,x_\perp){1\over \alpha-\ie}\big)\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(0){1\over \alpha-\ie}\big)\Gamma\Phi(x_\bu,x_\perp)\rangle_A
=~-{1\over\alpha_q}\!\int\!dx_\bu ~e^{-i\alpha x_\bu}\langle\bar\psi()\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_\star$).

Next, we will use QCD equation of motion to reduce quark-quark-gluon TMDs to leading-twist TMDs (cf. 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 \bar\psi(x_\bu,x_\perp)\slashed{p}_2\brA_i(x_\bu,x_\perp)\psi(0)\rangle_A
\label{tw3mael1}\\
=~-\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}\langle \bar\psi(x_\bu,x_\perp)\notA(x_\bu,x_\perp)\slashed{p}_2\gamma_i\psi(0)\rangle_A
\nonumber\\
=~\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\nonumber\\
\times~\big[\langle \hbsi(x_\bu,x_\perp)\slashed{k}_\perp\slashed{p}_2\gamma_i\hsi(0)\rangle_A
\langle \hbsi(x_\bu,x_\perp)\!\stackrel{\leftarrow}{\notD}_\perp\!\gamma_j\slashed{p}_2\gamma_i\hsi(0)\rangle_A\big].
\nonumber
\end{eqnarray}

Using QCD equations of motion 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 \hbsi(x_\bu,x_\perp)\slashed{k}_\perp\slashed{p}_2\gamma_i\hsi(0)\rangle_A
+\alpha_q\langle \hbsi(x_\bu,x_\perp)\slashed{p}_1\slashed{p}_2\gamma_i\hsi(0)\rangle_A\big]
\nonumber\\
=~\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\Big[
-k_i\langle \hbsi(x_\bu,x_\perp)\slashed{p}_2\hsi(0)\rangle_A
+~\alpha_q{s\over 2} \langle \hbsi(x_\bu,x_\perp)\gamma_i\hsi(0)\rangle_A
\nonumber\\
-~i\epsilon_{ij}k^j\langle \hbsi(x_\bu,x_\perp)\slashed{p}_2\gamma_5\hsi(0)\rangle_A
+i{s\over 2}\alpha\epsilon_{ij}
\langle \hbsi(x_\bu,x_\perp)\gamma^j\gamma_5\hsi(0)\rangle_A\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
[It can be demonstrated that $f_1(x,k_\perp^2)$, $f_\perp(x,k_\perp^2)$, and $f_\perp(x,k_\perp^2)$ have the same type of (logarithmic) behavior at low $x$.
Indeed, the low-$x$ behavior is determined by interaction with pomeron.  This interaction is specified by so-called impact factor
and it is easy to check that the impact factor for all three TMDs is similar.
, so
it should be neglected in our kinematical region $s\gg Q^2\gg q_\perp^2$  and we get

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

By complex conjugation

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

For the  corresponding antiquark distributions we get

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

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

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

Next, let us consider

\begin{eqnarray}
{g\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \hsi(x_\bu,x_\perp)\slashed{p}_2\notA(x_\bu,x_\perp)\hsi(0)\rangle_A
\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 \hsi(x_\bu,x_\perp)\slashed{k}_\perp\slashed{p}_2\hsi(0)\rangle_A
+i\langle \hsi(x_\bu,x_\perp)\stackrel{\leftarrow}{\notD}_\perp\slashed{p}_2\hsi(0)\rangle_A\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 \hbsi(x_\bu,x_\perp)\slashed{k}_\perp\slashed{p}_2\hsi(0)\rangle_A
+\alpha_q\langle \hbsi(x_\bu,x_\perp)\slashed{p}_1\slashed{p}_2\hsi(0)\rangle_A\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}
{g\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \hbsi^f(x_\bu,x_\perp)\slashed{p}_2\notA(x_\bu,x_\perp)\hsi^f(0)\rangle_A
~=~i{k_\perp^2\over m_N}h_{1f}^\perp(\alpha_q,k_\perp).
\nonumber\\
\label{9.56}
\end{eqnarray}

By complex conjugation we obtain

\begin{eqnarray}
{g\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \hbsi^f(x_\bu,x_\perp)\slashed{p}_2\notA(0)\hsi^f(0)\rangle_A
~=~i{k_\perp^2\over m_N}h_{1f}^\perp(\alpha_q,k_\perp).
\nonumber\\
\label{9.57}
\end{eqnarray}

For corresponding antiquark distributions one gets in a similar way

\begin{eqnarray}
{g\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \hbsi^f(0)\slashed{p}_2\notA(x_\bu,x_\perp)\hsi^f(x_\bu,x_\perp)\rangle_A
i{k_\perp^2\over m_N}\barh_{1f}^\perp(\alpha_q,k_\perp),
\nonumber\\
{g\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \hbsi^f(0)\slashed{p}_2\notA(0)\hsi^f(x_\bu,x_\perp)\rangle_A
i{k_\perp^2\over m_N}\barh_{1f}^\perp(\alpha_q,k_\perp).
\nonumber\\
\label{9.58}
\end{eqnarray}

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

For the Z-boson hadronic tensor we need these operators with extra $\gamma_5$. From formula (<ref>) we see that

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \bar\psi^f(x_\bu,x_\perp)\slashed{p}_2\brA_i(x_\bu,x_\perp)\gamma_5\psi^f(0)\rangle_A=-i\epsilon_{ij}k^jf_1^f(\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 \bar\psi^f(x_\bu,x_\perp)\brA_i(0)\slashed{p}_2\gamma_5\psi^f(0)\rangle_A=~i\epsilon_{ij}k^jf_1^f(\alpha_q,k_\perp)
\label{gamma5scancel}
\end{eqnarray}

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 \hsi(x_\bu,x_\perp)\notA(x_\bu,x_\perp)\slashed{p}_2\notA(0)\hsi(0)\rangle_A
\nonumber\\
=~{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \hsi(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 \hsi(x_\bu,x_\perp)\notA(x_\bu,x_\perp)\sigma_{\star i}\notA(0)\hsi(0)\rangle_A
\nonumber\\
=~{1\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \hsi(x_\bu,x_\perp)\Big(\slashed{k}_\perp
+i\stackrel{\leftarrow}{\notD}\big)\sigma_{\star 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}
{g\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \hsi(0)\notA(0)\slashed{p}_2\notA(x_\bu,x_\perp)\hsi(x_\bu,x_\perp)\rangle_A
\nonumber\\
=~-{k_\perp^2\over 16\pi^3}\barf_1(\alpha_q,k_\perp)~+~O(\alpha_q,\beta_q)
\nonumber\\
{g\over 8\pi^3s}\!\int\! dx_\bu dx_\perp~e^{-i\alpha_qx_\bu+i(k,x)_\perp}
\langle \hsi(x_\bu,x_\perp)\notA(x_\bu,x_\perp)\sigma_{\star i}\notA(0)\hsi(0)\rangle_A
\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

\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 \hbsi_f(x_\bu,x_\perp)\big[A_i(x)\sigma_{\star j}-\half g_{ij}A^k\sigma_{\star k}(x)\big]\hsi_f(0)\rangle_A~
\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 \hbsi_f(0)[A_i(0)\sigma_{\bu j}-\half g_{ij}A^k\sigma_{\bu k}(0)]\hsi_f(x_\bu,x_\perp)\rangle_A~
\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 element.

§.§ Matrix elements for exchange-type power corrections

We  parametrize “exchange” TMDs  as follows

\begin{eqnarray}
{1\over 4\pi^3s^2}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^{x_\bu}\! dx'_\bu\langle\bsi (x_\bu,x_\perp)\breF_{\star i}(0)\slashed{p}_2\psi(x'_\bu,x_\perp)\rangle_A
\nonumber\\
{1\over 4\pi^3s^2}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^{x_\bu}\! dx'_\bu~\langle \bsi (x_\bu,x_\perp)\slashed{p}_2\breF_{\star i}(0)\psi(x'_\bu,x_\perp)\rangle_A
\nonumber\\
{1\over 4\pi^3s^2}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^{x_\bu}\! dx'_\bu~\langle\bsi (x'_\bu,x_\perp)\slashed{p}_2\breF_{\star i}(0)\psi(x_\bu,x_\perp)\rangle_A
\nonumber\\
{1\over 4\pi^3s^2}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^{x_\bu}\! dx'_\bu~\langle\bsi (x'_\bu,x_\perp)\breF_{\star i}(0)\slashed{p}_2\psi(x_\bu,x_\perp)\rangle_A
\nonumber\\
\label{paramj}
\end{eqnarray}

where $\breF_{\mu\nu}~\equiv~F_{\mu\nu}-i\gamma_5\tilde{F}_{\mu\nu}$.
By complex conjugation we get

\begin{eqnarray}
{1\over 4\pi^3s^2}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^0\! dx'_\bu~\langle\bar\psi(x'_\bu,0_\perp)\slashed{p}_2\breF_{\star i}(x_\bu,x_\perp)\psi(0)\rangle_A
~=~k_i j_1^\ast(\alpha,k_\perp),
\nonumber\\
{1\over 4\pi^3s^2}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^0\! dx'_\bu~\langle\bar\psi(x'_\bu,0_\perp)\breF_{\star i}(x_\bu,x_\perp)\slashed{p}_2\psi(0)\rangle_A
\nonumber\\
{1\over 4\pi^3s^2}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^0\! dx'_\bu\langle\bar\psi(0)\breF_{\star i}(x_\bu,x_\perp)\slashed{p}_2\psi(x'_\bu,0_\perp)\rangle_A
~=~k_i \barj_1^\ast(\alpha,k_\perp),
\nonumber\\
{1\over 4\pi^3s^2}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\!\int_{-\infty}^0\! dx'_\bu\langle\bar\psi(0)\slashed{p}_2\breF_{\star i}(x_\bu,x_\perp))\psi(x'_\bu,0_\perp)\rangle_A
\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_\star$, and $\breF_{\star i}\leftrightarrow \breF_{\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_{\star 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_{\star i}(0)
\nonumber
\end{eqnarray*}

Let us now consider the  corresponding matrix elements from Sect. <ref>. As shown in Ref. [16] ,
Eqs. (<ref>) and (<ref>) lead to (here $x=(x_\bu,x_\perp)$)

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\bsi(x)
\bigg[\begin{array}{c}
\brA_i(0)\slashed{p}_2\\
\slashed{p}_2\brA_i(0)
\end{array}\bigg]
\Big[
\begin{array}{c}j_1(\alpha,k_\perp)\\
\end{array}\Big],
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\big(\bsi {1\over\alpha}\big)(x_\bu,x_\perp)
\bigg[
\begin{array}{c}
\slashed{p}_2\brA_i(0)\\
\brA_i(0)\slashed{p}_2
\end{array}\bigg]
\psi(x_\bu,x_\perp)\rangle_A~
\begin{array}{c}
\barj_1(\alpha,k_\perp)\\
\barj_2(\alpha,k_\perp)
\end{array}\Big]
\nonumber\\
{1\over 8\pi^3s}\!\int\! dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\bsi(0)
\bigg[\begin{array}{c}
\brA_i(x)\slashed{p}_2\\
\slashed{p}_2\brA_i(x)
\end{array}\bigg]
\begin{array}{c}
\barj_1^\ast(\alpha,k_\perp)\\
\barj_2^\ast(\alpha,k_\perp)
\end{array}\Big],
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\big(\bsi {1\over\alpha}\big)(0)
\bigg[\begin{array}{c}
\slashed{p}_2\brA_i(x)\\
\brA_i(x)\slashed{p}_2
\end{array}
\bigg]
\psi(0)\rangle_A
\Big[\begin{array}{c}j_1^\ast(\alpha,k_\perp)\\
\end{array}\Big]
\label{paramjse}
\end{eqnarray}

For the target matrix elements, we obtain ($x\equiv x_\star,x_\perp$)

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\star~e^{-i\beta x_\star+i(k,x)_\perp}\langle\bsi(x)
\bigg[
\begin{array}{c}
\breB_i(0)\slashed{p}_1\\
\slashed{p}_1\breB_i(0)
\end{array}\bigg]
\end{array}\Big],
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\star~e^{-i\alpha x_\star+i(k,x)_\perp}
\langle\big(\bsi {1\over\beta}\big)(x)
\bigg[
\begin{array}{c}
\slashed{p}_1\breB_i(0)\\
\breB_i(0)\slashed{p}_1
\end{array}\bigg]
\psi(x)\rangle_B
\begin{array}{c}
\barj_1(\beta,k_\perp)\\
\barj_2(\beta,k_\perp)
\end{array}\Big],
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\star~e^{-i\beta x_\star+i(k,x)_\perp}
\langle\bsi(0)
\bigg[
\begin{array}{c}
\breB_i(x)\slashed{p}_1\\
\slashed{p}_1\breB_i(x)
\end{array}\bigg]
\begin{array}{c}
\barj_1^\ast(\beta,k_\perp)\\
\barj_2^\ast(\beta,k_\perp)
\end{array}\Big],
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\star~e^{-i\beta x_\star+i(k,x)_\perp}
\langle\big(\bsi {1\over\beta}\big)(0)
\bigg[
\begin{array}{c}
\slashed{p}_1\breB_i(x)\\
\breB_i(x)\slashed{p}_1
\end{array}\Big]
\psi(0)\rangle_B
\end{array}\Big]
\label{paramjset}
\end{eqnarray}

We need also parametrization of matrix elements with extra $\gamma_5$. Since $\brA_i\gamma_5=i\epsilon_{ij}\brA^j$ and
$\breB_i\gamma_5=i\epsilon_{ij}\breB^j$ we get from Eqs. (<ref>), (<ref>)

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\bsi(x)
\bigg[\begin{array}{c}
\brA_i(0)\slashed{p}_2\\
\slashed{p}_2\brA_i(0)
\end{array}\bigg]
\Big[
\begin{array}{c}
\end{array}\Big],
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\big(\bsi {1\over\alpha}\big)(x)
\bigg[
\begin{array}{c}
\slashed{p}_2\brA_i(0)\\
\brA_i(0)\slashed{p}_2
\end{array}\bigg]\gamma_5
\psi(x)\rangle_A~
\begin{array}{c}
\barj_2(\alpha,k_\perp)
\end{array}\Big]
\nonumber\\
{1\over 8\pi^3s}\!\int\! dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\bsi(0)
\bigg[\begin{array}{c}
\brA_i(x)\slashed{p}_2\\
\slashed{p}_2\brA_i(x)
\end{array}\bigg]
\begin{array}{c}
\barj_1^\ast(\alpha,k_\perp)\\
\end{array}\Big],
\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\big(\bsi {1\over\alpha}\big)(0)
\bigg[\begin{array}{c}
\slashed{p}_2\brA_i(x)\\
\brA_i(x)\slashed{p}_2
\end{array}
\bigg]\gamma_5
\psi(0)\rangle_A
\Big[\begin{array}{c}
\end{array}\Big],
\nonumber
\end{eqnarray}

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\bsi(x)
\bigg[\begin{array}{c}
\brA_i(0)\slashed{p}_2\\
\slashed{p}_2\brA_i(0)
\end{array}\bigg]
\Big[
\begin{array}{c}
\end{array}\Big],
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\big(\bsi {1\over\alpha}\big)(x)
\bigg[
\begin{array}{c}
\slashed{p}_2\brA_i(0)\\
\brA_i(0)\slashed{p}_2
\end{array}\bigg]\gamma_5
\psi(x)\rangle_A~
\begin{array}{c}
\barj_2(\alpha,k_\perp)
\end{array}\Big]
\nonumber\\
{1\over 8\pi^3s}\!\int\! dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\bsi(0)
\bigg[\begin{array}{c}
\brA_i(x)\slashed{p}_2\\
\slashed{p}_2\brA_i(x)
\end{array}\bigg]
\begin{array}{c}
\barj_1^\ast(\alpha,k_\perp)\\
\end{array}\Big],
\label{paramjse5}\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\bu~e^{-i\alpha x_\bu+i(k,x)_\perp}
\langle\big(\bsi {1\over\alpha}\big)(0)
\bigg[\begin{array}{c}
\slashed{p}_2\brA_i(x)\\
\brA_i(x)\slashed{p}_2
\end{array}
\bigg]\gamma_5
\psi(0)\rangle_A
\Big[\begin{array}{c}
\end{array}\Big]
\nonumber
\end{eqnarray}

where $x\equiv x_\bu,x_\ast$,

\begin{eqnarray}
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\star~e^{-i\beta x_\star+i(k,x)_\perp}\langle\bsi(x)
\bigg[
\begin{array}{c}
\breB_i(0)\slashed{p}_1\\
\slashed{p}_1\breB_i(0)
\end{array}\bigg]
\end{array}\Big],
\label{paramjset5}\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\star~e^{-i\alpha x_\star+i(k,x)_\perp}
\langle\big(\bsi {1\over\beta}\big)(x)
\bigg[
\begin{array}{c}
\slashed{p}_1\breB_i(0)\\
\breB_i(0)\slashed{p}_1
\end{array}\bigg]\gamma_5
\psi(x)\rangle_B
\begin{array}{c}
\barj_2(\beta,k_\perp)
\end{array}\Big],
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\star~e^{-i\beta x_\star+i(k,x)_\perp}
\langle\bsi(0)
\bigg[
\begin{array}{c}
\breB_i(x)\slashed{p}_1\\
\slashed{p}_1\breB_i(x)
\end{array}\bigg]
\begin{array}{c}
\barj_1^\ast(\beta,k_\perp)\\
\end{array}\Big],
\nonumber\\
{1\over 8\pi^3s}\!\int\! d^2x_\perp dx_\star~e^{-i\beta x_\star+i(k,x)_\perp}
\langle\big(\bsi {1\over\beta}\big)(0)
\bigg[
\begin{array}{c}
\slashed{p}_1\breB_i(x)\\
\breB_i(x)\slashed{p}_1
\end{array}\bigg]\gamma_5
\psi(0)\rangle_B
\end{array}\Big]
\nonumber
\end{eqnarray}

where $x\equiv x_\star,x_\ast$

It is convenient to introduce the following combinations for the parametrization of the exchange-type power corrections

\begin{eqnarray}
(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
J^{1ff'}_{+-}(q,k_\perp)~=~(j_1+\barj_1)^f(\alpha_q,k_\perp)(j^\ast_1-\barj^\ast_1)^{f'}(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
J^{1ff'}_{-+}(q,k_\perp)~=~(j_1-\barj_1)^f(\alpha_q,k_\perp)(j^\ast_1+\barj^\ast_1)^{f'}(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
J^{2ff'}_{+-}(q,k_\perp)~=~(j_2+\barj_2)^f(\alpha_q,k_\perp)(j^\ast_2-\barj^\ast_2)^{f'}(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
J^{2ff'}_{-+}(q,k_\perp)~=~(j_2-\barj_2)^f(\alpha_q,k_\perp)(j^\ast_2+\barj^\ast_2)^{f'}(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\label{Js}
\end{eqnarray}

\begin{eqnarray}
(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
(\beta_q,(q-k)_\perp)~+~{\rm c.c.},
\nonumber\\
(\beta_q,(q-k)_\perp)~+~{\rm c.c.}
\nonumber\\
(\beta_q,(q-k)_\perp)~+~{\rm c.c.}
\label{Is}
\end{eqnarray}

§.§ Products of leptonic tensor and hadronic tensors' structures

The lepton momenta in the Collins-Soper frame are parametrized as

\begin{equation}
\hspace{0mm}
l~=~{Q\over 2}(1,\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta),~~~
l'~=~{Q\over 2}(1,-\sin\theta\cos\phi,-\sin\theta\sin\phi,-\cos\theta)
\end{equation}

The  products of the symmetric leptonic tensor with unit vectors $X$, $Y$, and $Z$ are

\begin{eqnarray}
L^{\mu\nu}\big({q_\mu q_\nu\over q^2}-g_{\mu\nu}\big)~=~2{(q\cdot l)(q\cdot l')\over q^2}+l\cdot l'
\nonumber\\
L^{\mu\nu}X_\mu X_\nu~=~2(l\cdot X)(l'\cdot X)+l\cdot l'~=~{Q^2\over 2}\big[-\sin^2\theta\cos^2\phi+1\big]
\nonumber\\
L^{\mu\nu}Z_\mu Z_\nu~=~2(l\cdot Z)(l'\cdot Z)+l\cdot l'~=~{Q^2\over 2}\big[-\cos^2\theta+1\big]
~=~{Q^2\over 2}\sin^2\theta
\nonumber\\
L^{\mu\nu}X_\mu Z_\nu~=~(l\cdot X)(l'\cdot Z)+(l'\cdot X)(l\cdot Z)
~=~-{Q^2\over 4}\sin2\theta\cos\phi
\nonumber\\
L^{\mu\nu}Y_\mu Z_\nu~=~(l\cdot Y)(l'\cdot Z)+(l'\cdot Y)(l\cdot Z)~=~-{Q^2\over 4}\sin 2\theta\sin\phi
\nonumber\\
L^{\mu\nu}Y_\mu X_\nu~=~(l\cdot Y)(l'\cdot X)+(l'\cdot Y)(l\cdot X)~=~-{Q^2\over 4}\sin^2\theta\sin2\phi
\label{WHdec}
\end{eqnarray}

and therefore

\begin{eqnarray}
{1\over Q^2}L^{\mu\nu}\Big(-\big(g_{\mu\nu}-{q_\mu q_\nu\over q^2}\big)(W_T+W_{\dd})
-2X_\mu X_\nu W_\dd
\label{WFdec}\\
+~Z_\mu Z_\nu(W_L-W_T-W_\dd)
-(X_\mu Z_\nu+X_\nu Z_\mu)W_\Delta\Big)
\nonumber\\
=~\half(W_T+W_L)\bigg[1+\cos^2\theta{W_T-W_L\over W_T+W_L}
+\sin2\theta\sin\phi{W_\Delta\over W_T+W_L}
+\sin^2\theta\cos2\phi {W_\dd\over W_T+W_L}\bigg]
\nonumber
\end{eqnarray}

We also need  products of antisymmetric leptonic tensor with hadron sctructures

\begin{eqnarray}
\epsilon_{\mu\nu\alpha\beta}l^\alpha {l'}^\beta\epsilon^{\mu\nu\lambda\rho}q_\lambda Z_\rho~
\\
-2[(l\cdot q)(l'\cdot Z)-(l'\cdot q)(l\cdot Z)]~=~Q[(l\cdot Z)-(l'\cdot Z)]~=~-Q^2\cos\theta,
\nonumber\\
\epsilon_{\mu\nu\alpha\beta}l^\alpha {l'}^\beta\epsilon^{\mu\nu\lambda\rho}q_\lambda X_\rho~
\nonumber\\
-2[(l\cdot q)(l'\cdot X)-(l'\cdot q)(l\cdot X)]~=~Q[(l\cdot X)-(l'\cdot X)]~=~-Q^2\sin\theta\cos\phi
\nonumber
\end{eqnarray}

and therefore ($Q_\parallel=\sqrt{Q^2+Q_\perp^2}$)

\begin{eqnarray}
\epsilon^{\mu\nu\alpha\beta}l_\alpha {l'}_\beta\epsilon_{\mu\nu\lambda\rho}q^\lambda
\Big[{Z^\rho\over Q_\parallel}\pizW_1^{\rmF f}(q)
+X^\rho{Q_\perp\over Q_\parallel Q}\pizW_2^{\rmF f}(q)\Big]
\nonumber\\
=~-Q^2\Big[\pizW_1^{\rmF f}(q)\cos\theta+{Q_\perp\over Q}\pizW_2^{\rmF f}(q)\sin\theta\cos\phi\Big]\Big[1+O\big({Q_\perp^2\over Q^2}\big)\Big].
\end{eqnarray}

[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]
John C. Collins and Davison E. Soper.
Angular Distribution of Dileptons in High-Energy Hadron Collisions.
Phys. Rev. D, 16:2219, 1977.

[3]
John C. Collins.
Simple Prediction of {QCD} for Angular Distribution of Dileptons
in Hadron Collisions.
Phys. Rev. Lett., 42:291, 1979.

[4]
E. Mirkes and J. Ohnemus.
$W$ and $Z$ polarization effects in hadronic collisions.
Phys. Rev. D, 50:5692–5703, 1994.

[5]
E. Mirkes and J. Ohnemus.
Angular distributions of Drell-Yan lepton pairs at the Tevatron:
Order $\alpha-s^{2}$ corrections and Monte Carlo studies.
Phys. Rev. D, 51:4891–4904, 1995.

[6]
Edmond L. Berger, Jian-Wei Qiu, and Ricardo A. Rodriguez-Pedraza.
Angular distribution of leptons from the decay of massive vector
Phys. Lett. B, 656:74–78, 2007.

[7]
Edmond L. Berger, Jian-Wei Qiu, and Ricardo A. Rodriguez-Pedraza.
Transverse momentum dependence of the angular distribution of the
Drell-Yan process.
Phys. Rev. D, 76:074006, 2007.

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

[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]
I. Balitsky and A. Tarasov.
Power corrections to TMD factorization for Z-boson production.
JHEP, 05:150, 2018.

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

[16]
Ian Balitsky.
Gauge-invariant TMD factorization for Drell-Yan hadronic tensor at
small x.
JHEP, 05:046, 2021.

[17]
Georges Aad et al.
Measurement of the angular coefficients in $Z$-boson events using
electron and muon pairs from data taken at $\sqrt{s}=8$ TeV with the ATLAS
JHEP, 08:159, 2016.

[18]
Vardan Khachatryan et al.
Angular coefficients of Z bosons produced in pp collisions at
$\sqrt{s}$ = 8 TeV and decaying to $\mu^+ \mu^-$ as a function of transverse
momentum and rapidity.
Phys. Lett. B, 750:154–175, 2015.

[19]
Markus A. Ebert, Johannes K. L. Michel, Iain W. Stewart, and Frank J. Tackmann.
Drell-Yan $q_{T}$ resummation of fiducial power corrections at
JHEP, 04:102, 2021.

[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]
R. D. Tangerman and P. J. Mulders.
Intrinsic transverse momentum and the polarized Drell-Yan process.
Phys. Rev., D51:3357–3372, 1995.

[22]
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.

[23]
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.

[24]
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.

[25]
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.

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

[27]
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.

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

[29]
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.

[30]
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.

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

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

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

[34]
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.

[35]
Yuri V. Kovchegov and Matthew D. Sievert.
Calculating TMDs of a Large Nucleus: Quasi-Classical Approximation
and Quantum Evolution.
Nucl. Phys., B903:164–203, 2016.

[14]