projectors in (4.1). In the rows and columns for $D$ and $27$, we have given
the numerical coefficients for $N=3$. The Collins-Soper kernels in channels
with antiquarks can be obtained from the above expressions by using the
relation (90) and the analogue of (4.3) for ${K}_{a_{1}a_{2}}$.

The matrices $M_{qg}$ and $M_{gg}$ are non-singular for generic values of the
functions $K_{y}$ and $K_{d}$. The eigenvalues of $M_{qg}$ can be calculated
analytically, but the expressions are rather lengthy and we refrain from
giving them here. For $N=3$, the eigenvalues of $M_{gg}$ read

$\displaystyle-3\bigl{(}K_{y}+{\textstyle\frac{1}{2}}K_{d}\bigr{)}\,,$
$\displaystyle-\tfrac{3}{2}K_{y}\,,$
$\displaystyle-\tfrac{3}{2}(K_{y}+K_{d})\,,$
$\displaystyle-(2K_{y}+K_{d})\pm\sqrt{(2K_{y}+K_{d})^{2}+3K_{d}^{2}}\,.$ (235)

####  7.2.1 Limit of small $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$

We now consider the limit where
$|\boldsymbol{z}_{1}|,|\boldsymbol{z}_{2}|\ll|\boldsymbol{y}|$. Expanding
(225), (231) and (233) in powers of $|\boldsymbol{z}_{1}|/|\boldsymbol{y}|$
and $|\boldsymbol{z}_{2}|/|\boldsymbol{y}|$ and comparing the leading terms
with the structure in (171), where the full kernel
${}^{RR^{\prime}}\\!K_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i})$
is split into three separate contributions, we can deduce the LO expressions
of the kernels ${}^{R}\mskip-1.5muK_{a}$ and ${}^{R}J$. We obtain

$\displaystyle{}^{1}\mskip-1.5muK_{q}(\boldsymbol{z};\mu)$
$\displaystyle=C_{F}\mskip 1.5muK(\boldsymbol{z};\mu)\,,$
$\displaystyle{}^{8}\mskip-1.5muK_{q}(\boldsymbol{z};\mu)$
$\displaystyle=-\frac{1}{2N}\mskip 1.5muK(\boldsymbol{z};\mu)\,,$
$\displaystyle{}^{1}\mskip-1.5muK_{g}(\boldsymbol{z};\mu)$
$\displaystyle=C_{A}\mskip 1.5muK(\boldsymbol{z};\mu)\,,$
$\displaystyle{}^{A}\mskip-1.5muK_{g}(\boldsymbol{z};\mu)$
$\displaystyle={}^{S}\mskip-1.5muK_{g}(\boldsymbol{z};\mu)=\frac{N}{2}\mskip
1.5muK(\boldsymbol{z};\mu)$ (236) and $\displaystyle{}^{1}\mskip-1.5muJ$
$\displaystyle=0\,,$
$\displaystyle{}^{8}\mskip-1.5muJ(\boldsymbol{y};\mu_{i})$
$\displaystyle=\frac{N}{2}\bigl{[}K(\boldsymbol{y};\mu_{1})+K(\boldsymbol{y};\mu_{2})\bigr{]}$
(237)

with $K(\boldsymbol{z};\mu)$ defined in (228). For the higher gluon
representations we find

$\displaystyle{}^{D}\mskip-1.5muK_{g}(\boldsymbol{z};\mu)$
$\displaystyle=0\,,$
$\displaystyle{}^{27}\mskip-1.5muK_{g}(\boldsymbol{z};\mu)$
$\displaystyle=-K(\boldsymbol{z};\mu)$ (238) and
$\displaystyle{}^{D}\mskip-1.5muJ(\boldsymbol{z};\mu_{i})$
$\displaystyle=3\bigl{[}K(\boldsymbol{y};\mu_{1})+K(\boldsymbol{y};\mu_{2})\bigr{]}\,,$
$\displaystyle{}^{27}\mskip-1.5muJ(\boldsymbol{z};\mu_{i})$
$\displaystyle=4\bigl{[}K(\boldsymbol{y};\mu_{1})+K(\boldsymbol{y};\mu_{2})\bigr{]}\,,$
(239)

where again we have set $N=3$. The expressions for $J$ are of course only
valid at perturbatively small $\boldsymbol{y}$.

Using the explicit form of $K(\boldsymbol{z};\mu)$, we obtain anomalous
dimensions

$\displaystyle{}^{1}\mskip-1.5mu\gamma_{K,q}$
$\displaystyle=2C_{F}\,\frac{\alpha_{s}}{\pi}\,,$
$\displaystyle{}^{8}\mskip-1.5mu\gamma_{K,q}$
$\displaystyle=-\frac{1}{N}\,\frac{\alpha_{s}}{\pi}\,,$
$\displaystyle{}^{1}\mskip-1.5mu\gamma_{J}$ $\displaystyle=0\,,$
$\displaystyle{}^{8}\mskip-1.5mu\gamma_{J}$
$\displaystyle=N\,\frac{\alpha_{s}}{\pi}\,.$ (240)

The expressions for other colour representations are readily obtained from
(7.2.1) and (237). We note that ${}^{8\,}\mskip-1.5mu\gamma_{K,q}$ and
${}^{27\,}\mskip-1.5mu\gamma_{K,g}$ are negative. Furthermore, we find that
${}^{R\,}\mskip-1.5mu\gamma_{J}>0$ for all $R\neq 1$.

For completeness we also give the anomalous dimensions $\gamma_{a}$, which can
be taken from the literature for single TMDs, e.g. from Echevarria:2016scs
(see our equation (221) for notation). The one-loop expressions are

$\displaystyle\gamma_{q}$ $\displaystyle=\frac{3}{2}\mskip 1.5muC_{F}\mskip
1.5mu\frac{\alpha_{s}}{\pi}\,,$ $\displaystyle\gamma_{g}$
$\displaystyle=\frac{\beta_{0}}{2}\,\frac{\alpha_{s}}{\pi}$ (241)

with $\beta_{0}=\frac{11}{3}\mskip 1.5muC_{A}-\frac{2}{3}\mskip 1.5mun_{F}$.
Finally, one has

$\displaystyle\gamma_{K,q}$ $\displaystyle=2C_{F}\mskip
1.5mu\frac{\alpha_{s}}{\pi}+\frac{C_{F}}{2}\mskip
1.5mu\biggl{(}\frac{67-3\pi^{2}}{9}\mskip 1.5muC_{A}-\frac{10}{9}\mskip
1.5mun_{F}\biggr{)}\biggl{(}\frac{\alpha_{s}}{\pi}\biggr{)}^{2}\,,$
$\displaystyle\gamma_{K,g}$ $\displaystyle=\frac{C_{A}}{C_{F}}\,\gamma_{K,q}$
(242)

for the cusp anomalous dimension at two-loop accuracy. Using the all-order
relation

$\displaystyle{}^{8}\mskip-1.5mu\gamma_{J}$
$\displaystyle=\frac{1}{2}\,\gamma_{K,g}\,,$ (243)

which follows from (185), we readily get the two-loop expression of
${}^{8}\mskip-1.5mu\gamma_{J}$ as well.

###  7.3 Matching coefficients for DTMDs

We now turn to the short-distance coefficients for matching DTMDs on DPDFs in
the large-$\boldsymbol{y}$ regime. In the colour singlet channel, they are
equal to the coefficients for matching TMDs on PDFs, which for most
polarisation combinations have been computed by several groups Collins:2011zzd
; Aybat:2011zv ; Bacchetta:2013pqa ; Echevarria:2015uaa ; Gutierrez-
Reyes:2017glx . In addition to extending these results to colour non-singlet
channels, we have computed the colour singlet coefficients for all
polarisation channels. In the following two subsections, we describe in detail
our method, which allows us to compute the one-loop matching coefficients from
real emission graphs alone and to handle gluon polarisation without using the
Levi-Civita tensor in $D=4-2\epsilon$ dimensions. Subtleties of
renormalisation in both the quark and gluon sector are discussed in section
7.3.3. We give our final results for the colour singlet coefficients in
section 7.3.4 and explain their extension to other colour channels in section
7.3.5.

####  7.3.1 General procedure

To compute the matching coefficients, we take matrix elements of the operator
relations in (160) and (161), using either the vacuum or a single parton as
external states. We define

$\displaystyle
2\pi\delta(p^{+}-p^{\prime+})\,2p^{+}\,{}^{R}\mskip-1.5mu\mathcal{M}_{ab}(x,\boldsymbol{z})$
$\displaystyle=\frac{1}{\mathcal{N}_{b}}\,\frac{1}{m(R)}\,\bigl{\langle}b,p^{\prime},r^{\prime}\big{|}\mskip
1.5mu{}^{R}O^{\mskip
1.5mu\underline{r}}_{a}(x,\boldsymbol{y},\boldsymbol{z})\mskip
1.5mu\big{|}b,p,r\bigr{\rangle}$
$\displaystyle=\varepsilon_{a}(R)\;\frac{\mathcal{N}_{a}}{\mathcal{N}_{b}}\,\frac{P_{R}^{\underline{r}\,\underline{s}}}{m(R)}\,\bigl{\langle}b,p^{\prime},r^{\prime}\big{|}\mskip
1.5muO^{\mskip 1.5mu\underline{s}}_{a}(x,\boldsymbol{y},\boldsymbol{z})\mskip
1.5mu\big{|}b,p,r\bigr{\rangle}\,,\hskip 40.00006pt$
$\displaystyle{}^{R}\mskip-1.5mu\mathcal{M}_{S,a}(\boldsymbol{z})$
$\displaystyle=\frac{P_{R}^{\underline{r}\,\underline{s}}}{m(R)}\,\bigl{\langle}0\big{|}\mskip
1.5muO_{S,a}^{\mskip
1.5mu\underline{r},\underline{s}}(\boldsymbol{y},\boldsymbol{z})\mskip
1.5mu\big{|}\mskip 1.5mu0\bigr{\rangle}\,,$ (244)

where $|\mskip 1.5mub,p,r\rangle$ is a parton state $b$ with momentum $p$ and
colour index $r$, and the colour projected operator ${}^{R}O^{\mskip
1.5mu\underline{r}}_{a}$ was introduced in (79). We recall that
$\underline{r}=(r,r^{\prime})$ and that repeated colour indices are summed
over. We gloss over the issue of parton polarisation for the time being. In
the colour singlet sector, one has $\varepsilon_{a}(1)/m(1)=1$ and the
contraction with
$P_{1}^{\underline{r}\,\underline{s}}\,\mathcal{N}_{a}/\mathcal{N}_{b}$ in
${}^{1}\mskip-1.5mu\mathcal{M}_{ab}$ implies a sum over the colour indices of
the operator and an average over the colour of the parton state. Thus,
${}^{1}\mskip-1.5mu\mathcal{M}_{ab}(x,\boldsymbol{z})$ is the unsubtracted TMD
for parton $a$ in target $b$ according to (3.1). Likewise,
${}^{1}\mskip-1.5mu\mathcal{M}_{S,a}(\boldsymbol{z})$ is the soft factor for
single TMDs. In analogy to (7.3.1) we define matrix elements
${}^{R}\mskip-1.5mu\mathcal{M}_{ab}(x)$ and
${}^{R}\mskip-1.5mu\mathcal{M}_{S,a}$ without argument $\boldsymbol{z}$ from
the operators with $\boldsymbol{z}=\boldsymbol{0}$; thus
${}^{1}\mskip-1.5mu\mathcal{M}_{ab}(x)$ is the PDF for parton $a$ in target
$b$. It is understood that all operators are renormalised, and we recall that
the operators at $\boldsymbol{z}=\boldsymbol{0}$ require additional
renormalisation compared with those at finite $\boldsymbol{z}$.

According to the operator matching equations (161) and (160), we thus have

$\displaystyle{}^{R}\mskip-1.5mu\mathcal{M}_{ac}(x,\boldsymbol{z})$
$\displaystyle=\sum_{b}{}^{R}C_{\text{us},ab}(x^{\prime},\boldsymbol{z})\,\underset{x}{\otimes}{}^{R}\mskip-1.5mu\mathcal{M}_{bc}(x^{\prime})\,,$
$\displaystyle{}^{R}\mskip-1.5mu\mathcal{M}_{S,a}(\boldsymbol{z})$
$\displaystyle={}^{R}C_{S,a}(\boldsymbol{z})\,{}^{R}\mskip-1.5mu\mathcal{M}_{S,a}\,.$
(245)

We now expand the quantities in these relations up to
$\mathcal{O}(\alpha_{s})$, denoting the zeroth and first orders with
superscripts $(0)$ and $(1)$, respectively. At tree level, the matrix elements
$\langle b,p^{\prime},r^{\prime}\mskip 1.5mu|\mskip 1.5muO^{\mskip
1.5mu\underline{s}}_{a}\mskip 1.5mu|b,p,r\rangle$ and $\langle 0|\mskip
1.5muO_{S,a}^{\mskip 1.5mu\underline{r},\underline{s}}\mskip 1.5mu|0\rangle$
have the colour structure $\delta_{rs}\,\delta_{r^{\prime}s^{\prime}}$, which
gives

$\displaystyle{}^{R}\mskip-1.5mu\mathcal{M}^{(0)}_{ab}(x,\boldsymbol{z})$
$\displaystyle={}^{R}\mskip-1.5mu\mathcal{M}^{(0)}_{ab}(x)=\delta_{ab}\,\delta(1-x)\,,$
$\displaystyle{}^{R}\mskip-1.5mu\mathcal{M}^{(0)}_{S,a}(\boldsymbol{z})$
$\displaystyle={}^{R}\mskip-1.5mu\mathcal{M}^{(0)}_{S,a}=1$ (246)

for the matrix elements on both sides of (7.3.1), which results in zeroth-
order matching coefficients

$\displaystyle{}^{R}C_{\text{us},ab}^{(0)}(x,\boldsymbol{z})$
$\displaystyle=\delta_{ab}\,\delta(1-x)\,,$
$\displaystyle{}^{R}C_{S,a}^{(0)}(\boldsymbol{z})$ $\displaystyle=1\,.$ (247)

Inserting this into the $\mathcal{O}(\alpha_{s})$ terms of (7.3.1), we obtain

$\displaystyle{}^{R}C_{\text{us},ab}^{(1)}(x,\boldsymbol{z})$
$\displaystyle={}^{R}\mskip-1.5mu\mathcal{M}^{(1)}_{ab}(x,\boldsymbol{z})-{}^{R}\mskip-1.5mu\mathcal{M}^{(1)}_{ab}(x)\,,$
$\displaystyle{}^{R}C_{S,a}^{(1)}(\boldsymbol{z})$
$\displaystyle={}^{R}\mskip-1.5mu\mathcal{M}^{(1)}_{S,a}(\boldsymbol{z})-{}^{R}\mskip-1.5mu\mathcal{M}^{(1)}_{S,a}\,.$
(248)

The one-loop matching coefficients are thus given as the difference of
$\mathcal{O}(\alpha_{s})$ matrix elements at finite and at zero
$\boldsymbol{z}$. Virtual graphs cancel in this difference, because they do
not depend on $\boldsymbol{z}$, so we can limit our calculation to real
graphs.

To obtain the matching coefficient ${}^{R}C_{ab}$ of the TMD, we need the
square root of ${}^{R}C_{S,a}$ according to (165). At
$\mathcal{O}(\alpha_{s})$, we have

$\displaystyle\sqrt{{}^{R}C_{S,a}^{(1)}(\boldsymbol{z};2Y)}$
$\displaystyle=1+\frac{1}{2}\,{}^{R}C_{S,a}^{(1)}(\boldsymbol{z};2Y)=1+{}^{R}C_{S,a}^{(1)}(\boldsymbol{z};Y)\,.$
(249)

In the second step, we have used that ${}^{R}C_{S,a}^{(1)}(\boldsymbol{z};Y)$
is linear in $Y$, which readily follows from its evolution equation (169) and
the fact that ${}^{R}\mskip-1.5muK_{a}=\mathcal{O}(\alpha_{s})$. Inserting
this into (165), we thus obtain at $\mathcal{O}(\alpha_{s})$

$\displaystyle{}^{R}C_{ab}(x,\boldsymbol{z};\zeta)=\delta_{ab}\,\delta(1-x)\,$
$\displaystyle\qquad+\lim_{Y_{L}\to-\infty}\Bigl{[}\mskip
1.5mu{}^{R}C^{(1)}_{\text{us},ab}(x,\boldsymbol{z};Y_{L})-\delta_{ab}\,\delta(1-x)\,{}^{R}C^{(1)}_{S,a}(\boldsymbol{z};Y_{C}-Y_{L})\mskip
1.5mu\Bigr{]}\,,$ (250)

where in the present section we use the definition (61) of $\zeta$, which
refers to a single parton with plus-momentum $xp^{+}$. We have omitted the
argument $\mu$ in all functions for brevity.

####  7.3.2 Calculation of gluon-gluon matching coefficients

We now describe in detail our procedure to compute the one-loop matching
coefficients in the pure gluon channel. Throughout this subsection, we
consider the colour singlet sector, $R=1$, and drop the corresponding index.
In turn, we now make gluon polarisation explicit. We take gluon operators
$O_{g}^{jj^{\prime}}(x,\boldsymbol{y},\boldsymbol{z})$ with open polarisation
indices, obtained from (21) and (17) by dropping the projection operator
$\Pi^{jj^{\prime}}_{a}$ in the latter equation. We also take open polarisation
indices $i$, $i^{\prime}$ for the gluons in the corresponding matrix element
(7.3.1), which thus carries four Lorentz indices,
$\mathcal{M}^{jj^{\prime}\\!,ii^{\prime}}_{gg}$. All of them are restricted to
be in the $D-2$ transverse dimensions. The subscripts $g$ refer of course only
to the parton type here, and not to its polarisation. For the gluon matrix
elements, we thus have

$\displaystyle\mathcal{M}^{jj^{\prime}\\!,ll^{\prime}}_{gg}(x,\boldsymbol{z})$
$\displaystyle=C^{jj^{\prime}\\!,ii^{\prime}}_{\text{us},gg}(x^{\prime},\boldsymbol{z})\,\underset{x}{\otimes}\mathcal{M}^{ii^{\prime}\\!,ll^{\prime}}_{gg}(x^{\prime})+\text{\\{quark-
gluon mixing terms\\}}$ (251)

and

$\displaystyle\mathcal{M}^{(0)\,jj^{\prime}\\!,ii}_{gg}(x,\boldsymbol{z})$
$\displaystyle=\mathcal{M}^{(0)\,jj^{\prime}\\!,ii^{\prime}}_{gg}(x)=C^{(0)\,jj^{\prime}\\!,ii}_{gg}(x,\boldsymbol{z})=\delta^{ij}\,\delta^{i^{\prime}j^{\prime}}\delta(1-x)\,,$
(252)

which is combined with the matching of the soft factor into

$\displaystyle C^{(1)\,jj^{\prime}\\!,ii}_{gg}(x,\boldsymbol{z};\zeta)$
$\displaystyle=\lim_{Y_{L}\to-\infty}\Bigl{[}\mskip
1.5muC^{(1)\,jj^{\prime}\\!,ii}_{\text{us},gg}(x,\boldsymbol{z};Y_{L})-\delta^{ij}\,\delta^{i^{\prime}j^{\prime}}\delta(1-x)\,C^{(1)}_{S,a}(\boldsymbol{z};Y_{C}-Y_{L})\mskip
1.5mu\Bigr{]}\,.$ (253)

The virtual graphs contributing to the one-loop matrix elements are shown in
figure 11. They are independent of $\boldsymbol{z}$ (the Fourier conjugated
distance to the transverse momentum carried by partons or eikonal lines) and
hence cancel in the matching coefficient. Therefore we only compute the real
graphs in figure 11. We omit so-called Wilson line self interactions, where a
gluon is exchanged between eikonal lines along the same direction $v$,
referring to chapter 13.7 of Collins:2011zzd  for further discussion.

Figure 10:  Virtual graphs contributing to the matrix elements
$\mathcal{M}^{(1)\,jj^{\prime}\\!,ii}_{gg}$ (a and b) and
$\mathcal{M}_{S,g}^{(1)}$ (c). Not shown are complex conjugate graphs and the
analogue of (a) with a quark instead of a gluon loop. All eikonal lines are in
the adjoint representation.

Figure 11:  Real graphs contributing to the matrix elements
$\mathcal{M}^{(1)\,jj^{\prime}\\!,ii}_{gg}$ (a and b) and
$\mathcal{M}_{S,g}^{(1)}$ (c). Not shown are the complex conjugates of graphs
(b) and (c).

We work in Feynman gauge and compute cut graphs (where the emitted gluon is
explicitly put on shell), using the Feynman rules given in appendix D. For
graph 11b, we then obtain

