that not all graphs are affected by double counting. The graphs in figure 5
for instance have no overlap with DPS, whereas there is overlap between SPS
and the SPS/DPS interference in figure 5b.

Figure 4:  Graphs that contribute to several terms in the cross section (5.3).
The rectangular boxes and oval blobs indicate the relevant hadronic matrix
elements, depending on whether the $g\to q\bar{q}$ splittings are collinear or
hard.

Figure 5:  (a) A graph for the SPS/DPS interference that does not have overlap
with DPS. (b) An SPS graph that has no overlap with DPS but with the SPS/DPS
interference.

A systematic formalism for separating the different contributions and for
adding them in the cross section without double counting has been presented in
Diehl:2017kgu . In the following, we briefly recapitulate the essentials of
this scheme for TMD factorisation, referring to that work for details. DTMDs
and other hadronic matrix elements are used in the transverse position (rather
than transverse momentum) representation. This is very convenient for treating
Collins-Soper evolution, which is then multiplicative (rather than involving
convolution integrals). The distinction between collinear and hard parton
splitting then corresponds to the transverse distance between the partons
being large or small, respectively.

Figure 6:  Transverse positions of the parton fields in parton splitting
contributions to (a) the quark-antiquark DTMD and (b) the twist-three TMD
$D_{q\bar{q}|g}$. Due to translation invariance, the latter depends only on
$\boldsymbol{y}_{+}$ and $\boldsymbol{Z}$. The positions in (b) correspond to
equation (3.9) in Diehl:2017kgu .

As can be seen in figure 6a, the relative distance between the two partons
with momentum fractions $x_{1}$ and $x_{2}$ in a DTMD is $\boldsymbol{y}_{+}$
for the partons in the amplitude and $\boldsymbol{y}_{-}$ for those in the
conjugate amplitude, with the distances defined in (120). The contribution
from small distances to the DPS cross section $\sigma_{\text{DPS}}$ is removed
by the factors $\Phi(\nu\boldsymbol{y}_{+})$ and $\Phi(\nu\boldsymbol{y}_{-})$
in (5), which we discuss further at the end of this section. Corresponding
factors are inserted in the interference
$\sigma_{\text{DPS/SPS}}+\sigma_{\text{SPS/DPS}}$ between SPS and DPS. The SPS
cross section $\sigma_{\text{SPS}}$ is not modified in this scheme and
computed as usual from TMDs and hard scattering cross sections, as described
for instance in Collins:2011zzd . The master formula for the overall cross
section then reads

$\displaystyle\sigma$
$\displaystyle=\sigma_{\text{DPS}\rule{0.0pt}{3.91806pt}}+\left[\mskip
1.5mu\sigma_{\text{DPS/SPS}}-\sigma_{\text{DPS},\,y_{-}\to
0\rule{0.0pt}{3.91806pt}}+\sigma_{\text{SPS/DPS}}-\sigma_{\text{DPS},\,y_{+}\to
0\rule{0.0pt}{3.91806pt}}\mskip 1.5mu\right]$ $\displaystyle\quad+\left[\mskip
1.5mu\sigma_{\text{SPS}\rule{0.0pt}{3.91806pt}}-\sigma_{\text{DPS/SPS},\,y_{+}\rightarrow
0}-\sigma_{\text{SPS/DPS},\,y_{-}\rightarrow
0}+\sigma_{\text{DPS},\,y_{\pm}\rightarrow 0\rule{0.0pt}{3.91806pt}}\mskip
1.5mu\right]\,,$ (147)

where all terms should be taken differential in
$dx_{1}\,dx_{2}\,d\bar{x}_{1}\,d\bar{x}_{2}\,d^{2}\boldsymbol{q}_{1}\,d^{2}\boldsymbol{q}_{2}$.
In this formula, the double counting problem between different contributions
is solved by subtraction terms, which we will discuss shortly.

Let us start with the interference between DPS and SPS. As the colour and
polarisation of the partons are not the main focus of the present section, we
will not display these quantum numbers explicitly. The following equations
hold for specific colour and helicity of each parton on the l.h.s. and with
appropriate sums on the r.h.s. The indices $\alpha_{i},\beta_{i}$ will just
denote the parton species. Up to kinematic and numeric factors, the
interference cross section is given by

$\displaystyle\frac{d\sigma_{\text{DPS/SPS}}}{dx_{1}\,dx_{2}\,d\bar{x}_{1}\,d\bar{x}_{2}\,d^{2}\boldsymbol{q}_{1}\,d^{2}\boldsymbol{q}_{2}}\,\propto\sum_{\genfrac{}{}{0.0pt}{1}{\alpha_{1},\alpha_{2},\alpha_{0}}{\beta_{1},\beta_{2},\beta_{0}}}\\!\\!\\!H_{\alpha_{1}\beta_{1}}\,H_{\alpha_{2}\beta_{2}}\,H^{*}_{\alpha_{0}\beta_{0}}\int
d^{2}\boldsymbol{Z}\,d^{2}\boldsymbol{y}_{+}\;e^{-i\mskip
1.5mu(\boldsymbol{q}_{1}+\boldsymbol{q}_{2})\boldsymbol{Z}}$
$\displaystyle\qquad\qquad\times e^{-i\mskip
1.5mu(\boldsymbol{q}_{1}-\boldsymbol{q}_{2})\boldsymbol{y}_{+}/2}\;\Phi(\nu\boldsymbol{y}_{+})\,D_{\beta_{1}\beta_{2}|\beta_{0}}(\bar{x}_{i},\boldsymbol{y}_{+},\boldsymbol{Z})\,D_{\alpha_{1}\alpha_{2}|\alpha_{0}}(x_{i},\boldsymbol{y}_{+},\boldsymbol{Z})$
(148)

for DPS in the amplitude, and by

$\displaystyle\frac{d\sigma_{\text{SPS/DPS}}}{dx_{1}\,dx_{2}\,d\bar{x}_{1}\,d\bar{x}_{2}\,d^{2}\boldsymbol{q}_{1}\,d^{2}\boldsymbol{q}_{2}}\,\propto\sum_{\genfrac{}{}{0.0pt}{1}{\alpha_{0},\alpha_{1},\alpha_{2}}{\beta_{0},\beta_{1},\beta_{2}}}\\!\\!\\!H_{\alpha_{0}\beta_{0}}\,H^{*}_{\alpha_{1}\beta_{1}}\,H^{*}_{\alpha_{2}\beta_{2}}\int
d^{2}\boldsymbol{Z}\,d^{2}\boldsymbol{y}_{-}\;e^{-i\mskip
1.5mu(\boldsymbol{q}_{1}+\boldsymbol{q}_{2})\boldsymbol{Z}}$
$\displaystyle\qquad\qquad\times e^{i\mskip
1.5mu(\boldsymbol{q}_{1}-\boldsymbol{q}_{2})\boldsymbol{y}_{-}/2}\;\Phi(\nu\boldsymbol{y}_{-})\,D_{\beta_{0}|\beta_{1}\beta_{2}}(\bar{x}_{i},\boldsymbol{y}_{-},\boldsymbol{Z})\,D_{\alpha_{0}|\alpha_{1}\alpha_{2}}(x_{i},\boldsymbol{y}_{-},\boldsymbol{Z})$
(149)

for DPS in the complex conjugate amplitude. Here
$D_{\alpha_{1}\alpha_{2}|\alpha_{0}}$ is a twist-three TMD with two partons
$\alpha_{1}\alpha_{2}$ in the amplitude and one parton $\alpha_{0}$ in the
conjugate amplitude. Its definition is similar to the one of DTMDs, with three
instead of four parton operators in the hadronic matrix element. The
transverse positions of the partons are as indicated in figure 6b. In
$D_{\alpha_{0}|\alpha_{1}\alpha_{2}}$, the roles of amplitude and conjugate
amplitude are interchanged. $H_{\alpha_{1}\beta_{1}}$ and
$H_{\alpha_{2}\beta_{2}}$ are the hard scattering amplitudes for the DPS
processes, whilst $H_{\alpha_{0}\alpha_{0}}$ is the amplitude for the hard SPS
interaction.

$\sigma_{\text{DPS/SPS}}$ receives a contribution from the graph in figure 5b
without the rectangular boxes. The quark loop in the conjugate amplitude
includes an integration over the full phase space, including the region where
the $g\to q\bar{q}$ splittings become collinear. By a Fourier transform, one
can rewrite the transverse part of the loop integration as an integral over
$\boldsymbol{y}_{-}$, after which the region of collinear splitting
corresponds to large $\boldsymbol{y}_{-}$. This contribution is already
included in $\sigma_{\text{DPS}}$, so that a subtraction term is necessary to
prevent double counting. It reads

$\displaystyle\frac{d\sigma_{\text{DPS},y_{-}\to
0}}{dx_{1}\,dx_{2}\,d\bar{x}_{1}\,d\bar{x}_{2}\,d^{2}\boldsymbol{q}_{1}\,d^{2}\boldsymbol{q}_{2}}\,\propto\,\sum_{\genfrac{}{}{0.0pt}{1}{\alpha_{1},\alpha_{2}}{\beta_{1},\beta_{2}}}H_{\alpha_{1}\beta_{1}}\,H_{\alpha_{2}\beta_{2}}\,H^{*}_{\alpha_{1}\beta_{1}}\,H^{*}_{\alpha_{2}\beta_{2}}\,$
$\displaystyle\qquad\qquad\times\int
d^{2}\boldsymbol{Z}\,d^{2}\boldsymbol{y}_{+}\,d^{2}\boldsymbol{y}_{-}\,\,e^{-i\mskip
1.5mu(\boldsymbol{q}_{1}+\boldsymbol{q}_{2})\boldsymbol{Z}-i\mskip
1.5mu(\boldsymbol{q}_{1}-\boldsymbol{q}_{2})(\boldsymbol{y}_{+}-\boldsymbol{y}_{-})/2}\;\Phi(\nu\boldsymbol{y}_{+})\,\Phi(\nu\boldsymbol{y}_{-})\,$
$\displaystyle\qquad\qquad\quad\times F_{\beta_{1}\beta_{2},\mskip
1.5muy_{-}\to
0}(\bar{x}_{i},\boldsymbol{y}_{\pm},\boldsymbol{Z})\,F_{\alpha_{1}\alpha_{2},\mskip
1.5muy_{-}\to 0}(x_{i},\boldsymbol{y}_{\pm},\boldsymbol{Z})\,.$ (150)

Compared with the DPS cross section given by (5) and (5), we have traded the
arguments $\boldsymbol{z}_{i}$, $\boldsymbol{y}$ of the DTMDs for
$\boldsymbol{y}_{\pm}$, $\boldsymbol{Z}$ and correspondingly changed the
integration according to (122). For simplicity we write $y_{-}\to 0$ instead
of $\boldsymbol{y}_{-}\to\boldsymbol{0}$ in subscripts. The DTMDs in (5.3)
correspond to the rectangular boxes in figure 5b and are given by

$\displaystyle F_{\alpha_{1}\alpha_{2},\,y_{-}\to
0}(x_{i},\boldsymbol{y}_{\pm},\boldsymbol{Z})$
$\displaystyle=\frac{\boldsymbol{y}_{-}^{l^{\prime}}}{\boldsymbol{y}_{-}^{2}}\,\bigl{[}U^{l^{\prime}}_{\alpha_{0}\to\alpha_{1}\alpha_{2}}(x_{i})\bigr{]}^{*}\,D_{\alpha_{1}\alpha_{2}|\alpha_{0}}(x_{i},\boldsymbol{y}_{+},\boldsymbol{Z})+\mathcal{O}(\alpha_{s}^{3/2})\,.$
(151)

This notation reflects that in the limit $\boldsymbol{y}_{-}\to\boldsymbol{0}$
the DTMD is dominated by the contribution where the two partons $\alpha_{1}$
and $\alpha_{2}$ in the conjugate amplitude originate from the splitting of a
single parton $\alpha_{0}$. The lowest order perturbative splitting kernel $U$
includes a factor of $\sqrt{\alpha_{s}\rule{0.0pt}{6.45831pt}}$. For given
$\alpha_{1}$ and $\alpha_{2}$, the parton species $\alpha_{0}$ is fixed at
this order, but its helicity must be summed over. As shown in Diehl:2017kgu ,
the combination
$\sigma_{\text{DPS}\rule{0.0pt}{3.91806pt}}+\sigma_{\text{DPS/SPS}}-\sigma_{\text{DPS},y_{-}\to
0\rule{0.0pt}{3.91806pt}}$ correctly treats the contribution of graph 5b to
the cross section, for small and for large $\boldsymbol{y}_{-}$, and without
double counting. We will see in more detail how this happens in section 6.5.
The other subtraction term in the first line of (5.3) is given by expressions
analogous to (5.3) and (151), with $\boldsymbol{y}_{-}$ replaced by
$\boldsymbol{y}_{+}$ and the roles of amplitude and conjugate amplitude being
interchanged.

The subtraction terms $\sigma_{\text{DPS/SPS},y_{+}\rightarrow 0}$ and
$\sigma_{\text{SPS/DPS},y_{-}\rightarrow 0}$ in the second line of (5.3) are
obtained from (5.3) and (5.3) by replacing all distributions $D$ with their
perturbative splitting approximation. For the case shown in figure 6b, this is
given by

$\displaystyle D_{\alpha_{1}\alpha_{2}|\alpha_{0},\,y_{+}\to
0}(x_{i},\boldsymbol{y}_{+},\boldsymbol{Z})$
$\displaystyle=\frac{\boldsymbol{y}_{+}^{l}}{\boldsymbol{y}_{+}^{2}}\;U_{\alpha_{0}\rightarrow\alpha_{1}\alpha_{2}}^{l}(x_{i})\,f_{\alpha_{0}}(x_{1}+x_{2},\boldsymbol{Z})+\mathcal{O}(\alpha_{s}^{3/2})\,,$
(152)

where on the r.h.s. we have a single parton TMD and the same splitting kernel
$U$ as in (151). An analogous expression holds for the limit
$\boldsymbol{y}_{-}\to\boldsymbol{0}$ of
$D_{\alpha_{0}|\alpha_{1}\alpha_{2}}$. The combination
$\sigma_{\text{DPS/SPS}}+\sigma_{\text{SPS}\rule{0.0pt}{3.91806pt}}-\sigma_{\text{DPS/SPS},y_{+}\rightarrow
0}$ correctly represents graph 5c in the region where the quark loop in the
conjugate amplitude is hard.

Finally, the term $\sigma_{\text{DPS},y_{\pm}\rightarrow 0}$ in (5.3) ensures
the correct treatment of the case when both quark loops in graph 5a are in the
collinear region. This term is given by the DPS cross section with each DTMD
approximated for the regime where both parton pairs are produced by
perturbative splitting:

$\displaystyle F_{\alpha_{1}\alpha_{2},\,y_{\pm}\to
0}(x_{i},\boldsymbol{y}_{\pm},\boldsymbol{Z})$
$\displaystyle=\frac{\boldsymbol{y}_{+}^{l}}{\boldsymbol{y}_{+}^{2}}\,\frac{\boldsymbol{y}_{-}^{l^{\prime}}}{\boldsymbol{y}_{-}^{2}}\,U^{l}_{\alpha_{0}\to\alpha_{1}\alpha_{2}}(x_{i})\,\bigl{[}U^{l^{\prime}}_{\alpha_{0}\to\alpha_{1}\alpha_{2}}(x_{i})\bigr{]}^{*}\,f_{\alpha_{0}}(x_{1}+x_{2},\boldsymbol{Z})+\mathcal{O}(\alpha_{s}^{2})\,.$
(153)

This is illustrated in figure 6a and can be obtained by inserting the
expression (152) into (151). Notice that
$\sigma_{\text{DPS},y_{\pm}\rightarrow 0}$ appears with a plus rather than a
minus sign in (5.3). This is a consequence of the recursive nature of double
counting subtractions, as explained in section 4.2 of Diehl:2017kgu .

To make contact with our notation in section 6.3.1, we rewrite (153) in terms
of a single splitting kernel as888The convention for $T^{jj^{\prime}}$ here is
the same as in Diehl:2017kgu  (arXiv version 2) and differs from the one in
Diehl:2011yj  as specified in section 7.4.

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2},\,y_{\pm}\to
0}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\frac{\boldsymbol{y}_{+}^{l}\mskip
1.5mu\boldsymbol{y}_{-}^{l^{\prime}}}{\boldsymbol{y}_{+}^{2}\mskip
1.5mu\boldsymbol{y}_{-}^{2}}\,\sum_{a_{0}}\frac{\alpha_{s}}{2\pi^{2}}\;{}^{R}T_{a_{0}\to
a_{1}a_{2}}^{ll^{\prime}}\biggl{(}\frac{x_{1}}{x_{1}+x_{2}}\biggr{)}\,\frac{f_{a_{0}}(x_{1}+x_{2},\boldsymbol{Z})}{x_{1}+x_{2}}+\mathcal{O}(\alpha_{s}^{2})\,,$
(154)

where we have restored the dependence on the colour representation and the
notation with labels $a_{i}$ that specify parton species and polarisation. For
transverse quark or linear gluon polarisation, the parton distributions and
the kernel $T_{a_{0}\to a_{1}a_{2}}$ carry additional transverse indices,
which are not displayed here (see section 7.4). The sum over $a_{0}$ reflects
the fact that there are TMDs for these polarisation states, even if the proton
is unpolarised.

We can now discuss the role of the function $\Phi$ introduced earlier. The
short-distance behaviour of the DTMDs and twist-three distributions, given in
(151) and (152), results in an integral

$\displaystyle\int\limits_{|\boldsymbol{y}_{+}|\ll
1/\Lambda}\frac{d^{2}\boldsymbol{y}_{+}}{\boldsymbol{y}_{+}^{2}}\;e^{-i\mskip
1.5mu(\boldsymbol{q}_{1}-\boldsymbol{q}_{2})\mskip
1.5mu\boldsymbol{y}_{+}/2}\;\Phi(\nu\boldsymbol{y}_{+})$ (155)

and its analogue for $\boldsymbol{y}_{-}$ in the cross section formulae for
DPS and for the SPS/DPS interference. The function $\Phi$ must satisfy
$\Phi(\boldsymbol{u})\to 0$ for $\boldsymbol{u}\to\boldsymbol{0}$ to ensure
that these integrals converge at small distances (rather than having a
logarithmic divergence). In order to keep the regions of large
$\boldsymbol{y}_{+}$ and $\boldsymbol{y}_{-}$ unaffected, one should
furthermore have $\Phi(\boldsymbol{u})\to 1$ for $|\boldsymbol{u}|\gg 1$. A
simple choice for $\Phi(\boldsymbol{u})$ is a step function in
$|\boldsymbol{u}|$, which corresponds to a hard cutoff on the
$\boldsymbol{y}_{\pm}$ integrals.

The integral (155) goes like $\log(\nu/q_{T})$ and thus depends on the
artificial parameter $\nu$ that controls which distances
$\boldsymbol{y}_{\pm}$ are included in what we _define_ to be DPS rather than
SPS. The $\nu$ dependence of the different terms in the cross section (5.3)
cancels to the perturbative order of the calculation, in close analogy with
the familiar case of renormalisation and factorisation scale dependence. As
discussed in Diehl:2017kgu , a suitable choice for $\nu$ is the lowest hard
scale, $\text{min}(Q_{1},Q_{2})$. With this choice, $\sigma_{\text{DPS}}$
contains a squared logarithm $\log^{2}(Q/q_{T})$, where we recall that $Q$
denotes the generic size of $Q_{1}$ and $Q_{2}$. After the subtractions in the
first line of (5.3), the SPS/DPS interference contains only a single
$\log(Q/q_{T})$ and the subtracted SPS term in the second line has no such
logarithm at the corresponding order in $\alpha_{s}$. This has an important
consequence. If one is satisfied with leading logarithmic accuracy, then the
terms in square brackets in (5.3) can be neglected at the order in
$\alpha_{s}$ where they have overlap with DPS, because at that order it is
$\sigma_{\text{DPS}}$ that has the highest power of $\log(Q/q_{T})$.

