DESY 21-148

UUITP-46/21

Double parton distributions out of bounds in colour space

M. Diehl${}^{\mskip 1.5mu1}$, J. R. Gaunt${}^{\mskip 1.5mu2}$, P.
Pichini${}^{\mskip 1.5mu1,3}$ and P. Plößl${}^{\mskip 1.5mu1}$

1 Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany
2 Department of Physics and Astronomy, University of Manchester, Manchester,
M13 9PL,
2 United Kingdom
3 Department of Physics and Astronomy, Uppsala University, Box 516, 75120
Uppsala, Sweden

We investigate the positivity of double parton distributions with a non-
trivial dependence on the parton colour. It turns out that positivity is not
preserved by leading-order evolution from lower to higher scales, in contrast
to the case in which parton colour is summed over. We also study the
positivity properties of the distributions at small distance between the two
partons, where they can be computed in terms of perturbative splitting kernels
and ordinary parton densities.

###### Contents

  1. 1 Introduction
  2. 2 Colour structure of DPDs
  3. 3 Collins-Soper evolution
  4. 4 DGLAP evolution
    1. 4.1 Evolution in the scale of one parton
    2. 4.2 Simultaneous evolution in all scales
  5. 5 DPDs from parton splitting at two-loop accuracy
  6. 6 Summary
  7. A Colour space projectors
  8. Acknowledgements

##  1 Introduction

Parton distribution functions and related quantities are crucial ingredients
for computing hadronic cross sections at high energies, and they are the main
quantities that describe the structure of hadrons at the level of quarks and
gluons. It is hence important to know and understand their general properties.
One of these properties is positivity. For ordinary parton distributions
(PDFs) this is just the statement $f_{a}(x)\geq 0$ if the parton $a$ and the
hadron are unpolarised. In the polarised case, this generalises to a set of
inequalities, namely the well-known Soffer bounds on polarised distributions
[1]. Corresponding bounds have been formulated for transverse-momentum
dependent distributions (TMDs) in [2], for impact parameter distributions in
[3], and for double parton distributions (DPDs) in [4]. The latter appear in
the description of double parton scattering and contain a wealth of
information about correlations between two partons in a hadron; for a recent
review we refer to the monograph [5]. DPDs have a non-trivial dependence not
only on the polarisation of the partons, but also on their colour, and
corresponding positivity bounds have been derived in [6].

Positivity bounds can be of considerable practical value. They may be used as
constraints in fits of PDFs, and in the context of spin physics, ansätze that
saturate certain bounds are often used to estimate the maximal allowed size of
spin asymmetries. A corresponding strategy for DPDs in spin or colour space
appears all the more attractive because our current knowledge of these
distributions is very incomplete.

Whether positivity bounds on parton distributions actually hold turns out to
be a non-trivial question. At leading order (LO) accuracy, the positivity of
PDFs can quite directly be deduced from the positivity of cross sections,
whilst the situation is more involved at next-to-leading order (NLO) in the
strong coupling [7, 8]. When derived along such lines, positivity holds for
renormalisation scales $\mu$ that are high enough for the approximations of
leading-twist dominance and of the perturbative expansion to be valid. To
formulate a corresponding approach for DPDs would be complicated due to the
large number of involved degrees of freedom (two pairs of partons in each
colliding hadron and two hard-scattering subprocesses), and we shall not
pursue such an avenue here.

Intuitively, the positivity bounds on parton distributions are a consequence
of their parton-model interpretation as number densities or linear
combinations of number densities. At a more formal level, one may use light-
cone quantisation and write appropriate linear combinations of distributions
as squared operator matrix elements that are summed over unobserved degrees of
freedom. Equivalently, one can represent the distributions in terms of light-
cone wave functions. This has for instance been done for ordinary PDFs in [9]
and for impact parameter distributions in [10]. A corresponding representation
holds for DPDs (see [11] for the momentum space version, which can readily be
adapted to the case of definite transverse parton position). A limitation of
this approach is that it does not account for the renormalisation of
ultraviolet divergences in the matrix elements (at least not in customary
schemes such as $\overline{\text{MS}}$) nor for subtleties related with Wilson
lines or the definition of light-cone gauge at infinite light-cone distances.
It has long been realised that ultraviolet subtractions can in principle
invalidate the positivity of distributions. Interestingly, a recent
investigation in a toy model has shown that this can happen even at scales
much above the intrinsic non-perturbative scale of the model [12].

An important result is that LO DGLAP evolution to higher scales conserves the
positivity of PDFs, both in the unpolarised sector [13, 14] and in the
polarised one [15, 16]. This means that if PDFs satisfy the positivity bounds
at a certain scale $\mu$, these bounds remain valid at higher scales when the
PDFs are evolved at leading order. Discussions for NLO evolution can be found
in [17, 18]. Conversely, experience shows that PDFs eventually turn negative
when evolved down to very low scales. Examples for this can for instance be
found in [19]. In [4] it was shown that positivity of spin dependent but
colour summed DPDs is conserved by LO DGLAP evolution to higher scales. The
first goal of the present paper is to investigate whether the same holds for
the bounds derived in [6] for DPDs with non-trivial colour dependence. In
addition to DGLAP evolution, we will also consider Collins-Soper evolution in
a rapidity variable, which appears when the parton colours are not summed
over. Following [6] we will limit ourselves to unpolarised partons throughout
this work.

For small transverse distances $y$ between the two partons, DPDs can be
computed in terms of ordinary parton densities and kernels for the
perturbative splitting of one parton into the two observed partons and (at
higher orders) additional unobserved ones. Using the results of the recent
two-loop calculation in [20], we can investigate to which extent positivity of
DPDs in colour space is realised in the small $y$ limit. This is the second
goal of our work.

This paper is organised as follows. In section 2, we specify the two
parametrisations for the colour structure of DPDs used in this work, specify
the property of positivity, and discuss the perturbative splitting mechanism
for DPDs at LO accuracy. In section 3, we analyse whether Collins-Soper
evolution from smaller to larger scales conserves positivity of DPDs, and in
section 4 we perform a corresponding analysis for LO DGLAP evolution. The
perturbative splitting mechanism at NLO accuracy is investigated in section 5.
Our results are summarised in section 6, and some technical formulae are
collected in an appendix.

##  2 Colour structure of DPDs

In this section, we discuss the general colour structure of quark and gluon
DPDs and state the hypothesis of positivity in colour space. Throughout this
work, we consider distributions for unpolarised partons.

As illustrated in figure 1, the colour structure of a DPD can be described in
terms of four colour indices, one for each parton field in its definition. We
generically write $F_{a_{1}a_{2}}^{r_{1}r_{1}^{\prime}\,r_{2}r_{2}^{\prime}}$,
where $a_{1}$ and $a_{2}$ denote the two parton flavours. The colour indices
are in the fundamental or adjoint representation as appropriate, with $r_{1}$
and $r_{2}$ referring to the partons in the amplitude of the scattering
process and $r_{1}^{\prime}$ and $r_{2}^{\prime}$ to the partons in the
conjugate amplitude. As described in [21], the definition of a DPD involves
the hadronic matrix element of two twist-two operators, as well as a soft
factor given as the matrix element of Wilson line operators in the vacuum.
Both matrix elements contain ultraviolet divergences that need to be
renormalised. One can take different renormalisation scales $\mu_{1}$ and
$\mu_{2}$ for the two partons, and the dependence of the DPD on these scales
is given by DGLAP equations, which will be discussed in section 4. The soft
factor removes rapidity divergences in the hadronic matrix element, in a
similar way as in the definition of transverse-momentum dependent
distributions [22]. This leads to a dependence of the DPD on a rapidity scale
$\zeta_{p}$, which is described by a Collins-Soper equation as discussed in
section 3.

(a) $F_{q\bar{q}}^{ii^{\prime}\mskip 1.5mujj^{\prime}}$ (b)
$F_{qg}^{ii^{\prime}\mskip 1.5muaa^{\prime}}$

Figure 1:  Assignment of colour labels for a quark-antiquark distribution (a)
and a quark-gluon distribution (b). The dashed vertical line indicates the
final state cut of the scattering process in which the distributions appear.

#### The $s$ and $t$ channel colour bases.

The four colour indices of a DPD must be coupled to an overall colour singlet.
This can be achieved by first coupling the colour of two parton pairs to an
irreducible representation and then coupling these two representations to an
overall singlet. Depending on the choice of parton pairs, we consider two
bases for the colour coupling.

In the $s$ channel basis, we pair the partons in the amplitude and in the
conjugate amplitude. The projection on irreducible representations can then be
written as

$\displaystyle{F}^{RR^{\mskip 1.5mu\prime}}_{a_{1}a_{2}}$
$\displaystyle=\frac{1}{m(R)}\;P_{\rule{0.0pt}{4.82224pt}\overline{R}\,\overline{R}{}^{\mskip
1.5mu\prime}}^{r_{1}r_{2}\,r_{1}^{\prime}r_{2}^{\prime}}\,F_{a_{1}a_{2}\phantom{R}}^{r_{1}r_{1}^{\prime}\,r_{2}r_{2}^{\prime}}\qquad\text{for
}(a_{1}a_{2})\neq(q\bar{q}),(\bar{q}q)$ (2.1)

with colour indices $r_{1},r_{2},r_{1}^{\prime},r_{2}^{\prime}$ in the
fundamental or adjoint representation as appropriate. The multiplicity $m(R)$
of the representation $R$ is $m(1)=1$, $m(8)=m(A)=m(S)=8$,
$m(\overline{3})=m(3)=3$, etc. The normalisation of ${F}^{RR^{\mskip
1.5mu\prime}}_{a_{1}a_{2}}$ in (2.1) follows the choice made in [6] and
corresponds to eq. (6) in [23]. We denote the conjugate of a representation
$R$ by $\overline{R}$, where it is understood that some representations like
the singlet and the octet are their own conjugate.

The matrix $P_{\mskip 1.5mu\overline{R}\,\overline{R}{}^{\mskip 1.5mu\prime}}$
in (2.1) projects the colour of the parton pair $a_{1}a_{2}$ on the
representation $R$ in the amplitude and on the representation $R^{\mskip
1.5mu\prime}$ in the conjugate amplitude. Its explicit form is given in the
appendix. Quark-antiquark distributions $F_{q\bar{q}}^{RR^{\mskip
1.5mu\prime}}$ are defined as in (2.1), but with transposed indices
$r_{2}^{\prime}\mskip 1.5mur_{1}^{\prime}$ in $P_{\mskip
1.5mu\overline{R}\,\overline{R}{}^{\mskip 1.5mu\prime}}$, which ensures that
covariant indices are always contracted with contravariant ones. Likewise, the
definition of $F_{\bar{q}q}^{RR^{\mskip 1.5mu\prime}}$ has transposed indices
$r_{2}\mskip 1.5mur_{1}$.

In the $t$ channel basis, we pair the partons with momentum fractions $x_{1}$
and $x_{2}$. Following [24, 21], we write

$\displaystyle{}^{R_{1}R_{2}}\mskip-1.5muF_{a_{1}a_{2}}$
$\displaystyle=\mathcal{N}_{a_{1}}\mathcal{N}_{a_{2}}\,\varepsilon(R_{1})\,\varepsilon(R_{2})\,\frac{1}{\sqrt{m(R_{1})}}\;P_{\rule{0.0pt}{4.82224pt}\overline{R}_{1}\overline{R}_{2}}^{r_{1}r_{1}^{\prime}\,r_{2}r_{2}^{\prime}}\,F_{a_{1}a_{2}\phantom{R}}^{r_{1}r_{1}^{\prime}\,r_{2}r_{2}^{\prime}}\qquad\text{for
}a_{1},a_{2}\neq\bar{q}$ (2.2)

with prefactors $\mathcal{N}_{q}=\sqrt{3}$, $\mathcal{N}_{g}=\sqrt{8}$,
$\varepsilon(A)=i$ and $\varepsilon(R)=1$ for $R\neq A$. The matrix $P_{\mskip
1.5mu\overline{R}_{1}\overline{R}_{2}}$ in (2.2) projects the colour indices
$r_{i}r_{i}^{\prime}$ of parton $a_{i}$ on the representation $R_{i}$ for
$i=1,2$. For distributions with antiquarks, one needs to transpose the
corresponding colour indices in $P_{\mskip
1.5mu\overline{R}_{1}\overline{R}_{2}}$. The correct ordering is hence
$r_{1}^{\prime}\mskip 1.5mur_{1}$ if $a_{1}=\bar{q}$ and $r_{2}^{\prime}\mskip
1.5mur_{2}$ if $a_{2}=\bar{q}$.

Combining the definition (2.1) with the completeness relation (A) and the
explicit form of the singlet projector $P_{11}$, one readily finds that

$\displaystyle{}^{11}\mskip-1.5muF_{a_{1}a_{2}}$
$\displaystyle=\sum_{R}m(R)\,F^{R\overline{R}}_{a_{1}a_{2}}$ (2.3)

for all parton combinations, where the sum runs over all relevant colour
representations $R$.

Throughout this work, we fix $N=3$ for the number of colours. The colour
factors used in later results thus have the values $C_{F}=4/3$, $C_{A}=3$, and
$T_{F}=1/2$. The accessible colour representations for different parton
combinations are given in table 1. We note that a quark and a gluon can couple
to $\overline{6}$ (see e.g. table 24 in [25]) rather than to $6$ as stated in
equation (8c) of [23]. As shown in [21], the $t$ channel distributions
${}^{R_{1}R_{2}}\mskip-1.5muF$ are real valued, except for the decuplet sector
in the pure gluon case, in which case one has $({}^{10\mskip
1.5mu\overline{10}}\mskip-1.5muF_{gg})^{*}={}^{\overline{10}\mskip
1.5mu10}\mskip-1.5muF_{gg}$. In the $s$ channel basis, this translates into
all distributions $F^{R_{1}R_{2}}$ being real, except for the mixed octet
combinations, where one finds $(F_{gg}^{AS})^{*}=F_{gg}^{SA}$.

$\begin{array}[]{c c c}&s\leavevmode\nobreak\
\text{channel}&t\leavevmode\nobreak\ \text{channel}\\\ \hline\cr
a_{1}a_{2}&RR^{\mskip 1.5mu\prime}&R_{1}R_{2}\\\ \hline\cr
q\bar{q}&11,88&11,88\\\ qq&\overline{3}\mskip 1.5mu3,6\mskip
1.5mu\overline{6}&11,88\\\ \bar{q}\bar{q}&3\mskip
1.5mu\overline{3},\overline{6}\mskip 1.5mu6&11,88\\\ gq&3\mskip
1.5mu\overline{3},\overline{6}\mskip 1.5mu6,15\mskip
1.5mu\overline{15}&11,S8,A8\\\ g\bar{q}&\overline{3}\mskip 1.5mu3,6\mskip
1.5mu\overline{6},\overline{15}\mskip 1.5mu15&11,S8,A8\\\ gg&\lx@intercol\hfil
11,SS,AA,S\mskip-1.5muA,AS,10\mskip 1.5mu\overline{10},\overline{10}\mskip
1.5mu10,27\mskip 1.5mu27\hfil\lx@intercol\\\ \hline\cr\end{array}$ Table 1:
Combinations of colour representations in the $s$ channel distributions
$F^{RR^{\mskip 1.5mu\prime}}_{a_{1}a_{2}}$ and the $t$ channel distributions
${}^{R_{1}R_{2}}\mskip-1.5muF_{a_{1}a_{2}}$. Two adjoint indices can couple to
a symmetric ($S$) or an antisymmetric ($A$) octet. For two-gluon
distributions, the colour combinations are identical in the $s$ and $t$
channels. Interchanging $a_{1}\leftrightarrow a_{2}$ implies interchanging
$R_{1}\leftrightarrow R_{2}$ while keeping $RR^{\mskip 1.5mu\prime}$
unchanged.