$\displaystyle
C_{\text{us},gg}^{(1)\,jj^{\prime},\,ii^{\prime}}(x,\boldsymbol{z};Y_{L}\to-\infty)\mskip
1.5mu\big{|}_{\text{fig.\leavevmode\nobreak\ \ref{fig:Cgg-
real}b}}=\frac{1}{xp^{+}}\int\\!\frac{d^{4-2\epsilon}k}{(2\pi)^{4-2\epsilon}}\,\bigl{(}e^{i\boldsymbol{k}\boldsymbol{z}}-1\bigr{)}\;2\pi\delta\big{(}(p-k)^{2}\big{)}\,\delta(k^{+}-xp^{+})\,$
$\displaystyle\qquad\times\frac{\delta_{aa^{\prime}}\delta_{de^{\prime}}}{N^{2}-1}\;\bigl{(}g\mu^{\epsilon}n^{\nu}f^{dbc}\bigr{)}\;\bigl{(}-g_{\nu\nu^{\prime}}\mskip
1.5mu\delta_{bb^{\prime}}\bigr{)}\;\frac{i\delta_{c^{\prime}d^{\prime}}\mskip
1.5mug_{\rho^{\prime}\tau^{\prime}}}{k^{2}-i0}\,\frac{i\delta_{ce}}{(p-k)\mskip
1.5mun-i0}$
$\displaystyle\qquad\times\Bigl{(}g\mu^{\epsilon}f^{a^{\prime}b^{\prime}c^{\prime}}\,\bigl{[}(p-2k)^{i^{\prime}}g^{\nu^{\prime}\rho^{\prime}}+(k+p)^{\nu^{\prime}}g^{i^{\prime}\rho^{\prime}}+(k-2p)^{\rho^{\prime}}g^{i^{\prime}\nu^{\prime}}\bigr{]}\Bigr{)}$
$\displaystyle\qquad\times\Bigl{(}-i\bigl{[}kn\,g^{ij}-p^{j}n^{i}\bigr{]}\,\delta_{ae}\Bigr{)}\,\Bigl{(}i\mskip
1.5mu\bigl{[}kn\,g^{\tau^{\prime}j^{\prime}}-k^{j^{\prime}}n^{\tau^{\prime}}\bigr{]}\,\delta_{d^{\prime}e^{\prime}}\Bigr{)}+C_{\text{us},\epsilon}\mskip
1.5mu\big{|}_{\text{fig.\leavevmode\nobreak\ \ref{fig:Cgg-real}b}}\,,$ (254)

where $n$ is the four-vector with $n^{-}=1$, $n^{+}=0$ and
$\boldsymbol{n}=\boldsymbol{0}$. The assignment of polarisation indices is
shown in figure 11a. All primed indices refer to the right of the final-state
cut, and $a,b,\ldots,e$ are colour octet indices. The factor
$e^{i\boldsymbol{k}\boldsymbol{z}}-1$ represents the difference between the
matrix element for nonzero and zero $\boldsymbol{z}$. In accordance with our
choice of light-cone coordinates, the components $p^{-}$ and $\boldsymbol{p}$
of the target momentum are zero. The factor $(xp^{+})^{-1}$ comes from the
definition (21) of the operator $O_{g}$. The ultraviolet counterterm
$C_{\text{us},\epsilon}$ will be discussed below.

In (254), we have taken the limit $Y_{L}\to-\infty$ by setting $v_{L}=n$ in
the eikonal propagator. (Only then can we use the Feynman rule given in figure
14). This gives a factor $p^{+}-k^{+}=p^{+}(1-x)$ in the denominator. The
resulting rapidity divergences in the convolution of $C_{\text{us},gg}$ with a
PDF or a DPDF are cancelled by corresponding divergences in the soft factor,
as we shall see shortly. In this sense, one should understand the limit
$Y_{L}\to-\infty$ in $C_{\text{us},gg}$ as being taken in the combination
(253).

Adding the complex conjugate of (254) and the expression for graph 11a, we get

$\displaystyle
C_{\text{us},gg}^{(1)\,jj^{\prime}\\!,ii^{\prime}}(x,\boldsymbol{z};Y_{L}\to-\infty)=\frac{\alpha_{s}\mskip
1.5muC_{A}}{\pi^{2}}\,(2\pi\mu)^{2\epsilon}\int\frac{d^{2-2\epsilon}\boldsymbol{k}}{\boldsymbol{k}^{2}}\,\bigl{(}e^{i\boldsymbol{k}\boldsymbol{z}}-1\bigr{)}$
$\displaystyle\qquad\times\Bigg{\\{}\bigg{[}\frac{x}{1-x}+\frac{(1-x)(1+x^{2})}{x}\bigg{]}\,\frac{1}{2(1-\epsilon)}\,\delta^{ii^{\prime}}\delta^{jj^{\prime}}$
$\displaystyle\qquad\quad+\biggl{[}\frac{x}{1-x}+\frac{2(1-x)}{1-\epsilon}\bigg{]}\,\frac{1}{2}\mskip
1.5mu\Big{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}-\delta^{ij^{\prime}}\delta^{i^{\prime}j}\Big{)}$
$\displaystyle\qquad\quad+\frac{x}{1-x}\;\frac{1}{2}\mskip
1.5mu\biggl{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}+\delta^{ij^{\prime}}\delta^{i^{\prime}j}-\frac{1}{1-\epsilon}\,\delta^{ii^{\prime}}\delta^{jj^{\prime}}\biggr{)}$
$\displaystyle\qquad\quad+x(1-x)\,\frac{k^{ii^{\prime}}}{\boldsymbol{k}^{2}}\,\delta^{jj^{\prime}}+\frac{1-x}{x}\,\delta^{ii^{\prime}}\,\frac{k^{jj^{\prime}}}{\boldsymbol{k}^{2}}$
$\displaystyle\qquad\quad+(1-x)\,\biggl{(}\delta^{ij}\frac{k^{i^{\prime}j^{\prime}}}{\boldsymbol{k}^{2}}+\delta^{i^{\prime}j^{\prime}}\,\frac{k^{ij}}{\boldsymbol{k}^{2}}-\delta^{ij^{\prime}}\,\frac{k^{i^{\prime}j}}{\boldsymbol{k}^{2}}-\delta^{i^{\prime}j}\,\frac{k^{ij^{\prime}}}{\boldsymbol{k}^{2}}\biggr{)}\Bigg{\\}}+C_{\text{us},\epsilon}\,,$
(255)

where we define the symmetric and traceless tensor

$\displaystyle v^{ij}$
$\displaystyle=\boldsymbol{v}^{i}\boldsymbol{v}^{j}-\frac{\boldsymbol{v}^{2}}{2(1-\epsilon)}\;\delta^{ij}$
(256)

for a given vector $\boldsymbol{v}$ in $2(1-\epsilon)$ transverse dimensions.
Using the integral relations in appendix E, we perform the transverse
integration and obtain

$\displaystyle
C_{\text{us},gg}^{(1)\,jj^{\prime}\\!,ii^{\prime}}(x,\boldsymbol{z};-\infty)=-\frac{\alpha_{s}\mskip
1.5muC_{A}}{\pi}\,(2\pi\mu)^{2\epsilon}\;\Gamma(1-\epsilon)\,\biggl{(}\frac{\boldsymbol{z}^{2}}{4\pi}\biggr{)}^{\epsilon}$
$\displaystyle\qquad\times\Bigg{\\{}\,\frac{1}{\epsilon}\,\biggl{[}\frac{x}{1-x}+\frac{(1-x)(1+x^{2})}{x}\biggr{]}\,\frac{1}{2(1-\epsilon)}\,\delta^{ii^{\prime}}\delta^{jj^{\prime}}$
$\displaystyle\qquad\quad+\frac{1}{\epsilon}\,\biggl{[}\frac{x}{1-x}+\frac{2(1-x)}{1-\epsilon}\biggr{]}\,\frac{1}{2}\,\Big{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}-\delta^{ij^{\prime}}\delta^{i^{\prime}j}\Big{)}$
$\displaystyle\qquad\quad+\frac{1}{\epsilon}\,\frac{x}{1-x}\;\frac{1}{2}\mskip
1.5mu\biggl{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}+\delta^{ij^{\prime}}\delta^{i^{\prime}j}-\frac{1}{1-\epsilon}\,\delta^{ii^{\prime}}\delta^{jj^{\prime}}\biggr{)}$
$\displaystyle\qquad\quad+x(1-x)\,\frac{z^{ii^{\prime}}}{\boldsymbol{z}^{2}}\,\delta^{jj^{\prime}}+\frac{1-x}{x}\,\delta^{ii^{\prime}}\,\frac{z^{jj^{\prime}}}{\boldsymbol{z}^{2}}$
$\displaystyle\qquad\quad+(1-x)\,\biggl{(}\delta^{ij}\,\frac{z^{i^{\prime}j^{\prime}}}{\boldsymbol{z}^{2}}+\delta^{i^{\prime}j^{\prime}}\,\frac{z^{ij}}{\boldsymbol{z}^{2}}-\delta^{ij^{\prime}}\,\frac{z^{i^{\prime}j}}{\boldsymbol{z}^{2}}-\delta^{i^{\prime}j}\,\frac{z^{ij^{\prime}}}{\boldsymbol{z}^{2}}\biggr{)}\Bigg{\\}}+C_{\text{us},\epsilon}\,.$
(257)

In (255) and (257) we have organised the tensor structure into terms that are
diagonal or antisymmetric or symmetric and traceless in the index pair
$jj^{\prime}$, for reasons that will become clear in section 7.3.3. We note
that the tensor in parentheses in the last line of (257) is nonzero in a
generic number $2-2\epsilon$ of transverse dimensions but vanishes for
$\epsilon=0$. This is easily seen by contracting this tensor with itself,
which gives an expression proportional to $\epsilon$.

We see that (257) has poles in $1/\epsilon$. They correspond to ultraviolet
divergences of the matrix element at $\boldsymbol{z}=\boldsymbol{0}$ and are
removed by the counterterm $C_{\text{us},\epsilon}$. If one splits the
$\boldsymbol{k}$ integration in (255) into the parts going with
$e^{i\boldsymbol{k}\boldsymbol{z}}$ and $-1$, then the $1/\epsilon$ poles
arise as infrared divergences of the term with
$e^{i\boldsymbol{k}\boldsymbol{z}}$, which is ultraviolet finite. The term
with $-1$ involves scaleless integrals, in which the ultraviolet and infrared
divergences add up to zero in dimensional regularisation. In the sum of terms,
the infrared divergences cancel and the ultraviolet ones remain. The
$\overline{\text{MS}}$ counterterm reads

$\displaystyle C_{\text{us},\epsilon}=\frac{\alpha_{s}\mskip
1.5muC_{A}}{\pi}\,\frac{S_{\epsilon}}{\epsilon}\;$
$\displaystyle\Bigg{\\{}\bigg{[}\frac{x}{1-x}+\frac{(1-x)(1+x^{2})}{x}\bigg{]}\,\frac{1}{2(1-\epsilon)}\,\delta^{ii^{\prime}}\delta^{jj^{\prime}}$
$\displaystyle+\bigg{[}\frac{x}{1-x}+2(1-x)\bigg{]}\,\frac{1}{2}\,\Big{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}-\delta^{ij^{\prime}}\delta^{i^{\prime}j}\Big{)}$
$\displaystyle+\frac{x}{1-x}\;\frac{1}{2}\,\biggl{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}+\delta^{ij^{\prime}}\delta^{i^{\prime}j}-\frac{1}{1-\epsilon}\,\delta^{ii^{\prime}}\delta^{jj^{\prime}}\biggr{)}\Bigg{\\}}\,,$
(258)

where

$\displaystyle S_{\epsilon}=\bigl{(}4\pi
e^{-\gamma_{E}}\bigr{)}^{\epsilon}\,.$ (259)

The origin of the factor $1/(1-\epsilon)$ multiplying
$\delta^{ii^{\prime}}\delta^{jj^{\prime}}$ is explained in section 7.3.3; note
that the naive implementation of $\overline{\text{MS}}$ counterterms as “the
residue of $1/\epsilon$ times $S_{\epsilon}$” would lead to an incorrect
result.

We now turn to the matching coefficient of the soft factor, which comes from
the graph in figure 11c and its complex conjugate. Taking again the limit
$Y_{L}\to-\infty$, we obtain

$\displaystyle
C_{S,g}^{(1)}(\boldsymbol{z};Y_{C}+\infty)=\int\frac{d^{4-2\epsilon}\ell}{(2\pi)^{4-2\epsilon}}\;\bigl{(}e^{-i\boldsymbol{\ell}\boldsymbol{z}}-1\bigr{)}\;2\pi\delta\bigl{(}\ell^{2}\bigr{)}\,\theta\bigl{(}\ell^{+}\bigr{)}\;\frac{\delta_{aa^{\prime}}\delta_{cc^{\prime}}}{N^{2}-1}$
$\displaystyle\qquad\times\frac{i}{\ell\mskip
1.5muv_{C}+i0}\;\bigl{(}-g\mu^{\epsilon}v_{C}^{\nu^{\prime}}f^{a^{\prime}b^{\prime}c^{\prime}}\bigr{)}\;\bigl{(}-g_{\nu\nu^{\prime}}\mskip
1.5mu\delta_{bb^{\prime}}\bigr{)}\;\bigl{(}g\mu^{\epsilon}n^{\nu}f^{cba}\bigr{)}\;\frac{i}{\ell\mskip
1.5mun-i0}+\text{c.c.}+C_{S,\epsilon}$
$\displaystyle\quad=\frac{\alpha_{s}\mskip
1.5muC_{A}}{4\pi^{2}}\,(2\pi\mu)^{2\epsilon}\int
d^{2-2\epsilon}\boldsymbol{\ell}\;\bigl{(}e^{-i\boldsymbol{\ell}\boldsymbol{z}}-1\bigr{)}\,\Bigg{[}\int_{0}^{\infty}\frac{d\ell^{+}}{\ell^{+}}\,\frac{2v_{C}^{+}}{2(\ell^{+})^{2}\,v_{C}^{-}+\boldsymbol{\ell}^{2}\,v_{C}^{+}-i0}$
$\displaystyle\hskip
180.00027pt-\frac{2}{\boldsymbol{\ell}^{2}}\int_{0}^{xp^{+}}\frac{d\ell^{+}}{\ell^{+}}+\frac{2}{\boldsymbol{\ell}^{2}}\int_{0}^{xp^{+}}\frac{d\ell^{+}}{\ell^{+}}\mskip
1.5mu\Bigg{]}+\text{c.c.}+C_{S,\epsilon}$
$\displaystyle\quad=-\frac{\alpha_{s}\mskip
1.5muC_{A}}{2\pi^{2}}\,(2\pi\mu)^{2\epsilon}\int\frac{d^{2-2\epsilon}\boldsymbol{\ell}}{\boldsymbol{\ell}^{2}}\;\bigl{(}e^{-i\boldsymbol{\ell}\boldsymbol{z}}-1\bigr{)}\,\Bigg{[}\log\frac{2(xp^{+})^{2}\,e^{-2Y_{C}}}{\boldsymbol{\ell}^{2}}-2\int_{0}^{xp^{+}}\frac{d\ell^{+}}{\ell^{+}}\mskip
1.5mu\Bigg{]}+C_{S,\epsilon}\,,$ (260)

where $C_{S,\epsilon}$ is the ultraviolet counterterm and “c.c.” denotes the
complex conjugate term (which results in an extra factor of $2$ in the last
line). The eikonal propagator depending on $n$ gives rise to the divergent
integral over $1/\ell^{+}$ in the third line. To combine this efficiently with
the divergence from $C_{\text{us},gg}$, we have subtracted and added a
divergent integral in the fourth line. The first two terms in square brackets
then combine to a finite integral, which gives the logarithm in the last line.
Having rewritten $v_{C}^{-}/v_{C}^{+}=-e^{-2Y_{C}}$, we recognise the variable
$\zeta$ from (61) appearing in this logarithm.

In the remaining divergent integral, we now change variables as
$\ell^{+}=p^{+}(1-x^{\prime})$, motivated by the corresponding form of the
gluon plus-momentum $p^{+}-k^{+}=p^{+}(1-x)$ crossing the cut in figure 11b.
Performing the Fourier transform using (344) and (345), we then obtain

$\displaystyle C_{S,g}^{(1)}(\boldsymbol{z};Y_{C}+\infty)$
$\displaystyle=-\frac{\alpha_{s}\mskip
1.5muC_{A}}{2\pi}\,(4\pi)^{\epsilon}\,\biggl{(}\frac{\mu^{2}\boldsymbol{z}^{2}}{4}\biggr{)}^{\epsilon}\;\frac{\Gamma(1-\epsilon)}{\epsilon}$
$\displaystyle\quad\times\Bigg{[}\log\frac{4}{\zeta\boldsymbol{z}^{2}}-\gamma_{E}+\psi(-\epsilon)+2\int_{1-x}^{1}\frac{dx^{\prime}}{1-x^{\prime}}\Bigg{]}+C_{S,\epsilon}\,,$
(261)

where $\psi(z)=\frac{d}{dz}\log\Gamma(z)$ is the digamma function. Its
expansion around $\epsilon=0$ gives $\psi(-\epsilon)=1/\epsilon+\ldots\,$, so
that overall we obtain a double pole in $\epsilon$. We recall that for a one-
loop quantity that has the form

$\displaystyle\frac{R_{-2}}{\epsilon^{2}}+\frac{R_{-1}}{\epsilon}+\mathcal{O}(\epsilon^{0})\,,$
(262)

the $\overline{\text{MS}}$-counterterm is given by

$\displaystyle-
S_{\epsilon}\left(\frac{R_{-2}}{\epsilon^{2}}+\frac{R_{-1}-S_{1}\mskip
1.5muR_{-2}}{\epsilon}\right)$ (263)

times the tree-level term $C_{S,g}^{(0)}=1$, where the coefficient $S_{1}$
comes from the Taylor expansion
$S_{\epsilon}=(4\pi\exp(-\gamma_{E}))^{\epsilon}=1+S_{1}\mskip
1.5mu\epsilon+\mathcal{O}(\epsilon^{2})$. The additional contribution
$S_{\epsilon}\mskip 1.5muS_{1}\mskip 1.5muR_{-2}/\epsilon$ compensates the
$1/\epsilon$ term that arises from Taylor expanding $S_{\epsilon}\mskip
1.5muR_{-2}/\epsilon^{2}$. Implementing this prescription, we obtain the final
form for the soft matching coefficient:

$\displaystyle C_{S,g}^{(1)}(\boldsymbol{z};Y_{C}+\infty)$
$\displaystyle=-\frac{\alpha_{s}\mskip
1.5muC_{A}}{2\pi}\,\bigg{[}-\frac{1}{2}L^{2}+L\,\log\frac{\mu^{2}}{\zeta}-\frac{\pi^{2}}{12}+2L\int_{1-x}^{1}\frac{dx^{\prime}}{1-x^{\prime}}\bigg{]}$
(264)

with

$\displaystyle L$
$\displaystyle=\log\frac{\mu^{2}\boldsymbol{z}^{2}}{b_{0}^{2}}\,,$ (265)

where $b_{0}$ is given in (180). We note in passing that $C_{S,q}$ is given by
the same expression, with the colour factor $C_{A}$ replaced by $C_{F}$.
Defining

$\displaystyle S_{L}$
$\displaystyle=-\frac{1}{2}L^{2}+L\,\log\frac{\mu^{2}}{\zeta}-\frac{\pi^{2}}{12}$
(266)

and combining $C_{\text{us},gg}$ with $C_{S,g}$ according to (253), we obtain
the full matching coefficient in tensor form:

$\displaystyle
C_{gg}^{(1)\,jj^{\prime}\\!,ii^{\prime}}(x,\boldsymbol{z};\mu,\zeta)$
$\displaystyle\quad\\!=\frac{\alpha_{s}\mskip
1.5muC_{A}}{2\pi}\,\biggl{\\{}-2L\,\bigg{[}\mskip
1.5mu\frac{x}{(1-x)_{+}}+\frac{(1-x)(1+x^{2})}{x}\bigg{]}+S_{L}\,\delta(1-x)\mskip
1.5mu\biggr{\\}}\,\frac{1}{2}\,\delta^{ii^{\prime}}\delta^{jj^{\prime}}$
$\displaystyle\qquad\\!+\frac{\alpha_{s}\mskip
1.5muC_{A}}{2\pi}\,\biggl{\\{}-2L\,\bigg{[}\mskip
1.5mu\frac{x}{(1-x)_{+}}+2(1-x)\bigg{]}-4(1-x)+S_{L}\,\delta(1-x)\mskip
1.5mu\biggr{\\}}\,\frac{1}{2}\,\Big{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}-\delta^{ij^{\prime}}\delta^{i^{\prime}j}\Big{)}$
$\displaystyle\qquad\\!+\frac{\alpha_{s}\mskip
1.5muC_{A}}{2\pi}\,\biggl{\\{}-2L\;\frac{x}{(1-x)_{+}}+S_{L}\,\delta(1-x)\mskip
1.5mu\biggr{\\}}\,\frac{1}{2}\,\Bigl{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}+\delta^{ij^{\prime}}\delta^{i^{\prime}j}-\delta^{ii^{\prime}}\delta^{jj^{\prime}}\Bigr{)}$
$\displaystyle\qquad\\!-\frac{\alpha_{s}\mskip
1.5muC_{A}}{2\pi}\;2x(1-x)\,\frac{z^{ii^{\prime}}}{\boldsymbol{z}^{2}}\,\delta^{jj^{\prime}}-\frac{\alpha_{s}\mskip
1.5muC_{A}}{2\pi}\;\frac{2(1-x)}{x}\,\delta^{ii^{\prime}}\,\frac{z^{jj^{\prime}}}{\boldsymbol{z}^{2}}\,.$
(267)