We should note that the structure of the higher-order terms in (151), (152)
and (153) is currently unknown. In particular, it cannot be excluded that
singularities will appear at points other than
$\boldsymbol{y}_{+}=\boldsymbol{0}$ and $\boldsymbol{y}_{-}=\boldsymbol{0}$.
The treatment of ultraviolet divergences and of the associated double counting
issues may thus be more involved beyond the leading order discussed here.

To conclude this section, we briefly discuss the scale evolution of the twist-
three TMDs $D$. As already noted, they are defined in terms of hadronic matrix
elements, with three parton fields and the corresponding Wilson lines instead
of the four fields we have in the unsubtracted DTMDs in (3.1). The three
parton fields are at different transverse positions (if all three positions
are equal, one has collinear twist-three distributions). It is easy to see
that the analysis of factorisation for the SPS/DPS interference will yield a
soft factor, with three products of Wilson lines associated with the three
parton fields in the hadronic matrix elements (rather than four products of
Wilson lines as in figure 1). One can readily adapt the discussion in sections
3.2 to 3.4 to this case and absorb the soft factor into the twist-three
distributions. The colour algebra is much simpler than in the DPS case, with a
single colour representation in the $q\bar{q}g$ channel, and two for the
channel with three gluons (one constructed with $f^{abc}$ and the other with
$d^{abc}$). The rapidity dependence of $D_{\alpha_{1}\alpha_{2}|\alpha_{0}}$
is finally given by a Collins-Soper equation as in (128) with a kernel
$K_{\alpha_{1}\alpha_{2}|\alpha_{0}}$. If one takes different renormalisation
scales $\mu_{1},\mu_{2}$ and $\mu_{0}$ for the three parton legs, then the
dependence of $D$ and $K$ on these scales is given by the analogues of (124)
and (129), with $\gamma_{K,a}$ replaced by $\gamma_{K,a}/2$ and $\gamma_{F,a}$
by $\gamma_{E,a}$ or $(\gamma_{E,a})^{*}$ as discussed at the end of section
3.4. The rapidity parameter $\zeta$ is defined by (125) and must be properly
rescaled in the argument of $\gamma_{F,a}$, namely by $x_{1}\zeta/x_{2}$,
$x_{2}\mskip 1.5mu\zeta/x_{1}$ and $(x_{1}+x_{2})^{2}\mskip
1.5mu\zeta/(x_{1}x_{2})$ for the partons with momentum fractions $x_{1}$,
$x_{2}$ and $x_{1}+x_{2}$, respectively.

##  6 Matching for small but perturbative transverse momenta

In the multi-scale regime $\Lambda\ll q_{T}\ll Q$, where the transverse
momenta $|\boldsymbol{q}_{i}|$ are small compared with $Q_{i}$ but large
compared with the scale $\Lambda$ of nonperturbative physics, the DPS cross
section, as well as its combination with SPS can be considerably simplified.
The reason is that the additional scale $q_{T}$ allows for more perturbative
calculations, which greatly enhances the predictive power of the theory.

We start our analysis by discussing the different transverse distance regions
in the DPS factorisation formula and find that two distinct regions of
$\boldsymbol{y}$ are relevant. The corresponding short-distance expansions of
DTMDs are derived in sections 6.2 and 6.3, respectively. In section 6.4 we
discuss the master formula (206) for combining the two types of expansion in
the cross section, where care must be taken to avoid double counting. After
this, we revisit the combination between DPS and SPS just discussed in section
5.3. In section 6.6 we go through the perturbative ingredients that are
necessary for a given resummation order as specified in table 3; their
availability determines the perturbative accuracy of cross section
computations.

###  6.1 Regions of transverse momenta and distances

In the DPS cross section (5), contributions from distances
$|\boldsymbol{z}_{i}|\gg 1/|\boldsymbol{q}_{i}|$ are suppressed by
oscillations of the Fourier exponent
$e^{-i\boldsymbol{q}_{1}\boldsymbol{z}_{1}-i\boldsymbol{q}_{2}\boldsymbol{z}_{2}}$.
If $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$ are sufficiently large, this
keeps $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$ in the region where
perturbation theory has predictive power for the dependence on these
variables.

Care is required if $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$ are
perturbatively large but $|\boldsymbol{q}_{1}+\boldsymbol{q}_{2}|$ is of order
$\Lambda$. According to (122) the sum
$|\boldsymbol{z}_{1}+\boldsymbol{z}_{2}|$ can then reach large values of order
$1/\Lambda$, as long as $|\boldsymbol{z}_{1}-\boldsymbol{z}_{2}|$ remains
small. The individual oscillations of
$e^{-i\boldsymbol{q}_{1}\boldsymbol{z}_{1}}$ and
$e^{-i\boldsymbol{q}_{2}\boldsymbol{z}_{2}}$ cancel in that case. Likewise,
$|\boldsymbol{z}_{1}-\boldsymbol{z}_{2}|$ can become as large as $1/\Lambda$
if $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$ are perturbatively large but
$|\boldsymbol{q}_{1}-\boldsymbol{q}_{2}|$ is of order $\Lambda$. The
perturbative splitting form (154) of DTMDs has a factorised dependence on the
variables $\boldsymbol{z}_{1}+\boldsymbol{z}_{2}$ and
$\boldsymbol{z}_{1}-\boldsymbol{z}_{2}$, and the short-distance expansion
(193) we will derive from it requires both distances to be small. We therefore
do not consider the particular phase space regions just mentioned and require
$|\boldsymbol{q}_{1}|$ and $|\boldsymbol{q}_{2}|$ as well as
$|\boldsymbol{q}_{1}+\boldsymbol{q}_{2}|$ and
$|\boldsymbol{q}_{1}-\boldsymbol{q}_{2}|$ to be much larger than $\Lambda$.
The dominant contribution to the $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$
integrals in the cross section formula (5) is then in the perturbative region.
Even with this requirement, one could still consider a region where
$|\boldsymbol{q}_{1}+\boldsymbol{q}_{2}|$ is much smaller than
$|\boldsymbol{q}_{1}|$ and $|\boldsymbol{q}_{2}|$. This multi-scale regime has
been discussed in section 5.2.1 of Diehl:2011yj  and in Blok:2011bu ;
Blok:2013bpa . We shall not investigate it in the present work.

Other observables in which sensitivity to large $\boldsymbol{z}_{1}$ and
$\boldsymbol{z}_{2}$ can be avoided are obtained by integrating the cross
section over one or both of $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$ up
to $q_{\text{max},1}$ and $q_{\text{max},2}$ respectively. In

$\displaystyle\int
d^{2}\boldsymbol{q}\;\theta\bigl{(}q_{\text{max}}^{2}-\boldsymbol{q}^{2}\bigr{)}\int\frac{d^{2}\boldsymbol{z}}{(2\pi)^{2}}\;e^{-i\boldsymbol{q}\boldsymbol{z}}\,W(\boldsymbol{z},\ldots)$
$\displaystyle=q_{\text{max}}\int\frac{d^{2}\boldsymbol{z}}{2\pi|\boldsymbol{z}|}\,J_{1}\bigl{(}q_{\text{max}}|\boldsymbol{z}|\bigr{)}\,W(\boldsymbol{z},\ldots)$
(156)

the region of $|\boldsymbol{z}|\gg 1/q_{\text{max}}$ is damped by
oscillations, as it is by the exponential factors in the differential cross
section. To avoid that the product of the two Bessel functions has a non-
oscillating component for $|\boldsymbol{z}_{1}|\approx|\boldsymbol{z}_{2}|$,
one must take $q_{\text{max},1}\neq q_{\text{max},2}$.

We henceforth assume that both $|\boldsymbol{z}_{1}|$ and
$|\boldsymbol{z}_{2}|$ are small, namely of order $1/q_{T}$. Since
$\boldsymbol{y}$ is integrated over without any Fourier exponent in (5) we
must still consider the two cases $|\boldsymbol{y}|\sim 1/\Lambda$ and
$|\boldsymbol{y}|\sim 1/q_{T}$, which are referred to as large
$\boldsymbol{y}$ and small $\boldsymbol{y}$ in the following. For even smaller
$|\boldsymbol{y}|\ll 1/q_{T}$, the approximations for DPS are not valid and
the subtraction formalism sketched in section 5.3 comes into play. The short-
distance expansion of DPDs in the two DPS regions just mentioned is quite
different and will now be discussed in turn.

###  6.2 The large-$\boldsymbol{y}$ region

For small $\boldsymbol{z}_{1}$, $\boldsymbol{z}_{2}$ and large
$\boldsymbol{y}$, the DTMD $F(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$ can be
computed in terms of DPDFs $F(x_{i},\boldsymbol{y})$, in close analogy to the
case of single parton distributions. This is commonly referred to as
“matching”. For DTMDs there are complications from colour, but our results
from section 4.4 lead to considerable simplifications. Generalising the
relations (9) and (2) between TMDs and PDFs, we can derive the general
matching equation (6.2.1) and the explicit solution (6.2.2) of the evolution
equations for DTMDs. This solution is promoted to the cross section level in
(6.2.2). We emphasise that the simple structure of (6.2.2) and (6.2.2) is due
to the additive form (171) of the Collins-Soper matrix kernel in the short-
distance limit considered here.

####  6.2.1 Short-distance expansion

We now derive the approximation of DTMDs in the region
$|\boldsymbol{z}_{1}|,|\boldsymbol{z}_{2}|\ll 1/\Lambda$, with
$\boldsymbol{y}$ kept fixed at a value such that
$|\boldsymbol{z}_{1}|,|\boldsymbol{z}_{2}|\ll|\boldsymbol{y}|$. Corrections to
this approximation are suppressed by powers of $|\boldsymbol{z}_{i}|\Lambda$
and of $|\boldsymbol{z}_{i}|/|\boldsymbol{y}|$. To make the analogy between
DTMDs and single parton TMDs transparent, we formulate the matching in terms
of operator product expansions around $\boldsymbol{z}_{i}=\boldsymbol{0}$.

Let us start with the expansion of the soft factor. We need the
small-$\boldsymbol{z}$ limit of the operators
$O_{S,a}(\boldsymbol{y},\boldsymbol{z})$ defined in (3.2), which contain two
Wilson line pairs separated by a transverse distance $\boldsymbol{z}$. The
colour structure of the short-distance expansion can be deduced from the
corresponding Feynman graphs. To minimise the number of eikonal propagators
that carry large momentum, interactions at the hard scale $1/|\boldsymbol{z}|$
must take place closest to the points where the left- and right-moving Wilson
lines meet, as shown in figure 7. In terms of the operators in (3.2) and
(3.2), we thus find the structure

$\displaystyle\bigl{[}{O}_{S,a}(\boldsymbol{y},\boldsymbol{z})\bigr{]}^{\underline{r},\underline{u}}$
$\displaystyle=C_{S,a}^{\underline{s}\,\underline{t}}(\boldsymbol{z})\;\bigl{[}O_{S,a}(\boldsymbol{y})\bigr{]}^{\underline{r},\underline{u}}_{\underline{s},\underline{t}}$
(157)

with a short-distance coefficient $C_{S,a}$ that can be computed in
perturbation theory. Here and in the following, the absence of the argument
$\boldsymbol{z}$ in operators or functions means that one has set
$\boldsymbol{z}=\boldsymbol{0}$. The analogue of this convention will be used
for arguments $\boldsymbol{z}_{i}$. Introducing colour projected coefficients
${}^{R}C_{S,a}=P_{R}^{\underline{s}\,\underline{t}}\,C_{S,a}^{\underline{s}\,\underline{t}}\big{/}m(R)$
and using (72), we obtain

$\displaystyle\bigl{[}O_{S,a}(\boldsymbol{y},\boldsymbol{z})\bigr{]}^{\underline{r},\underline{u}}$
$\displaystyle=\sum_{R}{}^{R}C_{S,a}(\boldsymbol{z})\,P_{R}^{\underline{s}\,\underline{t}}\,\bigl{[}O_{S,a}(\boldsymbol{y})\bigr{]}^{\underline{r},\underline{u}}_{\underline{s},\underline{t}}$
$\displaystyle=\sum_{R}{}^{R}C_{S,a}(\boldsymbol{z})\,P_{R}^{\underline{r}\,\underline{s}}\,P_{R}^{\underline{u}\,\underline{t}}\,\bigl{[}O_{S,a}(\boldsymbol{y})\bigr{]}^{\underline{s},\underline{t}}_{\underline{v},\underline{v}}\,,$
(158)

where the second step is completely analogous to the derivation of the
relation (4.4) between the extended soft factor and the usual one. Taking the
vacuum expectation value of
$O_{S,a_{1}}(\boldsymbol{y},\boldsymbol{z}_{1})\,O_{S,a_{2}}(\boldsymbol{0},\boldsymbol{z}_{2})$
with appropriate colour projections, we obtain

$\displaystyle{}^{RR^{\prime}}\\!S_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle={}^{R}C_{S,a_{1}}(\boldsymbol{z}_{1})\,{}^{R}C_{S,a_{2}}(\boldsymbol{z}_{2})\,{}^{RR}\mskip-1.5muS(\boldsymbol{y})\,\delta_{RR^{\prime}}\,.$
(159)

For computing the short-distance coefficients, it will be convenient to
contract (6.2.1) with $P_{R^{\prime}}^{\underline{r}\,\underline{u}}$, which
gives

$\displaystyle
P_{R}^{\underline{r}\,\underline{u}}\bigl{[}O_{S,a}(\boldsymbol{y},\boldsymbol{z})\bigr{]}^{\underline{r},\underline{u}}$
$\displaystyle={}^{R}C_{S,a}(\boldsymbol{z})\,P_{R}^{\underline{r}\,\underline{u}}\,\bigl{[}O_{S,a}(\boldsymbol{y})\bigr{]}^{\underline{r},\underline{u}}\,.$
(160)

Figure 7:  Interactions of Wilson lines in the operator
$O_{S,a}(\boldsymbol{y},\boldsymbol{z})$ in the limit of small
$\boldsymbol{z}$. Short-distance interactions are located in the box at the
centre, whereas gluons outside the box indicate possible interactions with
other Wilson lines at distances much larger than $\boldsymbol{z}$. In the soft
factor $S_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})$ these distances are
of order $\boldsymbol{y}$ for small $\boldsymbol{z}_{1}$ and
$\boldsymbol{z}_{2}$.

For approximating the unsubtracted collinear matrix element $F_{\text{us}}$,
we need the small-$\boldsymbol{z}$ limit of the parton operators in (21). We
use their projections (79) on definite colour representations, which cannot
mix with each other due to gauge invariance. The short-distance expansion
involves a convolution integral over momentum fractions and reads999The lower
limit of the convolution in (161) is $x^{\prime}=0$. When matrix elements
between hadrons with plus momentum $p^{+}$ are taken, the lower limit becomes
$x^{\prime}=x$, as in (10).

$\displaystyle{}^{R}O^{\mskip
1.5mu\underline{r}}_{a}(x,\boldsymbol{y},\boldsymbol{z})$
$\displaystyle=\sum_{b}{}^{R}C_{\text{us},ab}(x^{\prime},\boldsymbol{z})\,\underset{x}{\otimes}{}^{R}O^{\mskip
1.5mu\underline{r}}_{b}(x^{\prime},\boldsymbol{y})\,.$ (161)

If $a$ or $b$ refer to transverse quark or linear gluon polarisation, then the
associated operators carry additional transverse indices, and the same holds
of course for $C_{\text{us},ab}$. Taking matrix elements as in (80), we obtain
an approximation

$\displaystyle{}^{R}\mskip-1.5muF_{\text{us},a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\sum_{b_{1},b_{2}}{}^{R}C_{\text{us},a_{1}b_{1}}(x_{1}^{\prime},\boldsymbol{z}_{1})\underset{x_{1}}{\otimes}{}^{R}C_{\text{us},a_{2}b_{2}}(x_{2}^{\prime},\boldsymbol{z}_{2})\underset{x_{2}}{\otimes}{}^{R}\mskip-1.5muF_{\text{us},b_{1}b_{2}}(x_{i}^{\prime},\boldsymbol{y})$
(162)

of unsubtracted DTMDs in terms of unsubtracted DPDFs and perturbative
coefficients $C_{\text{us},ab}(x,\boldsymbol{z})$.

We now combine the soft factor with the unsubtracted DTMD. Using the relation
(137) between $s$ and $S$, we obtain

$\displaystyle{}^{RR^{\prime}}\mskip-1.5mus_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},Y)$
$\displaystyle=\sqrt{{}^{R}C_{S,a_{1}}(\boldsymbol{z}_{1};\mu_{1},2Y)}\,\sqrt{{}^{R}C_{S,a_{2}}(\boldsymbol{z}_{2};\mu_{2},2Y)}\;{}^{RR}\mskip-1.5mus(\boldsymbol{y};\mu_{i},Y)\,\delta_{RR^{\prime}}$
(163)

from (159), where we have restored the dependence on renormalisation scales
and rapidities. Combining this with (162) and using the definition (84) of
DPDs, we get

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},Y_{C})$
$\displaystyle=\sum_{b_{1},b_{2}}{}^{R}C_{a_{1}b_{1}}(x_{1}\\!\mskip-1.5mu\smash{{}^{\prime}},\boldsymbol{z}_{1};\mu_{1},Y_{C})\underset{x_{1}}{\otimes}{}^{R}C_{a_{2}b_{2}}(x_{2}\\!\mskip-1.5mu\smash{{}^{\prime}},\boldsymbol{z}_{2};\mu_{2},Y_{C})$
$\displaystyle\qquad\underset{x_{2}}{\otimes}\lim_{Y_{L}\to-\infty}{}^{R}\mskip-1.5muF_{\text{us},b_{1}b_{2}}(x_{i}\\!\mskip-1.5mu\smash{{}^{\prime}},\boldsymbol{y};\mu_{i},Y_{L})\big{/}{}^{RR}\mskip-1.5mus(\boldsymbol{y};\mu_{i},Y_{C}-Y_{L})$
(164)

for a DPD in a right moving proton, where the short-distance coefficients are
defined as

$\displaystyle{}^{R}C_{ab}(x,\boldsymbol{z};\mu,Y_{C})$
$\displaystyle=\lim_{Y_{L}\to-\infty}\frac{{}^{R}C_{\text{us},ab}(x,\boldsymbol{z};\mu,Y_{L})}{\sqrt{{}^{R}C_{S,a}(\boldsymbol{z};\mu,2Y_{C}-2Y_{L})\rule{0.0pt}{8.1805pt}}}\,.$
(165)

Switching variables from $Y_{C}$ to $\zeta$ and again using the definition
(84) of DPDs, we obtain the final form of the matching equation:

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},\zeta)$
$\displaystyle=\sum_{b_{1},b_{2}}{}^{R}C_{a_{1}b_{1}}(x_{1}\\!\mskip-1.5mu\smash{{}^{\prime}},\boldsymbol{z}_{1};\mu_{1},x_{1}\zeta/x_{2})$
$\displaystyle\qquad\quad\underset{x_{1}}{\otimes}{}^{R}C_{a_{2}b_{2}}(x_{2}\\!\mskip-1.5mu\smash{{}^{\prime}},\boldsymbol{z}_{2};\mu_{2},x_{2}\zeta/x_{1})\underset{x_{2}}{\otimes}{}^{R}\mskip-1.5muF_{b_{1}b_{2}}(x_{i}\\!\mskip-1.5mu\smash{{}^{\prime}},\boldsymbol{y};\mu_{i},\zeta)\,.$
(166)

The rescaling factors $x_{1}/x_{2}$ or $x_{2}/x_{1}$ of $\zeta$ in the short-
distance coefficients arise for the same reason as discussed after (125). Note
that they involve the parton momentum fractions $x_{1}$ or $x_{2}$ and not the
integration variables $x_{1}^{\prime}$ or $x_{2}^{\prime}$ of the convolution
(10).

