DESY 19-024, CERN-TH-2019-013

Two-loop splitting in double parton distributions

Markus Diehl1, Jonathan R. Gaunt2, Peter Plößl3, Andreas Schäfer3

1 Deutsches Elektronen-Synchroton DESY, 22603 Hamburg, Germany

2 CERN Theory Division, 1211 Geneva 23, Switzerland

3 Institut für Theoretische Physik, Universität Regensburg, 93040 Regensburg,
Germany

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>regensburg.de

July 18, 2019

## Abstract

Double parton distributions (DPDs) receive a short-distance contribution from
a single parton splitting to yield the two observed partons. We investigate
this mechanism at next-to-leading order (NLO) in perturbation theory.
Technically, we compute the two-loop matching of both the position and
momentum space DPDs onto ordinary PDFs. This also yields the $1\to 2$
splitting functions appearing in the evolution of momentum-space DPDs at NLO.
We give results for the unpolarised, colour-singlet DPDs in all partonic
channels. These quantities are required for calculations of double parton
scattering at full NLO. We discuss various kinematic limits of our results,
and we verify that the $1\to 2$ splitting functions are consistent with the
number and momentum sum rules for DPDs.

###### Contents

  1. 1 Introduction
  2. 2 General analysis: ultraviolet behaviour and scale dependence
    1. 2.1 Basic definitions and notation
    2. 2.2 Renormalisation and evolution of DPDs
    3. 2.3 Splitting kernels at higher order
      1. 2.3.1 Momentum space kernels
      2. 2.3.2 Position space kernels
  3. 3 Matching between momentum and position space DPDs
    1. 3.1 Matching at zero $\Delta$
    2. 3.2 Matching at nonzero $\Delta$
    3. 3.3 Scale independence of matching
  4. 4 Two-loop calculation
    1. 4.1 Channels and graphs
    2. 4.2 Performing the calculation
      1. 4.2.1 Real emission diagrams
      2. 4.2.2 Virtual loop diagrams
  5. 5 Results
    1. 5.1 Parton number and momentum sum rules
    2. 5.2 Active-spectator symmetry
    3. 5.3 Kinematic limits
      1. 5.3.1 Threshold limit: large $x_{1}+x_{2}$
      2. 5.3.2 Small $x_{1}+x_{2}$
      3. 5.3.3 Small $x_{1}$ or $x_{2}$
    4. 5.4 Numerical illustration
    5. 5.5 Comparison to other results in the literature
    6. 5.6 The $\mathcal{N}=1$ SUSY Relation
  6. 6 Conclusion

##  1 Introduction

Consider the integrated cross section for the production of some final state
with associated hard scale $Q$ at a hadron collider. As is well known, the
dominant mechanism for this is single parton scattering (SPS), where one
parton from each proton collides to form the state of interest. If the final
state can be divided into two subsets $A$ and $B$ with associated hard scales
$Q_{A}$ and $Q_{B}$, then another mechanism is two separate parton-parton
collisions, one yielding $A$ and the other yielding $B$. This is called double
parton scattering, or DPS. The integrated cross section for this mechanism is
suppressed compared to that for SPS by $\Lambda^{2}/Q^{2}$, where $\Lambda$
denotes a hadronic scale and $Q=\min(Q_{A},Q_{B})$. Despite this, there are
several motivations to study this process, and DPS occupies a special place
amongst the power suppressed corrections for the production of $A+B$. First,
DPS may compete with SPS when the latter is suppressed by small or multiple
coupling constants, with a well-known example being same-sign $W$ pair
production [1, 2, 3, 4, 5, 6]. When followed by leptonic $W$ decays, this is
an important process, since same-sign dilepton production is a search channel
for new physics (see e.g. [7, 8, 9]). Second, although DPS is a power
correction at the level of the integrated cross section, in differential cross
sections there can be regions in which DPS is competitive with SPS. In
particular, for small transverse momenta of $A$ and $B$, DPS and SPS
contribute at the same power [10, 11, 12, 13]. For several processes one also
finds that DPS is an important contribution at large rapidity separation of
$A$ and $B$ (see e.g. [14, 15, 7]). Third, the DPS cross section scales
roughly as the fourth power of a parton distribution, whilst SPS scales as the
second power (and other $\Lambda^{2}/Q^{2}$ corrections scale as the third
power). Thus for a given final state, as we increase the collider energy and
probe smaller parton momentum fractions $x$, DPS grows faster than the other
contributions, and hence is especially relevant at the LHC and future, higher
energy colliders. Finally, DPS is an interesting process to study in its own
right: it provides information on the correlation between partons in the
proton and thus reveals an aspect of hadron structure that remains
inaccessible in single parton distributions.

A variety of DPS processes have been investigated at the ISR [16], the CERN
SPS [17], the Tevatron [18, 19, 20, 21, 22, 23, 24, 25, 26], and the LHC [27,
28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 4] (see e.g. figure 4 of [35] and
figure 15 of [38] for overviews). These measurements cover final states with
jets, electroweak gauge bosons, heavy quarkonia, and open charm. The DPS
contribution to pair production of heavy gauge or Higgs bosons (such as
$W^{+}W^{-}$ or $H\mskip 1.5muZ$) is estimated to be small, but may well
become important in the high-luminosity phase of the LHC, given the
overarching goal to probe electroweak symmetry breaking with highest possible
precision. A detailed review of experimental, phenomenological, and
theoretical aspects of DPS can be found in [39].

Either one or both parton pairs in a DPS process can arise from a single
parton via a perturbative “$1\to 2$” branching process. At next-to-leading
order (NLO) and above, this branching may produce additional partons that are
emitted into the final state, in addition to the two “active” partons that
subsequently enter hard collisions. The process in which we have $1\to 2$
splittings in both protons overlaps with a loop correction to SPS, in the
sense that a given Feynman graph corresponds to either DPS or SPS depending on
the momenta of its internal lines. Likewise, the process where we have a $1\to
2$ splitting in just one proton overlaps with a “twist 2 $\times$ twist 4”
power-suppressed process. Developing a theoretical QCD framework for DPS that
includes $1\to 2$ splitting but avoids double counting between SPS and DPS
(and the twist 2 $\times$ twist 4 mechanism) is quite nontrivial. First
approaches to this problem are given in [13, 40, 41, 42]. More recently, some
of us developed a new formalism [43] that possesses several advantages over
the previous approaches. This involves a description of the DPS cross section
in terms of so-called position space double parton distributions (DPDs)
$F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})$, where $\boldsymbol{y}$ is the
separation in transverse space between the two partons.

One of the advantages of the approach in [43] is that it can be formulated at
all orders in $\alpha_{s}$, with higher-order contributions that can be
practicably computed. This opens the perspective to performing DPS
calculations at full next-to-leading order accuracy. The framework of [43] is
designed to make maximal re-use of quantities that are already computed for
SPS. The partonic scattering cross sections $\hat{\sigma}$ for DPS and for SPS
are identical. NLO effects in parton emissions from the incoming legs are
correctly described by including the known NLO corrections to the DGLAP
splitting kernels [44, 45, 46, 47, 48, 49, 50, 51, 52] in the corresponding
kernels of the “double DGLAP” equation, which describes the scale evolution of
$F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})$. A missing key ingredient,
specific to DPS, is the effect of the $1\to 2$ splitting at NLO. This is the
subject of the present paper. What is technically required in the notation of
[43] is the two-loop kernel $V_{a_{1}a_{2},a_{0}}(x_{1},x_{2},\boldsymbol{y})$
for the perturbative matching of the DPDs onto the PDFs $f_{a_{0}}(x)$ at
small $y$.

Why should one consider performing such calculations? Apart from unknown
higher-order corrections, a major source of uncertainty in DPS cross section
calculations is the lack of knowledge of the nonperturbative “input” DPDs.
Unlike ordinary parton distributions (PDFs), these are very poorly constrained
due to the limited experimental information on DPS. One can expect knowledge
of the DPDs to improve in future, as further information emerges from
experiment, and as this information is combined with theoretical constraints
[53, 54] and information emerging from lattice computations (see [55] for work
in this direction). However, it is unlikely that the situation for DPDs will
approach that for PDFs in the near-term future. Despite a considerable
uncertainty on the shape of the nonperturbative DPDs in general, there are a
number of motivations for analysing DPS at NLO accuracy. Concerning the
partonic cross sections $\hat{\sigma}$, it is well known from SPS that in many
cases their LO expressions give only a rough approximation, especially for
final states containing jets. On the side of DPDs, the $1\to 2$ splitting
mechanism requires only PDFs as a nonperturbative input and can thus be
brought under good theoretical control. Several studies have emphasised the
quantitative impact of this mechanism on DPS [56, 57, 58, 59, 60, 61, 62, 43].
It is thus natural to go beyond the LO approximation for this contribution and
compute the splitting at order $\alpha_{s}^{2}$. Moreover, for some parton
combinations, this is the first order that gives a nonzero result. An example
is the $u\bar{d}$ distribution, which is in particular relevant for same-sign
$W$ pair production. The recent study in [63] finds that the splitting
mechanism has some impact on observables even in such a situation. A general
aim of going beyond LO is to get a handle on the perturbative convergence. As
previously mentioned, favourable regions for DPS typically involve small $x$
values, and in this region there is potential that large logarithms
$\log(1/x)$ can have a significant impact, and even spoil the perturbative
convergence. By studying the effect of the NLO corrections on DPDs, and more
generally on DPS luminosities and cross sections, we can get an indication of
whether there are convergence issues caused (for example) by small $x$
logarithms. Finally, it is quite possible that by constructing appropriate
ratios of cross sections (involving, for example, the same parton flavours but
different scales), one may minimise the dependence on the nonperturbative
input and more directly test the perturbative part of the cross section. In
this case, NLO computations clearly should improve the overall accuracy.

In this paper we compute the matching coefficients
$V_{a_{1}a_{2},a_{0}}(x_{1},x_{2},\boldsymbol{y})$ at two loops for all
partonic channels (i.e. for all combinations of parton labels
$a_{0},a_{1},a_{2}$). We restrict ourselves to the matching coefficients for
the colour-singlet, unpolarised DPDs. These are expected to be the dominant
DPDs at small $x$ and high $Q$ [64, 65]; in fact they are the only DPDs
considered in most phenomenological analyses of DPS. The coefficients for the
colour interference and polarised DPDs will be presented in future work.

Aside from DPDs $F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})$ in position
space, one can also consider momentum space DPDs
$F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{\Delta})$. The variable
$\boldsymbol{\Delta}$ is Fourier conjugate to $\boldsymbol{y}$ and represents
the difference in transverse momentum between a parton entering the DPD
operator in the amplitude and the corresponding parton in the conjugate
amplitude. A description of DPS in terms of momentum space DPDs is considered
in a number of works [12, 13, 56, 66, 61]. Furthermore, the momentum-space
DPDs at $\boldsymbol{\Delta=0}$ are the objects that satisfy the number and
momentum sum rules of [53], as was discussed in [54]. Similarly to the
position space DPDs at small $\boldsymbol{y}$, momentum space DPDs at large
$\boldsymbol{\Delta}$ can be matched onto PDFs with perturbative kernels
$W_{a_{1}a_{2},a_{0}}(x_{1},x_{2},\boldsymbol{\Delta})$. For the momentum
space DPDs, the effect of $1\to 2$ splittings also shows up in the evolution
equations, which contain an inhomogeneous term where the PDFs are convolved
with “$1\to 2$ evolution kernels” $P_{a_{1}a_{2},a_{0}}(x_{1},x_{2})$ [67, 68,
66, 11]. In [54] it was shown that the kernels $P_{a_{1}a_{2},a_{0}}$ are
related by sum rules to the usual DGLAP splitting functions at all orders in
$\alpha_{s}$. Finally, one can introduce matching kernels
$U_{a_{1}a_{2},a_{0}}(x_{1},x_{2})$ that allow one to compute the
distributions $F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{\Delta=0})$ that
correspond to a given set of position space DPDs. This allows one to make use
of the sum rule constraints of [53] also if one calculates DPS cross sections
in terms of position space DPDs. The computation of the kernels
$U_{a_{1}a_{2},a_{0}}$ at the leading order level is given in section 7 of
[43] and will be extended to NLO in the present work.

This paper is organised as follows. In section 2 we define the kernels
$V_{a_{1}a_{2},a_{0}}$, $W_{a_{1}a_{2},a_{0}}$ and $P_{a_{1}a_{2},a_{0}}$,
analyse their renormalisation scale behaviour, and derive relations between
them. The general structure of the matching between momentum and position
space DPDs is discussed in section 3. In section 4 we give details of the two-
loop computation of the matching kernels. We present our results and discuss
various properties of them in section 5. A brief summary of this paper is
given in section 6.

##  2 General analysis: ultraviolet behaviour and scale dependence

In the next two subsections we briefly review the definitions of bare and of
renormalised DPDs, setting up our notation on the way. In section 2.3 we then
derive how renormalisation affects the splitting kernels for DPDs and their
logarithmic dependence on the renormalisation scale $\mu$. This will later
allow us to obtain renormalised kernels from the computation of bare Feynman
graphs for the relevant splitting processes.

###  2.1 Basic definitions and notation

To begin with let us recall that bare, i.e. unrenormalised, DPDs for partons
$a_{1}$ and $a_{2}$ are defined by matrix elements

$\displaystyle F_{B,a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})$
$\displaystyle=(x_{1}\mskip 1.5mup^{+})^{-n_{1}}\,(x_{2}\mskip
1.5mup^{+})^{-n_{2}}\;2p^{+}\int
dy^{-}\,\frac{dz^{-}_{1}}{2\pi}\,\frac{dz^{-}_{2}}{2\pi}\;e^{i\mskip
1.5mu(x_{1}z_{1}^{-}+x_{2}z_{2}^{-})\mskip 1.5mup^{+}}$
$\displaystyle\quad\times\langle\mskip 1.5mup\mskip
1.5mu|\,\mathcal{O}_{a_{1}}(y,z_{1})\,\mathcal{O}_{a_{2}}(0,z_{2})\,|\mskip
1.5mup\mskip
1.5mu\rangle\bigl{|}_{z_{1}^{+}=z_{2}^{+}=y_{\phantom{1}}^{+}=0\,,\boldsymbol{z}_{1}=\boldsymbol{z}_{2}=\boldsymbol{0}}\,,$
(1)

where $n_{i}=0$ for quarks and antiquarks and $n_{i}=1$ for gluons. We use
light-cone coordinates $v^{\pm}=(v^{0}\pm v^{3})/\sqrt{2}$ for any four-vector
$v^{\mu}$ and write its transverse part in boldface,
$\boldsymbol{v}=(v^{1},v^{2})$. It is understood that the transverse proton
momentum $\boldsymbol{p}$ is zero and that the proton polarisation is averaged
over. The twist-two operators

$\displaystyle\mathcal{O}_{q}(y,z)$
$\displaystyle=\frac{1}{2}\,\bar{q}\bigl{(}y-z/2\bigr{)}\mskip
1.5mu\gamma^{+}\mskip 1.5muq\bigl{(}y+z/2\bigr{)}\,,$
$\displaystyle\mathcal{O}_{\bar{q}}(y,z)$
$\displaystyle=-\frac{1}{2}\,\bar{q}\bigl{(}y+z/2\bigr{)}\mskip
1.5mu\gamma^{+}\mskip 1.5muq\bigl{(}y-z/2\bigr{)}\,,$
$\displaystyle\mathcal{O}_{g}(y,z)$
$\displaystyle=G^{+i}\bigl{(}y-z/2\bigr{)}\mskip
1.5muG^{+i}\bigl{(}y+z/2\bigr{)}$ (2)

for unpolarised partons are constructed from bare field operators and are
familiar from the definition of bare PDFs,

$\displaystyle f_{B,a}(x)$ $\displaystyle=(x\mskip
1.5mup^{+})^{-n}\int\frac{dz^{-}}{2\pi}\,e^{i\mskip
1.5muxz^{-}p^{+}}\langle\mskip 1.5mup\mskip
1.5mu|\,\mathcal{O}_{a}(0,z)\,|\mskip 1.5mup\mskip
1.5mu\rangle\bigl{|}_{z^{+}=0\,,\boldsymbol{z}=\boldsymbol{0}}\,.$ (3)

In the present work, we consider colour-singlet DPDs, so that the colour
indices of the quark or gluon fields in (2.1) are implicit and to be summed
over. The form (2.1) holds in light-cone gauge $A^{+}=0$, and in other gauges
Wilson lines need to be inserted between the two parton fields. The DPDs in
(2.1) depend on the transverse distance $\boldsymbol{y}$ between the two
partons. Distributions in transverse momentum representation are obtained by a
Fourier transform

$\displaystyle F_{B,a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{\Delta})$
$\displaystyle=\int
d^{2-2\epsilon}\boldsymbol{y}\;e^{i\boldsymbol{y}\boldsymbol{\Delta}}\,F_{B,a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y})\,,$
(4)

where it is understood that bare quantities are evaluated in $d=4-2\epsilon$
space-time dimensions. After ultraviolet renormalisation, which will be
discussed in section 2.2, one obtains finite parton distributions and can set
$d=4$. Due to their different ultraviolet behaviour, renormalised DPDs in
position and in momentum space are _not_ Fourier conjugate to each other,
unlike their bare counterparts.

In the limit of small $y=|\boldsymbol{y}|$ or large
$\Delta=|\boldsymbol{\Delta}|$, renormalised DPDs can be computed in terms of
renormalised PDFs and perturbative splitting kernels as

$\displaystyle F_{a_{1}a_{2}}(x_{1},x_{2},y,\mu)$
$\displaystyle=\frac{1}{\pi\mskip
1.5muy^{2}}\,\sum_{a_{0}}\int\limits_{x_{1}+x_{2}}^{1}\frac{dz}{z^{2}}\;V_{a_{1}a_{2},a_{0}}\biggl{(}\frac{x_{1}}{z},\frac{x_{2}}{z},a_{s}(\mu),\log\frac{\mu^{2}y^{2}}{b_{0}^{2}}\biggr{)}\,f_{a_{0}}(z,\mu)\,,$
$\displaystyle F_{a_{1}a_{2}}(x_{1},x_{2},\Delta,\mu)$
$\displaystyle=\sum_{a_{0}}\int\limits_{x_{1}+x_{2}}^{1}\frac{dz}{z^{2}}\;W_{a_{1}a_{2},a_{0}}\biggl{(}\frac{x_{1}}{z},\frac{x_{2}}{z},a_{s}(\mu),\log\frac{\mu^{2}}{\Delta^{2}}\biggr{)}\,f_{a_{0}}(z,\mu)\,,$
(5)

where $\mu$ is the renormalisation scale and we abbreviate

$\displaystyle a_{s}$ $\displaystyle=\frac{\alpha_{s}}{2\pi}\,,$
$\displaystyle b_{0}$ $\displaystyle=2e^{-\gamma}\,,$ (6)

with $\gamma$ being the Euler-Mascheroni constant. The general form of the
convolution integrals in (2.1) was already given in [43] and follows from
boost invariance and dimensional analysis. The splitting mechanism expressed
in (2.1) also induces an inhomogeneous term in the evolution equation for
momentum space DPDs, which reads [54]

$\displaystyle\frac{d}{d\ln\mu^{2}}\,F_{a_{1}a_{2}}(x_{1},x_{2},\Delta,\mu)$
$\displaystyle=\sum_{a_{0}}\int\limits_{x_{1}+x_{2}}^{1}\frac{dz}{z^{2}}\;P_{a_{1}a_{2},a_{0}}\biggl{(}\frac{x_{1}}{z},\frac{x_{2}}{z},a_{s}(\mu)\biggr{)}\,f_{a_{0}}(z,\mu)$
$\displaystyle\quad+\\{\text{homogeneous terms}\\}\,,$ (7)

where the homogeneous terms have the same form as in the evolution of
$F_{a_{1}a_{2}}(x_{1},x_{2},y,\mu)$ and are given in (22) below. At leading
order (LO) in $a_{s}$, the kernel $P_{a_{1}a_{2},a_{0}}$ can be expressed in
terms of the familiar LO DGLAP splitting functions, as specified in (37). The
evolution equation (2.1) then takes the form that has been derived long ago
[67, 68]. Furthermore, one has $V_{a_{1}a_{2},a_{0}}=P_{a_{1}a_{2},a_{0}}$ at
LO, as shown in [11]. In the present work we compute all kernels $V$, $W$ and
$P$ at next-to-leading order (NLO) in $a_{s}$.

We refer to $V_{a_{1}a_{2},a_{0}}$ and $W_{a_{1}a_{2},a_{0}}$ as “$1\to 2$
splitting kernels” and to $P_{a_{1}a_{2},a_{0}}$ as “$1\to 2$ evolution
kernel”, where “$1\to 2$” refers to having one initial parton $a_{0}$ that
splits into the two partons $a_{1},a_{2}$ specified in the DPD. Starting at
NLO, the splitting process can involve further unobserved partons in the final
state. At LO the flavour of $a_{0}$ is uniquely fixed for given $a_{1},a_{2}$,
so that the sum over $a_{0}$ can be omitted in (2.1) and (2.1), but this no
longer holds at NLO.

Let us now introduce some elements of notation that will be useful for the
calculations in the following sections.

##### Dimensional regularisation.

Working in $d=4-2\epsilon$ dimensions, we expand $\epsilon$ dependent
quantities $Q(\epsilon)$ in a Laurent series around $\epsilon=0$,

$\displaystyle Q(\epsilon)$ $\displaystyle=\sum_{k}\epsilon^{k}\mskip
1.5mu[Q]_{k}\,.$ (8)

In particular, $[Q]_{0}$ gives the finite piece and $[Q]_{-1}$ the residue of
the $1/\epsilon$ pole.

##### Convolution integrals.

