# Perturbative $T$-odd asymmetries in the Drell-Yan process revisited

Valery E. Lyubovitskij  0000-0001-7467-572X Institut für Theoretische Physik,
Universität Tübingen,
Kepler Center for Astro and Particle Physics,
Auf der Morgenstelle 14, D-72076 Tübingen, Germany  Departamento de Física y
Centro Científico Tecnológico de Valparaíso-CCTVal,
Universidad Técnica Federico Santa María, Casilla 110-V, Valparaíso, Chile
Millennium Institute for Subatomic Physics at the High-Energy Frontier
(SAPHIR) of ANID,
Fernández Concha 700, Santiago, Chile  Werner Vogelsang  0000-0002-4003-3099
Institut für Theoretische Physik, Universität Tübingen,
Kepler Center for Astro and Particle Physics,
Auf der Morgenstelle 14, D-72076 Tübingen, Germany  Fabian Wunder
0009-0007-4136-7844 Institut für Theoretische Physik, Universität Tübingen,
Kepler Center for Astro and Particle Physics,
Auf der Morgenstelle 14, D-72076 Tübingen, Germany  Alexey S. Zhevlakov
0000-0002-7775-5917 Millennium Institute for Subatomic Physics at the High-
Energy Frontier (SAPHIR) of ANID,
Fernández Concha 700, Santiago, Chile

###### Abstract

We calculate the perturbative $T$-odd contributions to the lepton angular
distribution in the Drell-Yan process. Using collinear factorization, we work
at the first order in QCD perturbation theory where these contributions
appear, ${\cal O}(\alpha_{s}^{2})$, and address both $W^{\pm}$ and
$\gamma/Z^{0}$ boson exchange. A major focus of our calculation is on the
regime where the boson’s transverse momentum $Q_{T}$ is much smaller than its
mass $Q$. We carefully expand our results up to next-to-next-to-leading power
in $Q_{T}/Q$. Our calculation provides a benchmark for studies of $T$-odd
contributions that employ transverse-momentum dependent parton distribution
functions. In the neutral-current case we compare our results for the $T$-odd
structure functions to available ATLAS data.

##  I Introduction

There is a long history of the study of $T$-odd asymmetries in QCD hard-
scattering processes, starting with the seminal papers DeRujula:1971nnp ;
DeRujula:1978bz ; Fabricius:1980wg ; Korner:1980np ; Hagiwara:1981qn ;
Hagiwara:1982cq  that showed how $T$-odd effects may be generated by
absorptive parts of QCD loop diagrams. $T$-odd behavior refers to non-
invariance of observables under so-called naive time-reversal, that is, under
reversal of momenta and spins without interchange of initial and final states.
As shown in the early papers, such behavior can occur even in theories that
are manifestly invariant under true time-reversal. Subsequent work Pire:1983tv
; Hagiwara:1984hi ; Bilal:1990wi ; Carlitz:1992fv ; Mirkes:1992hu ;
Brandenburg:1995nv ; Ahmed:1999ix ; Korner:2000zr ; Yokoya:2007xe ;
Hagiwara:2007sz ; Frederix:2014cba ; Gehrmann-DeRidder:2016cdi ; Benic:2019zvg
; Benic:2021gya ; Benic:2024fvk ; Abele:2022spu  explored $T$-odd QCD
phenomena in a wide range of scattering reactions, among them especially the
Drell-Yan process. In Ref. Hagiwara:1984hi  it was proposed to study $T$-odd
asymmetries appearing in the angular distributions of the charged leptons from
the decay of $W^{\pm}$ bosons produced with high transverse momentum at hadron
colliders. The asymmetries manifest themselves as terms proportional to
$\sin\phi$ or $\sin 2\phi$ in the lepton distribution, where $\phi$ is a
suitably defined azimuthal angle between the lepton plane and the hadron
plane. The $T$-odd part of the spin-averaged differential cross section for
$W^{\pm}$ production was expanded in Hagiwara:1984hi  in terms of three
structure functions which were computed to lowest order of perturbation
theory. The results were later obtained independently in Ref. Mirkes:1992hu
and extended to the case of longitudinal polarization of one of the initial
hadrons in Yokoya:2007xe ; Benic:2024fvk .

In parallel developments, it was realized that $T$-odd effects in QCD may also
arise in hadronic matrix elements, especially in parton distribution functions
(PDFs) Sivers:1989cc ; Boer:1997nt ; Brodsky:2002cx ; Collins:2002kn , where
they are associated with correlations among three-momenta, transverse momenta,
and polarizations of partons and hadrons and again generate azimuthal-angle
dependent terms. This has given rise to an intensive experimental program
aiming at the extraction of such $T$-odd transverse-momentum dependent parton
distributions (TMDs) (for a recent review, see Ref. Boussarie:2023izj ),
either in semi-inclusive lepton scattering (SIDIS) or via the Drell-Yan
process.

The precise connection between $T$-odd effects in perturbative (collinearly-
factorized) hard scattering on the one hand and $T$-odd TMDs on the other has
been an area of active research as well. This issue is important both
theoretically and for phenomenology, where it is central for the “matching” of
resummed calculations based on TMDs to fixed-order perturbation theory. While
much progress has been made for leading-twist observables Ji:2006ub ;
Ji:2006vf ; Gutierrez-Reyes:2017glx ; Gutierrez-Reyes:2018iod ;
Scimemi:2019gge ; Rein:2022odl , it was also realized that for many of the
azimuthal-angle dependent terms in the Drell-Yan and SIDIS cross sections –
both $T$-odd and $T$-even – this matching is non-trivial and will involve TMD
PDFs at next-to-leading power in the hard scale Bacchetta:2008xw ; Boer:2006eq
. Correspondingly, TMD factorization theorems at next-to-leading power were
developed in the literature Ebert:2021jhy ; Gamberg:2022lju ; Rodini:2022wic ;
Rodini:2023plb .

In the present paper, we advance this area of research by specifically
exploring the low-transverse momentum limit of the $T$-odd terms appearing in
the Drell-Yan hard-scattering calculation. The work is carried out in the
spirit of Ref. Abele:2022spu  that addressed $T$-odd effects in SIDIS, but
goes well beyond it in terms of calculational techniques. For our purpose, we
first perform an independent new analytical calculation of the lowest-order
$T$-odd terms in the Drell-Yan cross section, recovering results at this order
from the previous literature Hagiwara:1984hi ; Mirkes:1992hu , but also
extending them by presenting results for pure $Z$-boson exchange and
$\gamma$-$Z$ interference. As one application, we will compare our results to
available ATLAS data ATLAS:2016rnf  for the $T$-odd angular terms taken around
the $Z$ resonance. Our main focus, however, is to carefully expand the results
for low $Q_{T}^{2}/Q^{2}$, where $Q_{T}$ and $Q$ are the boson’s transverse
momentum and mass, respectively. We do this to first and second power in this
ratio, identifying logarithmic behavior as well. As a byproduct we also
uncover a novel simple relation between two of the $T$-odd structure functions
valid at leading power in $Q_{T}^{2}/Q^{2}$ for both partonic channels,
$q\bar{q}$ annihilation and $qg$ Compton scattering.

We hope that our explicit results will be useful in testing TMD factorization
at next-to-leading power and ultimately contribute to a better understanding
of the matching between the TMD and collinearly-factorized regimes. As has
been shown in Refs. Boer:1999mm ; Lu:2011mz ; Liu:2012fha , at leading power
$\gamma/Z$ interference generates a $\sin 2\phi$ azimuthal dependence in the
unpolarized Drell-Yan cross section, entering with the Boer-Mulders function
Boer:1997nt , while a term proportional to $\sin\phi$ is not generated.
Effects beyond leading power have been investigated in Ref. Ebert:2020dfc . We
will not address the potential ramifications of our results for TMDs in this
paper, but rather view our work as a complementary attempt to connect the
physics of small and large $Q_{T}$ and as providing (part of) a “library” for
hard-scattering functions at low transverse momenta. We stress that our
techniques for expanding the cross sections for low $Q_{T}^{2}/Q^{2}$ are
completely general and may also be used in a variety of other settings, such
as for $T$-even contributions, collisions of polarized hadrons, and so forth,
or perhaps even at the next order in perturbation theory, as available from
Gehrmann-DeRidder:2016cdi .

Our paper is organized as follows. In Sec. II we present the definition of the
structure functions parametrizing the lepton-angular distribution for the
Drell-Yan process and the main ingredients for the perturbative calculation of
the $T$-odd contributions. In Sec. III we collect the analytic results, and
subsequently in Sec. IV the small $Q_{T}$ expansion is performed. In Sec. V we
compare our results with the ATLAS data. Section VI concludes our paper. Some
calculational details are collected in the Appendices A–C, and the lengthy
results for the small-$Q_{T}$ expansion to next-next-to-leading-power are
presented in supplemental material to this article.

##  II $T$-odd structure of the Drell-Yan hadronic tensor

The hadronic tensor for the Drell-Yan process can be written in terms of nine
structure functions $W_{i}$. Using a basis of orthogonal unit vectors
$T^{\mu}=q^{\mu}/\sqrt{Q^{2}}=(1,0,0,0)$, $X^{\mu}=(0,1,0,0)$,
$Z^{\mu}=(0,0,0,1)$,
$Y^{\mu}=\epsilon^{\mu\nu\alpha\beta}T_{\nu}Z_{\alpha}X_{\beta}=(0,0,1,0)$,
constructed from the hadron and virtual boson momenta, it reads

$\displaystyle W^{\mu\nu}$ $\displaystyle=$
$\displaystyle(X^{\mu}X^{\nu}+Y^{\mu}Y^{\nu})W_{T}\,+\,i(X^{\mu}Y^{\nu}-Y^{\mu}X^{\nu})W_{T_{P}}\,+\,Z^{\mu}Z^{\nu}W_{L}$
(1) $\displaystyle+$
$\displaystyle(Y^{\mu}Y^{\nu}-X^{\mu}X^{\nu})W_{\Delta\Delta}\,-\,(X^{\mu}Y^{\nu}+Y^{\mu}X^{\nu})W_{\Delta\Delta_{P}}$
$\displaystyle-$
$\displaystyle(X^{\mu}Z^{\nu}+Z^{\mu}X^{\nu})W_{\Delta}\,-\,(Y^{\mu}Z^{\nu}+Z^{\mu}Y^{\nu})W_{\Delta_{P}}$
$\displaystyle+$ $\displaystyle
i(Z^{\mu}X^{\nu}-X^{\mu}Z^{\nu})W_{\nabla}\,+\,i(Y^{\mu}Z^{\nu}-Z^{\mu}Y^{\nu})W_{\nabla_{P}}\,,$

where $q$ is the momentum of the gauge boson $\gamma$, $W^{\pm}$, or $Z^{0}$,
with $q^{2}=Q^{2}$ its Minkowski momentum squared, and
$\epsilon^{\mu\nu\alpha\beta}$ is the four-dimensional Levi-Civita tensor
defined via ${\rm
tr}(\gamma^{5}\gamma^{\mu}\gamma^{\nu}\gamma^{\alpha}\gamma^{\beta})=4\,i\,\epsilon^{\mu\nu\alpha\beta}$,
with $\epsilon^{0123}=-\epsilon_{0123}=-1$.

The number of structure functions, $9=3\times 3$, is determined by the number
of possible helicity settings of the gauge boson in the amplitude and its
complex conjugate. In case of the purely weak Drell-Yan reactions or for
$\gamma$-$Z^{0}$ interference we have nine functions, while in case of the
purely electromagnetic Drell-Yan we have only four $T$-even structure
functions. In general, the Drell-Yan hadronic structure functions may be
classified as:

  * (a)

two transverse functions, the $P$-even $W_{T}$ and the $P$-odd $W_{T_{P}}$,

  * (b)

one longitudinal function $W_{L}$ which is $P$-even,

  * (c)

two transverse-transverse interference (double-spin-flip) functions, the
$P$-even $W_{\Delta\Delta}$ and the $P$-odd $W_{\Delta\Delta_{P}}$,

  * (d)

four transverse-longitudinal interference (single-spin-flip) functions, the
$P$-even $W_{\Delta}$, $W_{\nabla}$ and the $P$-odd $W_{\Delta_{P}}$,
$W_{\nabla_{P}}$.

The lepton angular distribution $dN/d\Omega$ is expanded in terms of the
hadronic structure functions as

$\displaystyle\frac{dN}{d\Omega}$ $\displaystyle=$
$\displaystyle\frac{3}{8\pi(2W_{T}+W_{L})}\,\biggl{[}g_{T}\,W_{T}\,+\,g_{L}\,W_{L}\,+\,g_{\Delta}\,W_{\Delta}\,+\,g_{\Delta\Delta}\,W_{\Delta\Delta}$
(2) $\displaystyle+$ $\displaystyle
g_{T_{P}}\,W_{T_{P}}\,+\,g_{\nabla_{P}}\,W_{\nabla_{P}}\,+\,g_{\nabla}\,W_{\nabla}\,+\,g_{\Delta\Delta_{P}}\,W_{\Delta\Delta_{P}}\,\,+\,g_{\Delta_{P}}\,W_{\Delta_{P}}\biggr{]}\,,$

where $g_{i}=g_{i}(\theta,\phi)$ denote the angular coefficients

$\displaystyle g_{{}_{T}}$ $\displaystyle=$ $\displaystyle
1+\cos^{2}\theta\,,\quad\hskip
17.07182ptg_{{}_{L}}\,=\,1-\cos^{2}\theta\,,\quad\hskip
12.80365ptg_{{}_{T_{P}}}\,=\,\cos\theta\,,$ $\displaystyle
g_{{}_{\Delta\Delta}}$ $\displaystyle=$ $\displaystyle\sin^{2}\theta\,\cos
2\phi\,,\quad\hskip 5.69046ptg_{{}_{\Delta}}\,=\,\sin
2\theta\,\cos\phi\,,\quad\hskip
4.26773ptg_{{}_{\nabla_{P}}}\,=\,\sin\theta\,\cos\phi\,,$ $\displaystyle
g_{{}_{\Delta\Delta_{P}}}$ $\displaystyle=$ $\displaystyle\sin^{2}\theta\,\sin
2\phi\,,\quad\hskip 1.13791ptg_{{}_{\Delta_{P}}}\,=\,\sin
2\theta\,\sin\phi\,,\quad\hskip
11.38092ptg_{{}_{\nabla}}\,=\,\sin\theta\,\sin\phi\,,$ (3)

with $\theta$ and $\phi$ being the polar and azimuthal angles of one of the
decay leptons in the center-of-mass system (c.m.s.) of the lepton pair. The
angle $\phi$ may be taken to define the orientation of the lepton plane with
respect to the hadron plane. In Fig. 1 we show the polar and the azimuthal
angles for the Drell-Yan process in the Collins-Soper frame.

The six angular coefficients $g_{i}$ ($i=T$, $L$, $\Delta\Delta$, $\Delta$,
$\Delta\Delta_{P}$, $\Delta_{P}$) in (II) are invariant under the $P$-parity
transformation $\theta\to\pi-\theta$ and $\phi\to\pi+\phi$, while the other
three coefficients $g_{i}$ ($i=T_{P}$, $\nabla$, $\nabla_{P}$) change their
sign in that case. Therefore, the six partial lepton angular distributions
$dN_{i}/d\Omega$ ($i$=$T$, $L$, $\Delta\Delta$, $\Delta$, $T_{P}$,
$\nabla_{P}$) are also $P$-invariant, whereas the other three distributions
$dN_{i}/d\Omega$ ($i$=$\Delta\Delta_{P}$, $\Delta_{P}$, $\nabla$) are $P$-odd
and also $T$-odd. As can be seen from Eq. (II), the latter distributions are
all proportional to either $\sin\phi$ or $\sin 2\phi$.

Figure 1: Definition of the polar and the azimuthal angles for the Drell-Yan
process in the Collins-Soper frame. The hadron plane is depicted in blue, the
lepton plane in red.

We note in passing that two other commonly employed, and equivalent,
parameterizations of the lepton angular distribution are Lam:1978pu ;
Collins:1977iv ; Mirkes:1992hu ; Boer:2006eq ; Berger:2007jw

$\displaystyle\frac{dN}{d\Omega}$ $\displaystyle=$
$\displaystyle\frac{3}{16\pi}\,\biggl{(}1+\cos^{2}\theta+\frac{A_{0}}{2}(1-3\cos^{2}\theta)+A_{1}\sin
2\theta\cos\phi+\frac{A_{2}}{2}\sin^{2}\theta\cos 2\phi$ (4) $\displaystyle+$
$\displaystyle A_{3}\sin\theta\cos\phi+A_{4}\cos\theta+A_{5}\sin^{2}\theta\sin
2\phi+A_{6}\sin 2\theta\sin\phi+A_{7}\sin\theta\sin\phi\biggr{)}\,,$

and

$\displaystyle\frac{dN}{d\Omega}$ $\displaystyle=$
$\displaystyle\frac{3}{4\pi}\,\frac{1}{\lambda+3}\,\biggl{(}1+\lambda\cos^{2}\theta+\mu\sin
2\theta\cos\phi+\frac{\nu}{2}\sin^{2}\theta\cos 2\phi$ (5) $\displaystyle+$
$\displaystyle\tau\sin\theta\cos\phi+\eta\cos\theta+\xi\sin^{2}\theta\sin
2\phi+\zeta\sin 2\theta\sin\phi+\chi\sin\theta\sin\phi\biggr{)}\,.$

The relations between the three sets of structure functions are recalled in
Appendix A.

The $T$-odd structure functions $W_{\nabla}$, $W_{\Delta\Delta_{P}}$, and
$W_{\Delta_{P}}$ are generated at ${\cal O}(\alpha_{s}^{2})$ in the strong
coupling constant $\alpha_{s}$ by the absorptive parts of parton scattering
amplitudes. The leading contributions arise from the interference of one-loop
and tree-level diagrams. The relevant channels are quark-antiquark
annihilation and quark-gluon Compton scattering. Their one-loop diagrams
providing an absorptive part for photon exchange in Drell-Yan are shown in
Figs. 3 and 3; the diagrams with $Z^{0}$ and $W^{\pm}$ bosons are generated
analogously. The ensuing $T$-odd effects were first studied in Ref.
Hagiwara:1984hi  and later recalculated in Mirkes:1992hu . Here we will
present an independent derivation that will allow us to explore the
low-$Q_{T}$ limit of the results.

Figure 2: One-loop diagrams for $q\bar{q}\to g\gamma$ that produce an
absorptive part.

Figure 3: One-loop diagrams for $qg\to q\gamma$ that produce an absorptive
part.

For our calculation of the $T$-odd structure functions we use a convenient
orthogonal basis of vectors $P,R,K$ Lyubovitskij:2021ges , defined by

$\displaystyle P^{\mu}$ $\displaystyle=$ $\displaystyle(p_{1}+p_{2})^{\mu}\,,$
$\displaystyle R^{\mu}$ $\displaystyle=$ $\displaystyle(p_{1}-p_{2})^{\mu}\,,$
$\displaystyle K^{\mu}$ $\displaystyle=$ $\displaystyle
k_{1}^{\mu}-P^{\mu}\,\frac{P\cdot k_{1}}{P^{2}}-R^{\mu}\,\frac{R\cdot
k_{1}}{R^{2}}=-q^{\mu}+P^{\mu}\,\frac{P\cdot q}{P^{2}}+R^{\mu}\,\frac{R\cdot
q}{R^{2}}\,,$ (6)

which obey the conditions

$\displaystyle P^{2}=-R^{2}=\hat{s}\,,\quad
K^{2}=-\frac{\hat{u}\hat{t}}{\hat{s}}\,,\quad P\cdot R=P\cdot K=R\cdot K=0\,.$
(7)

Here $p_{1}$, $p_{2}$, and $k_{1}$ are the momenta of the two initial partons
and the final-state parton, respectively, satisfying the momentum conservation
relation $p_{1}+p_{2}=k_{1}+q$. Furthermore, $\hat{s}=(p_{1}+p_{2})^{2}$,
$\hat{t}=(p_{1}-q)^{2}$, $\hat{u}=(p_{2}-q)^{2}$, with
$\hat{s}+\hat{t}+\hat{u}=Q^{2}$ the parton-level Mandelstam variables.

The $(P,R,K)$ and $(T,X,Y,Z)$ bases are related by

$\displaystyle X^{\mu}$ $\displaystyle=$
$\displaystyle\frac{T^{\mu}\,\sqrt{1+\rho^{2}}}{\rho}-\frac{P^{\mu}z_{12}^{+}+R^{\mu}z_{12}^{-}}{2Q\rho\sqrt{1+\rho^{2}}}\,=\,\frac{\rho\Big{(}P^{\mu}z_{12}^{+}+R^{\mu}z_{12}^{-}\Big{)}}{2Q\sqrt{1+\rho^{2}}}-\frac{K^{\mu}\,\sqrt{1+\rho^{2}}}{Q\rho}\,,$
$\displaystyle Z^{\mu}$ $\displaystyle=$
$\displaystyle\frac{P^{\mu}z_{12}^{-}+R^{\mu}z_{12}^{+}}{2Q\sqrt{1+\rho^{2}}}\,,$
$\displaystyle Y^{\mu}$ $\displaystyle=$ $\displaystyle-\,\epsilon^{\mu
PRK}\frac{z_{1}z_{2}}{Q^{3}\rho(1+\rho^{2})}\,,$ (8)

where $z_{12}^{\pm}=z_{1}\pm z_{2}$, $Q=\sqrt{Q^{2}}$, and $\epsilon^{\mu
PRK}=\epsilon^{\mu\nu\alpha\beta}\,P_{\nu}\,R_{\alpha}\,K_{\beta}$. We have
$z_{i}=x_{i}/\xi_{i}$, with the momentum fractions $\xi_{i}$ of the partons
defined by $p_{i}=\xi_{i}P_{i}$, and with the momentum fractions of the light-
cone components of the gauge boson,

$\displaystyle x_{1,2}=e^{\pm y}\sqrt{\frac{Q^{2}+Q_{T}^{2}}{s}}$ (9)

at non-zero $Q_{T}$. We also introduce the variables

$\displaystyle x_{1,2}^{0}=e^{\pm y}\,\frac{Q}{\sqrt{s}}$ (10)

relevant in the $Q_{T}=0$ limit.

The hadronic structure functions $W(x_{1},x_{2},\rho^{2})$ for the Drell-Yan
process with colliding hadrons $H_{1}$ and $H_{2}$ are related to the parton-
level structure functions $w^{ab}(x_{1},x_{2},\rho^{2})$ by the QCD collinear
factorization formula

$\displaystyle
W(x_{1},x_{2},\rho^{2})=\frac{1}{x_{1}x_{2}}\,\sum\limits_{a,b}\,\int\limits_{x_{1}}^{1}\,dz_{1}\,\int\limits_{x_{2}}^{1}\,dz_{2}\
w^{ab}(z_{1},z_{2},\rho^{2})\,f_{a/H_{1}}\Big{(}\frac{x_{1}}{z_{1}}\Big{)}\,f_{b/H_{2}}\Big{(}\frac{x_{2}}{z_{2}}\Big{)}\,,$
(11)

where $f_{i/H}(\xi)$ is the parton distribution function (PDF) describing the
$\xi$ distribution of partons of type $i$ in hadron $H$. We are suppressing
here the scale dependence of the PDFs.

We may project onto the parton-level $T$-odd structure functions in the
following way (here we drop parton labels):