Let us emphasise that (6.2.1) involves mixing between quark and gluon
distributions. Therefore, the combination of
${}^{RR}\mskip-1.5mus(\boldsymbol{y})$ and
${}^{R}\mskip-1.5muF_{\text{us},b_{1}b_{2}}(x_{i},\boldsymbol{y})$ into a DPDF
works only because ${}^{RR}\mskip-1.5muS(\boldsymbol{y})$ is the same for
quarks and gluons (for representations $R$ accessible to quarks).

Let us finally set $R=1$ and show that in the colour singlet channel, the
coefficient ${}^{1}C_{ab}(x,\boldsymbol{z})$ in (6.2.1) is identical to the
coefficient $C_{ab}(x,\boldsymbol{z})$ in the small-$\boldsymbol{z}$ expansion
(9) of a single parton TMD. The corresponding identity for the coefficient
${}^{1}C_{\text{us},ab}(x,\boldsymbol{z})$ readily follows from taking the
matrix element of the operator relation (161) between two proton states, which
is related to the unsubtracted TMD $f_{\text{us},a}$ via the first relation in
(3.1). For ${}^{1}C_{S,a}(\boldsymbol{z})$ the identity is obtained by taking
the vacuum expectation value of (160), which gives the soft factor
$S_{a}(\boldsymbol{z})$ relevant for constructing the TMD
$f_{a}(x,\boldsymbol{z})$ via (58) and (60). The short-distance coefficient
for $f_{a}(x,\boldsymbol{z})$ is then obtained as in (165).

Notice that the vacuum expectation value of the colour projected operator on
the r.h.s. of (160) gives unity for $R=1$, because the Wilson lines along both
$v_{L}$ and $v_{R}$ appear in the combination $WW^{\dagger}=\mathbb{1}$. This
reflects the fact that no soft factor appears in the definition of ordinary
PDFs. We thus find

$\displaystyle
S_{a}(\boldsymbol{z};\mu,Y)={}^{1}C_{S,a}(\boldsymbol{z};\mu,Y)$ (167)

for small $\boldsymbol{z}$, i.e. the soft factor for single TMD factorisation
is purely perturbative at small distances.

####  6.2.2 Evolution equations and their solution

Let us now establish the consequences of the short-distance expansion (6.2.1)
for scale and rapidity evolution. From (163) we deduce that

$\displaystyle\frac{\partial}{\partial
Y}\log\sqrt{{}^{R}C_{S,a_{1}}(\boldsymbol{z}_{1};\mu_{1},2Y)}+\frac{\partial}{\partial
Y}\log\sqrt{{}^{R}C_{S,a_{2}}(\boldsymbol{z}_{2};\mu_{2},2Y)}$
$\displaystyle\quad=\frac{\partial}{\partial
Y}\log{}^{RR}\mskip-1.5mus_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},Y)-\frac{\partial}{\partial
Y}\log{}^{RR}\mskip-1.5mus(\boldsymbol{y};\mu_{i},Y)$
$\displaystyle\quad={}^{RR}\mskip-1.5muK_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i})-{}^{R}\\!J(\boldsymbol{y};\mu_{i})\phantom{\frac{\partial}{\partial}}$
(168)

with the kernels ${}^{RR}\mskip-1.5muK_{a_{1}a_{2}}$ and ${}^{R}\\!J$ from
(35) and (139). The rapidity derivative of
${}^{R}C_{S,a}(\boldsymbol{z};\mu,2Y)$ must thus be independent of $Y$.
Defining

$\displaystyle{}^{R}\mskip-1.5muK_{a}(\boldsymbol{z};\mu)$
$\displaystyle=\frac{\partial}{\partial
Y}\log{}^{R}C_{S,a}(\boldsymbol{z};\mu,Y)\,,$ (169)

we obtain

$\displaystyle\frac{\partial}{\partial\log\zeta}\,{}^{R}C_{ab}(x,\boldsymbol{z};\mu,\zeta)$
$\displaystyle=\frac{1}{2}\,{}^{R}\mskip-1.5muK_{a}(\boldsymbol{z};\mu)\,{}^{R}\mskip-1.5muC_{ab}(x,\boldsymbol{z};\mu,\zeta)$
(170)

for the coefficients in (165). Inserting this into (6.2.1), we find

$\displaystyle{}^{RR^{\prime}}\\!K_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i})$
$\displaystyle=\delta_{RR^{\prime}}\,\bigl{[}{}^{R}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1};\mu_{1})+{}^{R}\mskip-1.5muK_{a_{2}}(\boldsymbol{z}_{2};\mu_{2})+{}^{R}\\!J(\boldsymbol{y};\mu_{i})\bigr{]}\,.$
(171)

This represents a significant simplification of the Collins-Soper equation for
DTMDs in the small $\boldsymbol{z}_{i}$ limit, both in the colour structure
and in the separation of the three distance variables.

Taking the derivative $\partial/\partial\log\mu_{1}$ of (171) and using the
renormalisation group equations (129) and (140) for
${}^{RR^{\prime}}\\!K_{a_{1}a_{2}}$ and ${}^{R}\\!J$, we deduce that

$\displaystyle\gamma_{K,a}(\mu)$
$\displaystyle={}^{R}\gamma_{K,a}(\mu)+{}^{R}\gamma_{J}(\mu)$ (172)

for all $R$, where we have introduced the anomalous dimension of the kernel
${}^{R}\mskip-1.5muK_{a}$,

$\displaystyle\frac{\partial}{\partial\log\mu}\,{}^{R}\mskip-1.5muK_{a}(\boldsymbol{z};\mu)$
$\displaystyle={}-{}^{R}\gamma_{K,a}(\mu)\,.$ (173)

In the colour singlet sector, we have ${}^{1}\mskip-1.5muJ=0$ and hence
${}^{1\,}\mskip-1.5mu\gamma_{J}=0$, which implies
${}^{1\,}\mskip-1.5mu\gamma_{K,a}=\gamma_{K,a}$. Indeed,
${}^{1}\mskip-1.5muK_{a}(\boldsymbol{z};\mu)$ is identical to the ordinary
Collins-Soper kernel $K_{a}(\boldsymbol{z};\mu)$ in (7), which is readily
shown by taking the rapidity derivative of (167).

Taking the $\mu_{1}$ derivative of the short-distance expansion (6.2.1) and
using the evolution equations (124) and (141) for DTMDs and DPDFs, we obtain

$\displaystyle\frac{\partial}{\partial\log\mu}{}^{R}C_{ac}(x,\boldsymbol{z};\mu,\zeta)$
$\displaystyle\qquad=\sum_{b}{}^{R}C_{ab}(x^{\prime},\boldsymbol{z};\mu,\zeta)\underset{x}{\otimes}\Bigl{[}\delta_{bc}\,\delta(1-x^{\prime})\,\gamma_{F,c}(\mu,\zeta)-2\,{}^{R}\mskip-1.5muP_{bc}(x^{\prime};\mu,\zeta)\Bigr{]}\,.$
(174)

Having computed the short-distance kernel ${}^{R}C_{ab}$ at a certain order in
$\alpha_{s}$, we can use this relation to reconstruct the evolution kernel
${}^{R}\mskip-1.5muP_{ab}$ at the same order. We will do this in section
7.3.5.

With the short-distance limit (171) of the Collins-Soper kernel, the general
solution (5.1) of the evolution equations for DTMDs takes the form

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},\zeta)$
$\displaystyle\quad=\exp\,\biggl{\\{}\int_{\mu_{01}}^{\mu_{1}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{1}}(\mu)-\gamma_{K,a_{1}}(\mu)\log\frac{\sqrt{x_{1}\zeta/x_{2}}}{\mu}\biggr{]}$
$\displaystyle\hskip
35.50008pt+\int_{\mu_{02}}^{\mu_{2}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{2}}(\mu)-\gamma_{K,a_{2}}(\mu)\log\frac{\sqrt{x_{2}\zeta/x_{1}}}{\mu}\biggr{]}$
$\displaystyle\hskip
35.50008pt+\Bigl{[}{}^{R}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1};\mu_{01})+{}^{R}\mskip-1.5muK_{a_{2}}(\boldsymbol{z}_{2};\mu_{02})+{}^{R}\\!J(\boldsymbol{y};\mu_{0i})\Bigr{]}\log\frac{\sqrt{\zeta}}{\sqrt{\zeta_{0}}}\biggr{\\}}\,{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{0i},\zeta_{0})\,,$
(175)

where different colour channels no longer mix with each other. Using (170), we
can rewrite the rapidity dependence of the short-distance kernels as

$\displaystyle{}^{R}C_{a_{1}b_{1}}(x_{1}^{\prime},\boldsymbol{z}_{1};\mu_{01},x_{1}\zeta_{0}/x_{2})$
$\displaystyle=\exp\,\biggl{[}{}^{R}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1};\mu_{01})\log\frac{\sqrt{x_{1}\zeta_{0}/x_{2}}}{\mu_{01}}\biggr{]}\,{}^{R}C_{a_{1}b_{1}}(x_{1}^{\prime},\boldsymbol{z}_{1};\mu_{01},\mu_{01}^{2})$
(176)

and likewise for the index $2$, which reduces the number of independent scales
in $C_{ab}$. Using the short-distance matching (6.2.1) for
${}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{0i},\zeta_{0})$
in (6.2.2), we obtain our master formula for DTMDs at small
$\boldsymbol{z}_{i}$ and large $\boldsymbol{y}$:

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},\zeta)$
$\displaystyle\quad=\exp\,\biggl{\\{}\int_{\mu_{01}}^{\mu_{1}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{1}}(\mu)-\gamma_{K,a_{1}}(\mu)\log\frac{\sqrt{x_{1}\zeta/x_{2}}}{\mu}\biggr{]}+{}^{R}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1};\mu_{01})\log\frac{\sqrt{x_{1}\zeta/x_{2}}}{\mu_{01}}$
$\displaystyle\qquad\quad\;\;+\int_{\mu_{02}}^{\mu_{2}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{2}}(\mu)-\gamma_{K,a_{2}}(\mu)\log\frac{\sqrt{x_{2}\mskip
1.5mu\zeta/x_{1}}}{\mu}\biggr{]}+{}^{R}\mskip-1.5muK_{a_{2}}(\boldsymbol{z}_{2};\mu_{02})\log\frac{\sqrt{x_{2}\mskip
1.5mu\zeta/x_{1}}}{\mu_{02}}$
$\displaystyle\qquad\quad\;\;+{}^{R}\\!J(\boldsymbol{y};\mu_{0i})\log\frac{\sqrt{\zeta}}{\sqrt{\zeta_{0}}}\biggr{\\}}$
$\displaystyle\quad\;\times\sum_{b_{1},b_{2}}{}^{R}C_{a_{1}b_{1}}(x_{1}^{\prime},\boldsymbol{z}_{1};\mu_{01},\mu_{01}^{2})\underset{x_{1}}{\otimes}{}^{R}C_{a_{2}b_{2}}(x_{2}^{\prime},\boldsymbol{z}_{2};\mu_{02},\mu_{02}^{2})\underset{x_{2}}{\otimes}{}^{R}\mskip-1.5muF_{b_{1}b_{2}}(x_{i}^{\prime},\boldsymbol{y};\mu_{0i},\zeta_{0}).$
(177)

In the colour singlet channel, the result (6.2.2) is a simple copy of its
analogue (2) for a single-parton TMD, with separate short-distance
coefficients $C_{ab}$ and Sudakov exponentials for each parton. For colour
non-singlet channels, the Collins-Soper kernels
${}^{R}\mskip-1.5muK_{a}(\boldsymbol{z})$ acquire a colour dependence, and
there is an additional term ${}^{R}\\!J(\boldsymbol{y})$ in the exponent. The
latter is also present in collinear DPS factorisation: in fact one simply has

$\displaystyle\exp\,\biggl{[}{}^{R}\\!J(\boldsymbol{y};\mu_{0i})\,\log\frac{\sqrt{\zeta}}{\sqrt{\zeta_{0}}}\,\biggr{]}\,{}^{R}\mskip-1.5muF_{b_{1}b_{2}}(x_{i}^{\prime},\boldsymbol{y};\mu_{0i},\zeta_{0})$
$\displaystyle={}^{R}\mskip-1.5muF_{b_{1}b_{2}}(x_{i}^{\prime},\boldsymbol{y};\mu_{0i},\zeta)$
(178)

according to (139). In this sense, the value of $\zeta_{0}$ in (6.2.2) is
irrelevant (contrary to the choice of $\mu_{0i}$, which appears in quantities
that are computed using fixed-order perturbation theory). The choice of
$\zeta_{0}$ does matter when this is the scale at which one formulates an
ansatz or model for the DPDFs (which is of course inevitable for concrete
calculations).

After inserting (6.2.2) and its counterpart for the left moving proton into
the definition (5) of $W$, we can combine the logarithms of $\zeta$ and
$\bar{\zeta}$ in the exponentials, as we did in the generic expression (5.1).
We thus obtain

$\displaystyle W_{\text{large $y$}}=\sum_{R}\eta_{a_{1}a_{2}}(R)$
$\displaystyle\qquad\times\exp\,\biggl{\\{}\int_{\mu_{01}}^{\mu_{1}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{1}}(\mu)-\gamma_{K,a_{1}}(\mu)\,\log\frac{{Q_{1}^{2}}}{\mu^{2}}\biggr{]}+{}^{R}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1};\mu_{01})\log\frac{{Q_{1}^{2}}}{\mu_{01}^{2}}$
$\displaystyle\hskip
44.50003pt+\int_{\mu_{02}}^{\mu_{2}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{2}}(\mu)-\gamma_{K,a_{2}}(\mu)\,\log\frac{{Q_{2}^{2}}}{\mu^{2}}\biggr{]}+{}^{R}\mskip-1.5muK_{a_{2}}(\boldsymbol{z}_{2};\mu_{02})\,\log\frac{{Q_{2}^{2}}}{\mu_{02}^{2}}\biggr{\\}}$
$\displaystyle\qquad\times\sum_{c_{1},c_{2},d_{1},d_{2}}{}^{R}C_{b_{1}d_{1}}(\bar{x}_{1}^{\prime},\boldsymbol{z}_{1};\mu_{01},\mu_{01}^{2})\underset{\bar{x}_{1}}{\otimes}{}^{R}C_{b_{2}d_{2}}(\bar{x}_{2}^{\prime},\boldsymbol{z}_{2};\mu_{02},\mu_{02}^{2})\,$
$\displaystyle\hskip
51.00014pt\underset{\bar{x}_{2}}{\otimes}{}^{R}C_{a_{1}c_{1}}(x_{1}^{\prime},\boldsymbol{z}_{1};\mu_{01},\mu_{01}^{2})\underset{x_{1}}{\otimes}{}^{R}C_{a_{2}c_{2}}(x_{2}^{\prime},\boldsymbol{z}_{2};\mu_{02},\mu_{02}^{2})$
$\displaystyle\qquad\quad\underset{x_{2}}{\otimes}\bigl{[}\Phi(\nu\boldsymbol{y})\bigr{]}^{2}\,\exp\,\biggl{[}{}^{R}\\!J(\boldsymbol{y};\mu_{0i})\mskip
1.5mu\log\frac{Q_{1}\mskip
1.5muQ_{2}}{\zeta_{0}}\,\biggr{]}\,{}^{R}\mskip-1.5muF_{d_{1}d_{2}}(\bar{x}_{i}^{\prime},\boldsymbol{y};\mu_{0i},\zeta_{0})\,{}^{R}\mskip-1.5muF_{c_{1}c_{2}}(x_{i}^{\prime},\boldsymbol{y};\mu_{0i},\zeta_{0})\,,$
(179)

which is the main result of this section. The product of regulator functions
$\Phi(\nu\boldsymbol{y}_{+})\mskip 1.5mu\Phi(\nu\boldsymbol{y}_{-})$ has been
approximated as appropriate for
$|\boldsymbol{z}_{1}|,|\boldsymbol{z}_{2}|\ll|\boldsymbol{y}|$. Note that the
dependence on $\boldsymbol{y}$, $\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$
is completely factorised in (6.2.2). The $\boldsymbol{y}$ integral can hence
be performed separately. For unpolarised or longitudinally polarised partons,
the short-distance coefficients $C$ are independent of the direction of
$\boldsymbol{z}_{1},\boldsymbol{z}_{2}$ by rotation invariance. The angular
part of the $\boldsymbol{z}_{i}$ integrations can then readily be performed
and turns the Fourier exponentials
$e^{-i\boldsymbol{q}{}_{i}\boldsymbol{z}_{i}}$ into Bessel functions. The
situation for transverse quark or linear gluon polarisation, where $C$ (as
well as $F$ and $\hat{\sigma}$) carries transverse indices, can be discussed
along similar lines.

Each of the two partonic cross sections $\hat{\sigma}_{i}$ in (5) and the four
short-distance kernels ${}^{R}C$ in (6.2.2) has an $\alpha_{s}$ expansion
starting at order $\alpha_{s}^{0}$. It is natural (although not mandatory) to
truncate their product to the highest order in $\alpha_{s}$ at which the
individual factors are computed.

To evaluate the kernels ${}^{R}\mskip-1.5muK$ and ${}^{R}C$ in (6.2.2) in
fixed-order perturbation theory, one should choose the scales $\mu_{0i}$ such
that no large logarithms appear at higher orders. A standard choice for single
hard scattering is to take $\mu_{0i}^{2}=b_{0}^{2}\mskip
1.5mu/\boldsymbol{z}_{i}^{2}$, which makes the one-loop expression of
${}^{R}\mskip-1.5muK(\boldsymbol{z}_{i};\mu_{0i})$ vanish. Here

$\displaystyle b_{0}$ $\displaystyle=2e^{-\gamma_{E}}$ (180)

with $\gamma_{E}$ being the Euler constant. While the cross section (5) is
dominated by $|\boldsymbol{z}_{i}|\sim 1/|\boldsymbol{q}_{i}|$ due to the
Fourier exponentials $e^{-i\boldsymbol{q}{}_{i}\boldsymbol{z}_{i}}$, one still
must integrate over the full range of these distances. To avoid evaluating the
DPDs at unreasonably small factorisation scales $\mu_{0i}$, one may modify the
above scale choice and take instead $\mu_{0i}=\mu_{z_{i}}$, where for any
transverse distance vector $\boldsymbol{b}$ we define

$\displaystyle\mu_{b}^{2}$
$\displaystyle=\frac{b_{0}^{2}}{\boldsymbol{b}^{*2}}\,,$
$\displaystyle\text{with \leavevmode\nobreak\
$\boldsymbol{b}^{*}(\boldsymbol{b})\xrightarrow[\boldsymbol{b}\to\boldsymbol{0}]{}\boldsymbol{b}$\leavevmode\nobreak\
and \leavevmode\nobreak\ $|\mskip
1.5mu\boldsymbol{b}^{*}(\boldsymbol{b})|\xrightarrow[|\boldsymbol{b}|\to\infty]{}b_{\text{max}}$}\,,$
(181)

where $b_{\text{max}}$ is chosen such that $\boldsymbol{b}^{*}$ remains in the
region where one trusts perturbative theory even when $\boldsymbol{b}$ becomes
large. A possible choice for $\boldsymbol{b}^{*}$ is

$\displaystyle\boldsymbol{b}^{*}$
$\displaystyle=\frac{\boldsymbol{b}}{\sqrt{1+\boldsymbol{b}^{2}/b_{\text{max}}^{2}}}\,,$