#### Density interpretation and positivity.

The parton model interpretation of DPDs can be obtained in the same way as for
single parton distributions by expressing the field operators in terms of
creation and annihilation operators in light-cone quantisation, neglecting all
complications from Wilson lines and from renormalisation. Details can for
instance be found in [22].

The $t$ channel distributions are normalised such that
${}^{11}\mskip-1.5muF_{a_{1}a_{2}}(x_{1},x_{2},y)$ is the probability density
for finding partons $a_{1}$ and $a_{2}$ with momentum fractions $x_{1}$ and
$x_{2}$ at a transverse distance $y$ from each other, with the density measure
being $dx_{1}\,dx_{2}\,d^{2}y$. The colours and polarisations of both partons
are summed over in ${}^{11}\mskip-1.5muF_{a_{1}a_{2}}$. Correspondingly, the
$s$ channel distribution $F_{a_{1}a_{2}}^{R\overline{R}}$ is the probability
density for finding the parton pair in one of the $m(R)$ states of the colour
representation $R$. This provides an intuitive interpretation of the relation
(2.3).

The positivity property for DPDs in full colour space is then the statement
that

$\displaystyle F^{R\overline{R}}_{a_{1}a_{2}}$ $\displaystyle\geq 0$
$\displaystyle\text{for all }R,$ (2.4)

which of course implies the weaker condition
${}^{11}\mskip-1.5muF_{a_{1}a_{2}}\geq 0$. Note that we define “positivity” as
including the value zero.

Note that in the pure gluon channel, the distributions in the $s$ channel
basis include the cases $F^{AS}_{gg}$ and $F^{S\mskip-1.5muA}_{gg}$, which
correspond not to densities but to interference terms in colour space (and
which may be complex valued as noted above). Given the large number of
accessible colour channels in that case, we will not consider two-gluon DPDs
in the remainder of this work, concentrating on the pure quark-antiquark
sector and on mixed quark-gluon or antiquark-gluon distributions.

#### Basis transformations.

Whilst the $s$ channel basis is natural for considering positivity, the
evolution of DPDs in the renormalisation and rapidity scales is much simpler
in the $t$ channel basis. We hence need the explicit transformations between
the two representations. In the pure quark sector, the transformations between
$s$ and $t$ channel bases read

$\displaystyle\begin{pmatrix}F_{qq}^{\overline{3}\mskip 1.5mu3}\\\\[3.00003pt]
F_{qq}^{6\mskip 1.5mu\overline{6}}\end{pmatrix}$
$\displaystyle=\mathbf{M}_{qq}\,\begin{pmatrix}{}^{11}\mskip-1.5muF_{qq}\\\\[3.99994pt]
{}^{88}\mskip-1.5muF_{qq}\end{pmatrix}\,,$
$\displaystyle\begin{pmatrix}F_{q\bar{q}}^{11}\\\\[3.00003pt]
F_{q\bar{q}}^{88}\end{pmatrix}$
$\displaystyle=\mathbf{M}_{q\bar{q}}\,\begin{pmatrix}{}^{11}\mskip-1.5muF_{q\bar{q}}\\\\[3.99994pt]
{}^{88}\mskip-1.5muF_{q\bar{q}}\end{pmatrix}\,,$ (2.5)

with

$\displaystyle\mathbf{M}_{qq}$
$\displaystyle=\frac{1}{9}\,\begin{pmatrix}1&-\sqrt{2}\,\\\\[1.99997pt]
1&\frac{1}{\sqrt{2}}\end{pmatrix}\,,$ $\displaystyle\mathbf{M}_{q\bar{q}}$
$\displaystyle=\frac{1}{9}\,\begin{pmatrix}1&\sqrt{8}\\\\[1.99997pt]
1&-\frac{1}{\sqrt{8}}\end{pmatrix}\,,$ $\displaystyle\mathbf{M}_{qq}^{-1}$
$\displaystyle=3\,\begin{pmatrix}1&2\\\\[1.99997pt]
-\sqrt{2}&\sqrt{2}\end{pmatrix}\,,$ $\displaystyle\mathbf{M}_{q\bar{q}}^{-1}$
$\displaystyle=\begin{pmatrix}1&8\\\\[1.99997pt]
\sqrt{8}&-\sqrt{8}\end{pmatrix}\,.$ (2.6)

The transformation matrix for two antiquarks is
$\mathbf{M}_{\bar{q}\bar{q}}=\mathbf{M}_{qq}$. For $gq$ and $g\bar{q}$
distributions, one has

$\displaystyle\begin{pmatrix}F_{gq}^{3\mskip 1.5mu\overline{3}}\\\\[3.00003pt]
F_{gq}^{\overline{6}\mskip 1.5mu6}\\\\[3.00003pt] F_{gq}^{15\mskip
1.5mu\overline{15}}\end{pmatrix}$
$\displaystyle=\mathbf{M}_{gq}\,\begin{pmatrix}{}^{11}\mskip-1.5muF_{gq}\\\\[3.99994pt]
{}^{S8}\mskip-1.5muF_{gq}\\\\[3.99994pt]
{}^{A8}\mskip-1.5muF_{gq}\end{pmatrix}\,,$
$\displaystyle\begin{pmatrix}F_{g\bar{q}}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{g\bar{q}}^{6\mskip 1.5mu\overline{6}}\\\\[3.00003pt]
F_{g\bar{q}}^{\overline{15}\mskip 1.5mu15}\end{pmatrix}$
$\displaystyle=\mathbf{M}_{g\bar{q}}\,\begin{pmatrix}{}^{11}\mskip-1.5muF_{\bar{q}g}\\\\[3.99994pt]
{}^{S8}\mskip-1.5muF_{g\bar{q}}\\\\[3.99994pt]
{}^{A8}\mskip-1.5muF_{g\bar{q}}\end{pmatrix}$ (2.7)

with

$\displaystyle\mathbf{M}_{gq}$
$\displaystyle=\frac{1}{24}\,\begin{pmatrix}1&\sqrt{\frac{5}{2}}&-\frac{3}{\sqrt{2}}\\\\[5.0pt]
1&-\sqrt{\frac{5}{2}}&-\frac{1}{\sqrt{2}}\\\\[5.0pt]
1&\sqrt{\frac{1}{10}}&\frac{1}{\sqrt{2}}\end{pmatrix}\,,$
$\displaystyle\mathbf{M}_{g\bar{q}}$
$\displaystyle=\frac{1}{24}\,\begin{pmatrix}1&\sqrt{\frac{5}{2}}&\frac{3}{\sqrt{2}}\\\\[5.0pt]
1&-\sqrt{\frac{5}{2}}&\frac{1}{\sqrt{2}}\\\\[5.0pt]
1&\sqrt{\frac{1}{10}}&-\frac{1}{\sqrt{2}}\end{pmatrix}\,,$
$\displaystyle\mathbf{M}_{gq}^{-1}$
$\displaystyle=3\,\begin{pmatrix}1&2&5\\\\[3.99994pt]
\sqrt{\frac{5}{2}}&-\sqrt{10}&\sqrt{\frac{5}{2}}\\\\[3.99994pt]
-\frac{3}{\sqrt{2}}&-\sqrt{2}&\frac{5}{\sqrt{2}}\end{pmatrix}\,,$
$\displaystyle\mathbf{M}_{g\bar{q}}^{-1}$
$\displaystyle=3\,\begin{pmatrix}1&2&5\\\\[3.99994pt]
\sqrt{\frac{5}{2}}&-\sqrt{10}&\sqrt{\frac{5}{2}}\\\\[3.99994pt]
\frac{3}{\sqrt{2}}&\sqrt{2}&-\frac{5}{\sqrt{2}}\end{pmatrix}\,.$ (2.8)

All transformation matrices fulfil the symmetry property
$\mathbf{M}_{a_{1}a_{2}}=\mathbf{M}_{a_{2}\mskip 1.5mua_{1}}$. The expressions
of $\mathbf{M}_{qq}$, $\mathbf{M}_{q\bar{q}}$, and $\mathbf{M}_{gq}$ were
already given in [6].

#### Perturbative splitting at leading order.

If the transverse distance $y$ between the two partons is small, DPDs can be
computed in terms of a perturbative splitting process and ordinary PDFs.
Example graphs for the perturbative splitting are shown in figure 2. This
mechanism is interesting in our context because it generates a non-trivial
colour dependence.

Let us take a closer look at the splitting process at one-loop order,
postponing the discussion of two-loop accuracy to section 5. The splitting
formula at leading order reads [24, 21]

$\displaystyle{}^{R_{1}R_{2}}\mskip-1.5muF_{a_{1}a_{2}}(x_{1},x_{2},y,\mu,\mu)$
$\displaystyle=\frac{a_{s}(\mu)}{\pi\mskip
1.5muy^{2}}\;{}^{R_{1}R_{2}}\mskip-1.5muV_{a_{1}a_{2},a_{0}}^{(1)}(u)\;\frac{f_{a_{0}}(x,\mu)}{x}\,,$
(2.9)

where $f_{a_{0}}(x,\mu)$ is the PDF for parton $a_{0}$ and we defined

$\displaystyle a_{s}(\mu)$ $\displaystyle=\frac{\alpha_{s}(\mu)}{2\pi}$ (2.10)

and

$\displaystyle x$ $\displaystyle=x_{1}+x_{2}\,,$ $\displaystyle u$
$\displaystyle=\frac{x_{1}}{x_{1}+x_{2}}\,.$ (2.11)

Note that (2.9) is an approximation for small $y$ and receives corrections
suppressed by a power of $a_{s}$ or $y^{2}\Lambda^{2}$, where $\Lambda$ is a
hadronic scale.

The splitting kernels ${}^{11}\mskip-1.5muV_{a_{1}a_{2},a_{0}}^{(1)}(z)$ for
the colour singlet channel are equal to the usual LO DGLAP splitting functions
without the distributional parts (plus prescription and delta function) at
$z=1$. One therefore has ${}^{11}\mskip-1.5muV_{a_{1}a_{2},a_{0}}^{(1)}(z)>0$.
This implies that ${}^{11}\mskip-1.5muF_{a_{1}a_{2}}>0$ at the scale $\mu$
where the LO splitting formula (2.9) is evaluated, provided of course that the
PDFs are positive. To avoid large $a_{s}^{2}$ corrections, one should take
$\mu\sim 1/y$.

The LO splitting kernels for other colour channels are proportional to
${}^{11}\mskip-1.5muV_{a_{1}a_{2},a_{0}}^{(1)}(z)$, and it is easy to
transform them to the $s$ channel colour basis. The result is

$\displaystyle F_{a_{1}a_{2}}^{R\overline{R}}$
$\displaystyle=\begin{cases}{}^{11}\mskip-1.5muF_{a_{1}a_{2}}\big{/}m(R_{0})&\text{for
}R=R_{0}\,,\\\ 0&\text{for }R\neq R_{0}\,,\end{cases}$ (2.12)

where $R_{0}$ is the colour representation of the initial parton $a_{0}$ of
the splitting process (with $R_{0}=A$ for $g\to gg$). The parton pair
$a_{1}a_{2}$ must be in the representation $R_{0}$ because the LO splitting
graphs are disconnected between the amplitude and conjugate amplitude (see
figure 2(a)). The result (2.12) then follows from the relation (2.3). With
${}^{11}\mskip-1.5muF_{a_{1}a_{2}}>0$ one hence finds positivity of DPDs in
colour space when taking the LO splitting approximation. We will see in
section 5 whether this still holds at NLO.

(a) (b)

(c)

Figure 2:  Graph for the perturbative splitting $g\to q\bar{q}$ in a DPD.
Panel (a) shows the LO graph, and panels (b) and (c) are examples of real and
virtual NLO graphs. The dashed vertical line indicates the final state cut of
the scattering process in which the distributions appear.

##  3 Collins-Soper evolution

In this section, we investigate how Collins-Soper evolution affects the
positivity of DPDs.

Collins-Soper evolution of DPDs does not mix different colour representations
in the $t$ channel basis, where one has

$\displaystyle\frac{\partial}{\partial\log\zeta_{p}}\,{}^{R_{1}R_{2}}F_{a_{1}a_{2}}(x_{1},x_{2},y,\mu_{1},\mu_{2},\zeta_{p})$
$\displaystyle=\frac{{}^{R_{1}}\mskip-1.5muJ(y,\mu_{1},\mu_{2})}{2}\;{}^{R_{1}R_{2}}F_{a_{1}a_{2}}(x_{1},x_{2},y,\mu_{1},\mu_{2},\zeta_{p})\,.$
(3.1)

Here we displayed all arguments of the functions. The Collins-Soper kernel
${}^{R_{1}}\mskip-1.5muJ$ depends only on the multiplicity of $R_{1}$ (which
is equal to the multiplicity of $R_{2}$) but not on the parton types. Note
that colour singlet distributions in the $t$ channel are $\zeta_{p}$
independent, i.e. ${}^{1}\mskip-1.5muJ=0$. For all parton combinations except
$gg$, the only non-trivial kernel needed is hence the one for the colour
octet. Remarkably, this kernel satisfies the exact relation [26]

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

where $K_{g}(y,\mu)$ is the Collins-Soper kernel for the evolution of single-
gluon TMDs.

Let us discuss the sign of the Collins-Soper kernel, which will be important
in the following. The renormalisation group equation for the Collins-Soper
kernel is solved by

$\displaystyle{}^{8}\mskip-1.5muJ(y,\mu_{1},\mu_{2})$
$\displaystyle={}^{8}\mskip-1.5muJ(y,\mu_{0},\mu_{0})-\int_{\mu_{0}}^{\mu_{1}}\frac{d\mu}{\mu}\;{}^{8}\mskip-1.5mu\gamma_{J}(\mu)-\int_{\mu_{0}}^{\mu_{2}}\frac{d\mu}{\mu}\;{}^{8}\mskip-1.5mu\gamma_{J}(\mu)$
(3.3)

with a positive anomalous dimension that is proportional to the cusp anomalous
dimension for adjoint Wilson lines:

$\displaystyle{}^{8}\mskip-1.5mu\gamma_{J}(\mu)$
$\displaystyle=a_{s}(\mu)\;{}^{8}\mskip-1.5mu\gamma_{J}^{(0)}+\mathcal{O}(a_{s}^{2})\,,$
$\displaystyle{}^{8}\mskip-1.5mu\gamma_{J}^{(0)}$ $\displaystyle=2C_{A}\,.$
(3.4)