$\displaystyle w_{\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(X^{\mu}Y^{\nu}+X^{\nu}Y^{\mu})\,w_{\mu\nu}$ (12)
$\displaystyle=$
$\displaystyle\frac{z_{1}z_{2}}{4Q^{4}\,(1+\rho^{2})^{3/2}}\,\biggl{[}\epsilon^{\mu
PRK}\Big{(}P^{\nu}z_{12}^{+}+R^{\nu}z_{12}^{-}\Big{)}+\epsilon^{\nu
PRK}\Big{(}P^{\mu}z_{12}^{+}+R^{\mu}z_{12}^{-}\Big{)}\biggr{]}\,w_{\mu\nu}\,,$
$\displaystyle w_{\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(Y^{\mu}Z^{\nu}+Y^{\nu}Z^{\mu})\,w_{\mu\nu}$ (13)
$\displaystyle=$
$\displaystyle\frac{z_{1}z_{2}}{4Q^{4}\,\rho\,(1+\rho^{2})^{3/2}}\,\biggl{[}\epsilon^{\mu
PRK}\Big{(}P^{\nu}z_{12}^{-}+R^{\nu}z_{12}^{+}\Big{)}+\epsilon^{\nu
PRK}\Big{(}P^{\mu}z_{12}^{-}+R^{\mu}z_{12}^{+}\Big{)}\biggr{]}\,w_{\mu\nu}\,,$
$\displaystyle w_{\nabla}$ $\displaystyle=$
$\displaystyle\frac{i}{2}(X^{\mu}Z^{\nu}-X^{\nu}Z^{\mu})\,w_{\mu\nu}=\frac{i\rho\,z_{1}z_{2}}{2Q^{2}\,(1+\rho^{2})}\,\Big{(}P^{\nu}R^{\mu}-P^{\mu}R^{\nu}\Big{)}\,w_{\mu\nu}\,.$
(14)

In the evaluation of the absorptive parts of the one loop diagrams we use the
following set of imaginary parts of scalar one-loop integrals
Lyubovitskij:2021ges ; Haug:2022hkr :

$\displaystyle{\rm Im}B_{0}(Q^{2})={\rm Im}B_{0}(\hat{s})=\pi\,,\ \ {\rm
Im}B_{0}(\hat{u})={\rm Im}B_{0}(\hat{t})=0\,,$ $\displaystyle{\rm
Im}C_{0}(\hat{s},0)=\frac{\pi}{\hat{s}}\biggl{(}\frac{1}{\bar{\epsilon}}-\log\frac{\hat{s}}{\mu^{2}}\biggr{)}\,,\
\ {\rm Im}C_{0}(\hat{u},0)={\rm Im}C_{0}(\hat{t},0)=0\,,$ $\displaystyle{\rm
Im}C_{0}(Q^{2},0)=\frac{\pi}{Q^{2}}\biggl{(}\frac{1}{\bar{\epsilon}}-\log\frac{Q^{2}}{\mu^{2}}\biggr{)}\,,\
\ {\rm
Im}C_{0}(Q^{2},\hat{s})=-\frac{\pi}{Q^{2}-\hat{s}}\,\log\frac{Q^{2}}{\hat{s}}\,,$
$\displaystyle{\rm
Im}C_{0}(Q^{2},\hat{u})=\frac{\pi}{Q^{2}-\hat{u}}\biggl{(}\frac{1}{\bar{\epsilon}}-\log\frac{Q^{2}}{\mu^{2}}\biggr{)}\,,\
\ {\rm
Im}C_{0}(Q^{2},\hat{t})=\frac{\pi}{Q^{2}-\hat{t}}\biggl{(}\frac{1}{\bar{\epsilon}}-\log\frac{Q^{2}}{\mu^{2}}\biggr{)}\,,$
$\displaystyle{\rm
Im}C_{0}(\hat{s},\hat{u})=\frac{\pi}{\hat{s}-\hat{u}}\biggl{(}\frac{1}{\bar{\epsilon}}-\log\frac{\hat{s}}{\mu^{2}}\biggr{)}\,,\
\ {\rm
Im}C_{0}(\hat{s},\hat{t})=\frac{\pi}{\hat{s}-\hat{t}}\biggl{(}\frac{1}{\bar{\epsilon}}-\log\frac{\hat{s}}{\mu^{2}}\biggr{)}\,,$
$\displaystyle{\rm
Im}D_{0}(Q^{2},\hat{s},\hat{u})=-\frac{2\pi}{\hat{s}\hat{u}}\log\frac{Q^{2}-\hat{u}}{Q^{2}}\,,\
\ {\rm
Im}D_{0}(Q^{2},\hat{s},\hat{t})=-\frac{2\pi}{\hat{s}\hat{t}}\log\frac{Q^{2}-\hat{t}}{Q^{2}}\,,$
$\displaystyle{\rm
Im}D_{0}(Q^{2},\hat{t},\hat{u})=\frac{2\pi}{\hat{u}\hat{t}}\biggl{(}-\frac{1}{\bar{\epsilon}}+\log\frac{Q^{2}}{\mu^{2}}-\log\frac{(Q^{2}-\hat{u})(Q^{2}-\hat{t})}{\hat{u}\hat{t}}\biggr{)}\,,$
(15)

where $B_{0}$ denotes the two-propagator (bubble) diagram, $C_{0}$ the three-
propagator (triangle) diagram, and $D_{0}$ the four-propagator (box) diagram.
We have used dimensional regularization with $D=4-2\epsilon$ space-time
dimensions; as usual
$1/\bar{\epsilon}=1/\epsilon+\log(4\pi)+\gamma_{\text{E}}$ with the Euler-
Mascheroni constant $\gamma_{E}$.

We note that in the cases of the electroweak contributions to
$w_{\Delta\Delta_{P}}$ and $w_{\Delta_{P}}$ we have to deal with an odd number
of $\gamma^{5}$ matrices in the relevant Dirac traces. We adopt the Larin
scheme for treating $\gamma_{5}$; details are discussed in Appendix B.

##  III Analytical results for the partonic $T$-odd structure functions

In this section we present our analytical results for the parton-level $T$-odd
structure functions $w_{\Delta\Delta_{P}}$, $w_{\Delta_{P}}$, and
$w_{\nabla}$. We first introduce some notation. We use the QCD color factors
$C_{F}=(N_{c}^{2}-1)/(2N_{c})=4/3$, $C_{A}=N_{c}=3$, $T_{F}=1/2$,
$C_{1}=C_{F}-N_{c}/2=-1/(2N_{c})=-1/6$, which at large $N_{c}$ scale as ${\cal
O}(N_{c})$, ${\cal O}(N_{c})$, ${\cal O}(1)$, ${\cal O}(1/N_{c})$,
respectively. Specifically, the color factors for $q\bar{q}$ annihilation and
$qg$ scattering are $C_{q\bar{q}}=C_{F}/N_{c}=(N_{c}^{2}-1)/(2N_{c}^{2})=4/9$
and $C_{qg}=T_{F}/N_{c}=1/(2N_{c})=1/6$. Furthermore, it is convenient to
introduce the coupling factors $g_{q\bar{q};i}=C_{q\bar{q}}\,g_{{\rm
EW};i}^{Z\gamma/W}\,e_{q}^{2}\,\alpha_{s}^{2}/(4\pi)$ and
$g_{qg;i}=C_{qg}\,g_{{\rm
EW};i}^{Z\gamma/W}\,e_{q}^{2}\,\alpha_{s}^{2}/(4\pi)$, where $i=1,2$. Here,
$e_{q}$ is the electric charge of a quark of flavor $q$. The electroweak
couplings $g_{{\rm EW};1}$ and $g_{{\rm EW};2}$, which incorporate the
products of couplings of the gauge bosons ($W^{\pm}$, $Z^{0}$, $\gamma$) with
quarks and leptons, are given by

$\displaystyle g_{\rm EW;1}^{Z\gamma}$ $\displaystyle=$ $\displaystyle
2\,g_{Zq}^{A}\,\Big{[}g_{Zq}^{V}\,\Big{(}(g_{Z\ell}^{V})^{2}+(g_{Z\ell}^{A})^{2}\Big{)}\,|D_{Z}(Q^{2})|^{2}+g_{Z\ell}^{V}\,\mathrm{Re}[D_{Z}(Q^{2})]\Big{]}\,,$
$\displaystyle g_{\rm EW;2}^{Z\gamma}$ $\displaystyle=$ $\displaystyle
2\,g_{Z\ell}^{A}\,\Big{[}g_{Z\ell}^{V}\,\Big{(}(g_{Zq}^{V})^{2}+(g_{Zq}^{A})^{2}\Big{)}\
|D_{Z}(Q^{2})|^{2}+g_{Zq}^{V}\,\mathrm{Re}[D_{Z}(Q^{2})]\Big{]}\,,$ (16)

in case of electrically neutral gauge bosons ($Z^{0}$, $\gamma$) and

$\displaystyle g_{\rm EW;1}^{W}$ $\displaystyle=$ $\displaystyle
2\,g_{Wqq^{\prime}}^{V}\,g_{Wqq^{\prime}}^{A}\,\Big{(}(g_{W\ell}^{V})^{2}+(g_{W\ell}^{A})^{2}\Big{)}|V_{qq^{\prime}}|^{2}\,|D_{W}(Q^{2})|^{2}\,,$
$\displaystyle g_{\rm EW;2}^{W}$ $\displaystyle=$ $\displaystyle
2\,g_{W\ell}^{V}\,g_{W\ell}^{A}\,\Big{(}(g_{Wqq^{\prime}}^{V})^{2}+(g_{Wqq^{\prime}}^{A})^{2}\Big{)}\,|V_{qq^{\prime}}|^{2}\,|D_{W}(Q^{2})|^{2}\,,$
(17)

in case of the $W^{\pm}$ gauge bosons, where

$\displaystyle g_{W\ell}^{V}$ $\displaystyle=$ $\displaystyle
g_{W\ell}^{A}\,=\,g_{Wqq^{\prime}}^{V}\,=\,g_{Wqq^{\prime}}^{A}\,=\,\frac{1}{2\sin\theta_{W}\sqrt{2}}\,,$
$\displaystyle g_{Z\ell}^{V}$ $\displaystyle=$
$\displaystyle-\frac{1-4\sin^{2}\theta_{W}}{2\sin 2\theta_{W}}\,,\quad\hskip
2.84544ptg_{Z\ell}^{A}\,=\,-\frac{1}{2\sin 2\theta_{W}}\,,$ $\displaystyle
g_{Zu}^{V}$ $\displaystyle=$
$\displaystyle\frac{1-8/3\sin^{2}\theta_{W}}{2e_{q}\sin 2\theta_{W}}\,,\quad
g_{Zd}^{V}\,=\,-\frac{1-4/3\sin^{2}\theta_{W}}{2e_{q}\sin 2\theta_{W}}\,,$
$\displaystyle g_{Zu}^{A}$ $\displaystyle=$ $\displaystyle\frac{1}{2e_{q}\sin
2\theta_{W}}\,,\quad\hskip 17.07182ptg_{Zd}^{A}\,=\,-\frac{1}{2e_{q}\sin
2\theta_{W}}\,.$ (18)

In the above expressions, $V_{qq^{\prime}}$ is the relevant element of the
Cabibbo-Kabayashi-Maskawa (CKM) matrix and $\theta_{W}$ is the Weinberg angle
measured to be $\sin^{2}\theta_{W}=0.23121$ ParticleDataGroup:2022pth .
Furthermore, $D_{G}(Q^{2})$ ($G=W^{\pm},Z^{0}$) denotes the product of the
Breit-Wigner propagator of a weak gauge boson and $Q^{2}$. Its real and
imaginary part are given by

$\displaystyle\mathrm{Re}[D_{G}(Q^{2})]$ $\displaystyle=$
$\displaystyle\frac{(Q^{2}-M_{G}^{2})Q^{2}}{(Q^{2}-M_{G}^{2})^{2}+M_{G}^{2}\Gamma_{G}^{2}}\,,$
$\displaystyle\mathrm{Im}[D_{G}(Q^{2})]$ $\displaystyle=$
$\displaystyle\frac{M_{G}\Gamma_{G}Q^{2}}{(Q^{2}-M_{G}^{2})^{2}+M_{G}^{2}\Gamma_{G}^{2}}\,.$
(19)

The masses $M_{G}$ and total widths $\Gamma_{G}$ of the bosons, taken from the
Particle Data Group ParticleDataGroup:2022pth , are $M_{W^{\pm}}=80.377\pm
0.012$ GeV, $M_{Z^{0}}=91.1876\pm 0.0021$ GeV, $\Gamma_{W^{\pm}}=2.085\pm
0.042$ GeV, and $\Gamma_{Z^{0}}=2.4955\pm 0.0023$ GeV. Note that in Eqs. (16)
and (17) the terms proportional to the squares of the Breit-Wigner propagators
($|D_{Z}(Q^{2})|^{2}$ and $|D_{W}(Q^{2})|^{2}$) correspond to the purely weak
$ZZ$ and $WW$ contributions to the couplings, while the terms proportional to
the real part of the Breit-Wigner $Z$ boson propagator
$\mathrm{Re}[D_{Z}(Q^{2})]$ correspond to $\gamma$-$Z$ interference.

As a final ingredient, we note that all one-loop partonic structure functions
contain a factor $\delta\big{(}(\hat{s}+\hat{t}+\hat{u}-Q^{2})/\hat{s}\big{)}$
arising from phase space and corresponding to the fact that the recoil in the
final state consists of a single massless parton. It is convenient to write

$\displaystyle
w^{ab}(z_{1},z_{2},\rho^{2})\,=\,\tilde{w}^{ab}(z_{1},z_{2},\rho^{2})\,\delta\big{(}(\hat{s}+\hat{t}+\hat{u}-Q^{2})/\hat{s}\big{)}\,.$
(20)

With this notation in place, we obtain the following partonic $T$-odd
structure functions. For the $q\bar{q}$ annihilation subprocess, we find

$\displaystyle\tilde{w}_{\Delta\Delta_{P}}^{q\bar{q}}$ $\displaystyle=$
$\displaystyle\frac{g_{q\bar{q};1}}{2}\,\sqrt{\frac{Q^{2}\hat{s}}{(Q^{2}-\hat{u})(Q^{2}-\hat{t})}}\,\biggl{[}-\frac{C_{F}}{2}\biggl{(}\frac{Q^{2}-\hat{t}}{Q^{2}-\hat{u}}+\frac{Q^{2}-\hat{u}}{Q^{2}-\hat{t}}\biggr{)}$
(21) $\displaystyle+$ $\displaystyle
C_{1}\biggl{(}\frac{Q^{2}-\hat{t}}{\hat{t}}\biggl{(}1-\frac{\hat{s}}{\hat{t}}\log\frac{Q^{2}-\hat{u}}{\hat{s}}\biggr{)}+\frac{Q^{2}-\hat{u}}{\hat{u}}\biggl{(}1-\frac{\hat{s}}{\hat{u}}\log\frac{Q^{2}-\hat{t}}{\hat{s}}\biggr{)}\biggr{]}\,,$
$\displaystyle\tilde{w}_{\Delta_{P}}^{q\bar{q}}$ $\displaystyle=$
$\displaystyle\frac{g_{q\bar{q};1}}{2}\,\,\frac{Q^{2}\hat{s}}{\sqrt{(Q^{2}-\hat{u})(Q^{2}-\hat{t}\,)\hat{u}\hat{t}}}\,\biggl{[}C_{F}\biggl{(}\frac{Q^{2}-\hat{t}}{Q^{2}-\hat{u}}-\frac{Q^{2}-\hat{u}}{Q^{2}-\hat{t}}\biggr{)}$
(22) $\displaystyle+$ $\displaystyle
C_{1}\biggl{(}\frac{Q^{2}-\hat{u}}{\hat{u}}\log\frac{Q^{2}-\hat{t}}{\hat{s}}-\frac{Q^{2}-\hat{t}}{\hat{t}}\log\frac{Q^{2}-\hat{u}}{\hat{s}}\biggr{)}\biggr{]}\,,$
$\displaystyle\tilde{w}_{\nabla}^{q\bar{q}}$ $\displaystyle=$ $\displaystyle
g_{q\bar{q};2}\,\sqrt{\frac{Q^{2}\hat{s}}{\hat{u}\hat{t}}}\,\biggl{[}\frac{C_{F}}{2}\frac{(2Q^{2}\hat{s}+\hat{u}\hat{t})(Q^{2}+\hat{s})(\hat{u}-\hat{t})}{(Q^{2}-\hat{u})^{2}(Q^{2}-\hat{t})^{2}}$
(23) $\displaystyle+$ $\displaystyle
C_{1}\biggl{(}-\frac{Q^{2}(\hat{u}-\hat{t})}{(Q^{2}-\hat{u})(Q^{2}-\hat{t})}+\frac{\hat{s}}{\hat{u}}\log\frac{Q^{2}-\hat{t}}{\hat{s}}-\frac{\hat{s}}{\hat{t}}\log\frac{Q^{2}-\hat{u}}{\hat{s}}\biggr{)}\biggr{]}\,.$

For $qg$ scattering we have

$\displaystyle\tilde{w}_{\Delta\Delta_{P}}^{qg}$ $\displaystyle=$
$\displaystyle\frac{g_{qg;1}}{2}\,\frac{\hat{u}}{\hat{s}}\,\sqrt{\frac{Q^{2}\hat{s}}{(Q^{2}-\hat{u})(Q^{2}-\hat{t})}}\,\biggl{[}\frac{C_{F}}{2}\frac{2Q^{2}+\hat{s}}{Q^{2}-\hat{t}}$
(24) $\displaystyle-$ $\displaystyle
C_{1}\biggl{(}\frac{Q^{2}}{Q^{2}-\hat{u}}+\frac{Q^{2}-\hat{u}}{\hat{s}}\log\frac{(Q^{2}-\hat{u})(Q^{2}-\hat{t})}{\hat{u}\hat{t}}+\frac{Q^{2}-\hat{t}}{\hat{t}}\biggl{(}1-\frac{Q^{2}-\hat{u}}{\hat{t}}\log\frac{Q^{2}-\hat{u}}{\hat{s}}\biggr{)}\biggr{)}\biggr{]}\,,$
$\displaystyle\tilde{w}_{\Delta_{P}}^{qg}$ $\displaystyle=$
$\displaystyle\frac{g_{qg;1}}{2}\,\frac{Q^{2}\hat{u}}{\sqrt{(Q^{2}-\hat{u})(Q^{2}-\hat{t})\hat{u}\hat{t}}}\,\biggl{[}\frac{C_{F}}{2}\frac{\hat{t}-2\hat{u}}{Q^{2}-\hat{t}}$
(25) $\displaystyle-$ $\displaystyle
C_{1}\biggl{(}\frac{\hat{u}}{Q^{2}-\hat{u}}+\frac{Q^{2}}{\hat{s}}-\frac{(Q^{2}-\hat{u})\hat{u}}{\hat{s}^{2}}\log\frac{(Q^{2}-\hat{u})(Q^{2}-\hat{t})}{\hat{u}\hat{t}}\biggr{)}\biggr{]}\,,$
$\displaystyle\tilde{w}_{\nabla}^{qg}$ $\displaystyle=$ $\displaystyle
g_{qg;2}\,\frac{\hat{u}}{\hat{s}}\,\sqrt{\frac{Q^{2}\hat{s}}{\hat{u}\hat{t}}}\,\biggl{[}-\frac{C_{F}}{2}\frac{2Q^{2}\hat{u}+\hat{s}\hat{t}}{(Q^{2}-\hat{t})^{2}}$
(26) $\displaystyle-$ $\displaystyle
C_{1}\biggl{(}\frac{2Q^{2}\hat{s}}{(Q^{2}-\hat{u})^{2}}-\frac{Q^{4}-\hat{u}\hat{t}}{(Q^{2}-\hat{u})(Q^{2}-\hat{t})}+\frac{\hat{u}}{\hat{t}}\log\frac{Q^{2}-\hat{u}}{\hat{s}}-\frac{\hat{u}}{\hat{s}}\log\frac{(Q^{2}-\hat{u})(Q^{2}-\hat{t})}{\hat{u}\hat{t}}\biggr{)}\biggr{]}\,.$

We note that in the large $N_{c}$ limit the terms proportional to $C_{1}$ in
the structure functions are suppressed by a factor $1/N_{c}^{2}$ relative to
those proportional to $C_{F}$ and $C_{A}$. The results for $g\bar{q}$
scattering are obtained from the ones for the $qg$ process by interchange of
momenta, $p_{1}\leftrightarrow p_{2}$, which corresponds to an interchange of
Mandelstam variables $u\leftrightarrow t$.

For calculating the hadronic structure functions via the factorization formula
(11) and subsequently investigating their small $Q_{T}$ behavior, it is
convenient to express our results in terms of the variables
$z_{1}=x_{1}/\xi_{1}$, $z_{2}=x_{2}/\xi_{2}$, and $\rho^{2}=Q^{2}/Q_{T}^{2}$.
Using

$\displaystyle\hat{s}=\frac{Q^{2}+Q_{T}^{2}}{z_{1}z_{2}}\,,\quad
Q^{2}-\hat{t}=\frac{Q^{2}+Q_{T}^{2}}{z_{1}}\,,\quad
Q^{2}-\hat{u}=\frac{Q^{2}+Q_{T}^{2}}{z_{2}}\,,$
$\displaystyle\frac{\hat{u}\hat{t}}{Q^{2}\hat{s}}=\rho^{2}\,,\quad\frac{(Q^{2}-\hat{u})(Q^{2}-\hat{t})}{Q^{2}\hat{s}}=1+\rho^{2}\,,$
(27)

one gets for $q\bar{q}$ annihilation

$\displaystyle\tilde{w}_{\Delta\Delta_{P}}^{q\bar{q}}$ $\displaystyle=$
$\displaystyle-\frac{g_{q\bar{q};1}}{4z_{1}z_{2}}\,\,\frac{1}{\sqrt{1+\rho^{2}}}\,\biggl{[}C_{A}\frac{z_{1}^{2}+z_{2}^{2}}{2}+C_{1}\Big{(}z_{1}^{2}F_{1}(z_{2})+z_{2}^{2}F_{1}(z_{1})\Big{)}\biggr{]}\,,$
(28) $\displaystyle\tilde{w}_{\Delta_{P}}^{q\bar{q}}$ $\displaystyle=$
$\displaystyle-\frac{g_{q\bar{q};1}}{2z_{1}z_{2}}\,\,\frac{1}{\rho\sqrt{1+\rho^{2}}}\biggl{[}C_{A}\frac{z_{1}^{2}-z_{2}^{2}}{2}+C_{1}\Big{(}z_{1}^{2}F_{2}(z_{2})-z_{2}^{2}F_{2}(z_{1})\Big{)}\biggr{]}\,,$
(29) $\displaystyle\tilde{w}_{\nabla}^{q\bar{q}}$ $\displaystyle=$
$\displaystyle-\frac{g_{q\bar{q};2}}{z_{1}z_{2}}\,\frac{1}{\rho}\biggl{[}\Big{(}C_{A}-\frac{\rho^{2}}{1+\rho^{2}}C_{F}\Big{)}\frac{z_{1}^{2}-z_{2}^{2}}{2}+C_{1}\Big{(}z_{1}F_{2}(z_{2})-z_{2}F_{2}(z_{1})\Big{)}\biggr{]}\,,$
(30)

and for $qg$ scattering

$\displaystyle\tilde{w}_{\Delta\Delta_{P}}^{qg}$ $\displaystyle=$
$\displaystyle-\frac{g_{qg;1}}{2}\,\frac{1-z_{2}}{z_{1}z_{2}}\,\frac{1}{\sqrt{1+\rho^{2}}}\,\biggl{[}\frac{C_{F}}{2}z_{1}\biggl{(}1+\frac{2z_{1}z_{2}}{1+\rho^{2}}\biggr{)}$
(31) $\displaystyle+$ $\displaystyle
C_{1}z_{1}z_{2}\biggl{(}\Big{(}F_{1}(z_{1})-\frac{1-\rho^{2}}{1+\rho^{2}}\Big{)}\frac{z_{2}}{2}+z_{1}\log\frac{\rho^{2}}{1+\rho^{2}}\biggr{)}\biggr{]}\,,$
$\displaystyle\tilde{w}_{\Delta_{P}}^{qg}$ $\displaystyle=$
$\displaystyle-\frac{g_{qg;1}}{2}\,\frac{1-z_{2}}{z_{1}z_{2}}\,\frac{1}{\rho\sqrt{1+\rho^{2}}}\,\biggl{[}\frac{C_{F}}{2}z_{1}(1+z_{1}-2z_{2})$
(32) $\displaystyle+$ $\displaystyle
C_{1}z_{2}\biggl{(}1-z_{2}-\frac{z_{1}^{2}z_{2}}{1+\rho^{2}}+z_{1}^{2}(1-z_{2})\log\frac{\rho^{2}}{1+\rho^{2}}\biggr{)}\biggr{]}\,,$
$\displaystyle\tilde{w}_{\nabla}^{qg}$ $\displaystyle=$ $\displaystyle-
g_{qg;2}\frac{1-z_{2}}{z_{1}z_{2}}\,\frac{1}{\rho}\biggl{[}C_{F}z_{1}\biggl{(}\frac{z_{1}(1-z_{2})}{1+\rho^{2}}+\frac{1-z_{1}}{2z_{2}}\biggr{)}$
(33) $\displaystyle+$ $\displaystyle
C_{1}z_{2}\biggl{(}\frac{z_{1}^{2}-2z_{2}}{1+\rho^{2}}+z_{1}(1-z_{2})\log\frac{\rho^{2}}{1+\rho^{2}}-(1-z_{2})F_{2}(z_{1})\biggr{)}\biggr{]}\,.$

Here, the functions $F_{1}$ and $F_{2}$ are defined as:

$\displaystyle F_{1}(z)$ $\displaystyle\equiv$
$\displaystyle\frac{1+z}{1-z}+\frac{2z\log(z)}{(1-z)^{2}}=2\sum\limits_{N=1}^{\infty}\,\frac{(1-z)^{N}}{(N+1)(N+2)}={\cal
O}(1-z)\,,$ $\displaystyle F_{2}(z)$ $\displaystyle\equiv$ $\displaystyle
1+\frac{z\log(z)}{1-z}=\frac{1-z}{2}\,\Big{(}1+F_{1}(z)\Big{)}=\sum\limits_{N=1}^{\infty}\,\frac{(1-z)^{N}}{N(N+1)}={\cal
O}(1-z)\,.$ (34)

They obey $F_{1}(1)=F_{2}(1)=0$ and $F_{1}(0)=F_{2}(0)=1$. For later
reference, we have also given their expansions around $z=1$.

We now have

$\displaystyle
W(x_{1},x_{2},\rho^{2})=\frac{1}{x_{1}x_{2}}\,\sum\limits_{a,b}\,\int\limits_{x_{1}}^{1}\,dz_{1}\,\int\limits_{x_{2}}^{1}\,dz_{2}\
\tilde{w}^{ab}(z_{1},z_{2},\rho^{2})\,\delta\left((1-z_{1})(1-z_{2})-\frac{\rho^{2}}{1+\rho^{2}}z_{1}z_{2}\right)\,f_{a/H_{1}}\Big{(}\frac{x_{1}}{z_{1}}\Big{)}\,f_{b/H_{2}}\Big{(}\frac{x_{2}}{z_{2}}\Big{)}\,.\quad$
(35)

This factorization is formally valid when $Q_{T}$ is of order $Q$, that is,
for $\Lambda_{\text{QCD}}\ll Q\sim Q_{T}$. For $Q_{T}\ll Q$ the appropriate
factorization formalism is TMD factorization. Here, we will take the collinear
factorization and extrapolate to small values of $Q_{T}$ by formally expanding
the result about $Q_{T}=0$. This will result in an expansion in powers of
$\rho^{2}$ which can be matched to TMD results in the region of intermediate
$Q_{T}$ for a smooth transition from the TMD to the collinear regime. We will
perform the expansion in $\rho^{2}$ beyond the leading power, which has been
the main focus in the existing literature, to provide information about which
higher-power corrections are accounted for in the collinear formalism.

##  IV Small-$Q_{T}$ expansion

In the expansion of hadronic structure functions of the form in Eq. (35) we
have three contributions. First, there is the direct dependence of the
partonic structure function on $Q_{T}$. Second, the phase space delta function
has non-trivial $Q_{T}$ dependence. Third, the variables $x_{1}$, $x_{2}$ have
implicit $Q_{T}$ dependence. The first type of contribution may be
straightforwardly taken into account by simple expansion of the partonic
structure functions. The second and third contributions require more
discussion.

The phase space delta function in (35) is well-known in the literature, see,
e.g., Refs. Meng:1995yn ; Boer:2006eq ; Berger:2007jw . Its expansion to
leading power in $\rho^{2}=Q_{T}^{2}/Q^{2}$ was also given in that reference
and reads

$\displaystyle\delta\\!\left((1-z_{1})(1-z_{2})-\frac{\rho^{2}}{1+\rho^{2}}z_{1}z_{2}\right)=\frac{\delta(1-z_{1})}{(1-z_{1})_{+}}+\frac{\delta(1-z_{2})}{(1-z_{2})_{+}}-\delta(1-z_{1})\delta(1-z_{2})\log\rho^{2}+\mathcal{O}(\rho^{2})$
$\displaystyle\,.$ (36)

Here the “plus”-distribution is defined by

$\int_{0}^{1}dz\,\frac{f(z)}{(1-z)_{+}}\equiv\int_{0}^{1}dz\,\frac{f(z)-f(1)}{1-z}\,,$
(37)

for a function $f$ that is regular at $z=1$. We note in passing that in Ref.
Gelfand  a general method for the expansion of distributions was developed,
based on Mellin integral techniques. Building on these ideas, we recently
formulated LVWZ:2023  an algorithm for the small-$Q_{T}$ expansion of singular
functions valid to arbitrary order of $\rho^{2}$ and arbitrary number of
radiated partons. This will be presented in a separate publication. Here we
are only concerned with the expansion of integrals containing the phase space
delta function in (36). For a general regular function $\varphi(z_{1},z_{2})$
such an integral can be expanded for small $Q_{T}$ including
$\mathcal{O}(\rho^{4},\rho^{4}\log\rho^{2})$ terms in the following way:

$\displaystyle I_{0}$ $\displaystyle\equiv$
$\displaystyle\int\limits_{x_{1}}^{1}dz_{1}\int\limits_{x_{2}}^{1}dz_{2}\,\delta\Big{(}(1-z_{1})(1-z_{2})-\frac{\rho^{2}}{1+\rho^{2}}z_{1}z_{2}\Big{)}\,\varphi(z_{1},z_{2})$
(38) $\displaystyle=$
$\displaystyle\int\limits_{x_{1}}^{1}dz_{1}\int\limits_{x_{2}}^{1}dz_{2}\,\biggl{(}\delta(1-z_{2})G_{1}(z_{1},z_{2})+\delta(1-z_{1})G_{1}(z_{2},z_{1})$
$\displaystyle+$
$\displaystyle\delta(1-z_{1})\delta(1-z_{2})G_{2}(z_{1},z_{2})\biggr{)}\,\varphi(z_{1},z_{2})+{\cal
O}(\rho^{6},\rho^{6}\log\rho^{2})\,,$

with

$\displaystyle G_{1}(z_{1},z_{2})$ $\displaystyle=$
$\displaystyle\frac{(1+\rho^{2})(1+\rho^{2}\partial_{z_{2}})+\rho^{4}\partial^{2}_{z_{2}}/2}{(1-z_{1})_{+}}$
$\displaystyle-$
$\displaystyle\frac{\rho^{2}(1+\rho^{2}+(1+3\rho^{2})\partial_{z_{2}}+\rho^{2}\partial^{2}_{z_{2}})}{(1-z_{1})^{2}_{+,1}}+\frac{\rho^{4}(1+2\partial_{z_{2}}+\partial^{2}_{z_{2}}/2)}{(1-z_{1})^{3}_{+,2}}\,,$
$\displaystyle G_{2}(z_{1},z_{2})$ $\displaystyle=$
$\displaystyle\rho^{2}\biggl{(}1+\rho^{2}\Big{(}\frac{1}{2}+\partial_{z_{1}}+\partial_{z_{2}}+\partial_{z_{1}z_{2}}^{2}\Big{)}\biggr{)}-\log\rho^{2}\biggl{(}(1+\rho^{2})(1+\rho^{2}(\partial_{z_{1}}+\partial_{z_{2}}+\partial^{2}_{z_{1}z_{2}}))$
(39) $\displaystyle+$
$\displaystyle\frac{\rho^{4}}{2}\,(\partial_{z_{1}}+\partial_{z_{2}})^{2}+\rho^{4}\partial^{2}_{z_{1}z_{2}}\,\Big{(}1+\partial_{z_{1}}+\partial_{z_{2}}+\frac{\partial^{2}_{z_{1}z_{2}}}{4}\Big{)}\biggr{)}\,.$

Here $1/(1-z)^{m}_{+,m-1}$ is a generalized plus distribution of power $m$,
defined by

$\int_{0}^{1}dz\frac{f(z)}{(1-z)^{m}_{+,m-1}}\equiv\int_{0}^{1}dz\frac{f(z)-\mathcal{T}^{m-1}_{z=1}f(z)}{(1-z)^{m}}\,,$
(40)

where $f(z)$ is again a sufficiently regular test function and
$\mathcal{T}^{m-1}_{z=1}f(z)$ denotes the Taylor polynomial of $f(z)$ about
$z=1$ to order $m-1$,

$\mathcal{T}^{m-1}_{z=1}f(z)=\sum_{k=0}^{m-1}\frac{(-1)^{k}f^{(k)}(1)}{k!}\,(1-z)^{k}.$
(41)

A lower integration bound of $x$ instead of zero introduces additional
boundary terms of the form

$\displaystyle\int\limits_{x}^{1}dz\,\frac{f(z)}{(1-z)^{m}_{+,m-1}}$
$\displaystyle=$
$\displaystyle\int\limits_{x}^{1}dz\,\biggl{[}\frac{1}{(1-z)^{m}_{+x,m-1}}+\delta(1-z)\,\log(1-x_{1})\,\frac{(-1)^{m-1}}{(m-1)!}\,\partial_{z}^{m-1}$
(42) $\displaystyle-$
$\displaystyle\delta(1-z)\,\sum\limits_{j=2}^{m}\,\frac{(-1)^{m-j}}{(j-1)\,(m-j)!}\,\biggl{(}\frac{1}{(1-x_{1})^{j-1}}-1\biggr{)}\,\partial_{z}^{m-j}\biggr{]}f(z_{1})\,,$

where $f(z)/(1-z)^{m}_{+x,m-1}$ is the generalized plus distribution defined
for an integral starting at a finite lower limit $x$, i.e.

$\int_{x}^{1}dz\frac{f(z)}{(1-z)^{m}_{+x,m-1}}\equiv\int_{x}^{1}dz\frac{f(z)-\mathcal{T}^{m-1}_{z=1}f(z)}{(1-z)^{m}}\,.$
(43)

Comparing Eq. (38) with (36) one can see that we reproduce the known leading
terms, while the terms of order $\rho^{2}$, $\rho^{2}\log\rho^{2}$,
$\rho^{4}$, and $\rho^{4}\log\rho^{2}$ are new.

Substituting the small-$Q_{T}$ expansion of the parton-level structure
functions $w^{ab}(z_{1},z_{2},\rho^{2})$ for the various partonic channels
into Eq. (11) we get for the contributions to the small-$Q_{T}$ expansion of
the hadronic structure function $W(x_{1},x_{2},\rho^{2})$:

$\displaystyle
W_{\delta-\textrm{fct.}}(x_{1},x_{2},\rho^{2})=W_{0}(x_{1},x_{2},L_{\rho})+\rho^{2}\,W_{1}(x_{1},x_{2},L_{\rho})+\rho^{4}\,W_{2}(x_{1},x_{2},L_{\rho})+{\cal
O}(\rho^{6})\,,$ (44)

were we have abbreviated

$\displaystyle L_{\rho}\equiv\log\rho^{2}\,.$ (45)

The expansion coefficients with $i=1,2,3$ have the structure

$\displaystyle W_{i}(x_{1},x_{2},L_{\rho})$ $\displaystyle=$
$\displaystyle\frac{1}{x_{1}x_{2}}\sum_{a,b}\,\biggl{[}R_{ab,i}(x_{1},x_{2},L_{\rho})\,f_{a/H_{1}}(x_{1})\,f_{b/H_{2}}(x_{2})$
(46) $\displaystyle+$ $\displaystyle\Big{(}P_{ba,i}\otimes
f_{b/H_{2}}\Big{)}(x_{2},x_{1},L_{\rho})\
f_{a/H_{1}}(x_{1})+\Big{(}P_{ab,i}\otimes
f_{a/H_{1}}\Big{)}(x_{1},x_{2},L_{\rho})\ f_{b/H_{2}}(x_{2})\biggr{]}\,,$

where

$\displaystyle\big{(}{\cal P}\otimes
f\big{)}(x,y,L_{\rho})=\int_{x}^{1}\,\frac{dz}{z}\,{\cal
P}(z,y,L_{\rho})\,f\Big{(}\frac{x}{z}\Big{)}$ (47)

denotes a generalized convolution, $R_{i}(x_{1},x_{2},L_{\rho})$,
$P_{ba,i}(z_{2},x_{1},L_{\rho})$, and $P_{ab,i}(z_{1},x_{2},L_{\rho})$ are
perturbative coefficient functions containing differential operators acting on
the PDFs $f_{a/H_{1}}(x_{1})$ and $f_{b/H_{2}}(x_{2})$. We note that the
generalized convolution (47) reverts to the ordinary one,

$\displaystyle\big{(}{\cal P}\otimes
f\big{)}(x)=\int_{x}^{1}\,\frac{dz}{z}\,{\cal
P}(z)\,f\Big{(}\frac{x}{z}\Big{)}\,,$ (48)

when ${\cal P}(z,y,L_{\rho})$ does not depend on $y$ and $L_{\rho}$. Details
are given in Appendix C. We stress that, as indicated in Eq. (44), the
functions $W_{i}$ may carry dependence on $\log\rho^{2}$, on top of the
overall power of $\rho$ that they multiply.

However, Equation (44) is not yet the complete expansion. As mentioned above,
we need to take into account that $x_{1}$ and $x_{2}$ are defined at finite
$Q_{T}$ (see Eq. (9)) and hence must also be expanded about their respective
values at $Q_{T}=0$, $x_{1}^{0}$ and $x_{2}^{0}$ in (10). Therefore, we
substitute $x_{i}=x_{i}^{0}\sqrt{1+\rho^{2}}$ as arguments of the structure
functions $W_{i}$ and perform the $\rho^{2}$ expansions of the latter. We now
present our final result for the full small-$Q_{T}$ expansion of the hadronic
structure functions, including the leading-power (LP) term $W^{\rm
LP}(x_{1}^{0},x_{2}^{0},L_{\rho})$, the next-to-leading-power (NLP) term
$W^{\rm NLP}(x_{1}^{0},x_{2}^{0},L_{\rho})$, and the next-next-to-leading-
power (NNLP)term $W^{\rm NNLP}(x_{1}^{0},x_{2}^{0},L_{\rho})$:

$\displaystyle W(x_{1},x_{2},\rho^{2})=W^{\rm
LP}(x_{1}^{0},x_{2}^{0},L_{\rho})+\rho^{2}\,W^{\rm
NLP}(x_{1}^{0},x_{2}^{0},L_{\rho})+\rho^{4}\,W^{\rm
NNLP}(x_{1}^{0},x_{2}^{0},L_{\rho})+{\cal O}(\rho^{6})\,,$ (49)

where

$\displaystyle W^{\rm LP}(x_{1}^{0},x_{2}^{0},L_{\rho})$ $\displaystyle=$
$\displaystyle W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})\,,$ (50) $\displaystyle
W^{\rm NLP}(x_{1}^{0},x_{2}^{0},L_{\rho})$ $\displaystyle=$ $\displaystyle
W_{1}(x_{1}^{0},x_{2}^{0},L_{\rho})+\frac{1}{2}\biggl{(}x_{1}^{0}\
\partial_{x_{1}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})+x_{2}^{0}\
\partial_{x_{2}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})\biggr{)}\,,$ (51)
$\displaystyle W^{\rm NNLP}(x_{1}^{0},x_{2}^{0},L_{\rho})$ $\displaystyle=$
$\displaystyle
W_{2}(x_{1}^{0},x_{2}^{0},L_{\rho})+\frac{1}{4}x_{1}^{0}x_{2}^{0}\
\partial_{x_{1}^{0}}\partial_{x_{2}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})$
(52) $\displaystyle-$ $\displaystyle\frac{1}{8}\biggl{(}x_{1}^{0}\
\partial_{x_{1}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})-4x_{1}^{0}\
\partial_{x_{1}^{0}}W_{1}(x_{1}^{0},x_{2}^{0},L_{\rho})-(x_{1}^{0})^{2}\
\partial^{2}_{x_{1}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})\biggr{)}$
$\displaystyle-$ $\displaystyle\frac{1}{8}\biggl{(}x_{2}^{0}\
\partial_{x_{2}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})-4x_{2}^{0}\
\partial_{x_{2}^{0}}W_{1}(x_{1}^{0},x_{2}^{0},L_{\rho})-(x_{2}^{0})^{2}\
\partial^{2}_{x_{2}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})\biggr{)}\,.$