Following [54], we use a compact notation for the different types of
convolution integrals that avoids giving explicit momentum fraction arguments.
In the following, let $D$ be a function of $x_{1},x_{2}$ whilst $A$, $B$, $C$
are functions of one momentum fraction only. We write

$\displaystyle A\underset{1}{\otimes}D$
$\displaystyle=\int\limits_{x_{1}}^{1}\frac{dz}{z}\;A(z)\,D\biggl{(}\frac{x_{1}}{z},x_{2}\biggr{)}\,,$
(9)

where the integration region is determined by the support properties of the
functions. If a one-variable function $A$ is a PDF then $A(x)=0$ for $x<0$ and
$x>1$, and if a two-variable function $D$ is a DPD, then $D(x_{1},x_{2})=0$
for $x_{1}<0$, $x_{2}<0$ and $x_{1}+x_{2}>1$. Corresponding support properties
hold for renormalisation factors and splitting kernels (which may include
delta functions and plus distributions at the endpoints of their support). In
analogy to (9), we define $A\underset{2}{\otimes}D$ and the combined
convolution
$A\underset{1}{\otimes}B\underset{2}{\otimes}D=A\underset{1}{\otimes}[B\underset{2}{\otimes}D]$.

In (2.1) and in (2.1) we encounter a second type of convolution,

$\displaystyle D\underset{12}{\otimes}A$
$\displaystyle=\int\limits_{x_{1}+x_{2}}^{1}\\!\\!\frac{dz}{z^{2}}\;D\biggl{(}\frac{x_{1}}{z},\frac{x_{2}}{z}\biggr{)}A(z)\,.$
(10)

Interchanging the order of integrations, it is straightforward to show that

$\displaystyle[A\underset{1}{\otimes}D]\underset{12}{\otimes}B$
$\displaystyle=A\underset{1}{\otimes}[D\underset{12}{\otimes}B]\,,$
$\displaystyle[A\underset{1}{\otimes}B\underset{2}{\otimes}D]\underset{12}{\otimes}C$
$\displaystyle=A\underset{1}{\otimes}B\underset{2}{\otimes}[D\underset{12}{\otimes}C]\,,$
(11)

so that we can write these convolutions without the brackets. We also have

$\displaystyle[D\underset{12}{\otimes}A]\underset{12}{\otimes}B$
$\displaystyle=D\underset{12}{\otimes}[A\otimes B]\,,\phantom{\int_{0}^{1}}$
(12)

where $A\otimes B$ is the usual convolution for functions of a single momentum
fraction. We will also write $D\underset{12}{\otimes}A\otimes B$ without
brackets. One can rewrite (10) in the form

$\displaystyle D\underset{12}{\otimes}A$
$\displaystyle=\frac{1}{x_{1}+x_{2}}\int\limits_{x_{1}+x_{2}}^{1}\\!\\!\\!dz\;D\biggl{(}z\mskip
1.5mu\frac{x_{1}}{x_{1}+x_{2}},z\mskip
1.5mu\frac{x_{2}}{x_{1}+x_{2}}\biggr{)}\,A\biggl{(}\frac{x_{1}+x_{2}}{z}\biggr{)}\,,$
(13)

which will be useful in section 5.

##### Inverse convolution (one variable).

We define $A^{-1}(x)$ as the solution of $A\otimes A^{-1}=\delta(1-x)$. We
will only use this inversion for functions $A(x)$ that have a series expansion
in $a_{s}$; one then obtains a unique inverse $A^{-1}(x)$ order by order in
perturbation theory. In particular, if $A(x)=1+a_{s}\mskip
1.5muA^{(1)}(x)+\mathcal{O}(a_{s}^{2})$, then $A^{-1}(x)=1-a_{s}\mskip
1.5muA^{(1)}(x)+\mathcal{O}(a_{s}^{2})$.

##### Parton indices.

The notation for convolutions can be extended to include appropriate summation
over indices denoting the parton type (quarks, antiquarks, gluons). For
quantities $A_{ab}(x),B_{ab}(x)$ and $f_{a}(x)$ depending on a single momentum
fraction, we write as usual

$\displaystyle[A\otimes B]_{ac}$ $\displaystyle=\sum_{b}A_{ab}\otimes
B_{bc}\,,$ $\displaystyle[A\otimes f]_{a}$
$\displaystyle=\sum_{b}A_{ab}\otimes f_{b}\,.$ (14)

Note that we do _not_ use the summation convention for parton indices, i.e.
sums are indicated explicitly when indices are given. For quantities
$D_{a_{1}a_{2},a_{0}}(x_{1},x_{2})$ and $F_{a_{1}a_{2}}(x_{1},x_{2})$
depending on two momentum fractions, we define

$\displaystyle[A\underset{1}{\otimes}D]_{a_{1}a_{2},a_{0}}$
$\displaystyle=\sum_{b}A_{a_{1}b}\underset{1}{\otimes}D_{b\mskip
1.5mua_{2},a_{0}}\,,$ $\displaystyle[A\underset{1}{\otimes}F]_{a_{1}a_{2}}$
$\displaystyle=\sum_{b}A_{a_{1}b}\underset{1}{\otimes}F_{b\mskip
1.5mua_{2}}\,,$ $\displaystyle[A\underset{2}{\otimes}D]_{a_{1}a_{2},a_{0}}$
$\displaystyle=\sum_{b}A_{a_{2}b}\underset{1}{\otimes}D_{a_{1}b,a_{0}}\,,$
$\displaystyle[A\underset{2}{\otimes}F]_{a_{1}a_{2}}$
$\displaystyle=\sum_{b}A_{a_{2}b}\underset{1}{\otimes}F_{a_{1}b}\,,$
$\displaystyle[D\underset{12}{\otimes}A]_{a_{1}a_{2},a_{0}}$
$\displaystyle=\sum_{b}D_{a_{1}a_{2},b}\underset{12}{\otimes}A_{b\mskip
1.5mua_{0}}\,,$ $\displaystyle[D\underset{12}{\otimes}f]_{a_{1}a_{2}}$
$\displaystyle=\sum_{b}D_{a_{1}a_{2},b}\underset{12}{\otimes}f_{b}\,.$ (15)

For the combination of convolutions in the first and second momentum fraction,
we then have

$\displaystyle[A\underset{1}{\otimes}B\underset{2}{\otimes}D]_{a_{1}a_{2},a_{0}}$
$\displaystyle=\sum_{b_{1},b_{2}}A_{a_{1}b_{1}}\underset{1}{\otimes}B_{a_{2}b_{2}}\underset{2}{\otimes}D_{b_{1}b_{2},a_{0}}\,.$
(16)

###  2.2 Renormalisation and evolution of DPDs

Let us recall the renormalisation of DPDs, which is well-known at LO and has
been formulated at arbitrary perturbative order in [54]. The twist-two
operators in (2.1) contain short-distance singularities that are renormalised
with a factor $Z(x,\mu)$ as111Our convention for renormalisation factors $Z$
of composite operators corresponds to the one in [69, 70]. Other authors, such
as the ones of [71], use $Z^{-1}$ instead.

$\displaystyle f(\mu)$ $\displaystyle=Z(\mu)\otimes f_{B}$ (17)

for ordinary PDFs. Correspondingly, $F_{B}(y)$ is renormalised by a factor $Z$
for each parton,

$\displaystyle F(y,\mu)$
$\displaystyle=Z(\mu)\underset{1}{\otimes}Z(\mu)\underset{2}{\otimes}F_{B}$
(18)

In addition to the divergences of the twist-two operators, the momentum-space
distributions $F_{B}(\Delta)$ have a $1\to 2$ splitting singularity, which
arises from integrating $F_{B}(y)$ over $\boldsymbol{y}$ in (4) because
$F_{B}(y)\propto y^{-2+2\epsilon}$ at small $y$. This singularity is
renormalised additively with a factor $Z_{s}(x_{1},x_{2},\mu)$ as

$\displaystyle F(\Delta,\mu)$
$\displaystyle=Z(\mu)\underset{1}{\otimes}Z(\mu)\underset{2}{\otimes}F_{B}(\Delta)+Z_{s}(\mu)\underset{12}{\otimes}f_{B}\,.$
(19)

Defining the ordinary DGLAP evolution kernel $P$ by

$\displaystyle\frac{d}{d\log\mu^{2}}\,Z(\mu)$ $\displaystyle=P(\mu)\otimes
Z(\mu)$ (20)

and using that bare matrix elements are $\mu$ independent, one obtains the
usual DGLAP equations for $f(\mu)$ by taking the $\mu$ derivative of (17). In
analogy, one obtains the homogeneous double DGLAP equation

$\displaystyle\frac{d}{d\log\mu^{2}}\,F(y,\mu)$
$\displaystyle=P(\mu)\underset{1}{\otimes}F(y,\mu)+P(\mu)\underset{2}{\otimes}F(y,\mu)$
(21)

from (18). Taking the $\mu$ derivative of (19), one obtains the inhomogeneous
double DGLAP equation

$\displaystyle\frac{d}{d\log\mu^{2}}\,F(\Delta,\mu)$
$\displaystyle=P(\mu)\underset{1}{\otimes}F(\Delta,\mu)+P(\mu)\underset{2}{\otimes}F(\Delta,\mu)+P_{s}(\mu)\underset{12}{\otimes}f(\mu)$
(22)

with the $1\to 2$ evolution kernel

$\displaystyle P_{s}(\mu)$
$\displaystyle=\biggl{[}\biggl{(}\frac{d}{d\log\mu^{2}}\,Z_{s}(\mu)\biggr{)}-P(\mu)\underset{1}{\otimes}Z_{s}(\mu)-P(\mu)\underset{2}{\otimes}Z_{s}(\mu)\biggr{]}\underset{12}{\otimes}Z^{-1}(\mu)\,,$
(23)

where $Z^{-1}$ was defined in the previous subsection. Making parton indices
explicit, we have DPDs $F_{a_{1}a_{2}}$ and PDFs $f_{a}$. Correspondingly, the
renormalisation factors and evolution kernels for PDFs read $Z_{ab}$ and
$P_{ab}$ as usual. For the quantities associated with $1\to 2$ splitting in
DPDs we write $Z_{a_{1}a_{2},a_{0}}$ and $P_{a_{1}a_{2},a_{0}}$, omitting the
subscript $s$ on $Z$ and $P$ when parton indices are written out.

###  2.3 Splitting kernels at higher order

Let us now see how renormalisation affects the $1\to 2$ splitting kernels $V$
and $W$. For unrenormalised DPDs at small $y$ or large $\Delta$, the splitting
mechanism gives factorisation formulae

$\displaystyle F_{B}(y)$ $\displaystyle=\frac{\Gamma(1-\epsilon)}{(\pi
y^{2})^{1-\epsilon}}\;V_{B}(y)\underset{12}{\otimes}f_{B}\,,$ $\displaystyle
F_{B}(\Delta)$ $\displaystyle=W_{B}(\Delta)\underset{12}{\otimes}f_{B}\,,$
(24)

with bare kernels $V_{B}$ and $W_{B}$, where the prefactor of $V_{B}(y)$ in
$F_{B}(y)$ is taken out for convenience.222Compared with [43], this prefactor
includes an additional $\Gamma(1-\epsilon)$. According to (17) to (19), we
then have

$\displaystyle F(y,\mu)$ $\displaystyle=\frac{\Gamma(1-\epsilon)}{(\pi
y^{2})^{1-\epsilon}}\;V(y,\mu)\underset{12}{\otimes}f(\mu)\,,$ $\displaystyle
F(\Delta,\mu)$ $\displaystyle=W(\Delta,\mu)\underset{12}{\otimes}f(\mu)$ (25)

with kernels

$\displaystyle V(y,\mu)$
$\displaystyle=Z(\mu)\underset{1}{\otimes}Z(\mu)\underset{2}{\otimes}V_{B}(y)\underset{12}{\otimes}Z^{-1}(\mu)\,,$
$\displaystyle W(\Delta,\mu)$
$\displaystyle=Z(\mu)\underset{1}{\otimes}Z(\mu)\underset{2}{\otimes}W_{B}(\Delta)\underset{12}{\otimes}Z^{-1}(\mu)+Z_{s}(\mu)\underset{12}{\otimes}Z^{-1}(\mu)$
(26)

for renormalised distributions.

#####  $\overline{\text{MS}}$ implementation and coupling renormalisation.

For the bare momentum space kernel, we have a perturbative expansion

$\displaystyle W_{B}(\Delta)$
$\displaystyle=\sum_{n}\biggl{(}\frac{\alpha_{0}}{2\pi}\biggr{)}^{n}\,W_{0}^{(n)}(\Delta)=\sum_{n}\biggl{(}\frac{\mu^{2\epsilon}}{S_{\epsilon}}\,\frac{\alpha_{s}(\mu)}{2\pi}\,Z_{\alpha}(\mu)\biggr{)}^{n}\,W_{0}^{(n)}(\Delta)$
$\displaystyle=\sum_{n}a_{s}^{n}(\mu)\,\biggl{(}\frac{\mu}{\Delta}\biggr{)}^{2\epsilon
n}Z^{n}_{\alpha}(\mu)\;W_{B}^{(n)}\,,$ (27)

where $\alpha_{0}$ is the bare coupling. In the last step, we used the
abbreviation (6) and introduced coefficients

$\displaystyle W_{B}^{(n)}(x_{1},x_{2},\epsilon)$
$\displaystyle=\Delta^{2\epsilon
n}\,S_{\epsilon}^{-n}\,W_{0}^{(n)}(x_{1},x_{2},\Delta,\epsilon)\,,$ (28)

which are $\Delta$ independent, because
$W_{0}^{(n)}(\Delta)\propto\Delta^{-2\epsilon n}$ for dimensional reasons. As
discussed in [54], we implement the $\overline{\text{MS}}$ prescription by
dividing bare graphs of order $n$ by $S_{\epsilon}^{n}$. The coupling
renormalisation constant $Z_{\alpha}$ is then $1$ plus a sum of pure poles in
$\epsilon$. The standard choice in the literature is $S_{\epsilon}=(4\pi
e^{-\gamma})^{\epsilon}$. An alternative choice
$S_{\epsilon}=(4\pi)^{\epsilon}/\Gamma(1-\epsilon)$ was proposed by Collins
[70], with the statement that for quantities for which ultraviolet divergences
give at most one power of $\epsilon^{-1}$ per loop, the two schemes give
identical results for renormalised quantities at $\epsilon=0$. We shall verify
this explicitly up to order $a_{s}^{2}$ in section 2.3.1.

Up to order $a_{s}$, we have

$\displaystyle Z_{\alpha}$
$\displaystyle=1-\frac{a_{s}}{\epsilon}\,\frac{\beta_{0}}{2}+\mathcal{O}(a_{s}^{2})\,,$
$\displaystyle\beta_{0}$ $\displaystyle=\frac{11}{3}\mskip
1.5muC_{A}-\frac{4}{3}\,T_{F}\mskip 1.5mun_{F}\,,$ (29)

where for brevity we do not display the $\mu$ dependence of $Z_{\alpha}$ and
$a_{s}$. Here $T_{F}=1/2$, and $n_{F}$ is the number of active quarks. This
gives

$\displaystyle\frac{da_{s}}{d\log\mu^{2}}$
$\displaystyle=\frac{\beta(a_{s})}{2\pi}=-\frac{\beta_{0}}{2}\,a_{s}^{2}+\mathcal{O}(a_{s}^{3})$
(30)

for the running coupling in $4$ dimensions, whereas the renormalisation group
derivative in $4-2\epsilon$ dimensions reads

$\displaystyle\frac{d}{d\log\mu^{2}}$
$\displaystyle=\frac{\partial}{\partial\log\mu^{2}}+\biggl{[}\frac{\beta(a_{s})}{2\pi}-\epsilon\mskip
1.5mua_{s}\biggr{]}\,\frac{\partial}{\partial a_{s}}\,.$ (31)

Expanding (2.3) up to order $a_{s}^{2}$, we obtain

$\displaystyle W_{B}(\Delta)$
$\displaystyle=a_{s}\,\biggl{(}\frac{\mu}{\Delta}\biggr{)}^{2\epsilon}W_{B}^{(1)}+a_{s}^{2}\,\biggl{(}\frac{\mu}{\Delta}\biggr{)}^{4\epsilon}W_{B}^{(2)}-a_{s}^{2}\,\biggl{(}\frac{\mu}{\Delta}\biggr{)}^{2\epsilon}\frac{\beta_{0}}{2\epsilon}\,W_{B}^{(1)}+\mathcal{O}(a_{s}^{3})$
(32)

for the bare $1\to 2$ splitting kernel.

##### Renormalisation factors and splitting functions.

The perturbative expansion of the renormalisation functions $Z$ and $Z_{s}$
reads

$\displaystyle Z(x,\epsilon,\mu)$
$\displaystyle=\delta(1-x)+\sum_{n=1}^{\infty}a_{s}^{n}(\mu)\,Z^{(n)}(x,\epsilon)\,,$
$\displaystyle Z_{s}(x_{1},x_{2},\epsilon,\mu)$
$\displaystyle=\sum_{n=1}^{\infty}a_{s}^{n}(\mu)\,Z_{s}^{(n)}(x_{1},x_{2},\epsilon)\,.$
(33)

Note that $Z_{s}$ has no tree level term. With the implementation of the
$\overline{\text{MS}}$ scheme discussed below (2.3), all counterterms are pure
poles in $\epsilon$, and we have

$\displaystyle\bigl{[}Z^{(1)}\bigr{]}_{-1}$ $\displaystyle=-P^{(0)}\,,$
$\displaystyle\bigl{[}Z^{(1)}_{s}\bigr{]}_{-1}$
$\displaystyle=-P_{s}^{(0)}\,,$
$\displaystyle\bigl{[}Z^{(2)}_{s}\bigr{]}_{-1}$
$\displaystyle=-P_{s}^{(1)}\big{/}\,2\,,$ (34)

where the evolution kernels are expanded as

$\displaystyle P(x,\mu)$
$\displaystyle=\sum_{n=0}^{\infty}a_{s}^{n+1}(\mu)\,P^{(n)}(x)\,,$
$\displaystyle P_{s}(x_{1},x_{2},\mu)$
$\displaystyle=\sum_{n=0}^{\infty}a_{s}^{n+1}(\mu)\,P_{s}^{(n)}(x_{1},x_{2})\,.$
(35)

The residues of $1/\epsilon$ in (34) are obtained by inserting (31) into (20)
and (23) and then comparing the $\mathcal{O}(\epsilon^{0})$ terms on both
sides, taking into account that the evolution kernels are independent of
$\epsilon$, see e.g. [54]. Comparing the $\mathcal{O}(\epsilon^{-1})$ terms,
one finds that $[Z]_{-2}$ and $[Z_{s}]_{-2}$ start at order $a_{s}^{2}$. In
particular, one has

$\displaystyle\bigl{[}Z_{s}^{(2)}\bigr{]}_{-2}$
$\displaystyle=\frac{1}{2}\,\biggl{(}\frac{\beta_{0}}{2}\,P_{s}^{(0)}+P^{(0)}\underset{1}{\otimes}P_{s}^{(0)}+P^{(0)}\underset{2}{\otimes}P_{s}^{(0)}+P_{s}^{(0)}\underset{12}{\otimes}P^{(0)}\biggr{)}\,.$
(36)

At LO, the inhomogeneous evolution kernel has the form

$\displaystyle P_{s}^{(0)}(x_{1},x_{2})$
$\displaystyle=\delta(1-x_{1}-x_{2})\,P_{s}^{(0)}(x_{1})$ (37)

with a kinematic constraint due to the fact that the splitting process gives
exactly two final state partons at this order. The one-variable kernel
$P_{s}(x)$ is equal to $P^{(0)}(x)$, except that it has no $\delta(1-x)$ terms
and no plus prescription on $1/(1-x)$ factors. The convolutions in (36) hence
turn into ordinary products:

$\displaystyle P^{(0)}\underset{1}{\otimes}P_{s}^{(0)}$
$\displaystyle=P^{(0)}\biggl{(}\frac{x_{1}}{1-x_{2}}\biggr{)}\,\frac{P_{s}^{(0)}(1-x_{2})}{1-x_{2}}\,,$
$\displaystyle P^{(0)}\underset{2}{\otimes}P_{s}^{(0)}$
$\displaystyle=P^{(0)}\biggl{(}\frac{x_{2}}{1-x_{1}}\biggr{)}\,\frac{P_{s}^{(0)}(x_{1})}{1-x_{1}}\,,$
$\displaystyle P_{s}^{(0)}\underset{12}{\otimes}P^{(0)}$
$\displaystyle=P_{s}^{(0)}\biggl{(}\frac{x_{1}}{x_{1}+x_{2}}\biggr{)}\,\frac{P^{(0)}(x_{1}+x_{2})}{x_{1}+x_{2}}\,.$
(38)

Note that for flavour diagonal transitions these expressions still contain
distributions, which will be made explicit in (5) and (5).

####  2.3.1 Momentum space kernels

We expand the renormalised kernel (2.3) in $4-2\epsilon$ dimensions as

$\displaystyle W(\Delta,\mu,\epsilon)$
$\displaystyle=\sum_{n=1}^{\infty}a_{s}^{n}(\mu)\,W^{(n)}(\Delta,\mu,\epsilon)\,,$
$\displaystyle W^{(n)}(\Delta,\mu,\epsilon)$
$\displaystyle=\sum_{m=0}^{n}\biggl{(}\frac{\mu}{\Delta}\biggr{)}^{2\epsilon
m}\,W^{(n,m)}(\epsilon)\,.$ (39)