At given $y$ one can hence always achieve a negative ${}^{8}\mskip-1.5muJ$ by
taking the scales $\mu_{1}$ and $\mu_{2}$ sufficiently high.

To make a more specific statement, we first consider small distances $y$,
where one can compute the kernel in perturbation theory and obtains

$\displaystyle{}^{8}\mskip-1.5muJ(y,\mu,\mu)$
$\displaystyle=-a_{s}(\mu)\;{}^{8}\mskip-1.5mu\gamma_{J}^{(0)}\,\log\frac{\mu^{2}\mskip
1.5muy^{2}}{b_{0}^{2}}+\mathcal{O}(a_{s}^{2})\,,$ (3.5)

where $b_{0}=2e^{-\gamma}\approx 1.12$ and $\gamma$ is the Euler-Mascheroni
constant. Bearing in mind that there are higher-order terms in (3.5), we see
that the transition from positive to negative ${}^{8}\mskip-1.5muJ$ happens at
$\mu$ around $b_{0}/y$, as long as $y$ remains in the perturbative regime.

Not much is known about ${}^{8}\mskip-1.5muJ$ or $K_{g}$ for $y$ in the
nonperturbative domain. The situation is different for the Collins-Soper
kernel $K_{q}$ for quark TMDs. Several phenomenological extractions find that
$K_{q}(y,\mu)$ is negative for large $y$, see e.g. [27] (figure 6), [28] and
[29] (figure 23). Furthermore a number of lattice determinations, covering a
distance range between about $0.1\operatorname{fm}$ and
$0.8\operatorname{fm}$, find that $K_{q}(y,\mu)<0$ at
$\mu=2\operatorname{GeV}$, see figure 7 in [30], figure 5 in [31], and figure
8 in [32].

We find it plausible to assume a qualitatively similar behaviour of
$K_{g}(y,\mu)$ and $K_{q}(y,\mu)$ as functions of $y$ (at perturbatively small
$y$, one actually has $K_{g}/K_{q}\approx C_{A}/C_{F}$ up to corrections of
order $\alpha_{s}^{4}$, see footnote 10 in [21]). Under this assumption, we
conclude that ${}^{8}\mskip-1.5muJ(y,\mu_{1},\mu_{2})<0$ for scales $\mu_{1}$
and $\mu_{2}$ sufficiently larger than $\max(b_{0}/y,2\operatorname{GeV})$.

#### Collins-Soper evolution in the $s$ channel.

In the $s$ channel, different colour representations mix under Collins-Soper
evolution. Starting with the $qq$ channel and using the basis transform given
in the previous section, we get

$\displaystyle\frac{\partial}{\partial\log\zeta_{p}}\,\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[1.99997pt] F_{qq}^{6\mskip 1.5mu\overline{6}}\end{pmatrix}$
$\displaystyle=\mathbf{M}_{qq}\;\frac{\partial}{\partial\log\zeta_{p}}\,\begin{pmatrix}{}^{11}\mskip-1.5muF_{qq}\\\\[1.99997pt]
{}^{88}\mskip-1.5muF_{qq}\end{pmatrix}=\mathbf{M}_{qq}\;\frac{1}{2}\begin{pmatrix}0&0\\\\[1.99997pt]
0&{}^{8}\mskip-1.5muJ\end{pmatrix}\,\begin{pmatrix}{}^{11}\mskip-1.5muF_{qq}\\\\[1.99997pt]
{}^{88}\mskip-1.5muF_{qq}\end{pmatrix}$
$\displaystyle=\frac{{}^{8}\mskip-1.5muJ}{2}\;\hat{\mathbf{J}}_{qq}\,\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[1.99997pt] F_{qq}^{6\mskip 1.5mu\overline{6}}\end{pmatrix}$ (3.6)

with

$\displaystyle\hat{\mathbf{J}}_{qq}=\mathbf{M}_{qq}\,\begin{pmatrix}0&0\\\\[1.00006pt]
0&1\end{pmatrix}\,\mathbf{M}_{qq}^{-1}$
$\displaystyle=\frac{1}{3}\,\begin{pmatrix}2&-2\\\\[1.00006pt]
-1&1\end{pmatrix}\,.$ (3.7)

Restoring all arguments, we find that the evolution equation is solved by

$\displaystyle\begin{pmatrix}F_{qq}^{\overline{3}\mskip 1.5mu3}\\\\[1.99997pt]
F_{qq}^{6\mskip
1.5mu\overline{6}}\end{pmatrix}(x_{1},x_{2},y,\mu_{1},\mu_{2},\zeta_{p})$
$\displaystyle=\mathbf{U}_{qq}(\alpha)\,\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[1.99997pt] F_{qq}^{6\mskip
1.5mu\overline{6}}\end{pmatrix}(x_{1},x_{2},y,\mu_{1},\mu_{2},\zeta_{0})$
(3.8)

with the matrix exponential

$\displaystyle\mathbf{U}_{qq}(\alpha)=\exp\bigl{(}-\alpha\,\hat{\mathbf{J}}_{qq}\bigr{)}$
$\displaystyle=\frac{1}{3}\,\begin{pmatrix}1+2e^{-\alpha}&2(1-e^{-\alpha})\\\\[1.00006pt]
1-e^{-\alpha}&2+e^{-\alpha}\end{pmatrix}\,,$ (3.9)

where we abbreviate

$\displaystyle\alpha$
$\displaystyle=-\frac{{}^{8}\mskip-1.5muJ(y,\mu_{1},\mu_{2})}{2}\,\log\frac{\zeta_{p}}{\zeta_{0}}\,.$
(3.10)

The Collins-Soper equation for the other parton combinations we consider has
the same form as (3) with appropriate changes in the colour labels and
matrices. In the quark-antiquark case, one has

$\displaystyle\hat{\mathbf{J}}_{q\bar{q}}=\mathbf{M}_{q\bar{q}}\,\begin{pmatrix}0&0\\\\[1.00006pt]
0&1\end{pmatrix}\,\mathbf{M}_{q\bar{q}}^{-1}$
$\displaystyle=\frac{1}{9}\,\begin{pmatrix}8&-8\\\\[1.00006pt]
-1&1\end{pmatrix}\,,$
$\displaystyle\mathbf{U}_{q\bar{q}}(\alpha)=\exp\bigl{(}-\alpha\,\hat{\mathbf{J}}_{q\bar{q}}\bigr{)}$
$\displaystyle=\frac{1}{9}\,\begin{pmatrix}1+8e^{-\alpha}&8(1-e^{-\alpha})\\\\[1.00006pt]
1-e^{-\alpha}&8+e^{-\alpha}\end{pmatrix}\,,$ (3.11)

and for $gq$ distributions we find

$\displaystyle\hat{\mathbf{J}}_{gq}=\mathbf{M}_{gq}\,\begin{pmatrix}0&0&0\\\\[1.00006pt]
0&1&0\\\\[1.00006pt] 0&0&1\end{pmatrix}\,\mathbf{M}_{gq}^{-1}$
$\displaystyle=\frac{1}{8}\,\begin{pmatrix}7&-2&-5\\\\[1.99997pt]
-1&6&-5\\\\[1.99997pt] -1&-2&3\end{pmatrix}\,,$
$\displaystyle\mathbf{U}_{gq}(\alpha)=\exp\bigl{(}-\alpha\,\hat{\mathbf{J}}_{gq}\bigr{)}$
$\displaystyle=\frac{1}{8}\,\begin{pmatrix}1+7e^{-\alpha}&2(1-e^{-\alpha})&5(1-e^{-\alpha})\\\\[1.00006pt]
1-e^{-\alpha}&2+6e^{-\alpha}&5(1-e^{-\alpha})\\\\[1.00006pt]
1-e^{-\alpha}&2(1-e^{-\alpha})&5+3e^{-\alpha}\,.\end{pmatrix}\,,$ (3.12)

where $\alpha$ is always given by (3.10). We furthermore get
$\hat{\mathbf{J}}_{\bar{q}\bar{q}}=\hat{\mathbf{J}}_{qq}$ and
$\hat{\mathbf{J}}_{g\bar{q}}=\hat{\mathbf{J}}_{gq}$ and corresponding
equalities for the evolution matrices $\mathbf{U}_{a_{1}a_{2}}$. The evolution
equations for distributions $F_{q_{1}q_{2}}$, $F_{q_{1}\bar{q}_{2}}$, etc.
with unequal flavours involve the same matrices as their counterparts for
equal flavours.

We see that for all parton combinations $a_{1}a_{2}$ except $gg$ (which we do
not consider) all elements of the evolution matrix $U_{a_{1}a_{2}}(\alpha)$
are positive for $\alpha>0$. This is the case for forward evolution
($\zeta_{p}>\zeta_{0}$), provided that ${}^{8}\mskip-1.5muJ<0$, which is the
case when $\mu_{1}$ and $\mu_{2}$ are sufficiently large. Under this
condition, Collins-Soper evolution to higher scales thus preserves positivity.

With the notation

$\displaystyle\bigl{(}\hat{\mathbf{J}}_{a_{1}a_{2}}\mskip
1.5muF_{a_{1}a_{2}}\bigr{)}{}^{R\overline{R}}$ $\displaystyle=\sum_{R^{\mskip
1.5mu\prime}}\bigl{(}\hat{\mathbf{J}}_{a_{1}a_{2}}\bigr{)}{}^{R\overline{R},\,R^{\mskip
1.5mu\prime}\mskip-1.5mu\overline{R}{}^{\mskip
1.5mu\prime}}\,F_{a_{1}a_{2}}^{R^{\mskip
1.5mu\prime}\mskip-1.5mu\overline{R}{}^{\mskip 1.5mu\prime}}$ (3.13)

the Collins-Soper equation in the $s$ channel basis reads

$\displaystyle\frac{\partial}{\partial\log\zeta_{p}}\mskip
1.5muF_{a_{1}a_{2}}^{R\overline{R}}$
$\displaystyle=\frac{{}^{8}\mskip-1.5muJ}{2}\,\bigl{(}\hat{\mathbf{J}}_{a_{1}a_{2}}\mskip
1.5muF_{a_{1}a_{2}}\bigr{)}{}^{R\overline{R}}_{\phantom{a_{1}}}$ (3.14)

for all parton combinations considered here. Using that
$(\hat{\mathbf{J}}_{a_{1}a_{2}})^{2}=\hat{\mathbf{J}}_{a_{1}a_{2}}$, we can
write its solution in the form

$\displaystyle F_{a_{1}a_{2}}^{R\overline{R}}(\zeta_{p})$
$\displaystyle=F_{a_{1}a_{2}}^{R\overline{R}}(\zeta_{0})+(e^{-\alpha}-1)\,\bigl{(}\hat{\mathbf{J}}_{a_{1}a_{2}}\mskip
1.5muF_{a_{1}a_{2}}(\zeta_{0})\bigr{)}{}^{R\overline{R}}\,.$ (3.15)

Writing the relation at (2.3) at rapidity scale $\zeta_{p}$ and using that
${}^{11}\mskip-1.5muF_{a_{1}a_{2}}$ is independent of $\zeta_{p}$, we get

$\displaystyle\sum_{R}m(R)\,F_{a_{1}a_{2}}^{R\overline{R}}(\zeta_{p})$
$\displaystyle={}^{11}\mskip-1.5muF_{a_{1}a_{2}}(\zeta_{0})\,.$ (3.16)

We can now discuss the behaviour of Collins-Soper evolution for large negative
$\alpha$, which is relevant when evolving backward with
${}^{8}\mskip-1.5muJ<0$, and when evolving forward at scales $\mu_{1}$ and
$\mu_{2}$ that are so low that ${}^{8}\mskip-1.5muJ>0$. In the regime where
the factor $e^{-\alpha}$ in (3.15) is much larger than $1$, the condition
(3.16) implies that the evolved distributions
$F_{a_{1}a_{2}}^{R\overline{R}}(\zeta_{p})$ must be positive in some colour
channels and negative in others. An exception to this statement is the case
where the initial conditions satisfy $(\hat{\mathbf{J}}_{a_{1}a_{2}}\mskip
1.5muF_{a_{1}a_{2}}(\zeta_{0}))^{R\overline{R}}=0$ for all $R$, so that all
distributions are independent of $\zeta_{p}$. In the $t$ channel basis, this
is tantamount to all distributions other than
${}^{11}\mskip-1.5muF_{a_{1}a_{2}}$ being zero. Apart from this special case,
evolution to large negative $\alpha$ always leads to a violation of
positivity.

##  4 DGLAP evolution

In this section, we investigate how leading-order DGLAP evolution affects the
positivity of DPDs. We limit ourselves to the evolution of two-quark and
quark-antiquark distributions, which are the simplest cases as far as mixing
and the number of colour channels are concerned. We first discuss evolution in
$\mu_{1}$ at fixed $\mu_{2}$ and $\zeta_{p}$ and then evolution in all three
scales simultaneously.

###  4.1 Evolution in the scale of one parton

Let us consider evolution in the renormalisation scale of one parton, which we
take to be the first one without loss of generality. For the parton
combinations of interest, the LO evolution equations in the $t$ channel basis
read

$\displaystyle\frac{\partial}{\partial\log\mu_{1}^{2}}\,{}^{R_{1}R_{2}}\mskip-1.5muF_{qa}(x_{1},x_{2},y,\mu_{1},\mu_{2},\zeta_{p})={}^{R_{1}R_{1}}\mskip-1.5muP_{qq}(z,\mu_{1},x_{1}^{2}\zeta_{p})\underset{x_{1}}{\otimes}{}^{R_{1}R_{2}}\mskip-1.5muF_{qa}(z,x_{2},y,\mu_{1},\mu_{2},\zeta_{p})$
$\displaystyle\qquad+\sum_{R^{\prime}}{}^{R_{1}R^{\mskip
1.5mu\prime}}\mskip-1.5muP_{qg}(z,\mu_{1},x_{1}^{2}\zeta_{p})\underset{x_{1}}{\otimes}{}^{R^{\mskip
1.5mu\prime}R_{2}}\mskip-1.5muF_{ga}(z,x_{2},y,\mu_{1},\mu_{2},\zeta_{p})\,,$
$\displaystyle\frac{\partial}{\partial\log\mu_{1}^{2}}\,{}^{R_{1}R_{2}}\mskip-1.5muF_{\bar{q}a}(x_{1},x_{2},y,\mu_{1},\mu_{2},\zeta_{p})={}^{R_{1}R_{1}}\mskip-1.5muP_{qq}(z,\mu_{1},x_{1}^{2}\zeta_{p})\underset{x_{1}}{\otimes}{}^{R_{1}R_{2}}\mskip-1.5muF_{\bar{q}a}(z,x_{2},y,\mu_{1},\mu_{2},\zeta_{p})$
$\displaystyle\qquad+\sum_{R^{\prime}}\varepsilon^{2}(R^{\prime})\,{}^{R_{1}R^{\mskip
1.5mu\prime}}\mskip-1.5muP_{qg}(z,\mu_{1},x_{1}^{2}\zeta_{p})\underset{x_{1}}{\otimes}{}^{R^{\mskip
1.5mu\prime}R_{2}}\mskip-1.5muF_{ga}(z,x_{2},y,\mu_{1},\mu_{2},\zeta_{p})\,,$
(4.1)