Here $\partial^{m}_{x_{1}}\partial^{n}_{x_{2}}W_{i}(x_{1},x_{2},L_{\rho})$
denotes the $m$-th partial derivative with respect to $x_{1}$ and the $n$-th
partial derivative with respect to $x_{2}$. The calculational techniques for
taking these derivatives are discussed in Appendix C.

Explicitly we obtain the following analytical results for the LP contributions
$W_{J}^{{\rm LP};ab}(x_{1}^{0},x_{2}^{0},L_{\rho})$ to the $T$-odd hadronic
structure functions (here $ab=q\bar{q},qg$ and
$J=\Delta\Delta_{P},\Delta_{P},\nabla$):

$\displaystyle W^{{\rm
LP};q\bar{q}}_{\Delta\Delta_{P}}(x_{1}^{0},x_{2}^{0},L_{\rho})$
$\displaystyle=$
$\displaystyle\frac{g_{q\bar{q};1}}{4x_{1}^{0}x_{2}^{0}}\,C_{A}\,\left(L_{\rho}+\frac{3}{2}\right)\,q_{1}(x_{1}^{0})\,\bar{q}_{2}(x_{2}^{0})$
(53) $\displaystyle-$
$\displaystyle\frac{g_{q\bar{q};1}}{4x_{1}^{0}x_{2}^{0}}\,\frac{C_{A}}{2C_{F}}\,\Big{[}q_{1}(x_{1}^{0})\,\Big{(}P_{qq}\otimes\bar{q}_{2}\Big{)}(x_{2}^{0})+\Big{(}P_{qq}\otimes
q_{1}\Big{)}(x_{1}^{0})\,\bar{q}_{2}(x_{2}^{0})\Big{]}$ $\displaystyle-$
$\displaystyle\frac{g_{q\bar{q};1}}{4x_{1}^{0}x_{2}^{0}}\,C_{1}\,\Big{[}q_{1}(x_{1}^{0})\,\Big{(}f_{1}\otimes\bar{q}_{2}\Big{)}(x_{2}^{0})+\Big{(}f_{1}\otimes
q_{1}\Big{)}(x_{1}^{0})\,\bar{q}_{2}(x_{2}^{0})\Big{]}\,,$ $\displaystyle
W^{{\rm LP};q\bar{q}}_{\nabla}(x_{1}^{0},x_{2}^{0},L_{\rho})$ $\displaystyle=$
$\displaystyle 2\,\beta\,W^{{\rm
LP};q\bar{q}}_{\Delta_{P}}(x_{1}^{0},x_{2}^{0})$ (54) $\displaystyle=$
$\displaystyle-\frac{g_{q\bar{q};1}}{\rho\,x_{1}^{0}x_{2}^{0}}\,\frac{C_{A}}{2C_{F}}\,\Big{[}q_{1}(x_{1}^{0})\,\Big{(}\tilde{P}_{qq}\otimes\bar{q}_{2}\Big{)}(x_{2}^{0})-\Big{(}\tilde{P}_{qq}\otimes
q_{1}\Big{)}(x_{1}^{0})\,\bar{q}_{2}(x_{2}^{0})\Big{]}$ $\displaystyle-$
$\displaystyle\frac{g_{q\bar{q};1}}{\rho\,x_{1}^{0}x_{2}^{0}}\,C_{1}\,\Big{[}q_{1}(x_{1}^{0})\,\Big{(}f_{2}\otimes\bar{q}_{2}\Big{)}(x_{2}^{0})-\Big{(}f_{2}\otimes
q_{1}\Big{)}(x_{1}^{0})\,\bar{q}_{2}(x_{2}^{0})\Big{]}\,,$ $\displaystyle
W^{{\rm LP};qg}_{\Delta\Delta_{P}}(x_{1}^{0},x_{2}^{0},L_{\rho})$
$\displaystyle=$
$\displaystyle-\frac{g_{qg;1}}{4x_{1}^{0}x_{2}^{0}}\,q(x_{1}^{0})\,\Big{(}P^{\prime}_{qg}\otimes
g\Big{)}(x_{2}^{0})\,,$ (55) $\displaystyle W^{{\rm
LP};qg}_{\nabla}(x_{1}^{0},x_{2}^{0},L_{\rho})$ $\displaystyle=$
$\displaystyle 2\,\beta\,W^{{\rm
LP};qg}_{\Delta_{P}}(x_{1}^{0},x_{2}^{0},L_{\rho})\,=\,-\frac{g_{qg;1}}{\rho\,x_{1}^{0}x_{2}^{0}}\,q(x_{1}^{0})\,\Big{(}P^{\prime\prime}_{qg}\otimes
g\Big{)}(x_{2}^{0})\,,$ (56)

where

$\displaystyle\beta$ $\displaystyle=$ $\displaystyle\frac{g_{\rm EW;1}}{g_{\rm
EW;2}}\,,$ $\displaystyle P_{qq}(z)$ $\displaystyle=$ $\displaystyle
C_{F}\biggl{[}\frac{1+z^{2}}{(1-z)_{+}}+\frac{3}{2}\delta(1-z)\biggr{]}\,,$
$\displaystyle\tilde{P}_{qq}(z)$ $\displaystyle=$ $\displaystyle
C_{F}(1+z)\,,$ $\displaystyle P^{\prime}_{qg}(z,L_{\rho})$ $\displaystyle=$
$\displaystyle C_{F}(1+2z)+C_{1}z\,(2L_{\rho}-z)\,,$ $\displaystyle
P^{\prime\prime}_{qg}(z,L_{\rho})$ $\displaystyle=$ $\displaystyle
C_{F}(1-z)+C_{1}z\,(L_{\rho}(1-z)+1-2z)\,,$ (57)

and

$\displaystyle f_{i}(z)=\frac{F_{i}(z)}{1-z},$ (58)