The coefficient $W^{(n,n)}$ equals the $n$th order term $\smash{W_{B}^{(n)}}$
of the bare kernel $W_{B}$ in (2.3), the coefficients $W^{(n,m)}$ with $0<m<n$
are products of lower order terms $\smash{W_{B}^{(m)}}$ with renormalisation
counterterms for the twist-two operators or the QCD coupling, whilst
$W^{(n,0)}$ comes from the counterterm for the splitting singularity.
Expanding $(\mu/\Delta)^{2\epsilon m}$ in $\epsilon$, we obtain the kernel in
$4$ dimensions,

$\displaystyle
W^{(n)}(\Delta,\mu)\underset{\text{def}}{=}W^{(n)}(\Delta,\mu,\epsilon=0)$
$\displaystyle=\sum_{k=0}^{n}\biggl{(}\log\frac{\mu^{2}}{\Delta^{2}}\biggr{)}^{k}\,W^{[n,k]}$
(40)

with

$\displaystyle W^{[n,0]}$
$\displaystyle=\sum_{m=0}^{n}\bigl{[}W^{(n,m)}\bigr{]}_{0}\,,$ $\displaystyle
W^{[n,k]}$
$\displaystyle=\sum_{m=1}^{n}\frac{m^{k}}{k!}\,\bigl{[}W^{(n,m)}\bigr{]}_{-k}\text{
\leavevmode\nobreak\ for\leavevmode\nobreak\ }k\geq 1\,.$ (41)

The absence of poles $1/\epsilon^{j}$ in $W^{(n)}$ yields the finiteness
conditions

$\displaystyle 0$
$\displaystyle=\sum_{m=0}^{n}\bigl{[}W^{(n,m)}\bigr{]}_{-j}\,,$ $\displaystyle
0$ $\displaystyle=\sum_{m=1}^{n}m^{k}\,\bigl{[}W^{(n,m)}\bigr{]}_{-k-j}\text{
\leavevmode\nobreak\ for\leavevmode\nobreak\ }k\geq 1\,,$ (42)

where $j\geq 1$ in all cases.

At LO we have

$\displaystyle W^{(1,1)}$ $\displaystyle=W_{B}^{(1)}\,,$ $\displaystyle
W^{(1,0)}$ $\displaystyle=Z_{s}^{(1)}=-\epsilon^{-1}\mskip
1.5muP_{s}^{(0)}\,,$ (43)

where in the last step we used (34). $W^{(1)}$ has at most single poles, and
the first condition in (42) implies

$\displaystyle\bigl{[}W_{B}^{(1)}\bigr{]}_{-1}=P_{s}^{(0)}\,.$ (44)

The coefficients of the kernel (40) in $4$ dimensions thus read

$\displaystyle W^{[1,0]}$
$\displaystyle=\bigl{[}W^{(1,1)}\bigr{]}_{0}=\bigl{[}W_{B}^{(1)}\bigr{]}_{0}\,,$
$\displaystyle W^{[1,1]}$
$\displaystyle=\bigl{[}W^{(1,1)}\bigr{]}_{-1}=P_{s}^{(0)}\,.$ (45)

At NLO we have

$\displaystyle W^{(2,2)}$ $\displaystyle=W_{B}^{(2)}\,,\phantom{\frac{1}{1}}$
$\displaystyle W^{(2,1)}$
$\displaystyle=-\frac{\beta_{0}}{2\epsilon}\,W_{B}^{(1)}+Z^{(1)}\underset{1}{\otimes}W_{B}^{(1)}+Z^{(1)}\underset{2}{\otimes}W_{B}^{(1)}-W_{B}^{(1)}\underset{12}{\otimes}Z^{(1)}$
$\displaystyle=-\frac{1}{\epsilon}\,\biggl{(}\frac{\beta_{0}}{2}\,W_{B}^{(1)}+P^{(0)}\underset{1}{\otimes}W_{B}^{(1)}+P^{(0)}\underset{2}{\otimes}W_{B}^{(1)}-W_{B}^{(1)}\underset{12}{\otimes}P^{(0)}\biggr{)}\,,$
$\displaystyle W^{(2,0)}$
$\displaystyle=Z_{s}^{(2)}-Z_{s}^{(1)}\underset{12}{\otimes}Z^{(1)}\phantom{\frac{1}{1}}$
$\displaystyle=-\frac{1}{2\epsilon}\,P_{s}^{(1)}+\frac{1}{\epsilon^{2}}\,\bigl{[}Z_{s}^{(2)}\bigr{]}_{-2}-\frac{1}{\epsilon^{2}}\,P_{s}^{(0)}\underset{12}{\otimes}P^{(0)}\,,$
(46)

where in the second steps we used (34) again. The finiteness conditions (42)
now read

$\displaystyle 0$
$\displaystyle=\bigl{[}W^{(2,2)}\bigr{]}_{-1}+\bigl{[}W^{(2,1)}\bigr{]}_{-1}+\bigl{[}W^{(2,0)}\bigr{]}_{-1}\,,$
$\displaystyle 0$
$\displaystyle=\bigl{[}W^{(2,2)}\bigr{]}_{-2}+\bigl{[}W^{(2,1)}\bigr{]}_{-2}+\bigl{[}W^{(2,0)}\bigr{]}_{-2}\,,$
$\displaystyle 0$
$\displaystyle=2\bigl{[}W^{(2,2)}\bigr{]}_{-2}+\bigl{[}W^{(2,1)}\bigr{]}_{-2}\,.$
(47)

Using (36) and the fact the double poles in $W^{(2,1)}$ are obtained by
replacing $W_{B}^{(1)}$ with $\epsilon^{-1}\mskip 1.5muP_{s}^{(0)}$, we get

$\displaystyle\bigl{[}W_{B}^{(2)}\bigr{]}_{-2}$
$\displaystyle=\frac{1}{2}\,\biggl{(}\frac{\beta_{0}}{2}\,P_{s}^{(0)}+P^{(0)}\underset{1}{\otimes}P_{s}^{(0)}+P^{(0)}\underset{2}{\otimes}P_{s}^{(0)}-P_{s}^{(0)}\underset{12}{\otimes}P^{(0)}\biggr{)}\,.$
(48)

This relation fixes the double pole of $W_{B}^{(2)}$ in terms of the LO
kernels and thus serves as a cross check of the two-loop calculation. For the
single poles, we obtain

$\displaystyle\bigl{[}W_{B}^{(2)}\bigr{]}_{-1}$
$\displaystyle=\frac{1}{2}\mskip
1.5muP_{s}^{(1)}+\biggl{(}\frac{\beta_{0}}{2}\,\bigl{[}W_{B}^{(1)}\bigr{]}_{0}+P^{(0)}\underset{1}{\otimes}\bigl{[}W_{B}^{(1)}\bigr{]}_{0}+P^{(0)}\underset{2}{\otimes}\bigl{[}W_{B}^{(1)}\bigr{]}_{0}-\bigl{[}W_{B}^{(1)}\bigr{]}_{0}\underset{12}{\otimes}P^{(0)}\biggr{)}\,,$
(49)

from which we can extract the NLO evolution kernel $P_{s}^{(1)}$. Collecting
the finite contributions, we get

$\displaystyle W^{[2,0]}$
$\displaystyle=\bigl{[}W_{B}^{(2)}\bigr{]}_{0}-\biggl{(}\frac{\beta_{0}}{2}\,\bigl{[}W_{B}^{(1)}\bigr{]}_{1}+P^{(0)}\underset{1}{\otimes}\bigl{[}W_{B}^{(1)}\bigr{]}_{1}+P^{(0)}\underset{2}{\otimes}\bigl{[}W_{B}^{(1)}\bigr{]}_{1}-\bigl{[}W_{B}^{(1)}\bigr{]}_{1}\underset{12}{\otimes}P^{(0)}\biggr{)}\,,$
$\displaystyle W^{[2,1]}$
$\displaystyle=P_{s}^{(1)}+\frac{\beta_{0}}{2}\,W^{[1,0]}+P^{(0)}\underset{1}{\otimes}W^{[1,0]}+P^{(0)}\underset{2}{\otimes}W^{[1,0]}-W^{[1,0]}\underset{12}{\otimes}P^{(0)}\,,$
$\displaystyle W^{[2,2]}$ $\displaystyle=\frac{1}{2}\mskip
1.5mu\biggl{(}\frac{\beta_{0}}{2}\,P_{s}^{(0)}+P^{(0)}\underset{1}{\otimes}P_{s}^{(0)}+P^{(0)}\underset{2}{\otimes}P_{s}^{(0)}-P_{s}^{(0)}\underset{12}{\otimes}P^{(0)}\biggr{)}\,,$
(50)

where we used (45) to replace $\bigl{[}W_{B}^{(1)}\bigr{]}_{0}$ by
$W^{[1,0]}$. Inserting (45) and (2.3.1) into the factorisation formula (25)
for $F(\Delta,\mu)$ and taking the $\mu$ derivative, one obtains the
inhomogeneous double DGLAP equation (22) up to order $a_{s}^{2}$ as a cross
check.

To evaluate (2.3.1) one needs $W_{B}^{(1)}$ up to $\mathcal{O}(\epsilon)$. The
full expressions for $W_{B}^{(1)}$ in $d=4-2\epsilon$ dimensions are easily
computed along the lines of section 5 in [11] and read

$\displaystyle W^{(1)}_{B}(x_{1},x_{2})$
$\displaystyle=\delta(1-x_{1}-x_{2})\,W^{(1)}_{B}(x_{1})$ (51) with
$\displaystyle\epsilon\mskip 1.5muR_{\epsilon}\,W^{(1)}_{B\mskip-1.5mu,\mskip
1.5mugg,g}(x)$ $\displaystyle=2C_{A}\,\biggl{[}\mskip
1.5mu\frac{x}{1-x}+\frac{1-x}{x}+x(1-x)\mskip 1.5mu\biggr{]}\,,$
$\displaystyle\epsilon\mskip 1.5muR_{\epsilon}\,W^{(1)}_{B\mskip-1.5mu,\mskip
1.5muq\bar{q},g}(x)$
$\displaystyle=T_{F}\,\frac{x^{2}+(1-x)^{2}-\epsilon}{1-\epsilon}\,,$
$\displaystyle\epsilon\mskip 1.5muR_{\epsilon}\,W^{(1)}_{B\mskip-1.5mu,\mskip
1.5muqg,q}(x)$ $\displaystyle=C_{F}\,\biggl{[}\mskip
1.5mu\frac{1+x^{2}}{1-x}-\epsilon(1-x)\mskip 1.5mu\biggr{]}\,.$ (52)

Here we introduced the factor

$\displaystyle R_{\epsilon}$
$\displaystyle=\frac{S_{\epsilon}}{(4\pi)^{\epsilon}}\,\frac{\Gamma(1-2\epsilon)}{\Gamma(1+\epsilon)\,\Gamma^{2}(1-\epsilon)}=1+\mathcal{O}(\epsilon^{2})$
(53)

with $S_{\epsilon}$ specified below (28). Since $W_{B}^{(1)}(x_{1},x_{2})$ has
the same kinematic constraint on $x_{1}$ and $x_{2}$ as $P_{s}^{(0)}$ in (37),
the convolutions in (2.3.1) are products as in (2.3), with $P_{s}^{(0)}$
replaced as appropriate.

The square of the tree-level graphs for the splitting $a_{0}\to a_{1}a_{2}$
that give rise to (2.3.1) appears in many higher-order calculations and has
been computed in many papers before. Note however that $W_{B}$ is computed for
transverse momenta of $a_{1}$ and $a_{2}$ that differ by
$\pm\boldsymbol{\Delta}$ between the amplitude and its conjugate (see section
4.1). This is typically not the case in other contexts. Nevertheless, the
expressions on the r.h.s. of (2.3.1) agree with the expressions for $P^{\mskip
1.5mun\neq 4}$ given in section 3 of [72].

##### Equivalence of $\overline{\text{MS}}$ scheme implementations.

Let us show that the two choices for $S_{\epsilon}$ in the implementation of
the $\overline{\text{MS}}$ scheme specified below (28) give the same kernels
$W^{[n,k]}$ in $4$ dimensions at LO and NLO. The choices coincide in the
constant and the first-order term $[S_{\epsilon}]_{\mskip 1.5mu1}$ of the
Taylor expansion of $S_{\epsilon}$ around $\epsilon=0$, but they differ in the
second-order term $[S_{\epsilon}]_{\mskip 1.5mu2}$. Because $W_{B}^{(1)}$ has
only a single pole in $\epsilon$, the renormalised LO kernels are identical
for the two choices, as are the associated renormalisation factors $Z^{(1)}$
and $Z_{s}^{(1)}$. Since $W_{B}^{(2)}$ contains at most double poles in
$\epsilon$, the NLO renormalisation factor $Z_{s}^{(2)}$ is the same for the
two choices as well. The only remaining dependence on $[S_{\epsilon}]_{2}$ in
the renormalised two loop kernel $W^{(2)}(\Delta,\mu)$ at $\epsilon=0$ can
thus come from the terms

$\displaystyle\bigl{[}W^{(2,2)}\bigr{]}_{-2}\,\bigl{[}S_{\epsilon}^{-2}\bigr{]}_{2}+\bigl{[}W^{(2,1)}\bigr{]}_{-2}\,\bigl{[}S_{\epsilon}^{-1}\bigr{]}_{2}$
$\displaystyle\quad=\Bigl{(}3\bigl{[}W^{(2,2)}\bigr{]}_{-2}+\bigl{[}W^{(2,1)}\bigr{]}_{-2}\mskip
1.5mu\Bigr{)}\;\Bigl{(}\bigl{[}S_{\epsilon}\bigr{]}_{1}\Bigr{)}^{2}-\Bigl{(}2\bigl{[}W^{(2,2)}\bigr{]}_{-2}+\bigl{[}W^{(2,1)}\bigr{]}_{-2}\mskip
1.5mu\Bigr{)}\;\bigl{[}S_{\epsilon}\bigr{]}_{2}\,,$ (54)

where the dependence on $S_{\epsilon}$ follows from (28) and is due to the
terms $W^{(2)}_{B}$ and $W^{(1)}_{B}$ in (2.3.1). Thanks to the last relation
in (2.3.1), the dependence on $[S_{\epsilon}]_{\mskip 1.5mu2}$ cancels in this
expression, which completes our argument.

####  2.3.2 Position space kernels

The bare kernels $V_{B}$ and $W_{B}$ are related by a Fourier transform
according to (4). Using equation (E.1) in [73], one gets for the $a_{s}^{n}$
term in (2.3)

$\displaystyle\int\frac{d^{2-2\epsilon}\boldsymbol{\Delta}}{(2\pi)^{2-2\epsilon}}\;e^{-i\boldsymbol{\Delta}\boldsymbol{y}}\,\biggl{(}\frac{\mu}{\Delta}\biggr{)}^{2\epsilon
n}W_{B}^{(n)}$ $\displaystyle=\frac{\Gamma(1-\epsilon)}{(\pi
y^{2})^{1-\epsilon}}\,\biggl{(}\frac{y\mu}{b_{0}}\biggr{)}^{2\epsilon
n}n\epsilon\,T_{\epsilon,n}\,W_{B}^{(n)}$ (55)

with $b_{0}$ given in (6) and

$\displaystyle T_{\epsilon,n}$ $\displaystyle=\frac{\Gamma(1-\epsilon-\epsilon
n)}{\Gamma(1+\epsilon
n)\,\Gamma(1-\epsilon)}\,e^{-2n\gamma\epsilon}=1+\zeta_{2}\,n\mskip
1.5mu\epsilon^{2}+\frac{\zeta_{3}}{3}\,(2n^{3}+3n^{2}+3n)\,\epsilon^{3}+\mathcal{O}(\epsilon^{4})\,,$
(56)

where $\zeta_{n}$ denotes the Riemann $\zeta$ function evaluated at integer
argument $n$. We define333Here and in (59) we use a convention that differs
from equations (3.15) and (3.16) in [43].

$\displaystyle V(y,\mu,\epsilon)$
$\displaystyle=\sum_{n=1}^{\infty}a_{s}^{n}(\mu)\,V^{(n)}(y,\mu,\epsilon)\,,$
$\displaystyle V^{(n)}(y,\mu,\epsilon)$
$\displaystyle=\sum_{m=1}^{n}\biggl{(}\frac{y\mu}{b_{0}}\biggr{)}^{2\epsilon
m}\,V^{(n,m)}(\epsilon)\,.$ (57)

In contrast to (39) there is no $m=0$ term here, because there is no splitting
singularity in $y$ space. By contrast, the counterterms for the twist-two
operators and the QCD coupling are common to momentum and position space. One
thus has

$\displaystyle V^{(n,m)}(\epsilon)$
$\displaystyle=m\epsilon\,T_{\epsilon,m}\,W^{(n,m)}(\epsilon)\,.$ (58)

For the kernel at $\epsilon=0$ we write

$\displaystyle
V^{(n)}(y,\mu)\underset{\text{def}}{=}V^{(n)}(y,\mu,\epsilon=0)$
$\displaystyle=\sum_{k=0}^{n-1}\biggl{(}\log\frac{y^{2}\mu^{2}}{b_{0}^{2}}\biggr{)}^{k}\,V^{[n,k]}\,.$
(59)

The coefficients $V^{[n,k]}$ and $V^{(n,m)}$ obey exactly the same relations
as $W^{[n,k]}$ and $W^{(n,m)}$ in (41) and (42), except that the sum over $m$
always starts at $m=1$.

At LO and NLO, the relation (58) implies

$\displaystyle\bigl{[}V^{(n,m)}\bigr{]}_{-k}$
$\displaystyle=\bigl{[}m\epsilon\,W^{(n,m)}\bigr{]}_{-k}=m\,\bigl{[}W^{(n,m)}\bigr{]}_{-k-1}$
$\displaystyle\text{ \leavevmode\nobreak\ for\leavevmode\nobreak\ }k\geq
0\text{ \leavevmode\nobreak\ and\leavevmode\nobreak\ }n=1,2\,.$ (60)

The kernel $V^{(1)}$ is finite and $V^{(2)}$ has at most single poles, so that
one has only one finiteness condition,

$\displaystyle 0$
$\displaystyle=\bigl{[}V^{(2,2)}\bigr{]}_{-1}+\bigl{[}V^{(2,1)}\bigr{]}_{-1}\,.$
(61)

This is satisfied thanks to (60) for $k=1$ and $n=2$, together with the third
relation in (2.3.1). For the finite part of the kernel, the relations (60)
imply

$\displaystyle V^{[n,k-1]}$
$\displaystyle=\sum_{m=1}^{n}\frac{m^{k}}{(k-1)!}\,\bigl{[}W^{(n,m)}\bigr{]}_{-k}=k\sum_{m=0}^{n}\frac{m^{k}}{k!}\,\bigl{[}W^{(n,m)}\bigr{]}_{-k}$
$\displaystyle=k\,W^{[n,k]}\hskip 50.00008pt\text{ \leavevmode\nobreak\
for\leavevmode\nobreak\ }1\leq k\leq n\text{ \leavevmode\nobreak\
and\leavevmode\nobreak\ }n=1,2\,,$ (62)

where we used the relation (41) and its analogue for $V^{[n,k]}$. This gives
the following simple relation between the kernels in 4 dimensions:

$\displaystyle V^{(n)}(y,\mu)$
$\displaystyle=\frac{\partial}{\partial\log\mu^{2}}\;W^{(n)}(\Delta=b_{0}/y,\mu)\hskip
30.00005pt\text{ for }n=1,2\,.$ (63)

Explicitly, we have

$\displaystyle V^{[1,0]}$
$\displaystyle=W^{[1,1]}=P_{s}^{(0)}\,,\phantom{\frac{1}{1}}$ $\displaystyle
V^{[2,0]}$
$\displaystyle=W^{[2,1]}=P_{s}^{(1)}+\frac{\beta_{0}}{2}\,\bigl{[}W_{B}^{(1)}\bigr{]}_{0}+P^{(0)}\underset{1}{\otimes}\bigl{[}W_{B}^{(1)}\bigr{]}_{0}+P^{(0)}\underset{2}{\otimes}\bigl{[}W_{B}^{(1)}\bigr{]}_{0}-\bigl{[}W_{B}^{(1)}\bigr{]}_{0}\underset{12}{\otimes}P^{(0)}\,,$
$\displaystyle V^{[2,1]}$
$\displaystyle=2W^{[2,2]}=\frac{\beta_{0}}{2}\,P_{s}^{(0)}+P^{(0)}\underset{1}{\otimes}P_{s}^{(0)}+P^{(0)}\underset{2}{\otimes}P_{s}^{(0)}-P_{s}^{(0)}\underset{12}{\otimes}P^{(0)}\,.$
(64)

Inserting this into the factorisation formula (25) for $F(y,\mu)$ and taking
the $\mu$ derivative, one correctly obtains the homogeneous double DGLAP
equation (21) up to order $a_{s}^{2}$.

##### Higher orders.

One may wonder whether (2.3.2) and thus (63) holds at all orders $n$. Consider
therefore $n=3$, where the $\mathcal{O}(\epsilon^{2})$ term of
$T_{\epsilon,n}$ contributes for the first time, namely to the finite part of
$V^{(3)}$ but not to its poles. With (56) and (60) one has

$\displaystyle V^{[3,0]}$
$\displaystyle=\sum_{m=1}^{3}\bigl{[}V^{(3,m)}\bigr{]}_{0}=\sum_{m=1}^{3}m\,\bigl{[}W^{(3,m)}\bigr{]}_{-1}+\zeta_{2}\,\sum_{m=1}^{3}m^{2}\,\bigl{[}W^{(3,m)}\bigr{]}_{-3}\,.$
(65)

The sum multiplying $\zeta_{2}$ is zero due to the constraint (42) for
$k=2,j=1$. We thus find that (2.3.2) remains valid for $n=3$. One order
higher, we have