which was proposed long ago Collins:1981va ; Collins:1984kg  and is
extensively used in TMD phenomenology, but other functional forms have also
been explored Bacchetta:2015ora . We note that it is not mandatory to take
factorisation scales proportional to inverse distances, referring to
Becher:2011xn ; Becher:2012yn  and DAlesio:2014mrz ; Echevarria:2015uaa  for
examples of scale setting in momentum space and to section 7.1.2 of
Ebert:2016gcn  for a brief synopsis. Most results in our paper do not depend
on choosing scales as in (181). An exception is the discussion of scale
setting in section 6.4 and appendix C, which needs to be adapted if a
different choice is made.

In the colour nonsinglet sector, the DPDs ${}^{R}\mskip-1.5muF$ in (6.2.2)
still involve the different scales $\mu_{0i}$ and $\zeta_{0}$. To disentangle
this dependence further, one may use (5.2) and (178) to rewrite
${}^{R}\mskip-1.5muF$ in terms of the distribution
${}^{R}\mskip-1.5mu\widehat{F}$, whose DGLAP evolution does not involve a
separate rapidity scale, and which is therefore better suited to make the
separation of scales explicit. At the cross section level, one can then
replace the last line of (6.2.2) by

$\displaystyle\bigl{[}\Phi(\nu\boldsymbol{y})\bigr{]}^{2}\,\exp\,\biggl{[}-\int_{\mu_{0}}^{\mu_{01}}\frac{d\mu}{\mu}\;{}^{R}\gamma_{J}(\mu)\,\log\frac{{Q_{1}^{2}}}{\mu^{2}}-\int_{\mu_{0}}^{\mu_{02}}\frac{d\mu}{\mu}\;{}^{R}\gamma_{J}(\mu)\,\log\frac{{Q_{2}^{2}}}{\mu^{2}}$
$\displaystyle\qquad+{}^{R}\\!J(\boldsymbol{y};\mu_{0},\mu_{0})\mskip
1.5mu\log\frac{Q_{1}\mskip
1.5muQ_{2}}{\zeta_{0}}\,\biggr{]}\,{}^{R}\mskip-1.5mu\widehat{F}_{d_{1}d_{2},\,\mu_{0},\zeta_{0}}(\bar{x}_{i}^{\prime},\boldsymbol{y};\mu_{0i})\,{}^{R}\mskip-1.5mu\widehat{F}_{c_{1}c_{2},\,\mu_{0},\zeta_{0}}(x_{i}^{\prime},\boldsymbol{y};\mu_{0i})\,,$
(182)

where we have explicitly indicated that the two renormalisation scales in
${}^{R}\\!J$ are taken equal. It is natural to take the starting scale of
evolution as $\mu_{0}=\sqrt{\zeta_{0}}$. A particular choice is
$\mu_{0}=\mu_{y}$, which for large $\boldsymbol{y}$ saturates at a hadronic
scale $\mu_{0}=b_{0}/b_{\text{max}}$. This choice ensures that the initial
conditions for DPD evolution do not involve widely different scales as
$\boldsymbol{y}$ becomes small.

The factor ${}^{R}\\!J(\boldsymbol{y};\mu_{0i})$ in (6.2.2) is a
nonperturbative function in the large-$\boldsymbol{y}$ region, but we can
ensure that it has the correct perturbative small-$\boldsymbol{y}$ behaviour
by copying part of the so-called $b^{*}$ trick formulated for ordinary TMDs in
Collins:1981va ; Collins:1984kg . We thus write

$\displaystyle{}^{R}\\!J(\boldsymbol{y};\mu_{01},\mu_{02})$
$\displaystyle=-\,{}^{R}\mskip-1.5mug_{J}(\boldsymbol{y})+{}^{R}\\!J(\boldsymbol{y}^{*}\mskip-1.5mu;\mu_{y},\mu_{y})-\int_{\mu_{y}}^{\mu_{01}}\frac{d\mu}{\mu}\;{}^{R}\gamma_{J}(\mu)-\int_{\mu_{y}}^{\mu_{02}}\frac{d\mu}{\mu}\;{}^{R}\gamma_{J}(\mu)\,,$
(183)

where the nonperturbative information is contained in

$\displaystyle{}^{R}\mskip-1.5mug_{J}(\boldsymbol{y})$
$\displaystyle={}^{R}\\!J(\boldsymbol{y}^{*}\mskip-1.5mu;\mu_{i})-{}^{R}\\!J(\boldsymbol{y};\mu_{i})\,.$
(184)

Because of (140), the dependence on the renormalisation scales drops out in
${}^{R}\mskip-1.5mug_{J}$, which satisfies
${}^{R}\mskip-1.5mug_{J}(\boldsymbol{0})=0$ by construction. The term
${}^{R}\\!J(\boldsymbol{y}^{*}\\!;\mu_{y},\mu_{y})$ in (183) can be evaluated
in fixed-order perturbation theory; according to the results in section 7.2 it
is zero up to terms of $\mathcal{O}(\alpha_{s}^{2})$. Of course, the form
(183) with $\mu_{0i}$ replaced by $\mu_{0}$ can be used in (6.2.2).

In the colour octet sector, the relation (110) implies that

$\displaystyle{}^{8}\mskip-1.5muJ(\boldsymbol{y};\mu,\mu)$
$\displaystyle=K_{g}(\boldsymbol{y};\mu)\,,$ (185)

where $K_{g}$ is the Collins Soper kernel for single gluon TMDs. While little
is known about this kernel for nonperturbative distances $\boldsymbol{y}$, one
may construct a model by connecting it with its counterpart $K_{q}$ in the
quark-antiquark channel, for which there is a considerable body of
phenomenology from single Drell-Yan production and from semi-inclusive DIS,
see e.g. Aidala:2014hva ; Collins:2014jpa  and references therein. A simple
possibility would be to assume Casimir scaling, $K_{g}/C_{A}=K_{q}/C_{F}$.
This scaling holds perturbatively up to $\mathcal{O}(\alpha_{s}^{3})$, as can
for instance be seen in appendix D of Echevarria:2016scs .101010Just recently
it has been found Moch:2017uml  that Casimir scaling for the cusp anomalous
dimension, $\gamma_{K,g}/C_{A}=\gamma_{K,q}/C_{F}$, is broken at
$\mathcal{O}(\alpha_{s}^{4})$. Equation (8) implies the breaking of Casimir
scaling for $K_{a}$ at the same order. For a related calculation see
Grozin:2017css .

####  6.2.3 Extrapolation to large $\boldsymbol{z}_{1}$ and
$\boldsymbol{z}_{2}$

If $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$ are not sufficiently large to
ensure dominance of small $\boldsymbol{z}_{i}$ in the TMD cross section, then
a more realistic description of the integrand for large $\boldsymbol{z}_{i}$
is necessary. Of course, a corresponding statement holds already for single
hard scattering. A widely used procedure in that case is the $b^{*}$ trick of
Collins:1981va ; Collins:1984kg . Let us briefly show how it can be adapted to
DPS. We recall the solution (131) of the rapidity evolution equation,

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},\zeta)$
$\displaystyle=\sum_{R^{\prime}}{}^{RR^{\prime}}\\!\exp\biggl{[}K_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i})\log\frac{\sqrt{\zeta}}{\sqrt{\zeta_{0}}}\,\biggr{]}{}^{R^{\prime}}\\!F_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},\zeta_{0})$

and introduce

$\displaystyle{}^{R}\mskip-1.5mug_{F,a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\zeta)$
$\displaystyle=\log\,\frac{{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}^{*}_{i},\boldsymbol{y};\mu_{i},\zeta)}{{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},\zeta)}\,,$
$\displaystyle{}^{RR^{\prime}}\\!g_{K,a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle={}^{RR^{\prime}}\\!K_{a_{1}a_{2}}(\boldsymbol{z}^{*}_{i},\boldsymbol{y};\mu_{i})-{}^{RR^{\prime}}\\!K_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i})$
(186)

with $\boldsymbol{z}^{*}$ defined by (181). Both $g_{F}$ and $g_{K}$ are
independent of $\mu_{i}$ according to (124) and (129), and both functions
depend on the parameter $b_{\text{max}}$ via $\boldsymbol{z}^{*}_{i}$. By
construction, both functions vanish at the point
$\boldsymbol{z}_{1}=\boldsymbol{z}_{2}=\boldsymbol{0}$. We can then write

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},\zeta)$
$\displaystyle=\sum_{R^{\prime}}{}^{RR^{\prime}}\\!\exp\biggl{[}{}-g_{K,a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})\log\frac{\sqrt{\zeta}}{\sqrt{\zeta_{0}}}\,+K_{a_{1}a_{2}}(\boldsymbol{z}_{i}^{*},\boldsymbol{y};\mu_{i})\log\frac{\sqrt{\zeta}}{\sqrt{\zeta_{0}}}\,\biggr{]}$
$\displaystyle\quad\times\exp\bigl{[}{}-{}^{R^{\prime}}\\!g_{F,a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\zeta_{0})\bigr{]}\,{}^{R^{\prime}}\\!F_{a_{1}a_{2}}(x_{i},\boldsymbol{z}^{*}_{i},\boldsymbol{y};\mu_{i},\zeta_{0})\,.$
(187)

In this expression, ${}^{RR^{\prime}}\\!K_{a_{1}a_{2}}$ and
${}^{R^{\prime}}\\!F_{a_{1}a_{2}}$ can be evaluated using the short-distance
expansion discussed in the previous subsection, provided that
$|\boldsymbol{z}_{1}^{*}|,|\boldsymbol{z}_{2}^{*}|\ll|\boldsymbol{y}|$. As a
consequence, ${}^{RR^{\prime}}\\!K_{a_{1}a_{2}}$ is diagonal in $R$ and
$R^{\prime}$. By contrast, the functions ${}^{RR^{\prime}}\\!g_{K,a_{1}a_{2}}$
and ${}^{R}\mskip-1.5mug_{F,a_{1}a_{2}}$ are entirely nonperturbative and need
to be modelled. This is a daunting task, because they depend on several
transverse variables, there are many functions in the different colour
channels and there is less guidance from data. Using the $b^{*}$ trick for
computing DPS processes with $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$ in
the nonperturbative region would therefore require additional theory input, or
strong simplifying assumptions.

The $b^{*}$ trick is not the only way to handle the region of large transverse
distances, and for TMD factorisation in SPS a variety of other methods have
been employed Laenen:2000de ; Kulesza:2002rh ; Bozzi:2005wk ; Qiu:2000hf ;
Becher:2011xn ; Becher:2012yn ; DAlesio:2014mrz ; Echevarria:2015uaa ;
Scimemi:2016ffw ; Scimemi:2017etj . It would be interesting to study if and
how they could be adapted to DPS. This is however beyond the scope of the
present work, where we focus on transverse momenta in the perturbative region.

###  6.3 The small-$\boldsymbol{y}$ region

We now consider the region where $|\boldsymbol{y}|$ is of the same order as
$|\boldsymbol{z}_{1}|$ and $|\boldsymbol{z}_{2}|$, with all distances being
small compared with $1/\Lambda$. Again, we will derive double parton analogues
of the relations (9) and (2) between TMDs and PDFs. They are given in (193),
(194) and (6.3.2) and involve collinear PDFs and twist-four distributions as
nonperturbative input, which also appear in the expansion (6.3.1) of DPDFs at
small $\boldsymbol{y}$. At the cross section level, we then obtain the
expression (6.3.2) in which all large logarithms are explicitly resummed. An
important finding is that the different combinations of collinear PDFs and
twist-four distributions have different power behaviour in the cross section,
specified in (195) and compared with the contribution from large
$\boldsymbol{y}$ in (196).

####  6.3.1 Short-distance expansion

Let us first establish that the soft factor, and hence the Collins-Soper
kernel ${}^{RR^{\prime}}\\!K_{a_{1}a_{2}}$, can be fully computed in
perturbation theory when all transverse distances become small. Repeating the
arguments that lead to (157), we obtain

$\displaystyle(S_{a_{1}a_{2}})^{\underline{r}_{1}\underline{r}_{2},\underline{u}_{1}\underline{u}_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=C_{S,a_{1}a_{2}\phantom{|}}^{\underline{s}_{1}\underline{s}_{2},\underline{t}_{1}\underline{t}_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})\,(S_{a_{1}a_{2}})^{\underline{r}_{1}\underline{r}_{2},\underline{u}_{1}\underline{u}_{2}}_{\underline{s}_{1}\underline{s}_{2},\underline{t}_{1}\underline{t}_{2}}(\boldsymbol{0})\,,$
(188)

which involves the extended soft factor (3.2) with all distance arguments set
to zero and a short-distance coefficient $C_{S,a_{1}a_{2}}$ describing
interactions between all Wilson lines. Inserting colour projectors as in (73),
we get

$\displaystyle(S_{a_{1}a_{2}})^{\underline{r}_{1}\underline{r}_{2},\underline{u}_{1}\underline{u}_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\sum_{RR^{\prime}}P_{R}^{\underline{v}_{1}\underline{v}_{2}}\,P_{R^{\prime}}^{\underline{w}_{1}\underline{w}_{2}}\,C_{S,a_{1}a_{2}\phantom{|}}^{\underline{v}_{1}\underline{v}_{2},\underline{w}_{1}\underline{w}_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})\,\frac{P_{R}^{\underline{s}_{1}\underline{s}_{2}}\,P_{R^{\prime}}^{\underline{t}_{1}\underline{t}_{2}}\,(S_{a_{1}a_{2}})^{\underline{r}_{1}\underline{r}_{2},\underline{u}_{1}\underline{u}_{2}}_{\underline{s}_{1}\underline{s}_{2},\underline{t}_{1}\underline{t}_{2}}(\boldsymbol{0})}{m(R)\,m(R^{\prime})}$
$\displaystyle=\sum_{RR^{\prime}}P_{R}^{\underline{v}_{1}\underline{v}_{2}}\,P_{R^{\prime}}^{\underline{w}_{1}\underline{w}_{2}}\,C_{S,a_{1}a_{2}\phantom{|}}^{\underline{v}_{1}\underline{v}_{2},\underline{w}_{1}\underline{w}_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})\,\frac{P_{R}^{\underline{r}_{1}\underline{s}_{1}}\,P_{R^{\prime}}^{\underline{u}_{1}\underline{t}_{1}}\,(S_{a_{1}a_{2}})^{\underline{s}_{1}\underline{r}_{2},\underline{t}_{1}\underline{u}_{2}}_{\underline{s}_{2}\underline{s}_{2},\underline{t}_{2}\underline{t}_{2}}(\boldsymbol{0})}{m(R)\,m(R^{\prime})}\,,$
(189)

where in the second step the identity (102) has been used for the indices
associated with parton $a_{1}$ in $S_{a_{1}a_{2}}$. The sum over
$\underline{s}_{2}$ on the r.h.s. now ties together pairs of conjugate Wilson
lines along the same direction and at the same transverse position, which
results in unit matrices according to
$W(\boldsymbol{\xi},v)W^{\dagger}(\boldsymbol{\xi},v)=\mathbb{1}$. For
$S_{qq}$ this is represented pictorially in figure 8. The same holds for the
sum over $\underline{t}_{2}$. The extended soft factor thus collapses to a
product of unit matrices, and we have

$\displaystyle(S_{a_{1}a_{2}})^{\underline{r}_{1}\underline{r}_{2},\underline{u}_{1}\underline{u}_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=\sum_{RR^{\prime}}P_{R}^{\underline{v}_{1}\underline{v}_{2}}\,P_{R^{\prime}}^{\underline{w}_{1}\underline{w}_{2}}\,C_{S,a_{1}a_{2}\phantom{|}}^{\underline{v}_{1}\underline{v}_{2},\underline{w}_{1}\underline{w}_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})\frac{P_{R}^{\underline{r}_{1}\underline{r}_{2}}\,P_{R^{\prime}}^{\underline{u}_{1}\underline{u}_{2}}}{m(R)\,m(R^{\prime})}\,.$
(190)

Projecting this on definite representations for
$(\underline{r}_{1}\underline{r}_{2})$ and for
$(\underline{u}_{1}\underline{u}_{2})$, we simply obtain

$\displaystyle{}^{RR^{\prime}}\\!S_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},Y)$
$\displaystyle={}^{RR^{\prime}}\\!C_{S,a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},Y)\,,$
(191)

where the projection ${}^{RR^{\prime}}\mskip-1.5muC$ is defined in analogy to
(81) and where we have restored all arguments. In analogy to the case (167)
for SPS, we thus find that when all transverse distances become small, the
soft factor for DPS is equal to its matching coefficient and thus can be
entirely computed in perturbation theory. The same then holds for the Collins-
Soper kernel ${}^{RR^{\prime}}\\!K$ and for the matrix ${}^{RR^{\prime}}\\!s$,
which can be constructed from ${}^{RR^{\prime}}\\!S$ as shown in appendix A.
An explicit check of this result is the fact that the one-loop expression for
${}^{RR^{\prime}}\\!K$ in section 7.2 is free of infrared divergences, as is
the two-loop soft factor for DPS computed in Vladimirov:2016qkd .

Figure 8:  Simplification in the extended soft factor $S_{qq}$ with all Wilson
lines at the same transverse position, using the identity in figure 2 and the
unitarity of Wilson lines. Analogous relations hold for indices in the adjoint
colour representation.

We now turn to the short-distance limit of the collinear matrix elements,
which can be discussed in close analogy to the small-$\boldsymbol{y}$ limit of
collinear DPDs in section 3.3 of Diehl:2017kgu . In the limit of small
$\boldsymbol{z}_{i}$ and $\boldsymbol{y}$, the unsubtracted distributions
${}^{R^{\prime}}\\!F_{\text{us}}$ are given by short-distance coefficients
times proton matrix elements of operators with all fields at transverse
position $\boldsymbol{0}$. The latter are given by collinear parton
distributions of different twist. Multiplying with the perturbative expression
of ${}^{RR^{\prime}}\\!s$, one obtains the expansion of ${}^{R}\mskip-1.5muF$.
Adapting the notation in Diehl:2017kgu , we write111111The distributions on
the r.h.s. of (192) are respectively denoted by $F_{\text{spl,pt}}$,
$F_{\text{int,pt}}$ and $F_{\text{tw3,pt}}$ in Diehl:2017kgu . For brevity we
omit the specification “pt” here.

$\displaystyle{}^{R}\mskip-1.5muF$
$\displaystyle={}^{R}\mskip-1.5muF_{\text{spl}}+{}^{R}\mskip-1.5muF_{\text{tw3}}+{}^{R}\mskip-1.5muF_{\text{int}}\,,$
(192)

where the short-distance expansion of the terms on the r.h.s. involves proton
matrix elements of operators with twist two, twist three and twist four,
respectively. Following the discussion in section 2.1 of Diehl:2017kgu , we
neglect the contribution with collinear twist-three distributions, which for
unpolarised protons are chiral odd. They have no dynamic cross talk with gluon
distributions and are therefore expected to be small compared with twist-two
and twist-four distributions at small $x$.

The splitting contribution $F_{\text{spl}}$ describes the case where a single
parton splits into partons $a_{1}$ and $a_{2}$, as shown in figure 6a at
lowest order in $\alpha_{s}$. It can be obtained from (154) by taking the
limit of small $\boldsymbol{Z}$. For the $\mathcal{O}(\alpha_{s})$ term, this
simply corresponds to replacing the TMD
$f_{a_{0}}(x_{1}+x_{2},\boldsymbol{Z})$ by the PDF $f_{a_{0}}(x_{1}+x_{2})$.
Since the proton is unpolarised, the parton $a_{0}$ in the PDF is unpolarised
as well. We thus obtain

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2},\,\text{spl}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu,\mu,\zeta)$
$\displaystyle=\frac{\boldsymbol{y}_{+}^{l}\mskip
1.5mu\boldsymbol{y}_{-}^{l^{\prime}}}{\boldsymbol{y}_{+}^{2}\mskip
1.5mu\boldsymbol{y}_{-}^{2}}\,\frac{\alpha_{s}(\mu)}{2\pi^{2}}\;{}^{R}T_{a_{0}\to
a_{1}a_{2}}^{ll^{\prime}}\biggl{(}\frac{x_{1}}{x_{1}+x_{2}}\biggr{)}\,\frac{f_{a_{0}}(x_{1}+x_{2};\mu)}{x_{1}+x_{2}}+\mathcal{O}(\alpha_{s}^{2})\,,$
(193)