where the second parton $a$ is a quark or an antiquark, and where
$R^{\prime}=1$ if $R_{1}=1$ and $R^{\prime}=A,S$ if $R_{1}=8$. In the second
equation we have made use of the charge conjugation relations between
$\bar{q}g$ and $qg$ splitting kernels; this gives a sign factor
$\varepsilon^{2}(A)=-1$ for the antisymmetric octet.

The evolution equations involve the Mellin convolution

$\displaystyle P(z,\ldots)\underset{x_{1}}{\otimes}F(z,x_{2},\ldots)$
$\displaystyle=\int_{x_{1}/(1-x_{2})}^{1}\frac{dz}{z}\,P(z,\ldots)\,F\biggl{(}\frac{x_{1}}{z},x_{2},\ldots\biggr{)}\,,$
(4.2)

whose lower integration boundary reflects the support region of DPDs in the
momentum fractions: $F(z_{1},z_{2},\ldots)$ is zero for $z_{1}+z_{2}>1$. The
evolution kernels can be written as

$\displaystyle{}^{RR^{\prime}}\\!P_{qb}(z,\mu,\zeta)$
$\displaystyle=a_{s}(\mu)\,c_{RR^{\prime}}\,\widetilde{P}_{qb}(z)+\Biggl{[}\frac{\gamma_{q}(\mu)}{2}+\frac{{}^{R}\mskip-1.5mu\gamma_{J}(\mu)}{4}\,\log\frac{\mu^{2}}{\zeta}\Biggr{]}\,\delta_{RR^{\mskip
1.5mu\prime}}\,\delta_{qb}\,\delta(1-z)+\mathcal{O}(a_{s}^{2})\,,$ (4.3)

where $b=q,g$, $\gamma_{q}(\mu)=3C_{F}\mskip
1.5mua_{s}(\mu)+\mathcal{O}(a_{s}^{2})$,
${}^{1}\mskip-1.5mu\gamma_{J}(\mu)=0$, and ${}^{8}\mskip-1.5mu\gamma_{J}(\mu)$
is given in (3.4). The $z$ dependent part of (4.3) involves the familiar
splitting functions

$\displaystyle\widetilde{P}_{qq}(z)$
$\displaystyle=C_{F}\,\frac{1+z^{2}}{(1-z)_{+}}\,,$
$\displaystyle\widetilde{P}_{qg}(z)$
$\displaystyle=T_{F}\,\bigl{[}z^{2}+(1-z)^{2}\bigr{]}$ (4.4)

and the colour factors

$\displaystyle c_{11}$ $\displaystyle=1\,,$ $\displaystyle c_{88}$
$\displaystyle=-1/8\,,$ $\displaystyle c_{8S}$ $\displaystyle=\sqrt{5}/4\,,$
$\displaystyle c_{8A}$ $\displaystyle=3/4\,.$ (4.5)

Whilst our analysis is limited to the LO approximation of the splitting
kernels, our results do not depend on whether one uses the LO or the NLO
approximation for the anomalous dimension ${}^{8}\mskip-1.5mu\gamma_{J}(\mu)$,
which is associated with Sudakov double logarithms. Taking such anomalous
dimensions at NLO corresponds to next-to-leading logarithmic (NLL)
approximation (see e.g. table 1 in [33] for single hard scattering and section
6.6 in [21] for DPS). If one takes ${}^{8}\mskip-1.5mu\gamma_{J}(\mu)$ at two-
loop accuracy, one may also want to use the two-loop rather than the one-loop
running of $\alpha_{s}(\mu)$. Our arguments in the present work do not depend
on that choice.

We will shortly need to know the sign of the Mellin convolutions (4.2). If
$F(z_{1},z_{2},\ldots)\geq 0$ for all $z_{1},z_{2}$, then we obviously have
$\widetilde{P}_{qg}(z,\ldots)\otimes F(z,x_{2},\ldots)\geq 0$. On the other
hand, the plus-prescription for $\widetilde{P}_{qq}$ involves a negative term
proportional to $F(x_{1},x_{2},\ldots)$:

$\displaystyle\widetilde{P}_{qq}(z)\underset{x_{1}}{\otimes}F(z,x_{2},\ldots)=C_{F}\lim_{\epsilon\to
0^{+}}\Biggl{[}\mskip 1.5mu$
$\displaystyle\int_{x_{1}/(1-x_{2})}^{1-\epsilon}\frac{dz}{z}\,\frac{1+z^{2}}{1-z}\,F\biggl{(}\frac{x_{1}}{z},x_{2},\ldots\biggr{)}$
$\displaystyle-F(x_{1},x_{2},\ldots)\int_{0}^{1-\epsilon}dz\;\frac{2}{1-z}\mskip
1.5mu\Biggr{]}\,.$ (4.6)

It can therefore have any sign, even if $F(z_{1},z_{2},\ldots)\geq 0$ for all
$z_{1},z_{2}$.

To illustrate this, let us consider the DPDs computed with the LO splitting
formula (2.9). We take the PDFs of the CT14lo PDF set [34], using the LHAPDF
interface [35] via ManeParse [36]. We evaluate the DPD at
$\mu=b_{0}/y=10\operatorname{GeV}$ and verify that the PDFs are positive at
that scale. For the strong coupling, we use the value
$\alpha_{s}(10\operatorname{GeV})=0.178$ provided by the PDF set. As shown in
figure 3, the convolution of $\widetilde{P}_{qq}$ with
${}^{11}\mskip-1.5muF_{q\bar{q}}$ so obtained is indeed negative in a large
region of the momentum fractions.

(a) $x_{2}=x_{1}$

(b) $x_{2}=0.1$

(c) $x_{2}=10^{-3}$

Figure 3:  The convolution
$x_{1}x_{2}\,P(z,\mu)\,\otimes_{x_{1}}\\!{}^{11}\mskip-1.5muF_{q\bar{q}}(z,x_{2},y,\mu,\mu)$
with $P(z,\mu)=a_{s}(\mu)\mskip 1.5mu\widetilde{P}_{qq}(z)$ and
${}^{11}\mskip-1.5muF_{q\bar{q}}$ computed from the LO splitting formula (2.9)
at $\mu=b_{0}/y=10\operatorname{GeV}$. A weighting factor $x_{1}x_{2}$ is
included to keep details visible at higher momentum fractions.

The evolution equations in the $s$ channel basis are derived along the same
lines as the Collins-Soper equation in (3). If the first parton is a quark, we
get

$\displaystyle\frac{\partial}{\partial\log\mu_{1}^{2}}\,\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{qq}^{6\mskip 1.5mu\overline{6}}\end{pmatrix}$
$\displaystyle=a_{s}(\mu_{1})\,\widetilde{P}_{qq}\underset{x_{1}}{\otimes}\widehat{\mathbf{P}}_{qq,q}\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{qq}^{6\mskip
1.5mu\overline{6}}\end{pmatrix}+a_{s}(\mu_{1})\,\widetilde{P}_{qg}\underset{x_{1}}{\otimes}\widehat{\mathbf{P}}_{qg,q}\begin{pmatrix}F_{gq}^{3\mskip
1.5mu\overline{3}}\\\\[3.00003pt] F_{gq}^{\overline{6}\mskip
1.5mu6}\\\\[3.00003pt] F_{gq}^{15\mskip 1.5mu\overline{15}}\end{pmatrix}$
$\displaystyle\quad+\frac{\gamma_{q}(\mu_{1})}{2}\,\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{qq}^{6\mskip
1.5mu\overline{6}}\end{pmatrix}+\frac{{}^{8}\mskip-1.5mu\gamma_{J}(\mu_{1})}{4}\,\log\frac{\mu_{1}^{2}}{x_{1}^{2}\zeta_{p}}\;\hat{\mathbf{J}}_{qq}\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{qq}^{6\mskip 1.5mu\overline{6}}\end{pmatrix}\,,$
(4.7)
$\displaystyle\frac{\partial}{\partial\log\mu_{1}^{2}}\,\begin{pmatrix}F_{q\bar{q}}^{11}\\\\[1.99997pt]
F_{q\bar{q}}^{88}\end{pmatrix}$
$\displaystyle=a_{s}(\mu_{1})\,\widetilde{P}_{qq}\underset{x_{1}}{\otimes}\widehat{\mathbf{P}}_{qq,\bar{q}}\begin{pmatrix}F_{q\bar{q}}^{11}\\\\[1.99997pt]
F_{q\bar{q}}^{88}\end{pmatrix}+a_{s}(\mu_{1})\,\widetilde{P}_{qg}\underset{x_{1}}{\otimes}\widehat{\mathbf{P}}_{qg,\bar{q}}\begin{pmatrix}F_{g\bar{q}}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{g\bar{q}}^{6\mskip 1.5mu\overline{6}}\\\\[3.00003pt]
F_{g\bar{q}}^{\overline{15}\mskip 1.5mu15}\end{pmatrix}$
$\displaystyle\quad+\frac{\gamma_{q}(\mu_{1})}{2}\,\begin{pmatrix}F_{q\bar{q}}^{11}\\\\[1.99997pt]
F_{q\bar{q}}^{88}\end{pmatrix}+\frac{{}^{8}\mskip-1.5mu\gamma_{J}(\mu_{1})}{4}\,\log\frac{\mu_{1}^{2}}{x_{1}^{2}\zeta_{p}}\;\hat{\mathbf{J}}_{q\bar{q}}\begin{pmatrix}F_{q\bar{q}}^{11}\\\\[1.99997pt]
F_{q\bar{q}}^{88}\end{pmatrix}\,,$ (4.8)

where the function arguments of $\widetilde{P}$ and $F$ are as those of $P$
and $F$ in (4.1). The colour mixing matrices read

$\displaystyle\widehat{\mathbf{P}}_{qq,q}$
$\displaystyle=\mathbf{M}_{qq}\,\begin{pmatrix}1&0\\\
0&-1/8\end{pmatrix}\mathbf{M}_{qq}^{-1}=\frac{1}{8}\,\begin{pmatrix}2&6\\\
3&5\end{pmatrix}\,,$ $\displaystyle\widehat{\mathbf{P}}_{qq,\bar{q}}$
$\displaystyle=\mathbf{M}_{q\bar{q}}\,\begin{pmatrix}1&0\\\
0&-1/8\end{pmatrix}\mathbf{M}_{q\bar{q}}^{-1}=\frac{1}{8}\,\begin{pmatrix}0&8\\\
1&7\end{pmatrix}\,,$ (4.9)

and

$\displaystyle\widehat{\mathbf{P}}_{qg,q}$
$\displaystyle=\mathbf{M}_{qq}\,\begin{pmatrix}1&0&0\\\
0&\sqrt{5}/4&3/4\end{pmatrix}\mathbf{M}_{gq}^{-1}=\frac{1}{6}\,\begin{pmatrix}4&12&0\\\
1&0&15\end{pmatrix}\,,$ $\displaystyle\widehat{\mathbf{P}}_{qg,\bar{q}}$
$\displaystyle=\mathbf{M}_{q\bar{q}}\,\begin{pmatrix}1&0&0\\\
0&\sqrt{5}/4&3/4\end{pmatrix}\mathbf{M}_{g\bar{q}}^{-1}=\frac{1}{24}\,\begin{pmatrix}64&0&0\\\
1&18&45\end{pmatrix}\,,$ (4.10)

whilst $\hat{\mathbf{J}}_{qq}$ and $\hat{\mathbf{J}}_{q\bar{q}}$ are given in
(3.7) and (3), respectively. With

$\displaystyle\widehat{\mathbf{P}}_{\bar{q}g,\bar{q}}$
$\displaystyle=\mathbf{M}_{\bar{q}\bar{q}}\,\begin{pmatrix}1&0&0\\\
0&\sqrt{5}/4&-3/4\end{pmatrix}\mathbf{M}_{g\bar{q}}^{-1}=\widehat{\mathbf{P}}_{qg,q}\,,$
$\displaystyle\widehat{\mathbf{P}}_{\bar{q}g,q}$
$\displaystyle=\mathbf{M}_{\bar{q}q}\,\begin{pmatrix}1&0&0\\\
0&\sqrt{5}/4&-3/4\end{pmatrix}\mathbf{M}_{gq}^{-1}=\widehat{\mathbf{P}}_{qg,\bar{q}}\,,$
(4.11)

we find that the evolution equations for $F_{\bar{q}\bar{q}}$ and
$F_{\bar{q}q}$ are respectively obtained from (4.1) and (4.1) by interchanging
$q\leftrightarrow\bar{q}$ in the DPDs and swapping their representation
labels, i.e. $F_{qq}^{\overline{3}\mskip 1.5mu3}\to
F_{\bar{q}\bar{q}}^{3\mskip 1.5mu\overline{3}}$, $F_{gq}^{3\mskip
1.5mu\overline{3}}\to F_{g\bar{q}}^{\overline{3}\mskip 1.5mu3}$, etc. The
splitting kernels, anomalous dimensions, and colour mixing matrices remain the
same.

The evolution equations have the same form for distributions with unequal
flavours, i.e. one may replace $F_{qq}\to F_{q_{1}q_{2}}$ and $F_{gq}\to
F_{gq_{2}}$ in (4.1), or $F_{q\bar{q}}\to F_{q_{1}\bar{q}_{2}}$ and
$F_{g\bar{q}}\to F_{g\bar{q}_{2}}$ in (4.1), whilst keeping the splitting
kernels and colour mixing matrices unchanged.

Before analysing the effect of the evolution equations (4.1) and (4.1) on
positivity, let us recall the situation for colour singlet distributions in
the $t$ channel basis. The LO evolution equation for
${}^{11}\mskip-1.5muF_{qa}$ reads

$\displaystyle\frac{\partial}{\partial\log\mu_{1}^{2}}\,{}^{11}\mskip-1.5muF_{qa}(x_{1},x_{2},\ldots)$
$\displaystyle=a_{s}(\mu_{1})\,\widetilde{P}_{qq}\underset{x_{1}}{\otimes}{}^{11}\mskip-1.5muF_{qa}+a_{s}(\mu_{1})\,\widetilde{P}_{qg}\underset{x_{1}}{\otimes}{}^{11}\mskip-1.5muF_{ga}+\frac{\gamma_{q}(\mu_{1})}{2}\,{}^{11}\mskip-1.5muF_{qa}\,,$
(4.12)