$\displaystyle V^{[4,0]}$
$\displaystyle=\sum_{m=1}^{4}\bigl{[}V^{(4,m)}\bigr{]}_{0}=\sum_{m=1}^{4}m\bigl{[}W^{(4,m)}\bigr{]}_{-1}+\zeta_{2}\,\sum_{m=1}^{4}m^{2}\,\bigl{[}W^{(4,m)}\bigr{]}_{-3}$
$\displaystyle\quad+\frac{\zeta_{3}}{3}\,\sum_{m=1}^{4}(2m^{4}+3m^{2}+3m)\,\bigl{[}W^{(4,m)}\bigr{]}_{-4}\,.$
(66)

According to (42), the sum over $m^{k}\,\bigl{[}W^{(4,m)}\bigr{]}_{-4}$ is
zero for $k=1,2,3$ but not for $k=4$. With (41) we get

$\displaystyle V^{[4,0]}$ $\displaystyle=W^{[4,1]}+16\mskip
1.5mu\zeta_{3}\,W^{[4,4]}$ (67)

and thus find that (2.3.2) is no longer valid for $n=4$.

##  3 Matching between momentum and position space DPDs

In the formalism developed in [43], DPS cross sections are calculated using
position space DPDs. On the other hand, one of the few theoretical constraints
we have on the form of DPDs is given by the sum rules introduced in [53],
which are formulated for momentum space DPDs at $\boldsymbol{\Delta=0}$. Due
to the splitting singularity, position and momentum space DPDs are defined
with a different ultraviolet renormalisation. As a consequence, they are not
simply related by a Fourier transform but by a short-distance matching
formula, which can be calculated in perturbation theory and was derived at LO
in [43]. We now generalise this matching to higher orders in $a_{s}$.

Consider a cut-off regularised momentum space DPD as defined in [43]:

$\displaystyle
F_{\Phi\mskip-1.5mu,\,a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{\Delta},\mu,\nu)$
$\displaystyle=\int
d^{2}\boldsymbol{y}\;e^{i\boldsymbol{y}\boldsymbol{\Delta}}\;\Phi(y\nu)\,F_{a_{1}a_{2}}(x_{1},x_{2},\boldsymbol{y},\mu)\,,$
(68)

where $\Phi$ is a function satisfying $\Phi(u)\to 1$ for $u\to\infty$ and
$\Phi(u)=O(u^{1+\delta})$ for $u\to 0$ with some $\delta>0$. This function
regulates the splitting divergence at small $y$ that one would encounter when
naively Fourier transforming the physical position space DPD in $d=4$
dimensions. Let us first consider the case where $F_{\Phi}$ is defined with a
hard cut-off $\Phi(u)=\Theta(u-b_{0})$. The difference between the
$\overline{\text{MS}}$ renormalised momentum space DPD (19) and (68) is then
given by

$\displaystyle F(\Delta,\mu)-F_{\Phi}(\Delta,\mu,\nu)=\lim_{\epsilon\to
0}\,\biggl{[}\,\int_{y<b_{0}/\nu}d^{2-2\epsilon}\boldsymbol{y}\;e^{i\boldsymbol{\Delta}\boldsymbol{y}}F(y,\epsilon)+Z_{s}(\epsilon)\underset{12}{\otimes}Z^{-1}(\epsilon)\otimes
f\,\biggr{]}$ $\displaystyle\qquad=\lim_{\epsilon\to
0}\,\biggl{[}\,\int_{y<b_{0}/\nu}d^{2-2\epsilon}\boldsymbol{y}\;F(y,\epsilon)+Z_{s}(\epsilon)\underset{12}{\otimes}Z^{-1}(\epsilon)\otimes
f\,\biggr{]}+\int_{y<b_{0}/\nu}d^{2}\boldsymbol{y}\;\bigl{[}e^{i\boldsymbol{\Delta}\boldsymbol{y}}-1\bigr{]}\,F(y)\,.$
(69)

Here we have explicitly written out the splitting counterterm from (19) and
indicated all $\epsilon$ dependence, using the shorthand notation
$F(y)=F(y,\epsilon=0)$ in the last term. For brevity, the dependence of
distributions and $Z$ factors on $\mu$ is not shown. The relation (3) allows
us to separate the matching into a part at $\Delta=0$ that involves DPDs in
$4-2\epsilon$ dimensions and a $\Delta$ dependent part that involves only DPDs
at the physical point $\epsilon=0$.

In following is understood that the cut-off scale $\nu$ is large enough to
justify replacing $F(y)$ with the perturbative splitting contribution (25), up
to corrections in powers of $\Lambda/\nu$. As discussed in section 3.3 of
[43], the corrections of order $\Lambda/\nu$ arise from twist-three
distributions, which are expected to be small at low $x_{1}$ and $x_{2}$. The
dominant power corrections are then of order $\Lambda^{2}/\nu^{2}$.

We now turn to the case where $\Phi$ is not a hard cut-off. The difference
between two momentum space DPDs defined with different regulator functions
$\Phi_{1}$ and $\Phi_{2}$ is given by

$\displaystyle F_{\Phi_{1}}-F_{\Phi_{2}}$ $\displaystyle=\int
d^{2}\boldsymbol{y}\;e^{i\boldsymbol{\Delta}\boldsymbol{y}}\,\bigl{[}\Phi_{1}(y\nu)-\Phi_{2}(y\nu)\bigr{]}\,F(y)\,.$
(70)

This involves only $F(y)$ at the physical point, because each regulator
$\Phi_{i}$ ensures that integrals over $y$ are finite in the ultraviolet. For
distances $y$ of hadronic size, the expression in square brackets vanishes
because $\Phi_{i}\to 1$ in that limit, so that one can replace $F(y)$ by its
perturbative splitting form (25). Using (59) one can thus express the
difference (70) in terms of $V^{[n,k]}\underset{12}{\otimes}f(\mu)$ and the
integrals

$\displaystyle\int\frac{d(y^{2})}{y^{2}}\,J_{0}(y\Delta)\,\biggl{(}\log\frac{y^{2}\mu^{2}}{b_{0}^{2}}\biggr{)}^{k}\,\bigl{[}\Phi_{1}(y\nu)-\Phi_{2}(y\nu)\bigr{]}\,,$
(71)

where the angular integration has been performed and has given rise to the
Bessel function $J_{0}$. It is thus straightforward to perform a conversion
between DPDs $F_{\Phi_{i}}$ defined with different regulator functions
$\Phi_{i}$. In the following we will therefore limit ourselves to the choice
$\Phi(u)=\Theta(u-b_{0})$.

###  3.1 Matching at zero $\Delta$

Let us first discuss matching at $\Delta=0$, which we wish to write as

$\displaystyle F(\Delta=0,\mu)-F_{\Phi}(\Delta=0,\mu,\nu)$
$\displaystyle=U(\mu,\nu)\underset{12}{\otimes}f(\mu)+\mathcal{O}(\Lambda/\nu)\,.$
(72)

To compute the matching kernel $U(x_{1},x_{2},\mu,\nu)$, we replace
$F(y,\epsilon)$ in the first term on the r.h.s. of (3) by its perturbative
splitting form (25). With the $y$ dependence in (57), this gives rise to the
integrals

$\displaystyle\int_{y<b_{0}/\nu}d^{2-2\epsilon}\boldsymbol{y}\;\frac{\Gamma(1-\epsilon)}{(\pi
y^{2})^{1-\epsilon}}\,\biggl{(}\frac{y\mu}{b_{0}}\biggr{)}^{2\epsilon m}$
$\displaystyle=\frac{1}{m\epsilon}\,\biggl{(}\frac{\mu}{\nu}\biggr{)}^{2\epsilon
m}\,.$ (73)

Expanding $U(x_{1},x_{2},\mu,\nu,\epsilon)$ as

$\displaystyle U(\mu,\nu,\epsilon)$
$\displaystyle=\sum_{n=1}^{\infty}a_{s}^{n}(\mu)\,U^{(n)}(\mu,\nu,\epsilon)\,,$
$\displaystyle U^{(n)}(\mu,\nu,\epsilon)$
$\displaystyle=\sum_{m=0}^{n}\biggl{(}\frac{\mu}{\nu}\biggr{)}^{2\epsilon
m}\,U^{(n,m)}(\epsilon)$ (74)

in analogy to (57), we thus find that the coefficients $U^{(n,m)}$,
$V^{(n,m)}$, and $W^{(n,m)}$ are related as

$\displaystyle U^{(n,m)}$
$\displaystyle=\frac{1}{m\epsilon}\,V^{(n,m)}=T_{\epsilon,m}\,W^{(n,m)}\qquad\text{
\leavevmode\nobreak\ for\leavevmode\nobreak\ }m\geq 1$ (75)

with $T_{\epsilon,m}$ from (56), where in the last step we have used the
relation (58) between $V^{(n,m)}$ and $W^{(n,m)}$. The $m=0$ coefficient
originates from the splitting counterterm in (3) and thus reads

$\displaystyle U^{(n,0)}$ $\displaystyle=W^{(n,0)}\,.$ (76)

In physical dimensions, i.e. at $\epsilon=0$, the matching kernel can again be
written as

$\displaystyle
U^{(n)}(\mu,\nu)\underset{\text{def}}{=}U^{(n)}(\mu,\nu,\epsilon=0)$
$\displaystyle=\sum_{m=0}^{n}\biggl{(}\log\frac{\mu^{2}}{\nu^{2}}\biggr{)}^{m}\,U^{[n,m]}\,,$
(77)

where $U^{[n,k]}$ and $U^{(n,m)}$ obey the same relations as $W^{[n,k]}$ and
$W^{(n,m)}$ in (41). The coefficients $U^{(n,m)}$ fulfil finiteness relations
analogous to those for $W^{(n,m)}$ in (42).

Since $T_{\epsilon,m}=1+\mathcal{O}(\epsilon^{2})$, the pole terms of
$U^{(n)}$ and $W^{(n)}$ coincide at LO and NLO, which ensures the validity of
the finiteness conditions for $U$ and fixes all coefficients of the logarithms
in (77). For the nonlogarithmic term we obtain

$\displaystyle U^{[n,0]}$
$\displaystyle=\sum_{m=0}^{n}\bigl{[}U^{(n,m)}\bigr{]}_{0}=\sum_{m=0}^{n}\bigl{[}W^{(n,m)}\bigr{]}_{0}+\zeta_{2}\,\sum_{m=1}^{n}m\mskip
1.5mu\bigl{[}W^{(n,m)}\bigr{]}_{-2}\text{ \leavevmode\nobreak\
for\leavevmode\nobreak\ }n=1,2,$ (78)

where in the second step we used (75), (76) and (56). The sum multiplying
$\zeta_{2}$ vanishes due to the constraint (42) for $k=1,j=1$. In total, we
thus have

$\displaystyle U^{[n,k]}$ $\displaystyle=W^{[n,k]}$ $\displaystyle\text{
\leavevmode\nobreak\ for\leavevmode\nobreak\ }0\leq k\leq n\text{
\leavevmode\nobreak\ and\leavevmode\nobreak\ }n=1,2.$ (79)

##### Higher orders.

Let us see whether the relation (79) remains valid at higher orders. To this
end, we compute the coefficient

$\displaystyle U^{[3,0]}$
$\displaystyle=\sum_{m=0}^{3}\bigl{[}U^{(3,m)}\bigr{]}_{0}$
$\displaystyle=\sum_{m=0}^{3}\bigl{[}W^{(3,m)}\bigr{]}_{0}+\zeta_{2}\,\sum_{m=1}^{3}m\mskip
1.5mu\bigl{[}W^{(3,m)}\bigr{]}_{-2}+\frac{\zeta_{3}}{3}\,\sum_{m=1}^{3}(2m^{3}+3m+3)\,\bigl{[}W^{(3,m)}\bigr{]}_{-3}$
$\displaystyle=W^{[3,0]}+4\mskip
1.5mu\zeta_{3}\,W^{[3,3]}\,,\phantom{\sum_{m=0}^{0}}$ (80)

where we have used (75) and (42). We thus find that (79) is no longer valid
for $n\geq 3$. We can, however, derive from (75) a relation between the
matching kernels and the position space kernels that is valid to all orders,
namely

$\displaystyle\bigl{[}U^{(n,m)}\bigr{]}_{-k-1}$
$\displaystyle=\frac{1}{m}\,\bigl{[}V^{(n,m)}\bigr{]}_{-k}$
$\displaystyle\text{ \leavevmode\nobreak\ for\leavevmode\nobreak\ }k\geq
0,\leavevmode\nobreak\ m\geq 1,\text{ \leavevmode\nobreak\ and
all\leavevmode\nobreak\ }n,$ (81)

In analogy to (2.3.2), we thus get

$\displaystyle U^{[n,k]}$ $\displaystyle=\frac{1}{k}\,V^{[n,k-1]}$
$\displaystyle\text{ \leavevmode\nobreak\ for\leavevmode\nobreak\ }1\leq k\leq
n\text{ \leavevmode\nobreak\ and all\leavevmode\nobreak\ }n$ (82)

which for $n=1,2$ is of course consistent with (2.3.2) and (79).

###  3.2 Matching at nonzero $\Delta$

For matching at nonzero $\Delta$, we need both terms on the r.h.s. of (3).
Using the perturbative splitting form of $F(y)$ (25) and the expansion (59) of
the kernel $V$, one readily obtains

$\displaystyle F(\Delta,\mu)-F_{\Phi}(\Delta,\mu,\nu)$
$\displaystyle=U(\mu,\nu)\underset{12}{\otimes}f(\mu)+\sum_{n=1}^{\infty}a_{s}^{n}(\mu)\sum_{k=0}^{n-1}I_{k}(\Delta,\nu,\mu)\,V^{[n,k]}\otimes
f(\mu)$ (83)

with corrections of $\mathcal{O}(\Lambda/\nu)$ and coefficients $I_{k}$ given
by the integrals

$\displaystyle I_{k}(\Delta,\nu,\mu)$
$\displaystyle=\int_{0}^{b_{0}^{2}/\nu^{2}}\frac{d(y^{2})}{y^{2}}\,\bigl{[}J_{0}(y\Delta)-1\bigr{]}\,\biggl{(}\log\frac{y^{2}\mu^{2}}{b_{0}^{2}}\biggr{)}^{k}\,,$
(84)

which arise after performing the angular part of the $y$ integration. We make
the $\mu$ dependence explicit by writing this as

$\displaystyle I_{k}(\Delta,\mu,\nu)$
$\displaystyle=\sum_{j=0}^{k}\binom{k}{j}\,I_{j}(\Delta,\nu,\nu)\,\biggl{(}\log\frac{\mu^{2}}{\nu^{2}}\biggr{)}^{j}\,,$
(85)

where

$\displaystyle I_{j}(\Delta,\nu,\nu)$
$\displaystyle=2\int_{0}^{a}\frac{dz}{z}\,\bigl{[}J_{0}(z)-1\bigr{]}\,\biggl{(}\log\frac{z^{2}}{a^{2}}\biggr{)}^{j}=\frac{1}{j+1}\int_{0}^{a}dz\,J_{1}(z)\,\biggl{(}\log\frac{z^{2}}{a^{2}}\biggr{)}^{j+1}$
(86)

with $a=b_{0}\mskip 1.5mu\Delta/\nu$. In the limits $\Delta\ll\nu$ and
$\nu\ll\Delta$, the above equations can be further simplified. In the former
case the limiting form of (85) is easily obtained by Taylor expanding the
Bessel function, yielding

$\displaystyle I_{k}(\Delta,\nu,\mu)$
$\displaystyle\underset{\Delta\ll\nu}{=}-\frac{b_{0}^{2}\mskip
1.5mu\Delta^{2}}{4\mu^{2}}\,\int_{0}^{\mu^{2}/\nu^{2}}dz\;\log^{k}z=-\frac{b_{0}^{2}\mskip
1.5mu\Delta^{2}}{4\nu^{2}}\,\Biggl{[}\biggl{(}\log\frac{\mu^{2}}{\nu^{2}}\biggr{)}^{k}+\sum_{j=0}^{k-1}c_{j}\,\biggl{(}\log\frac{\mu^{2}}{\nu^{2}}\biggr{)}^{j}\Biggr{]}$
(87)

with numerical coefficients $c_{j}$. The limiting behaviour for $\nu\ll\Delta$
is obtained by writing out the polynomial series of $\log^{\mskip
1.5muk+1}(z^{2}/a^{2})=\bigl{[}\mskip
1.5mu\log(\nu^{2}/\Delta^{2})+\log(z^{2}/b_{0}^{2})\mskip
1.5mu\bigr{]}^{\mskip 1.5muk+1}$ and extending the $z$ integration to infinity
in (86). This gives

$\displaystyle I_{j}(\Delta,\nu,\nu)$
$\displaystyle\underset{\nu\ll\Delta}{=}\,\frac{1}{j+1}\Biggl{[}\biggl{(}\log\frac{\nu^{2}}{\Delta^{2}}\biggr{)}^{j+1}+\sum_{i=0}^{j}d_{i}\,\biggl{(}\log\frac{\nu^{2}}{\Delta^{2}}\biggr{)}^{i}\,\Biggr{]}$
(88)

with numerical coefficients $d_{i}$. Plugging this back into (85) and
rearranging the binomial sum, we obtain the following expression for the large
$\Delta$ behaviour:

$\displaystyle I_{k}(\Delta,\nu,\mu)$
$\displaystyle\underset{\nu\ll\Delta}{=}\,\frac{1}{k+1}\Biggl{[}\biggl{(}\log\frac{\mu^{2}}{\Delta^{2}}\biggr{)}^{k+1}-\biggl{(}\log\frac{\mu^{2}}{\nu^{2}}\biggr{)}^{k+1}\,\Biggr{]}+\cdots\,,$
(89)

where the ellipsis denotes terms with fewer powers of logarithms.

To use the matching equation (83) truncated in $a_{s}$, one must avoid large
logarithms at higher orders. For the term with $U(\mu,\nu)$ this requires one
to take $\nu\sim\mu$. Satisfying this condition, one obtains large logarithms
in (89), so that the choice $\nu\ll\Delta$ has to be avoided. On the other
hand, no large logarithms appear for $\Delta\ll\nu$ according to (87). Fixed
order matching is thus possible both for $\Delta\sim\nu$ and for
$\Delta\ll\nu$. The latter choice is in fact unavoidable for matching at the
point $\Delta=0$.

For $j=0,1,2$ a closed form of $I_{j}$ can be given in terms of the
generalised hypergeometric functions ${}_{p}F_{q}$, namely

$\displaystyle I_{0}(\Delta,\nu,\nu)$
$\displaystyle=-\frac{a^{2}}{4}\,{}_{2}F_{3}\biggl{(}1,1;2,2,2;-\frac{a^{2}}{4}\biggr{)}\,,$
$\displaystyle I_{1}(\Delta,\nu,\nu)$
$\displaystyle=\frac{a^{2}}{4}\,{}_{3}F_{4}\biggl{(}1,1,1;2,2,2,2;-\frac{a^{2}}{4}\biggr{)}\,,$
$\displaystyle I_{2}(\Delta,\nu,\nu)$
$\displaystyle=-\frac{a^{2}}{2}\,{}_{4}F_{5}\biggl{(}1,1,1,1;2,2,2,2,2;-\frac{a^{2}}{4}\biggr{)}\,.$
(90)

###  3.3 Scale independence of matching

We now discuss the dependence of the matching equation (83) on the cut-off
scale $\nu$ and its behaviour under renormalisation group evolution. Let us
first verify explicitly that the $\nu$ dependence is the same on the left and
right hand sides. The logarithmic derivative of the l.h.s. is

$\displaystyle-\frac{d}{d\log\nu^{2}}\,F_{\Phi}(\Delta,\mu,\nu)$
$\displaystyle=-\frac{d}{d\log\nu^{2}}\,\pi\int_{b_{0}^{2}/\nu^{2}}^{\infty}dy^{2}\,J_{0}(y\Delta)\,F(y,\mu)$
$\displaystyle=-\frac{\pi\mskip
1.5mub_{0}^{2}}{\nu^{2}}\,J_{0}\biggl{(}\frac{b_{0}\mskip
1.5mu\Delta}{\nu}\biggr{)}\,F\biggl{(}y=\frac{b_{0}}{\nu},\mu\biggr{)}$
$\displaystyle=-J_{0}\biggl{(}\frac{b_{0}\mskip
1.5mu\Delta}{\nu}\biggr{)}\,\sum_{n=1}^{\infty}a_{s}^{n}\,\sum_{k=0}^{n-1}\biggl{(}\log\frac{\mu^{2}}{\nu^{2}}\biggr{)}^{k}\,V^{[n,k]}\underset{12}{\otimes}f(\mu)\,,$
(91)

where in the last step we used the perturbative splitting form of $F(y)$ in
(25) and (59), as is appropriate for $y=b_{0}/\nu$. For the r.h.s. of (83) one
gets

$\displaystyle\frac{d}{d\log\nu^{2}}\,\sum_{n=1}^{\infty}a_{s}^{n}\,\Biggl{\\{}\sum_{k=0}^{n}\biggl{(}\log\frac{\mu^{2}}{\nu^{2}}\biggr{)}^{k}\,U^{[n,k]}+\sum_{k=0}^{n-1}I_{k}(\Delta,\nu,\mu)\,V^{[n,k]}\mskip
1.5mu\Biggr{\\}}\underset{12}{\otimes}f(\mu)$
$\displaystyle\quad=-\sum_{n=1}^{\infty}a_{s}^{n}\,\sum_{k=0}^{n-1}\biggl{(}\log\frac{\mu^{2}}{\nu^{2}}\biggr{)}^{k}\Biggl{\\{}(k+1)\,U^{[n,k+1]}+\biggl{[}J_{0}\biggl{(}\frac{b_{0}\mskip
1.5mu\Delta}{\nu}\biggr{)}-1\biggr{]}\,V^{[n,k]}\mskip
1.5mu\Biggr{\\}}\underset{12}{\otimes}f(\mu)$ (92)