The rapidity divergences in the $C_{\text{us},gg}$ and $C_{S,g}$ have been
combined into the plus-distribution

$\displaystyle\frac{1}{(1-x)_{+}}$
$\displaystyle=\frac{1}{1-x}-\delta(1-x)\int_{0}^{1}\frac{dx^{\prime}}{1-x^{\prime}}\,,$
(268)

which gives a finite result in convolution integrals. The lower integration
limit of the $x^{\prime}$ integral has changed from $1-x$ in (264) to $0$ in
(268) because $C_{S,g}$ is multiplied with $\delta(1-x)$ in the combination
formula (253).

Note that the method we have used avoids introducing the antisymmetric tensor
$\epsilon_{\lambda\mu\nu\rho}$ in $D=4-2\epsilon$ dimensions. With the final
result (267) in the $D=4$ physical dimensions, it is straightforward to
perform a projection onto unpolarised or longitudinally polarised gluons. On
the operator side, one uses the spin projectors $\delta^{jj^{\prime}}$ or
$i\epsilon^{jj^{\prime}}$ from (18), and for the gluon states in the matrix
element, one contracts with $\delta^{ii^{\prime}}/2$ or
$-i\epsilon^{ii^{\prime}}/2$, which gives the average or half the difference
of the two helicity states, respectively.141414The sign difference for
longitudinal polarisation reflects that one has polarisation vectors
$\varepsilon$ for incoming gluons and $\varepsilon^{*}$ for outgoing ones. For
linear gluon polarisation, we see that $\tau^{jj^{\prime}\\!,kk^{\prime}}$ in
(18) projects on the symmetric and traceless part in the indices
$jj^{\prime}$, with the result depending on two transverse indices
$kk^{\prime}$. Matching coefficients involving linear gluon polarisation thus
carry indices, as already remarked earlier.

We will give the results of projecting (267) onto definite gluon polarisations
after discussing in more detail the renormalisation counterterms we have used
in our calculations.

####  7.3.3 Subtleties of renormalisation

As explained earlier, the renormalisation of TMDs involves only virtual
graphs, which drop out in the calculation of the one-loop matching
coefficients. It thus remains to discuss the renormalisation of the matrix
elements $\mathcal{M}^{jj^{\prime}\\!,ll^{\prime}}_{ab}(x)$ and
$\mathcal{M}_{S,a}$, which are respectively related with PDFs and with the
soft factor at $\boldsymbol{z}=\boldsymbol{0}$.

##### Gluon polarisation.

We begin with renormalisation of gluon PDFs, which gives the counterterm (258)
for the real graphs in figure 11. We recall that the renormalisation of a
scalar PDF operator (e.g. the one for unpolarised gluons) reads

$\displaystyle{O}(x)=Z(x^{\prime})\underset{x}{\otimes}{O}_{B}(x^{\prime})$
(269)

and gives

$\displaystyle\mathcal{M}^{(1)}(x)$
$\displaystyle=Z^{(0)}(x^{\prime})\underset{x}{\otimes}\mathcal{M}_{B}^{(1)}(x^{\prime})+Z^{(1)}(x^{\prime})\underset{x}{\otimes}\mathcal{M}_{B}^{(0)}(x^{\prime})\,,$
(270)

where $Z^{(0)}(x)=\delta(1-x)$ and the first-order term $Z^{(1)}(x)$ is given
by the $\overline{\text{MS}}$ prescription. Bare quantities are denoted with a
subscript $B$, and we ignore mixing with quark operators for the time being.
For the gluon operators $O^{jj^{\prime}}(x)$ with open indices, which we are
now interested in, the renormalisation factor $Z$ in (269) turns into a tensor
with four indices. As explained in section 6.5 of Collins:1984xc ,
renormalisation is significantly simplified by first splitting the tensor
operator into its antisymmetric, traceless symmetric diagonal terms, which in
$2-2\epsilon$ transverse dimensions reads

$\displaystyle O^{jj^{\prime}}$
$\displaystyle=\frac{1}{2}\left(O^{jj^{\prime}}-O^{j^{\prime}j}\right)+\frac{1}{2}\left(O^{jj^{\prime}}+O^{j^{\prime}j}+\frac{1}{1-\epsilon}\,\delta^{jj^{\prime}}O^{kk}\right)-\frac{1}{2(1-\epsilon)}\,\delta^{jj^{\prime}}O^{kk}\,.$
(271)

The three terms of this decomposition correspond to different irreducible
representations of the rotation group. Therefore, each of them is renormalised
by a scalar renormalisation factor, and there is no mixing between them. To
implement (269) in the tensor case, we project all tensors in $jj^{\prime}$ in
the same way, which in particular gives

$\displaystyle\delta^{ij}\delta^{i^{\prime}j^{\prime}}$
$\displaystyle=\frac{1}{2}\left(\delta^{ij}\delta^{i^{\prime}j^{\prime}}-\delta^{ij^{\prime}}\delta^{i^{\prime}j}\right)+\frac{1}{2}\left(\delta^{ij}\delta^{i^{\prime}j^{\prime}}+\delta^{ij^{\prime}}\delta^{i^{\prime}j}-\frac{1}{1-\epsilon}\,\delta^{ii^{\prime}}\delta^{jj^{\prime}}\right)+\frac{1}{2(1-\epsilon)}\,\delta^{ii^{\prime}}\delta^{jj^{\prime}}$
(272)

for the tree-level tensors in (252). This explains the appearance of the
factors $1/(1-\epsilon)$ in the counterterm (258). We note that the
coefficients of the diagonal, antisymmetric and traceless symmetric tensors in
that term are proportional to the well-known DGLAP kernels for the associated
gluon polarisation (apart from the $\delta(1-x)$ terms from virtual graphs),
which confirms the consistency of the procedure. Contracting the last term in
(272) with the projector $\delta^{jj^{\prime}}$ for unpolarised gluons and
taking the average over incoming gluon polarisations by contracting with
$\delta^{ii^{\prime}}\big{/}\bigl{(}2(1-\epsilon)\bigr{)}$, one obtains unity,
which corresponds to the correct normalisation of the tree-level matrix
element.

Notice also that the three terms in (272) have the same symmetry properties in
the index pairs $jj^{\prime}$ and $ii^{\prime}$, which reflects that there are
no transitions between states corresponding to different representations of
the rotation group.

##### Quark polarisation.

We now turn our attention to quark distributions and their renormalisation.
The graphs to be computed for the matching coefficients are as in figure 11,
with appropriate replacements of gluons by quarks and with eikonal lines in
the fundamental representation. The spinor indices at the top of the graphs
are to be contracted with the Dirac matrix of the relevant operator $O_{q}$,
$O_{\Delta q}$ or $O_{\delta q}$, and corresponding contractions are done for
the Dirac indices of the quark states in the matrix element.

As is well known, the treatment of quark polarisation in this context requires
a consistent definition of $\gamma_{5}$ and of the $\epsilon$ tensor in
dimensional regularisation. We use the scheme of ’t Hooft, Veltman and
Breitenlohner, Maison (HVBM) tHooft:1972tcz ; Breitenlohner:1977hr , where
$\gamma_{5}=\frac{i}{4!}\mskip 1.5mu\epsilon_{\lambda\mu\nu\rho}\mskip
1.5mu\gamma^{\lambda}\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}$. Here
$\epsilon_{\lambda\mu\nu\rho}$ is the usual antisymmetric tensor if all
indices are in the four physical dimensions (with $\epsilon_{0123}=+1$) and
zero otherwise. We also introduce its counterpart
$\epsilon^{ij}=\epsilon^{+-ij}$, where $i,j$ take values in the $D-2$
transverse dimensions, as well as the tensor $\bar{\delta}^{ij}$, which equals
$\delta^{ij}$ if $i,j=1,2$ and is zero for one or both of $i,j$ in the
unphysical dimensions.

For transverse quark polarisation, we avoid the use of $\gamma_{5}$ by taking
the Dirac matrix $\Gamma_{\delta q}^{j}$ in the operator
$\mathcal{O}_{\delta}$ as

$\displaystyle\Gamma_{\delta q}^{j}$ $\displaystyle=\frac{1}{2}\mskip
1.5mu\epsilon^{jj^{\prime}}\sigma^{j^{\prime}+}=\frac{i}{2}\mskip
1.5mu\epsilon^{jj^{\prime}}\gamma^{j^{\prime}}\gamma^{+}\,,$ (273)

which is equal to the form $\tfrac{i}{2}\mskip 1.5mu\sigma^{j+}\gamma_{5}$ in
(15) in the physical dimensions, $j=1,2$, whereas it differs for values of $j$
in the unphysical sector. The fermion trace in the matching coefficient
$C_{\delta q\mskip 1.5mu\delta q}$ then contains two matrices $\sigma^{j+}$
and $\sigma^{i-}$. This gives exactly the same result as one would obtain for
the trace with $\sigma^{j+}\gamma_{5}$ and $\sigma^{i-}\gamma_{5}$ using the
“naive dimensional regularisation” prescription, where one assumes that
$\gamma_{5}$ has square $1$ and anticommutes with all $\gamma^{\mu}$. That
prescription is often used in the context of transverse quark polarisation,
but we prefer the formulation with (273). Note that with “naive dimensional
regularisation” one cannot treat all polarised matching coefficients on an
equal footing, since it can only be used to compute traces with an even number
of $\gamma_{5}$ matrices. The matrix $\sigma^{i+}$ has also been used to
define quark transversity in the recent work Gutierrez-Reyes:2017glx , where
the matching coefficients for polarised TMDs were revisited in a systematic
way.

For longitudinal quark polarisation, the matrix $\gamma_{5}$ is unavoidable.
In the HVBM scheme, one can rewrite the Dirac matrix in the operator
$\Gamma_{\Delta q}$ as

$\displaystyle\Gamma_{\Delta q}$ $\displaystyle=\frac{1}{2}\mskip
1.5mu\gamma^{+}\gamma_{5}=-\frac{i}{4}\mskip 1.5mu\epsilon^{jj^{\prime}}\mskip
1.5mu\gamma^{[j^{\prime}}\gamma^{j]}\mskip 1.5mu\gamma^{+}\,,$ (274)

where we define $\gamma^{[i}\gamma^{j]}=\tfrac{1}{2}\mskip
1.5mu\bigl{(}\gamma^{i}\gamma^{j}-\gamma^{j}\gamma^{i}\bigr{)}$. A
corresponding replacement can be made for the matrix $\gamma^{-}\gamma_{5}$
that appears in the fermion trace for matrix elements with incoming
longitudinally polarised quarks. Postponing the contraction with the
$\epsilon$ tensor to the end of the calculation, one can hence represent
longitudinally polarised quarks by an operator $\mathcal{O}_{\Delta
q}^{jj^{\prime}}=-\frac{1}{4}\mskip 1.5mu\bar{q}\mskip
1.5mu\gamma^{[j^{\prime}}\gamma^{j]}\gamma^{+}\mskip 1.5muq$ with two
transverse indices. This makes it easy to see that one has mixing under
renormalisation with longitudinally polarised gluons, represented by the
operator $\mathcal{O}_{\Delta g}^{jj^{\prime}}=G^{+[j^{\prime}}\mskip
1.5muG^{+j]}$ that has the same transformation behaviour under rotations. For
all matching coefficients with longitudinal polarisation, ${\Delta q\Delta
q}$, ${\Delta g\Delta q}$, ${\Delta q\Delta g}$ and ${\Delta g\Delta g}$, one
then has a tensor matching coefficient $C_{ab}^{jj^{\prime},ii^{\prime}}$ that
is antisymmetric in $jj^{\prime}$. As the only terms with ultraviolet
divergences in a matching coefficient arise from the matrix elements at
$\boldsymbol{z}=\boldsymbol{0}$, they must be independent of $\boldsymbol{z}$.
This leaves
$\delta^{ij}\delta^{i^{\prime}j^{\prime}}-\delta^{ij^{\prime}}\delta^{i^{\prime}j}$
as the only possible tensor structure of the UV divergent part of the matching
coefficient. This is also the tensor structure of the tree-level matrix
element, which multiplies the one-loop counterterm according to (270). For the
renormalised matching coefficient one thus has

$\displaystyle C^{(1)\mskip 1.5mujj^{\prime},ii^{\prime}}$
$\displaystyle=\frac{1}{2}\mskip
1.5mu\bigl{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}-\delta^{ij^{\prime}}\delta^{i^{\prime}j}\bigr{)}\,B^{(1)}(x,\epsilon)-\frac{1}{2}\mskip
1.5mu\bigl{(}\delta^{ij}\delta^{i^{\prime}j^{\prime}}-\delta^{ij^{\prime}}\delta^{i^{\prime}j}\bigr{)}\,Z^{(1)}(x,\epsilon)$
$\displaystyle\quad+\\{\text{terms finite for $\epsilon\to 0$}\\}\,,$ (275)

where $B^{(1)}$ contains all ultraviolet divergent parts of the bare matrix
element and $Z^{(1)}$ is the relevant one-loop renormalisation factor. The
$\overline{\text{MS}}$ prescription fixes the latter to be $S_{\epsilon}$
times the pole part of $B^{(1)}$. The outcome of this discussion is that one
obtains the same counterterm (and hence the same renormalised result) when
using the tensor form (7.3.3) or when contracting it with

$\displaystyle\frac{1}{2}\,\epsilon^{jj^{\prime}}\epsilon^{ii^{\prime}}$
$\displaystyle=\frac{1}{2}\mskip
1.5mu\bigl{(}\bar{\delta}^{ij}\bar{\delta}^{i^{\prime}j^{\prime}}-\bar{\delta}^{ij^{\prime}}\bar{\delta}^{i^{\prime}j}\bigr{)}$
(276)

and working with scalar matching coefficients, which corresponds to using the
conventional operators ${\textstyle\frac{1}{2}}\mskip 1.5mu\bar{q}\mskip
1.5mu\gamma^{+}\gamma_{5}\mskip 1.5muq$ and
$i\epsilon^{jj^{\prime}}G^{+j^{\prime}}\mskip 1.5muG^{+j}$ in the HVBM scheme.
Working with (7.3.3) has the advantage that one does not need to distinguish
between physical and unphysical dimensions during the computation.

In the computations of Echevarria:2015uaa  and Gutierrez-Reyes:2017glx , the
tensor $\bar{\delta}$ on the r.h.s. of (276) was replaced by the full
transverse metric $\delta$ in $D-2$ dimensions, as a modification of the
proposal by Larin Larin:1993tq , where
$\epsilon^{\lambda\mu\nu\rho}\epsilon^{\lambda^{\prime}\mu^{\prime}\nu^{\prime}\rho^{\prime}}$
was replaced by products of $D$ dimensional instead of $4$ dimensional metric
tensors (cf. Moch:2015usa  for a recent discussion). We see that one obtains
the same result when one contracts (7.3.3) with
${\textstyle\frac{1}{2}}\bigl{(}{\delta}^{ij}{\delta}^{i^{\prime}j^{\prime}}-{\delta}^{ij^{\prime}}{\delta}^{i^{\prime}j}\bigr{)}$
instead of (276) and imposes $\overline{\text{MS}}$ subtraction. The essential
point for this is that $B^{(1)}$ and $Z^{(1)}$ have the _same_ $\epsilon$
dependent prefactor in this case. Since poles in $1/\epsilon$ cancel in
$B^{(1)}-Z^{(1)}$ by construction, only the $\epsilon\to 0$ limit of their
prefactor enters in the renormalised matching coefficient. The preceding
argument readily generalises to order $\alpha_{s}^{n}$, where one should
replace $B^{(1)}-Z^{(1)}$ with $\sum_{m}B^{(n-m)}\otimes Z^{(m)}$, where
$B^{(0)}=Z^{(0)}=\delta(1-x)$. If, however, the divergent part of a quantity
involves different tensor structures, then the procedure just described is no
longer guaranteed to work correctly.

We note that (274) is a special case of the identity

$\displaystyle\frac{1}{2}\mskip
1.5mu\bigl{(}\gamma_{\lambda}\gamma_{5}-\gamma_{5}\gamma_{\lambda}\bigr{)}$
$\displaystyle=\frac{i}{3!}\,\epsilon_{\lambda\mu\nu\rho}\,\gamma^{[\mu}\gamma^{\nu}\gamma^{\rho]}$
(277)

in $D$ dimensions, which has been used for a long time when discussing the
axial current in the HVBM scheme Collins:1984xc ; Larin:1993tq . The explicit
antisymmetrisation on the l.h.s. of (277), necessary to make the current
Hermitian, can be omitted in (274) because $\\{\gamma^{+},\gamma_{5}\\}=0$. As
shown in Collins:1984xc ; Larin:1993tq , the divergence of the flavour
nonsinglet axial current is nonzero if one uses $\overline{\text{MS}}$
renormalisation and the HVBM scheme, but a zero divergence can be achieved by
an additional finite renormalisation. The resulting nonsinglet current has
zero anomalous dimension. At order $\alpha_{s}$, the same finite
renormalisation achieves that the divergence of the flavour singlet axial
current is given by the Adler-Bell-Jackiw anomaly, as shown in Larin:1993tq
(at higher orders, the required renormalisation factors differ for the flavour
singlet and nonsinglet currents). In Vogelsang:1996im ; Matiounine:1998re
this finite renormalisation was extended to the nonlocal axial current
$\mathcal{O}_{\Delta q}$ that defines longitudinally polarised quark
distributions, see also Weber:1991wd ; Stratmann:1995fn . The finite
renormalisation can be written as

$\displaystyle O^{\text{NS}}_{\Delta q}$
$\displaystyle=Z_{5}^{\text{NS}}\otimes O^{\text{NS}}_{\Delta
q,\overline{\text{MS}}}\,,$ $\displaystyle O^{\text{S}}_{\Delta q}$
$\displaystyle=Z_{5}^{\text{S}}\otimes O^{\text{S}}_{\Delta
q,\overline{\text{MS}}}\,,$ $\displaystyle O_{\Delta g}$
$\displaystyle=O_{\Delta g,\overline{\text{MS}}}$ (278)

for the flavour singlet and nonsinglet operators in $x$ space, where the
currents on the r.h.s. are renormalised by the standard $\overline{\text{MS}}$
prescription. At order $\alpha_{s}$ one has

$\displaystyle Z_{5}^{\text{NS}}(x)$
$\displaystyle=Z_{5}^{\text{S}}(x)=1-\frac{\alpha_{s}}{2\pi}\,4C_{F}\mskip
1.5mu(1-x)\,;$ (279)

higher orders can be found in Moch:2014sna .151515Note that $Z_{5}$, called
$Z$ in Moch:2014sna , is different from $Z^{5}$ given in Matiounine:1998re .
In that work, the renormalisation of the bare current is written in terms of
the product $Z^{-1}Z^{5}$ of two renormalisation factors (or matrices in the
singlet sector), where $Z^{-1}$ does _not_ correspond to the
$\overline{\text{MS}}$ prescription. This finite renormalisation achieves in
particular that the DGLAP kernels $P_{qq}$ and $P_{\Delta q\Delta q}$ agree up
to order $\alpha_{s}^{2}$ Vogelsang:1996im ; Matiounine:1998re . It also
achieves that the hard-scattering coefficients of the Drell-Yan process in
collinear factorisation satisfy $W_{\Delta q\mskip
1.5mu\Delta\bar{q}}=-W_{q\bar{q}}$ at order $\alpha_{s}$ Weber:1991wd ;
Ravindran:2002na , where $W_{\Delta q\mskip 1.5mu\Delta\bar{q}}$ is for the
incoming quark and antiquark both being longitudinally polarised and
$W_{q\bar{q}}$ for both being unpolarised. In the results for the matching
kernels given in the next subsection, we have included the finite
renormalisation (278). At the order we are working at, this only affects the
coefficient $C_{\Delta q\Delta q}$.

#####  $\overline{\text{MS}}$ variants.

All results given in this paper so far refer to the prescription that a
renormalisation factor $Z$ has the form

$\displaystyle Z$ $\displaystyle=\sum_{n=0}^{\infty}(S_{\epsilon}\mskip
1.5mu\alpha_{s})^{n}\,Z^{(n)}\,,$ (280)

where $\alpha_{s}$ is the renormalised coupling and $Z^{(n)}$ for $n\geq 1$ is
a finite sum of poles in $\epsilon$, fixed uniquely by the requirement that
the renormalised quantity be finite at $\epsilon=0$. Equivalently, one may
express all quantities in terms of the rescaled coupling
$\bar{\alpha}_{s}=S_{\epsilon}\mskip 1.5mu\alpha_{s}$ and use minimal
subtraction, where all counterterms are sums of poles. The factor
$S_{\epsilon}$ is given in (259) and was introduced in Bardeen:1978yd  to
simplify the finite terms left after renormalisation.