where $a$ is a quark or an antiquark. If ${}^{11}\mskip-1.5muF_{qa}$ and
${}^{11}\mskip-1.5muF_{ga}$ are non-negative for all momentum fractions, the
terms involving $\widetilde{P}_{qg}$ or $\gamma_{q}$ in (4.12) are non-
negative as well and hence conserve positivity. As we have shown in figure 3,
the first term in (4.12) can be negative. However, the part that is
responsible for a negative sign is proportional to
${}^{11}\mskip-1.5muF_{qa}(x_{1},x_{2},\ldots)$ itself, as is easily seen in
(4.1). This negative contribution hence decreases in magnitude when
${}^{11}\mskip-1.5muF_{qa}(x_{1},x_{2},\ldots)$ approaches zero from above,
and closer inspection shows that it cannot lead to a violation of positivity.
This is shown in more detail in appendix B of [4]. Overall, positivity is
hence conserved for LO evolution of ${}^{11}\mskip-1.5muF_{qa}$ to higher
scales, and the same can be shown for all other parton combinations.

We now analyse the sign of the different terms on the r.h.s. of (4.1) and
(4.1) under the assumption that the $s$ channel distributions on the r.h.s.
are non-negative for all momentum fractions. For brevity, we write
$D_{1}F=\partial F/\partial\log\mu_{1}^{2}$.

  1. 1.

The scale variation $D_{1}F_{qa}^{R\overline{R}}$ receives contributions
$a_{R}\mskip 1.5mu\widetilde{P}_{qq}\otimes F_{qa}^{R\overline{R}}$ and
$b_{R}\mskip 1.5mu\widetilde{P}_{qq}\otimes F_{qa}^{R^{\mskip
1.5mu\prime}\mskip-1.5mu\overline{R}{}^{\mskip 1.5mu\prime}}$ with
coefficients $a_{R}\geq 0$, $b_{R}>0$, where $R^{\mskip 1.5mu\prime}\neq R$. A
negative contribution from the term with $a_{R}$ cannot lead to a violation of
positivity, as just discussed. However, the contribution from the term with
$b_{R}$ can remain large and negative even if $F_{qa}^{R\overline{R}}$
approaches zero. This term can therefore lead to a zero crossing of the
distribution as one evolves to higher scales.

  2. 2.

The matrices $\widehat{\mathbf{P}}_{qg,a}$ have no negative elements, so that
the terms with $\widetilde{P}_{qg}$ are all non-negative.

  3. 3.

The terms with $\gamma_{q}$ are all non-negative.

  4. 4.

The contributions with ${}^{8}\mskip-1.5mu\gamma_{J}$ to
$D_{1}F_{qq}^{\overline{3}3}$ and $D_{1}F_{qq}^{6\overline{6}}$ have opposite
sign, as well as those to $D_{1}F_{q\bar{q}}^{11}$ and
$D_{1}F_{q\bar{q}}^{88}$. This follows from the form of
$\hat{\mathbf{J}}_{qq}$ and $\hat{\mathbf{J}}_{q\bar{q}}$ in (3.7) and (3).
Whether these contributions can violate positivity depends on the sign of
$\log\mu_{1}^{2}/(x_{1}^{2}\mskip 1.5mu\zeta_{p})$.

It follows that LO DGLAP evolution of DPDs is not guaranteed to conserve their
positivity in colour space, in contrast to the colour singlet distributions
${}^{11}\mskip-1.5muF_{a_{1}a_{2}}$. This is one of our main results.

As a numerical illustration, let us take the initial conditions provided by
the perturbative splitting mechanism at LO. According to (2.12), we have
$F_{q\bar{q}}^{11}=F_{g\bar{q}}^{6\overline{6}}=F_{g\bar{q}}^{\overline{15}\mskip
1.5mu15}=0$ at the scale $\mu$ where the splitting formula is evaluated.
Inserting this in the evolution equation (4.1) and taking the LO approximation
(3.4) of ${}^{8}\mskip-1.5mu\gamma_{J}$, we obtain

$\displaystyle\frac{\partial}{\partial\log\mu_{1}^{2}}\,F_{q\bar{q}}^{11}(\mu_{1},\mu_{2})=a_{s}(\mu)\Biggl{[}\mskip
1.5mu$
$\displaystyle\widetilde{P}_{qq}\underset{x_{1}}{\otimes}F_{q\bar{q}}^{88}(\mu,\mu)+\frac{8}{3}\,\widetilde{P}_{qg}\underset{x_{1}}{\otimes}F_{g\bar{q}}^{\overline{3}\mskip
1.5mu3}(\mu,\mu)$
$\displaystyle-\frac{2}{9}\;{}^{8}\mskip-1.5mu\gamma_{J}^{(0)}\mskip
1.5mu\log\frac{\mu^{2}}{x_{1}^{2}\zeta_{p}}\,F_{q\bar{q}}^{88}(\mu,\mu)\mskip
1.5mu\Biggr{]}$ (4.13)

at the point $\mu_{1}=\mu_{2}=\mu$. We evaluate the r.h.s. numerically with
the same settings as in figure 3 and the choice

$\displaystyle\zeta_{p}=\mu^{2}/(x_{1}x_{2})\,,$ (4.14)

which is natural for the perturbative splitting mechanism [20]. With this
choice, the term with $\zeta_{p}$ is zero at $x_{1}=x_{2}$. At that point, a
negative value of (4.1) must hence be due to the term with
$\widetilde{P}_{qq}$. We see in figure 4 that there are regions in $x_{1}$ and
$x_{2}$ for which $D_{1}F_{q\bar{q}}^{11}<0$ at the scale $\mu$. With
$F_{q\bar{q}}^{11}=0$ at that scale, positivity is thus explicitly violated by
evolution to a higher scale for the first parton.

(a) $x_{2}=x_{1}$

(b) $x_{2}=0.1$

(c) $x_{2}=10^{-3}$

Figure 4:  The right-hand side of (4.1) with $q\bar{q}=u\bar{u}$, evaluated at
$\mu=b_{0}/y=10\operatorname{GeV}$ and $\zeta_{p}=\mu^{2}/(x_{1}x_{2})$ with
LO splitting DPDs. The individual terms going with $\widetilde{P}_{qq}$,
$\widetilde{P}_{qg}$, and $\gamma_{J}$ are shown as well. The PDFs used in the
splitting formula are specified in the text below (4.1). A weighting factor
$x_{1}x_{2}$ is included to keep details visible at higher momentum fractions.

###  4.2 Simultaneous evolution in all scales

When computing double parton scattering cross sections, the choice of scales
$\mu_{1}$, $\mu_{2}$ and $\zeta_{p}$ is driven by the kinematics of the
process. In particular, taking $\mu_{1}\neq\mu_{2}$ is natural for processes
with two hard scales of very different size. On the other hand, the simplest
setting for the physical interpretation of a DPD is with all relevant scales
set equal. We therefore investigate the evolution of DPDs with a common
renormalisation scale $\mu=\mu_{1}=\mu_{2}$ for both partons and the rapidity
scale given by $\zeta_{p}=\mu^{2}/(x_{1}x_{2})$ as in (4.14).

Let us briefly comment on the factor $x_{1}x_{2}$ in (4.14). As explained in
[21, 20], the definition of $\zeta_{p}$ involves the rapidity regulator and
the plus-momentum of the target proton. By contrast, $x_{1}x_{2}\mskip
1.5mu\zeta_{p}$ refers to the regulator and the plus-momenta of the two
extracted partons, and in this sense is more closely related to the scales
$\mu_{1}$ and $\mu_{2}$ that refer to the renormalisation of the operators
associated with the two partons. This motivates our choice (4.14), along with
the fact that $x_{1}x_{2}\mskip 1.5mu\zeta_{p}$ is the combination of
variables appearing at higher orders in the perturbative splitting formula for
DPDs, see (5) and (5.2).

Combining the DGLAP equations for the first and second parton with the
Collins-Soper equation (3) yields

$\displaystyle\frac{\partial}{\partial\log\mu^{2}}\,\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{qq}^{6\mskip 1.5mu\overline{6}}\end{pmatrix}$
$\displaystyle=a_{s}(\mu)\,\widetilde{P}_{qq}\underset{x_{1}}{\otimes}\widehat{\mathbf{P}}_{qq,q}\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{qq}^{6\mskip
1.5mu\overline{6}}\end{pmatrix}+a_{s}(\mu)\,\widetilde{P}_{qq}\underset{x_{2}}{\otimes}\widehat{\mathbf{P}}_{qq,q}\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{qq}^{6\mskip 1.5mu\overline{6}}\end{pmatrix}$
$\displaystyle\quad+a_{s}(\mu)\,\widetilde{P}_{qg}\underset{x_{1}}{\otimes}\widehat{\mathbf{P}}_{qg,q}\begin{pmatrix}F_{gq}^{3\mskip
1.5mu\overline{3}}\\\\[3.00003pt] F_{gq}^{\overline{6}\mskip
1.5mu6}\\\\[3.00003pt] F_{gq}^{15\mskip
1.5mu\overline{15}}\end{pmatrix}+a_{s}(\mu)\,\widetilde{P}_{qg}\underset{x_{2}}{\otimes}\widehat{\mathbf{P}}_{qg,q}\begin{pmatrix}F_{qg}^{3\mskip
1.5mu\overline{3}}\\\\[3.00003pt] F_{qg}^{\overline{6}\mskip
1.5mu6}\\\\[3.00003pt] F_{qg}^{15\mskip 1.5mu\overline{15}}\end{pmatrix}$
$\displaystyle\quad+\gamma_{q}(\mu)\,\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{qq}^{6\mskip
1.5mu\overline{6}}\end{pmatrix}+\frac{{}^{8}\mskip-1.5muJ(y,\mu,\mu)}{2}\;\hat{\mathbf{J}}_{qq}\,\begin{pmatrix}F_{qq}^{\overline{3}\mskip
1.5mu3}\\\\[3.00003pt] F_{qq}^{6\mskip 1.5mu\overline{6}}\end{pmatrix}\,.$
(4.15)

With $\zeta_{p}$ chosen as in (4.14), the term with
${}^{8}\mskip-1.5mu\gamma_{J}$ in (4.1) has cancelled against the
corresponding term from the evolution equation in $\mu_{2}$.

The sign of the r.h.s. of (4.2) can be analysed along the same lines as in the
previous subsection:

  1. 1.

The terms with $\widetilde{P}_{qq}$ can lead to a violation of positivity
because of the negative part in the plus-distribution and of the positive off-
diagonal entries in the matrix $\widehat{\mathbf{P}}_{qq,q}$.

  2. 2.

The terms with $\widetilde{P}_{qg}$ conserve positivity.

  3. 3.

The terms with $\gamma_{q}$ conserve positivity.

  4. 4.

If ${}^{8}\mskip-1.5muJ<0$ then the corresponding terms do not violate
positivity, because negative contributions are proportional to the DPD being
evolved, whilst the DPD in the other colour channel gives a positive
contribution. This is consistent with our findings for Collins-Soper evolution
in section 3.

The leading order expression (3.5) of ${}^{8}\mskip-1.5muJ(y,\mu,\mu)$
contains an explicit logarithm $\log(\mu^{2}y^{2})$. In the leading double
logarithmic approximation, one therefore has to keep only the last term in the
evolution equation (4.2), which then reduces to the LO Collins-Soper equation
with $\mu^{2}=x_{1}x_{2}\mskip 1.5mu\zeta_{p}$.

The evolution equation for $F_{q\bar{q}}$ has the same form as (4.2) and
involves the colour mixing matrices $\widehat{\mathbf{P}}_{qq,\bar{q}}$,
$\widehat{\mathbf{P}}_{qg,\bar{q}}$, and $\hat{\mathbf{J}}_{q\bar{q}}$. Its
discussion proceeds in full analogy.

In summary, we find that evolution to higher scales can violate positivity of
DPDs in colour space, both when one evolves in the renormalisation scale of
one parton and when one evolves in all scales simultaneously.

##  5 DPDs from parton splitting at two-loop accuracy

In section 2 we saw that the perturbative splitting mechanism gives DPDs that
satisfy positivity if the splitting is computed at LO, i.e. at one-loop
accuracy. This is not surprising, since the LO splitting formula in the $s$
channel basis can be written as a squared matrix element (in the mixed
representation of definite plus-momentum and transverse position for the
observed partons). Starting from two loops, the splitting formula has explicit
logarithms of the renormalisation scale $\mu$ and of the rapidity parameter
$\zeta_{p}$, which respectively result from subtractions for ultraviolet and
rapidity divergences. As a consequence of these subtractions, positivity is no
longer guaranteed. It is then natural to ask whether the resulting DPDs at
small $y$ violate positivity in colour space, and if so, by how much. We
address this question in the present section, using the results of the recent
two-loop calculation in [20].

The generalisation of the LO splitting formula (2.9) to higher orders has the
form

$\displaystyle{}^{R_{1}R_{2}}\mskip-1.5muF_{a_{1}a_{2}}(x_{1},x_{2},y,\mu,\mu,\zeta_{p})$
$\displaystyle\quad=\frac{1}{\pi\mskip
1.5muy^{2}}\,\sum_{a_{0}}\,{}^{R_{1}R_{2}}V_{a_{1}a_{2},a_{0}}\bigl{(}z_{1},z_{2},a_{s}(\mu),L_{y}(\mu),L_{\zeta}(\mu)\bigr{)}\underset{x}{\otimes}f_{a_{0}}(z,\mu)$
(5.1)

with

$\displaystyle L_{y}(\mu)$ $\displaystyle=\log\frac{\mu^{2}}{\mu_{y}^{2}}\,,$
$\displaystyle\mu_{y}$ $\displaystyle=\frac{b_{0}}{y}\,,$ $\displaystyle
L_{\zeta}(\mu)$ $\displaystyle=\log\frac{\mu^{2}}{x_{1}x_{2}\mskip
1.5mu\zeta_{p}}$ (5.2)

and the special convolution

$\displaystyle V(z_{1},z_{2})\underset{x}{\otimes}f(z)$
$\displaystyle=\frac{1}{x}\int_{x}^{1}dz\;V\bigl{(}uz,(1-u)z\mskip
1.5mu\bigr{)}\,f\biggl{(}\frac{x}{z}\biggr{)}\,.$ (5.3)

Here we write

$\displaystyle z_{1}$ $\displaystyle=uz\,,$ $\displaystyle z_{2}$
$\displaystyle=(1-u)\mskip 1.5muz$ (5.4) and recall that $\displaystyle x_{1}$
$\displaystyle=ux\,,$ $\displaystyle x_{2}$ $\displaystyle=(1-u)\mskip 1.5mux$
(5.5)

according to (2.11). The analogue of (5) for distributions in the $s$ channel
basis is readily obtained using the transformations in section 2. We limit our
attention to the quark-antiquark sector and consider the distributions

$\displaystyle F_{ud}^{\overline{3}3}$ $\displaystyle=4\mskip
1.5muF_{ud}^{6\overline{6}}\,,$ $\displaystyle F_{u\bar{d}}^{11}$
$\displaystyle=64\mskip 1.5muF_{u\bar{d}}^{88}$ (5.6)

and

$\displaystyle F_{uu}^{\overline{3}3}\,,\qquad F_{uu}^{6\overline{6}}\,,\qquad
F_{u\bar{u}}^{11}\,,\qquad F_{u\bar{u}}^{88}\,.$ (5.7)