where we have restored the dependence on the UV renormalisation scale, taken
equal for the partons with momentum fraction $x_{1}$ and $x_{2}$. A $\zeta$
dependence appears only at order $\alpha_{s}^{2}$. To avoid large logarithms
of the higher-order terms in the region of small
$|\boldsymbol{z}_{1}|\sim|\boldsymbol{z}_{2}|\sim|\boldsymbol{y}|$, any choice
where $1/\mu$ and $1/\sqrt{\zeta}$ are of the order of these distances will
do. A particular choice is given in (200) below.

The term $F_{\text{int}}$ in (192) is referred to as the “intrinsic”
contribution to the DPD and may be thought of as describing parton pairs
$a_{1}$, $a_{2}$ in the “intrinsic” proton wave function. It is the only
contribution starting at order $\alpha_{s}^{0}$ and reads

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2},\,\text{int}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu,\mu,\zeta)$
$\displaystyle={}^{R}G_{a_{1}a_{2}}(x_{1},x_{2},x_{2},x_{1};\mu)+\mathcal{O}(\alpha_{s})\,,$
(194)

where ${}^{R}G$ denotes a collinear twist-four distribution. The lowest-order
term is simply obtained by setting
$\boldsymbol{z}_{i}=\boldsymbol{y}=\boldsymbol{0}$ in the matrix element (3.1)
(this must be done before renormalisation and hence in $D=4-2\epsilon$
dimensions). The $\mathcal{O}(\alpha_{s})$ term involves short-distance
interactions and is briefly discussed in section 3.3 of Diehl:2017kgu .

Inserting (192) with (193) and (194) into the cross section formula (5), we
can derive the power behaviour of the contribution from the region
$|\boldsymbol{z}_{1}|\sim|\boldsymbol{z}_{2}|\sim|\boldsymbol{y}|\sim
1/q_{T}$, using that the twist-four distribution scales like
$G\sim\Lambda^{2}$. For the scaled DPS cross section we obtain

$\displaystyle
Q^{4}\,\frac{d\sigma_{\text{DPS}}}{d^{2}\boldsymbol{q}_{1}\,d^{2}\boldsymbol{q}_{2}}\bigg{|}_{\text{small
$y$}}$ $\displaystyle\sim\begin{cases}\alpha_{s}^{2}/q_{T}^{2}&\text{from
$F_{\text{spl}}\times F_{\text{spl}}$ \leavevmode\nobreak\ (1v1)}\\\
\alpha_{s}\mskip 1.5mu\Lambda^{2}/q_{T}^{4}&\text{from $F_{\text{spl}}\times
F_{\text{int}}$ \leavevmode\nobreak\ (2v1)}\\\
\Lambda^{4}/q_{T}^{6}&\text{from $F_{\text{int}}\times F_{\text{int}}$
\leavevmode\nobreak\ (2v2)}\\\ \end{cases}$ (195) where we have also indicated
the lowest order in $\alpha_{s}$ for each contribution. For the contribution
from large $|\boldsymbol{y}|\sim 1/\Lambda$ discussed in section 6.2, we have
$\displaystyle
Q^{4}\,\frac{d\sigma_{\text{DPS}}}{d^{2}\boldsymbol{q}_{1}\,d^{2}\boldsymbol{q}_{2}}\bigg{|}_{\text{large
$y$}}$ $\displaystyle\sim\Lambda^{2}/q_{T}^{4}\,,$ (196)

as follows from the small-distance expansion (6.2.1) and the power behaviour
$F\sim\Lambda^{2}$ of a DPDF at large $\boldsymbol{y}$. We note that (195) and
(196) do not hold in the kinematic region
$|\boldsymbol{q}_{1}+\boldsymbol{q}_{2}|\sim\Lambda$, which we have excluded
from our considerations, and which has a different power counting as seen from
table 1 in Diehl:2011yj . Notice also that some parton combinations, such as
$(a_{1}a_{2})=(d\bar{u})$, cannot be obtained by parton splitting at
$\mathcal{O}(\alpha_{s})$. In this case, $F_{\text{spl}}$ starts at
$\mathcal{O}(\alpha_{s}^{2})$ and the powers of $\alpha_{s}$ in (195) must be
adjusted accordingly.

The three contributions in (195) are referred to as 1v1, 2v1 and 2v2 terms as
indicated. The 1v1 (read “one versus one”) term is associated with graphs as
in figure 5a, where one parton in each proton initiates a sequence of short-
distance interactions (a splitting at scale $q_{T}$ followed by hard
scattering at scale $Q$). The 2v2 term corresponds to graphs where two
“intrinsic” partons in each proton directly enter the hard subprocesses, as
shown in figure 9a. In the 2v1 term one starts with two “intrinsic” partons in
one proton, whereas in the other proton one starts with one parton that splits
into two at scale $q_{T}$. This is shown in figure 9b. It is natural to
associate the term in (196) with 2v2 as well, since at large $\boldsymbol{y}$
one cannot identify any perturbative splitting contribution.

Figure 9:  Counterparts to the 1v1 graph in figure 5a. In the 2v2 graph (a)
there is no perturbative splitting, and in the 2v1 graph (b) one has a
perturbative splitting in one of the two protons. The box encloses the DTMD of
the upper proton and indicates that the splittings take place at transverse
distances of order $1/q_{T}$.

Among the terms in (195) and (196), only the term with $F_{\text{spl}}\times
F_{\text{spl}}$ is leading as far as powers of $\Lambda^{2}/q_{T}^{2}$ are
concerned. However, the other terms involve fewer powers in $\alpha_{s}$, and
they rise more strongly as the momentum fractions $x$ become small, because
one roughly expects the intrinsic part of the DPDF to grow like the square of
a single PDF in that limit. From this perspective, one may still discard the
term $F_{\text{int}}\times F_{\text{int}}$, which is power suppressed compared
with the large-$\boldsymbol{y}$ contribution but has the same power in
$\alpha_{s}$ and a similar small $x$ behaviour. However, the inclusion of this
term in the cross section will in general not require much additional effort,
and it renders the combination of the contributions from small and large
$\boldsymbol{y}$ more straightforward (see section 6.4).

The expansion of the DPDFs $F(x_{i},\boldsymbol{y})$ for small
$\boldsymbol{y}$ proceeds in exactly the same manner. It can be obtained from
(193) and (194) by setting
$\boldsymbol{z}_{1}=\boldsymbol{z}_{2}=\boldsymbol{0}$ before performing UV
renormalisation and setting $D=4$. This gives

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2},\,\text{spl}}(x_{i},\boldsymbol{y};\mu,\mu,\zeta)$
$\displaystyle=\frac{\boldsymbol{y}^{l}\mskip
1.5mu\boldsymbol{y}^{l^{\prime}}}{\boldsymbol{y}^{4}}\,\frac{\alpha_{s}(\mu)}{2\pi^{2}}\;{}^{R}T_{a_{0}\to
a_{1}a_{2}}^{ll^{\prime}}\biggl{(}\frac{x_{1}}{x_{1}+x_{2}}\biggr{)}\,\frac{f_{a_{0}}(x_{1}+x_{2};\mu)}{x_{1}+x_{2}}+\mathcal{O}(\alpha_{s}^{2})\,,$
$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2},\,\text{int}}(x_{i},\boldsymbol{y};\mu,\mu,\zeta)$
$\displaystyle={}^{R}G_{a_{1}a_{2}}(x_{1},x_{2},x_{2},x_{1};\mu)+\mathcal{O}(\alpha_{s})\,,$
(197)

where for $F_{\text{int}}$ the difference between (194) and (6.3.1) only
appears at $\mathcal{O}(\alpha_{s})$. We note that only the symmetric part of
$T_{a_{0}\to a_{1}a_{2}}^{ll^{\prime}}$ survives the contraction with
$\boldsymbol{y}^{l}\mskip 1.5mu\boldsymbol{y}^{l^{\prime}}$. As a consequence,
combinations of $a_{1},a_{2}$ with exactly one longitudinal polarisation give
zero, since one then has $T^{ll^{\prime}}\propto\epsilon^{ll^{\prime}}$ (see
section 7.4).

Little is known in phenomenology about twist-four distributions. The
computation using models of light-cone wave functions in Braun:2011aw  can
give guidance at large but not at small $x$. If one has a model for
${F}_{\text{int}}(x_{i},\boldsymbol{y})$ and is content with leading-order
accuracy in $\alpha_{s}$, one may combine (194) and (6.3.1) to approximate

$\displaystyle{}^{R}\mskip-1.5muF_{\text{int},\,a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu,\mu,\zeta)$
$\displaystyle={}^{R}\mskip-1.5muF_{\text{int},\,a_{1}a_{2}}(x_{i},\boldsymbol{y};\mu,\mu,\zeta)+\mathcal{O}(\alpha_{s})$
(198)

and use the intrinsic part of the DPDF on the r.h.s. in the
small-$\boldsymbol{y}$ contribution to the cross section, without further
approximation. The scales $\mu$ and $\zeta$ should of course be chosen
appropriately.

####  6.3.2 Evolved DPDs and cross section

Inserting the decomposition (192) of the DTMDs into (5.1) and neglecting
twist-three distributions, we obtain

$\displaystyle{}^{R}\mskip-1.5muF_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i},\zeta)$
$\displaystyle\quad=\exp\,\biggl{\\{}\int_{\mu_{0}}^{\mu_{1}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{1}}(\mu)-\gamma_{K,a_{1}}(\mu)\,\log\frac{\sqrt{x_{1}\zeta/x_{2}}}{\mu}\biggr{]}+{}^{1}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1};\mu_{0})\,\log\frac{\sqrt{\zeta}}{\sqrt{\zeta_{0}}}$
$\displaystyle\qquad\hskip
15.00002pt+\int_{\mu_{0}}^{\mu_{2}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{2}}(\mu)-\gamma_{K,a_{2}}(\mu)\,\log\frac{\sqrt{x_{2}\mskip
1.5mu\zeta/x_{1}}}{\mu}\biggr{]}+{}^{1}\mskip-1.5muK_{a_{2}}(\boldsymbol{z}_{2};\mu_{0})\,\log\frac{\sqrt{\zeta}}{\sqrt{\zeta_{0}}}\,\biggr{\\}}$
$\displaystyle\qquad\times\sum_{R^{\prime}}{}^{RR^{\prime}}\\!\exp\,\biggl{[}M_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})\mskip
1.5mu\log\frac{\sqrt{\zeta}}{\sqrt{\zeta_{0}}}\,\biggr{]}\,$
$\displaystyle\qquad\quad\times\Bigl{[}{}^{R^{\prime}}\\!F_{a_{1}a_{2},\,\text{spl}}({x}_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{0},\mu_{0},\zeta_{0})+{}^{R^{\prime}}\\!F_{a_{1}a_{2},\,\text{int}}({x}_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{0},\mu_{0},\zeta_{0})\Bigr{]}$
(199)

as our master formula for DTMDs at small $\boldsymbol{z}_{i}$ and small
$\boldsymbol{y}$. For the scales in the DTMDs we take

$\displaystyle\mu_{0}$ $\displaystyle=\mu_{Z}\,,$ $\displaystyle\zeta_{0}$
$\displaystyle=\mu_{Z}^{2}$ (200)

with $\mu_{Z}$ defined as in (181). We can then use the short-distance
approximations (193) and (194) for $F_{\text{spl}}$ and $F_{\text{int}}$. The
choice (200) will turn out to be well suited for combining the different
contributions to the overall cross section. In contrast to (6.2.2), we have
taken the two renormalisation scales $\mu_{01}$ and $\mu_{02}$ equal in
(6.3.2), given that we do not have an optimised short-distance expansion of
${F}$ for different scales. One can however optimise the resummation of
logarithms in the exponent by writing

$\displaystyle{}^{1}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1};\mu_{0})$
$\displaystyle={}^{1}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1};\mu_{01})+\int_{\mu_{0}}^{\mu_{01}}\frac{d\mu}{\mu}\,\gamma_{K,a_{1}}(\mu)$
(201)

and likewise for ${}^{1}\mskip-1.5muK_{a_{2}}$. A natural choice for the
additional scales is $\mu_{0i}=\mu_{z_{i}}$, as already suggested for (6.2.2).
Inserting the evolved DTMDs into (5.1), we have

$\displaystyle W_{\text{small $y$}}$
$\displaystyle=\exp\,\biggl{\\{}\int_{\mu_{0}}^{\mu_{1}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{1}}(\mu)-\gamma_{K,a_{1}}(\mu)\,\log\frac{{Q_{1}^{2}}}{\mu^{2}}\biggr{]}+{}^{1}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1},\mu_{0})\,\log\frac{Q_{1}\mskip
1.5muQ_{2}}{\zeta_{0}}$
$\displaystyle\qquad\;\;+\int_{\mu_{0}}^{\mu_{2}}\frac{d\mu}{\mu}\,\biggl{[}\gamma_{a_{2}}(\mu)-\gamma_{K,a_{2}}(\mu)\,\log\frac{{Q_{2}^{2}}}{\mu^{2}}\biggr{]}+{}^{1}\mskip-1.5muK_{a_{2}}(\boldsymbol{z}_{2},\mu_{0})\,\log\frac{Q_{1}\mskip
1.5muQ_{2}}{\zeta_{0}}\biggr{\\}}$
$\displaystyle\quad\times\sum_{RR^{\prime}}\Bigl{[}{}^{R}\mskip-1.5muF_{b_{1}b_{2},\,\text{spl}}(\bar{x}_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{0},\mu_{0},\zeta_{0})+{}^{R}\mskip-1.5muF_{b_{1}b_{2},\,\text{int}}(\bar{x}_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{0},\mu_{0},\zeta_{0})\Bigr{]}$
$\displaystyle\qquad\times\Phi(\nu\boldsymbol{y}_{+})\mskip
1.5mu\Phi(\nu\boldsymbol{y}_{-})\,\eta_{a_{1}a_{2}}(R)\;{}^{RR^{\prime}}\\!\exp\,\biggl{[}M_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})\log\frac{Q_{1}\mskip
1.5muQ_{2}}{\zeta_{0}}\,\biggr{]}\,$
$\displaystyle\qquad\times\Bigl{[}{}^{R^{\prime}}\\!F_{a_{1}a_{2},\,\text{spl}}({x}_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{0},\mu_{0},\zeta_{0})+{}^{R^{\prime}}\\!F_{a_{1}a_{2},\,\text{int}}({x}_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{0},\mu_{0},\zeta_{0})\Bigr{]}$
(202)

at the cross section level.

At present, one can only use the leading order (LO) expressions (193) and
(194) of the DTMD expansions, because the next-to-leading order (NLO) terms
have not been computed for either case. The leading $\mathcal{O}(\alpha_{s})$
term of $F_{\text{spl}}$ depends on $\boldsymbol{y}_{\pm}$ but not on
$\boldsymbol{Z}$, and the leading $\mathcal{O}(\alpha_{s}^{0})$ expression for
$F_{\text{int}}$ is independent of any transverse distance. A $\boldsymbol{Z}$
dependence via $\mu_{0}$ and $\zeta_{0}$ cancels in $W_{\text{small $y$}}$, up
to orders in $\alpha_{s}$ beyond the accuracy of the calculation. However,
$W_{\text{small $y$}}$ acquires a dependence on $\boldsymbol{Z}$ through the
$\boldsymbol{z}_{1}$ and $\boldsymbol{z}_{2}$ dependence of the Collins-Soper
kernels ${}^{1}\mskip-1.5muK_{a_{1}}(\boldsymbol{z}_{1};\mu_{0})$,
${}^{1}\mskip-1.5muK_{a_{2}}(\boldsymbol{z}_{2};\mu_{0})$ and
$M_{a_{1}a_{2}}(\boldsymbol{z}_{i},\boldsymbol{y})$. Since these kernels start
at $\mathcal{O}(\alpha_{s})$, the corresponding $\boldsymbol{Z}$ dependence
comes with a large logarithm $\log(Q_{1}Q_{2}/\zeta_{0})$ for each power of
$\alpha_{s}$ and is thus enhanced compared to the $\alpha_{s}$ corrections in
$F_{\text{spl}}$ and $F_{\text{int}}$, so that it is consistent to keep the
former while neglecting the latter. In this way, rapidity evolution provides a
nontrivial $\boldsymbol{Z}$ dependence and thus a nontrivial dependence on the
conjugate momentum $\boldsymbol{q}_{1}+\boldsymbol{q}_{2}$ in the cross
section, even if one uses the short-distance expansions (193) and (194) at
leading order. In physical terms, this is because the exponential in (6.3.2)
resums graphs for the soft factor that have a real gluon emission enhanced by
large rapidity logarithms, even though the graphs giving the leading-order
expansions of $F_{\text{spl}}$ and $F_{\text{int}}$ have no real gluon
emission into the final state.

###  6.4 Combining large and small $\boldsymbol{y}$

As we have seen in the previous subsections, the different regions of
$\boldsymbol{z}_{1},\boldsymbol{z}_{2}$ and $\boldsymbol{y}$ contributing to
the cross section involve different approximations, whose results are given in
(6.2.2) and (6.3.2). An important point is that the approximated expressions
give leading-power contributions not only in the regions for which they are
designed, but also outside these regions. With the short-distance behaviour
$F(x_{i},\boldsymbol{y})\sim 1/\boldsymbol{y}^{2}$ following from (6.3.1) we
see for instance that the integral of $W_{\text{large $y$}}$ over
$\boldsymbol{y}$ extends down to the smallest values allowed by the cutoff
function $\Phi(\nu\boldsymbol{y})$, where the approximations used to obtain
the expression are clearly invalid.

To deal with this problem we adapt the subtraction formalism discussed in
chapter 10 of Collins:2011zzd , which we briefly sketch now. The formalism is
formulated for a given Feynman graph $\Gamma$ in momentum space. It starts
with the smallest regions of loop momenta and works its way towards
increasingly larger regions. In this context, a region $r^{\prime}$ is called
smaller than $r$ (i.e. $r^{\prime}<r$) if hard momenta in $r$ are collinear or
soft in $r^{\prime}$, or if collinear momenta in $r^{\prime}$ are soft in $r$.
A more general definition of the relation $r^{\prime}<r$ is discussed in
chapters 5.4.1 and 10.1.3 of Collins:2011zzd . For each region $r$ that gives
a leading contribution to the cross section, a set of approximations $T_{r}$
is applied to the graph. In each approximated term $T_{r}\mskip 1.5mu\Gamma$
one still integrates over _all_ loop momenta, thus avoiding momentum cutoffs.
This leads to double counting of contributions from the overlap between
different regions, which are recursively removed by subtraction terms. The
full graph is then approximated by

$\displaystyle\Gamma$ $\displaystyle\approx\sum_{r}C_{r}\mskip 1.5mu\Gamma$
with $\displaystyle C_{r}\mskip 1.5mu\Gamma$ $\displaystyle=T_{r}\mskip
1.5mu\Gamma-\sum_{r^{\prime}<r}T_{r}\mskip 1.5muC_{r^{\prime}}\mskip
1.5mu\Gamma\,.$ (203)

$C_{r}\mskip 1.5mu\Gamma$ provides a valid approximation of $\Gamma$ in the
region $r$ and in all smaller regions $r^{\prime}<r$, thanks to the
subtraction terms $T_{r}\mskip 1.5muC_{r^{\prime}}\mskip 1.5mu\Gamma$, which
are obtained by applying the approximations for both $r$ and the smaller
regions. In the case where one has only two regions $r^{\prime}<r$, the
subtraction term is simply $T_{r}\mskip 1.5muT_{r^{\prime}}\mskip
1.5mu\Gamma$. In Collins:2011zzd  it was shown that if the approximation
$T_{r}$ is good up to power corrections in $\Lambda/Q$ in its design region
$r$, then $\sum_{r}C_{r}\mskip 1.5mu\Gamma$ is good with the same accuracy for
the full graph.