with $F_{i}$ as defined in Eq. (III). Note that the $f_{i}$ are regular
functions with $f_{1}(1)=1/3$, $f_{2}(1)=1/2$.

As shown in (56), there is an interesting relation between the structure
functions $W^{{\rm LP};ab}_{\nabla}(x_{1}^{0},x_{2}^{0},L_{\rho})$ and
$W^{{\rm LP};ab}_{\Delta_{P}}(x_{1}^{0},x_{2}^{0},L_{\rho})$:

$\displaystyle W^{{\rm
LP};ab}_{\nabla}(x_{1}^{0},x_{2}^{0},L_{\rho})=2\beta\,W^{{\rm
LP};ab}_{\Delta_{P}}(x_{1}^{0},x_{2}^{0},L_{\rho})\,,$ (59)

valid both for the $q\bar{q}$ and the $qg$ subprocess at leading power. The
NLP and NNLP contributions to the $T$-odd hadronic structure function are
listed in the supplemental material to this paper.

To illustrate the numerical behavior of these expansions, we consider the
$q\bar{q}$ contribution to the hadronic double-flip structure function,
$W^{q\bar{q}}_{\Delta\Delta_{P}}(x_{1},x_{2})$, as an example. In Fig. 4 we
compare the full expression without $Q_{T}$ expansion with the LP, NLP, and
NNLP results. Here we use the CTEQ 6.1M PDFs of Ref. Stump:2003yu , taken from
LHAPDF Buckley:2014ana , along with their ManeParse Clark:2016jgm  Mathematica
implementation. We choose $\sqrt{s}=8$ TeV and $Q=100$ GeV as representative
of the kinematics in the ATLAS measurements ATLAS:2016rnf . As one can see,
the LP piece describes the full result only at low $Q_{T}$ and rapidly departs
from it for $Q_{T}>10$ GeV or $\rho^{2}>0.01$. By contrast, already inclusion
of the NLP term leads to excellent agreement with the full result out to
$Q_{T}=40$ GeV ($\rho^{2}=0.16$), only marginally further improved by the NNLP
contribution. In particular, for $Q_{T}=20$ GeV, the LP result deviates from
the full one by about 20%, whereas at NNLP the relative deviation is only
$\sim 0.4\%$.

Figure 4: Comparison of the full analytical result for the quark-channel
contribution to $W_{\Delta\Delta_{P}}$ (black solid line, taken from Eq. (30))
with expansions to LP (dashed), NLP (dot-dashed), NNLP (blue solid) as given
in the supplemental material.

##  V Comparison to ATLAS data

As shown in Eq. (4), the Drell-Yan cross section can be expressed in terms of
eight angular coefficients $A_{i=0,\ldots,7}$. The relations of these
coefficients to the hadronic helicity structure functions are recalled in
Appendix A. Previous experimental and phenomenological studies mostly focused
on the first five $A_{i}$ coefficients CMS:2015cyj ; ATLAS:2016rnf ;
NA10:1986fgk ; LHCb:2022tbc ; CDF:2011ksg ; Lambertsen:2016wgj ; Gehrmann-
DeRidder:2016cdi ; Richter-Was:2016avq ; Mirkes:1990vn  which are related to
the $T$-even structure functions. The $T$-odd structure functions have
received less attention. Experimentally, it has not yet been possible to
measure the $T$-odd angular coefficients in $W$ boson production, but results
in neutral-current scattering in the vicinity of the $Z$-boson mass peak are
available from the ATLAS Collaboration ATLAS:2016rnf . Specifically, the
measurement was performed in the $Z$-boson invariant mass window of 80–100
GeV, as a function of $Q_{T}$, and also in three bins of rapidity $y$: (a)
$|y|<1,(b)1<|y|<2$, and (c) $2<|y|<3.5$. We note that near $Q=m_{Z}$ the
contribution by $\gamma$-$Z^{0}$ interference is suppressed relative to that
for pure $Z^{0}$ exchange. In the following, we compare our results for the
angular coefficients $A_{5},A_{6}$, and $A_{7}$ to the ATLAS data. Here we use
the full expressions at ${\cal O}(\alpha_{s}^{2})$ for the helicity structure
functions $W_{\Delta\Delta_{P}},W_{\Delta_{P}}$, and $W_{\nabla}$. In the
denominator of the coefficients, we use the ${\cal O}(\alpha_{s})$ expressions
for the transverse and longitudinal structure functions $W_{T}$ and $W_{L}$
(see details in Refs. Altarelli:1984pt ; Boer:2006eq ; Berger:2007jw ). This
approach thus consistently gives the leading contribution to $A_{5},A_{6}$,
and $A_{7}$, which is of order $\alpha_{s}$.

To begin with, we investigate the rapidity and $Q_{T}$ dependences of the
angular coefficients $A_{5}$, $A_{6}$, $A_{7}$ near the $Z^{0}$ pole. Fig. 5
shows the rapidity for fixed $Q_{T}$, while Fig. 7 presents the results as
functions of $Q_{T}$ in one of the three rapidity bins accessed by ATLAS. As
the plots show, the coefficients are overall small, reaching at most 0.1–0.2%
near $Q_{T}\sim m_{Z}$ or towards larger $|y|$. This finding does not really
come as a surprise: Small values of the $T$-odd $A_{5,6,7}$ coefficients have
been predicted in Refs. Mirkes:1992hu ; Hagiwara:1984hi  also for $W$ boson
production. In case of $Z$ boson production $A_{5,6,7}$ are further suppressed
because of the smallness of the corresponding weak couplings, relative to the
couplings appearing in $W_{T}$ and $W_{L}$.

Figure 5: $T$-odd angular coefficients $A_{5}$, $A_{6}$ and $A_{7}$ for
various transverse momenta $Q_{T}$ of the lepton pair, as functions of pair
rapidity at $Q\sim m_{Z}$ and $\sqrt{s}=8$ TeV.

The increase of $A_{5,6,7}$ at larger values of rapidity in Fig. 5 – which is
consistent with the leading-power relation between $W_{\nabla}$ and
$W_{\Delta_{P}}$ we found in Eq. (59) – follows the trend observed in Ref.
Hagiwara:1984hi  for $W$ boson production in $p\bar{p}$ collisions at
$\sqrt{s}=540$ GeV. Likewise, a similar dependence on rapidity was found for
the angular coefficients in Refs. ATLAS:2016rnf ; CMS:2015cyj . At fixed
rapidity, the $T$-odd coefficients are small for small $Q_{T}$ and then
increase, peaking when $Q_{T}$ is near the $Z$ mass (see Fig. 7). We also show
in the figure the individual contributions by $q\bar{q}$ and $qg$ scattering,
the latter dominating for all kinematics. Among the $T$-odd structure
functions, $W_{\Delta\Delta_{P}}$, being symmetric under interchange
$z_{1}\leftrightarrow z_{2}$, has the largest contribution from quark-
antiquark annihilation.

Figure 6: $Q_{T}$ dependence of the angular coefficients $A_{5}$, $A_{6}$, and
$A_{7}$ for the rapidity interval $1<|y|<2$ at $Q\sim m_{Z}$ and $\sqrt{s}=8$
TeV. We also show the individual $q\bar{q}$ and $qg$ contributions.

Figure 7: $Q_{T}$ dependence of the angular coefficients $A_{5}$, $A_{6}$, and
$A_{7}$ for rapidity $1<|y|<2$ at $\sqrt{s}=8$ TeV and for different $Q$.

Figure 8: Comparison of our result for $A_{5}$ to ATLAS data ATLAS:2016rnf  at
$Q\sim m_{Z}$. The black and red experimental points denote the unregularized
and regularized data, respectively, and show their statistical error. The left
panel shows the results for $1<|y|<2$, while the right is for $2<|y|<3.5$.
(The scale on the $y$-axis has been chosen for better visibility of the small
values of the coefficient; as a result, some data points fall outside the plot
range. The full data set is shown in the left lower inset in each plot.)

Figure 9: Same as Fig. 10, but for the coefficient $A_{6}$ for rapidity
$1<|y|<2$.

Figure 10: Same as Fig. 10, but for the coefficient $A_{7}$.

Figure 7 explores the role played by the pair mass $Q$ for the $Q_{T}$
distribution. We show results for $Q=80$ GeV, $Q=m_{Z}$, and $Q=100$ GeV,
which span the range of $Q$ used for the ATLAS measurements. As one can see,
$A_{5}$ and $A_{6}$ are rather insensitive to $Q$, whereas $A_{7}$ exhibits a
strong dependence, even turning negative at high $Q$.

We now turn to the actual comparison with the ATLAS data ATLAS:2016rnf . Of
their three rapidity bins ($|y|<1,1<|y|<2$, and $2<|y|<3.5$) we only use the
two with higher $|y|$ since, as we saw above, the angular coefficients are
very small for $|y|<1$. We note that ATLAS presents the data in two ways, as
an “unregularized” and a “regularized” set. The regularization smoothes the
data by correcting for bin migration. This procedure involves the use of
Monte-Carlo pseudo data which are at lowest-order accuracy. Details about the
data regularization method are presented in Appendices C and E of Ref.
ATLAS:2016rnf .

Figures 10-10 show the comparison. (Our definition for the coefficient $A_{7}$
differs from that in Ref. ATLAS:2016rnf  by a factor $2$, which has been
properly taken into account in Fig. 10.) Neither the regularized nor the
unregularized data are in particularly good agreement with our theoretical
predictions. At best, there is qualitative agreement in that the theoretical
results show positive values for all three $T$-odd angular coefficients, with
a similar trend in the data. In particular, for the bin $1<|y|<2$ one observes
a rise of $A_{5,6,7}$ with $Q_{T}$ up to about $Q_{T}\sim m_{Z}$, exactly as
predicted theoretically. Quantitatively, however, the regularized data – which
is the set primarily to be used for comparisons – shows overall much higher
coefficients than obtained in our calculation. We note that ATLAS used the
DYNNLO package Catani:2009sm  to obtain theoretical results for $A_{5,6,7}$.
DYNNLO predicts values of up to $0.005$ for the coefficients $A_{5,6,7}$.
However, the prediction of non-zero values is at the limit of sensitivity of
both the theoretical calculation and the data.

More robust conclusions regarding the agreement of data and theory for the
$T$-odd structure functions will only become possible after significant
improvements of both the experimental results and the theoretical description.
For the former, it will be important to decrease the statistical uncertainty
of the data and to reduce the reliance on the regularization procedure. On the
theoretical side, predictions beyond the leading order would be desirable in
order to reduce the theoretical uncertainty. We note that the next-order
corrections could be obtained from Ref. Gehrmann-DeRidder:2016cdi . As we
mentioned earlier, the numerators and denominators of the coefficients
$A_{5,6,7}$ enter with different powers of $\alpha_{s}$, which likely makes
them sensitive to scale variations and higher-order corrections. As a result,
it is hard to estimate the theoretical uncertainty of our predictions.

##  VI Conclusion

We have performed a detailed analysis of perturbative $T$-odd effects in
charged- and neutral-current Drell-Yan processes, taking into account
$W^{\pm}$ and $Z^{0}$ exchange, as well as $\gamma$-$Z^{0}$ interference. To
this end, we have computed the relevant $T$-odd structure functions for the
$q\bar{q}$ annihilation and $qg$ Compton channels at order ${\cal
O}(\alpha_{s}^{2})$, where they become nonvanishing thanks to absorptive
contributions to loop amplitudes. While the corresponding results are not new,
we have used them in novel ways. Foremost, we have presented a new formalism
to expand the results for low transverse momentum, or low $\rho=Q_{T}/Q$, with
the goal of facilitating comparisons to frameworks that analyze $T$-odd
effects in terms of TMDs, especially at non-leading power. Our new formalism
is completely general and can in principle be used to obtain expansions to
arbitrary order in $Q_{T}/Q$. As a proof of concept, we have applied it to the
$T$-odd structure functions and expanded them to order ${\cal O}(\rho^{4})$.
In doing so, we uncovered a new relation between two of the $T$-odd structure
functions, $W^{\text{LP}}_{\nabla}(x_{1},x_{2})$ and
$W^{\text{LP}}_{\Delta_{P}}(x_{1},x_{2})$, valid at leading power in the
small-$Q_{T}$ expansion. Although in the present paper we have not attempted
to connect our results to calculations based on TMD factorization, we think
that our paper has much to offer for such comparisons in the future.

We have also presented numerical results for the validity of the expansion in
$\rho^{2}$, and we have compared our full results for the $T$-odd structure
functions to available data from the ATLAS experiment. We found that in the
present situation it is impossible to draw any quantitative conclusions from
this comparison.

In the present paper we have restricted our analysis to the case of the
$T$-odd effects in the Drell-Yan process with unpolarized beams. Extensions to
the $T$-even sector and to polarized scattering will be natural extensions of
our work.

###### Acknowledgements.

We are grateful to Marc Schlegel for helpful discussions. A.S.Z. thanks
Tübingen University for warm hospitality during his visit during which the
main part of this paper was completed. This work was funded by BMBF (Germany)
“Verbundprojekt 05P2021 (ErUM-FSP T01) – Run 3 von ALICE am LHC: Perturbative
Berechnungen von Wirkungsquerschnitten für ALICE” (Förderkennzeichen:
05P21VTCAA), by ANID PIA/APOYO AFB220004 (Chile), by FONDECYT (Chile) under
Grants No. 1230160 and No. 1240066, and by ANID$-$Millennium
Program$-$ICN2019_044 (Chile).

##  Appendix A Relations among different sets of the structure functions

The three sets of structure functions $\\{A_{i}\\}$, $\\{W_{i}\\}$ and
$\\{\lambda,\mu,\nu,\ldots\\}$ are related as Lam:1978pu ; Collins:1977iv ;
Mirkes:1992hu ; Boer:2006eq ; Berger:2007jw

$\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\frac{W_{T}-W_{L}}{W_{T}+W_{L}}=\frac{2-3A_{0}}{2+A_{0}}\,,\quad\mu\,=\,\frac{W_{\Delta}}{W_{T}+W_{L}}=\frac{2A_{1}}{2+A_{0}}\,,\quad\nu\,=\,\frac{2W_{\Delta\Delta}}{W_{T}+W_{L}}=\frac{2A_{2}}{2+A_{0}}\,,$
$\displaystyle\tau$ $\displaystyle=$
$\displaystyle\frac{W_{\nabla_{P}}}{W_{T}+W_{L}}=\frac{2A_{3}}{2+A_{0}}\,,\quad\hskip
5.69046pt\eta\,=\,\frac{W_{T_{P}}}{W_{T}+W_{L}}=\frac{2A_{4}}{2+A_{0}}\,,\quad\hskip
2.84544pt\xi\,=\,\frac{W_{\Delta\Delta_{P}}}{W_{T}+W_{L}}=\frac{2A_{5}}{2+A_{0}}\,,$
$\displaystyle\zeta$ $\displaystyle=$
$\displaystyle\frac{W_{\Delta_{P}}}{W_{T}+W_{L}}=\frac{2A_{6}}{2+A_{0}}\,,\quad\hskip
4.26773pt\chi\,=\,\frac{W_{\nabla}}{W_{T}+W_{L}}=\frac{2A_{7}}{2+A_{0}}\,,$
(60)

and

$\displaystyle A_{0}$ $\displaystyle=$
$\displaystyle\frac{2W_{L}}{2W_{T}+W_{L}}=\frac{2(1-\lambda)}{3+\lambda}\,,\quad
A_{1}\,=\,\frac{2W_{\Delta}}{2W_{T}+W_{L}}=\frac{4\mu}{3+\lambda}\,,\quad
A_{2}\,=\,\frac{4W_{\Delta\Delta}}{2W_{T}+W_{L}}=\frac{4\nu}{3+\lambda}\,,$
$\displaystyle A_{3}$ $\displaystyle=$
$\displaystyle\frac{2W_{\nabla_{P}}}{2W_{T}+W_{L}}=\frac{4\tau}{3+\lambda}\,,\quad\hskip
11.38092ptA_{4}\,=\,\frac{2W_{T_{P}}}{2W_{T}+W_{L}}=\frac{4\eta}{3+\lambda}\,,\quad\hskip
1.42271ptA_{5}\,=\,\frac{2W_{\Delta\Delta_{P}}}{2W_{T}+W_{L}}=\frac{4\xi}{3+\lambda}\,,$
$\displaystyle A_{6}$ $\displaystyle=$
$\displaystyle\frac{2W_{\Delta_{P}}}{2W_{T}+W_{L}}=\frac{4\zeta}{3+\lambda}\,,\quad\hskip
11.38092ptA_{7}\,=\,\frac{2W_{\nabla}}{2W_{T}+W_{L}}=\frac{4\chi}{3+\lambda}\,.$
(61)

##  Appendix B Treatment of $\gamma^{5}$ matrix in calculation of structure
functions

In the calculation of the $T$-odd structure functions we have to deal with the
$\gamma^{5}$ matrix in dimensional regularization. We encounter two different
cases:

  * (a)

Dirac traces with an even number (in practice, two) of $\gamma^{5}$ matrices
in case of $w_{\nabla}$;

  * (b)

Dirac traces with an odd number (in practice, one) of $\gamma^{5}$ matrices in
case of $w_{\Delta\Delta_{P}}$ and $w_{\Delta_{P}}$.

For case (a), it is permissible even in dimensional regularization to anti-
commute the two $\gamma^{5}$ matrices towards each other and to use
$(\gamma^{5})^{2}=1$. As a result, the contributions of [vector $\otimes$
vector] and [axial-vector $\otimes$ axial-vector] couplings to $w_{\nabla}$
are identical.

In case (b) we use techniques established in the literature for the treatment
of $\gamma^{5}$ in dimensional regularization. In particular, for the axial-
vector spin matrix we use the Larin prescription Akyeampong:1973xi ;
Larin:1991tj ; Larin:1993tq ; Zijlstra:1992kj , expressing $\gamma^{5}$ as the
product of the four-dimensional Levi-Civita tensor
$\epsilon^{\mu\nu\alpha\beta}$ and three gamma matrices:

$\displaystyle\gamma^{\mu}\gamma^{5}=\frac{i}{6}\,\epsilon^{\mu\nu\alpha\beta}\,\gamma_{\nu}\,\gamma_{\alpha}\,\gamma_{\beta}\,.$
(62)

Next, in order to evaluate the structure functions $w_{\Delta\Delta_{P}}$ and
$w_{\Delta_{P}}$, we apply the method proposed in Ref. Zijlstra:1992kj  for
the contraction of two Levi-Civita tensors. One Levi-Civita tensor occurs in
the definition of the structure functions $w_{\Delta\Delta_{P}}$ and
$w_{\Delta_{P}}$ [see Eqs. (12) and (13)] and the other appears because of the
Larin substitution (62). In general, the product of two Levi-Civita tensors
must be evaluated in terms of the $D$-dimensional Kronecker tensor
$\delta^{\mu}_{\nu}$ in order to preserve Lorentz invariance Zijlstra:1992kj :

$\displaystyle\epsilon^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}\,\epsilon_{\nu_{1}\nu_{2}\nu_{3}\nu_{4}}=-{\rm
det}\left|\begin{array}[]{cccc}\delta^{\mu_{1}}_{\nu_{1}}&\delta^{\mu_{1}}_{\nu_{2}}&\delta^{\mu_{1}}_{\nu_{3}}&\delta^{\mu_{1}}_{\nu_{4}}\\\\[8.53581pt]
\delta^{\mu_{2}}_{\nu_{1}}&\delta^{\mu_{2}}_{\nu_{2}}&\delta^{\mu_{2}}_{\nu_{3}}&\delta^{\mu_{2}}_{\nu_{4}}\\\\[8.53581pt]
\delta^{\mu_{3}}_{\nu_{1}}&\delta^{\mu_{3}}_{\nu_{2}}&\delta^{\mu_{3}}_{\nu_{3}}&\delta^{\mu_{3}}_{\nu_{4}}\\\\[8.53581pt]
\delta^{\mu_{4}}_{\nu_{1}}&\delta^{\mu_{4}}_{\nu_{2}}&\delta^{\mu_{4}}_{\nu_{3}}&\delta^{\mu_{4}}_{\nu_{4}}\\\
\end{array}\right|\,.$ (67)

For our purposes we need to evaluate two types of contractions [see Eqs. (12)
and (13)],

$\displaystyle i\epsilon^{\mu
PRK}\,\gamma_{\mu}\gamma_{5}=-\frac{1}{6}\,\epsilon^{\mu
PRK}\,\epsilon_{\mu\rho\alpha\beta}\,\gamma^{\rho}\,\gamma^{\alpha}\,\gamma^{\beta}\,,$
(68)

and

$\displaystyle i\epsilon^{\mu
PRK}\,\gamma_{\nu}\gamma_{5}=-\frac{1}{6}\,\epsilon^{\mu
PRK}\,\epsilon_{\nu\rho\alpha\beta}\,\gamma^{\rho}\,\gamma^{\alpha}\,\gamma^{\beta}\,.$
(69)

Using Eq. (67) we get

$\displaystyle i\epsilon^{\mu
PRK}\,\gamma_{\mu}\gamma_{5}=(D-3)\,\not\\!P\not\\!R\not\\!K\,,$ (70)

and

$\displaystyle i\epsilon^{\mu
PRK}\,\gamma_{\nu}\gamma_{5}=\delta^{\mu}_{\nu;\perp}\,\not\\!P\not\\!R\not\\!K+\gamma^{\mu}_{\perp}\,\Big{(}P_{\nu}\not\\!K\not\\!R+R_{\nu}\not\\!P\not\\!K+K_{\nu}\not\\!R\not\\!P\Big{)}\,,$
(71)

where $D=4-2\epsilon$ and $\delta^{\mu}_{\nu;\perp}$ and
$\gamma^{\mu}_{\perp}$ are the perpendicular $D$-dimensional Kronecker tensor
and gamma matrices, respectively, defined as

$\displaystyle\delta^{\mu}_{\nu;\perp}$ $\displaystyle\equiv$
$\displaystyle\delta^{\mu}_{\nu}-\frac{P^{\mu}P_{\nu}}{P^{2}}-\frac{R^{\mu}R_{\nu}}{R^{2}}-\frac{K^{\mu}K_{\nu}}{K^{2}}\,,$
$\displaystyle\gamma^{\mu}_{\perp}$ $\displaystyle\equiv$
$\displaystyle\delta^{\mu}_{\nu;\perp}\,\gamma^{\nu}=\gamma^{\mu}-\frac{P^{\mu}\\!\\!\not\\!P}{P^{2}}-\frac{R^{\mu}\\!\\!\not\\!R}{R^{2}}-\frac{K^{\mu}\\!\\!\not\\!K}{K^{2}}\,,$
(72)

which are manifestly orthogonal to all momenta of the basis $(P,R,K)$. The
obey the following conditions:

$\displaystyle\delta^{\mu}_{\nu;\perp}\,\delta^{\nu}_{\mu;\perp}=\delta^{\mu}_{\nu;\perp}\,\delta^{\nu}_{\mu}=\delta^{\mu}_{\mu;\perp}=D-3\,,$
$\displaystyle\delta^{\mu}_{\nu;\perp}\,P_{\mu}=\delta^{\mu}_{\nu;\perp}\,R_{\mu}=\delta^{\mu}_{\nu;\perp}\,K_{\mu}=0\,.$
(73)

Similarly one can define the perpendicular $D$-dimensional metric tensor
$g^{\mu\nu}_{\perp}$ introduced in Ref. Lyubovitskij:2021ges ,