Among these, only $F_{u\bar{u}}^{88}$ is nonzero at order $a_{s}$, whilst all
others start at order $a_{s}^{2}$.

#### Parton combinations appearing first at two loops.

At order $a_{s}^{2}$ the distributions $F_{ud}$ and $F_{u\smash{\bar{d}}}$
receive a contribution only from the kernels $V_{qq^{\prime}\\!,q}$ or
$V_{q\bar{q}^{\prime}\\!,q}$, which respectively correspond to the graphs in
figure 5(a) and 5(d) and further graphs with identical topology. As a
consequence, distributions for different colour representations are
proportional to each other, as specified in equation (4.40) of [20]. This
leads to the relations in (5.6).

Furthermore, the splitting kernels for $F_{ud}$ and $F_{u\smash{\bar{d}}}$ are
proportional to each other, because
${}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}={}^{11}\mskip-1.5muV_{q\bar{q}^{\prime}\\!,q}$
at two-loop accuracy. Using the basis transform (2.5) and the proportionality
factors between $t$ channel singlet and octet kernels in equation (4.35) of
[20], one obtains

$\displaystyle F_{ud}^{\overline{3}3}$ $\displaystyle=\frac{1}{\pi\mskip
1.5muy^{2}}\Bigl{[}\mskip
1.5muV_{qq^{\prime}\\!,q}^{\overline{3}3}(z_{1},z_{2})\underset{x}{\otimes}f_{u}(z)+V_{qq^{\prime}\\!,q}^{\overline{3}3}(z_{2},z_{1})\underset{x}{\otimes}f_{d}(z)\mskip
1.5mu\Bigr{]}\,,$ $\displaystyle F_{u\bar{d}}^{11}$
$\displaystyle=\frac{1}{\pi\mskip 1.5muy^{2}}\Bigl{[}\mskip
1.5muV_{q\bar{q}^{\prime}\\!,q}^{11}(z_{1},z_{2})\underset{x}{\otimes}f_{u}(z)+V_{q\bar{q}^{\prime}\\!,q}^{11}(z_{2},z_{1})\underset{x}{\otimes}f_{\bar{d}}(z)\mskip
1.5mu\Bigr{]}$ (5.8)

with

$\displaystyle V_{qq^{\prime}\\!,q}^{\overline{3}3}$
$\displaystyle=\frac{2}{9}\,{}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}\,,$
$\displaystyle V_{q\bar{q}^{\prime}\\!,q}^{11}$
$\displaystyle=\frac{8}{9}\,{}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}\,.$ (5.9)

Here the convolution $V(z_{2},z_{1})\otimes f(z)$ is defined as in (5.3) with
$u$ and $1-u$ interchanged on the r.h.s. In the last term of (5) we used the
relation $V_{\bar{q}q^{\prime}\\!,\bar{q}}=V_{q\bar{q}^{\prime}\\!,q}$ from
charge conjugation invariance. Notice that the kernel
${}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}$ includes an overall factor
$a_{s}^{2}$.

(a) $V_{qq^{\prime}\\!,q}$ (b) $V_{qq,q}^{v}$

(c) $V_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}$ (d) $V_{q\bar{q}^{\prime}\\!,q}$

(e) $V_{q\bar{q},q}^{v}$

Figure 5:  Real graphs for the splitting $a_{0}\to a_{1}a_{2}$ in the pure
quark channels. The parton lines at the bottom of the graph correspond to
$a_{0}$, and those on top to $a_{1}$, $a_{2}$, $a_{2}$ and $a_{1}$ from left
to right. Additional graphs are obtained by interchanging
$a_{1}\leftrightarrow a_{2}$ and by reflection w.r.t. the final state cut.

The distributions $F_{uu}^{\overline{3}3}$ and $F_{uu}^{6\overline{6}}$
receive contributions from the kernels $V_{qq^{\prime}\\!,q}$ and
$V_{qq,q}^{v}$ with different weights and are therefore not proportional to
each other. Using equation (3.3) in [20], we find that the relevant kernels
read

$\displaystyle V_{qq,q}^{\overline{3}3}(z_{1},z_{2})$
$\displaystyle=\frac{4}{18}\,\Bigl{[}{}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}(z_{1},z_{2})+{}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}(z_{2},z_{1})\Bigr{]}+\frac{4}{6}\,{}^{11}\mskip-1.5muV_{qq,q}^{v}(z_{1},z_{2})\,,$
$\displaystyle V_{qq,q}^{6\overline{6}}(z_{1},z_{2})$
$\displaystyle=\frac{1}{18}\,\Bigl{[}{}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}(z_{1},z_{2})+{}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}(z_{2},z_{1})\Bigr{]}-\frac{1}{6}\,{}^{11}\mskip-1.5muV_{qq,q}^{v}(z_{1},z_{2})$
(5.10)

and are to be convolved with $f_{u}(z)$.

At order $a_{s}^{2}$, the kernel ${}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}$ and
hence all kernels in (5.9) and (5) require ultraviolet renormalisation but no
subtraction of rapidity divergences. As a consequence, they depend linearly on
the renormalisation group logarithm $L_{y}$. A natural choice of scale in the
fixed-order formula (5) is $\mu=\mu_{y}$, so that $L_{y}=0$. In a numerical
study of the two-gluon distributions ${}^{R\overline{R}}\mskip-1.5muF_{gg}$ in
[20], we found that the size of $a_{s}^{2}$ corrections relative to the
$a_{s}$ term is moderate for $\mu=\mu_{y}$ but grows substantially for
$\mu=\mu_{y}/2$ or $\mu=2\mu_{y}$. In the present work, we therefore consider
a smaller amount of variation and take $\mu=1.2\mskip 1.5mu\mu_{y}$ or
$\mu=\mu_{y}/1.2$ as alternative scales, which corresponds to $L_{y}\approx\pm
0.36$. We checked that the size of $a_{s}^{2}$ corrections in the two-gluon
sector remains moderate for these values. The DPDs depend on $\mu$ as
specified by the relevant DGLAP equations. When comparing the fixed-order
formula (5) for different $\mu$, one thus sees evolution effects truncated at
the lowest order in $a_{s}$.

(a) ${}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}$ at $\mu=\mu_{y}$ (b)
${}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}$ at $\mu=1.2\mskip 1.5mu\mu_{y}$

Figure 6:  The kernel ${}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}(z_{1},z_{2})$,
which appears in the DPDs for the parton combinations $ud$ and $u\bar{d}$ as
specified in (5.6), (5), and (5.9). A global factor $a_{s}^{2}$ has been
divided out. See (5.4) for the relation between $z_{1},z_{2}$ and $u,z$.

(a) $V_{qq,q}^{\overline{3}3}$ for $\mu=\mu_{y}$ (b)
$V_{qq,q}^{\overline{3}3}$ for $\mu=1.2\mskip 1.5mu\mu_{y}$

(c) $V_{qq,q}^{6\overline{6}}$ for $\mu=\mu_{y}$ (d)
$V_{qq,q}^{6\overline{6}}$ for $\mu=1.2\mskip 1.5mu\mu_{y}$

Figure 7:  As figure 7 but for the kernels in (5), which give the DPDs for the
parton combination $uu$. Since these kernels are symmetric under the
interchange $z_{1}\leftrightarrow z_{2}$ of the two momentum fractions, we can
limit ourselves to the region $u\leq 0.5$

In figures 7 and 7, we plot the kernels
${}^{11}\mskip-1.5muV_{qq^{\prime}\\!,q}$, $V_{qq,q}^{\overline{3}3}$, and
$V_{qq,q}^{6\overline{6}}$. We see that they are negative over wide ranges of
$z$ and $u$, although they are obtained from a sum of graphs that correspond
to the squared amplitude for $q\to qq^{\prime}\mskip-1.5mu\bar{q}^{\prime}$ or
$q\to qq\bar{q}$. This illustrates that the subtraction of ultraviolet
divergences can indeed lead to a negative result. We note that the kernels go
to large positive values for $z\to 1$, which can be traced back to terms
diverging like $\log 1/(1-z)$ in that limit. Finally, we observe that the
kernels increase with $\mu$ for all $z$ and $u$.

We now investigate to which extent the negative regions in the splitting
kernels lead to negative DPDs. To this end, we evaluate the splitting formula
for the distributions in (5.6) and (5.7) with the central PDFs from the
CT14nlo set, having checked that the PDFs are positive in the kinematics of
interest. As we did in section 4, we fix $y$ such that
$\mu_{y}=10\operatorname{GeV}$.

(a) $x_{2}=x_{1}$, $\mu=\mu_{y}/1.2$

(b) $x_{2}=x_{1}$, $\mu=\mu_{y}$

(c) $x_{2}=0.01$, $\mu=\mu_{y}/1.2$

(d) $x_{2}=0.01$, $\mu=\mu_{y}$

Figure 8:  The distributions $F_{ud}^{\overline{3}3}$ and $F_{u\bar{d}}^{11}$
computed with the two-loop splitting formula for
$\mu_{y}=10\operatorname{GeV}$ and two choices of $\mu$. The other two colour
channels are given in (5.6).

(a) $x_{2}=x_{1}$, $\mu=\mu_{y}$

(b) $x_{2}=x_{1}$, $\mu=1.2\mskip 1.5mu\mu_{y}$

(c) $x_{2}=0.01$, $\mu=\mu_{y}$

(d) $x_{2}=0.01$, $\mu=1.2\mskip 1.5mu\mu_{y}$

Figure 9:  The distributions $F_{uu}^{\overline{3}3}$ and
$F_{uu}^{6\overline{6}}$ computed with the two-loop splitting formula for
$\mu_{y}=10\operatorname{GeV}$. Note that the two values of $\mu$ are not the
same as in figure 8.

In figure 8 we show $F_{ud}^{\overline{3}3}$ and $F_{u\bar{d}}^{11}$ for
$\mu=\mu_{y}/1.2$ and $\mu=\mu_{y}$. The curves are scaled such that the
difference between them originates from the different PDFs in (5) and not from
the different normalisation of the kernels in (5.9). We observe a strong
effect of scale evolution: for $x_{1}\sim x_{2}$ the distributions are
negative at the smaller scale but positive at $\mu=\mu_{y}$ (and also at
$\mu=1.2\mskip 1.5mu\mu_{y}$, which is not shown in the figure). We also note
that the negative values in figure 8(c) are tiny compared with the size of the
same distributions at other values of $x_{1}$ for the same $x_{2}$.

The distributions $F_{uu}^{\overline{3}3}$ and $F_{uu}^{6\overline{6}}$ are
shown in figure 9, with a scaling factor such that the difference between the
curves is due to the contribution of ${}^{11}\mskip-1.5muV_{qq,q}^{v}$ to the
kernels in (5). The situation for $\mu=\mu_{y}/1.2$ (not shown in the figure)
is qualitatively similar to the one at $\mu=\mu_{y}$, where we find negative
values at $x_{1}\sim x_{2}$ for $F_{uu}^{\overline{3}3}$ but not for
$F_{uu}^{6\overline{6}}$. At the higher scale $\mu=1.2\mskip 1.5mu\mu_{y}$,
all values are positive. As in figure 8(c), the negative values in figure 9(c)
are tiny compared with the size of the distribution at other momentum
fractions. In this sense, the violations of positivity we have shown so far
may be regarded as minor.

#### Quark-antiquark distributions.

The last two distributions in (5.7) are for a quark and an antiquark of equal
flavour. $F_{u\bar{u}}^{11}$ receives contributions from all three kernels in
the bottom row of figure 5 and from the real two-loop graphs for the splitting
$g\to q\bar{q}$, an example of which is shown in figure 2(b).
$F_{u\bar{u}}^{88}$ receives contributions from the same graphs, from the LO
graph in figure 2(a), and from virtual two-loop graphs such as the one in
figure 2(c). The virtual graphs depend on the colour of the observed partons
in the same way as the LO graph and hence do not contribute to
$F_{u\bar{u}}^{11}$.

The distributions $F_{u\bar{u}}^{11}$ and $F_{u\bar{u}}^{88}$ require both
ultraviolet renormalisation and the subtraction of rapidity divergences. The
latter appears in the splitting process $g\to q\bar{q}$, whose kernels have a
more complicated structure than the ones considered so far. Up to order
$a_{s}^{2}$, they can be written as

$\displaystyle V_{q\bar{q},g}^{11}$
$\displaystyle=a_{s}^{2}\,D\;{}^{8}\mskip-1.5mu\gamma_{J}^{(0)}\,\delta(1-z)\,{}^{11}\mskip-1.5muV_{q\bar{q},g}^{(1)}(u)+a_{s}^{2}\,\Bigl{(}V_{q\bar{q},g}^{11\,[2,0]}(z,u)+L_{y}\,V_{q\bar{q},g}^{11\,[2,1]}(z,u)\Bigr{)}\,,$
$\displaystyle V_{q\bar{q},g}^{88}$
$\displaystyle=\frac{a_{s}}{8}\,\Bigl{(}1-a_{s}\mskip
1.5muD\;{}^{8}\mskip-1.5mu\gamma_{J}^{(0)}\Bigr{)}\,\delta(1-z)\,{}^{11}\mskip-1.5muV_{q\bar{q},g}^{(1)}(u)+a_{s}^{2}\,\Bigl{(}V_{q\bar{q},g}^{88\,[2,0]}(z,u)+L_{y}\,V_{q\bar{q},g}^{88\,[2,1]}(z,u)\Bigr{)}\,,$

where
${}^{11}\mskip-1.5muV_{q\bar{q},g}^{(1)}(u)=\bigl{[}u^{2}+(1-u)^{2}\mskip
1.5mu\bigr{]}\big{/}2$ is the LO splitting kernel appearing in (2.9) and

$\displaystyle D$ $\displaystyle=\frac{1}{36}\mskip
1.5mu\biggl{(}L_{y}^{2}-2L_{y}\mskip
1.5muL_{\zeta}+\frac{\pi^{2}}{6}\biggr{)}$ (5.12)

contains the double logarithms associated with rapidity divergences. We note
that (5.12) is obtained with the standard definition of the
$\overline{\text{MS}}$ scheme, and that the term $\pi^{2}/6$ is absent if one
instead uses the definition proposed by Collins in section 3.2.6 of [22]. The
coefficients $V^{[2,0]}$ and $V^{[2,1]}$ in (5) are smooth functions of $u$
but distributions in $z$. They consist of a regular part, a part proportional
to the plus-distribution $1/[1-z]_{+}$, and a term proportional to
$\delta(1-z)$. The regular part is a smooth function of $z$ but may have a
$\log(1-z)$ singularity for $z\to 1$, similar to what we saw for the pure
quark kernels in figures 7 and 7.

In figure 10 we show the different contributions to $F_{u\bar{u}}^{11}$ for
$x_{1}=x_{2}$. The regular part of the $g\to q\bar{q}$ kernel turns out to be
negative for all values of $z$ and $u$ and results in a large negative
contribution to the DPD. At $\mu=1.2\mskip 1.5mu\mu_{y}$, another negative
contribution comes from the part of the kernel that goes with the plus-
distribution $1/[1-z]_{+}$.