In Diehl:2017kgu , this formalism was applied to the double counting problem
between single and double parton scattering, leading to the formulation we
already presented in section 5.3. At variance to the original formulation in
Collins:2011zzd , momentum regions in (203) were replaced by regions of the
transverse distances $\boldsymbol{y}_{+}$ and $\boldsymbol{y}_{-}$, starting
from the region of large distances and going towards small ones in the
recursive construction of subtraction terms.

We now use the same procedure to remove the double counting between the
regions of “large” and “small” $\boldsymbol{y}$ as characterised in table 1.
The treatment of the regions where $|\boldsymbol{y}_{+}|$ or
$|\boldsymbol{y}_{-}|$ are of order $1/Q$ is postponed to the next subsection.
It is always understood that $|\boldsymbol{z}_{1}|$ and $|\boldsymbol{z}_{2}|$
are in the perturbative region, of order $1/q_{T}$, as discussed in section
6.1.

region | power counting | approximations
---|---|---
DPS, large $\boldsymbol{y}$ | $|\boldsymbol{y}|,|\boldsymbol{y}_{+}|,|\boldsymbol{y}_{-}|\sim 1/\Lambda$ | $|\boldsymbol{z}_{i}|\ll|\boldsymbol{y}|,1/\Lambda$
DPS, small $\boldsymbol{y}$ | $|\boldsymbol{y}|,|\boldsymbol{y}_{+}|,|\boldsymbol{y}_{-}|\sim 1/q_{T}$ | $|\boldsymbol{z}_{i}|,|\boldsymbol{y}|\ll 1/\Lambda$
SPS | $|\boldsymbol{y}_{+}|,|\boldsymbol{y}_{-}|\sim 1/Q$ | $|\boldsymbol{y}_{+}|,|\boldsymbol{y}_{-}|\ll|\boldsymbol{z}_{i}|\ll 1/\Lambda$
Table 1:  Regions of $\boldsymbol{y}$ discussed in the text. For power
counting, $|\boldsymbol{z}_{i}|\sim 1/q_{T}$ is always assumed.

According to the general construction, the subtraction term for large
$\boldsymbol{y}$ in the small-$\boldsymbol{y}$ term is obtained by applying
the approximations in the first and second row of table 1, i.e. by using both
$|\boldsymbol{z}_{i}|\ll|\boldsymbol{y}|$ and $|\boldsymbol{y}|\ll 1/\Lambda$.
For a given Feynman graph (and thus for any finite sum over graphs) the order
in which these approximations are made does not matter. This suggests that we
can start from either $W_{\text{large $y$}}$ or $W_{\text{small $y$}}$ and
then apply the approximation for the other region. At this point, we must
however note that the expressions of $W$ are not for fixed-order graphs but
contain logarithms resummed to all orders in the strong coupling. This is
reflected in the scale dependence of $\alpha_{s}$ and of the nonperturbative
functions. We must therefore make sure that the choice of scales in the
different terms is compatible with the way in which the subtraction formalism
works. In $W_{\text{small $y$}}$, this concerns the scales $\mu_{0}$ and
$\zeta_{0}$ at which one evaluates the short-distance expansions (193) and
(194) of DTMDs at fixed order in $\alpha_{s}$.

We now define the double counting subtraction term as

$\displaystyle W_{\text{sub}}$ $\displaystyle=W_{\text{small
$y$}}\,\bigl{|}_{\text{approx.\leavevmode\nobreak\ for}\leavevmode\nobreak\
|\boldsymbol{z}_{i}|\ll|\boldsymbol{y}|}$ (204)

with the small-$\boldsymbol{y}$ expression for $W$ given in (6.3.2). We should
take the limit $|\boldsymbol{z}_{1}|,|\boldsymbol{z}_{2}|\ll|\boldsymbol{y}|$
in all parts of that expression. In particular, this means replacing
$\boldsymbol{y}_{\pm}$ with $\boldsymbol{y}$ in $\Phi$ and in the leading-
order term of the short-distance expansion (193) of $F_{\text{spl}}$, and
replacing ${}^{RR^{\prime}}\\!M_{a_{1}a_{2}}$ with
$\delta_{RR^{\prime}}\,\bigl{[}{}^{R}\mskip-1.5muK_{a_{1}}+{}^{R}\mskip-1.5muK_{a_{2}}-{}^{1}\mskip-1.5muK_{a_{1}}-{}^{1}\mskip-1.5muK_{a_{2}}+{}^{R}\\!J\bigr{]}$
according to (133) and (171).

We note at this point that the integral of (6.3.2) over $\boldsymbol{y}$ is
actually divergent at large $\boldsymbol{y}$. Up to logarithmic corrections,
the 2v2 term $F_{\text{int}}\times F_{\text{int}}$ is constant and the 2v1
term $F_{\text{spl}}\times F_{\text{int}}$ goes like $1/\boldsymbol{y}^{2}$ in
that limit, and for $R=R^{\prime}=1$ there is no large-$\boldsymbol{y}$
suppression from the exponential since ${}^{1}\mskip-1.5muJ=0$. The
construction (204) ensures that the region $|\boldsymbol{y}|\gg 1/q_{T}$
cancels in $W_{\text{small $\boldsymbol{y}$}}-W_{\text{sub}}$. When evaluating
the cross section numerically, one may want to use an upper cutoff of order
$1/\Lambda$ on the integral over $|\boldsymbol{y|}$, thus avoiding to evaluate
terms in a region where they cancel anyway.

For the subtraction formalism to work, we need that

$\displaystyle W_{\text{sub}}$ $\displaystyle\approx W_{\text{large
$\boldsymbol{y}$}}$ for
$|\boldsymbol{y}|\sim|\boldsymbol{z}_{1}|,|\boldsymbol{z}_{2}|$ (205)

up to terms in $\alpha_{s}$ that are beyond the accuracy of the calculation.
In the specified region of $\boldsymbol{y}$, we can use the short-distance
limit of ${J}(\boldsymbol{y})$ and of the DPDFs in (6.2.2), the latter being
obtained by adding $F_{\text{spl}}$ and $F_{\text{int}}$ from (6.3.1). If one
retains only the leading-order terms for $F_{\text{spl}}+F_{\text{int}}$ on
both sides of (205), one finds manifest equality of the two expressions up to
differences due to scale dependence, which we will discuss shortly.

The full DPS cross section is then obtained from

$\displaystyle W_{\text{DPS}}(\nu)$ $\displaystyle=W_{\text{large
$y$}}(\nu^{\prime})-W_{\text{sub}}(\nu^{\prime})+W_{\text{small $y$}}(\nu)\,,$
(206)

where we make explicit the choice of cutoff parameters for the
$\boldsymbol{y}$ integration, taking $\nu^{\prime}\sim q_{T}$ and $\nu\sim Q$.
In the region $|\boldsymbol{y}|\gg|\boldsymbol{z}_{1}|,|\boldsymbol{z}_{2}|$
the last two terms cancel by virtue of (204), and one is left with the first
term, which is designed to give a correct approximation of the cross section
there. For $|\boldsymbol{y}|\sim|\boldsymbol{z}_{1}|,|\boldsymbol{z}_{2}|$,
the first and second terms cancel according to (205), and the third term gives
a correct approximation of the cross section. In this way, $W_{\text{DPS}}$
leads to a correct approximation of the cross section for $|\boldsymbol{y}|$
of order $1/q_{T}$ and larger.

We can take a more restrictive cutoff parameter $\nu^{\prime}\sim q_{T}$ in
the first two terms, because according to (205) the sum of these terms is
already suppressed for $|\boldsymbol{y}|\sim 1/q_{T}$, and the dependence on
the exact value of $\nu^{\prime}$ is thus beyond the accuracy of the
calculation. Taking $\nu^{\prime}$ lower than $Q$ is also more economical for
numerical calculations since it excludes a $\boldsymbol{y}$ region that is not
needed. By contrast, $\nu$ in the last term should be of order $Q$ (see the
next subsection).

Returning to the issue of scale dependence, we should first specify what is to
be used for the nonperturbative quantities in the expression (6.2.2) of
$W_{\text{large $y$}}$. For the Collins-Soper kernel we take the form in (183)
with $J(\boldsymbol{y}^{*};\mu_{y},\mu_{y})$ and $\gamma_{J}$ evaluated at
fixed perturbative order; this ensures the correct small-$\boldsymbol{y}$
limit. Using the DGLAP equations, the DPDFs are to be evolved from $\mu_{0i}$
to $\mu_{0}=\sqrt{\zeta_{0}}=\mu_{y}$, and it is at that scale that one makes
an ansatz which at small $\boldsymbol{y}$ tends to the fixed-order form
specified by (6.3.1).121212An explicit example for such an ansatz in the
colour singlet sector can be found in section 9.2 of Diehl:2017kgu . As usual,
this scale choice ensures that higher-order corrections are not accompanied by
large logarithms.

We now recall that in $W_{\text{small $y$}}$, and hence by construction also
in $W_{\text{sub}}$, we take $\mu_{0}=\sqrt{\zeta_{0}}=\mu_{Z}$ for the scales
at which we perform the short-distance expansion of the DTMDs. The difference
between the two sides in (205) is thus due to different scales of the short-
distance matching and to the different scales at which the DPDs and the
Collins-Soper kernels are evaluated in (6.2.2) and (6.3.2). However, all these
scales, $\mu_{y}$, $\mu_{z_{i}}$ and $\mu_{Z}$ are of the same order in the
region $|\boldsymbol{y}|\sim|\boldsymbol{z}_{i}|$ where we need (205). Using
the same methods as in section 6.2 of Diehl:2017kgu , one can show that the
approximation (205) holds up to higher-order terms in $\alpha_{s}$. These
terms are not accompanied by large logarithms and beyond the accuracy of the
calculation.

###  6.5 Combining DPS with SPS at short distances

In section 5.3 we recalled how to combine DPS, SPS and their interference in
TMD factorisation at generic values of $q_{T}$. We now derive the form of the
double counting subtraction terms for the case where $q_{T}\gg\Lambda$. We
restrict ourselves to the lowest order of $\alpha_{s}$ in the perturbative
short-distance coefficients. An extension to higher orders is beyond the scope
of the present work.

Let us start by discussing the relevant transverse distance scales and the
approximations they imply. We recall that the DPS cross section for
$q_{T}\gg\Lambda$ is dominated by distances
$|\boldsymbol{z}_{1}|,|\boldsymbol{z}_{2}|\ll 1/\Lambda$ thanks to the Fourier
exponentials in (5). In full analogy, the Fourier exponentials in (5.3) and
(5.3) ensure that $\sigma_{\text{DPS/SPS}}$ is dominated by
$|\boldsymbol{Z}|,|\boldsymbol{y}_{+}|\ll 1/\Lambda$ and
$\sigma_{\text{SPS/DPS}}$ by $|\boldsymbol{Z}|,|\boldsymbol{y}_{-}|\ll
1/\Lambda$. Likewise, $\sigma_{\text{SPS}}$ is dominated by
$|\boldsymbol{Z}|\ll 1/\Lambda$, where in the following we use the variable
$\boldsymbol{Z}$ instead of $\boldsymbol{z}$ in the SPS formula (1). According
to the formalism recalled in the previous subsection, the subtraction terms
for the overlap between double and single parton scattering are obtained by
applying the approximations appropriate for the SPS region to the DPS term in
the cross section. A corresponding prescription holds for the SPS/DPS
interference. The limits $\boldsymbol{y}_{+}\to\boldsymbol{0}$,
$\boldsymbol{y}_{-}\to\boldsymbol{0}$ or
$\boldsymbol{y}_{\pm}\to\boldsymbol{0}$ in the subtraction terms hence imply
that these distances are taken to be much smaller than $1/q_{T}$, because
$q_{T}$ is set to zero when one computes the hard scattering process for SPS.
We thus obtain the approximations given in table 2. The constraint on
$\boldsymbol{y}_{+}$ in the third row results from
$\boldsymbol{y}_{+}=\boldsymbol{y}_{-}+(\boldsymbol{z}_{1}-\boldsymbol{z}_{2})$
and the constraints on $\boldsymbol{y}_{-}$, $\boldsymbol{z}_{1}$ and
$\boldsymbol{z}_{2}$. We thus find that all terms in the overall cross section
(5.3) are dominated by distances in the perturbative region, except for
$\sigma_{\text{DPS}}$.

partial cross section | approximations
---|---
$\sigma_{\text{DPS/SPS}}$ | $|\boldsymbol{y}_{+}|\ll 1/\Lambda$
$\sigma_{\text{DPS/SPS},\,y_{+}\to 0}$ | $|\boldsymbol{y}_{+}|\ll 1/q_{T}$
$\sigma_{\text{DPS},\,y_{-}\to 0}$ |  $|\boldsymbol{y}_{-}|\ll 1/q_{T}$, $|\boldsymbol{y}_{+}|\ll 1/\Lambda$
$\sigma_{\text{DPS},\,y_{\pm}\to 0}$ | $|\boldsymbol{y}_{-}|,|\boldsymbol{y}_{+}|\ll 1/q_{T}$
Table 2:  Conditions on $\boldsymbol{y}_{\pm}$ for different terms in the
overall cross section (5.3) if $q_{T}\gg\Lambda$. Analogous relations hold for
$\sigma_{\text{SPS/DPS}}$ and the associated subtraction terms. In all cases
one has $|\boldsymbol{Z}|\ll 1/\Lambda$.

In all terms that are dominated by short distances, we can use the
perturbative splitting expressions (152) and (153) of the twist-three TMDs and
the DTMDs. Notice that the LO terms in these expressions have a factorised
dependence on $\boldsymbol{y}_{+}$, $\boldsymbol{y}_{-}$ and $\boldsymbol{Z}$,
so that they do not depend on the relative size of these distances (this will
no longer hold at NLO). Moreover, we can take the short-distance limit of the
TMDs $f_{a}(x,\boldsymbol{Z};\mu,\zeta)$. At LO accuracy we can simply replace
them with collinear PDFs $f_{a}(x;\mu)$, provided that $\mu^{2}\sim\zeta\sim
1/\boldsymbol{Z}^{2}$. Collinear PDFs are then the only nonperturbative input
necessary to compute the cross section, apart of course from the DPDs and
Collins-Soper kernels $J$ in $\sigma_{\text{DPS}}$. In a complete treatment,
one would also have collinear twist-three distributions, but we neglect these
as explained in section 5.3.

From the dependence on $\boldsymbol{y}_{\pm}$ of the splitting expressions
(152) and (153) and from the cross section formulae (1), (5.3) and (5.3) we
deduce that the SPS cross section, the SPS/DPS interference and all associated
subtraction terms have the same power behaviour as the 1v1 term of the DPS
cross section in (195), namely
$Q^{4}\,d\sigma/(d^{2}\boldsymbol{q}_{1}\,d^{2}\boldsymbol{q}_{2})\sim
1/q_{T}^{2}$.

We must now take a closer look at the scale dependence of the distributions.
Let us denote by $\mu_{h}$ the hard scale chosen in the computation of the SPS
amplitude, and by $\mu_{1}$ and $\mu_{2}$ the scales for the two DPS
amplitudes. To simplify the presentation, we first assume that all hard scales
are equal, $\mu=\mu_{h}=\mu_{1}=\mu_{2}$. The adaptation to the case of
different hard scales is discussed in appendix C. As we have explicitly shown
for single and double parton TMDs, the evolution in $\mu$ and $\zeta$ of these
distributions is multiplicative:

$\displaystyle f_{a}(x,\boldsymbol{z};\mu,\zeta)$
$\displaystyle=E_{2;a}(\boldsymbol{z};\mu,\zeta;\mu_{0},\zeta_{0})\,f_{a}(x,\boldsymbol{z};\mu_{0},\zeta_{0})\,,$
$\displaystyle
F_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu,\zeta)$
$\displaystyle=E_{4;a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu,\zeta;\mu_{0},\zeta_{0})\,F_{a_{1}a_{2}}(x_{i},\boldsymbol{z}_{i},\boldsymbol{y};\mu_{0},\zeta_{0})\,,$
(207)

where it is understood that $F$ is a vector and the evolution factor $E_{4}$ a
matrix in the space of colour representations $R$. From the discussion at the
end of section 5.3 if follows that for the twist-three TMDs, we have

$\displaystyle
D_{\alpha_{1}\alpha_{2}|\alpha_{0}}(x_{i},\boldsymbol{y}_{+},\boldsymbol{Z};\mu,\zeta)$
$\displaystyle=E_{3;\alpha_{1}\alpha_{2}|\alpha_{0}}(x_{i},\boldsymbol{y}_{+},\boldsymbol{Z};\mu,\zeta;\mu_{0},\zeta_{0})\,D_{\alpha_{1}\alpha_{2}|\alpha_{0}}(x_{i},\boldsymbol{y}_{+},\boldsymbol{Z};\mu_{0},\zeta_{0})\rule[-7.74998pt]{0.0pt}{0.0pt}$
(208)

and an analogous equation for $D_{\alpha_{0}|\alpha_{1}\alpha_{2}}$.

Let us now specify the evolution and scale dependence of the distributions
appearing in $\sigma_{\text{SPS}}$, the interference terms and
$\sigma_{\text{DPS}}$ in the region where all transverse distances are small
compared with $1/\Lambda$. In a schematic notation, we have

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

where $f(\mu_{Z})$ on the r.h.s. is a collinear PDF. We have specified only
the splitting part $F_{\text{spl}}$ of $F$, because the detailed form of
$F_{\text{int}}$ will not be needed. For brevity we have omitted parton
labels, as well as momentum fraction and transverse position arguments on the
r.h.s. The factors $\boldsymbol{y}_{+}^{l}/\boldsymbol{y}_{+}^{2}$ and
$\boldsymbol{y}_{-}^{l^{\prime}}/\boldsymbol{y}_{-}^{2}$ which accompany the
splitting kernels $U(\mu_{Z})$ and $U^{*}(\mu_{Z})$ have been omitted as well.
For brevity, we display neither the rapidity parameter $\zeta$ in the evolved
distributions nor the starting value $\zeta_{0}$. For the latter, we follow
(200) and take $\zeta_{0}=\mu_{Z}^{2}$. Since $\zeta$ is defined by (125), it
must be rescaled by $(x_{1}+x_{2})^{2}\mskip 1.5mu\zeta/(x_{1}x_{2})$ when
used in $f(x_{1}+x_{2},\boldsymbol{Z})$ and in $E_{2}$. The same forms as in
(6.5) are to be taken for the distributions in the left-moving proton, with
$x_{i}$ replaced by $\bar{x}_{i}$ and $\zeta$ by $\bar{\zeta}$.

The subtraction terms should be such that, up to higher orders in $\alpha_{s}$
and up to power corrections, the overall cross section (5.3) is given by

$\displaystyle d\sigma_{\text{DPS}}$ $\displaystyle\text{ \,for\,
}|\boldsymbol{y}_{+}|\sim|\boldsymbol{y}_{-}|\gg 1/\nu\,,$ $\displaystyle
d\sigma_{\text{SPS}}$ $\displaystyle\text{ \,for\,
}|\boldsymbol{y}_{+}|\sim|\boldsymbol{y}_{-}|\sim 1/\nu\,,$ $\displaystyle
d\sigma_{\text{DPS/SPS}}$ $\displaystyle\text{ \,for\,
}|\boldsymbol{y}_{+}|\gg|\boldsymbol{y}_{-}|\sim 1/\nu\,,$ $\displaystyle
d\sigma_{\text{SPS/DPS}}$ $\displaystyle\text{ \,for\,
}|\boldsymbol{y}_{-}|\gg|\boldsymbol{y}_{+}|\sim 1/\nu\,,$ (210)