$\displaystyle
g^{\mu\nu}_{\perp}=g^{\mu\nu}-\frac{P^{\mu}P^{\nu}}{P^{2}}-\frac{R^{\mu}R^{\nu}}{R^{2}}-\frac{K^{\mu}K^{\nu}}{K^{2}}\,.$
(74)

One should mention that the identities (69) and (71) are generalizations of
(68) and (70), respectively. In particular, Eq. (70) follows from (71) for
$\nu\to\mu$. In this limit, $\delta^{\mu}_{\mu;\perp}=D-3$, the second term on
the r.h.s. of Eq. (71) vanishes, and we arrive at the identity (70).

We also stress that using the orthogonal basis $(P,R,K)$ along with the
orthogonal metric tensors and gamma matrices turns out to be very useful and
economical in our analytical calculations. In addition, we note that the use
of identity (71) is further simplified in case of the evaluation of
$w_{\Delta\Delta_{P}}$ and $w_{\Delta_{P}}$. The reason is that the Levi-
Civita tensor $\epsilon^{\mu PRK}$ on the l.h.s. of Eq. (71) is accompanied by
the basis vector $P^{\nu}$ or $R^{\nu}$ [see Eqs. (12) and (13)]. Therefore,
the first Lorentz structure
$\delta^{\mu}_{\nu;\perp}\not\\!\\!P\not\\!\\!R\not\\!\\!K$ on the r.h.s. of
Eq. (71) vanishes thanks to
$P^{\nu}\,\delta^{\mu}_{\nu;\perp}=R^{\nu}\,\delta^{\mu}_{\nu;\perp}=0$. The
second Lorentz structure
$\gamma^{\mu}_{\perp}\,\Big{(}P_{\nu}\not\\!\\!K\not\\!\\!R+R_{\nu}\not\\!\\!P\not\\!\\!K+K_{\nu}\not\\!\\!R\not\\!\\!P\Big{)}$
on the r.h.s. of the identity (71) is also simplified because of the
orthogonality of the $(P,R,K)$ basis. In particular, depending on the
accompanying momentum $P^{\nu}$ or $R^{\nu}$, we deduce from the identity (71)
two simplified identities useful for the calculation of $w_{\Delta\Delta_{P}}$
and $w_{\Delta_{P}}$:

$\displaystyle i\epsilon^{\mu PRK}\,\not\\!P\gamma_{5}$ $\displaystyle=$
$\displaystyle\gamma^{\mu}_{\perp}\not\\!K\not\\!R\ P^{2}\,,$ $\displaystyle
i\epsilon^{\mu PRK}\,\not\\!R\gamma_{5}$ $\displaystyle=$
$\displaystyle\gamma^{\mu}_{\perp}\not\\!P\not\\!K\ R^{2}\,.$ (75)

##  Appendix C Calculational technique for the partial derivatives of the
structure functions $W_{i}(x_{1},x_{2})$

In this Appendix we discuss the calculational technique for the partial
derivatives of the structure functions $W_{i}(x_{1},x_{2})$ defined in Eq.
(46), which have the following form

$\displaystyle W_{i}(x_{1},x_{2},L_{\rho})$ $\displaystyle=$
$\displaystyle\frac{1}{x_{1}x_{2}}\sum_{a,b}\,\biggl{[}R_{ab,i}(x_{1},x_{2},L_{\rho})\,f_{a/H_{1}}(x_{1})\,f_{b/H_{2}}(x_{2})$
(76) $\displaystyle+$ $\displaystyle\Big{(}P_{ba,i}\otimes
f_{b/H_{2}}\Big{)}(x_{2},x_{1},L_{\rho})\
f_{a/H_{1}}(x_{1})+\Big{(}P_{ab,i}\otimes
f_{a/H_{1}}\Big{)}(x_{1},x_{2},L_{\rho})\ f_{b/H_{2}}(x_{2})\biggr{]}\,,$

where $R_{ab,i}(x_{1},x_{2},L_{\rho})$, $P_{ab,i}(z_{1},x_{2},L_{\rho})$, and
$P_{ba,i}(z_{2},x_{1},L_{\rho})$ are perturbative coefficient functions
containing differential operators acting on the PDFs $f_{a/H_{1}}(x_{1})$ and
$f_{b/H_{2}}(x_{2})$. In particular, the functions
$R_{ab,i}(x_{1},x_{2},L_{\rho})$, $P_{ab,i}(z_{1},x_{2},L_{\rho})$, and
$P_{ba,i}(z_{2},x_{1},L_{\rho})$ can be expanded using the set of the
differential operators to factorize the dependence on the variables $x_{1}$
and $x_{2}$, $z_{1}$ and $x_{2}$, $z_{2}$ and $x_{1}$, respectively (see
detailed discussion in the supplemental material). One can see that
$W_{i}(x_{1},x_{2},L_{\rho})$ is composed of three main terms, one term
proportional to the perturbative function $R_{i}(x_{1},x_{2},L_{\rho})$ times
the PDFs, and two terms each containing a single integral representation of
the convolution of the perturbative functions $P_{a,i}(z_{1},L_{\rho})$, and
$P_{b,i}(z_{2},L_{\rho})$ with the PDFs. Also the latter terms depends on the
second variable $x_{i}=x_{1}$ or $x_{2}$ via the simple form of a product of
$f(x_{i})$ with some polynomial in $x_{i}$.

The first term can be trivially differentiated with respect to $x_{1}$ and
$x_{2}$ to the desired order. In case of the other two terms one can
straightforwardly differentiate the non-integral terms. It remains to specify
how to take the partial derivative of the terms with integrals. Let us discuss
the treatment of such terms by considering the following integral:

$\displaystyle J(x)$ $\displaystyle=$ $\displaystyle\frac{I(x)}{x}\,,\qquad
I(x)\,=\,(P\otimes
f)(x)=\int\limits_{x}^{1}\frac{dz}{z}\,P(z)\,f\Big{(}\frac{x}{z}\Big{)}\,.$
(77)

The $n$-th derivative of the integral $J(x)$ can be taken using the binomial
formula

$\displaystyle\frac{\partial^{n}J(x)}{\partial x^{n}}$ $\displaystyle=$
$\displaystyle\sum\limits_{m=0}^{n}\,\frac{n!}{m!}\,\frac{(-1)^{n-m}}{x^{n-m+1}}\,\frac{\partial^{m}I(x)}{\partial
x^{m}}\,.$ (78)

Therefore, we only need to derive an analytical formula for the $n$-th
derivative of the integral $I(x)$, where $n$ is an arbitrary natural number.
In our derivation we will consider two possible choices for the function
$P(z)$:

  * (a)

$P(z)=R(z)$ is a regular function of the variable $z$,

  * (b)

$P(z)=[1/(1-z)^{m}]_{+,m-1}$ is the generalized plus distribution of power
$m$, defined in Eq. (40) in the main text.

We first consider the simpler case (a). Here, the first-order derivative of
$I(x)$ reads

$\displaystyle\frac{\partial I(x)}{\partial x}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
x}\,\int\limits_{x}^{1}\frac{dz}{z}\,R(z)\,f\Big{(}\frac{x}{z}\Big{)}=\frac{\partial}{\partial
x}\,\int\limits_{x}^{1}\frac{d\xi}{\xi}\,R(x/\xi)\,f(\xi)$ (79)
$\displaystyle=$
$\displaystyle-\int\limits_{x}^{1}\frac{d\xi}{\xi}\,\delta(\xi-x)\,R(x/\xi)\,f(\xi)+\int\limits_{x}^{1}\frac{d\xi}{\xi}\,\frac{\partial
R(x/\xi)}{\partial x}\,f(\xi)$ $\displaystyle=$
$\displaystyle-\int\limits_{x}^{1}\frac{d\xi}{\xi^{2}}\,\delta(1-x/\xi)\,R(x/\xi)\,f(\xi)+\int\limits_{x}^{1}\frac{d\xi}{\xi^{2}}\,\frac{\partial
R(x/\xi)}{\partial(x/\xi)}\,f(\xi)$ $\displaystyle=$
$\displaystyle-\frac{1}{x}\int\limits_{x}^{1}dz\,\delta(1-z)\,R(z)\,f\Big{(}\frac{x}{z}\Big{)}+\frac{1}{x}\int\limits_{x}^{1}dz\,\frac{\partial
R(z)}{\partial z}\,f\Big{(}\frac{x}{z}\Big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{x}\biggl{[}-R(1)\,f(x)+\int\limits_{x}^{1}dz\,R^{\prime}(z)\,f\Big{(}\frac{x}{z}\Big{)}\biggr{]}\,,$

where $R^{\prime}(z)=\partial R(z)/\partial z$.

One can prove by induction that the $n$-th derivative of $I(x)$ is given by

$\displaystyle\frac{\partial^{n}I(x)}{\partial
x^{n}}=-\sum\limits_{k=1}^{n}\,\biggl{(}\frac{f(x)}{x^{k}}\biggr{)}^{(n-k)}_{x}\,R^{(k-1)}(1)+\frac{1}{x^{n}}\,\int\limits_{x}^{1}dz\,z^{n-1}\,R^{(n)}(z)\,f\Big{(}\frac{x}{z}\Big{)}\,,$
(80)

where

$\displaystyle R^{(k)}(z)=\partial^{k}R(z)/\partial
z^{k}\,,\qquad(\cdots)^{(k)}_{x}=\partial^{k}(\cdots)/\partial x^{k}\,.$ (81)

For case (b) we recall the definition for the generalized plus distribution
when applied to a PDF $f(x/z)$:

$\displaystyle\int\limits_{x}^{1}dz\,\frac{f\Big{(}\frac{x}{z}\Big{)}}{(1-z)^{m}_{+,m-1}}$
$\displaystyle=$
$\displaystyle\int\limits_{x}^{1}dz\,\biggl{[}\frac{1}{(1-z)^{m}_{+x,m-1}}+\delta(1-z)\,\log(1-x)\,\frac{(-1)^{m-1}}{(m-1)!}\,\partial_{z}^{m-1}$
(82) $\displaystyle-$
$\displaystyle\delta(1-z)\,\sum\limits_{j=2}^{m}\,\frac{(-1)^{m-j}}{(j-1)\,(m-j)!}\,\biggl{(}\frac{1}{(1-x)^{j-1}}-1\biggr{)}\,\partial_{z}^{m-j}\biggr{]}f\Big{(}\frac{x}{z}\Big{)}\,.$

The $x$-dependence of $f(x/z)/(1-z)^{m}_{+x,m-1}$ induces a subtraction of the
$(m-1)$-th order Taylor polynomial $T[f(x/z)]^{m-1}_{z=1}$ evaluated at $z=1$,

$\displaystyle\frac{f\Big{(}\frac{x}{z}\Big{)}}{(1-z)^{m}_{+x,m-1}}$
$\displaystyle=$
$\displaystyle\frac{f\Big{(}\frac{x}{z}\Big{)}-T\Big{[}f\Big{(}\frac{x}{z}\Big{)}\Big{]}^{m-1}_{z=1}}{(1-z)^{m}}\,,$
$\displaystyle T\Big{[}f\Big{(}\frac{x}{z}\Big{)}\Big{]}^{m-1}_{z=1}$
$\displaystyle=$
$\displaystyle\sum\limits_{k=0}^{m-1}\frac{(z-1)^{k}}{k!}\,\partial_{z}^{k}f\Big{(}\frac{x}{z}\Big{)}\Big{|}_{z=1}$
(83) $\displaystyle=$ $\displaystyle
f(x)+\sum\limits_{k=1}^{m-1}\sum\limits_{\ell=1}^{k}\,C_{k-1}^{\ell-1}\,(1-z)^{k}\,\frac{x^{\ell}}{\ell!}\,f^{(\ell)}(x)\,,$

where $C_{m}^{\ell}=m!/(\ell!(m-\ell)!)$ is the binomial coefficient and
$f^{(\ell)}(x)=\partial^{\ell}f(x)/\partial x^{\ell}$. Note that the partial
derivatives $\delta(1-z)\partial_{z}^{m-j}f(x/z)$ with respect to $z$ in Eq.
(82) can be simplified and reduced to derivatives with respect to $x$ as

$\displaystyle\delta(1-z)\,\partial_{z}^{m}f\Big{(}\frac{x}{z}\Big{)}=\delta(1-z)\,\sum\limits_{\ell=1}^{m}(-1)^{m}C_{m}^{\ell}\,\frac{(m-1)!}{(\ell-1)!}\,x^{\ell}\,f^{(\ell)}(x)\,.$
(84)

In all cases the derivatives $\partial^{n}I(x)/\partial x^{n}$ can be easily
taken by changing the integration variable $z\to x/\xi$ and, after some
simplifications, returning back to the integration over $z$. We stress that
the choice $P(z)=R(z)[1/(1-z)^{m}]_{+,m-1}$, where $R(z)$ is a regular
function, can be reduced to case (b) by carrying out a Taylor expansion of
$R(z)$ around $z=1$:

$\displaystyle
R(z)=\sum\limits_{k=0}^{\infty}\,\frac{(z-1)^{k}}{k!}\,R^{(k)}(1)\,$ (85)

where $R^{(k)}(1)=(\partial^{k}R(z)/\partial z^{k})_{z=1}$, which results in
cancellation of the respective powers of $(1-z)$ between the numerator and the
denominator of the integrand of $I(x)$. In particular, the distribution
$P_{qq}(z)$ defined in Eq. (IV) contains the term $(1+z^{2})/(1-z)_{+}$, which
can be represented as the sum of a regular term corresponding to case (a) and
a single distribution $2/(1-z)_{+}$ corresponding to case (b),

$\displaystyle\frac{1+z^{2}}{(1-z)_{+}}=-(1+z)+\frac{2}{(1-z)_{+}}\,.$ (86)

In case (b) the $n$-th derivative of the integral $I(x)$ reads

$\displaystyle\frac{\partial^{n}I(x)}{\partial x^{n}}$ $\displaystyle=$
$\displaystyle-\lim\limits_{z\to
1}\,\sum\limits_{k=1}^{n}\,\left(\frac{1}{x^{k}}\,\frac{\partial^{k-1}}{\partial
z^{k-1}}\,\biggl{[}\frac{1}{(1-z)^{m}_{+x,m-1}}\biggr{]}\,z^{k-2}f\Big{(}\frac{x}{z}\Big{)}\right)^{(n-k)}_{x}$
(87) $\displaystyle+$
$\displaystyle\frac{1}{x^{n}}\,\int\limits_{x}^{1}dz\,\frac{\partial^{n}}{\partial
z^{n}}\,\biggl{[}\frac{1}{(1-z)^{m}_{+,m-1}}\biggr{]}\,z^{n-1}f\Big{(}\frac{x}{z}\Big{)}$
$\displaystyle=$ $\displaystyle-\lim\limits_{z\to
1}\,\sum\limits_{k=1}^{n}\,\left(\frac{(m+k-2)!}{x^{k}\,(m-1)!}\,\frac{z^{k-2}f\Big{(}\frac{x}{z}\Big{)}}{(1-z)^{m+k-1}_{+x,m+k-2}}\right)^{(n-k)}_{x}$
$\displaystyle+$
$\displaystyle\frac{(m+n-1)!}{x^{n}\,(m-1)!}\,\int\limits_{x}^{1}dz\,\,\frac{z^{n-1}f\Big{(}\frac{x}{z}\Big{)}}{(1-z)^{m+n}_{+,m+n-1}}$
$\displaystyle=$ $\displaystyle\lim\limits_{z\to
1}\,\sum\limits_{k=1}^{n}\,\left(\frac{(-1)^{m+k}}{x^{k}\,(m-1)!\,(m+k-1)}\,\frac{\partial^{m+k-1}}{\partial
z^{m+k-1}}\,\biggl{[}z^{k-2}\,f\Big{(}\frac{x}{z}\Big{)}\biggr{]}\right)^{(n-k)}_{x}$
$\displaystyle+$
$\displaystyle\frac{(m+n-1)!}{x^{n}\,(m-1)!}\,\int\limits_{x}^{1}dz\,\,\frac{z^{n-1}f\Big{(}\frac{x}{z}\Big{)}}{(1-z)^{m+n}_{+,m+n-1}}\,.$

## References

  * (1) A. De Rujula, J. M. Kaplan, and E. De Rafael, Nucl. Phys. B 35, 365 (1971).
  * (2) A. De Rujula, R. Petronzio, and B. E. Lautrup, Nucl. Phys. B 146, 50 (1978).
  * (3) K. Fabricius, I. Schmitt, G. Kramer, and G. Schierholz, Phys. Rev. Lett. 45, 867 (1980).
  * (4) J. G. Körner, G. Kramer, G. Schierholz, K. Fabricius, and I. Schmitt, Phys. Lett. B 94, 207 (1980).
  * (5) K. Hagiwara, K. I. Hikasa, and N. Kai, Phys. Rev. Lett. 47, 983 (1981).
  * (6) K. Hagiwara, K. I. Hikasa, and N. Kai, Phys. Rev. D 27, 84 (1983).
  * (7) B. Pire and J. P. Ralston, Phys. Rev. D 28, 260 (1983).
  * (8) K. Hagiwara, K. I. Hikasa, and N. Kai, Phys. Rev. Lett. 52, 1076 (1984).
  * (9) A. Bilal, E. Masso, and A. De Rujula, Nucl. Phys. B 355, 549 (1991).
  * (10) R. D. Carlitz and R. S. Willey, Phys. Rev. D 45, 2323 (1992).
  * (11) E. Mirkes, Nucl. Phys. B 387, 3 (1992).
  * (12) A. Brandenburg, L. J. Dixon, and Y. Shadmi, Phys. Rev. D 53, 1264 (1996) [arXiv:hep-ph/9505355 [hep-ph]].
  * (13) M. Ahmed and T. Gehrmann, Phys. Lett. B 465, 297 (1999) [arXiv:hep-ph/9906503 [hep-ph]].
  * (14) J. G. Körner, B. Melic, and Z. Merebashvili, Phys. Rev. D 62, 096011 (2000) [arXiv:hep-ph/0002302 [hep-ph]].
  * (15) H. Yokoya, Prog. Theor. Phys. 118, 371 (2007) [arXiv:0705.2481 [hep-ph]].
  * (16) K. Hagiwara, K. Mawatari, and H. Yokoya, JHEP 12, 041 (2007) [arXiv:0707.3194 [hep-ph]].
  * (17) R. Frederix, K. Hagiwara, T. Yamada, and H. Yokoya, Phys. Rev. Lett. 113, 152001 (2014) [arXiv:1407.1016 [hep-ph]].
  * (18) A. Gehrmann-De Ridder, T. Gehrmann, E. W. N. Glover, A. Huss, and T. A. Morgan, JHEP 07, 133 (2016) [arXiv:1605.04295 [hep-ph]].
  * (19) S. Benic, Y. Hatta, H. N. Li, and D. J. Yang, Phys. Rev. D 100, 094027 (2019) [arXiv:1909.10684 [hep-ph]].
  * (20) S. Benić, Y. Hatta, A. Kaushik, and H. N. Li, Phys. Rev. D 104, 094027 (2021) [arXiv:2109.05440 [hep-ph]].
  * (21) S. Benić, Y. Hatta, A. Kaushik, and H. N. Li, [arXiv:2402.02267 [hep-ph]].
  * (22) M. Abele, M. Aicher, F. Piacenza, A. Schäfer, and W. Vogelsang, Phys. Rev. D 106, 014020 (2022) [arXiv:2204.13967 [hep-ph]].
  * (23) D. W. Sivers, Phys. Rev. D 41, 83 (1990).
  * (24) D. Boer and P. J. Mulders, Phys. Rev. D 57, 5780 (1998) [arXiv:hep-ph/9711485 [hep-ph]].
  * (25) S. J. Brodsky, D. S. Hwang, and I. Schmidt, Phys. Lett. B 530, 99 (2002) [arXiv:hep-ph/0201296 [hep-ph]].
  * (26) J. C. Collins, Phys. Lett. B 536, 43 (2002) [arXiv:hep-ph/0204004 [hep-ph]].
  * (27) R. Boussarie, M. Burkardt, M. Constantinou, W. Detmold, M. Ebert, M. Engelhardt, S. Fleming, L. Gamberg, X. Ji, and Z. B. Kang, et al. [arXiv:2304.03302 [hep-ph]].
  * (28) X. Ji, J. W. Qiu, W. Vogelsang, and F. Yuan, Phys. Rev. Lett. 97, 082002 (2006) [arXiv:hep-ph/0602239 [hep-ph]].
  * (29) X. Ji, J. W. Qiu, W. Vogelsang, and F. Yuan, Phys. Rev. D 73, 094017 (2006) [arXiv:hep-ph/0604023 [hep-ph]].
  * (30) D. Gutiérrez-Reyes, I. Scimemi, and A. A. Vladimirov, Phys. Lett. B 769, 84 (2017) [arXiv:1702.06558 [hep-ph]].
  * (31) D. Gutierrez-Reyes, I. Scimemi, and A. Vladimirov, JHEP 07, 172 (2018) [arXiv:1805.07243 [hep-ph]].
  * (32) I. Scimemi, A. Tarasov, and A. Vladimirov, JHEP 05, 125 (2019) [arXiv:1901.04519 [hep-ph]].
  * (33) F. Rein, S. Rodini, A. Schäfer, and A. Vladimirov, JHEP 01, 116 (2023) [arXiv:2209.00962 [hep-ph]].
  * (34) A. Bacchetta, D. Boer, M. Diehl, and P. J. Mulders, JHEP 08, 023 (2008) [arXiv:0803.0227 [hep-ph]].
  * (35) D. Boer and W. Vogelsang, Phys. Rev. D 74, 014004 (2006) [hep-ph/0604177].
  * (36) S. Rodini and A. Vladimirov, JHEP 09, 117 (2023) [arXiv:2211.04494 [hep-ph]].
  * (37) L. Gamberg, Z. B. Kang, D. Y. Shao, J. Terry, and F. Zhao, [arXiv:2211.13209 [hep-ph]].
  * (38) M. A. Ebert, A. Gao, and I. W. Stewart, JHEP 06, 007 (2022) [erratum: JHEP 07, 096 (2023)] [arXiv:2112.07680 [hep-ph]].
  * (39) S. Rodini and A. Vladimirov, [arXiv:2306.09495 [hep-ph]].
  * (40) G. Aad et al. (ATLAS Collaboration), JHEP 08, 159 (2016) [arXiv:1606.00689 [hep-ex]].
  * (41) D. Boer, Phys. Rev. D 60, 014012 (1999) [arXiv:hep-ph/9902255 [hep-ph]].
  * (42) Z. Lu and I. Schmidt, Phys. Rev. D 84, 094002 (2011) [arXiv:1107.4693 [hep-ph]].
  * (43) T. Liu and B. Q. Ma, Eur. Phys. J. C 73, 2291 (2013) [arXiv:1201.2472 [hep-ph]].
  * (44) M. A. Ebert, J. K. L. Michel, I. W. Stewart, and F. J. Tackmann, JHEP 04, 102 (2021) [arXiv:2006.11382 [hep-ph]].
  * (45) E. L. Berger, J. -W. Qiu, and R. A. Rodriguez-Pedraza, Phys. Rev. D 76, 074006 (2007) [arXiv:0708.0578 [hep-ph]]; Phys. Lett. B 656, 74 (2007) [arXiv:0707.3150 [hep-ph]].
  * (46) C. S. Lam and W. -K. Tung, Phys. Rev. D 18, 2447 (1978).
  * (47) J. C. Collins and D. E. Soper, Phys. Rev. D 16, 2219 (1977).
  * (48) V. E. Lyubovitskij, F. Wunder, and A. S. Zhevlakov, JHEP 06, 066 (2021) [arXiv:2102.08943 [hep-ph]].
  * (49) J. Haug and F. Wunder, JHEP 02, 177 (2023) [arXiv:2211.14110 [hep-ph]].
  * (50) R. L. Workman et al. (Particle Data Group), PTEP 2022, 083C01 (2022).
  * (51) R. Meng, F. I. Olness, and D. E. Soper, Phys. Rev. D 54, 1919 (1996) [arXiv:hep-ph/9511311 [hep-ph]].
  * (52) I. M. Gelfand and G. E. Shilov, Generalized functions. Vol. 1. Properties and operations, Academic Press, New York, NY, 1964, 423 p.
  * (53) V. E. Lyubovitskij, W. Vogelsang, F. Wunder, and A. S. Zhevlakov, in preparation.
  * (54) D. Stump, J. Huston, J. Pumplin, W. K. Tung, H. L. Lai, S. Kuhlmann, and J. F. Owens, JHEP 10, 046 (2003) [arXiv:hep-ph/0303013 [hep-ph]].
  * (55) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, and G. Watt, Eur. Phys. J. C 75, 132 (2015) [arXiv:1412.7420 [hep-ph]].
  * (56) D. B. Clark, E. Godat, and F. I. Olness, Comput. Phys. Commun. 216, 126 (2017) [arXiv:1605.08012 [hep-ph]].
  * (57) M. Lambertsen and W. Vogelsang, Phys. Rev. D 93, 114013 (2016) [arXiv:1605.02625 [hep-ph]].
  * (58) V. Khachatryan et al. (CMS Collaboration), Phys. Lett. B 750, 154 (2015) [arXiv:1504.03512 [hep-ex]].
  * (59) S. Falciano et al. (NA10 Collaboration), Z. Phys. C 31, 513 (1986).
  * (60) R. Aaij et al. (LHCb Collaboration), Phys. Rev. Lett. 129, 091801 (2022) [arXiv:2203.01602 [hep-ex]].
  * (61) E. Richter-Was and Z. Was, Eur. Phys. J. C 77, 111 (2017) [arXiv:1609.02536 [hep-ph]].
  * (62) T. Aaltonen et al. (CDF Collaboration), Phys. Rev. Lett. 106, 241801 (2011) [arXiv:1103.5699 [hep-ex]].
  * (63) E. Mirkes, J. G. Korner, and G. A. Schuler, Phys. Lett. B 259, 151 (1991).
  * (64) G. Altarelli, R. K. Ellis, M. Greco, and G. Martinelli, Nucl. Phys. B 246, 12 (1984).
  * (65) S. Catani, L. Cieri, G. Ferrera, D. de Florian, and M. Grazzini, Phys. Rev. Lett. 103, 082001 (2009) [arXiv:0903.2120 [hep-ph]].
  * (66) D. A. Akyeampong and R. Delbourgo, Nuovo Cim. A 17, 578 (1973).
  * (67) S. A. Larin and J. A. M. Vermaseren, Phys. Lett. B 259, 345 (1991).
  * (68) S. A. Larin, Phys. Lett. B 303, 113 (1993) [arXiv:hep-ph/9302240 [hep-ph]].
  * (69) E. B. Zijlstra and W. L. van Neerven, Phys. Lett. B 297, 377 (1992).

## Supplementary Material: NLP and NNLP contributions to the $T$-odd hadronic
structure functions $W_{i}(x_{1},x_{2})$

Here we list the higher-order terms $W^{\rm
NLP}(x_{1}^{0},x_{2}^{0},L_{\rho})$ and $W^{\rm
NNLP}(x_{1}^{0},x_{2}^{0},,L_{\rho})$ in the small-$Q_{T}$ expansion,

$\displaystyle W^{\rm NLP}(x_{1}^{0},x_{2}^{0},L_{\rho})$ $\displaystyle=$
$\displaystyle
W_{1}(x_{1}^{0},x_{2}^{0},L_{\rho})+\frac{1}{2}\biggl{(}x_{1}^{0}\
\partial_{x_{1}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})+x_{2}^{0}\
\partial_{x_{2}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})\biggr{)}\,,$ (S1)
$\displaystyle W^{\rm NNLP}(x_{1}^{0},x_{2}^{0},L_{\rho})$ $\displaystyle=$
$\displaystyle
W_{2}(x_{1}^{0},x_{2}^{0},L_{\rho})+\frac{1}{4}x_{1}^{0}x_{2}^{0}\
\partial_{x_{1}^{0}}\partial_{x_{2}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})$
(S2) $\displaystyle-$ $\displaystyle\frac{1}{8}\biggl{(}x_{1}^{0}\
\partial_{x_{1}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})-4x_{1}^{0}\
\partial_{x_{1}^{0}}W_{1}(x_{1}^{0},x_{2}^{0},L_{\rho})-(x_{1}^{0})^{2}\
\partial^{2}_{x_{1}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})\biggr{)}$
$\displaystyle-$ $\displaystyle\frac{1}{8}\biggl{(}x_{2}^{0}\
\partial_{x_{2}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})-4x_{2}^{0}\
\partial_{x_{2}^{0}}W_{1}(x_{1}^{0},x_{2}^{0},L_{\rho})-(x_{2}^{0})^{2}\
\partial^{2}_{x_{2}^{0}}W_{0}(x_{1}^{0},x_{2}^{0},L_{\rho})\biggr{)}\,.$