A variant of the scheme has been proposed in chapter 3.2.6 of Collins:2011zzd
, where instead of $S_{\epsilon}$ in (259) one uses

$\displaystyle S_{\epsilon}^{\text{JCC}}$
$\displaystyle=\frac{(4\pi)^{\epsilon}}{\Gamma(1-\epsilon)}\,.$ (281)

The difference between the two versions is

$\displaystyle S_{\epsilon}-S_{\epsilon}^{\text{JCC}}$
$\displaystyle=\frac{\pi^{2}}{12}\,\epsilon^{2}+\mathcal{O}(\epsilon^{3})\,.$
(282)

For one-loop quantities with only a single pole in $\epsilon$, this is of no
consequence, but it does matter for quantities with a double pole. From (263)
we see that the counterterm is shifted by $\pi^{2}/12$ times the coefficient
of the double pole if one uses $S_{\epsilon}^{\text{JCC}}$ instead of
$S_{\epsilon}$. For the one-loop matching coefficients, this scheme change
affects the soft factor and simply removes the term $-\pi^{2}/12$ in the
quantity $S_{L}$ defined in (266).

Like any choice of renormalisation prescription, the choice of $S_{\epsilon}$
must drop out in physical quantities such as the overall cross section. At
$\mathcal{O}(\alpha_{s})$ this choice affects the matching coefficients
${}^{R}C_{ab}$, the hard-scattering coefficients $C_{H}^{2}$ (with squared
logarithms in (7.1) indicating the presence of double poles $1/\epsilon^{2}$
before renormalisation), as well as DPDFs in colour non-singlet channels and
DTMDs in any colour channel (via one-loop renormalisation of the soft factors
in their definitions). Collins-Soper and DGLAP kernels, as well as anomalous
dimensions, are not affected by the choice of $S_{\epsilon}$ at
$\mathcal{O}(\alpha_{s})$.

####  7.3.4 Matching coefficients for DTMDs on DPDFs

We now give our final results for the matching coefficients for DTMDs at large
$\boldsymbol{y}$ in the colour singlet sector, which are equal to the matching
coefficients for single parton TMDs. We have computed them not only with the
rapidity regulator used in this paper, but also with the so-called $\delta$
regulator specified in GarciaEchevarria:2011rb , and we find complete
agreement between the two versions.

In the gluon sector, we start from the tensor form (267) and project onto
unpolarised or longitudinally polarised gluons, which gives

$\displaystyle C_{gg}(x,\boldsymbol{z})=\delta(1-x)$
$\displaystyle-L\,\biggl{[}\,P_{gg}(x)-\frac{1}{2}\,\delta(1-x)\,\gamma_{g}\,\biggr{]}+\frac{\alpha_{s}}{2\pi}\,C_{A}\mskip
1.5muS_{L}\,\delta(1-x)\,,$ (283) $\displaystyle C_{\Delta g\Delta
g}(x,\boldsymbol{z})=\delta(1-x)$ $\displaystyle-L\,\biggl{[}\,P_{\Delta
g\Delta
g}(x)-\frac{1}{2}\,\delta(1-x)\,\gamma_{g}\,\biggr{]}+\frac{\alpha_{s}}{2\pi}\,C_{A}\mskip
1.5muS_{L}\,\delta(1-x)$
$\displaystyle-\frac{\alpha_{s}}{2\pi}\,4C_{A}\,(1-x)\,,$ (284)

where we recall that

$\displaystyle L$
$\displaystyle=\log\frac{\mu^{2}\boldsymbol{z}^{2}}{b_{0}^{2}}\,,$
$\displaystyle S_{L}$
$\displaystyle=-\frac{1}{2}L^{2}+L\,\log\frac{\mu^{2}}{\zeta}-\sigma\mskip
1.5mu\frac{\pi^{2}}{12}$

with $\sigma=1$ if one takes $S_{\epsilon}$ from (259) for the
$\overline{\text{MS}}$ counterterms, whereas $\sigma=0$ if one takes
$S_{\epsilon}^{\text{JCC}}$ from (281). The DGLAP splitting kernels
$P_{gg}(x)$ and $P_{\Delta g\Delta g}(x)$ include a term proportional to the
anomalous dimension $\gamma_{g}=\beta_{0}\mskip 1.5mu\alpha_{s}/(2\pi)$, which
comes from virtual graphs (see figure 11). This term is subtracted again in
the above matching coefficients, which only receive contributions from real
graphs as explained earlier. We list all kernels $P_{ab}$ at the end of this
subsection.

Matching coefficients involving linear gluon polarisation carry tensor
indices, and we extract from (267) the forms

$\displaystyle C_{\delta g\mskip 1.5mu\delta g}^{jj^{\prime}\\!,ii^{\prime}}$
$\displaystyle=\tau^{jj^{\prime}\\!,ii^{\prime}}\mskip 1.5muC_{\delta g\mskip
1.5mu\delta g}\,,\phantom{\biggl{(}\frac{1}{2}\biggr{)}}$ (285) $\displaystyle
C_{\delta g\mskip 1.5mug}^{jj^{\prime}}$
$\displaystyle=\biggl{(}\frac{\boldsymbol{z}^{j}\boldsymbol{z}^{j^{\prime}}}{\boldsymbol{z}^{2}}-\frac{1}{2}\,\delta^{jj^{\prime}}\biggr{)}\,C_{\delta
g\mskip 1.5mug}\,,$ (286) $\displaystyle C_{g\mskip 1.5mu\delta
g}^{ii^{\prime}}$
$\displaystyle=\biggl{(}\frac{2\boldsymbol{z}^{i}\boldsymbol{z}^{i^{\prime}}}{\boldsymbol{z}^{2}}-\delta^{ii^{\prime}}\biggr{)}\,C_{g\mskip
1.5mu\delta g}\,,$ (287)

where $\tau^{jj^{\prime}\\!,ii^{\prime}}$ is defined in (19) and the scalar
coefficient functions read

$\displaystyle C_{\delta g\mskip 1.5mu\delta g}(x,\boldsymbol{z})$
$\displaystyle=\delta(1-x)-L\,\biggl{[}\,P_{\delta g\mskip 1.5mu\delta
g}(x)-\frac{1}{2}\,\delta(1-x)\,\gamma_{g}\,\biggr{]}+\frac{\alpha_{s}}{2\pi}\,C_{A}\mskip
1.5muS_{L}\,\delta(1-x)\,,$ (288) $\displaystyle C_{\delta g\mskip
1.5mug}(x,\boldsymbol{z})$
$\displaystyle=-\frac{\alpha_{s}}{2\pi}\,2C_{A}\,\frac{1-x}{x}\,,$ (289)
$\displaystyle C_{g\mskip 1.5mu\delta g}(x,\boldsymbol{z})$
$\displaystyle=-\frac{\alpha_{s}}{2\pi}\,2C_{A}\,x(1-x)\,.$ (290)

The prefactor in (286) has been chosen such that $C_{\delta g\mskip 1.5mug}$
coincides with the corresponding matching coefficient defined for single gluon
TMDs in Echevarria:2015uaa ; Gutierrez-Reyes:2017glx .

In the quark sector, we obtain

$\displaystyle C_{qq}(x,\boldsymbol{z})$
$\displaystyle=\delta(1-x)-L\,\biggl{[}\,P_{qq}(x)-\frac{1}{2}\,\delta(1-x)\,\gamma_{q}\,\biggr{]}+\frac{\alpha_{s}}{2\pi}\,C_{F}\mskip
1.5muS_{L}\,\delta(1-x)$ $\displaystyle\qquad\qquad\hskip
8.50006pt+\frac{\alpha_{s}}{2\pi}\,C_{F}\,(1-x)\,,$ (291) $\displaystyle
C_{\Delta q\Delta q}(x,\boldsymbol{z})$
$\displaystyle=C_{qq}(x,\boldsymbol{z})\,,$ (292)

with $\gamma_{q}=3C_{F}\mskip 1.5mu\alpha_{s}/(2\pi)$. For transverse quark
polarisation, one has again a tensor valued coefficient

$\displaystyle C_{\delta q\mskip 1.5mu\delta q}^{j,i}$
$\displaystyle=\delta^{ji}_{\phantom{\delta q}}C_{\delta q\mskip 1.5mu\delta
q}^{\phantom{j}}$ (293)

with

$\displaystyle C_{\delta q\mskip 1.5mu\delta q}(x,\boldsymbol{z})=\delta(1-x)$
$\displaystyle-L\,\biggl{[}\,P_{\delta q\mskip 1.5mu\delta
q}(x)-\frac{1}{2}\,\delta(1-x)\,\gamma_{q}\,\biggr{]}+\frac{\alpha_{s}}{2\pi}\,C_{F}\mskip
1.5muS_{L}\,\delta(1-x)\,.$ (294)

We note that the tensor structure $\boldsymbol{z}^{j}\boldsymbol{z}^{i}$ would
be allowed in (293) but does not appear at $\mathcal{O}(\alpha_{s})$. For
transitions between quarks and gluons, we obtain

$\displaystyle C_{qg}(x,\boldsymbol{z})$
$\displaystyle=-L\,P_{qg}(x)+\frac{\alpha_{s}}{2\pi}\,2T_{F}\,x(1-x)\,,$ (295)
$\displaystyle C_{\Delta q\Delta g}(x,\boldsymbol{z})$
$\displaystyle=-L\,P_{\Delta q\Delta
g}(x)+\frac{\alpha_{s}}{2\pi}\,2T_{F}(1-x)$ (296)

and

$\displaystyle C_{gq}(x,\boldsymbol{z})$
$\displaystyle=-L\,P_{gq}(x)+\frac{\alpha_{s}}{2\pi}\,C_{F}\,x\,,$ (297)
$\displaystyle C_{\Delta g\Delta q}(x,\boldsymbol{z})$
$\displaystyle=-L\,P_{\Delta g\Delta
q}(x)-\frac{\alpha_{s}}{2\pi}\,2C_{F}(1-x)\,,$ (298)

where we have re-introduced the normalisation factor $T_{F}=1/2$. The
analogues of (289) and (290) for quark-gluon transitions are

$\displaystyle C_{\delta g\mskip 1.5muq}(x,\boldsymbol{z})$
$\displaystyle=-\frac{\alpha_{s}}{2\pi}\,2C_{F}\,\frac{1-x}{x}\,,$ (299)
$\displaystyle C_{q\mskip 1.5mu\delta g}(x,\boldsymbol{z})$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\;2T_{F}\,x(1-x)\,,$ (300)

where the scalar coefficients are defined as in (286) and (287) with the
subscript $g$ replaced by $q$ for the unpolarised parton.

DPDs associated with linear gluon or transverse quark polarisation are tensor
valued and can be decomposed into basis tensors that multiply scalar (or
pseudoscalar) distributions. We give these decompositions in appendix F and
list the resulting matching equations between DTMDs and DPDFs in appendix G.

Matching coefficients $C_{ab}$ with exactly one longitudinal polarisation
($\Delta g$ or $\Delta q$) are zero due to parity invariance. Coefficients
with exactly one transverse quark polarisation $\delta q$ vanish because the
relevant graphs involve the trace of an odd number of $\gamma$ matrices. Due
to charge conjugation invariance, the coefficients $C_{ab}$ remain the same if
one replaces all quark indices by antiquark ones. These symmetry relations
hold at all orders in $\alpha_{s}$. At one-loop level, there are no
transitions between quarks and antiquarks or between quarks of different
flavours. These appear starting at $\mathcal{O}(\alpha_{s}^{2})$.

Most of the above coefficients have been calculated in the literature before.
The coefficients in (283), (284), (289) and (297) to (299) were calculated in
Echevarria:2015uaa , using SCET and the $\delta$ regulator in the version of
GarciaEchevarria:2011rb . We find agreement between our results and those in
the arXiv version 5 of Echevarria:2015uaa . Our results in (291) to (296)
agree with Aybat:2011zv  (which contains unpolarised results only) and with
Bacchetta:2013pqa . The results in (283), (286), (297) and (299) agree with
Catani:2013tia  and Becher:2012yn . Apart from $C_{\delta g\mskip 1.5mu\delta
g}$, $C_{g\mskip 1.5mu\delta g}$ and $C_{q\mskip 1.5mu\delta g}$, all matching
coefficients were recently calculated in Gutierrez-Reyes:2017glx , using the
$\delta$ regulator briefly described in our appendix B. We agree with the
results in the arXiv version 2 of Gutierrez-Reyes:2017glx .

To the best of our knowledge, $C_{\delta g\mskip 1.5mu\delta g}$, $C_{g\mskip
1.5mu\delta g}$ and $C_{q\mskip 1.5mu\delta g}$ in (288), (290) and (300) have
not been given in the literature before. They are not relevant for single
parton TMDs in a nucleon, since the PDF for a linearly polarised gluon
vanishes in that case. This is because linear gluon polarisation corresponds
to the interference between gluons with helicity $+1$ and $-1$ in the
amplitude and its conjugate.

Let us finally list the leading-order DGLAP splitting functions, which were
first derived in Altarelli:1977zs ; Artru:1989zv . They are given by

$\displaystyle P_{qq}(x)$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\,C_{F}\,\frac{1+x^{2}}{(1-x)_{+}}+\frac{1}{2}\,\delta(1-x)\,\gamma_{q}\,,$
(301) $\displaystyle P_{\Delta q\Delta q}(x)$ $\displaystyle=P_{qq}(x),$ (302)
$\displaystyle P_{\delta q\mskip 1.5mu\delta q}(x)$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\,C_{F}\,\frac{2x}{(1-x)_{+}}+\frac{1}{2}\,\delta(1-x)\,\gamma_{q}$
(303)

for quark-quark transitions and by

$\displaystyle P_{gg}(x)$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\,2C_{A}\,\biggl{[}\,\frac{x}{(1-x)_{+}}+\frac{(1-x)(1+x^{2})}{x}\,\biggr{]}+\frac{1}{2}\,\delta(1-x)\,\gamma_{g}\,,$
(304) $\displaystyle P_{\Delta g\Delta g}(x)$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\,2C_{A}\,\biggl{[}\,\frac{x}{(1-x)_{+}}+2(1-x)\,\biggr{]}+\frac{1}{2}\,\delta(1-x)\,\gamma_{g}\,,$
(305) $\displaystyle P_{\delta g\mskip 1.5mu\delta g}(x)$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\,2C_{A}\,\frac{x}{(1-x)_{+}}+\frac{1}{2}\,\delta(1-x)\,\gamma_{g}$
(306)

for gluons. The splitting functions mixing quarks and gluons are

$\displaystyle P_{qg}(x)$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\,T_{F}\Big{[}x^{2}+(1-x)^{2}\Big{]}\,,$
$\displaystyle P_{gq}(x)$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\,C_{F}\,\frac{1+(1-x)^{2}}{x}\,,$ (307)
$\displaystyle P_{\Delta q\Delta g}(x)$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\,T_{F}\,\Big{[}x^{2}-(1-x)^{2}\Big{]}\,,$
$\displaystyle P_{\Delta g\Delta q}(x)$
$\displaystyle=\frac{\alpha_{s}}{2\pi}\,C_{F}\,\frac{1-(1-x)^{2}}{x}\,.$ (308)

The symmetry properties of the matching coefficients $C_{ab}$ detailed below
(300) also hold for the splitting functions $P_{ab}$. In addition, $P_{g\mskip
1.5mu\delta g}=P_{\delta g\mskip 1.5mug}=0$ because the collinear operators
$O_{g}(x)$ and $O_{\delta g}^{jj^{\prime}}(x)$ transform differently under
rotations and hence cannot mix under renormalisation.

####  7.3.5 Colour non-singlet channels

It is easy to derive the matching coefficients in colour non-singlet channels
from the results we have just presented. All real graphs in figure 11 have the
same colour factor, so that the colour projections in (7.3.1) affect all
graphs in the same way. The same holds for the graphs for quark-quark
transitions and for transitions between quarks and gluons (in the latter case,
the soft factor does not contribute to the one-loop matching coefficients and
one only needs the analogues of graphs 11a and b). We can therefore write

$\displaystyle{}^{R}C^{(1)}_{ab}(x,\boldsymbol{z};\mu,\zeta)$
$\displaystyle=c_{ab}(R)\;{}^{1}\mskip-1.5muC^{(1)}_{ab}(x,\boldsymbol{z};\mu,\zeta)\,.$
(309)

The ratios $c_{ab}(R)$ of colour factors between the representation $R$ and
the singlet channel can be determined from the basic ladder graphs, which have
the topology of figure 11a. They are readily obtained from the evolution
kernels given in section 5.1.3 in Diehl:2011yj  and read

$\displaystyle c_{qq}(8)$ $\displaystyle=-\frac{1}{N^{2}-1}\,,$ $\displaystyle
c_{gq}(A)$ $\displaystyle=c_{qg}(A)=\sqrt{\frac{N^{2}}{2(N^{2}-1)}}\,,$
$\displaystyle c_{gg}(A)$ $\displaystyle=c_{gg}(S)=\frac{1}{2}\,,$
$\displaystyle c_{gq}(S)$
$\displaystyle=c_{qg}(S)=\sqrt{\frac{N^{2}-4}{2(N^{2}-1)}}$ (310)

and for $N=3$

$\displaystyle c_{gg}(D)$ $\displaystyle=0\,,$ $\displaystyle c_{gg}(27)$
$\displaystyle=-1/3\,.$ (311)

For transitions involving antiquarks one has
$c_{\bar{q}\bar{q}}(8)=c_{qq}(8)$, $c_{g\bar{q}}(S)=c_{\bar{q}g}(S)=c_{qg}(S)$
and $c_{g\bar{q}}(A)=c_{\bar{q}g}(A)=-c_{qg}(A)$. The minus sign in the last
case reflects the fact that two gluons coupled to an antisymmetric octet have
negative charge parity.

We can now derive the explicit form of the LO DGLAP kernels in a general
colour representation $R$. For this, we use the evolution equation (6.2.2) for
the matching coefficients ${}^{R}C_{ab}$. At order $\alpha_{s}$, we can
replace ${}^{R}C_{ab}$ on the r.h.s. by its lowest-order value
$\delta_{ab}\,\delta(1-x^{\prime})$ and thus obtain

$\displaystyle{}^{R}\mskip-1.5muP_{ac}(x;\mu,\zeta)$
$\displaystyle=\frac{1}{2}\,\delta_{ac}\,\delta(1-x)\,\gamma_{F,c}(\mu,\zeta)-\frac{1}{2}\,\frac{\partial}{\partial\log\mu}{}^{R}C_{ac}(x,\boldsymbol{z};\mu,\zeta)\,.$
(312)

Combining this with (309), we readily find that the leading-order DGLAP
evolution kernels of collinear DPDs in colour non-singlet channels are given
by

$\displaystyle{}^{R}\mskip-1.5muP_{ab}(x;\mu,\zeta)$
$\displaystyle=\frac{1}{2}\,\delta_{ab}\,\delta(1-x)\,\gamma_{F,a}(\mu,\zeta)-\frac{c_{ab}(R)}{2}\,\frac{\partial}{\partial\log\mu}{}^{1}C_{ab}(x,\boldsymbol{z};\mu,\zeta)$
$\displaystyle=\frac{1-c_{aa}(R)}{2}\,\delta_{ab}\,\delta(1-x)\,\gamma_{F,a}(\mu,\zeta)+c_{ab}(R)\,{}^{1}\mskip-1.5muP_{ab}(x;\mu)$
$\displaystyle=-\frac{1}{2}\,\delta_{ab}\,\delta(1-x)\,{}^{R}\gamma_{J}(\mu)\,\log\frac{\sqrt{\zeta}}{\mu}$
$\displaystyle\quad+\frac{1-c_{aa}(R)}{2}\,\delta_{ab}\,\delta(1-x)\,\gamma_{a}(\mu)+c_{ab}(R)\,{}^{1}\mskip-1.5muP_{ab}(x;\mu)\,.$
(313)

They are hence easily recovered from their counterparts in the colour singlet
channel. In the last step we have brought the kernel into the form (143),
using the explicit $\zeta$ dependence (5) of $\gamma_{F,a}$, the general
relation (172) between anomalous dimensions, and the one-loop relation
${}^{R}\gamma_{K,a}=c_{aa}(R)\,{}^{1}\mskip-1.5mu\gamma_{K,a}$, which follows
from (309). We recall that the colour singlet kernels
${}^{1}\mskip-1.5muP_{aa}(x;\mu)$ contain a contribution
$\delta(1-x)\,\gamma_{a}/2$ from virtual graphs, see (301) to (306). The last
line of (7.3.5) simply reflects that these graphs have the same colour factor
for all channels $R$.

Using (7.3.5), one can rewrite the $\zeta$ independent distributions
$\widehat{F}$ introduced in (5.2) as