(a) $\mu=\mu_{y}$ (b) $\mu=1.2\mskip 1.5mu\mu_{y}$

Figure 10:  Individual contributions to $F_{u\bar{u}}^{11}$ at $x_{1}=x_{2}$
at two-loop level, evaluated with $x_{1}x_{2}\mskip 1.5mu\zeta_{p}=\mu^{2}$
for two values of $\mu$. The curves labelled “$g$” are for the splitting $g\to
q\bar{q}$, and the curve labelled “$q$” is for the sum of all splitting
contributions initiated by a quark or an antiquark. The labels “regular”,
“$[1-z]_{+}$”, and “$\delta(1-x)$” refer to the different parts of the kernels
$V^{11\,[2,0]}+L_{y}\,V^{11\,[2,1]}$ discussed below equation (5.12).

Adding up all contributions, one obtains the distribution shown in the two
upper rows of figure 11. The negative contributions dominate for $x_{1}=x_{2}$
up to a few $0.01$, whereas for larger momentum fractions the positive
contributions from quark or antiquark splitting gradually takes over. In the
two lower rows of figure 11, we see that the regions of negative
$F_{u\bar{u}}^{11}$ are centred around $x_{1}\sim x_{2}$, with values that are
small compared with the size of the distribution at other values of $x_{1}$ at
the same $x_{2}$. This is the same phenomenon that we observed earlier for
$F_{ud}$, $F_{u\smash{\bar{d}}}$, and $F_{uu}^{\overline{3}3}$. However, in
the present case the negative values around $x_{1}\sim x_{2}$ decrease with
$\mu$, so that the violation of positivity in $F_{u\bar{u}}^{11}$ becomes more
pronounced as $\mu$ becomes larger.

In figure 11 we also see that a mild variation of the rapidity parameter
$\zeta_{p}$ has a rather small effect in the kinematics considered here. The
choice of $\zeta_{p}$ is only relevant if $\mu\neq\mu_{y}$, because
$L_{\zeta}$ is multiplied by $L_{y}$ in (5.12).

Let us finally turn to the distribution $F_{u\bar{u}}^{88}$. Evaluated at one-
loop accuracy, this distribution is positive, so that negative values at two-
loop level can only appear when the contribution of order $a_{s}^{2}$ is
larger in size than the one of order $a_{s}$. In such a case, one may of
course worry whether the unknown contributions of yet higher orders will
change the sign of the distribution again. This situation is qualitatively
different from the one for the distributions discussed so far, where the terms
of order $a_{s}^{2}$ give the first non-vanishing contribution.

(a) $x_{2}=x_{1}$, $\mu=\mu_{y}$

(b) $x_{2}=x_{1}$, $\mu=1.2\mskip 1.5mu\mu_{y}$

(c) $x_{2}=0.01$, $\mu=\mu_{y}$

(d) $x_{2}=0.01$, $\mu=1.2\mskip 1.5mu\mu_{y}$

Figure 11:  The distribution $F_{u\bar{u}}^{11}$ computed with the two-loop
splitting formula for $\mu_{y}=10\operatorname{GeV}$ and two choices of $\mu$.

(a) $x_{2}=x_{1}$, $\mu=\mu_{y}$

(b) $x_{2}=x_{1}$, $\mu=1.2\mskip 1.5mu\mu_{y}$

(c) $x_{2}=0.01$, $\mu=\mu_{y}$

(d) $x_{2}=0.01$, $\mu=1.2\mskip 1.5mu\mu_{y}$

Figure 12:  The distribution $F_{u\bar{u}}^{88}$ computed for
$\mu_{y}=10\operatorname{GeV}$, either with one-loop (LO) or with the sum of
one- and two-loop splitting kernels (LO+NLO). All plots are for
$x_{1}x_{2}\mskip 1.5mu\zeta_{p}=\mu^{2}$.

(a) $x_{1}=0.1$, $\mu=1.2\mskip 1.5mu\mu_{y}$

(b) $x_{2}=10^{-3}$, $\mu=1.2\mskip 1.5mu\mu_{y}$

(c) $x_{1}=10^{-3}$, $\mu=\mu_{y}/1.2$ (d) detail of (c)

Figure 13:  Continuation of figure 12. Notice that $x_{1}$ and not $x_{2}$ is
fixed in panels (a), (c), and (d).

In figures 12 and 13 we show $F_{u\bar{u}}^{88}$ for different values of the
momentum fractions. We always take $x_{1}x_{2}\mskip 1.5mu\zeta_{p}=\mu^{2}$,
bearing in mind that according to (5) the effect of varying $\zeta_{p}$ is
eight times smaller for $F_{u\bar{u}}^{88}$ than it is for
$F_{u\bar{u}}^{11}$. As can be seen in figure 12, the difference between the
distributions at $\mu=\mu_{y}$ to $\mu=1.2\mskip 1.5mu\mu_{y}$ is rather
small, in contrast to what we found for the other distributions discussed so
far. This is not surprising, because for $F_{u\bar{u}}^{88}$ the relative
effect of changing the scale from $\mu_{a}$ to $\mu_{b}$ is of order
$a_{s}\log(\mu_{a}/\mu_{b})$, whereas it is of order $\log(\mu_{a}/\mu_{b})$
for distributions that receive their first nonzero contribution at two loops.

In all panels of figures 12 and 13, we show the result obtained with either
the one-loop kernel or with the sum of one- and two-loop kernels. In both
cases we take the same NLO PDFs, so that the difference between the LO and the
LO+NLO curves directly shows the impact of the two-loop kernels. We find that
the size of the $a_{s}^{2}$ corrections is often moderate but becomes large in
several kinematic situations.

  1. 1.

As discussed in section 4.3 of [20], the two-loop corrections for $g\to
q\bar{q}$ and $q\to q^{\prime}\mskip-1.5mu\bar{q}^{\prime}$ are enhanced at
small $x_{1}+x_{2}$. This is seen in the left panels of figures 12(a), 12(b),
13(b), and 13(c). As follows from equations (4.46) and (4.48) in [20], the
enhanced corrections provide a positive contribution to $F_{u\bar{u}}^{88}$.

An all-order resummation of the enhanced corrections using techniques from
small-$x$ factorisation may be possible, but details of this have not been
worked out.

  2. 2.

As explained in section 4.4 of [20], the splitting graph for $u\to u\bar{u}$
in figure 5(d) leads to a behaviour of the DPD like $1/(1-u)\approx
x_{1}/x_{2}$ for $x_{2}\ll x_{1}$, which is absent at order $a_{s}$. This
explains why in figure 13(a) the LO+NLO result for the scaled DPD
$x_{1}x_{2}\mskip 1.5muF_{u\bar{u}}^{88}$ goes to a finite value when
$x_{2}\ll x_{1}$, whereas the LO result goes to zero. It also explains the
huge relative NLO corrections seen in the right panel of figure 13(b). In the
latter case, the splitting process $u\to u\bar{u}$ is further enhanced by the
fact that the $u$ quark distribution becomes the dominant PDF with increasing
$x$.

The appearance of an additional power $1/(1-u)$ is unique for the step from LO
to NLO in this channel and will not repeat itself at yet higher orders.

  3. 3.

For $x_{2}\gg x_{1}$ we find a large negative two-loop contribution to
$F_{u\bar{u}}^{88}$, which is seen in the right panel of figure 13(c) and more
clearly in figure 13(d), where the LO+NLO result becomes just slightly
negative. This can be traced back to a negative term

$\displaystyle-\frac{5a_{s}^{2}}{24}\,\delta(1-z)\,{}^{11}\mskip-1.5muV_{q\bar{q},g}^{(1)}(u)\,\log^{2}\mskip-1.5muu$
(5.13)

in the kernel $V_{q\bar{q},g}^{88}$, which is enhanced by two powers of $\log
u$ compared with the LO expression.

It would require further analysis to understand whether this type of
enhancement repeats itself at yet higher orders and, if so, whether it can be
resummed to all orders. We therefore cannot say whether the negative values
seen in figure 13(d) will disappear when higher order terms are included.

Notice that figures 13(c) and figure 13(d) refer to the lowest scale
$\mu=\mu_{y}/1.2$ considered in this study. The corresponding plots for
$\mu=\mu_{y}$ or $1.2\mskip 1.5mu\mu_{y}$ show large negative NLO corrections
as well, but the LO+NLO curves do not become negative any more.

We remark in passing that the enhanced two-loop contributions described in
points 1 and 3 scale like the LO kernel for $g\to q\bar{q}$ and hence cancel
out in $F_{u\bar{u}}^{11}$. By contrast, the $1/(1-u)$ behaviour discussed in
point 2 appears both in $F_{u\bar{u}}^{88}$ and $F_{u\bar{u}}^{11}$.

#### Colour summed distributions

Let us finally consider the colour summed distributions
${}^{11}\mskip-1.5muF_{ud}$, ${}^{11}\mskip-1.5muF_{u\smash{\bar{d}}}$,
${}^{11}\mskip-1.5muF_{uu}$, and ${}^{11}\mskip-1.5muF_{u\bar{u}}$. The
relations (2.3) and (5.6) imply that ${}^{11}\mskip-1.5muF_{ud}$ and
${}^{11}\mskip-1.5muF_{u\smash{\bar{d}}}$ are negative whenever their $s$
channel counterparts are. By contrast, we find that
${}^{11}\mskip-1.5muF_{uu}$ and ${}^{11}\mskip-1.5muF_{u\bar{u}}$ remain
positive for the kinematic settings shown in figures 9 to 13 .

##  6 Summary

In the context of the parton model, DPDs in the $s$ channel colour basis are
probability densities for finding two partons in a definite colour state (with
specified longitudinal momenta and specified transverse distance from each
other). This leads to the expectation that $s$ channel DPDs for unpolarised
partons should be non-negative. If it is satisfied, this positivity property
provides valuable constraints on the colour dependence of DPDs, along with a
strategy to model them by saturating the positivity bounds at a certain scale.

In the present work, we investigate whether evolution of DPDs to higher scales
preserves the positivity property under the assumption that it holds at the
starting scale. We limit ourselves to unpolarised partons, and we exclude two-
gluon DPDs from our consideration because their colour structure is much more
involved than the one of pure quark or quark-gluon distributions.

DPDs in the $s$ channel colour basis are subject to Collins-Soper evolution in
the rapidity parameter $\zeta_{p}$, and this evolution includes mixing between
different colour channels. Provided that the renormalisation scales $\mu_{1}$
and $\mu_{2}$ associated with the two partons are large enough compared with
$1/y$, the Collins-Soper kernel ${}^{8}\mskip-1.5muJ(y,\mu_{1},\mu_{2})$ for
DPDs is negative. Under this condition, evolution to higher $\zeta_{p}$
preserves positivity. Conversely, backward evolution to sufficiently small
$\zeta_{p}$ eventually leads to negative distributions (except for the special
case in which all $t$ channel distributions other than
${}^{11}\mskip-1.5muF_{a_{1}a_{2}}$ are zero and hence all $s$ channel
distributions are independent of $\zeta_{p}$).

We next consider the DGLAP equations for evolution in one of the scales
$\mu_{1}$ or $\mu_{2}$, with the evolution kernels taken at LO. We find that
evolution to higher scales is not guaranteed to preserve positivity: there are
initial conditions that satisfy positivity but lead to negative $s$ channel
distributions at higher scales. This is due to the convolution of a plus
distribution in the evolution kernel with a DPD different from the one being
evolved. There is no contribution of this type in the LO evolution equations
for polarised colour summed DPDs ${}^{11}\mskip-1.5muF_{a_{1}a_{2}}$, which
conserve positivity in the same way as the LO evolution of polarised PDFs [4].
In a numerical illustration, we choose initial conditions where certain $s$
channel DPDs are zero and see that they turn negative at slightly higher
scales. We also study joint DGLAP and Collins-Soper evolution in the common
scale $\mu=\mu_{1}=\mu_{2}=\sqrt{x_{1}x_{2}\mskip 1.5mu\zeta_{p}}$ and find
that positivity is not preserved, for the same reasons as above.

At small inter-parton distance $y$, the initial conditions for DPD evolution
can be computed using the perturbative splitting mechanism, setting
$\mu=\mu_{1}=\mu_{2}\approx 1/y$ and using a fixed-order truncation of the DPD
splitting kernels. It is easy to see that at order $a_{s}$ one obtains DPDs
that satisfy positivity in colour space. At order $a_{s}^{2}$ this no longer
holds: in a numerical study we obtain negative values for colour channels in
which the distributions are zero at order $a_{s}$, and also for the
distribution $F_{u\bar{u}}^{88}$, which is nonzero at order $a_{s}$. Negative
values are also found for the colour summed distributions
${}^{11}\mskip-1.5muF_{ud}$ and ${}^{11}\mskip-1.5muF_{u\smash{\bar{d}}}$. In
several cases, the considered distributions have no $\zeta_{p}$ dependence at
order $a_{s}^{2}$, and negative values of them can be uniquely traced back to
the subtraction of ultraviolet divergences implied in the definition of twist-
two operators. The explicit form of this subtraction is given in section 2.6
of [20]. The negative values we find for the distributions are small compared
with the size of the same distributions at other values of the momentum
fractions $x_{1}$ and $x_{2}$. In this sense, the violations of positivity we
have seen may be regarded as “relatively small”. In view of this, we should
also caution that negative values obtained with the splitting formula at order
$a_{s}^{2}$ may turn into positive ones when yet higher orders are included.

Note that the violations of positivity just discussed refer to DPDs defined
with $\overline{\text{MS}}$ renormalisation of twist-two operators. By
contrast, the violation of positivity by forward DGLAP evolution described
earlier occurs at LO and is hence not specific to the $\overline{\text{MS}}$
scheme.

We conclude that the positivity of DPDs in full colour space cannot be taken
for granted and can be violated in physically realistic settings. Using
positivity as a guide for modelling DPDs at large $y$ may still be an option
when there is a lack of better information. It should however be done with due
caution, and one should check whether the chosen initial conditions give
positive distributions when evolved to higher scales.

##  Appendix A Colour space projectors

In this appendix, we list the colour space projectors that appear in the
definitions (2.1) and (2.2) of DPDs in the $s$ and $t$ channel bases. In the
pure quark and the mixed quark-gluon sector, we have