Also we have abbreviated $L_{\rho}\equiv\log\rho^{2}$.

The higher-order terms $W_{1;J}^{ab}(x_{1},x_{2},L_{\rho})$ and
$W_{2;J}^{ab}(x_{1},x_{2},L_{\rho})$ are given by

$\displaystyle W^{ab}_{k;J}(x_{1},x_{2},L_{\rho})$ $\displaystyle=$
$\displaystyle-\frac{g_{ab}}{4x_{1}x_{2}}\biggl{[}R_{ab,k}^{J}(x_{1},x_{2},L_{\rho})\,f_{a/H_{1}}(x_{1})\,f_{b/H_{2}}(x_{2})$
(S3) $\displaystyle+$ $\displaystyle\Big{(}P_{ba,k}^{J}\otimes
f_{b/H_{2}}\Big{)}(x_{2},x_{1},L_{\rho})\
f_{a/H_{1}}(x_{1})+\Big{(}P_{ab,k}^{J}\otimes
f_{a/H_{1}}\Big{)}(x_{1},x_{2},L_{\rho})\ f_{b/H_{2}}(x_{2})\biggr{]}\,,$

where $k=1,2$ and $J=\Delta\Delta_{P},\Delta_{P},\nabla$.

In the following, for convenience, we introduce functions $F_{i}(z)$ with
$i=1,2,3,4,5,6$. Functions $F_{i}(z)$ with $i=3,4,5,6$ are obtained from
functions $F_{1}(z)$ and $F_{2}(z)$ by subtraction of the leading and sub-
leading terms in an expansion in $(1-z)$. The main idea of such substracted
functions is to make the ratios $F_{i}(z)/(1-z)^{n}$ regular at $z=1$. All
functions $F_{i}(z)$ with $i=1,\ldots,6$ are defined as

$\displaystyle F_{1}(z)$ $\displaystyle=$
$\displaystyle\frac{1+z}{1-z}+\frac{2z\log(z)}{(1-z)^{2}}=2\sum\limits_{N=1}^{\infty}\,\frac{(1-z)^{N}}{(N+1)(N+2)}={\cal
O}(1-z)\,,$ $\displaystyle F_{2}(z)$ $\displaystyle=$ $\displaystyle
1+\frac{z\log(z)}{1-z}=\frac{1-z}{2}\,\Big{(}1+F_{1}(z)\Big{)}=\sum\limits_{N=1}^{\infty}\,\frac{(1-z)^{N}}{N(N+1)}={\cal
O}(1-z)\,,$ $\displaystyle F_{3}(z)$ $\displaystyle=$ $\displaystyle
F_{1}(z)-\frac{1-z}{3}=2\sum\limits_{N=2}^{\infty}\,\frac{(1-z)^{N}}{(N+1)(N+2)}={\cal
O}\Big{(}(1-z)^{2}\Big{)}\,,$ $\displaystyle F_{4}(z)$ $\displaystyle=$
$\displaystyle
F_{1}(z)-\frac{1-z}{3}-\frac{(1-z)^{2}}{6}=2\sum\limits_{N=3}^{\infty}\,\frac{(1-z)^{N}}{(N+1)(N+2)}={\cal
O}\Big{(}(1-z)^{3}\Big{)}\,,$ $\displaystyle F_{5}(z)$ $\displaystyle=$
$\displaystyle
F_{2}(z)-\frac{1-z}{2}=\frac{1-z}{2}F_{1}(z)=\sum\limits_{N=2}^{\infty}\,\frac{(1-z)^{N}}{N(N+1)}={\cal
O}\Big{(}(1-z)^{2}\Big{)}\,,$ $\displaystyle F_{6}(z)$ $\displaystyle=$
$\displaystyle
F_{2}(z)-\frac{1-z}{2}-\frac{(1-z)^{2}}{6}=\frac{1-z}{2}\Big{(}F_{1}(z)-\frac{1-z}{3}\Big{)}=\frac{1-z}{2}F_{3}(z)$
(S4) $\displaystyle=$
$\displaystyle\sum\limits_{N=3}^{\infty}\,\frac{(1-z)^{N}}{N(N+1)}={\cal
O}\Big{(}(1-z)^{3}\Big{)}\,.$

They obey the following normalization conditions:

$\displaystyle F_{i}(1)=0\,,\quad F_{1}(0)=F_{2}(0)=1\,,\quad
F_{3}(0)=\frac{2}{3}\,,\quad F_{4}(0)=F_{5}(0)=\frac{1}{2}\,,\quad
F_{6}(0)=\frac{1}{3}\,.$ (S5)

It is clear that the ratios

$\displaystyle f_{i_{1}}(z)=\frac{F_{i_{1}}(z)}{1-z}\,,\ i_{1}=1,2;\qquad
f_{i_{2}}(z)=\frac{F_{i_{2}}(z)}{(1-z)^{2}}\,,\ i_{2}=3,5;\qquad
f_{i_{3}}(z)=\frac{F_{i_{3}}(z)}{(1-z)^{3}}\,,\ i_{3}=4,6\,,$ (S6)

are regular at $z=1$. In the following, we express the NLP and NNLP
contributions to the hadronic structure functions in terms of the functions
$f_{i}(z)$.

The perturbative functions $R_{ab,k}^{J}(x_{1},x_{2},L_{\rho})$,
$P_{ab,k}^{J}(z_{1},x_{2},L_{\rho})$, and $P_{ba,k}^{J}(z_{2},x_{1},L_{\rho})$
with $k=1,2$ can be expanded as

$\displaystyle R_{ab,k}^{J}(x_{1},x_{2},L_{\rho})$ $\displaystyle=$
$\displaystyle\sum\limits_{N=1}^{9}\,R_{ab;N}^{k;J}(L_{\rho})\,T^{R}_{N}(x_{1},x_{2})\,,$
$\displaystyle P_{ab,k}^{J}(z_{1},x_{2},L_{\rho})$ $\displaystyle=$
$\displaystyle\sum\limits_{N=1}^{3}\,P_{a,1;N}^{k;J}(z_{1},L_{\rho})\,T^{P}_{N}(x_{2})\,,$
$\displaystyle P_{ba,k}^{J}(z_{2},x_{1},L_{\rho})$ $\displaystyle=$
$\displaystyle\sum\limits_{N=1}^{3}\,P_{ab,2;N}^{k;J}(z_{2},L_{\rho})\,T^{P}_{N}(x_{1})\,,$
(S7)

where

$\displaystyle T^{R}_{1}(x_{1},x_{2})=1\,,\hskip
71.13188ptT^{R}_{2}(x_{1},x_{2})=x_{1}\,\partial_{x_{1}}\,,\hskip
45.52458ptT^{R}_{3}(x_{1},x_{2})=x_{2}\,\partial_{x_{2}}\,,$ $\displaystyle
T^{R}_{4}(x_{1},x_{2})=x_{1}x_{2}\,\partial_{x_{1}x_{2}}^{2}\,,\hskip
31.2982ptT^{R}_{5}(x_{1},x_{2})=x_{1}^{2}\,\partial_{x_{1}}^{2}\,,\hskip
45.52458ptT^{R}_{6}(x_{1},x_{2})=x_{2}^{2}\,\partial_{x_{2}}^{2}\,,$
$\displaystyle
T^{R}_{7}(x_{1},x_{2})=x_{1}^{2}x_{2}\,\partial_{x_{1}}^{2}\partial_{x_{2}}\,,\hskip
25.6073ptT^{R}_{8}(x_{1},x_{2})=x_{1}x_{2}^{2}\,\partial_{x_{1}}\partial_{x_{2}}^{2}\,,\hskip
21.33955ptT^{R}_{9}(x_{1},x_{2})=x_{1}^{2}x_{2}^{2}\,\partial_{x_{1}}^{2}\partial_{x_{2}}^{2}$
(S8)

and

$\displaystyle T^{P}_{1}(x)=1\,,\qquad T^{P}_{2}(x)=x\,\partial_{x}\,,\qquad
T^{P}_{3}(x)=x^{2}\,\partial_{x}^{2}\,.$ (S9)

Below we list only the nonvanishing terms $R_{ab;N}^{k;J}$,
$P_{ab,1;N}^{k;J}(z_{1})$, $P_{ab,2;N}^{k;J}(z_{2})$ for the various structure
functions.

For the $\Delta\Delta_{P}$ structure function we have, in the $q\bar{q}$
annihilation channel:

$\displaystyle R_{qq;1}^{1;\Delta\Delta_{P}}$ $\displaystyle=$ $\displaystyle
C_{A}+\Big{(}\frac{C_{A}}{2}+\frac{4C_{1}}{3}\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qq;2}^{1;\Delta\Delta_{P}}$ $\displaystyle=$ $\displaystyle
R_{qq;3}^{1;\Delta\Delta_{P}}\,=\,\Big{(}C_{A}-\frac{C_{1}}{3}\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qq;4}^{1;\Delta\Delta_{P}}$ $\displaystyle=$ $\displaystyle-
C_{A}\,L_{\rho}\,.$ (S10)

$\displaystyle P_{qq,1;1}^{1;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle
P_{qq,2;1}^{1;\Delta\Delta_{P}}(z)\,=\,-\frac{C_{A}}{4}\,\frac{1+3z+z^{2}-z^{3}}{(1-z)^{2}_{+,1}}-\frac{C_{1}}{6}\,\frac{1+2z-3z^{2}-2z^{3}}{(1-z)^{2}_{+,1}}$
$\displaystyle-$ $\displaystyle\frac{C_{1}}{2}\,f_{3}(z)(1+3z)\,,$
$\displaystyle P_{qq,1;2}^{1;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle
P_{qq,2;2}^{1;\Delta\Delta_{P}}(z)\,=\,\frac{C_{A}}{2}\,\frac{z(1+z^{2})}{(1-z)^{2}_{+,1}}+\frac{C_{1}}{3}\,\,\frac{z}{(1-z)_{+}}$
(S11) $\displaystyle+$ $\displaystyle C_{1}\,f_{3}(z)\,z\,.$

$\displaystyle R_{qq;1}^{2;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\Big{(}C_{A}+\frac{4C_{1}}{3}\Big{)}-\Big{(}\frac{3C_{A}}{8}+C_{1}\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qq;2}^{2;\Delta\Delta_{P}}$ $\displaystyle=$ $\displaystyle
R_{qq;3}^{2;\Delta\Delta_{P}}\,=\,-\Big{(}C_{A}-\frac{C_{1}}{3}\Big{)}-\Big{(}\frac{3C_{A}}{2}+\frac{5C_{1}}{6}\Big{)}L_{\rho}\,,$
$\displaystyle R_{qq;4}^{2;\Delta\Delta_{P}}$ $\displaystyle=$ $\displaystyle
C_{A}+\Big{(}\frac{3C_{A}}{2}+\frac{4C_{1}}{3}\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qq;5}^{2;\Delta\Delta_{P}}$ $\displaystyle=$ $\displaystyle
R_{qq;6}^{2;\Delta\Delta_{P}}\,=\,-\Big{(}\frac{3C_{A}}{4}-\frac{C_{1}}{12}\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qq;7}^{2;\Delta\Delta_{P}}$ $\displaystyle=$ $\displaystyle
R_{qq;8}^{2;\Delta\Delta_{P}}\,=\,\Big{(}\frac{C_{A}}{2}-\frac{C_{1}}{6}\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qq;9}^{2;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{4}\,L_{\rho}\,.$ (S12) $\displaystyle
P_{qq,1;1}^{2;\Delta\Delta_{P}}(z)$ $\displaystyle=$ $\displaystyle
P_{qq,2;1}^{2;\Delta\Delta_{P}}(z)\,=\,\frac{3C_{A}}{16}\,\frac{(1+z)(1+5z-3z^{2}+z^{3})}{(1-z)^{3}_{+,2}}$
$\displaystyle+$
$\displaystyle\frac{C_{1}}{48}\,\frac{9+42z-60z^{2}-18z^{3}+19z^{4}}{(1-z)^{3}_{+,2}}+\frac{3C_{1}}{8}\,f_{4}(z)\,(1+6z+z^{2})\,,$
$\displaystyle P_{qq,1;2}^{2;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle
P_{qq,2;2}^{2;\Delta\Delta_{P}}(z)\,=\,-\frac{C_{A}}{4}\,\frac{z(3+z+3z^{2}-3z^{3})}{(1-z)^{3}_{+,2}}-\frac{C_{1}}{12}\frac{z(9-9z-z^{2}-3z^{3})}{(1-z)^{3}_{+,2}}$
$\displaystyle-$ $\displaystyle\frac{C_{1}}{2}\,f_{4}(z)\,z(3+z)\,,$
$\displaystyle P_{qq,1;3}^{2;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle
P_{qq,2;3}^{2;\Delta\Delta_{P}}(z)\,=\,\frac{C_{A}}{4}\,\frac{z^{2}(1+z^{2})}{(1-z)^{3}_{+,2}}+\frac{C_{1}}{12}\frac{z^{2}(3-z)}{(1-z)^{2}_{+,1}}+\frac{C_{1}}{2}\,f_{4}(z)\,z^{2}\,.$
(S13)

For the $\Delta_{P}$ structure function we have, in the $q\bar{q}$
annihilation channel:

$\displaystyle R_{qq;2}^{1;\Delta_{P}}$ $\displaystyle=$ $\displaystyle-
R_{qq;3}^{1;\Delta_{P}}\,=\,-\frac{L_{\rho}}{\rho}\,(2C_{A}+C_{1})\,.$ (S14)

$\displaystyle P_{qq,1;1}^{1;\Delta_{P}}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq,2;1}^{1;\Delta_{P}}(z)\,=\,\frac{C_{A}}{2\rho}\,\frac{1+3z-z^{2}+z^{3}}{(1-z)^{2}_{+,1}}+\frac{C_{1}}{2\rho}\,\frac{1+2z-3z^{2}+2z^{3}}{(1-z)^{2}_{+,1}}$
$\displaystyle+$ $\displaystyle\frac{C_{1}}{2\rho}\,f_{1}(z)\,(1+3z)\,,$
$\displaystyle P_{qq,1;2}^{1;\Delta_{P}}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq,2;2}^{1;\Delta_{P}}(z)\,=\,-\frac{C_{A}}{\rho}\,\frac{z(1+z)}{(1-z)_{+}}-\frac{C_{1}}{\rho}\,\frac{z}{(1-z)_{+}}-\frac{C_{1}}{\rho}\,f_{1}(z)\,z\,.$
(S15)

$\displaystyle R_{qq;2}^{2;\Delta_{P}}$ $\displaystyle=$ $\displaystyle-
R_{qq;3}^{2;\Delta_{P}}\,=\,\frac{2C_{A}+C_{1}}{\rho}\,+\frac{3}{\rho}\,\Big{(}C_{A}+\frac{19C_{1}}{18}\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qq;5}^{2;\Delta_{P}}$ $\displaystyle=$ $\displaystyle-
R_{qq;6}^{2;\Delta_{P}}\,=\,\frac{3}{2\rho}\,\Big{(}C_{A}+\frac{2C_{1}}{9}\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qq;7}^{2;\Delta_{P}}$ $\displaystyle=$ $\displaystyle-
R_{qq;8}^{2;\Delta_{P}}\,=\,-\frac{2C_{A}+C_{1}}{2\rho}\,L_{\rho}\,.$ (S16)

$\displaystyle P_{qq,1;1}^{2;\Delta_{P}}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq,2;1}^{2;\Delta_{P}}(z)\,=\,-\frac{3C_{A}}{8\rho}\,\frac{1+6z+2z^{3}-z^{4}}{(1-z)^{3}_{+,2}}$
$\displaystyle-$
$\displaystyle\frac{C_{1}}{24\rho}\frac{12+57z-75z^{2}+39z^{3}-17z^{4}}{(1-z)^{3}_{+,2}}-\frac{3C_{1}}{8\rho}\,f_{3}(z)\,(1+6z+z^{2})\,,$
$\displaystyle P_{qq,1;2}^{2;\Delta_{P}}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq,2;2}^{2;\Delta_{P}}(z)\,=\,\frac{C_{A}}{2\rho}\,\frac{z(3+z-3(1-z)z^{2})}{(1-z)^{3}_{+,2}}$
$\displaystyle+$
$\displaystyle\frac{C_{1}}{6\rho}\frac{z(12-11z-2z^{2}+7z^{3})}{(1-z)^{3}_{+,2}}+\frac{C_{1}}{2\rho}\,f_{3}(z)\,z(3+z)\,,$
$\displaystyle P_{qq,1;3}^{2;\Delta_{P}}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq,2;3}^{2;\Delta_{P}}(z)\,=\,-\frac{C_{A}}{2\rho}\,\frac{z^{2}(1+z)}{(1-z)^{2}_{+,1}}-\frac{C_{1}}{6\rho}\frac{z^{2}(4-z)}{(1-z)^{3}_{+,2}}-\frac{C_{1}}{2\rho}\,f_{3}(z)\,z^{2}\,.$
(S17)

For the $\nabla$ structure function we have, in the $q\bar{q}$ annihilation
channel:

$\displaystyle R_{qq;N}^{1;\nabla}=2R_{qq;N}^{1;\Delta_{P}}\,,\quad N=2,3\,.$
(S18)

$\displaystyle P_{qq,1;1}^{1;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,1;1}^{1;\Delta_{P}}(z)+\frac{C_{1}}{\rho}\,\Big{(}1-f_{1}(z)\,(1+z)\Big{)}\,,$
$\displaystyle P_{qq,1;2}^{1;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,1;2}^{1;\Delta_{P}}(z)\,,$ $\displaystyle P_{qq,2;N}^{1;\nabla}(z)$
$\displaystyle=$ $\displaystyle-P_{qq,1;N}^{1;\nabla}(z)\,,\quad N=1,2\,.$
(S19)

$\displaystyle R_{qq;2}^{2;\nabla}$ $\displaystyle=$ $\displaystyle-
R_{qq;3}^{2;\nabla}=2R_{qq;2}^{2;\Delta_{P}}+\frac{C_{1}}{3\rho}\,L_{\rho}\,,$
$\displaystyle R_{qq;5}^{2;\nabla}$ $\displaystyle=$ $\displaystyle-
R_{qq;6}^{2;\nabla}\,=\,2R_{qq;5}^{2;\Delta_{P}}\,=\,-2R_{qq;6}^{2;\Delta_{P}}\,,$
$\displaystyle R_{qq;7}^{2;\nabla}$ $\displaystyle=$ $\displaystyle-
R_{qq;8}^{2;\nabla}\,=\,2R_{qq;7}^{2;\Delta_{P}}\,=\,-2R_{qq;8}^{2;\Delta_{P}}\,.$
(S20)

$\displaystyle P_{qq,1;1}^{2;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,1;1}^{2;\Delta_{P}}(z)$ $\displaystyle-$
$\displaystyle\frac{C_{A}}{4\rho}\,(1+z)-\frac{C_{1}}{12\rho}\,\frac{12-13z-2z^{2}+7z^{3}}{(1-z)^{2}_{+,1}}$
$\displaystyle+$ $\displaystyle\frac{C_{1}}{4\rho}\,f_{3}(z)\,(1+3z)(3+z)\,,$
$\displaystyle P_{qq,1;2}^{2;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,1;2}^{2;\Delta_{P}}(z)$ $\displaystyle+$
$\displaystyle\frac{C_{1}}{3\rho}\,\frac{z(2-3z+z^{2})}{(1-z)^{2}_{+,1}}-\frac{C_{1}}{\rho}\,f_{3}(z)\,z(1+z)\,,$
$\displaystyle P_{qq,2;N}^{2;\nabla}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq,1;N}^{2;\nabla}(z)\,,\quad N=1,2\,.$ (S21)