$\displaystyle{}^{R}\mskip-1.5mu\widehat{F}_{a_{1}a_{2},\,\mu_{0},\zeta_{0}}(x_{i},\boldsymbol{y};\mu_{i})$
$\displaystyle=\exp\,\biggl{[}\bigl{(}1-c_{a_{1}a_{1}}(R)\bigr{)}\int_{\mu_{0}}^{\mu_{1}}\frac{d\mu}{\mu}\,\gamma_{a_{1}}(\mu)$
$\displaystyle\qquad\,+\bigl{(}1-c_{a_{2}a_{2}}(R)\bigr{)}\int_{\mu_{0}}^{\mu_{2}}\frac{d\mu}{\mu}\,\gamma_{a_{2}}(\mu)\mskip
1.5mu\biggr{]}\;{}^{R}\mskip-1.5mu\widetilde{F}_{a_{1}a_{2},\,\mu_{0},\zeta_{0}}(x_{i},\boldsymbol{y};\mu_{i})\,,$
(314)

where $\widetilde{F}$ satisfies the evolution equation (5.2) with kernels
$c_{a_{1}b_{1}}(R)\,{}^{1}\mskip-1.5muP_{a_{1}b_{1}}(x_{1}^{\prime};\mu_{1},\mu_{1}^{2})$
instead of
${}^{R}\mskip-1.5muP_{a_{1}b_{1}}(x_{1}^{\prime};\mu_{1},\mu_{1}^{2})$. One
can thus use numerical code for the one-loop evolution of colour singlet DPDs
by rescaling the evolution kernels.

###  7.4 Splitting kernels for DPDs

Both DTMDs and DPDFs can be matched on single-parton distributions at small
$\boldsymbol{y}$ as specified in (154), (193) and (6.3.1). At
$\mathcal{O}(\alpha_{s})$, all three matching equations involve the same
kernels ${}^{R}T_{a_{0}\to a_{1}a_{2}}^{ll^{\prime}}$, which were computed in
section 5.2.2 of Diehl:2011yj . We list them here in the notation of
Diehl:2017kgu  (arXiv version 2), where compared to Diehl:2011yj  the kernels
$T$ have the opposite order of indices and include a colour factor. Thus one
has for instance

$\displaystyle{}^{1}T^{ll^{\prime}}_{g\to q\bar{q}}\,\big{|}_{\text{here}}$
$\displaystyle=T_{F}\mskip 1.5muT^{l^{\prime}\mskip-1.5mul}_{g\to
q\bar{q}}\,\big{|}_{\text{Ref.\leavevmode\nobreak\
\smash{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Diehl:2011yj}{\@@citephrase{(}}{\@@citephrase{)}}}}}}\,,$
$\displaystyle{}^{1}T^{ll^{\prime}}_{q\to g\mskip
1.5mu\Delta{q}}\,\big{|}_{\text{here}}$ $\displaystyle=C_{F}\mskip
1.5muT^{l^{\prime}\mskip-1.5mul}_{q\to g\mskip
1.5mu\Delta{q}}\,\big{|}_{\text{Ref.\leavevmode\nobreak\
\smash{\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Diehl:2011yj}{\@@citephrase{(}}{\@@citephrase{)}}}}}}\,.$
(315)

The new assignment of indices is such that $l$ refers to the amplitude and
$l^{\prime}$ to its complex conjugate. We must note a mistake in equation
(5.62) of Diehl:2011yj , where the correct order of indices in the kernels $T$
and $U$ is $l^{\prime}\mskip-1.5mul$ and not $ll^{\prime}$.

We only give the kernels for the splitting of an unpolarised parton $a_{0}$
here. In the matching equation (154) for DTMDs on TMDs, one can also have
transverse quark or linear gluon polarisation for $a_{0}$, but this
possibility is absent if one matches on collinear PDFs in an unpolarised
proton.

We start with the kernels for DPDs in the colour singlet channel. For a gluon
splitting into a quark and an antiquark, one has

$\displaystyle{}^{1}T_{g\to q\bar{q}}^{ll^{\prime}}(u)$
$\displaystyle={}-{}^{1}T_{g\to\Delta q\Delta\bar{q}}^{ll^{\prime}}(u)$
$\displaystyle=T_{F}\mskip 1.5mu(u^{2}+\bar{u}^{2})\,\delta^{ll^{\prime}}\,,$
$\displaystyle{}^{1}T_{g\to\Delta q\mskip 1.5mu\bar{q}}^{ll^{\prime}}(u)$
$\displaystyle={}-{}^{1}T_{g\to q\Delta\bar{q}}^{ll^{\prime}}(u)$
$\displaystyle=-iT_{F}\mskip 1.5mu(u-\bar{u})\,\epsilon^{ll^{\prime}}\,,$
$\displaystyle\hskip 10.00002pt\Bigl{[}{}^{1}T_{g\to\delta q\mskip
1.5mu\delta\bar{q}}^{ll^{\prime}}(u)\Bigr{]}^{jj^{\prime}}$
$\displaystyle=-2T_{F}\mskip
1.5muu\bar{u}\;\delta^{ll^{\prime}}\delta^{jj^{\prime}}\,,$ (316)

where $\bar{u}=1-u$. For a quark splitting into a gluon and a quark, the
kernels are given by

$\displaystyle{}^{1}T_{q\to gq}^{ll^{\prime}}(u)$ $\displaystyle=C_{F}\mskip
1.5mu\frac{1+\bar{u}^{2}}{u}\,\delta^{ll^{\prime}}\,,$
$\displaystyle{}^{1}T_{q\to\Delta g\mskip 1.5muq}^{ll^{\prime}}(u)$
$\displaystyle=iC_{F}\mskip
1.5mu\frac{1+\bar{u}^{2}}{u}\,\epsilon^{ll^{\prime}}\,,$
$\displaystyle{}^{1}T_{q\to\Delta g\Delta q}^{ll^{\prime}}(u)$
$\displaystyle=C_{F}\mskip 1.5mu(1+\bar{u})\,\delta^{ll^{\prime}}\,,$
$\displaystyle{}^{1}T_{q\to g\Delta q}^{ll^{\prime}}(u)$
$\displaystyle=iC_{F}\mskip
1.5mu(1+\bar{u})\,\epsilon^{ll^{\prime}}\,,\phantom{\frac{1^{2}}{1}}$
$\displaystyle\left[{}^{1}T_{q\to\delta g\mskip
1.5muq}^{ll^{\prime}}(u)\right]^{jj^{\prime}}$ $\displaystyle=2C_{F}\mskip
1.5mu\frac{\bar{u}}{u}\,\tau^{ll^{\prime},jj^{\prime}}$ (317)

with $\tau^{ll^{\prime},jj^{\prime}}$ defined in (19). The kernels for an
antiquark splitting into a gluon and antiquark have the same form due to
charge conjugation invariance. Finally, for a gluon splitting into two gluons
one has

$\displaystyle{}^{1}T_{g\to gg}^{ll^{\prime}}(u)$
$\displaystyle=2C_{A}\left(\frac{u}{\bar{u}}+\frac{\bar{u}}{u}+u\bar{u}\right)\,\delta^{ll^{\prime}}\,,$
$\displaystyle{}^{1}T_{g\to\Delta g\Delta g}^{ll^{\prime}}(u)$
$\displaystyle=2C_{A}(2-u\bar{u})\,\delta^{ll^{\prime}}\,,\phantom{\frac{1}{2}}$
$\displaystyle{}^{1}T_{g\to g\Delta g}^{ll^{\prime}}(u)$
$\displaystyle=2iC_{A}\left(2\bar{u}+\frac{u}{\bar{u}}\right)\,\epsilon^{ll^{\prime}}\,,$
$\displaystyle\left[{}^{1}T_{g\to\delta g\mskip 1.5mu\delta
g}^{ll^{\prime}}(u)\right]^{jj^{\prime},kk^{\prime}}$
$\displaystyle=C_{A}\mskip
1.5muu\bar{u}\;\delta^{ll^{\prime}}\tau^{jj^{\prime},kk^{\prime}}\,,$
$\displaystyle\left[{}^{1}T_{g\to g\mskip 1.5mu\delta
g}^{ll^{\prime}}(u)\right]^{kk^{\prime}}$ $\displaystyle=2C_{A}\mskip
1.5mu\frac{u}{\bar{u}}\,\tau^{ll^{\prime},kk^{\prime}}\,.$ (318)

All other kernels that cannot be obtained from those above by interchanging
$a_{1}$ and $a_{2}$ or by changing quarks into antiquarks in (7.4) are zero.
In appendix H we give the matching equations that follow from these results
for the scalar and pseudoscalar DTMDs or DPDFs defined in appendix F.

Each of the above kernels corresponds to exactly one Feynman graph, so that
the results for other colour channels can be obtained by changing the overall
factor as

$\displaystyle{}^{R}T_{a_{0}\to a_{1}a_{2}}$ $\displaystyle=c_{a_{0}\to
a_{1}a_{2}}(R)\;{}^{1}T_{a_{0}\to a_{1}a_{2}}$ (319)

with

$\displaystyle c_{q\to gq}(A)$ $\displaystyle=-\frac{N}{\sqrt{2}}\,,$
$\displaystyle c_{q\to gq}(S)$ $\displaystyle=\sqrt{\frac{N^{2}-4}{2}}\,,$
$\displaystyle c_{g\to q\bar{q}}(8)$
$\displaystyle=-\frac{1}{\sqrt{N^{2}-1}}\,,$ $\displaystyle c_{g\to gg}(A)$
$\displaystyle=-\frac{\sqrt{N^{2}-1}}{2}\,,$ $\displaystyle c_{g\to gg}(S)$
$\displaystyle=-c_{g\to gg}(A)\,,$ $\displaystyle c_{g\to gg}(D)$
$\displaystyle=0\,,$ $\displaystyle c_{g\to gg}(27)$
$\displaystyle=-\sqrt{3}\,.\phantom{\sqrt{\frac{N^{2}-4}{2}}}$ (320)

For the colour representations $D$ and $27$, we have given the results for
$N=3$. Notice that in all cases except $c_{g\to q\bar{q}}(8)$ and $c_{g\to
gg}(D)$, the colour non-singlet channel is enhanced over the singlet one for
SU(3).

##  8 Summary

We have performed a systematic analysis of double parton scattering for
measured transverse momenta, extending the framework for factorisation and
resummation formulated in Collins:2011zzd  for single parton scattering with
colourless particles in the final state. A major challenge is the description
of soft gluon exchange and the associated rapidity evolution, which has a much
richer colour structure in DPS than in SPS. We handle this structure by
projecting suitable pairs of partons onto definite colour representations.

Based on perturbative results up to two loops Vladimirov:2016qkd , we have
proposed equation (37) for the rapidity dependence of the soft factor in DPS,
which is a Collins-Soper equation for matrices in colour space. Assuming the
validity of that equation, we have given a concise definition of DPDs in (3.3)
and derived their general properties. Applied to SPS, our construction
provides an alternative form for the square root factor in the definition of
single-parton TMDs by Collins Collins:2011zzd . Transverse momentum dependent
DPDs (DTMDs) follow the evolution equations (124) and (128), which can be
solved in closed analytic form as given in (5.1). The rapidity evolution of
DTMDs mixes different colour channels, and its solution involves a matrix
exponential in the space of colour representations. One-loop expressions of
the evolution kernels and anomalous dimensions are given in section 7.2 for
all parton combinations and colour channels.

In collinear factorisation, i.e. when the transverse momenta of the final-
state particles are integrated over, the soft factor for DPS simplifies
considerably: it becomes diagonal in colour representation space and
independent of the parton type. This can be shown using only colour algebra
and the fact that pairs of Wilson lines associated with initial and final
state partons cancel as $WW^{\dagger}=\mathbb{1}$. Important consequences of
this result are the simple structure (4.5) of the cross section for DPS
processes producing coloured particles, as well as the fact that evolution in
rapidity does not mix different colour channels for collinear DPDs (called
DPDFs here). The corresponding evolution equations can be found in (140) and
(141), and an analytic expression that makes the rapidity dependence explicit
is provided in (5.2). At one-loop accuracy, the evolution kernels in colour
non-singlet channels are related in a simple way to the usual DGLAP kernels,
as specified in (7.3.5).

If the size $q_{T}$ of the measured transverse momenta in a DPS process is
large compared with nonperturbative scales $\Lambda$ (but still small compared
with the scale $Q$ of the hard-scattering processes), one can use TMD
factorisation with DTMDs expressed in terms of perturbative kernels and
simpler hadronic matrix elements. There are two regimes for this short-
distance matching, depending on the relative size of the transverse distance
$\boldsymbol{y}$ between the two partons and the distances
$\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$ that are conjugate to the
measured transverse momenta in the cross section formula.

In the first regime, when $\boldsymbol{y}$ is much larger than
$\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$, the matching is very similar to
the familiar matching of single parton TMDs onto single parton PDFs. The
resulting expressions are given in (6.2.2) and (6.2.2). Apart from a
suppression factor for colour non-singlet DPDFs, which is controlled by their
rapidity evolution kernel ${}^{R}\mskip-1.5muJ(\boldsymbol{y})$, we obtain a
product of Sudakov exponentials and of matching kernels (one for each parton)
with DPDFs. At one-loop order, the matching kernels
${}^{R}\mskip-1.5muC(x,\boldsymbol{z})$ are due to ladder-type graphs and
coincide with the ones for single-parton TMDs in the colour singlet channel
$R=1$. Most of these kernels have been given in the literature before, but we
have recomputed them in section 7.3 and extended the calculation to provide
results for all combinations of parton type, polarisation and colour needed
for DPS. We use the HVBM scheme for defining $\gamma_{5}$ and the $\epsilon$
tensor in dimensional regularisation, but point out a way to avoid these
quantities before taking the number of dimensions to $D=4$.

In the second regime, when $\boldsymbol{y}$ is of the same size as
$\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$, one can match the DTMDs onto
collinear distributions of twist two (i.e. usual PDFs), twist three or twist
four. We neglect collinear twist-three distributions since in an unpolarised
proton they are chiral-odd and hence expected to play a minor role in the
regime of small parton momentum fractions typical of DPS processes. The
matching onto twist-four distributions is presently only known at lowest
order, where it is trivial and given in (194). The contribution that matches
onto ordinary PDFs corresponds to the splitting of one parton into two. At
leading order it has the form (193) with splitting kernels given in section
7.4. An important point is that in the DPS cross section, terms matching onto
collinear distributions of different twist come with different powers of
$\Lambda/q_{T}$ and of $\alpha_{s}$ as specified in (195).

To combine the expressions for matching at large or small $\boldsymbol{y}$, we
adapt the subtraction formalism of Collins:2011zzd . In this formalism, double
counting is avoided by a subtraction term in the cross section. This term can
be easily obtained from the small-$\boldsymbol{y}$ or the
large-$\boldsymbol{y}$ expression without further computation, but it requires
some care when choosing scales as discussed in section 6.4. The same
subtraction formalism is used to combine the cross sections for DPS, SPS and
their interference as specified in sections 5.3 and 6.5. In the regime where
$\Lambda\ll q_{T}\ll Q$ one needs a limited set of nonperturbative quantities
for computing the cross section: the DPDFs
${}^{R}\mskip-1.5muF(x_{i},\boldsymbol{y})$ and their rapidity evolution
kernel ${}^{R}\mskip-1.5muJ(\boldsymbol{y})$, collinear twist-four
distributions (which arise from DPDs in the short-distance limit) and ordinary
single-parton PDFs. The colour octet kernel
${}^{8}\mskip-1.5muJ(\boldsymbol{y})$ turns out to be equal to the Collins-
Soper kernel for single-gluon TMDs, providing a surprising connection between
two rather different areas of factorisation.

In the multi-scale regime $\Lambda\ll q_{T}\ll Q$ we thus have a significantly
increased predictive power of the theory. Rather than a huge set of
transverse-momentum dependent distributions, the TMD cross section involves
the same nonperturbative functions as the cross section integrated over
$q_{T}$. All other ingredients in the factorisation formula are of
perturbative nature. In most (although not all) cases they are known at least
to the first nontrivial order in $\alpha_{s}$. The results presented in this
work provide a starting point for phenomenological investigations, for
instance of electroweak diboson production with measured transverse boson
momenta. An important task will be to assess which types of correlations and
which regions of $\boldsymbol{y}$ are important in a given process and
kinematical setting.

Before concluding, we wish to discuss the question of scheme dependence in the
treatment of large gluon rapidities. Matrix elements of Wilson lines along
lightlike directions typically have rapidity divergences, which originate from
regions of gluon momenta with infinite rapidity. In the present paper, we
follow the approach of Collins Collins:2011zzd  and use Wilson line directions
away from the light cone to regulate these divergences. The regulated Wilson
lines then depend on rapidity variables. In the DPD construction laid out in
section 3.3, the parameter $\zeta$ is defined in terms of a central rapidity
$Y_{C}$ that specifies the range of gluon rapidities effectively included in
the DPD. Several other schemes for regulating rapidity divergences have been
proposed and used in the literature, see GarciaEchevarria:2011rb ;
Echevarria:2015byo ; Echevarria:2016scs ; Chiu:2011qc ; Chiu:2012ir ;
Li:2016axz  and Becher:2010tm . Quantities that depend on a regulator variable
will in general differ between schemes. By contrast, quantities that enter a
cross section formula in the same way in two schemes must be equal. This holds
in particular for TMDs, and hence also for the associated Collins-Soper
kernels and matching coefficients onto PDFs. Discussions of equivalence
between specific schemes can be found in Collins:2012uy ; Echevarria:2012js ;
Collins:2017oxh .

It follows that our results in sections 4 to 7 are independent of the rapidity
regulator scheme.161616A few equations in these sections are obviously
specific to the regulator we use, such as the definition (125) of $\zeta$ in
terms of $Y_{C}$. Of course, this only holds if one defines the variables
$\zeta$ and $\bar{\zeta}$ in a consistent way across different schemes, see
for instance appendix B here and section 6.1 in Vladimirov:2017ksc . We have
checked by explicit calculation that the one-loop matching coefficients
presented in section 7.3 remain the same when one uses the $\delta$ regulator
specified in GarciaEchevarria:2011rb . By contrast, our arguments in section 3
are specifically formulated for Wilson lines with finite rapidities and would
need to be adapted to other schemes, using an appropriate translation between
the regulator variables as given in appendix B for a particular case.

As already noted, we absorb the soft factor into the DPDs, corresponding to
what is done for SPS in several TMD factorisation approaches Collins:2011zzd ;
GarciaEchevarria:2011rb . A different route was taken in Manohar:2012jr ,
where the soft factor explicitly appears in the final DPS cross section
formula. Such an approach is often taken in the SCET literature in situations
when the soft factor can be computed in perturbation theory (e.g. for single
Drell-Yan production at large $q_{T}$). In DPS, however, one needs the soft
factor at large distances $\boldsymbol{y}$. Including it into the DPDs
therefore limits the number of nonperturbative functions required to compute
the cross section.

Given the results of the present paper and the arguments regarding Glauber
gluon exchange in Diehl:2015bca , we consider that TMD factorisation for DPS
is now established at the same level of rigour as for single hard scattering,
apart from the following issues. $(i)$: As already mentioned, there is no
complete all-order proof of the properties we assumed for the soft factor in
section 3.2, but significant process in this direction has been made in
Vladimirov:2017ksc . $(ii)$: There is no all-order proof of the nonabelian
Ward identities required for factorising soft gluon exchange between left and
right moving partons into the soft matrix. $(iii)$: As noted below (253), we
have omitted Wilson line self interactions in our calculation, following the
logic discussed for SPS in chapter 13.7 of Collins:2011zzd . A more thorough
understanding of this issue, in particular for DPS, has not yet been achieved.
Progress on any of these three points will constitute a further step towards a
rigorous treatment of double parton scattering in QCD.

###### Acknowledgements.

We gratefully acknowledge discussions with Alessandro Bacchetta, John Collins,
Tom van Daal, Miguel G. Echevarría, Jonathan Gaunt, Zhongbo Kang, Piet
Mulders, Riccardo Nagar, Alexei Prokudin, Ted Rogers, Maximilian Stahlhofen,
Iain Stewart, Frank Tackmann, Alexey Vladimirov, Werner Vogelsang and Wouter
Waalewijn. TK is supported by the European Community under the “Ideas” program
QWORK (contract 320389). He acknowledges the hospitality of the Munich
Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of
excellence “Origin and Structure of the Universe”. All figures were made using
JaxoDraw Binosi:2003yf ; Binosi:2008ig . For calculations we have used the
FeynCalc package Mertig:1990an ; Shtabovenko:2016sxi , the ColorMath package
Sjodahl:2012nk  and FORM Vermaseren:2000nd ; Kuipers:2012rf .

##  Appendix A Matrix manipulations for the soft factor

Referring to section 3.2, we now show that the set of properties 1a – 1c for
the soft factor is equivalent with properties 2a – 2b. As a corollary, we find
that the soft matrix $S(Y)$ is positive definite for arbitrary large $Y$ if
these properties hold.

##### Properties 2a – 2b follow from properties 1a – 1c.

The matrix $s(0)$ is nonsingular by assumption (property 1b), so that one can
define a matrix $\widehat{K}$ by