where it is understood that $\nu\sim Q$. In the subtraction terms,
renormalisation and rapidity scales should be chosen such that the
cancellations required to avoid double counting are not spoiled by large
logarithms at higher orders.

To achieve this, we take so-called profile scales, which vary smoothly between
$\mu_{Z}$ and $\mu$ as a function of the appropriate transverse distance. We
therefore introduce a profile function $p$ that satisfies

$\displaystyle p({u};\mu_{a},\mu_{b})$ $\displaystyle\approx\mu_{a}$ for
${u}\sim 1$, $\displaystyle p({u};\mu_{a},\mu_{b})$
$\displaystyle\approx\mu_{b}$ for ${u}\gg 1$; (211)

a concrete example is given below. For the subtraction terms between DPS and
the SPS/DPS interference, we now take

$\displaystyle\sigma_{\text{DPS/SPS},y_{+}\to 0}$ $\displaystyle:D\sim
E_{3}(\mu;\hat{\mu})\,U(\hat{\mu})\,E_{2}(\hat{\mu};\mu_{Z})\,f(\mu_{Z})$
$\displaystyle\mbox{with }\hat{\mu}=p\bigl{(}\nu\mskip
1.5mu|\boldsymbol{y}_{+}|;\mu_{Z},\mu\bigr{)}\,,$
$\displaystyle\sigma_{\text{SPS/DPS},y_{-}\to 0}$ $\displaystyle:D\sim
E_{3}(\mu;\hat{\mu})\,U^{*}(\hat{\mu})\,E_{2}(\hat{\mu};\mu_{Z})\,f(\mu_{Z})$
$\displaystyle\mbox{with }\hat{\mu}=p\bigl{(}\nu\mskip
1.5mu|\boldsymbol{y}_{-}|;\mu_{Z},\mu\bigr{)}\,.$ (212)

This corresponds to a two-step matching: at the scale $\mu_{Z}$ a twist-two
TMD is matched onto a PDF. The twist-two TMD is then evolved up to the scale
$\hat{\mu}$, where a twist-three TMD is matched onto the twist-two TMD using
the lowest order term in (152). The twist-three TMD is then evolved up to the
final scale $\mu$. It is understood that at the intermediate scale $\hat{\mu}$
the rapidity parameter is taken equal to

$\displaystyle\widehat{\zeta}$ $\displaystyle=p(u;\mu_{Z}^{2},\zeta)\,,$ (213)

where the first argument $u$ is the same as the first argument in $\hat{\mu}$.
We thus find that

$\displaystyle\sigma_{\text{DPS/SPS},y_{+}\to
0}\approx\begin{cases}\sigma_{\text{DPS/SPS}}&\mbox{for
$|\boldsymbol{y}_{+}|\sim 1/\nu$}\\\ \sigma_{\text{SPS}}&\mbox{for
$|\boldsymbol{y}_{+}|\gg 1/\nu$}\end{cases}$ (214)

where we have used that for distances $|\boldsymbol{y}_{+}|\gg 1/Q$ the hard
scattering amplitude for SPS can be written as

$\displaystyle H_{\alpha_{0}\beta_{0}}\big{|}_{y_{+}\gg 1/Q}$
$\displaystyle=\frac{\boldsymbol{y}_{+}^{k}}{\boldsymbol{y}_{+}^{2}}\,\frac{\boldsymbol{y}_{+}^{l}}{\boldsymbol{y}_{+}^{2}}\;U^{k}_{\alpha_{0}\to\alpha_{1}\alpha_{2}}\,U^{l}_{\beta_{0}\to\beta_{1}\beta_{2}}\,H_{\alpha_{1}\beta_{1}}\,H_{\alpha_{2}\beta_{2}}$
(215)

at leading order in $\alpha_{s}$. To make this relation explicit, one needs to
Fourier transform the appropriate transverse momentum in the SPS graph to
position space. For diagrams that have kinematic overlap between SPS and the
SPS/DPS interference, but not with DPS, such as the one in figure 5b, the
behaviour in (214) and in its analogue for $\sigma_{\text{SPS/DPS},y_{-}\to
0}$ solve the double counting problem completely. If there is also overlap
with the DPS region, as in the three graphs of figure 5, then we need the DPS
subtraction terms, where we choose scales as follows:

$\displaystyle\sigma_{\text{DPS},y_{-}\to 0}$ $\displaystyle:F\sim
E_{4}(\mu;\hat{\mu})\,U^{*}(\hat{\mu})\,E_{3}(\hat{\mu};\mu_{Z})\,U(\mu_{Z})\,f(\mu_{Z})\,,$
$\displaystyle\sigma_{\text{DPS},y_{+}\to 0}$ $\displaystyle:F\sim
E_{4}(\mu;\hat{\mu})\,U(\hat{\mu})\,E_{3}(\hat{\mu};\mu_{Z})\,U^{*}(\mu_{Z})\,f(\mu_{Z})\,,$
$\displaystyle\sigma_{\text{DPS},y_{\pm}\to 0}$ $\displaystyle:F\sim
E_{4}(\mu;\hat{\mu})\,U(\hat{\mu})\
U^{*}(\hat{\mu})\,E_{2}(\hat{\mu};\mu_{Z})\,f(\mu_{Z})$ (216) with
$\displaystyle\hat{\mu}$
$\displaystyle=p\bigl{(}\nu\min\\{|\boldsymbol{y}_{+}|,|\boldsymbol{y}_{-}|\\};\mu_{Z},\mu\bigr{)}\,.$
(217)

In all three terms we have a two-step matching, with a twist-three TMD at the
intermediate stage in the first two cases and with a twist-two TMD in the
third case. If both $|\boldsymbol{y}_{+}|$ and $|\boldsymbol{y}_{-}|$ are
large compared with $1/\nu$, we thus have

$\displaystyle\sigma_{\text{DPS},y_{-}\to 0}$
$\displaystyle\approx\sigma_{\text{DPS/SPS}}\,,$
$\displaystyle\sigma_{\text{DPS},y_{+}\to 0}$
$\displaystyle\approx\sigma_{\text{SPS/DPS}}\,,$
$\displaystyle\sigma_{\text{DPS},y_{\pm}\to 0}$
$\displaystyle\approx\sigma_{\text{SPS}}$ (218)

so that the sum in (5.3) leaves only $\sigma_{\text{DPS}}$, as it should. It
remains to discuss the regions where one or both of $|\boldsymbol{y}_{+}|$ and
$|\boldsymbol{y}_{-}|$ is of order $1/\nu$. In each of these regions,
$\sigma_{\text{DPS}}$ can be obtained from $W_{\text{small $y$}}$, because
$|\boldsymbol{y}|\ll 1/\Lambda$. Furthermore, the integral over
$W_{\text{small $y$}}$ in these regions is dominated by the 1v1 term, with the
2v1 and 2v2 terms being power suppressed by at least a relative factor
$\Lambda^{2}/(\nu q_{T})$. This follows from the dependence of
$F_{\text{spl}}$ and $F_{\text{int}}$ on $\boldsymbol{y}_{+}$ and
$\boldsymbol{y}_{-}$. Finally, one finds that all three terms in (6.5) are
approximately equal to the 1v1 term in $\sigma_{\text{DPS}}$, which is thus
cancelled in the combination (5.3). This leaves us with $\sigma_{\text{SPS}}$,
$\sigma_{\text{DPS/SPS}}$ or $\sigma_{\text{SPS/DPS}}$ as appropriate. Overall
we thus have achieved a correct description of the cross section in all
relevant regions of $\boldsymbol{y}_{+}$ and $\boldsymbol{y}_{-}$.

Profile scales have been used in other contexts, see e.g. Ligeti:2008ac ;
Abbate:2010xh . For our purpose here, a suitable profile $p$ for $\hat{\mu}$
is given in equation (6.33) of Diehl:2017kgu . An alternative choice is

$\displaystyle p(u;\mu_{a},\mu_{b})=\begin{cases}\mu_{a}&\mbox{ \,for\, }u\leq
u_{a}\\\ \bigl{[}1-\frac{1}{2}(1-\cos u^{\prime}\pi)\bigr{]}\mskip
1.5mu\mu_{a}+\frac{1}{2}(1-\cos u^{\prime}\pi)\,\mu_{b}&\mbox{ \,for\,
}u_{a}<u<u_{b}\\\ \mu_{b}&\mbox{ \,for\, }u\geq u_{b}\\\ \end{cases}$ (219)

with

$\displaystyle u^{\prime}=\frac{u-u_{a}}{u_{b}-u_{a}}\,.$ (220)

The transition points for the $u$ dependence should be such that $u_{a}\sim 1$
and $u_{b}\gg 1$. A natural choice is $u_{a}={\nu}/{\mu}$ and
$u_{b}={\nu}/{\mu_{Z}}$.

###  6.6 Perturbative accuracy

Equations (2), (6.2.2) and (6.3.2) express TMDs in terms of collinear
distributions and perturbative quantities, and they permit the resummation of
large logarithms in the cross section. In the following, we give a brief
overview of the order in $\alpha_{s}$ at which the different quantities are
available. We also take a closer look at the type of logarithms that are
resummed at all orders. To this end, we count all large scales
$\mu_{i},\zeta,\nu$ as order $Q$. Logarithms of $|\boldsymbol{z}_{i}|$ turn
into logarithms of $|\boldsymbol{q}_{i}|$ after Fourier transformation, and we
count $\mu_{i0}\sim 1/|\boldsymbol{z}_{i}|\sim q_{T}$.

An important point in this context is that a variety of schemes for TMD
factorisation and for the associated resummation are used in the literature.
Some care is needed when converting perturbative expressions from one scheme
to another. A systematic discussion of this issue is given in Collins:2017oxh
, along with a comprehensive list of higher-order results for the quark
channel in the scheme we use here (referred to as CSS2 in that paper).

#### SPS

The SPS cross section can be calculated within the single TMD formalism
outlined in section 2. To compute the spectrum in $\boldsymbol{q}_{1}$ and
$\boldsymbol{q}_{2}$, one needs the relevant parton-level cross sections, e.g.
for gauge boson pair production, as a function of the relative transverse
momentum $\boldsymbol{q}_{1}-\boldsymbol{q}_{2}$. Logarithms of the net
transverse momentum $\boldsymbol{q}_{1}+\boldsymbol{q}_{2}$ are resummed by
the exponential in (2).

Perturbative ingredients in the cross section are the anomalous dimensions
$\gamma_{K}$, $\gamma_{a}$, the Collins-Soper kernel ${}^{1}\mskip-1.5muK$,
the short-distance (or matching) coefficients ${}^{1}C$, and finally the
coefficient $C_{H}^{2}$ of the hard-scattering cross section, defined as
$\hat{\sigma}$ divided by its tree-level value.131313For definiteness we
include the colour-singlet labels on $K$ and $C$ here, although this is the
only colour channel that appears in SPS. In table 3 we specify different
levels of accuracy, following the scheme of Abbate:2010xh ; Berger:2010xi .
After expanding $\alpha_{s}(\mu^{\prime})$ in $\gamma_{a}(\mu^{\prime})$ and
$\gamma_{K}(\mu^{\prime})$ around $\alpha_{s}(\mu)$, the integral over
$\mu^{\prime}$ in (2) gives powers of $\log(\mu/\mu_{0})$. We thus find that
$\gamma_{K}$ goes with double logarithms $\alpha_{s}(Q)\log^{2}(Q/q_{T})$,
whilst $\gamma_{a}$ and ${}^{1}\mskip-1.5muK$ go with single logarithms
$\alpha_{s}(Q)\log(Q/q_{T})$. One therefore requires one order higher for
$\gamma_{K}$ than for $\gamma_{a}$ and ${}^{1}\mskip-1.5muK$ (see table 3).
Likewise, one finds that one needs the same perturbative order for the $\beta$
function as for $\gamma_{K}$. At accuracy NkLO one then has control over all
terms from $\alpha_{s}^{n}\log^{n+1}$ to $\alpha_{s}^{n+k}\log^{n+1}$ in the
exponent. $\gamma_{K}$, $\gamma_{a}$ and ${}^{1}\mskip-1.5muK$ all start at
$\mathcal{O}(\alpha_{s})$ and are known up to $\mathcal{O}(\alpha_{s}^{3})$,
both in the quark and in the gluon channel. A compilation of these results can
for instance be found in appendix D of Echevarria:2016scs . The notation in
that work is related to the one we are using by

$\displaystyle\gamma_{K,a}$ $\displaystyle=2\mskip
1.5mu\Gamma_{\text{cusp}}^{a}\,,$ $\displaystyle\gamma_{a}$
$\displaystyle={}-\gamma_{V}^{a}\,,$ $\displaystyle\gamma_{F,a}(\mu,\zeta)$
$\displaystyle=\gamma^{a}(\mu,\zeta)$ (221)

and

$\displaystyle{}^{1}\mskip-1.5muK_{a}(\boldsymbol{z},\mu)$
$\displaystyle=-2\mathcal{D}^{a}(\mu,\boldsymbol{b}_{T})$
$\displaystyle\text{with }\boldsymbol{z}=\boldsymbol{b}_{T}\,,$ (222)

where $a=q,g$ and it is understood that the strong coupling is taken at scale
$\mu$. Results for the quark channel are also given in Collins:2017oxh , where
the same notation as here is used (except that the Collins-Soper kernel is
denoted by $\tilde{K}$ instead of $K$).

| $\gamma_{K}$ | $\gamma_{a}$ |  ${}^{R}\mskip-1.5muK$, ${}^{RR^{\prime}}\\!M$ |  ${}^{R}C$, ${}^{RR^{\prime}}\mskip-1.5muC_{\text{int}}$ | $C_{H}^{2}\rule[-4.73611pt]{0.0pt}{4.30554pt}$ | ${}^{R}\mskip-1.5muP$
---|---|---|---|---|---|---
LL | $\alpha_{s}$ | $-$ | $-$ | $\alpha_{s}^{0}$ | $\alpha_{s}^{0}$ | $\alpha_{s}$
NLL | $\alpha_{s}^{2}$ | $\alpha_{s}$ | $\alpha_{s}$ | $\alpha_{s}^{0}$ | $\alpha_{s}^{0}$ | $\alpha_{s}$
NLL′ | $\alpha_{s}^{2}$ | $\alpha_{s}$ | $\alpha_{s}$ | $\alpha_{s}$ | $\alpha_{s}$ | $\alpha_{s}^{2}$
NNLL | $\alpha_{s}^{3}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{2}$ | $\alpha_{s}$ | $\alpha_{s}$ | $\alpha_{s}^{2}$
NNLL′ | $\alpha_{s}^{3}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{2}$ | $\alpha_{s}^{3}$
Table 3:  Different levels of accuracy and the associated orders of
$\alpha_{s}$ for the perturbative ingredients in the SPS cross section. The
order required for the $\beta$ function is the same as for $\gamma_{K}$.

At lowest order one has ${}^{1}C_{ab}=\delta_{ab}\,\delta(1-x)$ and
$C_{H}^{2}=1$. These coefficients are not associated with the resummation of
large logarithms, and conventionally one requires them at one order in
$\alpha_{s}$ lower than $\gamma_{a}$. This corresponds to counting
$\log(Q/q_{T})$ as order $1/\alpha_{s}$. Taking them at the same order as
$\gamma_{a}$ is indicated by a prime in the table. For unpolarised quarks or
gluons, ${}^{1}C$ has been computed up to $\mathcal{O}(\alpha_{s}^{2})$ in
Echevarria:2016scs  and Gehrmann:2014yya . With unpolarised protons this is
sufficient for $q\bar{q}$ annihilation processes, whereas for $gg$ initiated
processes one also requires the coefficient ${}^{1}C_{\delta g\mskip 1.5mug}$
for linear gluon polarisation in the TMDs. This is well known from the process
$gg\to H$ (see e.g. Echevarria:2015uaa  and references therein).

To the best of our knowledge, the hard-scattering coefficient $C_{H}^{2}$ has
not been given in the literature for final states relevant in our context,
such as $W^{+}W^{-}$, $ZH$, $HH$, etc. It should however be possible to
extract it from the virtual corrections to the relevant amplitudes, which are
known up to $\mathcal{O}(\alpha_{s}^{2})$ in some channels (for gauge boson
pair production see e.g. Gehrmann:2015ora ; Caola:2015ila  and references
therein). Without this additional effort, the SPS contribution to the cross
section can currently be evaluated at NLL accuracy.

There is one more type of large logarithms hidden in (2), namely the
logarithms of $q_{T}/\Lambda$ that are resummed by DGLAP evolution of the PDFs
on the r.h.s., from a typical starting scale of order $\Lambda$ up to
$\mu_{0}\sim 1/|\boldsymbol{z}|$. A customary choice is to take the DGLAP
kernels ${}^{1}P$ at one power in $\alpha_{s}$ higher than the matching
coefficients ${}^{1}C$, as shown in the table. This amounts to using LO
evolution together with LO matching coefficients, etc. The perturbative
accuracy of ${}^{1}C$ is then retained even if $\log(q_{T}/\Lambda)$ is as
large as $1/\alpha_{s}$. For unpolarised PDFs, the DGLAP kernels are known up
to $\mathcal{O}(\alpha_{s}^{3})$ Moch:2004pa ; Vogt:2004mw .

#### Interference between SPS and DPS

The perturbative ingredients required for computing the SPS/DPS interference
term, as laid out in section 6.5, are the Collins-Soper kernels
$K_{\alpha_{1}\alpha_{2}|\alpha_{0}}$ for twist-three TMDs, the splitting
kernels $U$ appearing in their short-distance approximation, and the
interference between SPS and DPS hard-scattering amplitudes for the process
under investigation. None of these have been calculated so far. It is easy to
compute the $\mathcal{O}(\alpha_{s})$ expressions of
$K_{\alpha_{1}\alpha_{2}|\alpha_{0}}$ and $U$ by adapting the corresponding
calculations for DPDs in Diehl:2011yj . Combining the known hard-scattering
amplitudes for the relevant SPS and DPS processes is a straightforward
exercise (but can be tedious, especially for heavy gauge bosons due to their
three helicities).

#### DPS: small-$\boldsymbol{y}$ contribution

Let us now discuss the perturbative ingredients in the small-$y$ expression
(6.3.2) for DPS. $\gamma_{K}$, $\gamma_{a}$ and ${}^{1}\mskip-1.5muK$ are the
same as in the SPS formula, and in addition one needs the matrix
${}^{RR^{\prime}}\\!M$, which starts at $\mathcal{O}(\alpha_{s})$ and goes
with single logarithms $\log(Q/q_{T})$. The expression for the matching
${}^{R}\mskip-1.5muF_{\text{spl}}$ on PDFs involves a coefficient
${}^{R}C_{\text{spl}}$, and the matching of ${}^{R}\mskip-1.5muF_{\text{int}}$
on twist-four distributions involves a coefficient
${}^{RR^{\prime}}\mskip-1.5muC_{\text{int}}$. The lowest order of
${C}_{\text{spl}}$ and ${C}_{\text{int}}$ is $\mathcal{O}(\alpha_{s})$ and
$\mathcal{O}(\alpha_{s}^{0})$, respectively — their explicit definitions are
not needed here. The hard-scattering subprocesses for DPS are simpler than for
SPS with the same final state. For Drell-Yan and for Higgs boson production
the hard-scattering coefficients $C_{H}^{2}$ are known up to
$\mathcal{O}(\alpha_{s}^{3})$. They can for instance be found in
Collins:2017oxh  for Drell-Yan production and in Ebert:2017uel  for both
channels.