We now present the results for the $qg$ channel. For the $\Delta\Delta_{P}$
structure function, we have:

$\displaystyle R_{qg;1}^{1;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle\Big{(}\frac{5C_{A}}{2}+4C_{1}\Big{(}\frac{11}{12}+L_{\rho}\Big{)}\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qg;2}^{1;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\Big{(}\frac{3C_{A}}{2}+2C_{1}(1+L_{\rho})\Big{)}\,L_{\rho}\,.$
(S22)

$\displaystyle P_{qg,1;1}^{1;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\frac{z^{2}(1+2z)}{(1-z)^{2}_{+,1}}+C_{1}\Big{(}\frac{1+5z}{3}+2zL_{\rho}\Big{)}\,\frac{z^{2}}{(1-z)^{2}_{+,1}}$
$\displaystyle+$ $\displaystyle C_{1}\,f_{3}(z)\,z^{2}\,,$ $\displaystyle
P_{qg,2;1}^{1;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{4}\,\frac{1+7z+2z^{2}}{(1-z)_{+}}-\frac{C_{1}}{6}\,\frac{3+33z-21z^{2}+7z^{3}+6z(1+3z)L_{\rho}}{(1-z)_{+}}\,,$
$\displaystyle P_{qg,2;2}^{1;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\,\frac{z(1+2z)}{(1-z)_{+}}+C_{1}\,\frac{z(1+2z-z^{2}+2zL_{\rho})}{(1-z)_{+}}\,.$
(S23)

$\displaystyle R_{qg;1}^{2;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{4}\,(10-3L_{\rho})-\frac{C_{1}}{3}\,(11+21L_{\rho}-18L_{\rho}^{2})\,,$
$\displaystyle R_{qg;2}^{2;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{4}\,(6-5L_{\rho})+C_{1}\,(2+2L_{\rho}-5L_{\rho}^{2})\,,$
$\displaystyle R_{qg;3}^{2;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\Big{(}\frac{7C_{A}}{2}+\frac{C_{1}}{2}\,(11+12L_{\rho})\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qg;4}^{2;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle\Big{(}\frac{5C_{A}}{2}+\frac{C_{1}}{3}\,(11+12L_{\rho})\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qg;5}^{2;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle\Big{(}\frac{5C_{A}}{4}+C_{1}\,(1+2L_{\rho})\Big{)}\,L_{\rho}\,,$
$\displaystyle R_{qg;7}^{2;\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\Big{(}\frac{3C_{A}}{4}+C_{1}\,(1+L_{\rho})\Big{)}\,L_{\rho}\,.$
(S24)

$\displaystyle P_{qg,1;1}^{2;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{4}\frac{z^{2}(3+9z-2z^{2})}{(1-z)^{3}_{+,2}}$
$\displaystyle+$
$\displaystyle\frac{C_{1}}{4}\Big{(}5-27z+15z^{2}-z^{3}-4z(3+z)L_{\rho}\Big{)}\,\frac{z^{2}}{(1-z)^{3}_{+,2}}$
$\displaystyle-$ $\displaystyle\frac{3C_{1}}{2}\,f_{4}(z)\,z^{2}(1+z)\,,$
$\displaystyle P_{qg,1;2}^{2;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\,\frac{z^{3}(1+2z)}{(1-z)^{3}_{+,2}}+\frac{C_{1}}{6}\frac{z^{3}(3+8z+z^{2}+12zL_{\rho})}{(1-z)^{3}_{+,2}}+C_{1}\,f_{4}(z)\,z^{3}\,,$
$\displaystyle P_{qg,2;1}^{2;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{16}\,\frac{3+36z+19z^{2}-2z^{3}}{(1-z)^{2}_{+,1}}$
$\displaystyle+$
$\displaystyle\frac{C_{1}}{8}\,\Big{(}3+52z-8z^{2}+4z^{3}-7z^{4}+6z(1+6z+z^{2})L_{\rho}\Big{)}\frac{1}{(1-z)^{2}_{+,1}}\,,$
$\displaystyle P_{qg,2;2}^{2;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{4}\,\frac{z(3+9z-2z^{2})}{(1-z)^{2}_{+,1}}$
$\displaystyle-$
$\displaystyle\frac{C_{1}}{6}\,\Big{(}9+39z-39z^{2}+13z^{3}+6z(3+z)L_{\rho}\Big{)}\,\frac{z}{(1-z)^{2}_{+,1}}\,,$
$\displaystyle P_{qg,2;3}^{2;\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{4}\,\frac{z^{2}(1+2z)}{(1-z)^{2}_{+,1}}+\frac{C_{1}}{2}\,\frac{z^{2}(1+2z-z^{2}+2zL_{\rho})}{(1-z)^{2}_{+,1}}\,.$
(S25)

For the $\Delta_{P}$ structure function we have, in the $qg$ channel:

$\displaystyle R_{qg;1}^{1;\Delta_{P}}$ $\displaystyle=$
$\displaystyle\frac{L_{\rho}}{2\rho}\,(C_{A}-6C_{1})\,,$ $\displaystyle
R_{qg;2}^{1;\Delta_{P}}$ $\displaystyle=$
$\displaystyle\frac{2L_{\rho}}{\rho}\,C_{1}\,.$ (S26)

$\displaystyle P_{qg,1;1}^{1;\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{2\rho}\,\frac{z^{2}}{(1-z)_{+}}-\frac{C_{1}}{\rho}\,\frac{z^{2}(1+z)}{(1-z)^{2}_{+,1}}\,,$
$\displaystyle P_{qg,2;1}^{1;\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{2\rho}\,\frac{1+z-z^{2}}{(1-z)_{+}}$
$\displaystyle-$
$\displaystyle\frac{C_{1}}{\rho}\,\frac{1+4z-10z^{2}+2z^{3}}{(1-z)_{+}}-\frac{C_{1}}{\rho}\,z(1+3z)L_{\rho}\,,$
$\displaystyle P_{qg,2;2}^{1;\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{\rho}\,z+\frac{2C_{1}}{\rho}\,\frac{z(1-2z^{2})}{(1-z)_{+}}+\frac{2C_{1}}{\rho}\,z^{2}\,L_{\rho}\,.$
(S27)

$\displaystyle R_{qg;1}^{2;\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{2\rho}\Big{(}1+\frac{3}{2}L_{\rho}\Big{)}+\frac{3C_{1}}{\rho}\Big{(}1-\frac{19}{6}L_{\rho}-2L_{\rho}^{2}\Big{)}\,,$
$\displaystyle R_{qg;2}^{2;\Delta_{P}}$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2\rho}L_{\rho}-\frac{2C_{1}}{\rho}(1-4L_{\rho}-2L_{\rho}^{2})\,,$
$\displaystyle R_{qg;3}^{2;\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{L_{\rho}}{\rho}(C_{A}-4C_{1})\,,$ $\displaystyle
R_{qg;4}^{2;\Delta_{P}}$ $\displaystyle=$
$\displaystyle\frac{L_{\rho}}{2\rho}\,(C_{A}-6C_{1})\,,$ $\displaystyle
R_{qg;5}^{2;\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{L_{\rho}}{2\rho}\,(C_{A}+2C_{1}(5+L_{\rho}))\,,$
$\displaystyle R_{qg;7}^{2;\Delta_{P}}$ $\displaystyle=$
$\displaystyle\frac{L_{\rho}}{\rho}\,C_{1}\,.$ (S28)

$\displaystyle P_{qg,1;1}^{2;\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{1}{4\rho}\,\Big{(}C_{A}(3+z^{2})+2C_{1}(7+10z+3z^{2}+4z^{2}L_{\rho})\Big{)}\,\frac{z^{2}}{(1-z)^{3}_{+,2}}\,,$
$\displaystyle P_{qg,1;2}^{2;\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{2\rho}\,\frac{z^{3}}{(1-z)^{2}_{+,1}}-\frac{C_{1}}{\rho}\,\frac{z^{3}(1+z)}{(1-z)^{3}_{+,2}}\,,$
$\displaystyle P_{qg,2;1}^{2;\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{8\rho}\,\frac{3+9z-5z^{2}+z^{3}}{(1-z)^{2}_{+,1}}$
$\displaystyle+$
$\displaystyle\frac{C_{1}}{4\rho}\,\frac{3+20z-37z^{2}-8z^{3}+6z^{4}}{(1-z)^{2}_{+,1}}+\frac{3C_{1}}{4\rho}\,\frac{z(1+6z+z^{2})}{(1-z)_{+}}\,L_{\rho}\,\,,$
$\displaystyle P_{qg,2;2}^{2;\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{2\rho}\,\frac{z(3-3z+z^{2})}{(1-z)^{2}_{+,1}}$
$\displaystyle-$
$\displaystyle\frac{C_{1}}{\rho}\,\frac{z(3+2z-14z^{2}+6z^{3})}{(1-z)^{2}_{+,1}}-\frac{C_{1}}{\rho}\,\frac{z^{2}(3+z)}{(1-z)_{+}}\,L_{\rho}\,,$
$\displaystyle P_{qg,2;3}^{2;\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2\rho}\,\frac{z^{2}}{(1-z)_{+}}+\frac{C_{1}}{\rho}\,\frac{z^{2}(1-2z^{2})}{(1-z)^{2}_{+,1}}+\frac{C_{1}}{\rho}\,\frac{z^{3}}{(1-z)_{+}}L_{\rho}\,.$
(S29)

Finally, for the $\nabla$ structure function we have, in the $qg$ channel:

$\displaystyle R_{qg;1}^{1;\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;1}^{1;\Delta_{P}}-\frac{2L_{\rho}}{\rho}\,(C_{A}-6C_{1})\,,$
$\displaystyle R_{qg;2}^{1;\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;2}^{1;\Delta_{P}}\,.$ (S30)

$\displaystyle P_{qg,1;1}^{1;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,1;1}^{1;\Delta_{P}}(z)+\frac{2}{\rho}\,\Big{(}C_{A}z-2C_{1}(2+z)\Big{)}\,\frac{z}{(1-z)_{+}}\,,$
$\displaystyle P_{qg,2;1}^{1;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;1}^{1;\Delta_{P}}(z)+\frac{1}{\rho}\,\frac{C_{A}-6C_{1}z}{(1-z)_{+}}+\frac{2C_{1}}{\rho}\,z(1+z)\,L_{\rho}\,,$
$\displaystyle P_{qg,2;2}^{1;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;2}^{1;\Delta_{P}}(z)\,.$ (S31)

$\displaystyle R_{qg;1}^{2;\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;1}^{2;\Delta_{P}}+\frac{C_{A}}{2\rho}(4-7L_{\rho})-\frac{C_{1}}{3\rho}\Big{(}36-47L_{\rho}-24L_{\rho}^{2}\Big{)}\Big{)}\,,$
$\displaystyle R_{qg;2}^{2;\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;2}^{2;\Delta_{P}}+\frac{L_{\rho}}{\rho}\,\Big{(}3C_{A}-2C_{1}(15+2L_{\rho})\Big{)}\,,$
$\displaystyle R_{qg;3}^{2;\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;3}^{2;\Delta_{P}}+\frac{4L_{\rho}}{\rho}\,(C_{A}-4C_{1})\,,$
$\displaystyle R_{qg;4}^{2;\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;4}^{2;\Delta_{P}}-\frac{2L_{\rho}}{\rho}\,(C_{A}-6C_{1})\,.$ (S32)

$\displaystyle P_{qg,1;1}^{2;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,1;1}^{2;\Delta_{P}}(z)$ $\displaystyle-$
$\displaystyle\frac{C_{A}}{2\rho}\,\frac{z^{2}(5+z)}{(1-z)^{2}_{+,1}}+\frac{C_{1}}{3\rho}\Big{(}48+37z+z^{2}(5+12L_{\rho})\Big{)}\,\frac{z}{(1-z)^{2}_{+,1}}$
$\displaystyle-$ $\displaystyle\frac{2C_{1}}{\rho}\,f_{3}(z)\,z^{2}\,,$
$\displaystyle P_{qg,1;2}^{2;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,1;2}^{2;\Delta_{P}}(z)+\frac{2}{\rho}\,\Big{(}C_{A}z-2C_{1}(2+z)\Big{)}\,\frac{z^{2}}{(1-z)^{2}_{+,1}}\,,$
$\displaystyle P_{qg,2;1}^{2;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,1;2}^{2;\Delta_{P}}(z)+\frac{C_{A}}{4\rho}\,\frac{1-5z-11z^{2}-z^{3}}{(1-z)^{2}_{+,1}}$
$\displaystyle+$
$\displaystyle\frac{C_{1}}{6\rho}\,\frac{3-12z+63z^{2}+92z^{3}-50z^{4}}{(1-z)^{2}_{+,1}}-\frac{C_{1}}{2\rho}\,\frac{z(1+3z)(3+z)}{(1-z)_{+}}\,L_{\rho}\,,$
$\displaystyle P_{qg,2;2}^{2;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;2}^{2;\Delta_{P}}(z)+\frac{1}{\rho}\,\Big{(}C_{A}-6C_{1}z)\,\frac{z^{2}(1+z)}{(1-z)^{2}_{+,1}}+\frac{2C_{1}}{\rho}\,\frac{z^{2}(1+z)}{(1-z)_{+}}\,L_{\rho}\,,$
$\displaystyle P_{qg,2;3}^{2;\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;3}^{2;\Delta_{P}}(z)\,.$ (S33)

The partial derivatives
$x_{1}^{m}x_{2}^{n}\,\partial^{m}_{x_{1}}\partial^{n}_{x_{2}}\
W^{ab}_{0;J}(x_{1},x_{2})$ and
$x_{1}^{m}x_{2}^{n}\,\partial^{m}_{x_{1}}\partial^{n}_{x_{2}}\
W^{ab}_{1;J}(x_{1},x_{2})$ are evaluated using the technique discussed in
Appendix C. To the accuracy we are working at, we need the derivatives
$\partial_{x_{i}}$, $\partial_{x_{1}}^{2}$ with $i=1,2$ and
$\partial_{x_{1}x_{2}}^{2}$ acting on the $W_{0}$ structure functions and
$\partial_{x_{i}}$ with $i=1,2$ acting on the $W_{1}$ structure functions.

These derivatives can be written as

$\displaystyle x_{i}\,\partial_{x_{i}}W_{k}=-W_{k}+V_{k}^{(i)}\,,$
$\displaystyle
x_{i}^{2}\,\partial_{x_{i}}^{2}W_{k}=-(1+3x_{i}\,\partial_{x_{i}})\,W_{k}+V_{k}^{(ii)}\,,$
$\displaystyle
x_{1}x_{2}\,\partial_{x_{1}x_{2}}^{2}W_{k}=-(1+x_{1}\,\partial_{x_{1}}+x_{2}\,\partial_{x_{2}})\,W_{k}+V_{k}^{(12)}\,,$
(S34)

where $V_{k}^{(L)}$ with $L=1,2,11,22,12$ are the specific expressions for the
corresponding derivatives, subprocess types, and type of structure function.
Following Eq. (S3) we write the $V_{k}^{(L)}$ as

$\displaystyle V^{ab;(L)}_{k;J}(x_{1},x_{2},L_{\rho})$ $\displaystyle=$
$\displaystyle-\frac{g_{ab}}{4x_{1}x_{2}}\biggl{[}R_{ab,k}^{(L);J}(x_{1},x_{2},L_{\rho})\,f_{a/H_{1}}(x_{1})\,f_{b/H_{2}}(x_{2})$
(S35) $\displaystyle+$ $\displaystyle\Big{(}P_{ba,k}^{(L);J}\otimes
f_{b/H_{2}}\Big{)}(x_{2},x_{1},L_{\rho})\ f_{a/H_{1}}(x_{1})$ $\displaystyle+$
$\displaystyle\Big{(}P_{ab,k}^{(L);J}\otimes
f_{a/H_{1}}\Big{)}(x_{1},x_{2},L_{\rho})\ f_{b/H_{2}}(x_{2})\biggr{]}\,,$

and expand the perturbative functions
$R_{ab,k}^{(L);J}(x_{1},x_{2},L_{\rho})$,
$P_{ab,k}^{(L);J}(z_{1},x_{2},L_{\rho})$, and
$P_{ba,k}^{(L);J}(z_{2},x_{1},L_{\rho})$ following Eq. (Supplementary
Material: NLP and NNLP contributions to the $T$-odd hadronic structure
functions $W_{i}(x_{1},x_{2})$):

$\displaystyle R_{ab,k}^{(L);J}(x_{1},x_{2},L_{\rho})$ $\displaystyle=$
$\displaystyle\sum\limits_{N=1}^{9}\,R_{ab;N}^{k(L);J}(L_{\rho})\,T^{R}_{N}(x_{1},x_{2})\,,$
$\displaystyle P_{ab,k}^{(L);J}(z_{1},x_{2},L_{\rho})$ $\displaystyle=$
$\displaystyle\sum\limits_{N=1}^{3}\,P_{ab,1;N}^{k(L);J}(z_{1},L_{\rho})\,T^{P}_{N}(x_{2})\,,$
$\displaystyle P_{ba,k}^{(L);J}(z_{2},x_{1},L_{\rho})$ $\displaystyle=$
$\displaystyle\sum\limits_{N=1}^{3}\,P_{ab,2;N}^{k(L);J}(z_{2},L_{\rho})\,T^{P}_{N}(x_{1})\,,$
(S36)

Below we again list only the nonvanishing $R_{ab;N}^{k(L);J}(L_{\rho})$,
$P_{ab,1;N}^{k(L);J}(z_{1},L_{\rho})$, $P_{ab,2;N}^{k(L);J}(z_{2},L_{\rho})$.

For the $\Delta\Delta_{P}$ structure function we have, in the $q\bar{q}$
annihilation channel:

$\displaystyle R_{qq;1}^{0(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle
R_{qq;1}^{0(2);\Delta\Delta_{P}}\,=\,R_{qq;2}^{0(12);\Delta\Delta_{P}}\,=\,R_{qq;3}^{0(12);\Delta\Delta_{P}}=-\frac{C_{1}}{3}\,,$
$\displaystyle R_{qq;2}^{0(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle R_{qq;3}^{0(2);\Delta\Delta_{P}}\,=\,-C_{A}\,(1+L_{\rho})\,,$
$\displaystyle R_{qq;1}^{0(11);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle
R_{qq;1}^{0(22);\Delta\Delta_{P}}\,=\,\frac{C_{A}}{2}+\frac{C_{1}}{6}\,,$
$\displaystyle R_{qq;2}^{0(11);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle
R_{qq;3}^{0(22);\Delta\Delta_{P}}\,=\,-C_{A}\,\Big{(}2+L_{\rho}\Big{)}-\frac{C_{1}}{3}\,,$
$\displaystyle R_{qq;4}^{0(12);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-C_{A}\,\Big{(}2+L_{\rho}\Big{)}\,,$ $\displaystyle
R_{qq;5}^{0(11);\Delta\Delta_{P}}$ $\displaystyle=$ $\displaystyle
R_{qq;6}^{0(22);\Delta\Delta_{P}}\,=\,-\frac{C_{A}}{2}\,(3+2L_{\rho})\,.$
(S37)

$\displaystyle P_{qq,1;1}^{0(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle P_{qq,2;1}^{0(2);\Delta\Delta_{P}}(z)\ =\
P_{qq,1;2}^{0(12);\Delta\Delta_{P}}(z)\ =\
P_{qq,2;2}^{0(12);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\frac{z(1+2z-z^{2})}{(1-z)^{2}_{+,1}}+C_{1}z\partial_{z}f_{1}(z)\,,$
$\displaystyle P_{qq,2;2}^{0(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle
P_{qq,1;2}^{0(2);\Delta\Delta_{P}}(z)\,=\,P_{qq,2;2}^{0(11);\Delta\Delta_{P}}(z)\,=\,P_{qq,2;3}^{0(11);\Delta\Delta_{P}}(z)\,=\,$
$\displaystyle=$ $\displaystyle
P_{qq,1;2}^{0(22);\Delta\Delta_{P}}(z)\,=\,P_{qq,1;3}^{0(22);\Delta\Delta_{P}}(z)\,=\,\frac{C_{A}}{2}\frac{1+z^{2}}{(1-z)_{+}}+C_{1}f_{1}(z)\,,$
$\displaystyle P_{qq,1;1}^{0(11);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle P_{qq,2;1}^{0(22);\Delta\Delta_{P}}(z)$ (S38) $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\frac{z(1+5z-3z^{2}+z^{3})}{(1-z)^{3}_{+,2}}+C_{1}\,\Big{(}z\partial_{z}+z^{2}\partial_{z}^{2}\Big{)}f_{1}(z)\,.$

$\displaystyle R_{qq;1}^{1(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle R_{qq;1}^{1(2);\Delta\Delta_{P}}=\frac{C_{1}}{6}\,,$
$\displaystyle R_{qq;2}^{1(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle
R_{qq;3}^{1(2);\Delta\Delta_{P}}\,=\,\frac{5}{2}C_{A}+C_{1}+L_{\rho}\Big{(}\frac{3}{2}C_{A}+C_{1}\Big{)}\,,$
$\displaystyle R_{qq;3}^{1(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle
R_{qq;2}^{1(2);\Delta\Delta_{P}}\,=\,-\Big{(}\frac{C_{A}}{2}+\frac{C_{1}}{6}\Big{)}\,,$
$\displaystyle R_{qq;4}^{1(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle
R_{qq;4}^{1(2);\Delta\Delta_{P}}\,=\,-\frac{C_{1}}{3}\Big{(}1+L_{\rho}\Big{)}\,,$
$\displaystyle R_{qq;5}^{1(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle
R_{qq;6}^{1(2);\Delta\Delta_{P}}\,=\,\Big{(}C_{A}-\frac{C_{1}}{3}\Big{)}\Big{(}\frac{1}{2}+L_{\rho}\Big{)}\,,$
$\displaystyle R_{qq;7}^{1(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle
R_{qq;8}^{1(2);\Delta\Delta_{P}}\,=\,-C_{A}\Big{(}\frac{1}{2}+L_{\rho}\Big{)}\,.$

$\displaystyle P_{qq,1;1}^{1(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle P_{qq,2;1}^{1(2);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{4}\frac{z(5+5z-3z^{2}+z^{3})}{(1-z)^{3}_{+,2}}-\frac{C_{1}}{3}\frac{z(1+z)(2-4z+z^{2})}{(1-z)^{3}_{+,2}}$
$\displaystyle-$
$\displaystyle\frac{C_{1}}{2}z\partial_{z}\Big{(}f_{3}(z)\,(1+3z)\Big{)}\,,$
$\displaystyle P_{qq,1;2}^{1(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle P_{qq,2;2}^{1(2);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\frac{z(1+z+3z^{2}-z^{3})}{(1-z)^{3}_{+,2}}+\frac{C_{1}}{3}\frac{z}{(1-z)^{2}_{+,1}}$
$\displaystyle+$ $\displaystyle
C_{1}z\partial_{z}\Big{(}f_{3}(z)\,z\Big{)}\,,$ $\displaystyle
P_{qq,2;2}^{1(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$ $\displaystyle
P_{qq,1;2}^{1(2);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{4}\frac{1+z+z^{2}-3z^{3}}{(1-z)^{2}_{+,1}}-\frac{C_{1}}{6}\frac{1-z^{2}-2z^{3}}{(1-z)^{2}_{+,1}}$
$\displaystyle-$ $\displaystyle\frac{C_{1}}{2}f_{3}(z)\,(1+z)\,,$
$\displaystyle P_{qq,2;3}^{1(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle P_{qq,1;3}^{1(2);\Delta\Delta_{P}}(z)$ (S40) $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\frac{z(1+z^{2})}{(1-z)^{2}_{+,1}}+\frac{C_{1}}{3}\frac{z}{(1-z)_{+}}$
$\displaystyle+$ $\displaystyle C_{1}f_{3}(z)\,z\,.$

For the $\Delta_{P}$ and $\nabla$ structure functions we have, in the
$q\bar{q}$ annihilation channel:

$\displaystyle R_{qq;1}^{0(1);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qq;1}^{0(1);\Delta_{p}}\,=\,-R_{qq;1}^{0(2);\nabla}\,=\,-2R_{qq;1}^{0(2);\Delta_{p}}$
$\displaystyle R_{qq;2}^{0(11);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qq;2}^{0(11);\Delta_{P}}\,=\,-R_{qq;3}^{0(22);\nabla}\,=\,-2R_{qq;3}^{0(22);\Delta_{P}}$
$\displaystyle R_{qq;3}^{0(12);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qq;3}^{0(12);\Delta_{P}}\,=\,-R_{qq;2}^{0(12);\nabla}\,=\,-2R_{qq;2}^{0(12);\Delta_{P}}$
$\displaystyle=$ $\displaystyle\frac{2}{\rho}\,\Big{(}2C_{A}+C_{1}\Big{)}\,,$
$\displaystyle R_{qq;1}^{0(11);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qq;1}^{0(11);\Delta_{P}}\,=\,-R_{qq;1}^{0(22);\nabla}\,=\,-2R_{qq;1}^{0(22);\Delta_{P}}$
(S41) $\displaystyle=$
$\displaystyle\frac{2}{\rho}\,\Big{(}C_{A}-\frac{C_{1}}{3}\Big{)}\,.$

$\displaystyle P_{qq,1;1}^{0(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,1;1}^{0(1);\Delta_{P}}(z)\,=\,-P_{qq,2;1}^{0(2);\nabla}(z)\,=\,-2P_{qq;2;1}^{0(2);\Delta_{P}}(z)$
$\displaystyle P_{qq,1;2}^{0(12);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,1;2}^{0(12);\Delta_{P}}(z)\,=\,-P_{qq,2;2}^{0(12);\nabla}(z)\,=\,-2P_{qq;2;2}^{0(12);\Delta_{P}}(z)$
$\displaystyle=$
$\displaystyle-\frac{2}{\rho}\,z\,\Big{(}C_{A}+2C_{1}\partial_{z}f_{2}(z)\Big{)}\,,$
$\displaystyle P_{qq,2;1}^{0(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,2;1}^{0(1);\Delta_{P}}(z)\,=\,-P_{qq,1;2}^{0(2);\nabla}(z)\,=\,-2P_{qq;1;2}^{0(2);\Delta_{P}}(z)$
$\displaystyle P_{qq,2;2}^{0(11);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,2;2}^{0(11);\Delta_{P}}(z)\,=\,-P_{qq,1;2}^{0(22);\nabla}(z)\,=\,-2P_{qq;1;2}^{0(22);\Delta_{P}}(z)$
$\displaystyle P_{qq,2;3}^{0(11);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,2;3}^{0(11);\Delta_{P}}(z)\,=\,-P_{qq,1;3}^{0(22);\nabla}(z)\,=\,-2P_{qq;1;3}^{0(22);\Delta_{P}}(z)$
$\displaystyle=$
$\displaystyle\frac{2}{\rho}\,\Big{(}C_{A}\,(1+z)+2C_{1}f(z)\Big{)}\,,$
$\displaystyle P_{qq,1;1}^{0(11);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qq,1;1}^{0(11);\Delta_{P}}(z)\,=\,-P_{qq,2;1}^{0(22);\nabla}(z)\,=\,-2P_{qq;2;1}^{0(22);\Delta_{P}}(z)$
(S42) $\displaystyle=$
$\displaystyle-\frac{2}{\rho}\,z\,\Big{(}C_{A}+2C_{1}(\partial_{z}+z\partial_{z}^{2})f_{2}(z)\Big{)}\,.$

$\displaystyle R_{qq;1}^{1(1);\nabla}$ $\displaystyle=$ $\displaystyle-
R_{qq;1}^{1(2);\nabla}\,=\,2R_{qq;1}^{1(1);\Delta_{P}}-\frac{C_{1}}{3\rho}\,=\,-\Big{(}2R_{qq;1}^{1(2);\Delta_{P}}-\frac{C_{1}}{3\rho}\Big{)}=-\frac{2}{\rho}\,\Big{(}C_{A}+\frac{7C_{1}}{3}\Big{)}\,,$
$\displaystyle R_{qq;2}^{1(1);\nabla}$ $\displaystyle=$ $\displaystyle-
R_{qq;3}^{1(2);\nabla}\,=\,2R_{qq;2}^{1(1);\Delta_{P}}\,=\,-2R_{qq;3}^{1(2);\Delta_{P}}$
$\displaystyle=$ $\displaystyle-
R_{qq;4}^{1(1);\nabla}\,=\,R_{qq;4}^{1(2);\nabla}\,=\,-2R_{qq;4}^{1(1);\Delta_{P}}\,=\,2R_{qq;4}^{1(2);\Delta_{P}}=-\frac{2(1+L_{\rho})}{\rho}\,\Big{(}2C_{A}+C_{1}\Big{)}\,,$
$\displaystyle R_{qq;3}^{1(1);\nabla}$ $\displaystyle=$ $\displaystyle-
R_{qq;2}^{1(2);\nabla}\,=\,2R_{qq;3}^{1(1);\Delta_{P}}\,=\,-2R_{qq;2}^{1(2);\Delta_{P}}\,=\,-\frac{2}{\rho}\,\Big{(}C_{A}+\frac{C_{1}}{3}\Big{)}\,,$
$\displaystyle R_{qq;5}^{1(1);\nabla}$ $\displaystyle=$ $\displaystyle-
R_{qq;6}^{1(2);\nabla}\,=\,2R_{qq;5}^{1(1);\Delta_{P}}\,=\,2R_{qq;6}^{1(2);\Delta_{P}}\,=\,-\frac{1+2L_{\rho}}{\rho}\,\Big{(}2C_{A}+C_{1}\Big{)}\,.$
(S43)

$\displaystyle P_{qq,1;1}^{1(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq;1}^{1(2);\nabla}(z)\,=\,2P_{qq,1;1}^{1(1);\Delta_{P}}(z)-\frac{C_{1}}{\rho}\,\partial_{z}\Big{(}f_{1}(z)\,(1+z)\Big{)}$
$\displaystyle=$
$\displaystyle-\Big{(}2P_{qq,1;1}^{1(2);\Delta_{P}}(z)-\frac{C_{1}}{\rho}\,\partial_{z}\Big{(}f_{1}(z)\,(1+z)\Big{)}\Big{)}$
$\displaystyle=$
$\displaystyle\frac{C_{A}}{\rho}\,\frac{z(5+z+3z^{2}-z^{3})}{(1-z)^{3}_{+,2}}+\frac{2C_{1}}{\rho}\,\frac{z(2-2z+3z^{2}-z^{3})}{(1-z)^{3}_{+,2}}$
$\displaystyle+$
$\displaystyle\frac{C_{1}}{\rho}\,\Big{(}2z\,\partial_{z}\Big{(}f_{1}(z)\,(1+3z)\Big{)}-\partial_{z}\Big{(}f_{1}(z)\,(1+z)\Big{)}\Big{)}\,,$
$\displaystyle P_{qq,1;2}^{1(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq,2;2}^{1(2);\nabla}(z)\,=\,2P_{qq,1;2}^{1(1);\Delta_{P}}(z)\,=\,-2P_{qq,122}^{1(2);\Delta_{P}}(z)$
$\displaystyle=$
$\displaystyle-\frac{2C_{A}}{\rho}\,\frac{z(1+2z-z^{2})}{(1-z)^{2}_{+,1}}-\frac{2C_{1}}{\rho}\,\frac{z}{(1-z)^{2}_{+,1}}-\frac{2C_{1}}{\rho}\,z\,\partial_{z}\Big{(}f_{1}(z)\,z\Big{)}\,,$
$\displaystyle P_{qq,2;2}^{1(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq,1;2}^{1(2);\nabla}(z)\,=\,2P_{qq,2;2}^{1(1);\Delta_{P}}(z)-\frac{C_{1}}{\rho}\,\Big{(}1-f_{1}(z)\,(1+z)\Big{)}$
$\displaystyle=$
$\displaystyle-\Big{(}2P_{qq,1;2}^{1(2);\Delta_{P}}(z)-\frac{C_{1}}{\rho}\,\Big{(}1-f_{1}(z)\,(1+z)\Big{)}\Big{)}$
$\displaystyle=$
$\displaystyle-\frac{C_{A}}{\rho}\,\frac{z(1+3z-z^{2}+z^{3})}{(1-z)^{2}_{+,1}}-\frac{C_{1}}{\rho}\,\frac{1-z+3z^{2}-3z^{3}+2z^{4}}{(1-z)^{2}_{+,1}}$
$\displaystyle+$ $\displaystyle\frac{C_{1}}{\rho}\,f_{1}(z)\,(1-3z^{2})\,,$
$\displaystyle P_{qq,2;3}^{1(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle-
P_{qq,1;3}^{1(2);\nabla}(z)\,=\,2P_{qq,2;3}^{1(1);\Delta_{P}}(z)\,=\,-2P_{qq,1;3}^{1(2);\Delta_{P}}(z)\,=\,$
(S44) $\displaystyle=$
$\displaystyle\frac{2C_{A}}{\rho}\,\frac{z(1+z)}{(1-z)_{+}}+\frac{2C_{1}}{\rho}\,\frac{z}{(1-z)_{+}}+\frac{2C_{1}}{\rho}\,z\,f_{1}(z)\,.$

For the $\Delta\Delta_{P}$ structure function we have, in the $qg$ channel:

$\displaystyle R_{qg;1}^{0(2);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle
R_{qg;3}^{0(22);\Delta\Delta_{P}}\,=\,R_{qg;2}^{0(12);\Delta\Delta_{P}}\,=\,-\frac{3C_{A}}{2}-2C_{1}(1+L_{\rho})\,,$
$\displaystyle R_{qg;1}^{0(22);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-C_{A}-2C_{1}L_{\rho}\,.$ (S45)

$\displaystyle P_{qg,2;2}^{0(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle
P_{qg,2;2}^{0(11);\Delta\Delta_{P}}(z)\,=\,P_{qg,2;3}^{0(11);\Delta\Delta_{P}}(z)$
$\displaystyle=$
$\displaystyle\frac{C_{A}}{2}(1+2z)+C_{1}(1+2z-z^{2}+2zL_{\rho})\,,$
$\displaystyle P_{qg,2;1}^{0(2);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle
P_{qg,2;2}^{0(12);\Delta\Delta_{P}}(z)\,=\,z\Big{(}C_{A}+2C_{1}(1-z+L_{\rho})\Big{)}\,,$
$\displaystyle P_{qg,2;1}^{0(22);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle z\Big{(}C_{A}+2C_{1}(1-2z+L_{\rho})\Big{)}\,.$ (S46)

$\displaystyle R_{qg;1}^{1(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-C_{A}-C_{1}\Big{(}\frac{11}{6}+2L_{\rho}\Big{)}\,,$
$\displaystyle R_{qg;2}^{1(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle\Big{(}C_{A}+C_{1}\,\Big{(}\frac{5}{3}+2L_{\rho}\Big{)}\Big{)}\,(1+L_{\rho})\,,$
$\displaystyle R_{qg;5}^{1(1);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\Big{(}\frac{3C_{A}}{4}+C_{1}\,(1+L_{\rho})\Big{)}\,(1+2L_{\rho})\,,$
$\displaystyle R_{qg;1}^{1(2);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{4}-C_{1}\,\Big{(}\frac{5}{3}-3L_{\rho}\Big{)}\,,$
$\displaystyle R_{qg;2}^{1(2);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle C_{A}+2C_{1}\,L_{\rho}\,,$ $\displaystyle
R_{qg;3}^{1(2);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle\frac{5C_{A}}{2}(1+L_{\rho})-4C_{1}\,(1-L_{\rho})\,\Big{(}\frac{11}{12}+L_{\rho}\Big{)}\,,$
$\displaystyle R_{qg;4}^{1(2);\Delta\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\Big{(}\frac{3C_{A}}{2}+2C_{1}\,(1+L_{\rho})\Big{)}\,(1+L_{\rho})\,.$
(S47)

$\displaystyle P_{qg,1;1}^{1(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle
C_{A}\frac{z^{2}(1+3z-z^{2})}{(1-z)^{3}_{+,2}}+\frac{C_{1}}{3}\frac{z^{2}(2+z(3-z)(5+6L_{\rho}))}{(1-z)^{3}_{+,2}}$
$\displaystyle+$ $\displaystyle
C_{1}z\,\partial_{z}\Big{(}f_{3}(z)\,z^{2}\Big{)}\,,$ $\displaystyle
P_{qg,2;2}^{1(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{4}\frac{1+5z-2z^{2}}{(1-z)_{+}}$ $\displaystyle-$
$\displaystyle\frac{C_{1}}{6}\frac{3+27z-33z^{2}+13z^{3}+6(1+z)L_{\rho}}{(1-z)_{+}}\,,$
$\displaystyle P_{qg,2;3}^{1(1);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\frac{z(1+2z)}{(1-z)_{+}}+C_{1}\frac{z(1+2z-z^{2}+2zL_{\rho})}{(1-z)_{+}}\,,$
$\displaystyle P_{qg,2;1}^{1(2);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{C_{A}}{2}\frac{z(4+2z-z^{2})}{(1-z)^{2}_{+,1}}$
$\displaystyle+$
$\displaystyle\frac{C_{1}}{3}\frac{z(18-21z+21z^{2}-7z^{3}+3(1+6z-3z^{2})L_{\rho})}{(1-z)^{2}_{+,1}}\,,$
$\displaystyle P_{qg,2;2}^{1(2);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\frac{z(1+4z-2z^{2})}{(1-z)^{2}_{+,1}}+C_{1}\frac{z(1+4z-5z^{2}+2z^{3}+2z(2-z)L_{\rho})}{(1-z)^{2}_{+,1}}\,,$
$\displaystyle P_{qg,1;2}^{1(2);\Delta\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle\frac{C_{A}}{2}\frac{z^{2}(1+2z)}{(1-z)^{2}_{+,1}}+\frac{C_{1}}{3}\frac{z^{2}(1+5z+6zL_{\rho})}{(1-z)^{2}_{+,1}}+C_{1}z^{2}f_{3}(z)\,.$
(S48)

Finally, for the $\Delta_{P}$ and $\nabla$ structure functions we have, in the
$qg$ channel:

$\displaystyle R_{qg;1}^{0(2);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;1}^{0(2);\Delta_{p}}\,=\,R_{qg;3}^{0(22);\nabla}\,=\,2R_{qg;3}^{0(22);\Delta_{p}}$
$\displaystyle=$ $\displaystyle
R_{qg;2}^{0(12);\nabla}\,=\,2R_{qg;2}^{0(12);\Delta_{p}}\,=\,\frac{4C_{1}}{\rho}\,,$
$\displaystyle R_{qg;1}^{0(22);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;1}^{0(22);\Delta_{p}}\,=\,\frac{2}{\rho}\Big{(}C_{A}+2C_{1}(4+L_{\rho})\Big{)}\,.$
(S49)

$\displaystyle P_{qg,2;2}^{0(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;2}^{0(1);\Delta_{P}}(z)\,=\,P_{qg,2;2}^{0(11);\nabla}(z)\,=\,2P_{qg,2;2}^{0(11);\Delta_{P}}(z)$
$\displaystyle=$ $\displaystyle
P_{qg,2;3}^{0(11);\nabla}(z)\,=\,2P_{qg;2;3}^{0(11);\Delta_{P}}(z)$
$\displaystyle=$
$\displaystyle\frac{2C_{A}}{\rho}(1-z)+\frac{4C_{1}}{\rho}\Big{(}1-2z^{2}+z(1-z)L_{\rho}\Big{)}\,,$
$\displaystyle P_{qg,2;1}^{0(2);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;1}^{0(2);\Delta_{P}}(z)\,=\,P_{qg,2;2}^{0(12);\nabla}(z)\,=\,2P_{qq,2;2}^{0(12);\Delta_{P}}(z)$
$\displaystyle=$
$\displaystyle-\frac{2C_{A}}{\rho}z-\frac{4C_{1}}{\rho}z\Big{(}4z-(1-2z)L_{\rho}\Big{)}\,,$
$\displaystyle P_{qg,2;1}^{0(22);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;1}^{0(22);\Delta_{P}}(z)\,=\,-\frac{2C_{A}}{\rho}z-\frac{4C_{1}}{\rho}z\Big{(}8z-(1-4z)L_{\rho}\Big{)}\,.$
(S50)

$\displaystyle R_{qg;1}^{1(1);\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{C_{A}-2C_{1}}{2\rho}\,,$ $\displaystyle
R_{qg;2}^{1(1);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;2}^{1(1);\Delta_{P}}-\frac{C_{A}-6C_{1}}{\rho}\,(1+2L_{\rho})$
$\displaystyle=$
$\displaystyle-\frac{C_{A}}{\rho}\,L_{\rho}+\frac{2C_{1}}{\rho}\,(2+5L_{\rho})\,,$
$\displaystyle R_{qg;5}^{1(1);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qq;5}^{1(1);\Delta_{P}}=\frac{2C_{1}}{\rho}\,(1+2L_{\rho})\,,$
$\displaystyle R_{qg;1}^{1(2);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;1}^{1(2);\Delta_{P}}-\frac{2C_{1}}{\rho}\,(3+2L_{\rho})$
$\displaystyle=$
$\displaystyle\frac{2C_{A}}{\rho}+\frac{4C_{1}}{\rho}\,(2+L_{\rho})\,,$
$\displaystyle R_{qg;2}^{1(2);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;2}^{1(2);\Delta_{P}}$ $\displaystyle=$
$\displaystyle-\frac{2C_{A}}{\rho}-\frac{4C_{1}}{\rho}\,(4+L_{\rho})\,,$
$\displaystyle R_{qg;3}^{1(2);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;3}^{1(2);\Delta_{P}}-\frac{C_{A}}{2\rho}\,(2+3L_{\rho})+\frac{3C_{1}}{\rho}\,(4+3L_{\rho})$
$\displaystyle=$
$\displaystyle-\frac{C_{A}}{2\rho}\,L_{\rho}+\frac{3C_{1}}{\rho}\,(2+L_{\rho})\,,$
$\displaystyle R_{qg;4}^{1(2);\nabla}$ $\displaystyle=$ $\displaystyle
2R_{qg;4}^{1(2);\Delta_{P}}-\frac{2C_{1}}{\rho}L_{\rho}\ =\
\frac{2C_{1}}{\rho}\,\Big{(}2+L_{\rho}\Big{)}\,.$ (S51)

$\displaystyle P_{qg,1;1}^{1(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,1;1}^{1(1);\Delta_{P}}(z)+\frac{2z}{\rho}\Big{(}C_{A}\frac{z(2-z)}{(1-z)^{2}_{+,1}}-4C_{1}\frac{1+2z}{(1-z)^{3}_{+,2}}\Big{)}$
$\displaystyle=$
$\displaystyle\frac{z}{\rho}\,\Big{(}C_{A}\frac{z(2-z)}{(1-z)^{2}_{+,1}}-2C_{1}\frac{4+10z+3z^{2}-z^{3}}{(1-z)^{3}_{+,2}}\Big{)}\,,$
$\displaystyle P_{qg,2;2}^{1(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;2}^{1(1);\Delta_{P}}(z)+\frac{z}{\rho}\Big{(}C_{A}\frac{1+3z}{(1-z)_{+}}-6C_{1}\frac{1+z}{(1-z)_{+}}+2C_{1}\,(1+z)\,L_{\rho}\Big{)}$
$\displaystyle=$
$\displaystyle-\frac{1}{\rho}\,\Big{(}C_{A}\frac{1-2z-2z^{2}}{(1-z)_{+}}+2C_{1}\frac{1+5z-7z^{2}+6z^{3}}{(1-z)_{+}}\Big{)}\,,$
$\displaystyle P_{qg,2;3}^{1(1);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;3}^{1(1);\Delta_{P}}(z)+\frac{2C_{A}}{\rho}\,\frac{z^{2}}{(1-z)_{+}}$
$\displaystyle=$
$\displaystyle\frac{2z}{\rho}\,\Big{(}\frac{C_{A}+2C_{1}(1-2z^{2})}{(1-z)_{+}}+2C_{1}\,z\,L_{\rho}\Big{)}\,,$
$\displaystyle P_{qg,2;1}^{1(2);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;1}^{1(1);\Delta_{P}}(z)-\frac{2z}{\rho}C_{1}\,\Big{(}\frac{3+6z-3z^{2}}{(1-z)^{2}_{+,1}}-(1+2z)\,L_{\rho}\Big{)}$
$\displaystyle=$
$\displaystyle-\frac{z}{\rho}\,\Big{(}C_{A}\frac{2-2z+z^{2}}{(1-z)^{2}_{+,1}}+2C_{1}\frac{(2-z)(4-5z+4z^{2})}{(1-z)^{2}_{+,1}}+8C_{1}zL_{\rho}\Big{)}\,,$
$\displaystyle P_{qg,1;2}^{1(2);\Delta_{P}}(z)$ $\displaystyle=$
$\displaystyle-\frac{z^{2}}{2\rho}\,\Big{(}\frac{C_{A}}{(1-z)_{+}}+2C_{1}\frac{1+z}{(1-z)^{2}_{+,1}}\Big{)}\,,$
$\displaystyle P_{qg,2;2}^{1(2);\nabla}(z)$ $\displaystyle=$ $\displaystyle
2P_{qg,2;2}^{1(1);\Delta_{P}}(z)+4C_{1}z^{2}(2-3z)\,L_{\rho}$ (S52)
$\displaystyle=$
$\displaystyle\frac{2z}{\rho}\,\Big{(}C_{A}+2C_{1}\frac{1-6z^{2}+4z^{3}}{(1-z)^{2}_{+,1}}+4C_{1}\,z\,L_{\rho}\Big{)}\,.$