$\displaystyle\widehat{K}s(0)$ $\displaystyle=s(0)K\,.$ (321)

The differential equation (35) (property 1a) is thus solved by

$\displaystyle s(Y)$ $\displaystyle=s(0)\mskip
1.5mue^{YK}=e^{Y\widehat{K}}\mskip 1.5mus(0)$ (322)

for all $Y$. We now define

$\displaystyle\widehat{S}(Y)$ $\displaystyle=s(Y/2)\,s^{\dagger}(Y/2)\,.$
(323)

Using (36) (property 1c) for $Y_{0}=Y/2$, we see that $\widehat{S}(Y)$
approximates $S(Y)$ for $Y\gg 1$.

For any complex vector $a$ one has

$\displaystyle a^{\dagger}\,\widehat{S}(Y)\mskip 1.5mua$
$\displaystyle=b^{\dagger}\mskip 1.5mub\,,$ (324)

with $b=s^{\dagger}(Y/2)\,a$. The expression in (324) is obviously positive or
zero, and since $s^{\dagger}(Y/2)$ is nonsingular by assumption (property 1b),
it is zero only for $a=0$. Therefore $\widehat{S}(Y)$ is positive definite for
all $Y$. This implies property 2b, and it also implies that $S(Y)$ is positive
for all $Y$ where it can be approximated by $\widehat{S}(Y)$.

Multiplying (321) with $s^{\dagger}(0)$ on the right, and using that $K$ is
Hermitian, we obtain

$\displaystyle\widehat{K}\widehat{S}(0)=\widehat{S}(0)\widehat{K}^{\dagger}\,.$
(325)

We thus finally have

$\displaystyle\widehat{S}(Y)$ $\displaystyle=e^{Y\widehat{K}/2}\mskip
1.5mu\widehat{S}(0)e^{Y\widehat{K}^{\dagger}/2}=e^{Y\widehat{K}}\widehat{S}(0)\,,$
(326)

where we used (322) and (323) in the first step and (325) in the second one.
This implies property 2a for $Y\gg 1$, where $\widehat{S}(Y)$ approximates
$S(Y)$.

##### Properties 1a – 1c follow from properties 2a – 2b.

We start by defining the matrix $\widehat{S}(Y)$ for all $Y$ by

$\displaystyle\widehat{S}(Y)$
$\displaystyle=\lim_{Y_{0}\to\infty}e^{(Y-Y_{0})\widehat{K}}S(Y_{0})\,.$ (327)

This obviously satisfies

$\displaystyle\frac{\partial}{\partial Y}\,\widehat{S}(Y)$
$\displaystyle=\widehat{K}\mskip 1.5mu\widehat{S}(Y)$ (328)

exactly. For $Y,Y_{0}\gg 1$ one has $S(Y)=e^{(Y-Y_{0})\widehat{K}}S(Y_{0})$ by
virtue of (37) (property 2a), which implies $\widehat{S}(Y)=S(Y)$ for $Y\gg
1$.

Because $S(Y)$ is Hermitian according to (3), it follows from (37) that
$\widehat{K}S(Y)=S(Y)\mskip 1.5mu\widehat{K}^{\dagger}$ for $Y\gg 1$. Using
this, one deduces from (327) that $\widehat{S}(Y)$ is Hermitian as well. With
(328) we then have

$\displaystyle\widehat{K}\widehat{S}(Y)=\widehat{S}(Y)\widehat{K}^{\dagger}$
(329)

for all $Y$. The solution of (328) can thus be written as

$\displaystyle\widehat{S}(Y)$ $\displaystyle=e^{(Y-Y_{1})\widehat{K}}\mskip
1.5mu\widehat{S}(Y_{1})=e^{(Y-Y_{1}-Y_{0})\widehat{K}}\mskip
1.5mu\widehat{S}(Y_{1})\mskip 1.5mue^{Y_{0}\widehat{K}^{\dagger}}$ (330)

for arbitrary $Y_{0}$ and $Y_{1}$.

We now recall the definition of the square root ${M}^{1/2}$ of a positive
definite matrix $M$. One can write $M=V^{\dagger}DV$, where $V$ is unitary and
$D$ is diagonal with positive entries $d_{ii}$. One defines ${D}^{1/2}$ as the
diagonal matrix with entries $\sqrt{d_{ii}}$, and furthermore
$\sqrt{M}=V^{\dagger}{D}^{1/2}V$. If $M$ is real then $V$ can be taken as an
orthogonal matrix, so that ${M}^{1/2}$ is real as well. Since
$\widehat{S}(Y_{1})$ is positive definite for some $Y_{1}$ (property 2b), one
can use (330) with $Y_{0}=(Y-Y_{1})/2$ and write
$\widehat{S}(Y)=t(Y)\,t^{\dagger}(Y)$, where
$t(Y)=e^{(Y-Y_{1})\widehat{K}/2}\mskip
1.5mu\bigl{[}\widehat{S}(Y_{1})\bigr{]}{}^{1/2}$ is a nonsingular matrix. With
the same argument as given below (323), it then follows that $\widehat{S}(Y)$
is positive for all $Y$.

We can now define

$\displaystyle s(Y)$ $\displaystyle=e^{Y\widehat{K}}\mskip
1.5mu\bigl{[}\widehat{S}(0)\bigr{]}{}^{1/2}\mskip 1.5muU\,,$ (331)

where $U$ is a unitary matrix. $s(Y)$ is nonsingular (property 1b) and
satisfies the decomposition (36) (property 1c) by virtue of (330) with
$Y_{1}=0$. Defining

$\displaystyle K$ $\displaystyle=U^{\dagger}\mskip
1.5mu\bigl{[}\widehat{S}(0)\bigr{]}{}^{-1/2}\widehat{K}\,\bigl{[}\widehat{S}(0)\bigr{]}{}^{1/2}\mskip
1.5muU\,,$ (332)

one readily finds that $K^{\dagger}=K$ and

$\displaystyle s(Y)$
$\displaystyle=\bigl{[}\widehat{S}(0)\bigr{]}{}^{1/2}\mskip 1.5muU\mskip
1.5mue^{YK}\,,$ (333)

which implies (35) (property 1a). The definitions (331) and (332) are
consistent with (321) above. We note that a more general definition of $s(Y)$
and $K$ is obtained by replacing $\widehat{S}(0)\to\widehat{S}(Y_{1})$ in
(331) to (333), as well as $Y\to(Y-Y_{1}/2)$ in the exponentials there.

##### Restrictions on $U$.

In principle, the matrix $U$ can be chosen freely as a smooth function of
$\boldsymbol{z}_{1}$, $\boldsymbol{z}_{2}$ and $\boldsymbol{y}$, but
independent of $Y$ and of the renormalisation scale $\mu$. We restrict
ourselves to choices such that $U=1$ for $\boldsymbol{y}=\boldsymbol{0}$.
Since $S(Y)$ is diagonal in that case (see section 4.4), the same then holds
for $\widehat{S}(Y)$, $s(Y)$ and for $K=\widehat{K}$.

We furthermore impose a number of restrictions on $U$ in different partonic
channels. We thus demand that

$\displaystyle U_{qq}=$ $\displaystyle U_{\bar{q}\bar{q}}^{*}\,,$
$\displaystyle U_{q\bar{q}}$ $\displaystyle=U_{\bar{q}q}^{*}\,,$
$\displaystyle U_{qg}$ $\displaystyle=U_{\bar{q}g}^{*}\,,$ $\displaystyle
U_{gg}$ $\displaystyle=U_{gg}^{*}\,,$ (334)

which is necessary to obtain the corresponding relations (41) for
$s_{a_{1}a_{2}}$ from those for $S_{a_{1}a_{2}}$. As discussed in section 4.3,
we also require that ${}^{RR^{\prime}}\mskip-1.5muU$, defined in analogy to
(81), is a real valued matrix in representation space and satisfies

$\displaystyle{}^{RR^{\prime}}\mskip-1.5muU_{qq}(\boldsymbol{z}_{1},\boldsymbol{z}_{2},\boldsymbol{y})$
$\displaystyle={}^{RR^{\prime}}\mskip-1.5muU_{q\bar{q}}(\boldsymbol{z}_{1},-\boldsymbol{z}_{2},\boldsymbol{y})$
$\displaystyle={}^{RR^{\prime}}\mskip-1.5muU_{\bar{q}q}(-\boldsymbol{z}_{1},\boldsymbol{z}_{2},\boldsymbol{y})={}^{RR^{\prime}}\mskip-1.5muU_{\bar{q}\bar{q}}(-\boldsymbol{z}_{1},-\boldsymbol{z}_{2},\boldsymbol{y})\,.$
(335)

We note that the trivial choice $U_{a_{1}a_{2}}=\mathbb{1}$ satisfies all
requirements just listed and gives
${}^{RR^{\prime}}\mskip-1.5muU_{a_{1}a_{2}}=\delta_{RR^{\prime}}$.

##  Appendix B Comparison with the $\delta$ regulator scheme

As emphasised in section 8, our treatment of soft and collinear factors uses
Wilson line directions away from the light cone to regulate rapidity
divergences, as was done by Collins in Collins:2011zzd . The regulated Wilson
lines then depend on large but finite rapidities.

The two-loop calculation of the DPS soft factor in Vladimirov:2016qkd  was
performed using the $\delta$ regulator scheme specified in Echevarria:2015byo
; Echevarria:2016scs . Given the particular relevance of this calculation for
our assumptions on the soft factor in section 3.2, we now compare the
variables in that scheme with the ones in the scheme of Collins. For
simplicity we perform the comparison for the SPS soft factor. It is easy to
extend the following arguments to DPS, including the appropriate rescaling of
$\zeta$, see (61) and (125). We note that there are two earlier variants of
the $\delta$ regulator, described in GarciaEchevarria:2011rb  and
Echevarria:2012js , and compared with the Collins scheme in Collins:2012uy
and Echevarria:2012js , respectively.

In the $\delta$ regulator scheme of Echevarria:2015byo ; Echevarria:2016scs ,
Wilson lines are taken along lightlike paths but modified by an exponential
damping at large distances $z^{-}$ ($z^{+}$). This damping is controlled by a
parameter $\delta^{+}$ ($\delta^{-}$), which transforms like the plus (minus)
component of a vector under longitudinal boosts. One has a correspondence of
variables

$\displaystyle Y_{R}$
$\displaystyle\,\leftrightarrow\,\log\frac{\mu}{\delta^{-}}\,,$ $\displaystyle
Y_{L}$ $\displaystyle\,\leftrightarrow\,\log\frac{\delta^{+}}{\mu}\,.$ (336)

Due to boost invariance, the soft factor $S$ only depends on
$Y_{R}-Y_{L}\leftrightarrow-\log(\delta^{+}\delta^{-}/\mu^{2})$. The
correspondence indicated by $\leftrightarrow$ is valid when taking
derivatives, so that $\mu$ in (336) could be replaced by another quantity that
has dimension of mass and is boost invariant. In our formalism, the soft
factor can be split into

$\displaystyle S(Y_{R}-Y_{L})$
$\displaystyle=\sqrt{S(2Y_{R}-2Y_{C})}\,\sqrt{S(2Y_{C}-2Y_{L})}\,,$ (337)

and the same splitting is performed with the $\delta$ regulator, with the
following correspondence of variables:

$\displaystyle 2(Y_{R}-Y_{C})$
$\displaystyle\,\leftrightarrow\,\log\frac{\mu^{2}}{(\delta^{-})^{2}}\,\frac{2(\bar{k}^{-})^{2}}{\bar{\zeta}\rule{0.0pt}{9.47217pt}}\,,$
$\displaystyle 2(Y_{C}-Y_{L})$
$\displaystyle\,\leftrightarrow\,\log\frac{\mu^{2}}{(\delta^{+})^{2}}\,\frac{2(k^{+})^{2}}{\zeta}\,.$
(338)

Here the momentum components $k^{+}$ and $\bar{k}^{-}$ correspond to the two
partons initiating a hard scattering, and we have normalised the rapidity
parameters as $\zeta\bar{\zeta}=(2k^{+}\bar{k}^{-})^{2}$, in accordance to the
convention (61) used in the present paper.171717The normalisation of $\zeta$
and $\bar{\zeta}$ is not relevant for the Collins Soper equation, but it does
matter for the argument of the anomalous dimensions $\gamma_{F,a}$. On both
sides of (338), the two expressions add up to twice the argument of the
original soft factor in the relevant formalism. After the square roots of the
soft factor are combined with unsubtracted collinear matrix elements (see
(60)), one removes the rapidity regulator in either scheme, taking
$\delta^{-}\to 0$ and $\delta^{+}\to 0$ or $Y_{R}\to\infty$ and
$Y_{L}\to-\infty$. The resulting distributions depend on the rapidity variable
$\zeta$ or $\bar{\zeta}$.

##  Appendix C Combining SPS and DPS for different hard scales

In this appendix, we discuss the generalisation of the scale setting discussed
in section 6.5 to the case where the scales of the two DPS processes
($\mu_{1}$ and $\mu_{2}$) and of the SPS process ($\mu_{h}$) are different. In
the four contributions (6.5) to the cross section, we now take