using the integral representation (84) of $I_{k}$. With the relation (82)
between $V$ and $U$, we see that both sides of the matching equation have the
same $\nu$ dependence at each fixed order in $a_{s}$, as it should be. This
statement holds up to power corrections in $\Lambda/\nu$, due to replacing
$F(y)$ with its perturbative splitting form.

Let us now consider the renormalisation scale dependence of (83). Using the
appropriate – homogeneous or inhomogeneous – double DGLAP equations, one
obtains for the l.h.s.

$\displaystyle\frac{d}{d\log\mu^{2}}\,\Bigl{[}F(\Delta,\mu)-F_{\Phi}(\Delta,\mu,\nu)\Bigr{]}=P_{s}\underset{12}{\otimes}f(\mu)$
$\displaystyle\qquad+P\underset{1}{\otimes}\Bigl{[}F(\Delta,\mu)-F_{\Phi}(\Delta,\mu,\nu)\Bigr{]}+P\underset{2}{\otimes}\Bigl{[}F(\Delta,\mu)-F_{\Phi}(\Delta,\mu,\nu)\Bigr{]}\,,$
(93)

where the term with $P_{s}$ arises from the $\overline{\text{MS}}$ momentum
space DPD. On the r.h.s. of (3.3) one can replace the difference of DPDs with
the matching expression (83). Comparing this to the renormalisation group
derivative of the r.h.s. of (83), one finds that the matching equation holds
at all $\mu$ if

  * •

$U(\mu,\nu)\underset{12}{\otimes}f(\mu)$ obeys the inhomogeneous double DGLAP
equation. At LO and NLO, this can be explicitly verified using the evolution
equation for $W(\mu,\nu)\underset{12}{\otimes}f(\mu)$ and the equality of $W$
and $U$ stated in (79),

  * •

the $\Delta$ dependent term

$\displaystyle\sum_{n=1}^{\infty}a_{s}^{n}(\mu)\sum_{k=0}^{n-1}I_{k}(\Delta,\nu,\mu)\,V^{[n,k]}\underset{12}{\otimes}f(\mu)$
(94)

satisfies the homogeneous double DGLAP equation. This follows from the fact
that the splitting form

$\displaystyle\sum_{n=1}^{\infty}a_{s}^{n}(\mu)\sum_{k=0}^{n-1}\biggl{(}\log\frac{y^{2}\mu^{2}}{b_{0}^{2}}\biggr{)}^{k}\,V^{[n,k]}\underset{12}{\otimes}f(\mu)$
(95)

of $F(y,\mu)$ satisfies the homogeneous equation and that
$I_{k}(\Delta,\nu,\mu)$ and $\log^{k}(y^{2}\mu^{2}/b_{0}^{2})$ satisfy the
same differential equation in $\mu$.

As is the case for the perturbative splitting forms of $F(\Delta,\mu)$ and
$F(y,\mu)$, the double DGLAP equations are fulfilled only up to order
$a_{s}^{n}$ if the perturbative series for the splitting or matching kernel is
truncated at order $a_{s}^{n}$.

##  4 Two-loop calculation

The aim of this work is to obtain the kernels $V$, $W$, $P_{s}$, and $U$
introduced in sections 2 and 3 at NLO accuracy. To this end it suffices to
calculate only the momentum space kernel $W$, since from this $V$, $P_{s}$,
and $U$ can readily be obtained as shown in the sections just mentioned. In
this section, we describe the computation of the NLO momentum space kernel
$W_{a_{1}a_{2},b}^{(2)}(\boldsymbol{\Delta})$ in some detail.

To begin with, we apply the factorisation formula (24) to the bare DPD
$F_{B\mskip-1.5mu,\mskip 1.5mua_{1}a_{2}/a_{0}}(\boldsymbol{\Delta})$ of
partons $a_{1}$ and $a_{2}$ in a parton $a_{0}$. Expanding the formula in
$a_{s}$, we obtain for the second order term

$\displaystyle F^{(2)}_{B\mskip-1.5mu,\mskip
1.5mua_{1}a_{2}/a_{0}}(\boldsymbol{\Delta})$
$\displaystyle=\sum_{b}\Bigl{[}W_{B\mskip-1.5mu,\mskip
1.5mua_{1}a_{2},b}^{(2)}(\boldsymbol{\Delta})\underset{12}{\otimes}f^{(0)}_{B\mskip-1.5mu,\mskip
1.5mub/a_{0}}+W_{B\mskip-1.5mu,\mskip
1.5mua_{1}a_{2},b}^{(1)}(\boldsymbol{\Delta})\underset{12}{\otimes}f^{(1)}_{B\mskip-1.5mu,\mskip
1.5mub/a_{0}}\Bigr{]}$ $\displaystyle=W_{B\mskip-1.5mu,\mskip
1.5mua_{1}a_{2},a_{0}}^{(2)}(\boldsymbol{\Delta})\,,$ (96)

where $f_{B\mskip-1.5mu,\mskip 1.5mub/a_{0}}$ denotes the bare PDF of parton
$b$ in parton $a_{0}$. In the second step we used the tree-level expression
$f^{(0)}_{B\mskip-1.5mu,\mskip 1.5mub/a_{0}}(x)=\delta_{b\mskip
1.5mua_{0}}\,\delta(1-x)$ and the fact that $f^{(1)}_{B\mskip-1.5mu,\mskip
1.5mub/a_{0}}=0$ because the corresponding loop integrals do not depend on any
dimensionful scale. We thus obtain the bare two-loop kernel directly from the
two-loop graphs for the bare DPD of partons $a_{1}$ and $a_{2}$ in an on-shell
parton $a_{0}$. We compute with massless quarks and gluons, so that the
relevant graphs have both infrared and ultraviolet divergences, both being
treated by dimensional regularisation.

###  4.1 Channels and graphs

Let us analyse at which order in $a_{s}$ a kernel $W_{a_{1}a_{2},a_{0}}$ is
first nonzero. The conservation of quark and antiquark flavour implies that
the splitting $a_{0}\to a_{1}a_{2}$ starts at

  1. 1.

LO for $g\to gg$, $g\to q\bar{q}$ and $q\to qg$,

  2. 2.

NLO for $g\to qg$, $q\to gg$ and $q_{j}\to q_{j}q_{k}$, $q_{j}\to
q_{j}\bar{q}_{k}$, $q_{j}\to q_{k}\bar{q}_{k}$,

where $j$ and $k$ label quark flavours and may be equal or different. We will
refer to these as “LO channels” and “NLO channels”, respectively. Further
channels at the same orders are obtained by interchanging partons $a_{1}$ and
$a_{2}$ or by charge conjugation, where the kernels for charge conjugated
channels are identical. All other splitting processes start either at NNLO or
N3LO.

Real emission graphs (or “real graphs” for short) are shown in figures 1 and 2
for LO and NLO channels, respectively. Further graphs for the same channels
are obtained by complex conjugation and by interchanging the lines for partons
$a_{1}$ and $a_{2}$ if $a_{1}=a_{2}$. The unobserved parton radiated into the
final state is uniquely determined for given $a_{0}$, $a_{1}$, and $a_{2}$: it
is a gluon for the LO channels and a quark or antiquark for the NLO channels.

For the LO channels, we have also virtual loop graphs (or “virtual graphs” for
short). They are obtained by dressing the LO splitting graphs with one vertex
correction or with one propagator correction for the lines of parton $a_{1}$
or $a_{2}$ (propagator corrections for the massless on-shell parton $a_{0}$
are zero in dimensional regularisation). This is depicted in figure 3.

At this point a brief comment is in order about the flavour structure for
channels where $a_{0},a_{1}$ and $a_{2}$ are only quarks or antiquarks. In
graphs 2(k), 2(l) and 2(n) there are separate fermion lines for quarks $q$ and
$q^{\prime}$, which may or may not have the same flavour. We denote the
corresponding kernels $W_{a_{1}a_{2},a_{0}}$ as specified in the figure. By
contrast, graphs 2(m) and 2(o) are of a valence type in the sense that they
involve a single quark flavour; we denote the corresponding kernels with a
superscript $v$. The complete kernels for a specific flavour transition are
thus obtained as

$\displaystyle W_{q_{j}\bar{q}_{k},q_{i}}$ $\displaystyle=\delta_{jk}\mskip
1.5muW_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}+\delta_{ij}\mskip
1.5muW_{q\bar{q}^{\prime}\\!,q}+\delta_{ij}\mskip
1.5mu\delta_{jk}\,W_{q\bar{q},q}^{v}\,,$ $\displaystyle W_{q_{j}q_{k},q_{j}}$
$\displaystyle=\delta_{jk}\mskip
1.5muW_{q^{\prime}\\!q,q}+\phantom{\delta_{ij}\,}W_{qq^{\prime}\\!,q}+\phantom{\delta_{ij}\,}\delta_{jk}\mskip
1.5muW_{qq,q}^{v}\,,$ (97)

where $W_{q^{\prime}\\!q,q}(x_{1},x_{2})=W_{qq^{\prime}\\!,q}(x_{2},x_{1})$.

(a) LD  (b) UD

(c) UND  (d) T2B

(e) (f)

(g) (h)

(i) (j)

(k) (l)

(m) (n)

(o) (p)

(q)

Figure 1:  Real graphs for LO 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. The topologies “LD”, “UD” etc. are explained in
section 4.2.

(a) (b)

(c) (d)

(e) (f)

(g) (h)

(i) (j)

(k) $W_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}$ (l) $W_{q\bar{q}^{\prime}\\!,q}$

(m) $W_{q\bar{q},q}^{v}$ (n) $W_{qq^{\prime}\\!,q}$

(o) $W_{qq,q}^{v}$

Figure 2:  Real graphs for NLO channels. The association of parton lines is as
in figure 1.

(a) (b)

(c) (d)

(e) (f)

(g) (h)

(i) (j)

(k) (l)

(m)

Figure 3:  Virtual graphs for the amplitude in LO channels. The tree graphs in
the complex conjugate amplitude are not shown.

##### Feynman vs. light-cone gauge.

We performed two independent calculations of the graphs making up the bare
partonic DPDs – one in covariant Feynman gauge, and the other in $A^{+}=0$
light-cone gauge. For cut graphs in Feynman gauge, we use the Feynman rules
given in appendix D of [73]. Not given there are the expressions for the four-
gluon vertex. The expression to be used to the left of the cut can readily be
found in textbooks (see e.g. section 3.2.4 in [70]), whereas to the right of
the cut one needs to take the complex conjugate expression (just as for the
quark-gluon vertex).

In addition to the graphs shown in figures 1, 2, and 3, one gets in Feynman
gauge also graphs with Wilson lines attaching to the active partons, as
already mentioned below (3). These Wilson line graphs can be obtained from the
graphs without Wilson lines following the prescription illustrated in figure
4. Generally speaking, Wilson lines appear in two cases. The first is when an
active parton, $a_{1}$ or $a_{2}$, is a gluon emerging from a three-gluon
vertex on one of the upper legs, as for instance in figures 1(a), 1(b), and
3(a). Such graphs give rise to two kinds of Wilson line graphs as shown in
figure 4(a) and 4(b). In the second case, the active parton is a quark or
antiquark originating from a quark-gluon vertex on one of the upper legs, as
is e.g. the case in figures 1(h), 1(l), 3(f), and 3(g). This leads to only one
Wilson line graph with a gluon attaching to a quark Wilson line, visualised in
figure 4(c) and 4(d).

Finally, in Feynman gauge one has to include the corresponding Fadeev-Popov
ghost version for all graphs containing closed gluon loops, such as in figure
3(a) and 3(d).

(a) (b)

(c) (d)

Figure 4:  Rules for obtaining the Wilson line graphs needed in Feynman gauge
from graphs without Wilson lines. In each of the four panels, the top left
parton is an active one ($a_{1}$ or $a_{2}$). These rules apply to both real
and virtual graphs, and corresponding rules hold on the right of the final
state cut. Note that the three-gluon vertex in the upper panels gives rise to
two different graphs.

##### Kinematics.

Figure 5:  Assignment of transverse momenta in of a real graph. Note that
compared with the symmetric assignment in [43] we have shifted the integration
variables $\boldsymbol{k}_{1}$ and $\boldsymbol{k}_{2}$ such that
$\boldsymbol{\Delta}$ appears only on the r.h.s. of the cut.

For real graphs, both in LO and NLO channels, the general kinematic structure
is illustrated in figure 5. We work in a frame where the incoming parton
$a_{0}$ has plus momentum $p^{+}$ and zero transverse and minus momentum. The
plus momenta $x_{1}\mskip 1.5mup^{+}$ and $x_{2}\mskip 1.5mup^{+}$ of the
active partons $a_{1}$ and $a_{2}$ are equal in the amplitude and its
conjugate, whereas their transverse momenta differ by $\pm\boldsymbol{\Delta}$
as shown in figure 5. Minus momenta are assigned in the same way as the
transverse ones. The parton $a_{3}$ that goes into the final state has plus
momentum $x_{3}\mskip 1.5mup^{+}$, and its momentum is uniquely fixed by the
momenta of the active partons $a_{1}$ and $a_{2}$, i.e.

$\displaystyle x_{3}$ $\displaystyle=1-x_{1}-x_{2}\,,$ $\displaystyle
k_{3}^{-}=-k_{1}^{-}-k_{2}^{-}\,,$
$\displaystyle\boldsymbol{k}_{3}=-\boldsymbol{k}_{1}-\boldsymbol{k}_{2}\,.$
(98)

To obtain collinear momentum space DPDs, we have to integrate the expressions
corresponding to the real graphs in figures 1 and 2 over $k_{1}^{-}$,
$k_{2}^{-}$, and $\Delta^{-}$, which corresponds to having
$z_{1}^{+}=z_{2}^{+}=y^{+}=0$ in the matrix element (2.1). As we are
interested in collinear DPDs, we also have to integrate over the transverse
momenta $\boldsymbol{k}_{1}$ and $\boldsymbol{k}_{2}$, which is tantamount to
the condition $\boldsymbol{z}_{1}=\boldsymbol{z}_{2}=\boldsymbol{0}$ in (2.1).

(a) vertex correction

(b) propagator correction

Figure 6:  Assignment of momenta in generic virtual graphs to the left of the
cut. An analogous assignment holds for loops to the right of the cut, with
$k_{1}$ and $k_{2}$ shifted by $+\Delta$ and $-\Delta$, respectively.

The kinematic structure of the virtual graphs we encounter in the LO channels
is exemplified in figure 6. Due to momentum conservation we now find that
$k_{2}^{-}=-k_{1}^{-}$ and $\boldsymbol{k}_{2}=-\boldsymbol{k}_{1}$, so that
the integrations over $k_{2}^{-}$ and $\boldsymbol{k}_{2}$ become redundant.
However, there now is an additional loop momentum $\ell$ over which we have to
integrate.

###  4.2 Performing the calculation

We now move on to outline the details of our calculation in Feynman and light-
cone gauge. The initial computation of the diagram expressions is done using
FORM [74] or FeynCalc [75].444In practice, FORM was used in the light-cone
gauge calculation, whilst FeynCalc was used for the Feynman gauge calculation.
Of course, either code could have been used for both computations. One then
needs to perform the integrations over phase space and loop momenta. As
mentioned before, there are two types of graph contributing to the bare
partonic DPD at NLO: real emission graphs in which one parton is emitted into
the final state, and virtual loop graphs with a loop in the amplitude or its
conjugate and no emitted particles in the final state. We discuss the
methodology for performing the integrations for the two cases in turn.

####  4.2.1 Real emission diagrams

For the real emission diagrams we have four different topologies. These are
exemplified by figure 1(a), 1(b), 1(c), 1(d), and we refer to them as LD
(lower diagonal), UD (upper diagonal), UND (upper non-diagonal) and T2B (top
to bottom). The four-gluon graphs in figure 1(e), 1(f), and 1(g) can be
considered to fit into the T2B topology, because all the denominator factors
in these diagrams are covered by those in the T2B topology.

For each graph, we begin by performing the integration over minus components
of the loop and phase space momenta. This is most easily done by changing
variables as follows:

$\displaystyle K^{-}$ $\displaystyle=k_{1}^{-}+k_{2}^{-}\,,$ $\displaystyle
k^{-}$ $\displaystyle=(k_{1}^{-}-k_{2}^{-})/2\,,$ $\displaystyle k^{\prime-}$
$\displaystyle=(k_{1}^{-}-k_{2}^{-})/2+\Delta^{-}\,.$ (99)

The integral over $K^{-}$ is performed using the on-shell delta function,
whilst the integrals over $k^{-}$ and $k^{\prime-}$ can be performed using
complex contour integration. For all topologies, one is able to close the
integration contour on only one pole, both in the $k^{-}$ and in the
$k^{\prime-}$ integration. We note in passing that the variable transformation
(99) is also useful in the general proofs of [76, 54].

As all plus components are fixed by the operators in (2.1), it remains to
integrate over $\boldsymbol{k}_{1}$ and $\boldsymbol{k}_{2}$ at fixed
$\boldsymbol{\Delta}$. Every graph may be written in a form where its
denominator factors are contained within the following set (this may require
some relabelling of integration momenta, and use of the fact that the
integrals are invariant under $\boldsymbol{\Delta}\to-\boldsymbol{\Delta}$):

$\displaystyle D_{1}$
$\displaystyle=\dfrac{(\boldsymbol{k}_{1}+\boldsymbol{\Delta})^{2}}{x_{1}}+\dfrac{(\boldsymbol{k}_{2}-\boldsymbol{\Delta})^{2}}{x_{2}}+\dfrac{(\boldsymbol{k}_{1}+\boldsymbol{k}_{2})^{2}}{x_{3}}\,,$
$\displaystyle D_{2}$
$\displaystyle=\dfrac{\boldsymbol{k}_{1}^{2}}{x_{1}}+\dfrac{\boldsymbol{k}_{2}^{2}}{x_{2}}+\dfrac{(\boldsymbol{k}_{1}+\boldsymbol{k}_{2})^{2}}{x_{3}}\,,$
$\displaystyle D_{3}$
$\displaystyle=(\boldsymbol{k}_{1}+\boldsymbol{\Delta})^{2}\,,$ $\displaystyle
D_{4}$ $\displaystyle=\boldsymbol{k}_{2}^{2}\,,\phantom{\frac{x_{1}}{x_{1}}}$
$\displaystyle\tilde{D}_{4}$ $\displaystyle=\boldsymbol{k}_{1}^{2}\,,$
$\displaystyle\tilde{D}_{5}$
$\displaystyle=(\boldsymbol{k}_{1}+\boldsymbol{k}_{2})^{2}\,.$ (100)

Graphs of all topologies have denominators $D_{1}$ and $D_{2}$. In addition LD
graphs have $\tilde{D}_{5}$, UD graphs have $D_{3}$ and $\tilde{D}_{4}$, T2B
graphs have $D_{3}$ and $\tilde{D}_{5}$, and UND graphs have $D_{3}$ and
$D_{4}$.

We express all numerator factors in the graphs in terms of the $D_{i}$ and
$\tilde{D}_{i}$, and then use integration by parts identities [77, 78], as
implemented in LiteRed [79], to express the results for all graphs in terms of
a small set of master integrals. For graphs with the topology UND the master
integrals have the general form

$\displaystyle I_{1}(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4})$
$\displaystyle=\int\dfrac{d^{d-2}\boldsymbol{k}_{1}\,d^{d-2}\boldsymbol{k}_{2}}{\prod_{i=1..4}D_{i}^{\alpha_{i}}}\,.$
(101)

Specifically, we need the integrals

$\displaystyle I_{1}(1,1,0,0)\,,\quad I_{1}(0,1,1,0)\,,\quad
I_{1}(1,1,1,0)\,,$ (102) $\displaystyle I_{1}(1,0,1,1)\,,\quad
I_{1}(1,1,1,1)\,,\quad I_{1}(2,1,1,1)\,,$

as well as integrals related to these by the simple transformation
$x_{1}\leftrightarrow x_{2}$.

For the remaining topologies, the master integrals are of the form

$\displaystyle I_{2}(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5})$
$\displaystyle=\int\dfrac{d^{d-2}\boldsymbol{k}_{1}\,d^{d-2}\boldsymbol{k}_{2}}{\prod_{i=1...3}D_{i}^{\alpha_{i}}\;\prod_{i=4...5}\tilde{D}_{i}^{\alpha_{i}}}\,,$
(103)

and specifically we need

$\displaystyle I_{2}(0,1,1,0,1)\,,\quad I_{2}(1,1,1,1,0)\,,$ (104)

as well as the integrals on the first line of (102) and integrals related to
these by $x_{1}\leftrightarrow x_{2}$.

We apply the method of differential equations [80, 81, 82, 83] to compute
these master integrals. The integrals depend on the external variables
$x_{1}$, $x_{2}$ and $\Delta^{2}$. Since $\Delta^{2}$ is the only dimensionful
external quantity, each master integral only depends on $\Delta^{2}$ via an
overall prefactor of a power of $\Delta^{2}$ that is easily determined by
dimensional analysis. The differential equations in $\Delta^{2}$ thus only
provide trivial information and are not used. We compute the master integrals
by solving the differential equations in $x_{1}$. For a particular master
integral, solving the equation determines the integral up to a “constant” of
integration, which here is actually a function of $x_{2}$. We determine this
constant of integration by computing the leading behaviour of the master
integrals in the vicinity of the line $x_{3}=0$. As a cross-check, we then
verify that the obtained solution satisfies the differential equation in
$x_{2}$.

The system of differential equations has an almost entirely triangular matrix
structure. For the $I_{1}$ integrals the system schematically looks as
follows:

$\displaystyle\frac{\partial}{\partial x_{1}}\begin{bmatrix}I_{1}(1,1,0,0)\\\
I_{1}(0,1,1,0)\\\ I_{1}(1,1,1,0)\\\ I_{1}(1,0,1,1)\\\ I_{1}(1,1,1,1)\\\
I_{1}(2,1,1,1)\end{bmatrix}=\begin{bmatrix}\blacksquare&0&0&0&0&0\\\
0&\blacksquare&0&0&0&0\\\ \blacklozenge&\blacklozenge&\blacksquare&0&0&0\\\
0&\blacklozenge&0&\blacksquare&0&0\\\
\blacklozenge&\blacklozenge&\blacklozenge&\blacklozenge&\blacksquare&\blacksquare\\\
\blacklozenge&\blacklozenge&\blacklozenge&\blacklozenge&\blacksquare&\blacksquare\end{bmatrix}\begin{bmatrix}I_{1}(1,1,0,0)\\\
I_{1}(0,1,1,0)\\\ I_{1}(1,1,1,0)\\\ I_{1}(1,0,1,1)\\\ I_{1}(1,1,1,1)\\\
I_{1}(2,1,1,1)\end{bmatrix}\,,$ (105)

where black squares $\blacksquare$ denote entries of the form $c(x_{1},x_{2})$
while black diamonds $\blacklozenge$ denote entries of the form
$c_{1}(x_{1},x_{2})+c_{2}(x_{1},x_{2})\,P_{x_{1}\leftrightarrow x_{2}}$. Here
$P_{x_{1}\leftrightarrow x_{2}}$ is the operator transforming a given master
integral $I_{1}(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4})$ to the
corresponding master integral with $x_{1}\leftrightarrow x_{2}$. This enables
us to solve most of the system using elementary methods for first-order
differential equations and forward substitution, working downwards from
$I_{1}(1,1,0,0)$ and $I_{1}(0,1,1,0)$. The two-denominator integrals can be
solved exactly in $d=4-2\epsilon$ dimensions, whilst for the remainder we
obtain a series in $\epsilon$ up to the required order. The only instance in
which we have a coupled system of differential equations is for
$I_{1}(1,1,1,1)$ and $I_{1}(2,1,1,1)$. Here we transform the $2\times 2$
system to the canonical basis [84] using Fuchsia [85], which allows us to
obtain a solution for these two integrals as a series in $\epsilon$. The
system of differential equations for the $I_{2}$ integrals is fully
triangular, so that application of elementary methods plus forward
substitution suffices in this case.

The leading behaviour of the integrals for $x_{3}\to 0$ is computed using the
method of regions [86]. Let us introduce the scaling parameter $\lambda\ll 1$,
and say that $x_{3}\sim\lambda$. Then for all master integrals the following
region gives a leading contribution in $\lambda$:

$\displaystyle R_{1}:$ $\displaystyle x_{3}$ $\displaystyle\sim\lambda\,,$
$\displaystyle x_{1},x_{2}$ $\displaystyle\sim 1\,,$
$\displaystyle\boldsymbol{k}_{1}^{2},\boldsymbol{k}_{2}^{2},(\boldsymbol{k}_{1}+\boldsymbol{\Delta})^{2}$
$\displaystyle\sim\Delta^{2}\,,$
$\displaystyle(\boldsymbol{k}_{1}+\boldsymbol{k}_{2})^{2}$
$\displaystyle\sim\lambda\mskip 1.5mu\Delta^{2}\,.$ (106)

For $I_{1}(1,0,1,1)$ only, we identify a further leading region, namely

$\displaystyle R_{2}:$ $\displaystyle x_{3}$ $\displaystyle\sim\lambda\,,$
$\displaystyle x_{1},x_{2}$ $\displaystyle\sim 1\,,$
$\displaystyle\boldsymbol{k}_{1}^{2},\boldsymbol{k}_{2}^{2},(\boldsymbol{k}_{1}+\boldsymbol{\Delta})^{2},(\boldsymbol{k}_{1}+\boldsymbol{k}_{2})^{2}$
$\displaystyle\sim\Delta^{2}\,.\hskip 48.50012pt$ (107)

For each region, we use the appropriate scaling and drop terms in the
denominator that are subleading in $\lambda$ (so that the result is
homogeneous in $\lambda$). Following this approximation, every master integral
has a sufficiently simple form to be solved to all orders in $\epsilon$ by the
method of Feynman parameters. Then one adds together the contributions from
the leading regions to obtain the leading behaviour in the limit $x_{3}\to 0$.
For any master integral, the region $R_{1}$ gives a non-integer power of
$x_{3}$ for $\epsilon\neq 0$, namely $x_{3}^{1-\epsilon}$. This is because all
denominators behave like $\lambda^{0}$, whilst the phase space contributes
$\lambda^{1-\epsilon}$. For $I_{1}(1,0,1,1)$ the region $R_{2}$ gives an
integer power of $x_{3}$, namely $x_{3}^{1}$, because $D_{1}$ behaves like
$\lambda^{-1}$, whilst all other denominators and the phase space behave like
$\lambda^{0}$.

For each master integral, we checked our analytic results against a numerical
computation of the integral performed by FIESTA 2 [87] at 10 randomly chosen
$(x_{1},x_{2})$ points. The results were found to agree within the precision
of the numerical computation.

When computing the expressions for the graphs, careful consideration of the
limit $x_{3}\to 0$ is needed, because some graphs are singular in this limit.
Fortunately, we know the full (all-order in $\epsilon$) behaviour of the
graphs for $x_{3}\to 0$ from the boundary condition computation for the master
integrals. For graphs associated with the non-UND topologies, this singularity
is always regulated by the dimensional regularisation parameter $\epsilon$,
since the non-UND master integrals only have the leading region $R_{1}$ for
$x_{3}\to 0$. We make a distributional expansion of the relevant
$x_{3}^{-1-n\epsilon}$ factors as follows:

$\theta(x_{3})\,x_{3}^{-1-n\epsilon}=-\dfrac{1}{n\epsilon}\,\delta(x_{3})+\mathcal{L}_{0}(x_{3})-n\epsilon\mskip
1.5mu\mathcal{L}_{1}(x_{3})+\dfrac{n^{2}\epsilon^{2}}{2}\,\mathcal{L}_{2}(x_{3})+\mathcal{O}(\epsilon^{3})\,,$
(108)

where we denote the plus distributions as

$\mathcal{L}_{n}(x)=\left[\dfrac{\theta(x)\,\ln^{n}(x)}{x}\right]_{+}\,.$
(109)

Useful relations for these distributions are given in [88]. Our final results
turn out to involve only $\mathcal{L}_{0}$ and not $\mathcal{L}_{n}$ with
$n>0$.

For graphs associated with the UND topologies, specifically in the
contribution from $I_{1}(1,0,1,1)$, one can have terms which just go as
$1/x_{3}$ and are not regulated by $\epsilon$. This corresponds to a rapidity
divergence and requires an additional regulator. In covariant gauge, the
rapidity regulator is typically inserted via some modification of the Wilson
lines, whilst in light-cone gauge it corresponds to a modification of the
gluon propagator. As noted in section 3.5 of [11] or in [64], DPDs in the
colour singlet channel do not suffer from rapidity divergences.555This is
analogous to what happens in the single scattering sector, where for the
collinear PDFs we have no rapidity divergences [89]. Therefore, the
divergences associated with $I_{1}(1,0,1,1)$ must cancel when we sum over
graphs, including the virtual loop graphs.

There are various possible choices for the rapidity regulator, such as the
analytic regulator of [90], the exponential regulator of [91], the “pure
rapidity regulator” of [92], the $\delta$ regulator of [93], the CMU $\eta$
regulator of [94, 95], or the use of Wilson lines tilted away from the light
cone as described in [70]. We tried each of the latter three options, and we
find that when consistently implementing the regulator in all real emission
and virtual loop graphs, the sum over graphs yields the same rapidity-
divergence-free result for all regulators. Note that we always remove the
rapidity regulator before performing ultraviolet renormalisation and setting
$\epsilon\to 0$, as is prescribed when using such regulators. Note also that
for each regulator we have a distributional expansion analogous to (108). We
will give more detail in a forthcoming publication devoted to two-loop results
for the colour non-singlet channels. In this case, the rapidity divergences
only cancel when one combines the bare DPDs with the appropriate soft factor.

####  4.2.2 Virtual loop diagrams

Just as for the real graphs, the integrations over minus momenta (including
$\ell^{-}$) in virtual graphs are easily done using complex contour
integration. The integrations over the transverse momenta $\boldsymbol{k}_{1}$
and $\boldsymbol{\ell}$ were performed using different methods in the Feynman
gauge and light-cone gauge calculations. In the former case we again used
integration by parts reduction, as we did for the real emission graphs. One
can even re-use the same master integrals, namely $I_{1}(0,1,1,0)$ and
$I_{1}(0,1,1,1)$, where the latter is obtained from $I_{1}(1,0,1,1)$ by
interchanging $x_{1}\leftrightarrow x_{2}$. In the calculation using light-
cone gauge, the adopted strategy closely followed that of [52] – indeed we
were able to directly re-use many of the results given in Appendix A of that
work.

Finally, we have to integrate $\ell^{+}$ over a finite range. This finite
range is a result of the earlier integrations over minus components: when
$\ell^{+}$ is outside a certain range, all poles move into one half-plane for
one of these integrations, and the result is zero. The $\ell^{+}$ integration
may be performed using standard integration techniques. However, this again
requires a consistent treatment of endpoint singularities. We have to perform
a distributional expansion analogous to (108) in cases where these
singularities are regulated by $\epsilon$, and when this is not the case we
have to insert a rapidity regulator. In light-cone gauge, the graphs
containing such a rapidity divergence are vertex correction graphs in which
the line with momentum $\ell$ is a gluon, such as figure 3(f). In Feynman
gauge the rapidity divergence is contained in the corresponding Wilson line
graphs. Here we insert the same regulator as used in the real emission
calculations.

##  5 Results

In this section, we present our results for the kernels $P_{s}$ and $W$ up to
NLO and discuss some of their properties. Results for the kernels $V$ and $U$
can readily be obtained from the relations (2.3.2) and (79).

##### Support and singularity structure

Let us briefly discuss the singularity structure of the NLO results. This
differs between the LO and NLO channels defined at the beginning of section
4.1. The kernels for LO channels have the general form

$\displaystyle K(x_{1},x_{2})$
$\displaystyle=K_{\text{reg}}(x_{1},x_{2})+\frac{K_{p}(x_{1},x_{2})}{[\mskip
1.5mu1-x_{1}-x_{2}\mskip
1.5mu]_{+}}+K_{\delta}(x_{1})\;\delta(1-x_{1}-x_{2})\,,$ (110)

where $K$ is either $W^{(2)}$ or $P_{s}^{(1)}$. The kernels for NLO channels
contain only a term $K_{\text{reg}}$. The functions $K_{\text{reg}}$ and
$K_{p}$ are regular in the region defined by the conditions $0<x_{1}$,
$0<x_{2}$ and $x_{1}+x_{2}\leq 1$, and $K_{\delta}$ is regular in the region
$0<x_{1}<1$. Note that these regions exclude the points $x_{1}=1$ and
$x_{2}=1$, where some kernels have power-law divergences of the form
$(1-x_{1})^{-n_{1}}(1-x_{2})^{-n_{2}}$ with integers $n_{1}$ and $n_{2}$.
These points are not reached in any physical double parton scattering process,
where both $x_{1}$ and $x_{2}$ must be strictly positive. We note that the
decomposition (110) is invariant under the simultaneous replacement

$\displaystyle K_{\text{reg}}(x_{1},x_{2})$ $\displaystyle\to
K_{\text{reg}}(x_{1},x_{2})+\varphi(x_{1},x_{2})\,,$ $\displaystyle
K_{p}(x_{1},x_{2})$ $\displaystyle\to
K_{p}(x_{1},x_{2})-(1-x_{1}-x_{2})\,\varphi(x_{1},x_{2})\,,$ (111)

where $\varphi(x_{1},x_{2})$ is a regular function. This leaves a freedom of
choice in the form of $K_{\text{reg}}$ and $K_{p}$ for a given kernel.

Inserting the form (110) into the convolution with a PDF and using the
representation (13), one obtains one-variable distributions $1/[1-z]_{+}$ and
$\delta(1-z)$, so that

$\displaystyle\Bigl{[}\mskip
1.5muK\underset{12}{\otimes}f\,\Bigr{]}(x_{1},x_{2})$
$\displaystyle=\frac{1}{x_{1}+x_{2}}\int\limits_{x_{1}+x_{2}}^{1}\\!\\!dz\,\Biggl{\\{}\mskip
1.5muK_{\text{reg}}(zu,z\bar{u})+\frac{K_{p}\bigl{(}zu,z\bar{u}\bigr{)}}{[\mskip
1.5mu1-z\mskip 1.5mu]_{+}}\mskip
1.5mu\Biggr{\\}}\;f\biggl{(}\frac{x_{1}+x_{2}}{z}\biggr{)}$
$\displaystyle\quad+K_{\delta}(u)\;\frac{f(x_{1}+x_{2})}{x_{1}+x_{2}}$ (112)

with

$\displaystyle u$ $\displaystyle=\frac{x_{1}}{x_{1}+x_{2}}\,,$
$\displaystyle\bar{u}$ $\displaystyle=1-u\,.$ (113)

The term with $K_{\delta}$ in (5) has the same structure as the convolution of
the LO kernels $W^{(1)}$ or $P_{s}^{(0)}$ with $f$.

Before presenting the kernels in detail, we note that

$\displaystyle W_{q\bar{q},g}(x_{1},x_{2})$
$\displaystyle=W_{\bar{q}q,g}(x_{1},x_{2})\,,$ $\displaystyle
W_{gq,g}(x_{1},x_{2})$ $\displaystyle=W_{g\bar{q},g}(x_{1},x_{2})\,,$
$\displaystyle W_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}(x_{1},x_{2})$
$\displaystyle=W_{\bar{q}^{\prime}\\!q^{\prime}\\!,q}(x_{1},x_{2})\,,$
$\displaystyle W_{qq^{\prime}\\!,q}(x_{1},x_{2})$
$\displaystyle=W_{q\bar{q}^{\prime}\\!,q}(x_{1},x_{2})\,.$ (114)

Analogous relations hold for the kernels $P_{a_{1}a_{2},a_{0}}$. The relations
in the first line follows from charge conjugation, whereas those in the second
line are obtained by reversing the quark line associated with $q^{\prime}$
(but not the one associated with $q$). The channels on the r.h.s. of (5) are
not discussed further in the following. We also recall that kernels for
specific flavour transitions are obtained from (4.1) and its analogue for
$P_{a_{1}a_{2},a_{0}}$, and that kernels for channels with an initial
$\bar{q}$ are equal to the kernels for the charge conjugate channels with an
initial $q$.

Useful building blocks for presenting the explicit kernels are the functions

$\displaystyle p_{gg}(x)$
$\displaystyle=\frac{x}{1-x}+\frac{1-x}{x}+x(1-x)\,,$ $\displaystyle
p_{qq}(x)$ $\displaystyle=\frac{1+x^{2}}{1-x}\,,$ $\displaystyle p_{qg}(x)$
$\displaystyle=x^{2}+(1-x)^{2}\,,$ $\displaystyle p_{gq}(x)$
$\displaystyle=\frac{1+(1-x)^{2}}{x}\,,$ (115)

which are proportional to the LO DGLAP evolution kernels away from the
endpoint $x=1$. The LO kernels $P_{s}^{(0)}$ and $W^{[1,0]}$ read

$\displaystyle P^{(0)}_{gg,g}(x)$ $\displaystyle=2C_{A}\,p_{gg}(x)$
$\displaystyle W^{[1,0]}_{gg,g}(x)$ $\displaystyle=0\,,$ $\displaystyle
P^{(0)}_{q\bar{q},g}(x)$ $\displaystyle=T_{F}\,p_{qg}(x)\,,$ $\displaystyle
W^{[1,0]}_{q\bar{q},g}(x)$ $\displaystyle=-2\mskip 1.5muT_{F}\mskip
1.5mux(1-x)\,$ $\displaystyle P^{(0)}_{qg,q}(x)$
$\displaystyle=C_{F}\,p_{qq}(x)\,,$ $\displaystyle W^{[1,0]}_{qg,q}(x)$
$\displaystyle=-C_{F}\,(1-x)\,.$ (116)

We find it useful to express the NLO kernels in terms of the overcomplete set
of momentum fractions $x_{1},x_{2}$ and $x_{3}=1-x_{1}-x_{2}$. This makes
symmetry properties more transparent and often allows for more compact
expressions. Furthermore we will use the notation $\bar{x}_{i}=1-x_{i}$. Let
us first discuss the evolution kernels $P_{s}^{(1)}$.

#####  $1\to 2$ evolution kernels.

Starting with the pure gluon channel, we introduce auxiliary functions

$\displaystyle
R_{gg,g}(x_{1},x_{2},x_{3})=C_{A}^{2}\;\Biggl{\\{}\frac{\log(1+x_{1}/x_{3})+\log(1+x_{1}/x_{2})}{x_{1}}\;\frac{x_{2}\,p_{gg}(x_{2})}{\bar{x}_{3}}$
$\displaystyle\qquad+2\,\biggl{[}\frac{x_{2}^{2}\,(3+4x_{2})}{\bar{x}_{1}^{4}}-\frac{x_{2}\,(4+9x_{2}+x_{2}^{2})}{\bar{x}_{1}^{3}}+\frac{1+6x_{2}+2x_{2}^{2}}{\bar{x}_{1}^{2}}-\frac{1+x_{2}}{\bar{x}_{1}}+\frac{1}{2}\biggr{]}$
$\displaystyle\qquad+\log
x_{1}\,\biggl{[}\frac{4x_{1}^{2}}{\bar{x}_{3}^{4}}-\frac{4x_{1}(1+x_{1})}{\bar{x}_{3}^{3}}+\frac{4(2+x_{1}+2x_{1}^{2})}{\bar{x}_{3}^{2}}-\frac{3-3x_{1}+x_{1}^{2}+x_{1}^{3}}{\bar{x}_{1}^{4}}\,x_{2}^{2}$
$\displaystyle\qquad\qquad+\frac{10-25x_{1}+12x_{1}^{2}-3x_{1}^{3}}{2\bar{x}_{1}^{2}}-\frac{2+7x_{1}-3x_{1}^{2}}{\bar{x}_{2}}-\frac{2}{\bar{x}_{1}\bar{x}_{2}\bar{x}_{3}}\mskip
1.5mu\biggr{]}$
$\displaystyle\qquad-\frac{\log\bar{x}_{1}}{x_{1}}\,\biggl{[}\frac{2+3x_{1}-9x_{1}^{2}+9x_{1}^{3}-x_{1}^{4}}{\bar{x}_{1}^{4}}\,x_{2}^{2}+\frac{6+4x_{1}-5x_{1}^{2}+10x_{1}^{3}-3x_{1}^{4}}{2\bar{x}_{1}^{2}}$
$\displaystyle\qquad\qquad-\frac{4-6x_{1}+9x_{1}^{2}-6x_{1}^{3}+3x_{1}^{4}}{\bar{x}_{1}\bar{x}_{2}}\mskip
1.5mu\biggr{]}\Biggr{\\}}\,,$ $\displaystyle
D_{gg,g}(x_{1})={}-C_{A}^{2}\;\Biggl{\\{}\frac{3}{4}+2\mskip
1.5mup_{gg}(x_{1})\,\biggl{[}\log
x_{1}\,\log\bar{x}_{1}+\frac{\pi^{2}}{6}-\frac{2}{3}\biggr{]}\Biggr{\\}}$
$\displaystyle\qquad+\beta_{0}\,C_{A}\,\frac{40-27x_{1}-20x_{1}^{3}}{12x_{1}}\,,$
(117)

in terms of which the $1\to 2$ evolution kernel can be written as

$\displaystyle
P^{(1)}_{gg,g}(x_{1},x_{2},x_{3})=R(x_{1},x_{2},x_{3})+R(x_{2},x_{1},x_{3})+R(x_{3},x_{2},x_{1})$
$\displaystyle\quad+R(x_{1},x_{3},x_{2})+R(x_{3},x_{1},x_{2})+R(x_{2},x_{3},x_{1})+\bigl{[}D(x_{1})+D(x_{2})\bigr{]}\,\delta(x_{3})\,.$
(118)

Here we dropped the subscripts on $R$ and $D$ for brevity. If we replace
$x_{3}=1-x_{1}-x_{2}$ and express $R$ as a function of only $x_{1}$ and
$x_{2}$, we have

$\displaystyle
P^{(1)}_{gg,g}(x_{1},x_{2})=R(x_{1},x_{2})+R(x_{2},x_{1})+R(1-x_{1}-x_{2},x_{2})+R(x_{1},1-x_{1}-x_{2})$
$\displaystyle\quad+R(1-x_{1}-x_{2},x_{1})+R(x_{2},1-x_{1}-x_{2})+\bigl{[}D(x_{1})+D(x_{2})\bigr{]}\,\delta(1-x_{1}-x_{2})\,.$
(119)

We see in (5) that, apart from the distribution term involving
$\delta(x_{3})$, the kernel is fully symmetric in the momentum fractions of
the three final-state gluons. The symmetry under the interchange
$x_{1}\leftrightarrow x_{2}$ is trivial, but the symmetry between an observed
gluon ($a_{1}$ or $a_{2}$) and the unobserved one ($a_{3}$) is not. We discuss
this “active-spectator symmetry” (where we refer to $a_{1},a_{2}$ as active
partons and to $a_{3}$ as the spectator) in section 5.2.

For the $g\to q\bar{q}$ channel, we define