The evolution kernel ${}^{RR^{\prime}}\\!M$ is known at
$\mathcal{O}(\alpha_{s}^{2})$ in the two-quark channel Vladimirov:2016qkd ;
its $\mathcal{O}(\alpha_{s})$ expressions for all channels are given in
section 7.2. The matching coefficients ${}^{R}C_{\text{spl}}$ have been
calculated at $\mathcal{O}(\alpha_{s})$ in Diehl:2011yj , and the related
kernels ${}^{R}T$ are compiled in section 7.4 here. For
${}^{RR^{\prime}}\mskip-1.5muC_{\text{int}}$, only the trivial order
$\alpha_{s}^{0}$ is currently known, which corresponds to matching
$F_{\text{int}}=G$. The kernels for the DGLAP evolution of the distributions
$G$ are known at $\mathcal{O}(\alpha_{s})$ and given in Bukhvostov:1985rn ;
for a more general discussion of twist-four evolution we refer to Braun:2009vc
. At a given level of accuracy, their order should be the same as the one of
${}^{R}\mskip-1.5muP$ in table 3.

When denoting the level of perturbative accuracy, one might think of combining
equal orders in $\alpha_{s}$ for $C_{\text{spl}}$ and $C_{\text{int}}$. On the
other hand, $C_{\text{spl}}$ starts one order higher in $\alpha_{s}$ than
$C_{\text{int}}$, and moreover the different combinations of $F_{\text{spl}}$
and $F_{\text{int}}$ in the cross section have different power behaviour in
$\Lambda^{2}/q_{T}^{2}$ according to (195). We thus find it more natural to
take one power in $\alpha_{s}$ more for $C_{\text{spl}}$ than for
$C_{\text{int}}$. With this naming convention, the small-$y$ expression of DPS
can currently be evaluated to NLL accuracy. Missing ingredients for achieving
NNLL in pure quark channels or NLL′ in channels with gluons are the matching
coefficients $C_{\text{spl}}$ and $C_{\text{int}}$ for DTMDs, as well as the
two-loop DGLAP kernels for twist-four evolution.

#### DPS: large-$\boldsymbol{y}$ contribution

We now turn our attention to the large-$y$ expression (6.2.2) of DPS. The
quantities $\gamma_{K}$, $\gamma_{a}$ and $C_{H}^{2}$ are the same as in the
small-$y$ expressions of DPS. The resummation of single logarithms
$\log(Q/q_{T})$ requires the Collins-Soper kernels ${}^{R}\mskip-1.5muK$ in
the different colour channels, which can be obtained from the perturbative
expression of ${}^{RR^{\prime}}\\!M(\boldsymbol{z}_{i},\boldsymbol{y})$ by
taking the limit $|\boldsymbol{z}_{i}|\ll|\boldsymbol{y}|$ and using that in
this limit one has the structure (171). In the two-quark channel one can thus
obtain the $\mathcal{O}(\alpha_{s}^{2})$ expression for ${}^{R}\mskip-1.5muK$
from the results of Vladimirov:2016qkd ; in all other channels one has the
one-loop expressions, which we list in section 7.2. For the matching
coefficients ${}^{R}C$ we will derive $\mathcal{O}(\alpha_{s})$ results in all
colour and polarisation channels in the next section; most of them have been
given in the literature before. One can thus evaluate gauge boson pair
production with NNLL accuracy, whereas for DPS processes with one or two Higgs
bosons in the final state, NLL′ accuracy is currently achievable.

Let us finally take a closer look at the evolution of the DPDFs, which
generates logarithms involving a hadronic starting scale $\Lambda$. The
situation is now more involved than the one discussed for SPS above. In
analogy to usual DGLAP evolution, there are single logarithms
$\log(q_{T}/\Lambda)$, whose resummation requires the $\zeta$ independent part
of the evolution kernels ${}^{R}\mskip-1.5muP$. We will give
$\mathcal{O}(\alpha_{s})$ expressions for these kernels in all polarisation
channels in section 7.3.5. The $\zeta$ dependence of the DPDFs, which has no
counterpart for ordinary parton distributions, is made explicit in (6.2.2).
Expanding $\alpha_{s}(\mu)$ in the anomalous dimension ${}^{R}\gamma_{J}(\mu)$
around a fixed scale, we can perform the integral and obtain double logarithms
$\log^{2}(Q/\Lambda)-\log^{2}(Q/q_{T})=\log(\Lambda
q_{T}/Q^{2})\,\log(\Lambda/q_{T})$ in the exponential. If both $\log(Q/q_{T})$
and $\log(q_{T}/\Lambda)$ are counted as order $1/\alpha_{s}$, one therefore
needs ${}^{R}\gamma_{J}$ at the same order as the cusp anomalous dimension
$\gamma_{K}$. For smaller values of $\log(q_{T}/\Lambda)$ one may instead be
content with ${}^{R}\gamma_{J}$ at the same order as $\gamma_{a}$ and
${}^{R}\mskip-1.5muK$. From the two-loop results in Vladimirov:2016qkd  one
can extract the $\mathcal{O}(\alpha_{s}^{2})$ term of
${}^{8}\mskip-1.5mu\gamma_{J}$, whereas for higher colour representations we
currently only know the $\mathcal{O}(\alpha_{s})$ results given in section
7.2. We note that in (6.2.2) there are also single logarithms $\log(Q/q_{T})$
multiplied with ${}^{R}\\!J(\boldsymbol{y})$, which at large $\boldsymbol{y}$
is beyond the control of perturbation theory.

##  7 One-loop results

In this section, we present the perturbative ingredients necessary for
evaluating DPS at NLL accuracy (see table 3). Some of the expressions are
quoted from the literature to make the presentation self-contained, whereas
for evolution kernels and matching coefficients we verify previous and derive
new results.

In section 7.1 we discuss the hard-scattering cross sections for Drell-Yan and
for Higgs production via gluon fusion, and in section 7.2 we give the one-loop
expressions of the Collins-Soper kernels ${}^{RR^{\prime}}\\!K_{a_{1}a_{2}}$
for all parton and colour channels. Section 7.3 is devoted to the
$\mathcal{O}(\alpha_{s})$ matching coefficients ${}^{R}C_{ab}$ of DTMDs in the
large-$\boldsymbol{y}$ region. After explaining a number of computational
details, we present results for all polarisation and colour combinations in
sections 7.3.4 and 7.3.5. As a corollary we obtain the one-loop DGLAP kernels
(7.3.5) for colour non-singlet DPDFs. In section 7.4 we list the
$\mathcal{O}(\alpha_{s})$ splitting kernels ${}^{R}T_{a_{0}\to a_{1}a_{2}}$ of
DTMDs in the small-$\boldsymbol{y}$ region, adapting the results of
Diehl:2011yj  to the notation used here.

###  7.1 TMD hard-scattering cross sections

We define the coefficients $C_{H}^{2}$ of the hard-scattering cross section as

$\displaystyle\hat{\sigma}(Q^{2},\mu^{2})=\hat{\sigma}_{0}(Q^{2})\,C_{H}^{2}(Q^{2},\mu^{2})$
(223)

for a given partonic channel, where $\hat{\sigma}_{0}$ is the cross section at
lowest order in $\alpha_{s}$. The expressions for $\hat{\sigma}_{0}$ can be
found in many places, for instance in Kasemets:2012pr  for the production of
an electroweak gauge boson $V$ and its subsequent decay into a lepton pair,
and in Echevarria:2015uaa  for Higgs production via $gg$ fusion in the limit
of large top mass. The hard-scattering coefficients read

$\displaystyle C_{H,\,q\bar{q}\to V}^{2}(Q^{2},\mu^{2})$
$\displaystyle=1+\frac{\alpha_{s}(\mu)}{2\pi}\,C_{F}\,\biggl{[}-\ln^{2}\frac{Q^{2}}{\mu^{2}}+3\ln\frac{Q^{2}}{\mu^{2}}-8+\frac{7\pi^{2}}{6}\mskip
1.5mu\biggr{]}\,,$ $\displaystyle C_{H,\,gg\to H}^{2}(Q^{2},\mu^{2})$
$\displaystyle=\biggl{(}1+\frac{\alpha_{s}(\mu)}{2\pi}\,C_{A}\,\biggr{[}-\ln^{2}\frac{Q^{2}}{\mu^{2}}+\frac{7\pi^{2}}{6}\mskip
1.5mu\biggr{]}\biggr{)}\,\biggl{(}1+\frac{\alpha_{s}(\mu)}{2\pi}\,\bigl{[}5C_{A}-3C_{F}\bigr{]}\biggr{)}$
(224)

and are known up to $\mathcal{O}(\alpha_{s}^{3})$, see for instance appendix A
of Ebert:2017uel . The colour factors are $C_{F}=(N^{2}-1)/(2N)$ and $C_{A}=N$
as usual. The first factor in $C_{H,\,gg\to H}^{2}$ is for two-gluon fusion
via the local current $F^{\mu\nu\mskip 1.5mua}F_{\mu\nu}^{a}$ and the second
factor for the coupling of this current to the Higgs boson via a top quark
loop. Throughout this section we drop the explicit indication of higher-order
corrections in equations, omitting a term $\\!{}+\mathcal{O}(\alpha_{s}^{2})$
in the first line of (7.1) and so forth.

An important property of the coefficients $C_{H}^{2}$ for Drell-Yan and Higgs
production is their independence of the incoming parton polarisations. This
holds at all orders in $\alpha_{s}$. For $gg\to H$ it follows from the fact
that angular momentum conservation only allows for two non-vanishing helicity
amplitudes, which are related by parity invariance. Therefore, all
$\alpha_{s}$ corrections can be expressed by a single hard-scattering
coefficient. The argument holds for the production of any scalar or
pseudoscalar boson via $gg$ fusion.

For $q\bar{q}$ annihilation into a boson via the vector current, one
additionally needs the conservation of quark chirality, which together with
angular momentum and parity conservation leaves a single independent helicity
amplitude. To appeal to quark chirality conservation in a calculation using
dimensional regularisation, one can adapt the argument given in Collins:1999un
. As long as one works in $D=4-2\epsilon$ dimensions, one can avoid specifying
the polarisation of the incoming partons and instead compute the hard
scattering as a matrix in Dirac space. Chirality conservation then technically
means that the hard-scattering amplitude contains an odd number of
$\gamma_{\mu}$ matrices at any perturbative order. The same holds for the
terms that must be subtracted to remove ultraviolet and infrared divergences.
After these subtractions have been done, one can take $D=4$ and use the
relation between quark helicity and chirality.

The extension of this argument to $q\bar{q}$ annihilation via the axial vector
current is simple in $D=4$ dimensions: anticommuting the $\gamma_{5}$ matrix
from the current, one finds $A_{\mu}=V_{\mu}\gamma_{5}$ for the Dirac matrix
of the hard-scattering amplitudes of the two currents (with their divergences
subtracted). Their independent helicity amplitudes are hence the same (or
opposite in sign). To complete the argument, one needs to discuss the
implementation of the electroweak axial current in dimensional regularisation,
which we shall not do here.

###  7.2 Collins-Soper kernels and anomalous dimensions

In the limit of small $\boldsymbol{z}_{1}$, $\boldsymbol{z}_{1}$ and
$\boldsymbol{y}$, the soft matrix for DPS can be directly computed in
perturbation theory, as we derived in (191). Setting $U=\mathbb{1}$ in (332),
we simply have $K=\widehat{K}=dS/dY$ at $\mathcal{O}(\alpha_{s})$. The one-
loop calculation of $S$ in the $qq$ channel is presented in detail in section
3.3.2 of Diehl:2011yj . We can rewrite the results there in such a way that we
have a separate dependence on the two scales $\mu_{1}$ and $\mu_{2}$.
Splitting off a diagonal part according to (133), we get

$\displaystyle K_{qq}$ $\displaystyle=C_{F}\mskip
1.5mu\bigl{[}K(\boldsymbol{z}_{1};\mu_{1})+K(\boldsymbol{z}_{2};\mu_{2})\bigr{]}\begin{pmatrix}1\leavevmode\nobreak\
&0\\\ 0\leavevmode\nobreak\ &1\end{pmatrix}+M_{qq}$ (225)

with

$\displaystyle M_{qq}$
$\displaystyle=\begin{pmatrix}0&\frac{\sqrt{N^{2}-1}}{2N}\,K_{d}\\\\[3.00003pt]
\frac{\sqrt{N^{2}-1}}{2N}\,K_{d}\leavevmode\nobreak\ &{}-\frac{N}{2}\mskip
1.5muK_{y}-\frac{1}{N}\mskip 1.5muK_{d}\,\end{pmatrix}\,,$ (226)

where the matrices are in colour representation space, with $R=(1,8)$. We
recall that $M_{qq}$ is renormalisation scale independent. The scale
independent combinations of kernels read

$\displaystyle K_{y}(\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=K(\boldsymbol{z}_{1};\mu)+K(\boldsymbol{z}_{2};\mu)-K\Bigl{(}\boldsymbol{y}+{\textstyle\frac{1}{2}}(\boldsymbol{z}_{1}+\boldsymbol{z}_{2});\mu\Bigr{)}-K\Bigl{(}\boldsymbol{y}-{\textstyle\frac{1}{2}}(\boldsymbol{z}_{1}+\boldsymbol{z}_{2});\mu\Bigr{)}\,,$
$\displaystyle K_{d}(\boldsymbol{z}_{i},\boldsymbol{y})$
$\displaystyle=K\Bigl{(}\boldsymbol{y}+{\textstyle\frac{1}{2}}(\boldsymbol{z}_{1}+\boldsymbol{z}_{2});\mu\Bigr{)}+K\Bigl{(}\boldsymbol{y}-{\textstyle\frac{1}{2}}(\boldsymbol{z}_{1}+\boldsymbol{z}_{2});\mu\Bigr{)}$
$\displaystyle\quad-K\Bigl{(}\boldsymbol{y}+{\textstyle\frac{1}{2}}(\boldsymbol{z}_{1}-\boldsymbol{z}_{2});\mu\Bigr{)}-K\Bigl{(}\boldsymbol{y}-{\textstyle\frac{1}{2}}(\boldsymbol{z}_{1}-\boldsymbol{z}_{2});\mu\Bigr{)}\,,$
(227)

where

$\displaystyle K(\boldsymbol{z};\mu)$
$\displaystyle=-\frac{\alpha_{s}(\mu)}{\pi}\,\log\frac{\boldsymbol{z}^{2}\mu^{2}}{b_{0}^{2}}$
(228)

is the one-loop Collins-Soper kernel for single parton TMDs, with the colour
factor removed. We recall that $b_{0}=2e^{-\gamma_{E}}$. We note that the two-
loop result in section 4 of Vladimirov:2016qkd  has the form given by (225),
(226) and (7.2) if one replaces $C_{F}\mskip 1.5muK(\boldsymbol{z};\mu)$ with
the two-loop Collins-Soper kernel for a single quark TMD.

In the solution (5.1) of the evolution equations, one needs the matrix
exponential of $M_{qq}$ times $L=\log\sqrt{\zeta/\zeta_{0}}$, which reads (see
section 3.4.2 of Diehl:2011yj )

$\displaystyle\exp[\mskip 1.5muLM_{qq}\mskip 1.5mu]$
$\displaystyle=\frac{1}{d_{+}-d_{-}}\begin{pmatrix}d_{+}e^{Ld_{-}}-d_{-}e^{Ld_{+}}\leavevmode\nobreak\
&\frac{\sqrt{N^{2}-1}}{2N}\,K_{d}\,\bigl{(}e^{Ld_{+}}-e^{Ld_{-}}\big{)}\,\\\\[3.00003pt]
\frac{\sqrt{N^{2}-1}}{2N}\,K_{d}\,\bigl{(}e^{Ld_{+}}-e^{Ld_{-}}\big{)}\leavevmode\nobreak\
&d_{+}e^{Ld_{+}}-d_{-}e^{Ld_{-}}\end{pmatrix}$ (229)

with

$\displaystyle d_{\pm}$
$\displaystyle=\frac{1}{2N}\Biggl{[}-\frac{N^{2}}{2}K_{y}-K_{d}\pm\sqrt{\biggl{(}\frac{N^{2}}{2}K_{y}+K_{d}\biggr{)}^{2}+(N^{2}-1)K_{d}^{2}}\;\Biggr{]}\,.$
(230)

An interesting situation arises if $|K_{d}|\ll\sqrt{N^{2}-1}\;|K_{y}|$, which
holds for instance in the large-$N$ limit. As shown in section 3.4.2 of
Diehl:2011yj , one then finds that ${}^{8}\mskip-1.5muF(\zeta)$ is suppressed
compared with ${}^{1}\mskip-1.5muF(\zeta)$ for $\zeta\gg\zeta_{0}$. Note that
the above condition always holds in the limit
$\boldsymbol{z}_{1}=\boldsymbol{z}_{2}=\boldsymbol{0}$, where $K_{d}$ tends to
zero.

It is easy to extend the above results to other parton channels, because the
one-loop graphs remain the same up to colour factors. For the quark-gluon
channels, we obtain

$\displaystyle K_{qg}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i})$
$\displaystyle=\bigl{[}\mskip 1.5muC_{F}\mskip
1.5muK(\boldsymbol{z}_{1};\mu_{1})+C_{A}\mskip
1.5muK(\boldsymbol{z}_{2};\mu_{2})\mskip 1.5mu\bigr{]}\mathbb{1}+M_{qg}$ (231)

with

$\displaystyle M_{qg}$
$\displaystyle=\begin{pmatrix}0&\frac{1}{\sqrt{2}}\mskip
1.5muK_{d}&0\\\\[3.99994pt] \frac{1}{\sqrt{2}}\mskip
1.5muK_{d}\leavevmode\nobreak\ &-\frac{N}{2}\mskip
1.5mu(K_{y}+{\textstyle\frac{1}{2}}K_{d})&\frac{\sqrt{N^{2}-4}}{4}\mskip
1.5muK_{d}\\\\[3.99994pt] 0&\frac{\sqrt{N^{2}-4}}{4}\mskip
1.5muK_{d}&-\frac{N}{2}\mskip
1.5mu(K_{y}+{\textstyle\frac{1}{2}}K_{d})\,\end{pmatrix}\,.$ (232)

Here the space of colour representations is $R=(1,A,S)$, corresponding to the
projectors in (4.1), and $\mathbb{1}$ is the unit matrix in that space. For
the two-gluon channel we get the kernel

$\displaystyle K_{gg}(\boldsymbol{z}_{i},\boldsymbol{y};\mu_{i})$
$\displaystyle=C_{A}\mskip
1.5mu\bigl{[}K(\boldsymbol{z}_{1};\mu_{1})+K(\boldsymbol{z}_{2};\mu_{2})\bigr{]}\mathbb{1}+M_{gg}$
(233)

with

$\displaystyle M_{gg}$
$\displaystyle=\begin{pmatrix}0&\frac{N}{\sqrt{N^{2}-1}}\mskip
1.5muK_{d}&0&0&0\\\\[3.00003pt] \frac{N}{\sqrt{N^{2}-1}}\mskip
1.5muK_{d}\leavevmode\nobreak\ \leavevmode\nobreak\ &-\frac{N}{2}\mskip
1.5mu(K_{y}+{\textstyle\frac{1}{2}}K_{d})&\frac{N}{4}\mskip
1.5muK_{d}&0&\sqrt{\frac{3}{8}}\mskip 1.5muK_{d}\\\\[3.00003pt]
0&\frac{N}{4}\mskip 1.5muK_{d}&-\frac{N}{2}\mskip
1.5mu(K_{y}+{\textstyle\frac{1}{2}}K_{d})&-\frac{3}{\sqrt{10}}\mskip
1.5muK_{d}&0\\\\[3.00003pt] 0&0&-\frac{3}{\sqrt{10}}\mskip
1.5muK_{d}&-3(K_{y}+{\textstyle\frac{1}{2}}K_{d})&-\sqrt{\frac{3}{5}}\mskip
1.5muK_{d}\\\\[3.00003pt] 0&\sqrt{\frac{3}{8}}\mskip
1.5muK_{d}&0&-\sqrt{\frac{3}{5}}\mskip
1.5muK_{d}&-4(K_{y}+{\textstyle\frac{1}{2}}K_{d})\,\end{pmatrix}$ (234)

in the space of colour representations $R=(1,A,S,D,27)$, corresponding to the