$\displaystyle\sigma_{\text{SPS}}:$ $\displaystyle
f(x_{1}+x_{2},\boldsymbol{Z})$ $\displaystyle\sim
E_{2}(\mu_{h};\mu_{Z})\,f(\mu_{Z})\,,$ $\displaystyle\sigma_{\text{DPS/SPS}}:$
$\displaystyle D(x_{i},\boldsymbol{y}_{+},\boldsymbol{Z})$ $\displaystyle\sim
E_{3}(\mu_{1},\mu_{2}|\mu_{h};\mu_{Z})\,U(\mu_{Z})\,f(\mu_{Z})\,,$
$\displaystyle\sigma_{\text{SPS/DPS}}:$ $\displaystyle
D(x_{i},\boldsymbol{y}_{-},\boldsymbol{Z})$ $\displaystyle\sim
E_{3}(\mu_{h}|\mu_{1},\mu_{2};\mu_{Z})\,U^{*}(\mu_{Z})\,f(\mu_{Z})\,,$
$\displaystyle\sigma_{\text{DPS}}:$ $\displaystyle
F(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$ $\displaystyle\sim
E_{4}(\mu_{1},\mu_{2}|\mu_{1},\mu_{2};\mu_{Z})\,U(\mu_{Z})\,U^{*}(\mu_{Z})\,f(\mu_{Z})+F_{\text{int}}(\mu_{1},\mu_{2})\,.$
(339)

In $E_{3}$ and $E_{4}$, we now allow different scales for all three or four
parton legs (see the remarks at the ends of sections 3.4 and 5.3). Scale
arguments in the amplitude and in its complex conjugate are separated by a
vertical bar “$|$”. If only a single scale is given, such as $\mu_{Z}$ in (C),
then all scales are taken equal. In the subtraction terms for the SPS/DPS
interference, we now take

$\displaystyle\sigma_{\text{DPS/SPS},y_{+}\to 0}:D$ $\displaystyle\sim
E_{3}(\bar{\mu}_{1},\bar{\mu}_{2}|\mu_{h};\hat{\mu}_{h})\,U(\hat{\mu}_{h})\,E_{2}(\hat{\mu}_{h};\mu_{Z})\,f(\mu_{Z})$
$\displaystyle\qquad\mbox{with\leavevmode\nobreak\ }\hat{\mu}_{h}=p(\nu\mskip
1.5mu|\boldsymbol{y}_{+}|;\mu_{Z},\mu_{h})\mbox{ \leavevmode\nobreak\
and\leavevmode\nobreak\ }\bar{\mu}_{i}=p(\nu\mskip
1.5mu|\boldsymbol{y}_{+}|;\mu_{i},\mu_{h})\,,$
$\displaystyle\sigma_{\text{SPS/DPS},y_{-}\to 0}:D$ $\displaystyle\sim
E_{3}(\mu_{h}|\bar{\mu}_{1},\bar{\mu}_{2};\hat{\mu}_{h})\,U^{*}(\hat{\mu}_{h})\,E_{2}(\hat{\mu}_{h};\mu_{Z})\,f(\mu_{Z})$
$\displaystyle\qquad\mbox{with\leavevmode\nobreak\ }\hat{\mu}_{h}=p(\nu\mskip
1.5mu|\boldsymbol{y}_{-}|;\mu_{Z},\mu_{h})\mbox{ \leavevmode\nobreak\
and\leavevmode\nobreak\ }\bar{\mu}_{i}=p(\nu\mskip
1.5mu|\boldsymbol{y}_{-}|;\mu_{i},\mu_{h})$ (340)

instead of (6.5). The profile scales $\hat{\mu}$ are the same as in the
single-scale case of section 6.5, whilst $\bar{\mu}$ interpolates between the
different hard scales. The relation (215) between the hard-scattering
amplitudes for SPS and DPS holds for equal scales of $\alpha_{s}$ in all its
terms. In all subtraction terms, the scale of $\alpha_{s}$ in the DPS
amplitudes $H_{\alpha_{1}\beta_{1}}$ and $H_{\alpha_{2}\beta_{2}}$ should be
taken as $\bar{\mu}_{i}$. This is for instance relevant for $gg\to H$, which
involves the strong coupling already at leading order. With the choices just
specified, one finds that the limiting relation (214) for the subtraction term
still holds.

Generalising (6.5), the subtraction terms for DPS can be taken as follows:

$\displaystyle\sigma_{\text{DPS},y_{-}\to 0}:$ $\displaystyle F$
$\displaystyle\sim
E_{4}(\bar{\mu}_{1},\bar{\mu}_{2}|\mu_{1},\mu_{2};\hat{\mu}_{1},\hat{\mu}_{2}|\hat{\mu}_{h},\hat{\mu}_{h})\,U^{*}(\hat{\mu}_{h})\,E_{3}(\hat{\mu}_{1},\hat{\mu}_{2}|\hat{\mu}_{h};\mu_{Z})\,U(\mu_{Z})\,f(\mu_{Z})\,,$
$\displaystyle\sigma_{\text{DPS},y_{+}\to 0}:$ $\displaystyle F$
$\displaystyle\sim
E_{4}(\mu_{1},\mu_{2}|\bar{\mu}_{1},\bar{\mu}_{2};\hat{\mu}_{h},\hat{\mu}_{h}|\hat{\mu}_{1},\hat{\mu}_{2})\,U(\hat{\mu}_{h})\,E_{3}(\hat{\mu}_{h}|\hat{\mu}_{1},\hat{\mu}_{2};\mu_{Z})\,U^{*}(\mu_{Z})\,f(\mu_{Z})\,,$
$\displaystyle\sigma_{\text{DPS},y_{\pm}\to 0}:$ $\displaystyle F$
$\displaystyle\sim
E_{4}(\bar{\mu}_{1},\bar{\mu}_{2}|\bar{\mu}_{1},\bar{\mu}_{2};\hat{\mu}_{h})\,U(\hat{\mu}_{h})\
U^{*}(\hat{\mu}_{h})\,E_{2}(\hat{\mu}_{h};\mu_{Z})\,f(\mu_{Z})$ (341)

with

$\displaystyle\hat{\mu}$
$\displaystyle=p\bigl{(}\nu\min\\{|\boldsymbol{y}_{+}|,|\boldsymbol{y}_{-}|\\};\mu_{Z},\mu_{h}\bigr{)}$
and $\displaystyle\bar{\mu}_{i}$
$\displaystyle=p\bigl{(}\nu\min\\{|\boldsymbol{y}_{+}|,|\boldsymbol{y}_{-}|\\};\mu_{i},\mu_{h}\bigr{)}\,.$
(342)

For $|\boldsymbol{y}_{+}|\gg 1/\nu$ and $|\boldsymbol{y}_{-}|\gg 1/\nu$ one
finds the limiting behaviour in (218), whereas all three terms in (C) tend to
the 1v1 contribution to DPS if one or both of the transverse distances is of
order $1/\nu$. Overall, we thus find that the subtraction formalism can be
adapted to work in the multi-scale case.

##  Appendix D Feynman rules

In this appendix we give the Feynman rules used in our calculation of the
matching coefficients in section 7.3. We compute cut graphs and thus need in
particular the rules for propagators and vertices to the right of the cut
(corresponding to the complex conjugate amplitude). This requires special care
for the three-gluon vertex. A derivation of these rules can for instance be
found in Belitsky:1997ay  for the QCD Lagrangian (but not for eikonal lines).
The rules given here correspond to the conventions in the chapters 3 and 7 of
Collins:2011zzd . In particular, the coupling constant is defined such that
the covariant derivative reads
$D_{\mu}=\partial_{\mu}+igA^{a}{}_{\\!\\!\\!\mu}\,t^{a}$.

The Feynman rules for the propagators and vertices arising from the QCD
Lagrangian in Feynman gauge can be found in figure 12, where we use the
notation

$\displaystyle V^{\mu\nu\rho}(p,q,r)$
$\displaystyle=(p-q)^{\rho}g^{\mu\nu}+(q-r)^{\mu}g^{\nu\rho}+(r-p)^{\nu}g^{\rho\mu}\,.$
(343)

The four-gluon vertex is not given here since it does not appear in our
calculations.

Figure 12:  Feynman rules for propagators and vertices in cut graphs. We use
Feynman gauge and set the quark mass to zero. $V^{\mu\nu\rho}(p,q,r)$ is
defined in (343). With these rules for cut quark and antiquark lines, there is
no minus sign for a closed fermion loop going across the cut.

Figure 13:  Feynman rules involving eikonal lines along the direction $v$. The
colour indices $r$ and $s$ refer to either the fundamental or the adjoint
representation, whereas $j$ and $k$ are colour triplet and $a$, $b$ and $c$
are colour octet indices.

Figure 14:  Feynman rule corresponding to the operator in a gluon distribution
when an eikonal line along $n$ is attached to the gluon. The Lorentz index $j$
is transverse w.r.t. the lightlike direction $n$. The rule for the graph
without the eikonal line is obtained by setting $p=q$.

The Feynman rules for the propagators and couplings involving eikonal lines
are given in figure 14. They arise from the expansion of Wilson line
operators, given by (16) and its analogue for the adjoint representation. We
use there the notation of Diehl:2011yj , where open and closed circles at the
ends of eikonal lines (in addition to arrows on and above them) were
introduced as a way to make the correspondence between graphs and mathematical
expressions unique. Let us briefly explain this.

First of all, the full circle indicates the (relative) past and the open
circle the (relative) future time direction when considering the path of the
Wilson line in space-time. This determines the sign of $i0$ in the eikonal
propagators. Secondly, the arrow _above_ the eikonal line propagator
corresponds to the momentum flow, just as for quark and gluon propagators, and
is necessary because the eikonal propagator is linear in the momentum carried
by the line. Finally, the arrow _on_ the eikonal line fixes the overall sign
of the propagator and determines the order of contraction of colour indices,
with matrix multiplication going against the direction of the arrow. This is
the same convention as for Dirac indices on fermion lines. For eikonal lines
in the fundamental colour representation, this arrow indicates the canonical
order of multiplication, with expressions involving Gell-Mann matrices $t_{a}$
rather than their transpose $(t_{a})^{T}$. For adjoint eikonal lines, the
order in which the indices are arranged can be chosen freely (writing e.g.
$f_{abc}$ or $-f_{acb}$).

The vertices for a gluon coupling to a fundamental or an adjoint eikonal line
are given in the third and fourth line of figure 14, respectively. In the last
row, the relative sign between the vertices to the left and the right of the
cut can be understood from the relation $f_{cab}=-i(T^{a})_{cb}$, where
$T^{a}$ acts as a generator of the colour group. The arrow on the line is
required to distinguish the role of the indices $c$ and $b$, given that
$f_{cab}$ is antisymmetric in its indices.181818We note that in figure 7.10 of
Collins:2011zzd  a minus sign is missing in the coupling of a gluon to an
adjoint eikonal line on the right of the cut. Otherwise, we agree with the
Feynman rules given there.

In figure 14 we give the Feynman rule for the operator in a gluon distribution
when an eikonal line is attached to the gluon. This rule was derived in
section 7.6 of Collins:2011zzd , see also Collins:1981uw ; Collins:2008sg . It
can be obtained by expanding the operator $n_{\mu}\mskip 1.5muF^{\mu
j}(z)\,W(z,n)$ in $g$. The Lorentz index $j$ in the gluon field strength is
restricted to be transverse to the lightlike vector $n$. Note that this rule
corresponds to the gluon operator (17) only if the Wilson line is along
$v_{L}=n$.

##  Appendix E Useful integrals

In this appendix, we give a number of Fourier integrals that are useful for
the calculations in section 7.3. We start with

$\displaystyle\int\frac{d^{2-2\epsilon}\boldsymbol{k}}{\left(\boldsymbol{k}^{2}\right)^{\alpha}}\,e^{i\boldsymbol{k}\boldsymbol{z}}$
$\displaystyle=\left(\frac{\boldsymbol{z}^{2}}{4\pi}\right)^{\epsilon+\alpha-1}\frac{\pi^{\alpha}\mskip
1.5mu\Gamma(1-\epsilon-\alpha)}{\Gamma(\alpha)},$ (344)
$\displaystyle\int\frac{d^{2-2\epsilon}\boldsymbol{k}}{\left(\boldsymbol{k}^{2}\right)^{\alpha}}\,e^{i\boldsymbol{k}\boldsymbol{z}}\,\log(\boldsymbol{k}^{2})$
$\displaystyle=\left(\frac{\boldsymbol{z}^{2}}{4\pi}\right)^{\epsilon+\alpha-1}\frac{\pi^{\alpha}\mskip
1.5mu\Gamma(1-\epsilon-\alpha)}{\Gamma(\alpha)}$
$\displaystyle\qquad\times\biggl{[}\mskip
1.5mu\log\frac{4}{\boldsymbol{z}^{2}}+\psi(\alpha)+\psi(1-\alpha-\epsilon)\mskip
1.5mu\biggr{]}\,,$ (345)

where $\psi(x)=\frac{d}{dz}\log\Gamma(z)$ is the digamma function. One can
obtain (345) by differentiation of (344) with respect to $\alpha$. A useful
relation for integrals involving tensors is given by

$\displaystyle\int\,\frac{d^{2-2\epsilon}{\boldsymbol{k}}}{\left(\boldsymbol{k}^{2}\right)^{\alpha}}\,e^{i{\boldsymbol{k}}\boldsymbol{z}}\,\boldsymbol{k}^{j}\boldsymbol{k}^{j^{\prime}}$
$\displaystyle=-\frac{\pi^{\alpha-1}\,\Gamma(3-\alpha-\epsilon)}{\Gamma(\alpha)}\left(\frac{\boldsymbol{z}^{2}}{4\pi}\right)^{\epsilon+\alpha-2}\biggl{[}\frac{\boldsymbol{z}^{j}\boldsymbol{z}^{j^{\prime}}}{\boldsymbol{z}^{2}}+\frac{\delta^{jj^{\prime}}}{2(\epsilon+\alpha-2)}\biggr{]}\,.$
(346)

We obtain (346) by rewriting $\boldsymbol{k}^{j}\boldsymbol{k}^{j^{\prime}}$
as
$-\partial^{2}/(\partial\boldsymbol{z}^{j}\partial\boldsymbol{z}^{j^{\prime}})$
on the l.h.s. of (346) and then performing these derivatives on the r.h.s. of
(344). From the above equations, it follows that

$\displaystyle\int\frac{d^{2-2\epsilon}\boldsymbol{k}}{\boldsymbol{k}^{2}}\,e^{i\boldsymbol{k}\boldsymbol{z}}\,\biggl{[}\frac{k^{jj^{\prime}}}{\boldsymbol{k}^{2}}+f(\epsilon)\,\delta^{jj^{\prime}}\biggr{]}$
$\displaystyle=-\pi\Gamma(1-\epsilon)\left(\frac{\boldsymbol{z}^{2}}{4\pi}\right)^{\epsilon}\biggl{[}\frac{z^{jj^{\prime}}}{\boldsymbol{z}^{2}}+\frac{f(\epsilon)}{\epsilon}\,\delta^{jj^{\prime}}\biggr{]}\,,$
(347)

where the traceless symmetric tensors $k^{jj^{\prime}}$ and $z^{jj^{\prime}}$
are defined by (256) and $f(\epsilon)$ is an arbitrary function. The integral
over $k^{jj^{\prime}}/\boldsymbol{k}^{4}$ is thus finite, whereas the one over
$f(\epsilon)/\boldsymbol{k}^{2}$ has an ultraviolet pole if $f(\epsilon)$ is
finite at $\epsilon=0$.

##  Appendix F Tensor decomposition of DPDs

In this appendix, we give the tensor decompositions of DTMDs and DPDFs in
terms of scalar and pseudoscalar functions. We emphasise that these
decompositions are only meant to be complete in the $2$ physical transverse
dimensions. In calculations that have to be done in $2-2\epsilon$ transverse
dimensions, one must either avoid the decompositions below (as we do) or
extend them appropriately.

DTMDs depend on the three transverse vectors $\boldsymbol{y}$,
$\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$, but one can construct the
tensors necessary for their decomposition from $\boldsymbol{y}$ and

$\displaystyle\tilde{\boldsymbol{y}}^{j}=\epsilon^{jj^{\prime}}\boldsymbol{y}^{j^{\prime}}$
(348)

together with the invariants $\delta^{jj^{\prime}}$ and
$\epsilon^{jj^{\prime}}$. This is because in 2 dimensions there are only two
linearly independent vectors. The choice of the above two vectors results in
the most straightforward correspondence between the DTMDs and the DPDFs, where
only the vectors $\boldsymbol{y}$ and $\tilde{\boldsymbol{y}}$ can be used.
The decomposition of the double quark distributions Diehl:2011yj  reads

$\displaystyle F_{qq}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=f_{qq}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})\,,$
$\displaystyle F_{\Delta q\Delta q}$ $\displaystyle=f_{\Delta q\Delta q}\,,$
$\displaystyle F_{\Delta q\mskip 1.5muq}$ $\displaystyle=g_{\Delta q\mskip
1.5muq}\,,$ $\displaystyle F_{q\Delta q}$ $\displaystyle=g_{q\Delta q}\,,$
$\displaystyle F_{\delta q\mskip 1.5muq}^{j}$
$\displaystyle=\tilde{\boldsymbol{y}}^{j}Mf_{\delta q\mskip
1.5muq}+\boldsymbol{y}^{j}Mg_{\delta q\mskip 1.5muq}\,,$ $\displaystyle
F_{q\mskip 1.5mu\delta q}^{j}$
$\displaystyle=\tilde{\boldsymbol{y}}^{j}Mf_{q\mskip 1.5mu\delta
q}+\boldsymbol{y}^{j}Mg_{q\mskip 1.5mu\delta q}\,,$ $\displaystyle F_{\delta
q\mskip 1.5mu\delta q}^{jj^{\prime}}$
$\displaystyle=\delta^{jj^{\prime}}f_{\delta q\mskip 1.5mu\delta
q}+2\tau^{jj^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}M^{2}f_{\delta q\mskip
1.5mu\delta q}^{t}$
$\displaystyle\quad+2\tau^{jj^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}}M^{2}g_{\delta
q\mskip 1.5mu\delta
q}^{s}+(\boldsymbol{y}^{j}\tilde{\boldsymbol{y}}^{j^{\prime}}-\tilde{\boldsymbol{y}}^{j}\boldsymbol{y}^{j^{\prime}})\,M^{2}g_{\delta
q\mskip 1.5mu\delta q}^{a}\,,$ (349)

where $M$ is the proton mass and the tensor
$\tau^{jj^{\prime}\\!,kk^{\prime}}$ is defined in (19). Note that one has
$\tau^{jj^{\prime}\\!,kk^{\prime}}\tau^{kk^{\prime}\\!,ll^{\prime}}=\tau^{jj^{\prime}\\!,ll^{\prime}}$.
We employ a shorthand notation where vectors $\boldsymbol{y}$ or
$\tilde{\boldsymbol{y}}$ appearing as an index of $\tau$ denote contraction,
i.e.
$\tau^{jj^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}=\tau^{jj^{\prime}\\!,kk^{\prime}}\,\boldsymbol{y}^{k}\boldsymbol{y}^{k^{\prime}}$
etc. $f$ and $g$ denote scalar and pseudo-scalar functions respectively. The
pseudoscalar functions are absent in the decomposition of DPDFs because one
cannot construct a pseudoscalar from $\boldsymbol{y}$ alone. Decompositions
analogous to (F) hold for quark-antiquark distributions and for distributions
of two antiquarks.

For quark-gluon distributions we write

$\displaystyle F_{qg}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=f_{qg}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})\,,$
$\displaystyle F_{\Delta q\Delta g}$ $\displaystyle=f_{\Delta q\Delta g}\,,$
$\displaystyle F_{\Delta q\mskip 1.5mug}$ $\displaystyle=g_{\Delta q\mskip
1.5mug}\,,$ $\displaystyle F_{q\Delta g}$ $\displaystyle=g_{q\Delta g}\,,$
$\displaystyle F_{\delta q\mskip 1.5mug}^{j}$
$\displaystyle=\tilde{\boldsymbol{y}}^{j}Mf_{\delta q\mskip
1.5mug}+\boldsymbol{y}^{j}Mg_{\delta q\mskip 1.5mug}\,,$ $\displaystyle
F_{q\mskip 1.5mu\delta g}^{jj^{\prime}}$
$\displaystyle=\tau^{jj^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}M^{2}f_{q\mskip
1.5mu\delta
g}+\tau^{jj^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}}M^{2}g_{q\mskip
1.5mu\delta g}\,,$ $\displaystyle F_{\delta q\mskip 1.5mu\delta
g}^{j,kk^{\prime}}$
$\displaystyle=-\tau^{\tilde{\boldsymbol{y}}j,kk^{\prime}}Mf_{\delta q\mskip
1.5mu\delta
g}-(\tilde{\boldsymbol{y}}^{j}\tau^{kk^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}+\boldsymbol{y}^{j}\tau^{kk^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}})\,M^{3}f_{\delta
q\mskip 1.5mu\delta g}^{t}$
$\displaystyle\quad-\tau^{\boldsymbol{y}j,kk^{\prime}}Mg_{\delta q\mskip
1.5mu\delta
g}-(\boldsymbol{y}^{j}\tau^{kk^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}-\tilde{\boldsymbol{y}}^{j}\tau^{kk^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}})\,M^{3}g_{\delta
q\mskip 1.5mu\delta g}^{t}\,.$ (350)

Analogous decompositions hold when the quark is replaced by an antiquark. Two-
gluon distributions are decomposed as

$\displaystyle F_{gg}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=f_{gg}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})\,,$
$\displaystyle F_{\Delta g\Delta g}$ $\displaystyle=f_{\Delta g\Delta g}\,,$
$\displaystyle F_{\Delta g\mskip 1.5mug}$ $\displaystyle=g_{\Delta g\mskip
1.5mug}\,,$ $\displaystyle F_{g\Delta g}$ $\displaystyle=g_{g\Delta g}\,,$
$\displaystyle F_{\delta g\mskip 1.5mug}^{jj^{\prime}}$
$\displaystyle=\tau^{jj^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}M^{2}f_{\delta
g\mskip
1.5mug}+\tau^{jj^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}}M^{2}g_{\delta
g\mskip 1.5mug}\,,$ $\displaystyle F_{g\mskip 1.5mu\delta g}^{jj^{\prime}}$
$\displaystyle=\tau^{jj^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}M^{2}f_{g\mskip
1.5mu\delta
g}+\tau^{jj^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}}M^{2}g_{g\mskip
1.5mu\delta g}\,,$ $\displaystyle F_{\delta g\mskip 1.5mu\delta
g}^{jj^{\prime},kk^{\prime}}$
$\displaystyle=\tau^{jj^{\prime}\\!,kk^{\prime}}f_{\delta g\mskip 1.5mu\delta
g}/2+(\tau^{jj^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}}\tau^{kk^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}}-\tau^{jj^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}\tau^{kk^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}})\,M^{4}f_{\delta
g\mskip 1.5mu\delta g}^{t}$
$\displaystyle\quad+(\tau^{jj^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}}\tau^{kk^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}+\tau^{jj^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}\tau^{kk^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}})\,M^{4}g_{\delta
g\mskip 1.5mu\delta g}^{s}$
$\displaystyle\quad+(\tau^{jj^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}}\tau^{kk^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}-\tau^{jj^{\prime}\\!,\boldsymbol{y}\boldsymbol{y}}\tau^{kk^{\prime}\\!,\boldsymbol{y}\tilde{\boldsymbol{y}}})\,M^{4}g_{\delta
g\mskip 1.5mu\delta g}^{a}\,.$ (351)

We note that in the decomposition of each DTMD, the tensors multiplying
different terms are orthogonal to each other.

The decompositions for the nonzero DPDFs in terms of real-valued scalar
functions have already been given in Diehl:2011yj ; Diehl:2013mla . They can
be directly obtained from the DTMD decompositions above by setting all pseudo-
scalar functions (i.e. all $g$’s) equal to zero (and removing the arguments
$\boldsymbol{z}_{i}$ in $F$ and $f$).

##  Appendix G Matching coefficients for (pseudo)scalar DPDs

Using the decompositions in the previous subsection and the results for the
one-loop matching coefficients given in section 7.3.4, we can perform the
necessary tensor contractions and obtain a full set of matching equations at
the level of scalar and pseudoscalar distributions. We emphasise that this is
done in the two physical space-time dimensions, since all intermediate
calculations requiring dimensional regularisation have been completed at this
point. For definiteness, let us recall that the tensor-valued coefficients in
section 7.3.4 appear in the matching relations as

$\displaystyle{}^{R}\mskip-1.5muF^{jj^{\prime}\\!,kk^{\prime}}_{\delta g\mskip
1.5mu\delta g}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\zeta)$
$\displaystyle={}^{R\,}\mskip-1.5muC^{jj^{\prime}\\!,ll^{\prime}}_{\delta
g\mskip 1.5mu\delta
g}(x_{1}^{\prime},\boldsymbol{z}_{1};x_{1}\zeta/x_{2})\underset{x_{1}}{\otimes}{}^{R\,}\mskip-1.5muC^{kk^{\prime}\\!,mm^{\prime}}_{\delta
g\mskip 1.5mu\delta g}(x_{2}^{\prime},\boldsymbol{z}_{2};x_{2}\mskip
1.5mu\zeta/x_{1})$
$\displaystyle\quad\underset{x_{2}}{\otimes}{}^{R}\mskip-1.5muF^{ll^{\prime}\\!,mm^{\prime}}_{\delta
g\mskip 1.5mu\delta g}(x_{i}^{\prime},\boldsymbol{y};\zeta)+\\{\text{quark-
gluon mixing terms}\\}$ (352)

for the distribution $F_{\delta g\mskip 1.5mu\delta g}$ and in an analogous
manner in other relations involving DPDs that carry polarisation indices. In
an expansion up to $\mathcal{O}(\alpha_{s})$, one keeps the one-loop terms
only in one of the two matching coefficients.

In the following we give the terms with matching coefficients for the parton
with momentum fraction $x_{1}$; the terms for the parton with momentum
fraction $x_{2}$ are fully analogous. We omit colour labels $R$ for brevity.
With the tensor structure of the kernels in section 7.3.4 we then obtain for
two-quark distributions

$\displaystyle f_{qq}$
$\displaystyle=C_{qq}\underset{x_{1}}{\otimes}f_{qq}+C_{qg}\underset{x_{1}}{\otimes}f_{gq}-\boldsymbol{y}^{2}M^{2}\omega_{1}C_{q\mskip
1.5mu\delta g}\underset{x_{1}}{\otimes}f_{\delta g\mskip 1.5muq}\,,$
$\displaystyle f_{\Delta q\Delta q}$ $\displaystyle=C_{\Delta q\Delta
q}\underset{x_{1}}{\otimes}f_{\Delta q\Delta q}+C_{\Delta q\Delta
g}\underset{x_{1}}{\otimes}f_{\Delta g\Delta q}\,,$ $\displaystyle g_{\Delta
q\mskip 1.5muq}$ $\displaystyle=g_{q\Delta
q}=0\,,\phantom{\underset{x_{1}}{\otimes}}$
$\displaystyle|\boldsymbol{y}|Mf_{\delta q\mskip 1.5muq}$
$\displaystyle=|\boldsymbol{y}|MC_{\delta q\mskip 1.5mu\delta
q}\underset{x_{1}}{\otimes}f_{\delta q\mskip 1.5muq}\,,$
$\displaystyle|\boldsymbol{y}|Mg_{\delta q\mskip 1.5muq}$
$\displaystyle=0\,,\phantom{\underset{x_{1}}{\otimes}}$
$\displaystyle|\boldsymbol{y}|Mf_{q\mskip 1.5mu\delta q}$
$\displaystyle=|\boldsymbol{y}|MC_{qq}\underset{x_{1}}{\otimes}f_{q\mskip
1.5mu\delta q}+|\boldsymbol{y}|MC_{qg}\underset{x_{1}}{\otimes}f_{g\mskip
1.5mu\delta q}+|\boldsymbol{y}|M\omega_{1}C_{q\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}(f_{\delta g\mskip 1.5mu\delta
q}+\boldsymbol{y}^{2}M^{2}f_{\delta g\mskip 1.5mu\delta q}^{t})\,,$
$\displaystyle|\boldsymbol{y}|Mg_{q\mskip 1.5mu\delta q}$
$\displaystyle=-|\boldsymbol{y}|M\tilde{\omega}_{1}C_{q\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}(f_{\delta g\mskip 1.5mu\delta
q}-\boldsymbol{y}^{2}M^{2}f_{\delta g\mskip 1.5mu\delta q}^{t})\,,$
$\displaystyle f_{\delta q\mskip 1.5mu\delta q}$ $\displaystyle=C_{\delta
q\mskip 1.5mu\delta q}\underset{x_{1}}{\otimes}f_{\delta q\mskip 1.5mu\delta
q}\,,$ $\displaystyle\boldsymbol{y}^{2}M^{2}f_{\delta q\mskip 1.5mu\delta
q}^{t}$ $\displaystyle=\boldsymbol{y}^{2}M^{2}C_{\delta q\mskip 1.5mu\delta
q}\underset{x_{1}}{\otimes}f_{\delta q\mskip 1.5mu\delta q}^{t}\,,$
$\displaystyle g_{\delta q\mskip 1.5mu\delta q}^{s}$ $\displaystyle=g_{\delta
q\mskip 1.5mu\delta q}^{a}=0\,,$ (353)