$\displaystyle R_{q\bar{q},g}(x_{1},x_{2},x_{3})=C_{A}\mskip
1.5muT_{F}\,\Biggl{\\{}\frac{\log(1+x_{3}/x_{2})}{x_{3}}\;\frac{p_{qg}(x_{1})+p_{qg}(x_{2})}{2}$
$\displaystyle\qquad-
x_{1}x_{2}\,\frac{3-10x_{3}+3x_{3}^{2}}{\bar{x}_{3}^{4}}+\frac{1-6x_{3}+5x_{3}^{2}-2x_{3}^{3}}{2\bar{x}_{3}^{2}}$
$\displaystyle\qquad-2\mskip 1.5mu\log
x_{1}\,\biggl{[}\frac{2x_{1}^{2}}{\bar{x}_{3}^{4}}-\frac{2x_{1}(1+x_{1})}{\bar{x}_{3}^{3}}+\frac{1+2x_{1}+4x_{1}^{2}}{\bar{x}_{3}^{2}}-\frac{1+6x_{1}-2x_{1}^{2}}{2\bar{x}_{3}}+1-x_{1}-\frac{x_{3}}{2}\biggr{]}$
$\displaystyle\qquad+2\log\bar{x}_{1}\,\biggl{[}\frac{1-2x_{1}+2x_{1}^{2}}{2\bar{x}_{3}}+1-x_{1}-\frac{x_{3}}{2}\biggr{]}$
$\displaystyle\qquad-2\mskip 1.5mu\log
x_{3}\,\biggl{[}x_{1}x_{2}\,\frac{2-3x_{3}+2x_{3}^{2}}{\bar{x}_{3}^{4}}-\frac{1-x_{3}+x_{3}^{2}}{2\bar{x}_{3}^{2}}\biggr{]}$
$\displaystyle\qquad-\frac{2\mskip
1.5mu\log\bar{x}_{3}}{x_{3}}\,\biggl{[}x_{1}x_{2}\,\frac{1+x_{3}-3x_{3}^{2}+3x_{3}^{3}}{\bar{x}_{3}^{4}}-\frac{1+x_{3}-2x_{3}^{2}+2x_{3}^{3}}{2\bar{x}_{3}^{2}}\biggr{]}\Biggr{\\}}$
$\displaystyle\quad+C_{F}\mskip
1.5muT_{F}\,\Biggl{\\{}-\biggl{[}\frac{3x_{3}}{\bar{x}_{1}^{2}}-\frac{1+5x_{3}}{\bar{x}_{1}}+3+2x_{3}\biggr{]}+2\mskip
1.5mu\log
x_{1}\,\biggl{[}\frac{x_{3}}{2\bar{x}_{1}^{2}}-\frac{1+x_{3}}{\bar{x}_{1}}-\frac{x_{3}}{2\bar{x}_{2}^{2}}+\frac{1+x_{3}}{\bar{x}_{2}}\biggr{]}$
$\displaystyle\qquad+2\mskip 1.5mu(2\mskip 1.5mu\log\bar{x}_{1}-\log
x_{3})\,\biggl{[}\frac{x_{3}}{2\bar{x}_{1}^{2}}-\frac{1+x_{3}}{\bar{x}_{1}}+1\biggr{]}\Biggr{\\}}\,,$
$\displaystyle S_{q\bar{q},g}(x_{1})={}-2\mskip 1.5mu(C_{A}-2C_{F})\mskip
1.5muT_{F}\,x_{1}\bar{x}_{1}\,,$ $\displaystyle
D_{q\bar{q},g}(x_{1})=C_{A}\mskip
1.5muT_{F}\,\Bigl{\\{}x_{1}^{2}-\log^{2}\\!x_{1}\;p_{qq}(x_{1})\Bigr{\\}}+C_{F}\mskip
1.5muT_{F}\,\Biggl{\\{}\frac{3}{2}-2x_{1}\bar{x}_{1}\bigl{(}\log
x_{1}+\log\bar{x}_{1}\bigr{)}$
$\displaystyle\qquad+\biggl{[}\log^{2}\\!x_{1}-\log
x_{1}\,\log\bar{x}_{1}-\frac{\pi^{2}}{6}+\frac{3}{2}\biggr{]}\,p_{qg}(x_{1})\Biggr{\\}}$
(120)

and have an evolution kernel

$\displaystyle P^{(1)}_{q\bar{q},g}(x_{1},x_{2},x_{3})$
$\displaystyle=R_{q\bar{q},g}(x_{1},x_{2},x_{3})+R_{q\bar{q},g}(x_{2},x_{1},x_{3})+\frac{S_{q\bar{q},g}(x_{1})+S_{q\bar{q},g}(x_{2})}{[\mskip
1.5mux_{3}\mskip 1.5mu]_{+}}$
$\displaystyle\quad+\bigl{[}D_{q\bar{q},g}(x_{1})+D_{q\bar{q},g}(x_{2})\bigr{]}\,\delta(x_{3})\,.$
(121)

The evolution kernel for the channel $g\to qg$ can be expressed by the same
functions as

$\displaystyle P^{(1)}_{qg,g}(x_{1},x_{2},x_{3})$
$\displaystyle=R_{q\bar{q},g}(x_{1},x_{3},x_{2})+R_{q\bar{q},g}(x_{3},x_{1},x_{2})+\frac{S_{q\bar{q},g}(x_{1})+S_{q\bar{q},g}(x_{3})}{x_{2}}\,.$
(122)

As we did in the pure gluon channel, we observe active-spectator symmetry for
the splitting process $g\to q\bar{q}g$. We obtain the kernel
$P^{(1)}_{qg,g}(x_{1},x_{2},x_{3})$ from
$P^{(1)}_{q\bar{q},g}(x_{1},x_{2},x_{3})$ if we drop the distributions terms,
omitting $\delta(x_{3})$ and replacing $1/[x_{3}]_{+}$ with $1/x_{3}$, and
then interchange $x_{2}\leftrightarrow x_{3}$.

Turning to the channel $q\to qg$, we define

$\displaystyle R_{qg,q}(x_{1},x_{2},x_{3})=C_{A}\mskip
1.5muC_{F}\,\Biggl{\\{}-\frac{\log(1+x_{3}/x_{1})+\log(1-x_{3}/\bar{x}_{1})}{x_{3}}\;\frac{x_{2}\,p_{gq}(x_{2})}{\bar{x}_{1}}$
$\displaystyle\qquad+2\,\biggl{[}\frac{3x_{1}}{\bar{x}_{1}^{4}}\,x_{3}^{2}-\frac{x_{1}}{\bar{x}_{1}^{2}}-\frac{3}{4}\biggr{]}-\bigl{(}\log
x_{1}+2\mskip 1.5mu\log\bar{x}_{1}\bigr{)}\,\biggl{[}x_{1}\mskip
1.5mu\frac{2x_{2}^{2}}{\bar{x}_{1}^{4}}+\frac{3+x_{2}^{2}}{\bar{x}_{1}^{2}}-\frac{1}{\bar{x}_{1}}\biggr{]}$
$\displaystyle\qquad+2\,\log
x_{2}\,\biggl{[}\frac{2x_{2}^{2}}{\bar{x}_{1}^{4}}-\frac{2x_{2}(1+x_{2})}{\bar{x}_{1}^{3}}+\frac{4+2x_{2}+x_{2}^{2}}{\bar{x}_{1}^{2}}-\frac{2+x_{2}}{\bar{x}_{1}}+\frac{1}{2}\biggr{]}\Biggr{\\}}$
$\displaystyle\quad+C_{F}^{2}\,\Biggl{\\{}\frac{2\mskip
1.5mu\log(1+x_{3}/x_{1})}{x_{3}}\;\frac{x_{2}\,p_{gq}(x_{2})}{\bar{x}_{1}}-2\,\biggl{[}\frac{x_{3}}{\bar{x}_{2}^{2}}+\frac{1+x_{3}}{2\bar{x}_{2}}-1\biggr{]}$
$\displaystyle\qquad-\log
x_{1}\,\biggl{[}\frac{x_{2}}{\bar{x}_{3}^{2}}-\frac{2-x_{2}}{\bar{x}_{3}}-\frac{4}{\bar{x}_{1}}+1\biggr{]}-\log
x_{2}\,\biggl{[}\frac{x_{2}}{\bar{x}_{3}^{2}}-\frac{2-x_{2}}{\bar{x}_{3}}-\frac{2-x_{2}}{\bar{x}_{2}^{2}}\,x_{3}+\frac{2}{\bar{x}_{2}}\biggr{]}$
$\displaystyle\qquad+\frac{2\mskip
1.5mu\log\bar{x}_{2}}{x_{2}}\,\biggl{[}\frac{x_{3}}{\bar{x}_{2}^{2}}-\frac{2}{\bar{x}_{2}}+\frac{x_{2}\,p_{gq}(x_{2})}{\bar{x}_{1}}\biggr{]}\Biggr{\\}}\,,$
$\displaystyle S_{qg,q}(x_{1})=2\mskip 1.5muC_{A}\mskip
1.5muC_{F}\,\bar{x}_{1}\,,$ $\displaystyle D_{qg,q}(x_{1})=C_{A}\mskip
1.5muC_{F}\,\Biggl{\\{}1-2\bar{x}_{1}\log\bar{x}_{1}+\biggl{[}\log^{2}\\!x_{1}-2\log
x_{1}\mskip
1.5mu\log\bar{x}_{1}-2\operatorname{Li}_{2}(x_{1})+\frac{4}{3}\biggr{]}\,p_{qq}(x_{1})\Biggr{\\}}$
$\displaystyle\qquad-C_{F}^{2}\,\Biggl{\\{}1+2\bar{x}_{1}\log
x_{1}+\biggl{[}\log^{2}\\!x_{1}+3\log
x_{1}-2\operatorname{Li}_{2}(x_{1})+\frac{\pi^{2}}{3}\biggr{]}\,p_{qq}(x_{1})\Biggr{\\}}$
$\displaystyle\qquad+\beta_{0}\mskip
1.5muC_{F}\,\Biggl{\\{}\bar{x}_{1}+\biggl{[}\log
x_{1}+\frac{5}{3}\biggr{]}\,p_{qq}(x_{1})\Biggr{\\}}$ (123)

and have an evolution kernel

$\displaystyle P^{(1)}_{qg,q}(x_{1},x_{2},x_{3})$
$\displaystyle=R_{qg,q}(x_{1},x_{2},x_{3})+R_{qg,q}(x_{1},x_{3},x_{2})+\frac{S_{qg,q}(x_{1})}{[\mskip
1.5mux_{3}\mskip 1.5mu]_{+}}+\frac{S_{qg,q}(x_{1})}{x_{2}}$
$\displaystyle\quad+D_{qg,q}(x_{1})\,\delta(x_{3})\,.$ (124)

In terms of the same functions, the evolution kernel for $q\to gg$ reads

$\displaystyle P^{(1)}_{gg,q}(x_{1},x_{2},x_{3})$
$\displaystyle=R_{qg,q}(x_{3},x_{2},x_{1})+R_{qg,q}(x_{3},x_{1},x_{2})+\frac{S_{qg,q}(x_{3})}{x_{1}}+\frac{S_{qg,q}(x_{3})}{x_{2}}\,.$
(125)

This is symmetric in the arguments $x_{1}$ and $x_{2}$, as it must be. Again
we observe active-spectator symmetry, this time under the exchange
$x_{1}\leftrightarrow x_{3}$.

For quark-antiquark transitions, we have

$\displaystyle
P^{(1)}_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}(x_{1},x_{2},x_{3})$
$\displaystyle=R_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}(x_{1},x_{2},x_{3})+R_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}(x_{2},x_{1},x_{3})$
(126)

with

$\displaystyle R_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}={}-C_{F}\mskip
1.5muT_{F}\,$
$\displaystyle\biggl{\\{}x_{1}x_{2}\,\frac{1-6x_{3}+x_{3}^{2}}{\bar{x}_{3}^{4}}+\frac{x_{3}}{\bar{x}_{3}^{2}}$
$\displaystyle+\bigl{(}2\mskip 1.5mu\log x_{1}-\log x_{3}-2\mskip
1.5mu\log\bar{x}_{3}\bigr{)}\,\frac{x_{1}^{2}+x_{2}^{2}}{2\bar{x}_{3}^{3}}\,p_{qq}(x_{3})\biggr{\\}}\,,$
(127)

and for the quark valence kernels we have

$\displaystyle P^{v\,(1)}_{qq,q}(x_{1},x_{2},x_{3})$
$\displaystyle=R_{qq,q}^{v}(x_{1},x_{2},x_{3})+R_{qq,q}^{v}(x_{2},x_{1},x_{3})$
(128)

with

$\displaystyle R_{qq,q}^{v}(x_{1},x_{2},x_{3})=-C_{F}\,(C_{A}-2C_{F})\,$
$\displaystyle\biggl{\\{}\frac{2x_{3}}{\bar{x}_{1}^{2}}-\frac{1+x_{3}}{\bar{x}_{1}}+\bigl{(}2\mskip
1.5mu\log\bar{x}_{1}-\log
x_{3}\bigr{)}\,\frac{1+x_{3}^{2}}{2\bar{x}_{1}\bar{x}_{2}}\biggr{\\}}\,.$
(129)

The evolution kernels for the remaining channels again obey active-spectator
symmetry and are given by

$\displaystyle P^{(1)}_{qq^{\prime}\\!,q}(x_{1},x_{2},x_{3})$
$\displaystyle=P^{(1)}_{q^{\prime}\\!q,q}(x_{2},x_{1},x_{3})=P^{(1)}_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}(x_{2},x_{3},x_{1})$
(130) and $\displaystyle P^{v\,(1)}_{q\bar{q},q}(x_{1},x_{2},x_{3})$
$\displaystyle=P^{v\,(1)}_{qq,q}(x_{1},x_{3},x_{2})\,.$ (131)

##### Terms originating from LO kernels

As specified in (2.3.1), the kernels $W^{[2,1]}$ and $W^{[2,2]}$ can be
constructed from $P_{s}^{(1)}$ and from additional terms that originate from
LO kernels. For LO channels, we find that these additional terms have the form

$\displaystyle W_{a_{1}a_{2},a_{0}}$ $\displaystyle:\leavevmode\nobreak\
\frac{\beta_{0}}{2}\,K_{a_{1}a_{2},a_{0}}+P^{(0)}_{a_{1}a_{1}}\underset{1}{\otimes}K_{a_{1}a_{2},a_{0}}+P^{(0)}_{a_{2}a_{2}}\underset{2}{\otimes}K_{a_{1}a_{2},a_{0}}-K_{a_{1}a_{2},a_{0}}\underset{12}{\otimes}P^{(0)}_{a_{0}a_{0}}\,,$
(132)

where $K=W^{[1,0]}$ or $P_{s}^{(0)}$. Note that repeated parton indices on the
r.h.s. are _not_ summed over. The flavour diagonal DGLAP kernels
$P^{(0)}_{a_{i}a_{i}}$ in (132) contain distribution terms, which we wish to
make explicit. To this end, we write the kernels in the form

$\displaystyle P_{gg}^{(0)}(x)$
$\displaystyle=P_{gg,\text{reg}}^{(0)}(x)+\frac{2C_{A}}{[\mskip 1.5mu1-x\mskip
1.5mu]_{+}}+\frac{\beta_{0}}{2}\,\delta(1-x)\,,$ $\displaystyle
P_{qq}^{(0)}(x)$
$\displaystyle=P_{qq,\text{reg}}^{(0)}(x)+\frac{2C_{F}}{[\mskip 1.5mu1-x\mskip
1.5mu]_{+}}+\frac{3}{2}\,C_{F}\,\delta(1-x)$ (133)

in analogy to (110), where the regular terms read

$\displaystyle P_{gg,\text{reg}}^{(0)}(x)$
$\displaystyle=2C_{A}\biggl{[}\frac{1-x}{x}+x(1-x)-1\biggr{]}\,,$
$\displaystyle P_{qq,\text{reg}}^{(0)}(x)$ $\displaystyle=-C_{F}\mskip
1.5mu(1+x)\,.$ (134)

We then have convolutions

$\displaystyle P_{gg}^{(0)}\underset{1}{\otimes}K$
$\displaystyle=P_{gg,\text{reg}}^{(0)}\biggl{(}\frac{x_{1}}{\bar{x}_{2}}\biggr{)}\,\frac{K(\bar{x}_{2})}{\bar{x}_{2}}+\frac{2C_{A}}{[\mskip
1.5mux_{3}\mskip
1.5mu]_{+}}\,K(\bar{x}_{2})+\frac{1}{2}\,\bigl{(}\beta_{0}-4C_{A}\mskip
1.5mu\log x_{1}\bigr{)}\,K(x_{1})\,\delta(x_{3})\,,$ $\displaystyle
P_{gg}^{(0)}\underset{2}{\otimes}K$
$\displaystyle=P_{gg,\text{reg}}^{(0)}\biggl{(}\frac{x_{2}}{\bar{x}_{1}}\biggr{)}\,\frac{K(x_{1})}{\bar{x}_{1}}+\frac{2C_{A}}{[\mskip
1.5mux_{3}\mskip
1.5mu]_{+}}\,K(x_{1})+\frac{1}{2}\,\bigl{(}\beta_{0}-4C_{A}\mskip
1.5mu\log\bar{x}_{1}\bigr{)}\,K(x_{1})\,\delta(x_{3})\,,$ $\displaystyle
K\underset{12}{\otimes}P_{gg}^{(0)}$
$\displaystyle=K\biggl{(}\frac{x_{1}}{\bar{x}_{3}}\biggr{)}\,\frac{P_{gg,\text{reg}}^{(0)}(\bar{x}_{3})}{\bar{x}_{3}}+\frac{2C_{A}}{[\mskip
1.5mux_{3}\mskip
1.5mu]_{+}}\;\frac{1}{\bar{x}_{3}}\,K\biggl{(}\frac{x_{1}}{\bar{x}_{3}}\biggr{)}+\frac{\beta_{0}}{2}\,K(x_{1})\,\delta(x_{3})\,,$
(135) and $\displaystyle P_{qq}^{(0)}\underset{1}{\otimes}K$
$\displaystyle=P_{qq,\text{reg}}^{(0)}\biggl{(}\frac{x_{1}}{\bar{x}_{2}}\biggr{)}\,\frac{K(\bar{x}_{2})}{\bar{x}_{2}}+\frac{2C_{F}}{[\mskip
1.5mux_{3}\mskip
1.5mu]_{+}}\,K(\bar{x}_{2})+\frac{1}{2}\,C_{F}\,\bigl{(}3-4\log
x_{1}\bigr{)}\,K(x_{1})\,\delta(x_{3})\,,$ $\displaystyle
P_{qq}^{(0)}\underset{2}{\otimes}K$
$\displaystyle=P_{qq,\text{reg}}^{(0)}\biggl{(}\frac{x_{2}}{\bar{x}_{1}}\biggr{)}\,\frac{K(x_{1})}{\bar{x}_{1}}+\frac{2C_{F}}{[\mskip
1.5mux_{3}\mskip
1.5mu]_{+}}\,K(x_{1})+\frac{1}{2}\,C_{F}\,\bigl{(}3-4\log\bar{x}_{1}\bigr{)}\,K(x_{1})\,\delta(x_{3})\,,$
$\displaystyle K\underset{12}{\otimes}P_{qq}^{(0)}$
$\displaystyle=K\biggl{(}\frac{x_{1}}{\bar{x}_{3}}\biggr{)}\,\frac{P_{qq,\text{reg}}^{(0)}(\bar{x}_{3})}{\bar{x}_{3}}+\frac{2C_{F}}{[\mskip
1.5mux_{3}\mskip
1.5mu]_{+}}\;\frac{1}{\bar{x}_{3}}\,K\biggl{(}\frac{x_{1}}{\bar{x}_{3}}\biggr{)}+\frac{3}{2}\,C_{F}\,K(x_{1})\,\delta(x_{3})\,,$
(136)

where we have omitted the parton labels on $K$ for brevity. To obtain these
relations, we used equation (B13) in [88] to bring all plus distribution terms
into the form $1/[x_{3}]_{+}$.

The convolution terms that appear in the different NLO channels are

$\displaystyle W_{qg,g}$ $\displaystyle:\leavevmode\nobreak\
P^{(0)}_{qg}\underset{1}{\otimes}K_{gg,g}+P^{(0)}_{gq}\underset{2}{\otimes}K_{q\bar{q},g}-K_{qg,q}\underset{12}{\otimes}P^{(0)}_{qg}\,,$
$\displaystyle W_{gg,q}$ $\displaystyle:\leavevmode\nobreak\
P^{(0)}_{gq}\underset{1}{\otimes}K_{qg,q}+P^{(0)}_{gq}\underset{2}{\otimes}K_{gq,q}-K_{gg,g}\underset{12}{\otimes}P^{(0)}_{gq}\,,$
(137)

and

$\displaystyle W_{q_{j}\bar{q}_{k},q_{i}}$ $\displaystyle:\leavevmode\nobreak\
\phantom{+P^{(0)}_{gq}\underset{1}{\otimes}K_{qg,q}\;\;}\delta_{ij}\,P^{(0)}_{qg}\underset{2}{\otimes}K_{qg,q}-\delta_{jk}\mskip
1.5muK_{q\bar{q},g}\underset{12}{\otimes}P^{(0)}_{gq}\,,$ $\displaystyle
W_{q_{j}q_{k},q_{j}}$ $\displaystyle:\leavevmode\nobreak\
\delta_{jk}\,P^{(0)}_{qg}\underset{1}{\otimes}K_{gq,q}+P^{(0)}_{qg}\underset{2}{\otimes}K_{qg,q}\,.$
(138)

In all these cases, only the flavour non-diagonal DGLAP kernels