$\displaystyle P_{\rule{0.0pt}{4.82224pt}11}^{ii^{\prime}\mskip
1.5mujj^{\prime}}$
$\displaystyle=\frac{1}{3}\,\delta_{i^{\prime}i}\,\delta_{j^{\prime}\mskip-1.5muj}\,,$
$\displaystyle P_{\rule{0.0pt}{4.82224pt}88}^{ii^{\prime}\mskip
1.5mujj^{\prime}}$ $\displaystyle=2\mskip
1.5mut^{a}_{i^{\prime}i}\,t^{a}_{j^{\prime}\mskip-1.5muj}\,,$ (A.1)
$\displaystyle P_{\rule{0.0pt}{4.82224pt}11}^{aa^{\prime}\mskip
1.5muii^{\prime}}$
$\displaystyle=\frac{1}{\sqrt{24}}\,\delta_{i^{\prime}\mskip-1.5mui}\,\delta^{aa^{\prime}}\,,$
$\displaystyle P_{\rule{0.0pt}{4.82224pt}S8}^{aa^{\prime}\mskip
1.5muii^{\prime}}$
$\displaystyle=\sqrt{\frac{6}{5}}\,d^{aa^{\prime}c}\,t^{c}_{i^{\prime}\mskip-1.5mui}\,,$
$\displaystyle P_{\rule{0.0pt}{4.82224pt}A8}^{aa^{\prime}\mskip
1.5muii^{\prime}}$
$\displaystyle=\sqrt{\frac{2}{3}}\,f^{aa^{\prime}c}\,t^{c}_{i^{\prime}\mskip-1.5mui}\,,$
(A.2)

and

$\displaystyle P_{\rule{0.0pt}{4.82224pt}3\mskip
1.5mu\overline{3}}^{ij\,i^{\prime}j^{\prime}}$
$\displaystyle=\frac{1}{2}\,\bigl{(}\delta_{i^{\prime}i}\delta_{j^{\prime}j}-\delta_{i^{\prime}j}\delta_{j^{\prime}i}\bigr{)}\,,$
$\displaystyle P_{\rule{0.0pt}{4.82224pt}\overline{6}\mskip
1.5mu6}^{ij\,i^{\prime}j^{\prime}}$
$\displaystyle=\frac{1}{2}\,\bigl{(}\delta_{i^{\prime}i}\delta_{j^{\prime}j}+\delta_{i^{\prime}j}\delta_{j^{\prime}i}\bigr{)}\,,$
(A.3) $\displaystyle P_{\rule{0.0pt}{4.82224pt}\overline{3}\mskip
1.5mu3}^{ia\,i^{\prime}\mskip-1.5mua^{\prime}}$
$\displaystyle=\frac{3}{4}\,\bigl{(}t^{a^{\prime}}t^{a}\bigr{)}_{i^{\prime}\mskip-1.5mui}\,,$
$\displaystyle P_{\rule{0.0pt}{4.82224pt}6\mskip
1.5mu\overline{6}}^{ia\,i^{\prime}\mskip-1.5mua^{\prime}}$
$\displaystyle=\frac{1}{2}\,\delta_{i^{\prime}\mskip-1.5mui}\,\delta^{aa^{\prime}}-\bigl{(}t^{a}t^{a^{\prime}}\bigr{)}_{i^{\prime}\mskip-1.5mui}-\frac{1}{2}\,\bigl{(}t^{a^{\prime}}t^{a}\bigr{)}_{i^{\prime}\mskip-1.5mui}\,,$
$\displaystyle
P_{\rule{0.0pt}{4.82224pt}\overline{15}\,15}^{ia\,i^{\prime}\mskip-1.5mua^{\prime}}$
$\displaystyle=\frac{1}{2}\,\delta_{i^{\prime}\mskip-1.5mui}\,\delta^{aa^{\prime}}+\bigl{(}t^{a}t^{a^{\prime}}\bigr{)}_{i^{\prime}\mskip-1.5mui}-\frac{1}{4}\,\bigl{(}t^{a^{\prime}}t^{a}\bigr{)}_{i^{\prime}\mskip-1.5mui}\,.$
(A.4)

For completeness, we also give the projectors for the pure gluon sector,
although they are not used in the present work.

$\displaystyle P_{\rule{0.0pt}{4.82224pt}11}^{aa^{\prime}\,bb^{\prime}}$
$\displaystyle=\frac{1}{8}\,\delta^{aa^{\prime}}\delta^{bb^{\prime}}\,,$
$\displaystyle P_{\rule{0.0pt}{4.82224pt}SS}^{aa^{\prime}\,bb^{\prime}}$
$\displaystyle=\frac{3}{5}\,d^{aa^{\prime}c}d^{bb^{\prime}c}\,,\qquad\qquad
P_{\rule{0.0pt}{4.82224pt}AA}^{aa^{\prime}\,bb^{\prime}}=\frac{1}{3}\,f^{aa^{\prime}c}f^{bb^{\prime}c}\,,\qquad\qquad
P_{\rule{0.0pt}{4.82224pt}S\mskip-1.5muA}^{aa^{\prime}\,bb^{\prime}}=\frac{1}{\sqrt{5}}\,d^{aa^{\prime}c}f^{bb^{\prime}c}\,,$
$\displaystyle P_{10\mskip 1.5mu\overline{10}}^{aa^{\prime}\,bb^{\prime}}$
$\displaystyle=\frac{1}{4}\,\bigl{(}\delta^{ab}\delta^{a^{\prime}b^{\prime}}-\delta^{ab^{\prime}}\delta^{a^{\prime}b}\bigr{)}-\frac{1}{2}\,P_{AA}^{aa^{\prime}\,bb^{\prime}}-\frac{i}{4}\,\bigl{(}d^{abc}f^{a^{\prime}b^{\prime}c}+f^{abc}d^{a^{\prime}b^{\prime}c}\bigr{)}\,,$
$\displaystyle P_{\rule{0.0pt}{4.82224pt}27\mskip
1.5mu27}^{aa^{\prime}\,bb^{\prime}}$
$\displaystyle=\frac{1}{2}\,\bigl{(}\delta^{ab}\delta^{a^{\prime}b^{\prime}}+\delta^{ab^{\prime}}\delta^{a^{\prime}b}\bigr{)}-P_{SS}^{aa^{\prime}\,bb^{\prime}}-P_{11}^{aa^{\prime}\,bb^{\prime}}\,.$
(A.5)

In all cases, we have set the number of colours to $N=3$. Further projectors
are obtained by exchanging the representation labels and the corresponding
indices:

$\displaystyle P_{\rule{0.0pt}{4.82224pt}RR^{\mskip
1.5mu\prime}}^{r_{1}r_{2}\,r_{1}^{\prime}r_{2}^{\prime}}$
$\displaystyle=P_{\rule{0.0pt}{4.82224pt}R^{\mskip
1.5mu\prime}\mskip-1.5muR}^{r_{1}^{\prime}r_{2}^{\prime}\,r_{1}r_{2}}\,.$
(A.6)

One readily verifies the completeness relations

$\displaystyle\sum_{R}P_{\rule{0.0pt}{4.82224pt}R\overline{R}}^{r_{1}r_{2}\,r_{2}^{\prime}r_{1}^{\prime}}$
$\displaystyle=\delta_{r_{1}r_{1}^{\prime}}\,\delta_{r_{2}r_{2}^{\prime}}$
$\displaystyle\text{for $P_{R\overline{R}}$ in \eqref{qq-pro}}\,,$
$\displaystyle\sum_{R}P_{\rule{0.0pt}{4.82224pt}R\overline{R}}^{r_{1}r_{2}\,r_{1}^{\prime}r_{2}^{\prime}}$
$\displaystyle=\delta_{r_{1}r_{1}^{\prime}}\,\delta_{r_{2}r_{2}^{\prime}}$
$\displaystyle\text{for $P_{R\overline{R}}$ in \eqref{qq-pro-s}, \eqref{qg-
pro-s}, \eqref{gg-pro}}\,,$ (A.7)

where in each case the sum runs over all available representations. The
multiplicity of a representation $R$ can be computed from the trace

$\displaystyle m(R)$
$\displaystyle=\begin{cases}P_{R\overline{R}}^{r_{1}r_{2}\,r_{2}r_{1}}&\text{for
$P_{R\overline{R}}$ in \eqref{qq-pro},}\\\\[6.00006pt]
P_{R\overline{R}}^{r_{1}r_{2}\,r_{1}r_{2}}&\text{for $P_{R\overline{R}}$ in
\eqref{qq-pro-s}, \eqref{qg-pro-s}, \eqref{gg-pro}.}\end{cases}$ (A.8)

## Acknowledgements

This work was in part supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) – Research Unit FOR 2926, grant number 409651613.
The work of JRG is supported by the Royal Society through Grant URF\R1\201500.

## References

  * [1] J. Soffer, _Positivity constraints for spin dependent parton distributions_ , _Phys. Rev. Lett._ 74 (1995) 1292 [hep-ph/9409254].
  * [2] A. Bacchetta, M. Boglione, A. Henneman and P. J. Mulders, _Bounds on transverse momentum dependent distribution and fragmentation functions_ , _Phys. Rev. Lett._ 85 (2000) 712 [hep-ph/9912490].
  * [3] M. Diehl and P. Hägler, _Spin densities in the transverse plane and generalized transversity distributions_ , _Eur. Phys. J. C_ 44 (2005) 87 [hep-ph/0504175].
  * [4] M. Diehl and T. Kasemets, _Positivity bounds on double parton distributions_ , _JHEP_ 05 (2013) 150 [1303.0842].
  * [5] P. Bartalini and J. R. Gaunt, eds., _Multiple Parton Interactions at the LHC_ , vol. 29. WSP, 2019, 10.1142/10646.
  * [6] T. Kasemets and P. J. Mulders, _Constraining double parton correlations and interferences_ , _Phys. Rev. D_ 91 (2015) 014015 [1411.0726].
  * [7] G. Altarelli, S. Forte and G. Ridolfi, _On positivity of parton distributions_ , _Nucl. Phys. B_ 534 (1998) 277 [hep-ph/9806345].
  * [8] A. Candido, S. Forte and F. Hekhorn, _Can $\overline{\mathrm{MS}}$ parton distributions be negative?_, _JHEP_ 11 (2020) 129 [2006.07377].
  * [9] S. J. Brodsky and G. P. Lepage, _Exclusive Processes in Quantum Chromodynamics_ , _Adv. Ser. Direct. High Energy Phys._ 5 (1989) 93.
  * [10] M. Diehl, _Generalized parton distributions in impact parameter space_ , _Eur. Phys. J. C_ 25 (2002) 223 [hep-ph/0205208].
  * [11] J. R. Gaunt, _Double parton scattering in proton-proton collisions_ , Ph.D. thesis, University of Cambridge, 2012.  https://doi.org/10.17863/CAM.16589.
  * [12] T. Rogers, _Renormalization schemes for parton densities_ , talk at the QCD Evolution Workshop 2021, https://indico.bnl.gov/event/6803/contributions/49158.
  * [13] L. Durand and W. Putikka, _Probabilistic Derivation of Parton Splitting Functions_ , _Phys. Rev. D_ 36 (1987) 2840.
  * [14] J. C. Collins and J.-w. Qiu, _A New Derivation of the Altarelli-Parisi Equations_ , _Phys. Rev. D_ 39 (1989) 1398.
  * [15] V. Barone, _On the QCD evolution of the transversity distribution_ , _Phys. Lett. B_ 409 (1997) 499 [hep-ph/9703343].
  * [16] C. Bourrely, J. Soffer and O. V. Teryaev, _The $Q^{2}$ evolution of Soffer inequality_, _Phys. Lett. B_ 420 (1998) 375 [hep-ph/9710224].
  * [17] W. Vogelsang, _Next-to-leading order evolution of transversity distributions and Soffer’s inequality_ , _Phys. Rev. D_ 57 (1998) 1886 [hep-ph/9706511].
  * [18] O. Martin, A. Schäfer, M. Stratmann and W. Vogelsang, _Soffer’s inequality and the transversely polarized Drell-Yan process at next-to-leading order_ , _Phys. Rev. D_ 57 (1998) 3084 [hep-ph/9710300].
  * [19] M. Diehl and P. Stienemeier, _Gluons and sea quarks in the proton at low scales_ , _Eur. Phys. J. Plus_ 135 (2020) 211 [1904.10722].
  * [20] M. Diehl, J. R. Gaunt and P. Plößl, _Two-loop splitting in double parton distributions: the colour non-singlet case_ , _JHEP_ 08 (2021) 040 [2105.08425].
  * [21] M. G. A. Buffing, M. Diehl and T. Kasemets, _Transverse momentum in double parton scattering: factorisation, evolution and matching_ , _JHEP_ 01 (2018) 044 [1708.03528].
  * [22] J. Collins, _Foundations of perturbative QCD_ , vol. 32. Cambridge University Press, 2013.
  * [23] M. Mekhfi, _Correlations in Color and Spin in Multiparton Processes_ , _Phys. Rev. D_ 32 (1985) 2380.
  * [24] M. Diehl, D. Ostermeier and A. Schäfer, _Elements of a theory for multiparton interactions in QCD_ , _JHEP_ 03 (2012) 089 [1111.0910].
  * [25] R. Slansky, _Group Theory for Unified Model Building_ , _Phys. Rept._ 79 (1981) 1.
  * [26] A. Vladimirov, _Soft factors for double parton scattering at NNLO_ , _JHEP_ 12 (2016) 038 [1608.04920].
  * [27] J. Collins and T. Rogers, _Understanding the large-distance behavior of transverse-momentum-dependent parton densities and the Collins-Soper evolution kernel_ , _Phys. Rev. D_ 91 (2015) 074020 [1412.3820].
  * [28] A. Bacchetta, V. Bertone, C. Bissolotti, G. Bozzi, F. Delcarro, F. Piacenza et al., _Transverse-momentum-dependent parton distributions up to N 3LL from Drell-Yan data_, _JHEP_ 07 (2020) 117 [1912.07550].
  * [29] I. Scimemi and A. Vladimirov, _Non-perturbative structure of semi-inclusive deep-inelastic and Drell-Yan scattering at small transverse momentum_ , _JHEP_ 06 (2020) 137 [1912.06532].
  * [30] P. Shanahan, M. Wagman and Y. Zhao, _Collins-Soper kernel for TMD evolution from lattice QCD_ , _Phys. Rev. D_ 102 (2020) 014511 [2003.06063].
  * [31] Lattice Parton collaboration, Q.-A. Zhang et al., _Lattice QCD Calculations of Transverse-Momentum-Dependent Soft Function through Large-Momentum Effective Theory_ , _Phys. Rev. Lett._ 125 (2020) 192001 [2005.14572].
  * [32] M. Schlemmer, A. Vladimirov, C. Zimmermann, M. Engelhardt and A. Schäfer, _Determination of the Collins-Soper Kernel from Lattice QCD_ , _JHEP_ 08 (2021) 004 [2103.16991].
  * [33] C. F. Berger, C. Marcantonini, I. W. Stewart, F. J. Tackmann and W. J. Waalewijn, _Higgs Production with a Central Jet Veto at NNLL+NNLO_ , _JHEP_ 04 (2011) 092 [1012.4480].
  * [34] S. Dulat, T.-J. Hou, J. Gao, M. Guzzi, J. Huston, P. Nadolsky et al., _New parton distribution functions from a global analysis of quantum chromodynamics_ , _Phys. Rev. D_ 93 (2016) 033006 [1506.07443].
  * [35] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _Eur. Phys. J. C_ 75 (2015) 132 [1412.7420].
  * [36] D. B. Clark, E. Godat and F. I. Olness, _ManeParse : A Mathematica reader for Parton Distribution Functions_ , _Comput. Phys. Commun._ 216 (2017) 126 [1605.08012].