where we have defined

$\displaystyle\omega_{1}$
$\displaystyle=-\frac{2\tau^{\boldsymbol{y}\boldsymbol{y},\boldsymbol{z}_{1}\\!\boldsymbol{z}_{1}}}{\boldsymbol{y}^{2}\mskip
1.5mu\boldsymbol{z}_{1}^{2}}=-\cos(2\varphi_{1})\,,$
$\displaystyle\tilde{\omega}_{1}$
$\displaystyle=-\frac{2\tau^{\boldsymbol{y}\tilde{\boldsymbol{y}},\boldsymbol{z}_{1}\\!\boldsymbol{z}_{1}}}{\boldsymbol{y}^{2}\mskip
1.5mu\boldsymbol{z}_{1}^{2}}=\sin(2\varphi_{1})$ (354)

with $\tilde{\boldsymbol{y}}$ given in (348). Here $\varphi_{1}$ is the angle
between $\boldsymbol{y}$ and $\boldsymbol{z}_{1}$ in the transverse plane,
oriented such that
$\epsilon^{jj^{\prime}}_{\phantom{1}}\boldsymbol{y}^{j}_{\phantom{1}}\boldsymbol{z}^{j^{\prime}}_{1}=|\boldsymbol{y}|\mskip
1.5mu|\boldsymbol{z}_{1}|\sin\varphi_{1}$. For mixed quark and gluon
distributions the matching reads

$\displaystyle f_{gq}$
$\displaystyle=C_{gg}\underset{x_{1}}{\otimes}f_{gq}-\boldsymbol{y}^{2}M^{2}\omega_{1}C_{g\mskip
1.5mu\delta g}\underset{x_{1}}{\otimes}f_{\delta g\mskip
1.5muq}+C_{gq}\underset{x_{1}}{\otimes}f_{qq}\,,$ $\displaystyle f_{\Delta
g\Delta q}$ $\displaystyle=C_{\Delta g\Delta
g}\underset{x_{1}}{\otimes}f_{\Delta g\Delta q}+C_{\Delta g\Delta
q}\underset{x_{1}}{\otimes}f_{\Delta q\Delta q}\,,$ $\displaystyle g_{\Delta
g\mskip 1.5muq}$ $\displaystyle=g_{g\Delta
q}=0\,,\phantom{\underset{x_{1}}{\otimes}}$
$\displaystyle\boldsymbol{y}^{2}M^{2}f_{\delta g\mskip 1.5muq}$
$\displaystyle=\boldsymbol{y}^{2}M^{2}C_{\delta g\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}f_{\delta g\mskip 1.5muq}-\omega_{1}C_{\delta
g\mskip 1.5mug}\underset{x_{1}}{\otimes}f_{gq}-\omega_{1}C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{qq}\,,$
$\displaystyle\boldsymbol{y}^{2}M^{2}g_{\delta g\mskip 1.5muq}$
$\displaystyle=-\tilde{\omega}_{1}C_{\delta g\mskip
1.5mug}\underset{x_{1}}{\otimes}f_{gq}-\tilde{\omega}_{1}C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{qq}\,,$
$\displaystyle|\boldsymbol{y}|Mf_{g\mskip 1.5mu\delta q}$
$\displaystyle=|\boldsymbol{y}|MC_{gg}\underset{x_{1}}{\otimes}f_{g\mskip
1.5mu\delta q}+|\boldsymbol{y}|M\omega_{1}C_{g\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}(f_{\delta g\mskip 1.5mu\delta
q}+\boldsymbol{y}^{2}M^{2}f_{\delta g\mskip 1.5mu\delta
q}^{t})+|\boldsymbol{y}|MC_{gq}\underset{x_{1}}{\otimes}f_{q\mskip 1.5mu\delta
q}\,,$ $\displaystyle|\boldsymbol{y}|Mg_{g\mskip 1.5mu\delta q}$
$\displaystyle=-|\boldsymbol{y}|M\tilde{\omega}_{1}C_{g\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}(f_{\delta g\mskip 1.5mu\delta
q}-\boldsymbol{y}^{2}M^{2}f_{\delta g\mskip 1.5mu\delta q}^{t})\,,$
$\displaystyle 2|\boldsymbol{y}|Mf_{\delta g\mskip 1.5mu\delta q}$
$\displaystyle=2|\boldsymbol{y}|MC_{\delta g\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}f_{\delta g\mskip 1.5mu\delta
q}+|\boldsymbol{y}|M\omega_{1}C_{\delta g\mskip
1.5mug}\underset{x_{1}}{\otimes}f_{g\mskip 1.5mu\delta
q}+|\boldsymbol{y}|M\omega_{1}C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{q\mskip 1.5mu\delta q}\,,$ $\displaystyle
2|\boldsymbol{y}|^{3}M^{3}f_{\delta g\mskip 1.5mu\delta q}^{t}$
$\displaystyle=2|\boldsymbol{y}|^{3}M^{3}C_{\delta g\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}f_{\delta g\mskip 1.5mu\delta
q}^{t}+|\boldsymbol{y}|M\omega_{1}C_{\delta g\mskip
1.5mug}\underset{x_{1}}{\otimes}f_{g\mskip 1.5mu\delta
q}+|\boldsymbol{y}|M\omega_{1}C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{q\mskip 1.5mu\delta q}\,,$ $\displaystyle
2|\boldsymbol{y}|Mg_{\delta g\mskip 1.5mu\delta q}$
$\displaystyle=-|\boldsymbol{y}|M\tilde{\omega}_{1}(C_{\delta g\mskip
1.5mug}\underset{x_{1}}{\otimes}f_{g\mskip 1.5mu\delta q}+C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{q\mskip 1.5mu\delta q})\,,$ $\displaystyle
2|\boldsymbol{y}|^{3}M^{3}g_{\delta g\mskip 1.5mu\delta q}^{t}$
$\displaystyle=-|\boldsymbol{y}|M\tilde{\omega}_{1}(C_{\delta g\mskip
1.5mug}\underset{x_{1}}{\otimes}f_{g\mskip 1.5mu\delta q}+C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{q\mskip 1.5mu\delta q})$ (355)

and

$\displaystyle f_{qg}$
$\displaystyle=C_{qq}\underset{x_{1}}{\otimes}f_{qg}+C_{qg}\underset{x_{1}}{\otimes}f_{gg}-\boldsymbol{y}^{2}M^{2}\omega_{1}C_{q\mskip
1.5mu\delta g}\underset{x_{1}}{\otimes}f_{\delta g\mskip 1.5mug}\,,$
$\displaystyle f_{\Delta q\Delta g}$ $\displaystyle=C_{\Delta q\Delta
q}\underset{x_{1}}{\otimes}f_{\Delta q\Delta g}+C_{\Delta q\Delta
g}\underset{x_{1}}{\otimes}f_{\Delta g\Delta g}\,,$ $\displaystyle g_{\Delta
q\mskip 1.5mug}$ $\displaystyle=g_{q\Delta
g}=0\,,\phantom{\underset{x_{1}}{\otimes}}$
$\displaystyle|\boldsymbol{y}|Mf_{\delta q\mskip 1.5mug}$
$\displaystyle=|\boldsymbol{y}|MC_{\delta q\mskip 1.5mu\delta
q}\underset{x_{1}}{\otimes}f_{\delta q\mskip 1.5mug}\,,$
$\displaystyle|\boldsymbol{y}|Mg_{\delta q\mskip 1.5mug}$
$\displaystyle=0\,,\phantom{\underset{x_{1}}{\otimes}}$
$\displaystyle\boldsymbol{y}^{2}M^{2}f_{q\mskip 1.5mu\delta g}$
$\displaystyle=\boldsymbol{y}^{2}M^{2}C_{qq}\underset{x_{1}}{\otimes}f_{q\mskip
1.5mu\delta
g}+\boldsymbol{y}^{2}M^{2}C_{qg}\underset{x_{1}}{\otimes}f_{g\mskip
1.5mu\delta g}-\omega_{1}C_{q\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}(f_{\delta g\mskip 1.5mu\delta
g}-\boldsymbol{y}^{4}M^{4}f_{\delta g\mskip 1.5mu\delta g}^{t})\,,$
$\displaystyle\boldsymbol{y}^{2}M^{2}g_{q\mskip 1.5mu\delta g}$
$\displaystyle=-\tilde{\omega}_{1}C_{q\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}(f_{\delta g\mskip 1.5mu\delta
g}+\boldsymbol{y}^{4}M^{4}f_{\delta g\mskip 1.5mu\delta g}^{t})\,,$
$\displaystyle|\boldsymbol{y}|Mf_{\delta q\mskip 1.5mu\delta g}$
$\displaystyle=|\boldsymbol{y}|MC_{\delta q\mskip 1.5mu\delta
q}\underset{x_{1}}{\otimes}f_{\delta q\mskip 1.5mu\delta g}\,,$
$\displaystyle|\boldsymbol{y}^{3}|M^{3}f_{\delta q\mskip 1.5mu\delta g}^{t}$
$\displaystyle=|\boldsymbol{y}|^{3}M^{3}C_{\delta q\mskip 1.5mu\delta
q}\underset{x_{1}}{\otimes}f_{\delta q\mskip 1.5mu\delta g}^{t}\,,$
$\displaystyle g_{\delta q\mskip 1.5mu\delta g}$ $\displaystyle=g_{\delta
q\mskip 1.5mu\delta g}^{t}=0\,.$ (356)

For two-gluon distributions we have

$\displaystyle f_{gg}$
$\displaystyle=C_{gg}\underset{x_{1}}{\otimes}f_{gg}-\boldsymbol{y}^{2}M^{2}\omega_{1}C_{g\mskip
1.5mu\delta g}\underset{x_{1}}{\otimes}f_{\delta g\mskip
1.5mug}+C_{gq}\underset{x_{1}}{\otimes}f_{qg}\,,$ $\displaystyle f_{\Delta
g\Delta g}$ $\displaystyle=C_{\Delta g\Delta
g}\underset{x_{1}}{\otimes}f_{\Delta g\Delta g}+C_{\Delta g\Delta
q}\underset{x_{1}}{\otimes}f_{\Delta q\Delta g}\,,$ $\displaystyle g_{\Delta
g\mskip 1.5mug}$ $\displaystyle=g_{g\Delta
g}=0\,,\phantom{\underset{x_{1}}{\otimes}}$
$\displaystyle\boldsymbol{y}^{2}M^{2}f_{\delta g\mskip 1.5mug}$
$\displaystyle=\boldsymbol{y}^{2}M^{2}C_{\delta g\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}f_{\delta g\mskip 1.5mug}-\omega_{1}C_{\delta
g\mskip 1.5mug}\underset{x_{1}}{\otimes}f_{gg}-\omega_{1}C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{qg}\,,$
$\displaystyle\boldsymbol{y}^{2}M^{2}g_{\delta g\mskip 1.5mug}$
$\displaystyle=-\tilde{\omega}_{1}C_{\delta g\mskip
1.5mug}\underset{x_{1}}{\otimes}f_{gg}-\tilde{\omega}_{1}C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{qg}\,,$
$\displaystyle\boldsymbol{y}^{2}M^{2}f_{g\mskip 1.5mu\delta g}$
$\displaystyle=\boldsymbol{y}^{2}M^{2}C_{gg}\underset{x_{1}}{\otimes}f_{g\mskip
1.5mu\delta g}-\omega_{1}C_{g\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}(f_{\delta g\mskip 1.5mu\delta
g}-\boldsymbol{y}^{4}M^{4}f_{\delta g\mskip 1.5mu\delta g}^{t})\,,$
$\displaystyle\boldsymbol{y}^{2}M^{2}g_{g\mskip 1.5mu\delta g}$
$\displaystyle=-\tilde{\omega}_{1}C_{g\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}(f_{\delta g\mskip 1.5mu\delta
g}+\boldsymbol{y}^{4}M^{4}f_{\delta g\mskip 1.5mu\delta g}^{t})\,,$
$\displaystyle 2f_{\delta g\mskip 1.5mu\delta g}$ $\displaystyle=2C_{\delta
g\mskip 1.5mu\delta g}\underset{x_{1}}{\otimes}f_{\delta g\mskip 1.5mu\delta
g}-\boldsymbol{y}^{2}M^{2}\omega_{1}C_{\delta g\mskip
1.5mug}\underset{x_{1}}{\otimes}f_{g\mskip 1.5mu\delta
g}-\boldsymbol{y}^{2}M^{2}\omega_{1}C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{q\mskip 1.5mu\delta g}\,,$ $\displaystyle
2\boldsymbol{y}^{4}M^{4}f_{\delta g\mskip 1.5mu\delta g}^{t}$
$\displaystyle=2\boldsymbol{y}^{4}M^{4}C_{\delta g\mskip 1.5mu\delta
g}\underset{x_{1}}{\otimes}f_{\delta g\mskip 1.5mu\delta
g}^{t}+\boldsymbol{y}^{2}M^{2}\omega_{1}C_{\delta g\mskip
1.5mug}\underset{x_{1}}{\otimes}f_{g\mskip 1.5mu\delta
g}+\boldsymbol{y}^{2}M^{2}\omega_{1}C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{q\mskip 1.5mu\delta g}\,,$ $\displaystyle
2\boldsymbol{y}^{4}M^{4}g_{\delta g\mskip 1.5mu\delta g}^{s}$
$\displaystyle=-\boldsymbol{y}^{2}M^{2}\tilde{\omega}_{1}C_{\delta g\mskip
1.5mug}\underset{x_{1}}{\otimes}f_{g\mskip 1.5mu\delta
g}-\boldsymbol{y}^{2}M^{2}\tilde{\omega}_{1}C_{\delta g\mskip
1.5muq}\underset{x_{1}}{\otimes}f_{q\mskip 1.5mu\delta g}\,,$ $\displaystyle
g_{\delta g\mskip 1.5mu\delta g}^{a}$ $\displaystyle=g_{\delta g\mskip
1.5mu\delta g}^{s}\,.$ (357)

Further relations are obtained by replacing quark by antiquark labels.

Let us comment on the $\boldsymbol{y}$ dependent factors in the above
relations. The factors of $|\boldsymbol{y}M|^{n}$ are such that they can be
completely absorbed into the distributions (such rescaled distributions were
denoted by $h$ in Diehl:2013mla ). The factors $\omega_{1}$ and
$\tilde{\omega}_{1}$ appear in the matching between DTMDs and DPDFs carrying
different polarisations; the azimuthal dependence they provide is required by
the conservation of angular momentum along the $z$ axis. The factor
$\tilde{\omega}_{1}$ appears in the matching between pseudoscalar and scalar
functions and ensures that the r.h.s. of the matching equations is odd under
parity.

##  Appendix H Splitting kernels for (pseudo)scalar DPDs

In this appendix, we give the matching of scalar and pseudoscalar DTMDs onto
PDFs in the region of small $\boldsymbol{y}$ and $\boldsymbol{z}_{i}$. The
following results readily follow from the splitting kernels given in section
7.4 and the tensor decompositions in appendix F. The scalar matching kernels
${}^{R}T_{a_{0}\rightarrow a_{1}a_{2}}$ used in the following are obtained
from those in section 7.4 by removing all tensors $\delta$, $\epsilon$ and
$\tau$, so that one has ${}^{1}T_{g\to q\bar{q}}=T_{F}\mskip
1.5mu(u^{2}+\bar{u}^{2})$, ${}^{1}T_{g\to\Delta q\mskip
1.5mu\bar{q}}=-iT_{F}\mskip 1.5mu(u-\bar{u})$, ${}^{1}T_{g\to\delta q\mskip
1.5mu\delta\bar{q}}=-2T_{F}\mskip 1.5muu\bar{u}$, ${}^{1}T_{q\to\delta g\mskip
1.5muq}=2C_{F}\mskip 1.5mu\bar{u}/u$, etc. The only exception to this rule is
that we define ${}^{1}T_{g\to\delta g\mskip 1.5mu\delta g}=2C_{A}\mskip
1.5muu\bar{u}$. Then the coefficients ${}^{1}T_{a_{0}\to a_{1}a_{2}}$ used
here agree with the coefficients $T_{a_{0}\to a_{1}a_{2}}$ in appendix B of
Diehl:2014vaa  for all channels. The kernels for other colour representations
are obtained from (319).

The matching between (pseudo)scalar DTMDs and PDFs then reads

$\displaystyle{}^{R}\mskip-1.5muf_{\text{spl},\,q\bar{q}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\frac{\boldsymbol{y}_{+}\boldsymbol{y}_{-}}{\boldsymbol{y}_{+}^{2}\mskip
1.5mu\boldsymbol{y}_{-}^{2}}\,\frac{\alpha_{s}}{2\pi^{2}}\,{}^{R\,}\mskip-1.5muT_{g\to
q\bar{q}}\,\frac{f_{g}(x_{1}+x_{2})}{x_{1}+x_{2}}\,,$
$\displaystyle{}^{R}\mskip-1.5muf_{\text{spl},\,\Delta
q\Delta\bar{q}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\frac{\boldsymbol{y}_{+}\boldsymbol{y}_{-}}{\boldsymbol{y}_{+}^{2}\mskip
1.5mu\boldsymbol{y}_{-}^{2}}\,\frac{\alpha_{s}}{2\pi^{2}}\,{}^{R\,}\mskip-1.5muT_{g\to\Delta
q\mskip 1.5mu\Delta\bar{q}}\,\frac{f_{g}(x_{1}+x_{2})}{x_{1}+x_{2}}\,,$
$\displaystyle{}^{R}\mskip-1.5mug_{\text{spl},\,\Delta q\mskip
1.5mu\bar{q}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\frac{\tilde{\boldsymbol{y}}(\boldsymbol{z}_{1}-\boldsymbol{z}_{2})}{\boldsymbol{y}_{+}^{2}\mskip
1.5mu\boldsymbol{y}_{-}^{2}}\,\frac{\alpha_{s}}{2\pi^{2}}\,{}^{R\,}\mskip-1.5muT_{g\to\Delta
q\mskip 1.5mu\bar{q}}\,\frac{f_{g}(x_{1}+x_{2})}{x_{1}+x_{2}}\,,$
$\displaystyle{}^{R}\mskip-1.5mug_{\text{spl},\,q\mskip
1.5mu\Delta\bar{q}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\frac{\tilde{\boldsymbol{y}}(\boldsymbol{z}_{1}-\boldsymbol{z}_{2})}{\boldsymbol{y}_{+}^{2}\mskip
1.5mu\boldsymbol{y}_{-}^{2}}\,\frac{\alpha_{s}}{2\pi^{2}}\,{}^{R\,}\mskip-1.5muT_{g\to
q\mskip 1.5mu\Delta\bar{q}}\,\frac{f_{g}(x_{1}+x_{2})}{x_{1}+x_{2}}\,,$
$\displaystyle{}^{R}\mskip-1.5muf_{\text{spl},\,\delta q\mskip
1.5mu\delta\bar{q}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\frac{\boldsymbol{y}_{+}\boldsymbol{y}_{-}}{\boldsymbol{y}_{+}^{2}\mskip
1.5mu\boldsymbol{y}_{-}^{2}}\,\frac{\alpha_{s}}{2\pi^{2}}\,{}^{R\,}\mskip-1.5muT_{g\to\delta
q\mskip 1.5mu\delta\bar{q}}\,\frac{f_{g}(x_{1}+x_{2})}{x_{1}+x_{2}}$ (358)

for quark-antiquark DTMDs. For splitting into gluon-quark distributions we
have

$\displaystyle{}^{R}\mskip-1.5muf_{\text{spl},\,gq}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\frac{\boldsymbol{y}_{+}\boldsymbol{y}_{-}}{\boldsymbol{y}_{+}^{2}\mskip
1.5mu\boldsymbol{y}_{-}^{2}}\,\frac{\alpha_{s}}{2\pi^{2}}\,{}^{R\,}\mskip-1.5muT_{q\to
gq}\,\frac{f_{q}(x_{1}+x_{2})}{x_{1}+x_{2}}\,,$
$\displaystyle{}^{R}\mskip-1.5muf_{\text{spl},\,\Delta g\Delta
q}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\frac{\boldsymbol{y}_{+}\boldsymbol{y}_{-}}{\boldsymbol{y}_{+}^{2}\mskip
1.5mu\boldsymbol{y}_{-}^{2}}\,\frac{\alpha_{s}}{2\pi^{2}}\,{}^{R\,}\mskip-1.5muT_{q\to\Delta
g\Delta q}\,\frac{f_{q}(x_{1}+x_{2})}{x_{1}+x_{2}}\,,$