$\displaystyle P^{(0)}_{qg}(x)$ $\displaystyle=P_{q\bar{q},g}^{(0)}(x)\,,$
$\displaystyle P^{(0)}_{gq}(x)$ $\displaystyle=P_{gq,q}^{(0)}(x)$ (139)

appear, after we used charge conjugation symmetry to replace
$P^{(0)}_{\bar{q}g}$ and $P^{(0)}_{g\bar{q}}$.

We note that for all LO and NLO channels $a_{0}\to a_{1}a_{2}$, there is
exactly one parton combination in each type of convolution term. This will no
longer hold for kernels at order $\alpha_{s}^{3}$, where for given $a_{0}$,
$a_{1}$, and $a_{2}$ there is more than one possibility for the spectator
partons.

In contrast to what we observed for the kernels $P_{s}^{(1)}$, active-
spectator symmetry does _not_ hold for the LO induced terms just discussed,
and therefore this symmetry does not hold for $W^{[2,1]}$ and $W^{[2,2]}$
either. This is most easily seen for $W_{gg,g}^{[2,2]}$. Apart from
distribution terms, all functions in the convolutions (2.3) are equal to
$2C_{A}\mskip 1.5mup_{gg}$ in that case. The exchange $x_{1}\leftrightarrow
x_{3}$ then leaves the first term in (2.3) invariant, whilst interchanging the
second and third terms. Because the latter enter $W_{gg,g}^{[2,2]}$ with
opposite signs in (2.3.1), active-spectator symmetry is broken.

##### Nonlogarithmic terms.

The full expressions of the kernels $W^{[2,0]}$, which are not accompanied by
logarithms in $W^{(2)}$, are rather lengthy, and we only discuss their general
features here. Their full expressions are given in the ancillary files
associated with this paper on https://arxiv.org. We note that the functions
$W^{[2,0]}$ are not needed for the position space matching kernels $V^{(2)}$,
but they do appear in the matching kernels $U^{(2)}$ between momentum and
position space DPDs.

The general form of the $W^{[2,0]}$ coefficients can be written as

$\displaystyle W_{gg,g}^{[2,0]}$
$\displaystyle:C_{A}^{2}\biggl{[}R_{gg,g}^{A}+\delta(x_{3})\,D_{gg,g}^{A}\biggr{]}+\beta_{0}\,C_{A}\,\delta(x_{3})\,D_{gg,g}^{\beta}\,,$
$\displaystyle W_{q\bar{q},g}^{[2,0]}$ $\displaystyle:C_{A}\mskip
1.5muT_{F}\mskip
1.5mu\Biggl{\\{}R_{q\bar{q},g}^{A}+\delta(x_{3})\,D_{q\bar{q},g}^{A}+\frac{S_{q\bar{q},g}^{A}}{[x_{3}]_{+}}\Biggr{\\}}+C_{F}\mskip
1.5muT_{F}\mskip
1.5mu\Biggl{\\{}R_{q\bar{q},g}^{F}+\delta(x_{3})\,D_{q\bar{q},g}^{F}+\frac{S_{q\bar{q},g}^{F}}{[x_{3}]_{+}}\Biggr{\\}}\,,$
$\displaystyle W_{qg,q}^{[2,0]}$ $\displaystyle:C_{A}\mskip
1.5muC_{F}\biggl{[}R_{qg,q}^{A}+\delta(x_{3})\,D_{qg,q}^{A}\biggr{]}+C_{F}^{2}\biggl{[}R_{qg,q}^{F}+\delta(x_{3})\,D_{qg,q}^{F}\biggr{]}+\beta_{0}\,C_{F}\,\delta(x_{3})\,D_{qg,q}^{\beta}\,,$
$\displaystyle W_{qg,g}^{[2,0]}$ $\displaystyle:C_{A}\mskip
1.5muT_{F}\,R_{qg,g}^{A}+C_{F}\mskip 1.5muT_{F}\,R_{qg,g}^{F}\,,$
$\displaystyle W_{gg,q}^{[2,0]}$ $\displaystyle:C_{A}\mskip
1.5muC_{F}\,R_{gg,q}^{A}+C_{F}^{2}\,R_{gg,q}^{F}\,,$ (140)

and

$\displaystyle W_{qq^{\prime}\\!,q}^{[2,0]}$ $\displaystyle:C_{F}\mskip
1.5muT_{F}\,R_{qq^{\prime}\\!,q}^{F}\,,$ $\displaystyle
W_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}^{[2,0]}$ $\displaystyle:C_{F}\mskip
1.5muT_{F}\,R_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}^{F}\,,$ $\displaystyle
W_{qq,q}^{v\,[2,0]}$ $\displaystyle:C_{F}\mskip
1.5mu(C_{A}-2C_{F})\,R_{qq,q}^{F}\,,$ $\displaystyle
W_{q\bar{q},q}^{v\,[2,0]}$ $\displaystyle:C_{F}\mskip
1.5mu(C_{A}-2C_{F})\,R_{q\bar{q},q}^{F}\,,$ (141)

where the functions $R_{a_{1}a_{2},a_{0}}^{c}$ with $c=A,F,\beta$ are regular
as specified below (110). They are independent of colour factors and of the
number $n_{F}$ of active quark flavours. We see that the colour structure of
$W^{[2,0]}_{a_{1}a_{2},a_{0}}$ is the same as for the corresponding kernel
$P_{a_{1}a_{2},a_{0}}^{(1)}$.

The regular parts $R^{c}$ are by far the most lengthy ones, containing
rational functions as well as the product of rational functions with
logarithms, with products of two logarithms, and with dilogarithms. The
logarithms have arguments $x_{1}$, $\bar{x}_{1}$, $x_{2}$, $\bar{x}_{2}$,
$x_{3}$, or $\bar{x}_{3}$. The arguments of the dilogarithms can be reduced to
the set $x_{1},\,x_{2},\,\bar{x}_{3},\,x_{1}/\bar{x}_{2},\,x_{2}/\bar{x}_{1}$,
and $-x_{2}/x_{1}$, using the relations given in section II.B of [96]. Note
that with these arguments, all logarithms and dilogarithms are real valued. In
the denominators of rational functions, the highest powers of $x_{i}$ and
$\bar{x}_{i}$ (with $i=1,2,3$) are the same as in the evolution kernels
$P_{s}^{(1)}$ given earlier.

The plus distribution parts $S^{c}$ of the $W^{[2,0]}$ coefficients are either
zero or simple polynomials in $x_{1}$ and $\bar{x}_{1}$. They determine the
leading behaviour in the limit $x_{1}+x_{2}\to 1$ and will be explicitly given
in (5.3.1) below.

The delta distribution parts $D^{c}$ are in general less lengthy than the
regular parts, but they involve rational functions and the product of rational
functions with trilogarithms, with products of logarithms and dilogarithms,
with dilogarithms, and with products of up to three logarithms. The arguments
of all logarithms and polylogarithms are only $x_{1}$ or $\bar{x}_{1}$. The
rational functions often coincide with the corresponding functions in (5), or
at least they have the same denominators as these functions.

###  5.1 Parton number and momentum sum rules

Sum rules for momentum space DPDs at $\Delta=0$ were proposed in [53]. They
are among the few theoretical constraints on DPDs that we have. A detailed
proof of these sum rules was given in [54], where it was also shown that the
renormalisation scale independence of the DPD sum rules implies corresponding
sum rules for the $1\to 2$ evolution kernels. These read

$\displaystyle\int\limits_{0}^{1-x_{1}}\\!\\!dx_{2}\;\bigl{[}P_{a_{1}q,a_{0}}(x_{1},x_{2})-P_{a_{1}\bar{q},a_{0}}(x_{1},x_{2})\bigr{]}$
$\displaystyle=\bigl{(}\delta_{a_{1}\bar{q}}-\delta_{a_{1}q}-\delta_{a_{0}\mskip
1.5mu\bar{q}}+\delta_{a_{0}\mskip 1.5muq}\bigr{)}\,P_{a_{1}a_{0}}(x_{1})$
(142) and
$\displaystyle\sum_{a_{2}}\int\limits_{0}^{1-x_{1}}\\!\\!dx_{2}\,x_{2}\;P_{a_{1}a_{2},a_{0}}(x_{1},x_{2})$
$\displaystyle=(1-x_{1})\,P_{a_{1}a_{0}}(x_{1})\,.$ (143)

Note that the distribution terms in the DGLAP evolution kernels on the r.h.s.
cancel. In the number sum rule (142), such terms only appear if $a_{0}=a_{1}$,
in which case the sum of Kronecker deltas gives zero. In the momentum sum rule
(143), distribution terms are removed by the prefactor $1-x_{1}$.

For the LO kernels, the sum rules (142) and (143) are readily verified using
the relation (37) and the list of possible transitions $a_{0}\to a_{1}a_{2}$.
At NLO, however, the sum rules provide nontrivial relations between the
evolution kernels $P_{s}^{(1)}$ and the DGLAP splitting functions, which are
well known at that order [44, 45, 46, 47, 48, 49, 50, 51, 52]. These relations
serve as a valuable cross check of our results.

##### Number sum rule.

Let us first consider the number sum rule (142) for the different parton
combinations. Starting with the ones that involve $1\to 2$ evolution kernels
in LO channels, we have

$\displaystyle\int dx_{2}\;P_{q\bar{q},g}^{(1)}$
$\displaystyle=P_{qg}^{(1)}\,,$ $\displaystyle\int dx_{2}\;P_{gq,q}^{(1)}$
$\displaystyle=P_{gq}^{(1)}\,,$ (144)

where the arguments of the kernels and integration boundaries are as in (142).
Here and in the following we omit NLO kernels that vanish identically, such as
$P_{qq,g}^{(1)}$ and $P^{(1)}_{g\bar{q},q}$. Turning to the NLO channels, we
use the notation

$\displaystyle P_{q_{i}q_{k}}$
$\displaystyle=\delta_{ik}P_{qq}^{V}+P_{qq}^{S}\,,$ $\displaystyle
P_{\bar{q}_{i}{q}_{k}}$
$\displaystyle=\delta_{ik}P_{\bar{q}{q}}^{V}+P_{\bar{q}{q}}^{S}$ (145)

for DGLAP kernels [97] and note that at NLO the relation

$\displaystyle P_{qq}^{S}$ $\displaystyle=P_{\bar{q}{q}}^{S}$ (146)

holds for the flavour singlet parts. Taking linear combinations of the number
sum rules for the different transitions $q_{j}\to q_{j}q_{k}$ and $q_{i}\to
q_{j}\bar{q}_{k}$ and using the symmetry relations (5), we obtain

$\displaystyle\int dx_{2}\;P_{qq,q}^{(1)\,v}$ $\displaystyle=\int
dx_{2}\;P_{q\bar{q},q}^{(1)\,v}\,,$ $\displaystyle\phantom{=}\;\mskip
1.5mu\int dx_{2}\;P_{\bar{q}q,q}^{(1)\,v}$
$\displaystyle=2\,P_{\bar{q}q}^{V\,(1)}\,,$ $\displaystyle\int
dx_{2}\;P_{q^{\prime}\\!q,q}^{(1)}$ $\displaystyle=\int
dx_{2}\;P_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}^{(1)}$
$\displaystyle=P_{qq}^{S\,(1)}\,.$ (147)

We checked explicitly that our results for the $1\to 2$ kernels fulfil all of
the above sum rules.

##### Momentum sum rule.

For NLO kernels, one finds the momentum sum rules

$\displaystyle\int
dx_{2}\;x_{2}\,\bigl{[}P_{gg,g}^{(1)}+2n_{f}P_{gq,g}^{(1)}\bigr{]}$
$\displaystyle=(1-x_{1})P_{gg}^{(1)}\,,$ $\displaystyle\int
dx_{2}\;x_{2}\,\bigl{[}P_{gq,q}^{(1)}+P_{gg,q}^{(1)}\bigr{]}$
$\displaystyle=(1-x_{1})P_{gq}^{(1)}\,,$ $\displaystyle\int
dx_{2}\;x_{2}\,\bigl{[}P_{q\bar{q},g}^{(1)}+P_{qg,g}^{(1)}\bigr{]}$
$\displaystyle=(1-x_{1})P_{qg}^{(1)}\,,$ $\displaystyle\int
dx_{2}\;x_{2}\,\bigl{[}P_{qg,q}^{(1)}+P_{qq,q}^{(1)\,v}+P_{q\bar{q},q}^{(1)\,v}+2n_{f}P_{qq^{\prime}\\!,q}^{(1)}\bigr{]}$
$\displaystyle=(1-x_{1})P_{qq}^{V\,(1)}\,,$ $\displaystyle\int
dx_{2}\;x_{2}\,P_{\bar{q}q,q}^{(1)\,v}$
$\displaystyle=(1-x_{1})P_{\bar{q}q}^{V\,(1)}\,,$ $\displaystyle\int
dx_{2}\;x_{2}\,\bigl{[}P_{q^{\prime}\\!q,q}^{(1)}+P_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}^{(1)}\bigr{]}$
$\displaystyle=(1-x_{1})P_{qq}^{S\,(1)}\,,$ (148)

all of which are fulfilled for the $1\to 2$ evolution kernels we have
calculated.

###  5.2 Active-spectator symmetry

In this subsection, we consider kernels away from the kinematic point
$x_{3}=0$, so that only real emission graphs contribute and distribution terms
do not contribute. We find active-spectator symmetry to hold for the evolution
kernels $P_{s}^{(1)}$ in all parton channels, but not for the any of the
kernels $W^{[2,k]}$ that make up the matching coefficient $W^{(2)}$.

The lack of this symmetry for $W^{(2)}$ is not surprising, because the
renormalisation (2.3) of $W$ clearly treats observed and spectator partons in
an asymmetric way. Furthermore, already in bare graphs the active partons are
singled out because their transverse momenta differ by
$\pm\boldsymbol{\Delta}$ in the amplitude and its conjugate. This is in fact
the only difference between active partons and spectators. To see this, we
consider the bare DPDs $F_{B}$ in light-cone perturbation theory, which can be
obtained from covariant perturbation theory by integrating over all minus
momenta (see e.g. chapter 7.2.3 of [70] and references therein). Both active
and spectator partons are then on their mass shell by construction. As shown
in section 5 of [54], for ${\Delta}={0}$ the numerator factors in light-cone
perturbation theory are also identical for active and spectator partons in
bare DPDs, provided of course that one considers unpolarised partons.

To understand why active-spectator symmetry holds for
$P^{(1)}_{a_{1}a_{2},a_{0}}(x_{1},x_{2},x_{3})$, we note that for $x_{3}>0$
these evolution kernels are associated with the ultraviolet divergences of
two-loop graphs for the splitting process $a_{0}\to a_{1}a_{2}a_{3}$, where
$a_{3}$ is the spectator parton. More precisely, these kernels arise from
kinematic configurations in which all three final state partons have large
transverse momenta. The finite value of ${\Delta}$ can be neglected in that
region, so that one should have active-spectator symmetry according to our
arguments in the previous paragraph. Note that divergences from configurations
in which only two of the three final state partons have transverse momenta in
the ultraviolet are not associated with $P_{s}^{(1)}$, but with the lower-
order splitting kernels $P_{s}^{(0)}$ or $P^{(0)}_{\phantom{s}}$.

Let us give another reason why the finite transverse momentum ${\Delta}$ in
our calculation should play no role for the kernels $P_{s}$. We obtained
$P_{s}^{(1)}$ from the computation of the splitting kernels $W(\Delta)$, where
${\Delta}$ plays the role of a hard scale and cannot be set to zero. On the
other hand, the evolution kernels $P_{s}$ appear in the evolution equation of
the DPDs $F(\Delta)$, which is valid at $\Delta=0$. One could hence compute
$P_{s}^{(1)}$ from graphs with $\Delta=0$, provided that one has a
dimensionful infrared regulator for separating infrared and ultraviolet poles
in dimensional regularisation.

###  5.3 Kinematic limits

In this subsection, we analyse the behaviour of the kernels $P_{s}^{(1)}$ and
$W^{(2)}$ in various kinematic limits.

####  5.3.1 Threshold limit: large $x_{1}+x_{2}$

To begin with, we consider the convolution (5) at the parton-level threshold
$x_{1}+x_{2}\to 1$, i.e. $x_{3}\to 0$. Writing out the plus-distribution, we
readily obtain

$\displaystyle K\underset{12}{\otimes}f$ $\displaystyle\underset{x_{3}\to
0}{=}{}-\bigl{[}\mskip 1.5muK_{p}(x_{1},x_{2})\,\log
x_{3}-K_{\delta}(x_{1})\mskip
1.5mu\bigr{]}\;\frac{f(x_{1}+x_{2})}{x_{1}+x_{2}}\,.$ (149)

A logarithm in $x_{3}$ is thus generated by the plus distribution term of the
kernel $K$, which hence dominates the threshold behaviour of the convolution.
Since the kernels do not contain distributions $[\mskip
1.5mux_{3}^{-1}\log^{k}x_{3}\mskip 1.5mu]_{+}$ with $k>0$, no higher powers of
$\log x_{3}$ are generated in the convolution (5). This situation is analogous
to the convolution of ordinary DGLAP splitting functions with PDFs, where no
such distributions appear even at NNLO [98, 99].

The plus distribution parts of the LO channel kernels are given by

$\displaystyle P_{q\bar{q},g;\,p}^{(1)}$ $\displaystyle=-2\mskip
1.5mu(C_{A}-2C_{F})\,T_{F}\,(x_{1}\bar{x}_{1}+x_{2}\mskip
1.5mu\bar{x}_{2})\,,$ $\displaystyle P_{qg,q;\,p}^{(1)}$
$\displaystyle=2C_{A}\mskip 1.5muC_{F}\,\bar{x}_{1}\,,$ (150)

and

$\displaystyle W^{[2,2]}_{gg,g;\,p}$
$\displaystyle=C_{A}^{2}\,\bigl{[}p_{gg}(x_{1})+p_{gg}(x_{2})\bigr{]}\,,$
$\displaystyle W^{[2,2]}_{q\bar{q},g;\,p}$
$\displaystyle=-\frac{(C_{A}-2C_{F})\,T_{F}}{2}\,\bigl{[}p_{qg}(x_{1})+p_{qg}(x_{2})\bigr{]}\,,$
$\displaystyle W^{[2,0]}_{q\bar{q},g;\,p}$ $\displaystyle=\frac{1}{2}\mskip
1.5muP_{q\bar{q},g;\,p}^{(1)}\,,$ $\displaystyle W^{[2,2]}_{qg,q;\,p}$
$\displaystyle=C_{A}\mskip 1.5muC_{F}\,p_{qq}(x_{1})\,,$ (151)

where we omitted coefficients $P_{a_{1}a_{2},a_{0};\,p}^{(1)}$ and
$W_{a_{1}a_{2},a_{0};\,p}^{[2,k]}$ that are exactly zero.

The kernels $P_{s}^{(1)}$ and $W^{(2)}$ for NLO channels have no plus
distribution terms. Their leading threshold behaviour is $\log^{n}x_{3}$ with
$n=0,1,2$, which gives a convolution $K\otimes f$ of order $x_{3}\mskip
1.5mu\log^{n}x_{3}$. Specifically, we have

$\displaystyle P_{gg,g}^{(1)}$ $\displaystyle\sim\mathcal{O}(\log x_{3})\,,$
$\displaystyle P_{qg,g}^{(1)}$ $\displaystyle\sim\mathcal{O}(\log x_{3})\,,$
$\displaystyle P_{gg,q}^{(1)}$ $\displaystyle\sim\mathcal{O}(\log x_{3})\,,$
$\displaystyle P_{qq^{\prime}\\!,q}^{(1)}$ $\displaystyle\sim\mathcal{O}(\log
x_{3})\,,$ $\displaystyle P_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}^{(1)}$
$\displaystyle\sim\mathcal{O}(\log x_{3})\,,$ $\displaystyle
P_{qq,q}^{v\,(1)}$ $\displaystyle\sim\mathcal{O}(\log x_{3})\,,$
$\displaystyle P_{q\bar{q},q}^{v\,(1)}$ $\displaystyle\sim\mathcal{O}(1)\,,$
(152)

and

$\displaystyle W^{[2,0]}_{qg,g}$
$\displaystyle\sim\mathcal{O}(\log^{2}x_{3})\,,$ $\displaystyle
W^{[2,1]}_{qg,g}$ $\displaystyle\sim\mathcal{O}(\log x_{3})\,,$ $\displaystyle
W^{[2,2]}_{qg,g}$ $\displaystyle\sim\mathcal{O}(1)\,,$ $\displaystyle
W^{[2,0]}_{gg,q}$ $\displaystyle\sim\mathcal{O}(\log^{2}x_{3})\,,$
$\displaystyle W^{[2,1]}_{gg,q}$ $\displaystyle\sim\mathcal{O}(\log x_{3})\,,$
$\displaystyle W^{[2,2]}_{gg,q}$ $\displaystyle\sim\mathcal{O}(1)\,,$
$\displaystyle W^{[2,0]}_{qq^{\prime}\\!,q}$
$\displaystyle\sim\mathcal{O}(\log^{2}x_{3})\,,$ $\displaystyle
W^{[2,1]}_{qq^{\prime}\\!,q}$ $\displaystyle\sim\mathcal{O}(\log x_{3})\,,$
$\displaystyle W^{[2,2]}_{qq^{\prime}\\!,q}$
$\displaystyle\sim\mathcal{O}(1)\,,$ $\displaystyle
W^{[2,0]}_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}$
$\displaystyle\sim\mathcal{O}(\log^{2}x_{3})\,,$ $\displaystyle
W^{[2,1]}_{q^{\prime}\\!\bar{q}^{\prime}\\!,q}$
$\displaystyle\sim\mathcal{O}(\log x_{3})\,,$ $\displaystyle