# Evolution and interpolation of double parton distributions using Chebyshev
grids

Markus Diehl  Riccardo Nagar  Peter Plößl  and Frank J. Tackmann

(v1: May 8, 2023
v2: June 26, 2023)

###### Abstract

Double parton distributions are the nonperturbative ingredients needed for
computing double parton scattering processes in hadron-hadron collisions. They
describe a variety of correlations between two partons in a hadron and depend
on a large number of variables, including two independent renormalization
scales. This makes it challenging to compute their scale evolution with
satisfactory numerical accuracy while keeping computational costs at a
manageable level. We show that this problem can be solved using interpolation
on Chebyshev grids, extending the methods we previously developed for ordinary
single-parton distributions. Using an implementation of these methods in the
C++ library ChiliPDF, we study for the first time the evolution of double
parton distributions beyond leading order in perturbation theory.

##  1 Introduction

Double parton scattering (DPS) is a mechanism in hadron-hadron collisions in
which two partons from each hadron initiate two separate hard-scattering
processes. The Tevatron and the LHC have provided evidence for this mechanism
for a wide range of final states with heavy flavors, jets, photons, or
electroweak gauge bosons [1, 2, 3, 4, 5, 6]. While often suppressed compared
with single-parton scattering, DPS can be significant and even dominate in
specific final states or kinematic regions. A prominent example is like-sign
$W$ production [5, 7, 8, 9, 10, 11, 12, 13], which (via its leptonic decay)
also provides a background to searches for physics beyond the Standard Model
[14, 15]. Significant progress in the theory description of DPS has been made
in the last decade [16, 17, 18, 19, 20, 21, 22, 23, 24], and there is a
sustained interest in phenomenological aspects (see for instance the recent
study in [25] and the references therein). A detailed account of theoretical
and experimental aspects of DPS is given in ref. [26].

Double parton distributions (DPDs) are the nonperturbative quantities needed
to compute DPS cross sections. They quantify a variety of correlations between
two partons in the proton [19, 20, 27] and thus reveal aspects of hadron
structure that are not accessible in ordinary parton distributions (PDFs). Our
knowledge of DPDs is still quite limited, although there is considerable
activity in devising theory guided ansätze [28, 29, 30, 31] and in computing
DPDs using lattice QCD [32] or quark models [33, 34, 35, 36, 37, 38, 39, 40].

Reflecting the multitude of information they contain, DPDs depend on a large
number of variables, namely the momentum fractions $x_{1}$ and $x_{2}$ of the
two partons, the transverse distance $y$ between them, and the renormalization
scales $\mu_{1}$ and $\mu_{2}$ associated with each parton (corresponding to
the factorization scales in the two hard-scattering processes that the partons
initiate). Their scale dependence is well understood (see eq. (2.1) below),
but the numerical delivery of evolved DPDs is a computational challenge. We
believe that this represents one of the bottlenecks for using realistic forms
of DPDs in predictions, either at the analytical level or in the form of Monte
Carlo event generators [41, 42, 43, 44, 45]. In the present work, we present a
method that addresses this challenge, allowing for high numerical accuracy at
a moderate computational cost.

We recall that the evolution equations for PDFs can be solved numerically
either by discretization in the momentum fraction $x$ or in Mellin space, and
that public codes are available for both options, see for instance refs. [46,
47, 48, 49, 50] and [51, 52, 53]. For DPDs, even the simplest realistic
initial conditions for evolution have a correlated dependence on $x_{1}$ and
$x_{2}$ (see eqs. (7) and (8) below), which precludes the analytic computation
of complex Mellin moments in these two variables. As far as we can see, this
makes a discretization in $x_{1}$ and $x_{2}$ inevitable.

The main reason why the handling of evolved distributions is much more
demanding for DPDs than for PDFs is the sheer amount of computer memory needed
to store the former. As an example, consider interpolation grids for $x_{1}$
and $x_{2}$ with $64$ points (for PDFs, this is at the lower end of typical
grid sizes in the LHAPDF interface [54]), and grids with $48$ points for $y$,
$\mu_{1}$, and $\mu_{2}$. This corresponds to $64^{2}\times 48^{3}\approx
4.5\times 10^{8}$ real numbers and thus (with 8 bytes for a double precision
floating-point number) to $3.4\operatorname{GiB}$ per parton flavor
combination. A full unpolarized DPD with $5$ active quark flavors has
$11^{2}=121$ parton flavor combinations. Getting evolved DPDs from pre-
computed grids would thus require an enormous amount of memory (with an added
penalty of accessing this memory for interpolation in the five variables). The
path we have chosen instead is to store a DPD only for a pair of initial
scales, and to evolve “on the fly” to the desired values of $\mu_{1}$ and
$\mu_{2}$. In the above example, this reduces the memory imprint to
$64^{2}\times 48$ points and thus to $1.5\operatorname{MiB}$ per flavor
combination.

We will see in subsection 3.4 that the computational effort for evolving a DPD
with $p_{x}$ grid points in $x_{1}$ and $x_{2}$ scales like $p_{x}^{3}$,
whilst the scaling with the number $p_{y}$ of points in $y$ is at most linear.
It is therefore of great importance to limit the size of interpolation grids
for the momentum fractions, without sacrificing numerical accuracy. In our
previous work [55] we found that a highly accurate interpolation and evolution
of PDFs in $x$ space is possible using Chebyshev grids. In the present paper,
we show how to adapt this approach to the case of DPDs. We find that with
$p_{x}\sim 54$ and $p_{y}\sim 40$ to $56$ points, one can achieve high
accuracy for the interpolation and integration of DPDs with $x_{1},x_{2}\geq
10^{-5}$ and $y\geq(7\,\mathrm{TeV})^{-1}$.

We have implemented our approach in the C++ library ChiliPDF,111Chebyshev
Interpolation Library for PDFs which is under development and which will be
made public in the future. Our implementation demonstrates that the methods we
will describe do work in practice. It also allows us to explore a number of
physics aspects. So far, DPD evolution has been studied only with DGLAP
kernels at leading order (LO). We extend this to next-to-leading and next-to-
next-to-leading order (NLO and NNLO). Moreover, our implementation includes
the change in the number of active quark flavors at specified matching scales,
and access to distributions with different scales ($\mu_{1}$ and $\mu_{2}$)
and flavor numbers ($n_{1}$ and $n_{2}$) for the two partons. Last but not
least, we can evolve DPDs for polarized partons, which may for instance have
measurable impact on observables in like-sign $W$ production according to the
studies in refs. [12, 13].

This paper is organized as follows. In section 2.1 we recall some basics about
DPDs and set up the corresponding notation. Using our implementation for a
physics study, we compare in section 2.2 the evolution of DPDs at LO, NLO, and
NNLO. In section 3, we summarize the method we developed in ref. [55] for
interpolating and evolving PDFs and describe its extension to DPDs. We study
the joint interpolation in the two momentum fractions $x_{1}$ and $x_{2}$ in
section 4, and the evolution in the two scales $\mu_{1}$ and $\mu_{2}$ in
section 5. Section 6 is devoted to the interpolation in the distance $y$
between the partons, from some minimum value to infinity. In section 7 we
cross check our implementation against the DPD evolution code DOVE, which was
originally introduced in ref. [28] and further developed in subsequent work.
We conclude in section 8. In appendix A we prove an important statement about
the independence of DPD evolution and flavor matching on the path taken in the
$(\mu_{1},\mu_{2})$ plane.

##  2 Basics and main results

###  2.1 Theory framework

To begin with, let us recall some properties of DPDs that are relevant to this
work. A DPD $F^{n_{1},n_{2}}_{a_{1}a_{2}}(x_{1},x_{2},y;\mu_{1},\mu_{2})$
depends on the longitudinal momentum fractions $x_{1}$ and $x_{2}$ of the two
partons, on their distance $y$ in the transverse plane, and on the
factorization scales $\mu_{1}$ and $\mu_{2}$ associated with each parton. The
labels $a_{1}$ and $a_{2}$ specify the flavor and polarization of each parton.
DPDs involving transverse quark or linear gluon polarization have open Lorenz
indices and depend on the difference $\vec{y}$ of the transverse parton
positions rather than on the length $y$ of this vector. As specified in
section 2 of ref. [56], they can be decomposed into scalar distributions that
depend on $y$. The evolution and matching equations discussed below hold
separately for these scalar distributions. Throughout this work, we consider
only DPDs in the color-singlet channel, i.e. where the color is summed over
separately for each of the two extracted partons.

Notice that we allow different numbers $n_{1}$ and $n_{2}$ of active quark
flavors for the two partons, which is useful when the hard-scattering
processes initiated by $a_{1}$ and $a_{2}$ take place at very different
scales. Technically, $n_{1}$ and $n_{2}$ are specified by the renormalisation
prescription for the twist-two operator pertaining to parton $a_{1}$ and
$a_{2}$, respectively.

The scale dependence of a DPD is given by the evolution equations

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\ln\mu_{1}^{2}}\,F^{n_{1},n_{2}}_{a_{1}a_{2}}(x_{1},x_{2},y;\mu_{1},\mu_{2})$
$\displaystyle=\sum_{b_{1}}\Bigl{[}P_{a_{1}b_{1}}^{n_{1}}(\mu_{1})\underset{1}{\otimes}F^{n_{1},n_{2}}_{b_{1}a_{2}}(x_{2},y;\mu_{1},\mu_{2})\Bigr{]}(x_{1})\,,$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\ln\mu_{2}^{2}}\,F^{n_{1},n_{2}}_{a_{1}a_{2}}(x_{1},x_{2},y;\mu_{1},\mu_{2})$
$\displaystyle=\sum_{b_{2}}\Bigr{[}P_{a_{2}b_{2}}^{n_{2}}(\mu_{2})\underset{2}{\otimes}F^{n_{1},n_{2}}_{a_{1}b_{2}}(x_{1},y;\mu_{1},\mu_{2})\Bigr{]}(x_{2})$
(1)

with separate Mellin convolutions

$\displaystyle\Bigr{[}P(\mu_{1})\underset{1}{\otimes}F(x_{2},y;\mu_{1},\mu_{2})\Bigr{]}(x_{1})$
$\displaystyle=\int_{x_{1}}^{1}\frac{\mathrm{d}z}{z}\;P(z;\mu_{1})\,F\Bigl{(}\frac{x_{1}}{z},x_{2},y;\mu_{1},\mu_{2}\Bigr{)}\,,$
$\displaystyle\Bigr{[}P(\mu_{2})\underset{2}{\otimes}F(x_{1},y;\mu_{1},\mu_{2})\Bigr{]}(x_{2})$
$\displaystyle=\int_{x_{2}}^{1}\frac{\mathrm{d}z}{z}\;P(z;\mu_{2})\,F\Bigl{(}x_{1},\frac{x_{2}}{z},y;\mu_{1},\mu_{2}\Bigr{)}$
(2)

for the two momentum fractions. Although we define each convolution integral
in eq. (2.1) with the upper integration limit equal to $1$, the effective
integration range is reduced by the support property

$\displaystyle F^{n_{1},n_{2}}_{a_{1}a_{2}}(x_{1},x_{2},y;\mu_{1},\mu_{2})$
$\displaystyle=0$ $\displaystyle\text{for }x_{1}+x_{2}>1$ (3)

of DPDs, which is conserved by the evolution equations (2.1). The evolution
kernels $P^{n}_{ab}(z;\mu)$ in eq. (2.1) are identical to the familiar DGLAP
kernels for the evolution of PDFs with $n$ active quark flavors. Note also
that the distance $y$ plays no active role in the scale evolution, so that
each “slice” of a DPD at fixed $y$ evolves by itself.

We remark in passing that the homogeneous evolution equations (2.1) hold for
DPDs depending on the interparton distance $y$. An inhomogeneous term appears
if one integrates the distributions over $\vec{y}$ or performs a Fourier
transform from $\vec{y}$ to its conjugate transverse momentum [19, 24]. This
term is closely related with the splitting contribution (7) discussed below
and has been extensively studied in the earlier literature [57, 58, 59, 28,
60]. Throughout this work, we will exclusively work with DPDs that depend on
$y$ and evolve as given in eq. (2.1).

The transition from $n-1$ to $n$ active flavors in a DPD is described by
matching equations

$\displaystyle F_{a_{1}a_{2}}^{n_{1},n_{2}}(x_{1},x_{2},y;\mu_{1},\mu_{2})$
$\displaystyle=\sum_{b_{1}}\Bigl{[}A_{a_{1}b_{1}}^{n_{1}}(m_{n_{1}};\mu_{1})\underset{1}{\otimes}F_{b_{1}a_{2}}^{n_{1}-1,n_{2}}(x_{2},y;\mu_{1},\mu_{2})\Bigr{]}(x_{1})\,,$
$\displaystyle F_{a_{1}a_{2}}^{n_{1},n_{2}}(x_{1},x_{2},y;\mu_{1},\mu_{2})$
$\displaystyle=\sum_{b_{2}}\bigl{[}A_{a_{2}b_{2}}^{n_{2}}(m_{n_{2}};\mu_{2})\underset{2}{\otimes}F_{a_{1}b_{2}}^{n_{1},n_{2}-1}(x_{1},y;\mu_{1},\mu_{2})\Bigr{]}(x_{2})\,,$
(4)

where $m_{n}$ is the mass of the quark with flavor number $n$ and the matching
kernels $A^{n}_{ab}(z,m_{n};\mu)$ are the same as the ones for PDFs [61]:

$\displaystyle f_{a}^{n}(x;\mu)$
$\displaystyle=\sum_{b}\Bigl{[}A_{ab}^{n}(m_{n};\mu)\underset{}{\otimes}f_{b}^{n-1}(\mu)\Bigr{]}(x)\,.$
(5)

Eq. (2.1) is a coupled set of differential equations in $\mu_{1}$ and
$\mu_{2}$. In appendix A we show that with an initial condition at some point
$(\mu_{01},\mu_{02})$ one obtains a unique DPD regardless of the evolution
path taken in the $(\mu_{1},\mu_{2})$ plane. We will also show that the result
of successive flavor matching does not depend on the order in which the flavor
thresholds for the two partons are crossed. Note that the path independence of
evolution and flavor matching in the $(\mu_{1},\mu_{2})$ plane is exact even
if the perturbative expansions of the DGLAP and matching kernels are
truncated. This is in contrast to the independence on the scale at which the
flavor matching is performed, which only holds to the perturbative order of
the truncation.

The following evolution and matching kernels are currently implemented in
ChiliPDF and used in the present work:

  * •

DGLAP evolution for unpolarized and longitudinally polarized partons up to
NNLO (order $\alpha_{s}^{3}$). For the NNLO kernels we use the approximate
parameterized forms given in refs. [62, 63] and refs. [64, 65]. The kernels in
these references were confirmed by the calculations in refs. [66, 67].

  * •

DGLAP evolution up to NLO for transversely polarized quarks and linearly
polarized gluons, with the kernels given in refs. [68, 69].

  * •

Flavor matching for unpolarized partons up to NNLO, i.e. order
$\alpha_{s}^{2}$. The corresponding kernels were computed independently in
refs. [61] and [70, 71, 72], and we verified that they agree with each other.

  * •

Flavor matching for polarized partons up to NLO. At order $\alpha_{s}$, the
matching kernels are directly proportional to the DGLAP kernels for the same
polarization.

We now take a closer look at the dependence of DPDs on the kinematic variables
$x_{1},x_{2}$, and $y$. The $y$ dependence is not directly observable, because
DPDs enter cross sections (differential or integrated) via so-called double
parton luminosities

$\displaystyle\mathcal{L}_{a_{1}a_{2},b_{1}b_{2}}^{n_{1},n_{2}}(x_{1},x_{2},\bar{x}_{1},\bar{x}_{2};\mu_{1},\mu_{2},y_{\min})$
$\displaystyle\quad=\int\mathrm{d}^{2}\vec{y}\;\,\theta(y-y_{\min})\;F_{a_{1}a_{2}}^{n_{1},n_{2}}(x_{1},x_{2},y;\mu_{1},\mu_{2})\,F_{b_{1}b_{2}}^{n_{1},n_{2}}(\bar{x}_{1},\bar{x}_{2},y;\mu_{1},\mu_{2})\,,$
(6)

where $y$ is integrated over with a lower cutoff $y_{\min}$ of order
$1/\min(\mu_{1},\mu_{2})$. The origin of this cutoff is explained in ref.
[24], where it is also shown how the dependence on this cutoff is removed when
the cross sections for single and double parton scattering are combined.

The form in eq. (2.1) holds for unpolarized and longitudinally polarized
partons; for transverse quark or linear gluon distributions additional tensors
depending on $\vec{y}$ appear under the integral.

At sufficiently small distances $y$, the dominant contribution to a DPD arises
from the perturbative splitting of a single parton into the two observed
partons $a_{1}$ and $a_{2}$. At leading order in the coupling, this mechanism
gives [19, 24]

$\displaystyle F_{a_{1}a_{2}}^{\text{spl}}(x_{1},x_{2},y;\mu_{y},\mu_{y})$
$\displaystyle=\frac{1}{\pi
y^{2}}\,\frac{\alpha_{s}(\mu_{y})}{2\pi}\,V^{(1)}_{a_{1}a_{2},a_{0}}\biggl{(}\frac{x_{1}}{x_{1}+x_{2}}\biggr{)}\,\frac{f_{a_{0}}(x_{1}+x_{2};\mu_{y})}{x_{1}+x_{2}}\,,$
(7)

where the scale $\mu_{y}$ should be taken of order $1/y$ so as to avoid large
logarithms $\ln(y\mu)$ in higher-order corrections. In addition, there is an
“intrinsic” two-parton contribution, which lacks the $1/y^{2}$ enhancement of
the splitting part at small $y$ but grows more strongly at small $x_{1}$ and
$x_{2}$ for gluons and sea quarks.

At distances $y$ of nonperturbative size, our knowledge of DPDs is rather
limited, and in practice one needs to make a model ansatz. Where needed, we
will take recourse to the model used in refs. [24, 30], where a DPD is written
as a sum $F^{\text{int}}+F^{\text{spl}}$ of intrinsic and splitting parts at
all $y$. At a starting scale $\mu_{0}$, the intrinsic part of the DPD is
assumed to be proportional to the product of two PDFs:

$\displaystyle F_{a_{1}a_{2}}^{\text{int}}(x_{1},x_{2},y;\mu_{0},\mu_{0})$
$\displaystyle=\frac{1}{4\pi
h_{a_{1}a_{2}}}\,\exp\biggl{[}-\frac{y^{2}}{4h_{a_{1}a_{2}}}\biggr{]}\,\rho(x_{1},x_{2})\,f_{a_{1}}(x_{1};\mu_{0})\,f_{a_{2}}(x_{2};\mu_{0})\;\,,$
(8)

with a phase space factor

$\displaystyle\rho(x_{1},x_{2})$
$\displaystyle=\frac{(1-x_{1}-x_{2})^{r}}{(1-x_{1})^{r}(1-x_{2})^{r}}$ (9)

to ensure that the DPD smoothly goes to zero at the kinematic boundary
$x_{1}+x_{2}\to 1$. Other models in the literature take only a power
$(1-x_{1}-x_{2})^{r}$; for the motivation of the form (9) we point to section
3.2 in ref. [28].

The splitting term in the model of refs. [24, 30] is given by the perturbative
splitting expression (7) multiplied with the same exponential factor
$\exp\bigl{[}-y^{2}/(4h_{a_{1}a_{2}})\bigr{]}$ as in the intrinsic part. This
retains the correct small-$y$ limit whilst providing a more realistic decrease
at large $y$.

###  2.2 Quantitative impact of higher orders in DPD evolution

With the methods presented in this work, we are able to evolve DPDs at
different perturbative orders, for all relevant polarization combinations and
with different scales and active flavor numbers of the two partons. As a
demonstration, we now compare DPDs evolved at LO, NLO, and NNLO from a common
starting condition.

As starting condition we take the perturbative splitting form (7) with
$n_{1}=n_{2}=3$ active flavors. The splitting formula is evaluated at
$\mu_{y}=2\,\mathrm{GeV}$ and $y\approx 0.561\,\mathrm{GeV}^{-1}$, which
corresponds to the scale choice $\mu_{y}=b_{0}/y$ with
$b_{0}=2e^{-\gamma_{E}}\approx 1.12$, where $\gamma_{E}$ is the Euler-
Mascheroni constant. This choice is suggested by the form of the NLO
corrections to the splitting formula (7) and was found to keep these
corrections at a moderate size [73, 74]. For the PDFs in the DPD splitting
formula, we take the default NNLO set of the MSHT20 parameterization [75].222
We use the central PDF from the MSHT20nnlo_as118 set obtained via LHAPDF [54].
The ChiliPDF library includes an interface for importing data from LHAPDF. The
parameterization assumes quark masses $m_{c}=1.4\,\mathrm{GeV}$,
$m_{b}=4.75\,\mathrm{GeV}$ and a five-flavor coupling
$\alpha_{s}(m_{Z})=0.118$.

We evolve these initial conditions to $\mu_{1}=9\,\mathrm{GeV}$ with $n_{1}=3$
for the first parton and to $\mu_{2}=m_{W}$ with $n_{2}=5$ for the second one.
This is a setting relevant for the production of a $J/\mskip-1.5mu\Psi$ and a
$W$, a channel for which double parton scattering has in fact been observed
experimentally [76, 4]. Flavor matching for the second parton is performed at
the charm and bottom quark masses, with kernels of the same order (LO, NLO, or
NNLO) as the evolution kernels. An exception is longitudinally polarized
evolution at NNLO, where we use NLO flavor matching.

(a) $gg$ (b) $u\bar{u}$

(c) $ug$ (d) $sc$

Figure -179: Unpolarized DPDs at $x_{2}=5\times 10^{-3}$, evolved from
$(\mu_{1},\mu_{2})=(2\,\mathrm{GeV},2\,\mathrm{GeV})$ to
$(9\,\mathrm{GeV},m_{W})$ at different perturbative orders. The starting
conditions and the procedure of flavor matching are specified in the text. The
ratio shown in the small panels is taken with respect to the LO result.

We show the effect of the higher-order DGLAP evolution for a selection of
unpolarized DPDs in figure -179. The distributions are given in the top panel
of each plot, and their ratios with respect to the leading-order distribution
in the lower panel. At small momentum fractions, higher-order evolution has an
effect of order $10\%$ on $F_{gg}$, whereas for the other flavor combinations
shown in the figure, the effect is of order $1$. In all cases, we find that
the change from LO to NLO is much larger than the rather moderate change from
NLO to NNLO.

(a) $\Delta g\Delta g$ (b) $\Delta u\Delta\bar{u}$

(c) $\delta g\delta g$ (d) $\delta u\delta\bar{u}$

Figure -178: As figure -179 but for polarized DPDs. Panels (a) and (b) are for
longitudinal, panel (c) for linear, and panel (d) for transverse polarization.

A selection of polarized DPDs is shown in figure -178. For linear gluon and
transverse quark polarization, we have only implemented evolution at LO and
NLO. We find that the NLO corrections to evolution introduce effects of order
$10\%$ for $\Delta u\Delta\bar{u}$ and $\delta u\delta\bar{u}$ in different
regions of $x_{1}$, where as for $\Delta g\Delta g$ they reach $70\%$ at small
$x_{1}$. For $\delta g\delta g$, one obtains a huge relative difference
between NLO and LO evolution in the limit $x_{1}\ll x_{2}$. This is because
for $z\to 0$ the splitting function $P_{\delta g\delta g}(z)$ vanishes like
$z$ at LO but grows like $1/z$ at NLO [69]. At both orders, the evolved
distribution at $x_{1}\ll x_{2}$ remains however small compared with its peak
around $x_{1}\sim x_{2}$.

As is the case for PDFs, DPDs are not observables, and distributions evolved
at different orders are combined with parton-level cross sections computed at
different orders. Nevertheless, the strong evolution effects we observed in
several channels and kinematic regions suggest that NLO corrections in double
parton scattering can be of substantial size.

In figure -177, we illustrate the scale dependence of DPDs, using evolution at
NNLO and flavor matching at the highest available order (NNLO for unpolarized
partons and NLO for polarized ones). We show the distributions at the starting
point $(\mu_{1},\mu_{2})=(2\,\mathrm{GeV},2\,\mathrm{GeV})$, and after
evolution to different combinations of the scales $\mu_{i}=9\,\mathrm{GeV}$ or
$\mu_{i}=m_{W}$ for the two partons, with associated flavor numbers $n_{i}=3$
in the first and $n_{i}=5$ in the second case. For unpolarized partons, we
find that the effect of evolution from $\mu_{1}=\mu_{2}=2\,\mathrm{GeV}$ to
$9\,\mathrm{GeV}$ is huge at small $x$, turning a shape that is flat or
increasing with $x$ at the low scale into a shape decreasing with $x$ at the
higher scale. Subsequent evolution to $m_{W}$ for one or both partons has a
significant, although less dramatic, effect. For polarized partons, evolution
is generally weaker, and we find DPDs increasing with $x$ at small $x$ even
when both partons are taken at the high scale $m_{W}$.

(a) $gg$ (b) $u\bar{u}$

(c) $\Delta g\Delta g$ (d) $\Delta u\Delta\bar{u}$

Figure -177:  Comparison of DPDs at different scales. The starting conditions
at $(\mu_{1},\mu_{2})=(2\,\mathrm{GeV},2\,\mathrm{GeV})$ and
$(n_{1},n_{2})=(3,3)$ are the same as those used for the previous two figures.
At higher scales, the number of active flavors is $n_{i}=3$ for
$\mu_{i}=9\,\mathrm{GeV}$ and $n_{i}=5$ for $\mu_{i}=m_{W}$. Evolution and
flavor matching are carried out at the highest available order.

##  3 Interpolation and evolution of DPDs in ChiliPDF

In this section we sketch the mathematical framework underlying ChiliPDF,
summarizing the relevant methods presented in ref. [55] and extending them
from PDFs to DPDs. In particular, we give explicit formulae for interpolation
and integration, and for the discretization of the DGLAP evolution equations.
A detailed account of Chebyshev interpolation and its applications can be
found in ref. [77].333The first chapters of this book are available on
https://people.maths.ox.ac.uk/trefethen/ATAP.

###  3.1 Chebyshev interpolation

We start by reviewing Chebyshev interpolation of a function $f(t)$ on the
interval $t\in[-1,1]$. The function is discretized on a grid consisting of the
so-called _Chebyshev points_ , which for given $N$ read

$\displaystyle t_{i}$ $\displaystyle=\cos\Bigl{(}\frac{i\pi}{N}\mskip
1.5mu\Bigr{)}\qquad\text{with}\quad i=0,\ldots,N\,.$ (10)

These are the points where the $N^{\mathrm{th}}$ Chebyshev polynomial of the
first kind, $T_{N}(t)$, assumes its maxima $+1$ and minima $-1$. They form a
descending series from $t_{0}=1$ to $t_{N}=-1$ and satisfy the symmetry
property $t_{N-i}=-t_{i}$. We call the set of Chebyshev points a _Chebyshev
grid_.

A sufficiently smooth function $f(t)$ on $t\in[-1,1]$ can be approximated by
the _Chebyshev interpolant_ $p_{N}(t)$. This is a series of Chebyshev
polynomials $T_{j}(t)$ with $j=0,\ldots,N$ whose expansion coefficients are
such that one has $p_{N}(t_{i})=f(t_{i})$ at the Chebyshev points. For details
see equations (2.6) to (2.8) in ref. [55].

#### Barycentric formula.

A simple and efficient way to compute the Chebyshev interpolant without having
to evaluate Chebyshev polynomials is given by the _barycentric formula_

$\displaystyle p_{N}(t)$ $\displaystyle=\sum_{i=0}^{N}f(t_{i})\,b_{i}(t)$ (11)

with the barycentric basis functions

$\displaystyle b_{i}(t)$
$\displaystyle=\beta_{i}\,\frac{(-1)^{i}}{t-t_{i}}\Bigg{/}\sum_{j=0}^{N}\beta_{j}\,\frac{(-1)^{j}}{t-t_{j}}\,,$
(12)

where $\beta_{0}=\beta_{N}=1/2$ and $\beta_{i}=1$ otherwise. Even though it is
not directly evident from eq. (12), the $b_{i}(t)$ are polynomials of order
$N$. The number of operations for evaluating the barycentric formula scales
linearly with $N$. The barycentric formula is found to be numerically stable
in the interpolation interval (however, this is _not_ the case for
extrapolating the function $f(t)$ outside this interval, see chapter 5 of ref.
[77]).

#### Accuracy estimate.

A computationally inexpensive estimate for the accuracy of Chebyshev
interpolation is obtained by interpolating $f(t)$ on the Chebyshev grid
_without_ the end points $t_{0}$ and $t_{N}$. The resulting interpolant
$q_{N-2}(t)$ is a polynomial of order $N-2$ and can be computed with an
adapted form of the barycentric formula, see equations (2.23) to (2.27) in
ref. [55]. Comparing $p_{N}(t)$ with $q_{N-2}(t)$ yields an estimate for the
interpolation accuracy, i.e. for the difference between $p_{N}(t)$ and $f(t)$.
This estimate was found to work well for PDF interpolation [55], and we will
investigate in later chapters how reliable it is in the case of DPDs.

#### Integration.

Using the expansion of $f$ in terms of Chebyshev polynomials $T_{i}(t)$ and
the analytic form for integrals of these polynomials, one readily obtains the
integration rule

$\displaystyle\int_{-1}^{1}\mathrm{d}t\,f(t)$
$\displaystyle\approx\sum_{i=0}^{N}\,w_{i}\mskip 1.5muf(t_{i})\;$ (13)

with weights

$\displaystyle w_{i}$
$\displaystyle=\frac{4\beta_{i}}{N}\,\sum_{\genfrac{}{}{0.0pt}{}{j=0}{\text{even}}}^{N}\,\beta_{j}\,\frac{\cos(j\mskip
1.5mu\theta_{i})}{1-j^{2}}\,.$ (14)

This integration rule is known as _Clenshaw-Curtis quadrature_. For a detailed
discussion of its accuracy (and comparison with Gauss quadrature), we refer to
chapter 19 of ref. [77] and to section 2 of ref. [55].

An estimate for the integration accuracy can be obtained by using a quadrature
rule on the Chebyshev grid without its end points, in full analogy to the
procedure described above for interpolation. The result is known as _Fejér’s
second rule_. Details are given in equations (2.31) to (2.33) of ref. [55].

###  3.2 Interpolation strategy

In this subsection, we consider the scales $\mu_{1}$ and $\mu_{2}$ to be fixed
and omit them as arguments for brevity. We discuss the discretization of DPDs
in the three variables $x_{1}$, $x_{2}$, and $y$, which we generically denote
by $z$. Depending on the typical dependence of a DPD on $z$, different
variable transformations $z\to u$ are used to map an interval
$[z_{\min},z_{\max}]$ onto a finite interval $[u_{\min},u_{\max}]$. For the
momentum fractions $x_{1}$ and $x_{2}$ we always use the transformation

$\displaystyle u(x_{i})$ $\displaystyle=\ln x_{i}$ $\displaystyle\text{for
}i=1,2,$ (15)

which corresponds to the choice for interpolating PDFs in ref. [55]. For both
$x_{1}$ and $x_{2}$, the interpolation interval has a fixed upper limit
$x_{\max}=1$, whereas the lower limits may be equal or different (but must be
nonzero). For the $y$ dependence (where $y_{\max}$ may be infinite or finite)
we use a number of different transformations, which are presented in
subsection 6.1.

For reasons given below, the $z$ intervals are usually split into a few
subintervals. On each subinterval, we perform a linear transformation from $u$
to a variable $t\in[-1,1]$, which is used for Chebyshev interpolation as
described in subsection 3.1.

Consider a particular subgrid with $N+1$ points $u_{0},\ldots,u_{N}$, which is
mapped by a linear transform onto the Chebyshev grid $t_{0},\ldots,t_{N}$
given in eq. (10). The corresponding grid points in the physical variable are
then given by the inverse variable transformation $z_{i}=z(u_{i})$. Similarly
to the PDF case, it is advantageous to interpolate DPDs scaled by their
momentum fractions,

$\displaystyle\widetilde{F}(x_{1},x_{2},y)$ $\displaystyle=x_{1}\mskip
1.5mux_{2}\mskip 1.5muF(x_{1},x_{2},y)\,,$ (16)

because this leads to less steep functions in $x_{1}$ and $x_{2}$. The
interpolation in one variable is given by the barycentric formula

$\displaystyle\widetilde{F}(z,\ldots)$
$\displaystyle\approx\sum_{i=0}^{N}\widetilde{F}_{i}(\ldots)\;b_{i}(u(z))$
$\displaystyle\text{for~{}}z_{0}\leq z\leq z_{N}\,,$ (17)

where the ellipsis denotes the other variables and

$\displaystyle\widetilde{F}_{i}(\ldots)$
$\displaystyle=\widetilde{F}(z_{i},\ldots)\,,$ $\displaystyle b_{i}(u)$
$\displaystyle=\beta_{i}\,\frac{(-1)^{i}}{u-u_{i}}\Bigg{/}\sum_{j=0}^{N}\beta_{j}\,\frac{(-1)^{j}}{u-u_{j}}\,.$
(18)

The interpolation of the full DPD is obtained by discretization in both
momentum fractions and the interparton distance and reads

$\displaystyle\widetilde{F}(x_{1},x_{2},y)$
$\displaystyle\approx\sum_{i=0}^{N_{x\mskip-1.5mu,1}}\sum_{j=0}^{N_{x\mskip-1.5mu,2}}\sum_{k=0}^{N_{y}}\widetilde{F}_{ijk}\;b_{i}\bigl{(}\ln(x_{1})\bigr{)}\,b_{j}\bigl{(}\ln(x_{2})\bigr{)}\,b_{k}\bigl{(}u(y)\bigr{)}\,,$
(19)

where $N_{x\mskip-1.5mu,1}$, $N_{x\mskip-1.5mu,2}$, and $N_{y}$ specify the
polynomial order of interpolation in each variable, and

$\displaystyle\widetilde{F}_{ijk}$ $\displaystyle=\widetilde{F}(x_{1,i},\mskip
1.5mux_{2,j},\mskip 1.5muy_{k})$ (20)

denotes the values of $\widetilde{F}$ on the interpolation grid. Eq. (18)
reflects that the form of the barycentric basis functions (12) remains
unchanged under a linear transform of the interpolation variable. Analogous
formulae can be used to interpolate functions derived from DPDs, such as the
Mellin convolutions of DPDs with an integral kernel (see subsection 3.3).

Note that we interpolate $\widetilde{F}(x_{1},x_{2},y)$ on a full rectangle in
the $(x_{1},x_{2})$ plane, given by $x_{1}\in[x_{1,\min},1]$ and
$x_{2}\in[x_{2,\min},1]$. This includes the region $x_{1}+x_{2}>1$ where the
DPD is zero. In this way, we can use the same interpolation grid in $x_{1}$
for all values of $x_{2}$ and vice versa, which in turn permits a simple
implementation of DGLAP evolution for the two partons. The implication of this
choice on the interpolation accuracy is discussed in section 4.

Chebyshev interpolation is “global” in the sense that the interpolant at a
given point $z$ depends on the function values at _all_ points $z_{i}$ in the
interpolation interval. For grids with many points, this can lead to an
accumulation of rounding errors, especially in regions of $z$ where the
function is much smaller than elsewhere in the interval. To limit this effect,
it is advantageous to split the full $z$ domain into several subintervals and
to use the barycentric formula (19) separately on each subinterval. In ref.
[55] this was found to substantially improve the interpolation accuracy for
PDFs, provided that the number of points in each subinterval does not become
too small. We find that the same is true for the interpolation of DPDs in
$x_{1}$ and $x_{2}$, as well as in $y$.

To specify such a composite grid with $k$ individual subgrids, we use the
notation

$\displaystyle[z_{0},\mskip 1.5muz_{1},\mskip 1.5mu\ldots,\mskip
1.5muz_{k}]_{(p_{1},\mskip 1.5mup_{2},\mskip 1.5mu\ldots,\mskip
1.5mup_{k})}\,,$ (21)

where $z_{i}$ are the subinterval boundaries and $p_{i}=N_{i}+1$ is the number
of Chebyshev points in subgrid number $i$. We refer to this as a
$(p_{1},p_{2},\ldots,p_{k})$-point grid. Note that neighboring subgrids share
their end points, such that the total number of grid points is
$p_{\text{total}}=\sum_{i=1}^{k}p_{i}-(k-1)$.

Splitting the interpolation interval in $y$ has the additional benefit that
one can use different variable transformations $u(y)$ for regions with
different characteristic $y$ dependence, such as the power law behavior from
perturbative splitting and an exponential decrease at large distances. Details
will be discussed in subsection 6.1.

###  3.3 Mellin convolution

Consider now the convolution of a DPD with an integral kernel, which may be a
DGLAP evolution kernel, a flavor matching kernel, or a hard-scattering
coefficient. The general expression to be computed reads

$\displaystyle\bigl{(}K\underset{}{\otimes}\widetilde{F}\,\bigr{)}(x)$
$\displaystyle=\int_{x}^{1}\frac{\mathrm{d}z}{z}\,K(z)\,\widetilde{F}\left(\frac{x}{z}\right)\,,$
(22)

where $x$ stands for $x_{1}$ or $x_{2}$ and the dependence of the rescaled DPD
on all other variables is omitted for now. The kernel $K(z)$ is appropriately
rescaled, given that $g_{1}=g_{2}\otimes g_{3}$ implies $h_{1}=h_{2}\otimes
h_{3}$ with $h_{k}(z)=z\mskip 1.5mug_{k}(z)$.

To discretize eq. (22), we consider for simplicity a single Chebyshev grid
with $p_{x}$ points for $x\in[x_{\min},1]$. The result of the convolution
$(K\otimes\widetilde{F})(x)$ is a function defined in the same $x$ interval as
$\widetilde{F}(x)$ and can hence be interpolated on the same grid. It is thus
sufficient to evaluate eq. (22) at the grid points $x_{i}$,

$\displaystyle\bigl{(}K\underset{}{\otimes}\widetilde{F}\,\bigr{)}(x_{i})$
$\displaystyle=\int_{x_{i}}^{1}\\!\frac{\mathrm{d}{z}}{z}\,K(z)\,\widetilde{F}\left(\frac{x_{i}}{z}\right)\approx\int_{x_{i}}^{1}\frac{\mathrm{d}{z}}{z}\,K(z)\,\sum_{i^{\prime}=0}^{p_{x}-1}\widetilde{F}(x_{i^{\prime}})\;b_{i^{\prime}}\\!\left(\ln\frac{x_{i}}{z}\right)\,,$
(23)

where in the second step we have approximated $\widetilde{F}(x_{i}/z)$ by its
interpolated form. Here and in the following, we use the number of grid points
$p_{x}$ rather than the polynomial order $N_{x}=p_{x}-1$ to specify summation
ranges. At the grid points, the convolution is thus obtained by a matrix
multiplication

$\displaystyle\bigl{(}K\underset{}{\otimes}\widetilde{F}\,\bigr{)}(x_{i})$
$\displaystyle\approx\sum_{i^{\prime}=0}^{p_{x}-1}K_{ii^{\prime}}\,\widetilde{F}(x_{i^{\prime}})\,,$
(24)

where the _kernel matrix_

$\displaystyle K_{ii^{\prime}}$
$\displaystyle=\int_{x_{i}}^{1}\frac{\mathrm{d}{z}}{z}\,K(z)\;b_{i^{\prime}}\\!\left(\ln\frac{x_{i}}{z}\right)$
(25)

can be computed using standard numerical integration techniques (we use an
adaptive Gauss-Kronrod routine). If $K(z)$ contains plus- or $\delta$
distributions, this requires rewriting eq. (25) as specified in equations
(4.5) to (4.10) of ref. [55].

It is straightforward to generalize the preceding discussion to the case where
the interval $[x_{\min},1]$ is split into several subintervals with
corresponding subgrids. The relation (24) remains valid, with an appropriately
generalized definition of the kernel matrix and with $p_{x}$ being the _total_
number of grid points in $x$.

The kernel matrix method can be readily used for flavor matching of DPDs,
discretizing the Mellin convolutions with the matching kernels in eq. (2.1).
It is also a crucial part in DPD evolution, which is presented next.

###  3.4 DGLAP evolution

We now discuss the evolution of DPDs, starting with evolution in $\mu_{1}$ at
a fixed value of $\mu_{2}$. The evolution equation for the rescaled DPD reads

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\ln\mu_{1}^{2}}\,\widetilde{F}(x_{1},x_{2},y;\mu_{1},\mu_{2})$
$\displaystyle=\Bigl{[}\widetilde{P}(\mu_{1})\underset{1}{\otimes}\widetilde{F}(x_{2},y;\mu_{1},\mu_{2})\Bigr{]}(x_{1})\,,$
(26)

where $\widetilde{P}(z;\mu_{1})=z\mskip 1.5muP(z;\mu_{1})$ is the rescaled
splitting function and we omit parton labels and the associated sums for the
time being. After discretization on one or several subgrids with a total
number of $p_{x\mskip-1.5mu,1}$, $p_{x\mskip-1.5mu,2}$, $p_{y}$ points in
$x_{1}$, $x_{2}$, $y$, the integro-differential equation (26) turns into a
linear system of ordinary differential equations in $\mu_{1}$,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\ln\mu_{1}^{2}}\,\widetilde{F}_{ijk}(\mu_{1},\mu_{2})$
$\displaystyle=\sum_{i^{\prime}=0}^{p_{x\mskip-1.5mu,1}-1}\,\widetilde{P}_{ii^{\prime}}(\mu_{1})\,\widetilde{F}_{i^{\prime}\mskip-1.5mujk}(\mu_{1},\mu_{2})\,,$
(27)

where $\widetilde{F}_{ijk}(\mu_{1},\mu_{2})$ is given in eq. (20) and
$\widetilde{P}_{ii^{\prime}}$ is the matrix for the kernel $\widetilde{P}(z)$
as defined in eq. (25) or its generalization to several subgrids. At NNLO
accuracy, we have

$\displaystyle\widetilde{P}_{ii^{\prime}}(\mu)$
$\displaystyle=\sum_{k=0}^{2}\biggl{[}\frac{\alpha_{s}(\mu)}{4\pi}\biggr{]}^{k+1}\;\widetilde{P}_{ii^{\prime}}^{(k)}\,,$
(28)

and after a change of variables from the scales $\mu_{i}$ to the evolution
times

$\displaystyle t_{i}$ $\displaystyle=-\ln\bigl{(}\alpha_{s}(\mu_{i})\bigr{)}$
$\displaystyle\text{for }i=1,2,$ (29)

we obtain

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t_{1}}\,\widetilde{F}_{ijk}(t_{1},t_{2})$
$\displaystyle=\sum_{i^{\prime}=0}^{p_{x\mskip-1.5mu,1}-1}\,K_{ii^{\prime}}(t_{1})\,\widetilde{F}_{i^{\prime}\mskip-1.5mujk}(t_{1},t_{2})$
(30)

with a new kernel matrix given by

$\displaystyle K_{ii^{\prime}}(t)$ $\displaystyle=\biggl{[}\mskip
1.5mu\beta_{0}+\frac{e^{-t}}{4\pi}\,\beta_{1}+\frac{e^{-2t}}{(4\pi)^{2}}\,\beta_{2}\,\biggr{]}^{-1}\;\biggl{[}\mskip
1.5mu\widetilde{P}_{ii^{\prime}}^{(0)}+\frac{e^{-t}}{4\pi}\,\widetilde{P}_{ii^{\prime}}^{(1)}+\frac{e^{-2t}}{(4\pi)^{2}}\,\widetilde{P}_{ii^{\prime}}^{(2)}\,\biggr{]}\,,$
(31)

where $\beta_{k}$ are the expansion coefficients of the QCD $\beta$ function
(see equation (5.4) in [55] for our normalization conventions).

To solve eq. (30) we use a Runge-Kutta algorithm with a fixed maximal step
size in $t_{1}$. Notice that the $t_{1}$ dependence of
$K_{ii^{\prime}}(t_{1})$ starts only at NLO, which leads to a more uniform
accuracy of the algorithm compared with the direct solution of eq. (27), where
$\widetilde{P}_{ii^{\prime}}(\mu_{1})$ changes with $\ln\mu_{1}$ already at
LO. As documented in appendix A of ref. [55], we find that the use of a
higher-order Runge-Kutta algorithm — in particular the $8^{\text{th}}$ order
algorithm of Dormand and Prince [78] — greatly improves the accuracy of
evolution for a given number $N_{\text{RK}}$ of intermediate time values at
which the r.h.s. of the differential equation is evaluated for evolution from
one value of $t_{1}$ to another.

The Runge-Kutta algorithm is also used to compute the transformation (29),
solving the differential equation for the running of $\alpha_{s}^{-1}$ as a
function of $\ln\mu$.

Flavor labels and sums have to be added in eq. (30) as in the original version
(2.1). The mixing between parton flavors under evolution is the same as for
PDFs, and correspondingly we solve the evolution equations in the same basis
that is used for PDFs in ref. [55]. For each of the two partons, we thus form
the linear combinations

$\displaystyle\Sigma^{\pm}=\sum_{i}q_{i}^{\pm},$ $\displaystyle
u^{\pm}-d^{\pm},$ $\displaystyle d^{\pm}-s^{\pm},$ $\displaystyle
s^{\pm}-c^{\pm},$ $\displaystyle c^{\pm}-b^{\pm},$ $\displaystyle
b^{\pm}-t^{\pm}\,,$ (32)

where $q^{\pm}=q\pm\bar{q}$. The combination $\Sigma^{-}$ and all flavor
differences then evolve by themselves, and mixing is reduced to the
combination $\Sigma^{+}$ and the gluon.

The system (30) needs to be solved separately for each combination of the
indices $j,k$ that correspond to the “inactive” variables $x_{2}$ and $y$
during evolution in $t_{1}$, and for each of the $2n_{2}+1$ active flavors of
the second parton. This gives a large number of
$p_{x\mskip-1.5mu,2}\;p_{y}\,(2n_{2}+1)$ combinations — typically several
thousand. It is much more efficient to apply the Runge-Kutta algorithm not to
each of these combinations, but to a set of $p_{x\mskip-1.5mu,1}$ basis
vectors $e_{i}^{j}$ that generates all possible initial conditions on the grid
for $x_{1}$. With the choice of basis $e_{i}^{\mskip 1.5muj}=\delta_{ij}$,
evolution from $t_{01}$ to $t_{1}$ requires solving the system

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t_{1}}\,U_{ij}(t_{1},t_{01})$
$\displaystyle=\sum_{i^{\prime}=0}^{p_{x\mskip-1.5mu,1}-1}K_{ii^{\prime}}(t_{1})\,U_{i^{\prime}j}(t_{1},t_{01})\,,$
(33)

with the boundary condition $U_{ij}(t_{01},t_{01})=\delta_{ij}$. The solution
of eq. (30) is then obtained by a matrix multiplication

$\displaystyle\widetilde{F}_{ijk}(t_{1},t_{2})$
$\displaystyle=\sum_{i^{\prime}=0}^{p_{x\mskip-1.5mu,1}-1}U_{ii^{\prime}}(t_{1},t_{01})\,\widetilde{F}_{i^{\prime}jk}(t_{01},t_{2})\,,$
(34)

where the evolution matrix $U_{ii^{\prime}}(t_{1},t_{01})$ acts as the
discretized version of the Green function for the original evolution equation
(26). With flavor labels reinstated, eq. (33) is solved separately for the
evolution kernels of the different channels (i.e. for $q_{i}^{+}-q_{j}^{+}$,
$q_{i}^{-}-q_{j}^{-}$, $\Sigma^{-}$, and the coupled system of $\Sigma^{+}$
and the gluon).

This approach is also attractive in the PDF case when a large set of different
PDFs is evolved simultaneously, and it is for instance used in the PDF
evolution codes APFEL [49] and EKO [53], with $U_{ii^{\prime}}(t_{1},t_{01})$
being called an evolution operator and an evolution kernel operator,
respectively. In ChiliPDF, the evolution of a DPD in one scale is implemented
in the same way as the evolution of a set of PDFs.

The number of operations for evaluating eq. (34) for all indices $i,j,k$ and
flavors for the second parton scales like
$(p_{x\mskip-1.5mu,1})^{2}\;p_{x\mskip-1.5mu,2}\;p_{y}\,(2n_{2}+1)$, and the
number of operations for solving the evolution equation (33) scales like
$(p_{x\mskip-1.5mu,1})^{3}\,N_{\text{RK}}$ with $N_{\text{RK}}$ introduced
below eq. (31). The number of grid points for the interpolation in the
momentum fractions is therefore a major factor determining the computational
speed of DPD evolution.

The preceding discussion is trivially extended to evolution in the scale
$\mu_{2}$ at fixed $\mu_{1}$. To evolve from $(\mu_{01},\mu_{02})$ to
$(\mu_{1},\mu_{2})$, we first evolve from $(\mu_{01},\mu_{02})$ to
$(\mu_{1},\mu_{02})$ and from there to $(\mu_{1},\mu_{2})$. We verified that
evolving to $(\mu_{01},\mu_{2})$ in the first step leads to the same final
result with very high numerical accuracy (see also subsection 5.1).

Let us briefly sketch how the initial conditions for DPD evolution are set up
in ChiliPDF. In a first step, the values of
$F_{a_{1}a_{2}}\bigl{(}x_{1},x_{2},y;\mu(y),\mu(y)\bigr{)}$ for a selected
number $n_{0}$ of active flavors are computed on the grid points in $x_{1}$,
$x_{2}$ and $y$, either from pre-defined functions for the splitting form (7),
or from functions implementing a product ansatz like the one in (8), or from
functions entirely provided by the user. As is adequate for the physics, the
input scale $\mu(y)$ may depend on $y$ as specified by the user. For each grid
point $y_{k}$, the discretized DPD is then evolved to a pair of scales
$(\mu_{0},\mu_{0})$ that is independent of $y$ and stored as an initial
condition. After this, all evolution calls can use the same set of evolution
matrices for all grid points $y_{k}$, which significantly reduces computing
time. One may also add two DPDs (after evolution to a common scale pair),
which is for instance needed for calculating $F^{\text{int}}+F^{\text{spl}}$.

In the procedure just described, the flavor number $n_{0}$ is equal for both
partons. DPDs with larger values of $n_{1}$ or $n_{2}$ are then obtained by
flavor matching, with initial conditions at the pair of scales
$(\mu_{1},\mu_{2})$ where the matching is carried out. More detail is given in
subsection 5.1.

##  4 Interpolation of DPDs in both momentum fractions

As shown in the previous section, our approach to interpolating DPDs in
$x_{1}$ and $x_{2}$ is a straightforward two-dimensional generalization of the
method presented in ref. [55] for the interpolation of PDFs in $x$. In the
present section, we study the performance of this approach. Let us point out
some particularities of the DPD case:

  1. 1.

To limit computation time, interpolation grids should not have too many
points, while still giving numerically reliable results. For PDFs one can
typically afford denser grids, and the demands on interpolation accuracy are
often higher than for DPDs.

  2. 2.

The DPD is a non-analytic function of $x_{1}$ at the point $x_{1}=1-x_{2}$
because it vanishes identically above that value. Since we interpolate on the
same grid in $x_{1}$ for all values of $x_{2}$, this singular point is
generally in the middle of an interpolation interval. Polynomial interpolation
is not very accurate in the vicinity of such a point, but we will see that
this effect is in general relatively mild.

  3. 3.

The perturbative splitting mechanism gives a dependence on $x_{1}$ that is
qualitatively different from the $x$ dependence of a PDF. In some channels,
the splitting contribution (7) goes to zero for $x_{1}\ll x_{2}$ and thus
results in a _strong_ rise of the DPD as a function of $x_{1}$.

Of course, the last two points also apply to interpolation in $x_{2}$ at a
given $x_{1}$.

In the following, we will study the interpolation accuracy for DPDs that have
a known analytic expression. To this end, we simplify the forms (8) and (7) to

$\displaystyle\widehat{F}_{a_{1}a_{2}}^{\mskip 1.5mu\text{int}\mskip
1.5mu(r)}(x_{1},x_{2})$
$\displaystyle=\frac{(1-x_{1}-x_{2})^{r}}{(1-x_{1})^{r}(1-x_{2})^{r}}\;f_{a_{1}}(x_{1})\,f_{a_{2}}(x_{2})$
(35) and $\displaystyle\widehat{F}_{a_{1}a_{2}}^{\mskip
1.5mu\text{spl}}(x_{1},x_{2})$
$\displaystyle=V^{(1)}_{a_{1}a_{2},a_{0}}\biggl{(}\frac{x_{1}}{x_{1}+x_{2}}\biggr{)}\,\frac{f_{a_{0}}(x_{1}+x_{2})}{x_{1}+x_{2}}\,,$
(36)

respectively. For the PDFs, we take the simple parameterizations of the Les
Houches benchmark PDFs given in section 1.32 of ref. [79] and section 4.4 of
ref. [80].

Corresponding to the typical use case, we take the same interpolation grids
for $x_{1}$ and $x_{2}$ (a situation in which different grids are useful is
encountered in section 7). We consider three composite grids of the form

$\displaystyle[10^{-5},5\times 10^{-3},0.5,1]_{(p,\,p,\,p)}$
$\displaystyle\text{with }p=\begin{cases}12&\text{(coarse $x$ grid)}\\\
16&\text{(medium $x$ grid)}\\\ 24&\text{(fine $x$ grid)}\end{cases}$ (37)

where the shorthand notation from eq. (21) has been used. We find that for the
given $x$ range, three subgrids with the above interval limits generally give
the best accuracy for a given total number of points. As a measure for the
quality of interpolation, we use the relative accuracy defined as

relative accuracy $\displaystyle=\bigl{|}\text{interpolated
result}/\text{exact result}-1\bigr{|}\,.$ (38)

In the following, this interpolation accuracy will be shown as a function of
$x_{1}$ for fixed values of $x_{2}$. The latter are _not_ grid points, so that
interpolation is performed for both momentum fractions.

###  4.1 Interpolation of the product form

We start with the product form (35) for the three values $r=0,1,2$ of the
power in the phase space factor. The case $r=0$ corresponds to a “naive”
product form with a discontinuity at $x_{1}+x_{2}=1$ that is not very
plausible on physics grounds. In fact, evolution to higher scales quickly
changes such a discontinuity to a steep but continuous decrease. We
nevertheless include the case $r=0$ here, in order to see how our method
performs in this limiting case. The following plots show the $u\bar{u}$
distribution; we verified that for other parton combinations the interpolation
accuracy is not significantly better or worse.

(a) $r=0,~{}x_{2}=10^{-4}$ (b) $r=0,~{}x_{2}=0.3$

(c) $r=1,~{}x_{2}=10^{-4}$ (d) $r=1,~{}x_{2}=0.3$

(e) $r=2,~{}x_{2}=10^{-4}$ (f) $r=2,~{}x_{2}=0.3$

Figure -176: Relative accuracy (38) for interpolation in $x_{1}$ and $x_{2}$,
evaluated for the product form $\widehat{F}_{u\bar{u}}^{\mskip
1.5mu\text{int}\mskip 1.5mu(r)}(x_{1},x_{2})$ of eq. (35) with $r=0,1,2$. The
interpolation grids are specified in eq. (37). Here and in similar plots, the
exact result that is being interpolated is shown in black.

In the left panels of figure -176 we see that all three grids yield an
excellent interpolation accuracy at small momentum fractions. Not
surprisingly, there is little dependence on the power $r$ in this region,
given that the phase space factor is very close to $1$ for $x_{1},x_{2}\ll 1$.
In the right panels of the figure, we see that a degradation of the relative
accuracy sets in as $x_{1}$ approaches its kinematic limit $1-x_{2}=0.7$, and
that the degradation is stronger for smaller $r$. Nevertheless, the accuracy
is better than 1% as long as the distance $1-(x_{1}+x_{2})$ from the kinematic
threshold is above $0.15$ for $r=1,2$ and above $0.2$ for $r=0$.

(a) $r=0,~{}x_{2}=0.3$ (b) $r=0,~{}x_{2}=0.3$

(c) $r=1,~{}x_{2}=0.3$ (d) $r=1,~{}x_{2}=0.3$

(e) $r=2,~{}x_{2}=0.3$ (f) $r=2,~{}x_{2}=0.3$

Figure -175: Interpolation of DPDs at large momentum fractions. The DPDs
considered are the same as in the right panels of figure -176. The plots on
the left compare the exact and interpolated forms of the DPDs. The plots on
the right show the difference between the interpolated and the exact values,
along with the exact form. Notice that the plots in the top row have a
different $x_{1}$ range than the remaining ones.

The plots in figure -175 show the $x_{1}$ region around the point where the
DPD has a non-analytic behavior, which is at $x_{1}=0.7$ in this example. In
these plots, we include the region $x_{1}>1-x_{2}$, where the interpolated DPD
oscillates around its true value zero. Note that this region is included when
one computes a Mellin convolution for the DPD with the kernel matrix method
described in subsection 3.3. However, the effect of oscillations around the
true value will partially cancel out in convolution integrals.

We see that for $r=0$, one obtains a rather poor approximation of the DPD in
the range $x_{1}\in[0.6,0.8]$ with the coarse or medium grid. In situations
where this region is deemed important, one should take a dense grid and check
the reliability of the results by comparing with an even denser grid. Of
course, this is only needed in the relevant subinterval (i.e. for
$x_{1}\in[0.5,1]$ in our example) and hence entails only a moderate increase
in the total number of grid points. For $r=1$, where the DPD is continuous but
has a discontinuous first derivative at $x_{1}=0.7$, we see that the fine grid
performs very well, and for $r=2$ already the medium grid yields a good
approximation of the DPD.

###  4.2 Interpolation of the splitting form

We now turn to the interpolation of the splitting form. We see in eq. (36)
that it decreases almost as fast as a PDF for $x_{1}+x_{2}\to 1$, so that
according to the results of the previous subsection one does not expect
particular problems with interpolation in that region. This is confirmed by
our numerical studies.

(a) $u\bar{u},~{}x_{2}=10^{-4}$ (b) $u\bar{u},~{}x_{2}=0.3$

(c) $\delta g\delta g,~{}x_{2}=10^{-4}$ (d) $\delta g\delta g,~{}x_{2}=0.3$

Figure -174: Relative accuracy (38) for interpolation in $x_{1}$ and $x_{2}$,
evaluated for the splitting form $\widehat{F}_{a_{1}a_{2}}^{\mskip
1.5mu\text{spl}}(x_{1},x_{2})$ of eq. (36). The interpolation grids are
specified in eq. (37). The plots in the top row are for an unpolarized
$u\bar{u}$ pair, and those in the bottom row for a pair of linearly polarized
gluons. Notice that in panel (d) a logarithmic scale for $x_{1}$ is used to
plot the DPD at $x_{2}=0.3$.

The shape of the splitting form is strongly influenced by the end-point
behavior of the splitting kernel $V^{(1)}_{a_{1}a_{2},a_{0}}$, and in the
following we show results for two qualitatively different cases:

$\displaystyle V^{(1)}_{u\bar{u},g}(z)$
$\displaystyle=T_{F}\bigl{[}z^{2}+(1-z)^{2}\bigr{]}\,,$ $\displaystyle
V^{(1)}_{\delta g\delta g,g}(z)$ $\displaystyle=2C_{A}\,z(1-z)\,.$ (39)

Whereas the first kernel is finite for $z\to 0$ and $z\to 1$, the second one
goes to zero in both limits. This implies a strong decrease of
$\widehat{F}_{\delta g\delta g}^{\mskip 1.5mu\text{spl}}(x_{1},x_{2})$ for
both $x_{1}\ll x_{2}$ and $x_{1}\gg x_{2}$. Other splitting kernels diverge at
the end points — namely when an unpolarized gluon produced by the splitting
becomes soft — and we verified that the interpolation accuracy for these
channels is not worse than the one for $u\bar{u}$.

We see in the top row of figure -174 that already the medium grid yields a
very good interpolation as long as the momentum fractions are not too close to
the kinematic limit. At small $x_{1}$ and $x_{2}$, the accuracy is somewhat
less good than for the product form in figure -176, which can be understood
because the splitting form has more structure in that region. This is also
seen for the interpolation of $\widehat{F}_{\delta g\delta g}^{\mskip
1.5mu\text{spl}}$ at small momentum fractions in panel (c), where the accuracy
of the medium grid is at the permille level. In panel (d) we see that the
linear decrease of the distribution for $x_{1}\ll x_{2}$ requires the fine
grid for a satisfactory interpolation at the smallest $x_{1}$, where the DPD
is many orders of magnitude smaller than its maximum value. On the other hand,
the medium grid still performs well for $x_{1}>10^{-4}$. We note that the
coarse grid does not give a good accuracy in any of the cases considered here,
and it is entirely unreliable in the case of panel (d). We conclude that at at
least 16 points per subgrid should be taken, unless one considers only
unpolarized DPDs and has only very moderate accuracy requirements.

###  4.3 Estimating the interpolating accuracy

Let us finally study the reliability of the accuracy estimate for
interpolation described in subsection 3.1, given by

relative accuracy estimate $\displaystyle=\biggl{|}\,\frac{\text{value
interpolated without end points}}{\text{value interpolated with end
points}}-1\,\biggr{|}\,.$ (40)

(a) $\widehat{F}_{u\bar{u}}^{\mskip 1.5mu\text{int}\mskip 1.5mu(0)}$ at
$x_{2}=0.3$ (b) $\widehat{F}_{u\bar{u}}^{\mskip 1.5mu\text{int}\mskip
1.5mu(1)}$ at $x_{2}=0.3$

(c) $\widehat{F}_{u\bar{u}}^{\mskip 1.5mu\text{int}\mskip 1.5mu(2)}$ at
$x_{2}=0.3$ (d) $\widehat{F}_{u\bar{u}}^{\mskip 1.5mu\text{spl}}$ at
$x_{2}=10^{-4}$

Figure -173: True interpolation accuracy in $x_{1}$ and $x_{2}$ compared with
its estimate (40). All values refer to the medium $x$ grid in eq. (37). The
DPDs used are specified in eqs. (35) and (36).

In figure -173, we compare this estimate with the true relative accuracy (38)
for the medium $x$ grid in eq. (37). We see that the estimate is typically not
more than two orders of magnitude away from the true accuracy (except around
points where the error estimate or the true error has a zero crossing).
Somewhat surprisingly, the estimate almost coincides with the true accuracy in
some cases, as seen in figure -173(a). Closer inspection reveals that the
estimate tends to be above the true accuracy close to subinterval end points;
this is expected because in that region interpolation without the end points
is particularly bad by construction. We find a qualitatively similar picture
for the accuracy estimate in other kinematic settings, and for the other two
grids in eq. (37). Although not perfectly accurate in all regions, the method
based on eq. (40) provides a good and easy-to-compute indicator for the
quality of interpolation.

##  5 Evolution and flavor matching of DPDs

In this section, we demonstrate the accuracy of DGLAP evolution and flavor
matching for DPDs with the methods presented in subsections 3.3 and 3.4. As a
default setting for DPD evolution, we take a maximal step size of $h=0.22$ for
the evolution time (29) in the Runge-Kutta method. We use the $8^{\text{th}}$
order algorithm of Dormand and Prince [78], where each Runge-Kutta step is
divided into 13 substeps for which the r.h.s. of the evolution equation is
evaluated.

As initial conditions we take the splitting form (7) evaluated with the Les
Houches benchmark PDFs of refs. [79, 80]. Following the settings in these
references, we take

$\displaystyle\mu_{0}$ $\displaystyle=\sqrt{2}\,\mathrm{GeV}$ (41)

as initial scale with $n=3$ active flavors for both partons. The evolution
accuracy is independent of $y$, and we fix this variable to the value
$y=b_{0}/\mu_{0}\approx 0.794\,\mathrm{GeV}^{-1}$. Unless specified otherwise,
flavor matching is carried out at the canonical scales $\mu_{q}=m_{q}$, with
the quark masses $m_{c}=\sqrt{2}\,\mathrm{GeV}$, $m_{b}=4.5\,\mathrm{GeV}$,
and $m_{t}=175\,\mathrm{GeV}$. The strong coupling at the initial scale is
taken to be $\alpha_{s}(\mu_{0})=0.35$. Throughout this section, we consider
unpolarized partons and use evolution and flavor matching at NNLO. We specify
the scales and the number of active flavors for the two partons in the form
$(\mu_{1},\mu_{2})$ and $(n_{1},n_{2})$.

Guided by the studies in section 4, we interpolate both $x_{1}$ and $x_{2}$ on
the grid

$\displaystyle[10^{-5},5\times 10^{-3},0.5,1]_{(16,16,24)}\,.$ (42)

This corresponds to the medium grid in eq. (37), except for the third
subinterval, where we afford a finer grid for improved interpolation accuracy
at large momentum fractions.

###  5.1 Path independence of evolution and flavor matching

For a given heavy quark flavor $q$, we apply the matching equations (2.1) at a
fixed scale $\mu_{q}$, so that in the $(\mu_{1},\mu_{2})$ plane we match at
the line $(\mu_{q},\mu_{2})$ for the first parton and at the line
$(\mu_{1},\mu_{q})$ for the second one. We show in appendix A that, under
these conditions, the result of evolution and flavor matching between two
points in the $(\mu_{1},\mu_{2})$ plane is independent of the evolution path,
including the order in which the different flavor matching steps are carried
out.

Figure -172: Two paths for evolution and flavor matching in the
$(\mu_{1},\mu_{2})$ plane. Next to each line segment, we specify the numbers
$(n_{1},n_{2})$ of active flavors for which the DPD is evolved. Flavor
matching is performed at the points marked by dots. The path in black is the
ChiliPDF default for matching from $(3,3)$ to $(6,6)$ and for subsequent
evolution up to the final scale pair.

We verify this numerically by comparing the results of evolving and matching
from $(1\,\mathrm{GeV},1\,\mathrm{GeV})$ with $(3,3)$ flavors to
$(500\,\mathrm{GeV},1000\,\mathrm{GeV})$ with $(6,6)$ flavors along the two
paths shown in figure -172. The path in black is the default taken in ChiliPDF
when matching from 3 to 6 flavors for both partons and then evolving to the
final pair of scales. The staircase-like path in red can be selected by
subsequent evolution calls for each line segment. Each point marked by a dot
in the figure corresponds to the initial condition of the DPD with given
flavor numbers $(n_{1},n_{2})$ in ChiliPDF. For any scale pair
$(\mu_{1},\mu_{2})$, the DPD with these flavor numbers is computed by
evolution from this initial condition.

The results of evolving and matching along the two specified paths are in
excellent agreement, as illustrated in figure -171 for the $t\bar{t}$
distribution, for which the differences are largest among all flavor
combinations. The agreement deteriorates close to the kinematic limit
$x_{1}=1-x_{2}$ but remains small compared with the interpolation accuracy in
that region. The choice for a particular evolution path in our implementation
is therefore of no consequence for the accuracy.

(a) $t\bar{t},~{}x_{2}=10^{-4}$ (b) $t\bar{t},~{}x_{2}=0.3$

Figure -171: Relative difference between the $t\bar{t}$ distribution obtained
by evolving and matching along one of the two paths in figure -172. Both paths
lead from $(1\,\mathrm{GeV},1\,\mathrm{GeV})$ with $(3,3)$ flavors to
$(500\,\mathrm{GeV},1000\,\mathrm{GeV})$ with $(6,6)$ flavors. The
distributions at $(1\,\mathrm{GeV},1\,\mathrm{GeV})$ are obtained by evolving
backwards from the starting conditions at $(\mu_{0},\mu_{0})$ specified at the
beginning of the section.

###  5.2 Evolution and flavor matching accuracy

We now study evolution and flavor matching to higher scales, using the initial
conditions specified at the beginning of this section. To assess the accuracy,
we use three settings:

  1. 1.

our default setting with the grid $[10^{-5},5\times
10^{-3},0.5,1]_{(16,16,24)}$ and a maximal step size $h=0.22$ for the Runge-
Kutta algorithm,

  2. 2.

a finer grid $[10^{-5},5\times 10^{-3},0.5,1]_{(40,40,40)}$ with the same
maximal step size $h=0.22$,

  3. 3.

the same finer grid with a smaller maximal step size $h=0.0088$.

The difference between settings 1 and 2 is taken as an estimate of the error
from discretization in $x_{1}$ and $x_{2}$ with our default grid, and the
difference between settings 2 and 3 is taken as an estimate of the error due
to the Runge-Kutta algorithm.

These errors are shown in figure -170 for evolution and matching to
$(100\,\mathrm{GeV},10\,\mathrm{TeV})$ with $(5,6)$ active flavors. The plots
are for selected parton combinations; we verified that the errors for the
other combinations are of similar size or even smaller. We see that the Runge-
Kutta errors are always negligible compared with the errors from
discretization. The latter are of similar size as (or even smaller than) the
interpolation errors for the initial conditions, which are shown for
$u\bar{u}$ in figure -174. Evolution to high scales does therefore _not_
degrade the overall accuracy of DPDs with our method. Only close to the
kinematic limit does the relative error from discretization become large. With
the chosen grid, it remains below 1% for $1-(x_{1}+x_{2})$ above $0.2$.

(a) $u\bar{u},~{}x_{2}=10^{-4}$ (b) $dg,~{}x_{2}=10^{-4}$

(c) $d\bar{s},~{}x_{2}=0.3$ (d) $g\bar{u},~{}x_{2}=0.3$

Figure -170: Relative Runge-Kutta and discretization errors for DPD evolution
and flavor matching from scales $(\mu_{0},\mu_{0})$ with $(3,3)$ flavors to
$(100\,\mathrm{GeV},10\,\mathrm{TeV})$ with $(5,6)$ flavors. The initial
conditions are specified at the beginning of the section.

In figure -169, we show the errors for evolution and matching to $(1.01\mskip
1.5mu\mu_{c},1.01\mskip 1.5mu\mu_{b})$ with $(4,5)$ flavors, i.e. to scales
just slightly above the matching scales for the final flavor numbers. The
accuracy remains excellent, even for DPDs with charm as first or bottom as
second parton, which are still tiny at the chosen scales.

(a) $\bar{c}b,~{}x_{2}=10^{-4}$ (b) $s\bar{s},~{}x_{2}=10^{-4}$

(c) $gg,~{}x_{2}=0.3$ (d) $\bar{u}b,~{}x_{2}=0.3$

Figure -169: As figure -170, but for different flavor combinations evolved and
matched to scales $(1.01\mskip 1.5mu\mu_{c},1.01\mskip 1.5mu\mu_{b})$ with
$(4,5)$ flavors. A dashed black line indicates that the interpolated function
is negative.

###  5.3 Backward evolution accuracy

A stringent way to test the accuracy of an evolution algorithm is backward
evolution from large to small scales. Whilst small inaccuracies of
distributions at a given scale are diminished when evolving to higher scales,
they are amplified by backward evolution. In the following, we adapt the
backward evolution test performed for PDFs in ref. [55] to the case of DPDs.

We start with the three-flavor input DPD specified at the beginning of this
section, match to $4$ flavors at the input scales $(\mu_{0},\mu_{0})$, and
then evolve to $(\mu_{\text{low}},\mu_{\text{low}})$ with
$\mu_{\text{low}}=m_{b}/2=2.25\,\mathrm{GeV}$. At this point, we match from
$4$ to $5$ flavors for both partons. Since the matching scale differs from
$m_{b}$, all DPDs are nonzero at that point, including the ones with bottom
quarks or antiquarks. We take these distributions as starting conditions for
evolution with $5$ flavors, evolving

  1. 1.

from $(\mu_{\text{low}},\mu_{\text{low}})$ to
$(\mu_{\text{high}},\mu_{\text{high}})$ with
$\mu_{\text{high}}=1\,\mathrm{TeV}$ in step 1,

  2. 2.

from $(\mu_{\text{high}},\mu_{\text{high}})$ back to
$(\mu_{\text{low}},\mu_{\text{low}})$ in step 2.

We use our default value of $h=0.22$ for the maximal Runge-Kutta step size.

The relative difference between the input to the first step and the output of
the second step is sensitive to the accuracy of forward _and_ backward
evolution. In figure -168, this quantity is shown for selected parton
combinations; the accuracy is not worse for other combinations. We find that
backward evolution is very precise, with the relative errors in our exercise
being below $10^{-6}$ over most of the phase space. This also holds for
distributions containing bottom quarks or antiquarks, which have zero
crossings at the scale $\mu_{\text{low}}$ where the ratio shown in the figure
is evaluated.

(a) $x_{2}=10^{-4}$ (b) $x_{2}=0.3$

Figure -168: Relative accuracy of evolving DPDs from
$\mu_{\text{low}}=m_{b}/2$ to $\mu_{\text{high}}=1\,\mathrm{TeV}$ and back to
$\mu_{\text{low}}$ for both partons. Details are given in the text.

We note that backward evolution is inherently more demanding for two scales
than for a single scale. Indeed, step 2 of our study corresponds to backward
evolution first from $(\mu_{\text{high}},\mu_{\text{high}})$ to
$(\mu_{\text{low}},\mu_{\text{high}})$ and then from
$(\mu_{\text{low}},\mu_{\text{high}})$ to
$(\mu_{\text{low}},\mu_{\text{low}})$. Inaccuracies accumulated in the first
substep are amplified in the second one. We repeated the above exercise for
evolution in only one scale, evolving from
$(\mu_{\text{low}},\mu_{\text{low}})$ up to
$(\mu_{\text{high}},\mu_{\text{low}})$ and then down to
$(\mu_{\text{low}},\mu_{\text{low}})$. We find that indeed the accuracy is
even higher in that case, by one or two orders of magnitude in some kinematic
regions.

##  6 Interpolation of DPDs in the distance between partons

In the discussion so far, the transverse distance $y$ between partons did not
play an active role, given that evolution and flavor matching of DPDs are
performed for each discretized value of $y$ individually, without any cross
talk between different values of $y$. Once a DPD is computed for selected
scales $(\mu_{1},\mu_{2})$ and flavor numbers $(n_{1},n_{2})$, two different
tasks are concerned with the $y$ dependence:

  1. 1.

interpolation of the DPD for specific values of $y$, in order to study how
evolution and matching have modified the $y$ dependence of the initial
conditions. In the study described in ref. [81], visible effects of evolution
on the large-$y$ dependence of DPDs were reported, and in ref. [24] it was
found that evolution can substantially flatten the $y^{-2}$ behavior generated
by the perturbative splitting mechanism at small distances. Numerically
reliable interpolation of DPDs over a wide range of $y$ is hence important for
understanding the dynamics of evolution.

  2. 2.

integration over $y$ from a minimum value $y_{\min}$ up to infinity. First and
foremost, this is needed for computing the double parton luminosities (2.1)
that enter in cross sections of double parton scattering processes. It is also
needed for evaluating sum rules (see section 7), which DPDs must satisfy and
which provide valuable constraints for their modeling. Note that one
integrates the product of two DPDs when computing double parton luminosities
and single DPDs when evaluating the sum rules.

Both tasks can be carried out using Chebyshev grids, as described in
subsection 3.1. To this end, one first needs to find variable transformations
that map the characteristic $y$ dependence of DPDs to a dependence in the
variable $u$ that is well suited for interpolation by high-order polynomials.
In the next two subsections, we show how this can be achieved. After that, we
will demonstrate the resulting accuracy of interpolation and integration.

###  6.1 Characteristic $y$ dependence and associated variable transformations

It is natural to distinguish the region of small $y$, where perturbative
mechanisms are at work from the domain of large $y$.

As mentioned earlier and discussed at length in ref. [24], DPDs at small $y$
can be represented as the superposition of a part due to perturbative
splitting, with a power behavior like $y^{-2}$, and an “intrinsic” part that
behaves like $y^{0}$ and is related to twist-four distributions. Both power
laws are modified by logarithms from loop corrections. The relative weight of
the two parts strongly depends on the parton combination and on the momentum
fractions $x_{1}$ and $x_{2}$ in the DPD. As already mentioned, the $y^{-2}$
behavior can be substantially flattened by evolution to high scales. In the
small-$y$ region, the transformation $u(y)$ should therefore be suited to
handle the interpolation of a power behavior $y^{-q}$ with $0\leq q\leq 2$.

A complication arises from the fact that for certain parton combinations, such
as $uu$ or $u\bar{d}$, the leading-order splitting contribution (7) is zero.
Quark-gluon mixing results in nonzero values for these combinations when
evolving from the scale $\mu_{y}\sim 1/y$ at which the splitting formula is
evaluated to the final scale $\mu$. Hence, the splitting part of the DPD at
scale $\mu$ has a zero around $y\sim 1/\mu$. The position of this zero is
shifted when the intrinsic part is added.

At large distances, one expects a gradual decrease of DPDs governed by some
characteristic nonperturbative mass scale. This is often modeled by a
Gaussian, $F\sim\exp(-M^{2}y^{2})$, but one may also argue in favor of a
Yukawa-like decrease, $F\sim\exp(-My)$. The characteristic mass scale $M$ may
depend on the parton combination and on $x_{1}$ and $x_{2}$ [81, 41], and it
can slowly change under evolution.

We now discuss variable transformations that we found to work well for the
types of $y$ dependence just discussed. An overview of the transformations and
their properties is given in tables 2 and 2. By convention, all
transformations are such that $u(y)$ has a positive derivative in the domain
where it is used. In the case of interpolation, the function $f(y)$ in table 2
represents a DPD at fixed momentum fractions or scales. For sum rule integrals
$f(y)$ is $2\pi y$ times a DPD (see eq. (55)), and for double parton
luminosities $f(y)$ is $2\pi y$ times the product of two DPDs (see eq. (52)).
The overall normalization of $f(y)$ is irrelevant for the interpolation or
integration accuracy and therefore set to an arbitrary value here.

no. | ${}-u(y)$ | $y(u)$ | $y(-1)$ | $\mathrm{d}y/\mathrm{d}u$
---|---|---|---|---
1 | $y^{-\alpha}$ | $|u|^{-1/\alpha}$ | $1$ | $\frac{1}{\alpha}\,|u|^{\mskip 1.5mu-(1+\alpha)/\alpha}$
2 | $\exp\Bigl{[}-\frac{1}{4}\mskip 1.5mu(m^{2}y^{2}+my)\Bigr{]}$ | $\frac{1}{2m}\,\bigl{(}\sqrt{16L(u)+1}-1\bigr{)}$ | $0$ | $\frac{4}{m}\,|u|^{-1}\,\bigl{(}16L(u)+1\bigr{)}^{-1/2}$
3 | $\exp\Bigl{[}-\frac{1}{4}\mskip 1.5mumy\Bigr{]}$ | $\frac{4}{m}\mskip 1.5muL(u)$ | $0$ | $\frac{4}{m}\mskip 1.5mu|u|^{-1}$
4 | $\exp\Bigl{[}1-\sqrt{1+my/2}\Bigr{]}$ | $\frac{2}{m}\,\bigl{(}L^{2}(u)+2L(u)\bigr{)}$ | $0$ | $\frac{4}{m}\,|u|^{-1}\,\bigl{(}L(u)+1\bigr{)}$

Table 1: Variable transformations suitable for interpolation in $y$. The
parameters $\alpha$ and $m$ must be positive. All transformation satisfy
$u\leq 0$, $y(0)=\infty$, and $\mathrm{d}y/\mathrm{d}u>0$ for $u\in[-1,0]$. We
abbreviate $L(u)=\ln\bigl{(}1/|u|\bigr{)}$.

no. | $f(y)$ | $f\bigl{(}y(u)\bigr{)}$ |  $\int^{\infty}\mathrm{d}y\,f(y)$ requires
---|---|---|---
1 | $y^{-\beta}$ | $|u|^{\mskip 1.5mu\beta/\alpha}$ | $\alpha<\beta-1$
2 | $\exp(-M^{2}y^{2})$ | $|u|^{\mskip 1.5mu(2M/m)^{2}}\,\exp\Bigl{[}\frac{M^{2}}{m}\mskip 1.5muy(u)\Bigr{]}$ | $m<2M$
3 | $\exp(-M^{2}y^{2})$ | $|u|^{\mskip 1.5muL(u)\;(4M/m)^{2}}$ | no condition
3 | $\exp(-My)$ | $|u|^{\mskip 1.5mu4M/m}$ | $m<4M$
4 | $\exp(-My)$ | $|u|^{\mskip 1.5muL(u)\;2M/m}\;|u|^{\mskip 1.5mu4M/m}$ | no condition

Table 2: Behavior of the variable transformations in table 2 for selected
functions $f(y)$. In the third column we used $\exp\bigl{[}-c\mskip
1.5muL(u)\bigr{]}=|u|^{c}$. The condition in the last column ensures that
$\mathrm{d}y/\mathrm{d}u\,f\bigl{(}y(u)\bigr{)}\to 0$ for $u\to 0$. It is
understood that the parameters $M$ and $m$ are always positive.

The last three transformations in the table are designed for functions with a
decrease like $\exp(-M^{2}y^{2})$ or $\exp(-My)$ at large $y$, possibly
modified by a power law in $y$. They map the physical interval
$y\in[0,\infty]$ onto $u\in[-1,0]$ and depend on a mass parameter $m$, which
is normalized such that $\mathrm{d}y/\mathrm{d}u=4/m$ at $y=0$. The integral

$\displaystyle\int_{y_{\min}}^{\infty}\\!\\!\mathrm{d}y\;f(y)$
$\displaystyle=\int_{u_{\min}}^{0}\\!\\!\mathrm{d}u\;\frac{\mathrm{d}y}{\mathrm{d}u}\,f\bigl{(}y(u)\bigr{)}$
(43)

is computed with the Clenshaw-Curtis quadrature rule, which includes the
interval boundary $u=0$ where the Jacobian $\mathrm{d}y/\mathrm{d}u$ is
infinite. The upper limit on $m$ specified in the last column of table 2
ensures that the product $\mathrm{d}y/\mathrm{d}u\,f\bigl{(}y(u)\bigr{)}$ is
zero at the point $u=0$, so that this point does not contribute to the
Clenshaw-Curtis rule. We note that the value of the mass scale $M$ in the
integrand depends on whether one integrates a product of two DPDs or a single
DPD.

Inverse power law:

Transformation no. 1 in the table is $u(y)=-y^{-\alpha}$. We will show that
with $\alpha\sim 0.2$ this is well suited for the range of power laws
$y^{-\beta}$ expected for DPDs at small distances. The more general form
$u(y)=-(y+a)^{-\alpha}$ can be used down to $y=0$ if $a>0$, but we will not
need this here.

Gaussian:

Transformation no. 2 in the table has a Gaussian falloff at large $y$ and is
designed for functions with a Gaussian behavior, $f(y)\sim\exp(-M^{2}y^{2})$.
The transformed function $f\bigl{(}y(u)\bigr{)}$ approximately behaves like a
power of $|u|$ for $u\to 0$, and we find that good interpolation and
integration accuracy is obtained for $m\sim M$.

Exponential:

Transformation no. 3 in the table is an exponential in $y$. It is suitable for
both $f(y)\sim\exp(-My)$ and $f(y)\sim\exp(-M^{2}y^{2})$. In the first case,
the transformed function $f\bigl{(}y(u)\bigr{)}$ behaves like a power of $|u|$
for $u\to 0$, whereas in the second case it falls off faster than any power of
$|u|$. In both cases, values $m\sim M$ yield good interpolation and
integration accuracy.

Exponential with square root:

Transformation no. 4 in the table is designed for functions
$f(y)\sim\exp(-My)$ with a Yukawa-type decrease. The transformed function
$f\bigl{(}y(u)\bigr{)}$ decreases faster than any power of $|u|$ for $u\to 0$,
and good accuracy is obtained with $m\sim M$.

When handling DPDs with different flavor combinations and in different
kinematic regions, the mass scale $M$ of $f(y)$ in the preceding discussion
does not have a unique value. $M$ should instead be understood as an average,
or in the case of the bounds in the last column of table 2 as the lower limit
of the relevant values.

###  6.2 Composite grids

Given the different behavior of DPDs at small and large $y$, it is natural to
use interpolation grids with different transformations $u(y)$. At small $y$,
there is an additional physics reason for using several subgrids. Not only the
scale $\mu_{y}\sim 1/y$ but also the number $n$ of active flavors for which
the splitting DPD is naturally initialized depends on $y$. Since the value of
a DPD in general changes with $n$, one should use a unique value of $n$ in a
subgrid, so as to avoid interpolating discontinuous functions. As discussed at
length in ref. [82], it is appropriate to evaluate the perturbative splitting
formula (7) with $n$ active flavors for $\gamma\mskip
1.5mum_{n}<\mu_{y}<\gamma\mskip 1.5mum_{n+1}$ with $\gamma\sim 1$. We
therefore place subgrid boundaries in $y$ at values around the inverse heavy
quark masses.

In the following studies, we will assume a $y$ dependence

$\displaystyle F(y;\kappa)$ $\displaystyle=y^{-2\kappa}\,\exp(-M^{2}y^{2})\,,$
(44)

for a single DPD, with $\kappa=0$ for the intrinsic and $\kappa=1$ for the
splitting part. We take $M=0.21\,\mathrm{GeV}$, which roughly corresponds to
the Gaussian parameter in quark-gluon DPDs of the model used in refs. [24,
30]. As for the function $f(y)$ in table 2, the overall normalization of
$F(y;\kappa)$ is irrelevant for the interpolation accuracy.

The smallest value of $y$ needed in a given physics setting corresponds to the
inverse of the largest hard scale that will be considered. In the following,
we will take $y_{\min}=1/(7\,\mathrm{TeV})$, which covers the production of
very heavy states at the LHC. For the interpolation of eq. (44) and for the
computation of the associated double parton luminosities, we use the following
combination of grids and transformations:

$\displaystyle\bigl{[}\mskip
1.5mu1/(7\,\mathrm{TeV}),1/m_{b},1/m_{c}\,\bigr{]}_{p_{1},\,p_{2}}$ inverse
power law transformation with $\alpha=0.2$, (45) $\displaystyle\bigl{[}\mskip
1.5mu1/m_{c},\infty\,\bigr{]}_{p_{3}}$ Gaussian transformation with $m=2M$,
(46)

where the “inverse power law” and the “Gaussian” transformation are
respectively given in entries no. 1 and 2 of table 2. We find that taking
$m=2M$ results in a somewhat better overall accuracy than the rule-of-thumb
value $m\sim M$ given in the previous subsection. We will consider three
different settings for the number of points in each subgrid:

$\displaystyle(p_{1},p_{2},p_{3})$
$\displaystyle=\begin{cases}(12,6,12)&\text{(coarse $y$ grid)}\\\
(16,8,16)&\text{(medium $y$ grid)}\\\ (24,12,24)&\text{(fine $y$
grid)}\end{cases}$ (47)

Notice the relatively small number of points in the central grid, which covers
the narrow $y$ region in which one would initialize the splitting part of the
DPD with $4$ active flavors. For the quark masses, we use the values specified
at the beginning of section 5, which results in subgrid boundaries at
$1/m_{b}\approx 0.22\,\mathrm{GeV}^{-1}$ and $1/m_{c}\approx
0.71\,\mathrm{GeV}^{-1}$.

###  6.3 Accuracy of interpolation in $y$

In the top panels of figure -167, we show the interpolation accuracy for the
form (44) with $\kappa=0$ and $\kappa=1$. In the first and second subgrid, a
very accurate interpolation is obtained even with the coarse grid setting. In
the third subgrid, the accuracy degrades rather quickly for large distances,
but for the medium grid it remains below 1% for $y$ up to about
$4\,\mathrm{GeV}^{-1}\approx 0.8\operatorname{fm}$. The fine grid yields an
acceptable accuracy at even larger distances.

(a) $\exp(-M^{2}y^{2})$ (b) $y^{-2}\exp(-M^{2}y^{2})$

(c) $F(y;y_{0})$ with $y_{0}=10^{-2}$ (d) $F(y;y_{0})$ with $y_{0}=10^{-3}$

Figure -167: Relative interpolation accuracy (38) on the grids specified in
eqs. (45), (46), and (47). The interpolated functions are $F(y;\kappa)$ in the
upper panels and $F(y;y_{0})$ in the lower ones, see eqs. (44) and (48). A
dashed black line indicates that the interpolated function is negative. Here
and in all following plots, the mass parameter in the functions is
$M=0.21\,\mathrm{GeV}$.

Before discussing this issue in more detail, we consider the interpolation of
a DPD with a zero crossing in $y$, which arises for certain parton
combinations as explained above. As an example for this case, we take the
simple form

$\displaystyle F(y;y_{0})$ $\displaystyle=(y-y_{0})\exp(-M^{2}y^{2})\,,$ (48)

with different values of $y_{0}$ and the same mass parameter
$M=0.21\,\mathrm{GeV}$ as before. The resulting accuracy is shown in the
bottom panels of figure -167. We find that the coarse grid performs poorly,
whereas with the medium grid the relative accuracy stays below 1% except in
the immediate vicinity of the zero crossing. Very high accuracy is obtained
with the fine grid.

We note that for the intermediate $y$ grid from $1/m_{b}$ to $1/m_{c}$, the
setting with $p_{2}=8$ grid points is sufficient to obtain a relative
interpolation accuracy of $10^{-6}$ or better in all cases considered here. In
each of the grids for smaller and at larger $y$, one may take either $16$ or
$24$ points, depending on how stringent the accuracy requirements are. This
gives $40$, $48$, or $56$ for the total number of points.444Here we count the
common end points of neighboring subgrids twice, since different subgrids may
contain DPDs with different flavor numbers.

(a) $\exp(-M^{2}y^{2})$ (b) $y^{-2}\exp(-M^{2}y^{2})$

Figure -166: Relative interpolation accuracy at large $y$ for the grid in eq.
(46) with $p_{3}=24$ points (fine grid) and for the composite grids specified
by eqs. (49), (50), and (51) (finer and finest grid). The interpolated
functions are the same as in figures -167(a) and -167(b), respectively.

With the grids considered so far, the interpolation accuracy at very large $y$
becomes poor mainly because there are too few grid points in that region. For
many purposes, the behavior of the distributions in the extreme tail region is
not of particular interest, and it contributes only little to integrals over
$y$. Nevertheless, we wish to show that with an appropriate choice of subgrids
and grid transformations, one can also accurately describe this region. To
this end, we replace the large-$y$ grid in eq. (46) with

$\displaystyle\bigl{[}1/m_{c},3\,\mathrm{GeV}^{-1},10\,\mathrm{GeV}^{-1}\bigr{]}_{p_{a},\,p_{b}}$
Gaussian transformation with $m=\sqrt{2}M$, (49)
$\displaystyle\bigl{[}10\,\mathrm{GeV}^{-1},\infty\bigr{]}_{p_{c}}$ Gaussian
transformation with $m=2M$, (50)

where we consider two settings for the number of points:

$\displaystyle(p_{a},p_{b},p_{c})$
$\displaystyle=\begin{cases}(12,12,8)&\text{(finer $y$ grid)}\\\
(16,16,8)&\text{(finest $y$ grid})\end{cases}$ (51)

In figure -166, the resulting interpolation accuracy is compared with the
previous setting of eq. (46) with $p_{3}=24$ points. We see that both the
finer and the finest grid extend the domain of very accurate interpolation up
to $y\approx 10\,\mathrm{GeV}^{-1}\approx 2\operatorname{fm}$, with only a
moderate increase in the number of grid points.

(a) $y^{-2}\exp(-M^{2}y^{2})$ (b) $F(y;y_{0})$ with $y_{0}=10^{-3}$

Figure -165: True and estimated relative interpolation accuracy in $y$, as
defined in eqs. (38) and (40), respectively. The curves are for the fine grid
setting in eq. (47).

In analogy to subsection 4.3, we can estimate the interpolation accuracy in
$y$ by comparing the results of interpolation with and without subinterval end
points. We find that this method works very well, as shown by the examples in
figure -165. The biggest discrepancies between true and estimated errors
appear near subgrid boundaries, and in the central grid, where the omission of
the end points decreases the polynomial order of interpolation from 11 to 9.

###  6.4 Accuracy of integration over $y$

Integration over $y$ from a lower cutoff $y_{\min}$ to infinity is needed to
evaluate the double parton luminosities (2.1) and DPD sum rules. If $y_{\min}$
coincides with the end point of a subinterval of the interpolation grid, one
can directly use Clenshaw-Curtis quadrature to compute the integral.
Otherwise, we perform a number of intermediate steps:

  1. 1.

identify the subinterval $[y_{a},y_{b}]$ that contains $y_{\min}$,

  2. 2.

define a reduced subgrid on the interval $[y_{\min},y_{b}]$ with the same
number of points and the same variable transformation as the original subgrid,

  3. 3.

compute the required DPDs on the reduced subgrid by interpolation on
$[y_{a},y_{b}]$.

The integral is then computed using Clenshaw-Curtis quadrature on the reduced
subgrid and on the subgrids for $y\geq y_{b}$.

Let us investigate the accuracy of this procedure for the double parton
luminosities

$\displaystyle\mathcal{L}(y_{\min},\kappa_{1}+\kappa_{2})$
$\displaystyle=2\pi\\!\int\limits_{y_{\min}}^{\infty}\\!\mathrm{d}y\;y\,F(y;\kappa_{1})\,F(y;\kappa_{2})=2\pi\\!\int\limits_{y_{\min}}^{\infty}\\!\frac{\mathrm{d}y\;y}{y^{\mskip
1.5mu2(\kappa_{1}+\kappa_{2})}}\,\exp(-2M^{2}y^{2})$ (52)

that correspond to the simplified shapes of DPDs in eq. (44). We consider the
cases $\kappa_{1}+\kappa_{2}=0,1,2$, which respectively correspond to taking
the intrinsic parts of both DPDs, the intrinsic part of one and the splitting
part of the other DPD, or the splitting parts of both DPDs. Although only the
sum of intrinsic and splitting parts enters in physical cross sections, it is
useful to investigate these different combinations separately, as was for
instance done in the studies in refs. [24, 82].

(a) $\exp(-2M^{2}y^{2})$ (b) $\exp(-2M^{2}y^{2})$

(c) $y^{-2}\exp(-2M^{2}y^{2})$ (d) $y^{-2}\exp(-2M^{2}y^{2})$

(e) $y^{-4}\exp(-2M^{2}y^{2})$ (f) $y^{-4}\exp(-2M^{2}y^{2})$

Figure -164: Left: Relative accuracy of computing the integrals (52) with the
interpolation grids specified by eqs. (45), (46), and (47). Right: True and
estimated integration accuracy for the fine grid setting of eq. (47). The
functions in the legends are the product $F(y;\kappa_{1})\,F(y;\kappa_{2})$ of
DPDs in the luminosity integral.

The integration accuracy can be computed using the exact expression of (52) in
terms of an incomplete gamma function. The result is shown in the left panels
of figure -164 for the interpolation grids introduced in subsection 6.2. We
see that for $y_{\min}$ below $1\,\mathrm{GeV}^{-1}$, good accuracy is
achieved even with the coarse grid. For cutoff values above
$1\,\mathrm{GeV}^{-1}$, the accuracy quickly becomes unsatisfactory unless the
fine grid is used. This reflects the degrading interpolation accuracy in that
region, which was shown in the upper panels of figure -167. In physics
applications, $y_{\min}$ is of the order of an inverse hard scale, so that
large values of the cutoff will only be needed in special circumstances, for
instance if one is interested in the partial contribution of large $y$ to an
integral. We see that the fine grid setting is useful in such a case.

In the right panels of figure -164, we compare the true integration accuracy
with its estimate obtained from comparing the values obtained with Clenshaw-
Curtis quadrature and with Fejér’s second rule (see subsection 3.1). We see
that for $y_{\min}$ below $1\,\mathrm{GeV}^{-1}$, one obtains an overestimate
of the true error by one or two orders of magnitude. This is qualitatively
similar to what we found for integrals of PDFs over $x$ in section 3.4 of ref.
[55].

##  7 Cross check: comparison with DOVE

As a cross check for our implementation of DPD evolution in ChiliPDF, we
performed a comparison with DOVE, a code that was introduced in ref. [28] and
further developed or adapted in subsequent work in refs. [81, 24, 44, 45,
30].555Results obtained with an early version of DOVE are given on the website
https://gsdpdf.hepforge.org. The comparison is done with the version that was
used in ref. [30] to compute the Gaunt-Stirling sum rules for DPDs. In the
following, we show the momentum sum rule for a $u$ quark and the number sum
rule for the parton combination $gd_{v}$, which read

$\displaystyle\sum_{a_{2}}\int\limits_{0}^{1-x_{1}}\mathrm{d}x_{2}\;x_{2}\,F_{u\mskip
1.5mua_{2}}^{\scalebox{0.6}{$\overline{\text{MS}}$}}(x_{1},x_{2};\mu)$
$\displaystyle=(1-x_{1})\,f_{u}(x_{1};\mu)\,,$ (53)
$\displaystyle\int\limits_{0}^{1-x_{1}}\mathrm{d}x_{2}\;F_{g\mskip
1.5mud_{v}}^{\scalebox{0.6}{$\overline{\text{MS}}$}}(x_{1},x_{2};\mu)$
$\displaystyle=f_{g}(x_{1};\mu)\,,$ (54)

where $F_{g\mskip 1.5mu\smash{d_{v}}}=F_{g\smash{d}}-F_{g\smash{\bar{d}}}$.
The expressions for general parton combinations are given in equations (3.1)
and (3.2) of ref. [28]. As shown in ref. [83], the sum rules hold for the
distributions

$\displaystyle
F_{a_{1}a_{2}}^{\scalebox{0.6}{$\overline{\text{MS}}$}}(x_{1},x_{2};\mu)$
$\displaystyle=2\pi\\!\\!\int\limits_{b_{0}/\nu}^{\infty}\\!\mathrm{d}y\;y\,F_{a_{1}a_{2}}(x_{1},x_{2},y;\mu,\mu)+F_{a_{1}a_{2}}^{\text{match}}(x_{1},x_{2};\nu,\mu)\,,$
(55)

which are integrated over $y$ with an $\overline{\text{MS}}$ renormalization
prescription for the $y^{-2}$ singularity due to perturbative splitting. The
matching term $F_{a_{1}a_{2}}^{\text{match}}(x_{1},x_{2};\nu,\mu)$ relates
this prescription with the cutoff in $y$ used on the r.h.s. of eq. (55). It
can be computed in terms of a PDF and a matching kernel in a similar fashion
as the splitting part in eq. (7), and its LO expression is given in equation
(18) of ref. [30]. The dependence on the matching scale $\nu$ cancels between
the two terms in (55) up to power corrections in $1/\nu^{2}$. We recall that
$b_{0}=2e^{-\gamma_{E}}\approx 1.12$.

The different roles played by $x_{1}$ and $x_{2}$ in the sum rules (53) and
(54) motivate us to take different grids for the two momentum fractions,

$\displaystyle x_{1}:$ $\displaystyle\qquad\bigl{[}5\times 10^{-5},5\times
10^{-3},0.5,0.9,1.0\bigr{]}_{(12,12,16,8)}\,,$ $\displaystyle x_{2}:$
$\displaystyle\qquad\bigl{[}5\times 10^{-5},5\times
10^{-3},0.5,1.0\bigr{]}_{(12,16,16)}\,.$ (56)

The first grid has an additional small subgrid that leads to a better
interpolation accuracy at large $x_{1}$, which is not needed for integration
over $x_{2}$. The sum rule integrals over $x_{2}$ are computed using Clenshaw-
Curtis quadrature over the full $x_{2}$ grid. This implies that the integrals
are computed with a lower cutoff $x_{\min}=5\times 10^{-5}$. The same cutoff
is used in the computation with DOVE.

We find that the integration over $y$ in the sum rules can be computed with
the rather small grids

$\displaystyle\bigl{[}b_{0}/(170\,\mathrm{GeV}),b_{0}/(1\,\mathrm{GeV})\bigr{]}_{(12)}$
$\displaystyle\text{inverse power law transformation with $\alpha=0.2$}\,,$
(57) $\displaystyle\bigl{[}b_{0}/(1\,\mathrm{GeV}),\infty\bigr{]}_{(8)}$
$\displaystyle\text{Gaussian transformation with $m=M$}\,,$ (58)

where $M=1\big{/}\sqrt{4h_{gg}}\approx 0.23\,\mathrm{GeV}$ is the Gaussian
mass parameter for gluons in the initial conditions of the model (see eq.
(8)).

We take the DPD model specified in section 3 of ref. [30] and evolve to
$\mu_{1}=\mu_{2}=144.6\,\mathrm{GeV}$ with $n=3$ active flavors for each
parton. At that scale, the momentum and number sum rules are evaluated with
$\nu=144.6\,\mathrm{GeV}$. These $\mu$ and $\nu$ values correspond to grid
points in our setup of DOVE and thus avoid additional interpolation. The
detailed grid settings and integration procedure in DOVE are as specified in
section 4 of ref. [30].

(a) $u$ momentum sum rule  (b) $gd_{v}$ number sum rule

Figure -163: Relative differences between the left-hand sides of the DPD sum
rules (53) and (54), evaluated using DOVE and ChiliPDF. Details are given in
the text.

In figure -163, we show the relative deviation between the left-hand sides of
the sum rules (53) and (54) evaluated in ChiliPDF and in DOVE. The numerical
differences increase with $x_{1}$ and reach several percent at $x_{1}=0.9$.
For sum rules with other parton combinations, we find similar or smaller
discrepancies. We also compared the sum rules evaluated at smaller values of
$\mu_{1}=\mu_{2}$ and $\nu$ and found somewhat smaller discrepancies than at
high scales.

An error estimate for the accuracy of evolution with DOVE in its original
version is given in figure 9 of [28]. Given that this estimate was carried out
with finer grids than the ones used here, and given that relatively simple
integration rules are used in the evaluation of the sum rules with DOVE, we
find the level of agreement between the two codes satisfactory for the purpose
of cross-validation.

##  8 Conclusions

We have shown that discretization on Chebyshev grids allows for a precise and
economical representation of DPDs as functions of the momentum fractions
$x_{1}$, $x_{2}$ and the distance $y$ between the partons. This generalizes
the methods we developed for PDFs in ref. [55]. We use grids for $x_{1}$ that
are independent of $x_{2}$ and vice versa. This gives the highest degree of
simplicity and decoupling between operations that involve only one momentum
fraction, which is the case for the Mellin convolutions that appear in DGLAP
evolution and in flavor matching. This method incurs interpolation errors
close to the kinematic boundary $x_{1}+x_{2}=1$ of the DPDs, which we find
manageable. We investigated the interpolation accuracy for typical shapes of
DPDs, both for the input distributions and after evolution. We find that grids
with $p_{x}=54$ points for $x\geq 10^{-5}$ yield a very accurate
representation of DPDs, with relative errors that are much smaller than
$10^{-3}$ in large parts of the phase space and remain below $1\%$ for
$x_{1}+x_{2}<0.8$ (see figures -170 and -169).

For the dependence on $y$, Chebyshev interpolation is used after variable
transformations that are adapted to the functional form and the typical mass
scale for the decrease of the DPDs at large distance. We find that with
$p_{y}\sim 48$ grid points, very good accuracy can be achieved for
interpolation and integration.

To solve the DGLAP equations for DPDs, we use a high-order Runge-Kutta
algorithm for computing the evolution matrices in eqs. (33) and (34), which
are the discretized version of the Green functions for the evolution
equations. This provides very high numerical accuracy for both forward and
backward evolution, with errors well below the ones due to discretization in
$x_{1}$ and $x_{2}$. With our algorithm, it is natural to evolve DPDs
separately in the scales $\mu_{1}$ and $\mu_{2}$ associated with the two
partons. We note that the computational effort scales linearly with $p_{y}$
but with the third power of $p_{x}$.

In our setup, a DPD with given flavor numbers $(n_{1},n_{2})$ is specified
numerically by its values at all grid points in $(x_{1},x_{2},y)$ and at one
pair $(\mu_{1},\mu_{2})$ of scales. Its values at other scales are computed
“on the fly” by solving the evolution equations. DPDs at other values of
$x_{1}$, $x_{2}$, and $y$ are computed using the barycentric interpolation
formula, which is fast and numerically stable. As discussed in the
introduction, we avoid grids in the renormalization scales, which would
require a huge amount of computer memory. This is a notable difference between
our approach and the widely used method to compute PDFs from grids in $x$ and
$\mu$, for instance via the LHAPDF interface [54], and it reflects the notable
difference in complexity between the distribution functions for one or for two
partons inside a proton.

With our present implementation of the above methods in ChiliPDF, the
evolution of a DPD with $n=5$ active flavors (and thus 121 parton flavor
combinations) from $\mu_{1}=\mu_{2}=15\,\mathrm{GeV}$ to $150\,\mathrm{GeV}$
takes 1 to 2 seconds666Timings were obtained on different computers with
Intel® Core™ i5, Intel® Core™ i7, or AMD EPYC™ 7513 32-Core processors. for
$p_{x}=54$ grid points in $x_{1}$ and $x_{2}$ and $p_{y}=48$ grid points in
$y$. With grids of this size, the largest part of computing time is spent on
multiplying the evolution matrices with the initial conditions, so that there
is only a weak dependence of the timing on the perturbative order of evolution
and on the initial and final scales. Since the current code is based on a
simple implementation of linear algebra and matrix multiplication, we expect
that significant gains can still be made with a dedicated performance tuning.

We used our code to investigate the impact of higher orders in the DGLAP
kernels on the evolution of DPDs. For some parton combinations, the change
from LO to NLO evolution produces $\mathcal{O}(1)$ effects at small momentum
fractions, whereas we find only moderate differences when changing from NLO to
NNLO kernels.

The numerical delivery of realistic DPDs has been a bottleneck in practice so
far. The presented methods and their implementation provide a practical tool
for computing evolved DPDs from user-given starting conditions, on a par with
what has long been a standard for PDFs. This increases our ability to use the
predictive power of perturbation theory and thus presents an important
ingredient to making theory predictions for double parton scattering more
realistic.

The evolution equations considered in this work apply to DPDs that are summed
separately over the color of each parton. DPDs with color correlations between
the two partons follow a different evolution pattern, with kernels that have
recently been computed at NLO [84]. The methods described in the present paper
can be adapted to this case, as will be reported in a future publication.

## Acknowledgments

We thank Florian Fabry, Oskar Grocholski, and Mees van Kampen for their
contributions to ChiliPDF. It is a pleasure to thank Jonathan Gaunt for
valuable discussions and input.

This work is in part supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) – grant number 409651613 (Research Unit FOR 2926)
and grant number 491245950. The work of RN is supported by the ERC Starting
Grant REINVENT-714788.

##  Appendix A Path independence of evolution and flavor matching

In this appendix we show that the set (2.1) of DGLAP equations has a unique
solution for a given initial condition at a point in the $(\mu_{1},\mu_{2})$
plane. In other words, the solution is independent of the path in the
$(\mu_{1},\mu_{2})$ plane that is taken to evolve from the initial condition
to the final scales. To establish this, it is enough to show that the second
derivatives in the renormalisation scales commute, i.e. that

$\displaystyle\frac{\partial}{\partial L_{1}}\,\frac{\partial}{\partial
L_{2}}\,F^{n_{1},n_{2}}_{a_{1}a_{2}}(L_{1},L_{2})$
$\displaystyle=\frac{\partial}{\partial
L_{1}}\biggl{(}\sum_{b_{2}}P^{n_{2}}_{a_{2}b_{2}}(L_{2})\underset{2}{\otimes}F^{n_{1},n_{2}}_{a_{1}b_{2}}(L_{1},L_{2})\biggr{)}$
$\displaystyle=\sum_{b_{2},b_{1}}P^{n_{2}}_{a_{2}b_{2}}(L_{2})\underset{2}{\otimes}P^{n_{1}}_{a_{1}b_{1}}(L_{1})\underset{1}{\otimes}F^{n_{1},n_{2}}_{b_{1}b_{2}}(L_{1},L_{2})$
(59) is equal to $\displaystyle\frac{\partial}{\partial
L_{2}}\,\frac{\partial}{\partial
L_{1}}\,F^{n_{1},n_{2}}_{a_{1}a_{2}}(L_{1},L_{2})$
$\displaystyle=\frac{\partial}{\partial
L_{2}}\biggl{(}\sum_{b_{1}}P^{n_{1}}_{a_{1}b_{1}}(L_{1})\underset{1}{\otimes}F^{n_{1},n_{2}}_{b_{1}a_{2}}(L_{1},L_{2})\biggr{)}$
$\displaystyle=\sum_{b_{1},b_{2}}P^{n_{1}}_{a_{1}b_{1}}(L_{1})\underset{1}{\otimes}P^{n_{2}}_{a_{2}b_{2}}(L_{2})\underset{2}{\otimes}F^{n_{1},n_{2}}_{b_{1}b_{2}}(L_{1},L_{2})\,.$
(60)

Here and in the following, we change evolution variables from $\mu_{i}$ to

$\displaystyle L_{i}$ $\displaystyle=2\ln\mu_{i}$ (61)

and drop the arguments $x_{1},x_{2}$, and $y$ of the DPDs for brevity. For
evolution at fixed $(n_{1},n_{2})$, i.e. without flavor matching, the desired
equality follows from the fact that the convolution integrals w.r.t. the first
and second momentum fraction of the DPD commute with each other. This can
readily be shown from their definition (2.1). Note that it is essential for
the argument that the evolution kernel for one parton does not depend on the
scale of the other parton.

We now extend this argument, showing that one obtains the same result when
flavor matching along different paths in the $(n_{1},n_{2})$ plane. Let us see
how eqs. (A) and (A) have to be modified to include the effects from flavor
matching. We first consider the case of PDFs, where the notation is less
involved. For a set of given matching scales $\tilde{\mu}_{n}$ for the
transition from $n-1$ to $n$ active flavors, we define the “canonical” number
of flavors at given $L=2\ln\mu$ as

$\displaystyle n(L)$ $\displaystyle=\begin{cases}3&\text{ for
$L<\widetilde{L}_{4}$}\\\ 4&\text{ for
$\widetilde{L}_{4}<L<\widetilde{L}_{5}$}\\\ 5&\text{ for
$\widetilde{L}_{5}<L<\widetilde{L}_{6}$}\\\ 6&\text{ for
$L>\widetilde{L}_{6}$}\end{cases}$ (62)

where $\widetilde{L}_{n}=2\ln\tilde{\mu}_{n}$. We further define PDFs

$\displaystyle f_{a}(L)$ $\displaystyle=f_{a}^{n(L)}(L)\,,$ (63)

with the canonical number of flavors at a given renormalisation scale. From
these distributions, one can obtain $f^{n}_{a}(L)$ for flavor numbers
different from the canonical values by evolution.777This implies that our
argument is _not_ tied to the use of a specific “variable flavor number
scheme”.

At the matching scales, the distributions in eq. (63) can have
discontinuities, which complicates their description via differential
equations. We therefore introduce a smoothened version of $f_{a}(L)$. Taking a
single matching step from $n-1$ to $n$ flavors at
$\widetilde{L}=2\ln\tilde{\mu}$ for simplicity, one can define a smooth
approximation $\theta_{\epsilon}(L)$ of the step function that takes the
values

$\displaystyle\theta_{\epsilon}(t)$ $\displaystyle=\begin{cases}0&\text{ for
$t<-\epsilon$}\\\ 1&\text{ for $t>\epsilon$}\end{cases}$ (64)

and is continuous and monotonic in the intermediate region
$-\epsilon<t<\epsilon$. Its derivative
$\delta_{\epsilon}(t)=\mathrm{d}\theta_{\epsilon}(t)/\mathrm{d}t$ is then a
smooth approximation of the Dirac distribution. Define now a PDF
$f^{\epsilon}_{a}(L)$ that evolves according to

$\displaystyle\frac{\mathrm{d}f^{\epsilon}_{a}(L)}{\mathrm{d}L}$
$\displaystyle=\sum_{b}P_{ab}^{n_{\epsilon}(L)}(L)\underset{}{\otimes}f^{\epsilon}_{b}(L)$
$\displaystyle\quad+\delta_{\epsilon}(L-\widetilde{L})\biggl{[}\sum_{b}A_{ab}^{n}(m;\widetilde{L})\underset{}{\otimes}f^{\epsilon}_{b}(\widetilde{L}-\epsilon)-f^{\epsilon}_{a}(\widetilde{L}-\epsilon)\biggr{]}\,,$
(65)

with

$\displaystyle n_{\epsilon}(L)=n-1+\theta_{\epsilon}(L-\widetilde{L})$ (66)

interpolating smoothly between $n-1$ and $n$. The DGLAP kernels for the
evolution of $f^{\epsilon}_{a}(L)$ are thus evaluated with flavor number $n-1$
for $L\leq\widetilde{L}-\epsilon$, with flavor number $n$ for
$L\geq\widetilde{L}+\epsilon$, and with fractional flavor number in the
intermediate $L$ region.888PDFs and DGLAP kernels for inactive flavors are set
to zero, so that the range of flavor sums does not explicitly depend on $n$.
One furthermore has

$\displaystyle f^{\epsilon}_{a}(\widetilde{L}+\epsilon)$
$\displaystyle=f^{\epsilon}_{a}(\widetilde{L}-\epsilon)+\int_{\widetilde{L}-\epsilon}^{\widetilde{L}+\epsilon}\mathrm{d}L\,\frac{\mathrm{d}f^{\epsilon}_{a}(L)}{\mathrm{d}L}$
$\displaystyle=f^{\epsilon}_{a}(\widetilde{L}-\epsilon)+\mathcal{O}(\epsilon)+\int_{\widetilde{L}-\epsilon}^{\widetilde{L}+\epsilon}\mathrm{d}L\;\delta_{\epsilon}(L-\widetilde{L})$
$\displaystyle\qquad\times\biggl{[}\sum_{b}A_{ab}^{n}(m;\widetilde{L})\underset{}{\otimes}f^{\epsilon}_{b}(\widetilde{L}-\epsilon)-f^{\epsilon}_{a}(\widetilde{L}-\epsilon)\biggr{]}$
$\displaystyle=\sum_{b}A_{ab}^{n}(m;\widetilde{L})\underset{}{\otimes}f^{\epsilon}_{b}(\widetilde{L}-\epsilon)+\mathcal{O}(\epsilon)\,,$
(67)

where the $\mathcal{O}(\epsilon)$ contribution comes from the integral over
the term with $P^{n_{\epsilon}(L)}$ in eq. (A). In the limit $\epsilon\to 0$,
this tends to the flavor matching prescription given in eq. (5).

It is straightforward to extend this construction from one to several matching
steps and from PDFs to DPDs. In analogy to eq. (63), DPDs with canonical
flavor numbers are defined as

$\displaystyle F_{a_{1}a_{2}}(L_{1},L_{2})$
$\displaystyle=F_{a_{1}a_{2}}^{n(L_{1}),\,n(L_{2})}(L_{1},L_{2})\,,$ (68)

and their smoothened version $F^{\epsilon}_{a_{1}a_{2}}(L_{1},L_{2})$ evolves
according to

$\displaystyle\frac{\partial}{\partial
L_{1}}\,F^{\epsilon}_{a_{1}a_{2}}(L_{1},L_{2})$
$\displaystyle=\sum_{b_{1}}P^{n_{\epsilon}(L_{1})}_{a_{1}b_{1}}(L_{1})\underset{1}{\otimes}F^{\epsilon}_{b_{1}a_{2}}(L_{1},L_{2})$
$\displaystyle\quad+\sum_{n}\delta_{\epsilon}(L_{1}-\widetilde{L}_{n})\sum_{b_{1}}\Delta
A^{n}_{a_{1}b_{1}}\underset{1}{\otimes}F^{\epsilon}_{b_{1}a_{2}}(\widetilde{L}_{n}-\epsilon,L_{2})\,,$
$\displaystyle\frac{\partial}{\partial
L_{2}}\,F^{\epsilon}_{a_{1}a_{2}}(L_{1},L_{2})$
$\displaystyle=\sum_{b_{2}}P^{n_{\epsilon}(L_{2})}_{a_{2}b_{2}}(L_{2})\underset{2}{\otimes}F^{\epsilon}_{a_{1}b_{2}}(L_{1},L_{2})$
$\displaystyle\quad+\sum_{n}\delta_{\epsilon}(L_{2}-\widetilde{L}_{n})\sum_{b_{2}}\Delta
A^{n}_{a_{2}b_{2}}\underset{2}{\otimes}F^{\epsilon}_{a_{1}b_{2}}(L_{1},\widetilde{L}_{n}-\epsilon)\,,$
(69)

where we abbreviated

$\displaystyle\Delta A^{n}_{ab}(z)$
$\displaystyle=A^{n}_{ab}(z,m_{n};\widetilde{L}_{n})-\delta_{ab}\,\delta(1-z)\,.$
(70)

As at the beginning of this section, one can then show that

$\displaystyle\frac{\partial}{\partial L_{1}}\,\frac{\partial}{\partial
L_{2}}\,F^{\epsilon}_{a_{1}a_{2}}(L_{1},L_{2})$
$\displaystyle=\frac{\partial}{\partial L_{2}}\,\frac{\partial}{\partial
L_{1}}\,F^{\epsilon}_{a_{1}a_{2}}(L_{1},L_{2})\,,$ (71)

so that with an initial condition given in the $(L_{1},L_{2})$ plane, eq. (A)
has a unique solution independent of the evolution path. Taking the limit
$\epsilon\to 0$, the desired result is obtained.

## References

  * [1] CDF collaboration, F. Abe et al., _Double parton scattering in $\bar{p}p$ collisions at $\sqrt{s}$ = 1.8 TeV_, _Phys. Rev. D_ 56 (1997) 3811.
  * [2] D0 collaboration, V. M. Abazov et al., _Study of double parton interactions in diphoton + dijet events in $p\bar{p}$ collisions at $\sqrt{s}$ = 1.96 TeV_, _Phys. Rev. D_ 93 (2016) 052008 [1512.05291].
  * [3] LHCb collaboration, R. Aaij et al., _Measurement of the J/ $\psi$ pair production cross-section in pp collisions at $\sqrt{s}$ = 13 TeV_, _JHEP_ 06 (2017) 047 [1612.07451].
  * [4] ATLAS collaboration, _Measurement of J/ $\psi$ production in association with a W± boson with pp data at 8 TeV_, _JHEP_ 01 (2020) 095 [1909.13626].
  * [5] CMS collaboration, _Evidence for WW production from double-parton interactions in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$_, _Eur. Phys. J. C_ 80 (2020) 41 [1909.06265].
  * [6] CMS collaboration, _Measurement of double-parton scattering in inclusive production of four jets with low transverse momentum in proton-proton collisions at $\sqrt{s}$ = 13 TeV_, _JHEP_ 01 (2022) 177 [2109.13822].
  * [7] CMS collaboration, _Observation of same-sign WW production from double parton scattering in proton-proton collisions at $\sqrt{s}$ = 13 TeV_, 2206.02681.
  * [8] A. Kulesza and W. Stirling, _Like sign $W$ boson production at the LHC as a probe of double parton scattering_, _Phys. Lett. B_ 475 (2000) 168 [hep-ph/9912232].
  * [9] J. R. Gaunt, C.-H. Kom, A. Kulesza and W. Stirling, _Same-sign W pair production as a probe of double parton scattering at the LHC_ , _Eur. Phys. J. C_ 69 (2010) 53 [1003.3953].
  * [10] K. Golec-Biernat and E. Lewandowska, _Electroweak boson production in double parton scattering_ , _Phys. Rev. D_ 90 (2014) 094032 [1407.4038].
  * [11] F. A. Ceccopieri, M. Rinaldi and S. Scopetta, _Parton correlations in same-sign $W$ pair production via double parton scattering at the LHC_, _Phys. Rev. D_ 95 (2017) 114030 [1702.05363].
  * [12] S. Cotogno, T. Kasemets and M. Myska, _Spin on same-sign $W$-boson pair production_, _Phys. Rev. D_ 100 (2019) 011503 [1809.09024].
  * [13] S. Cotogno, T. Kasemets and M. Myska, _Confronting same-sign W-boson production with parton correlations_ , _JHEP_ 10 (2020) 214 [2003.03347].
  * [14] CMS collaboration, _Search for new physics in same-sign dilepton events in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$_, _Eur. Phys. J. C_ 76 (2016) 439 [1605.03171].
  * [15] CMS collaboration, _Search for top quark partners with charge 5/3 in the same-sign dilepton and single-lepton final states in proton-proton collisions at $\sqrt{s}=13$ TeV_, _JHEP_ 03 (2019) 082 [1810.03188].
  * [16] B. Blok, Y. Dokshitzer, L. Frankfurt and M. Strikman, _The Four jet production at LHC and Tevatron in QCD_ , _Phys. Rev. D_ 83 (2011) 071501 [1009.2714].
  * [17] J. R. Gaunt and W. Stirling, _Double Parton Scattering Singularity in One-Loop Integrals_ , _JHEP_ 06 (2011) 048 [1103.1888].
  * [18] M. G. Ryskin and A. M. Snigirev, _A Fresh look at double parton scattering_ , _Phys. Rev._ D83 (2011) 114047 [1103.3495].
  * [19] M. Diehl, D. Ostermeier and A. Schäfer, _Elements of a theory for multiparton interactions in QCD_ , _JHEP_ 03 (2012) 089 [1111.0910].
  * [20] A. V. Manohar and W. J. Waalewijn, _A QCD Analysis of Double Parton Scattering: Color Correlations, Interference Effects and Evolution_ , _Phys. Rev. D_ 85 (2012) 114009 [1202.3794].
  * [21] M. G. Ryskin and A. M. Snigirev, _Double parton scattering in double logarithm approximation of perturbative QCD_ , _Phys. Rev._ D86 (2012) 014018 [1203.2330].
  * [22] J. R. Gaunt, _Single Perturbative Splitting Diagrams in Double Parton Scattering_ , _JHEP_ 01 (2013) 042 [1207.0480].
  * [23] B. Blok, Y. Dokshitzer, L. Frankfurt and M. Strikman, _Perturbative QCD correlations in multi-parton collisions_ , _Eur. Phys. J. C_ 74 (2014) 2926 [1306.3763].
  * [24] M. Diehl, J. R. Gaunt and K. Schönwald, _Double hard scattering without double counting_ , _JHEP_ 06 (2017) 083 [1702.06486].
  * [25] O. Fedkevych and A. Kulesza, _Double parton scattering in four-jet production in proton-proton collisions at the LHC_ , _Phys. Rev. D_ 104 (2021) 054021 [2008.08347].
  * [26] P. Bartalini and J. R. Gaunt, eds., _Multiple Parton Interactions at the LHC_ , vol. 29. World Scientific Publishing, 2019, 10.1142/10646.
  * [27] T. Kasemets and S. Scopetta, _Parton correlations in double parton scattering_ , _Adv. Ser. Direct. High Energy Phys._ 29 (2018) 49 [1712.02884].
  * [28] J. R. Gaunt and W. J. Stirling, _Double Parton Distributions Incorporating Perturbative QCD Evolution and Momentum and Quark Number Sum Rules_ , _JHEP_ 03 (2010) 005 [0910.4347].
  * [29] K. Golec-Biernat, E. Lewandowska, M. Serino, Z. Snyder and A. M. Stasto, _Constraining the double gluon distribution by the single gluon distribution_ , _Phys. Lett. B_ 750 (2015) 559 [1507.08583].
  * [30] M. Diehl, J. R. Gaunt, D. M. Lang, P. Plößl and A. Schäfer, _Sum rule improved double parton distributions in position space_ , _Eur. Phys. J. C_ 80 (2020) 468 [2001.10428].
  * [31] K. Golec-Biernat and A. M. Staśto, _Momentum sum rule and factorization of double parton distributions_ , _Phys. Rev. D_ 107 (2023) 054020 [2212.02289].
  * [32] G. S. Bali, M. Diehl, B. Gläßle, A. Schäfer and C. Zimmermann, _Double parton distributions in the nucleon from lattice QCD_ , _JHEP_ 09 (2021) 106 [2106.03451].
  * [33] H.-M. Chang, A. V. Manohar and W. J. Waalewijn, _Double Parton Correlations in the Bag Model_ , _Phys. Rev. D_ 87 (2013) 034009 [1211.3132].
  * [34] M. Rinaldi, S. Scopetta and V. Vento, _Double parton correlations in constituent quark models_ , _Phys. Rev. D_ 87 (2013) 114021 [1302.6462].
  * [35] W. Broniowski and E. Ruiz Arriola, _Valence double parton distributions of the nucleon in a simple model_ , _Few Body Syst._ 55 (2014) 381 [1310.8419].
  * [36] M. Rinaldi, S. Scopetta, M. Traini and V. Vento, _Double parton correlations and constituent quark models: a Light Front approach to the valence sector_ , _JHEP_ 12 (2014) 028 [1409.1500].
  * [37] W. Broniowski, E. Ruiz Arriola and K. Golec-Biernat, _Generalized Valon Model for Double Parton Distributions_ , _Few Body Syst._ 57 (2016) 405 [1602.00254].
  * [38] T. Kasemets and A. Mukherjee, _Quark-gluon double parton distributions in the light-front dressed quark model_ , _Phys. Rev. D_ 94 (2016) 074029 [1606.05686].
  * [39] M. Rinaldi, S. Scopetta, M. C. Traini and V. Vento, _Correlations in Double Parton Distributions: Perturbative and Non-Perturbative effects_ , _JHEP_ 10 (2016) 063 [1608.02521].
  * [40] M. Rinaldi and F. A. Ceccopieri, _Relativistic effects in model calculations of double parton distribution function_ , _Phys. Rev. D_ 95 (2017) 034040 [1611.04793].
  * [41] R. Corke and T. Sjostrand, _Multiparton Interactions with an x-dependent Proton Size_ , _JHEP_ 05 (2011) 009 [1101.5953].
  * [42] B. Blok and P. Gunnellini, _Dynamical approach to MPI in W+dijet and Z+dijet production within the PYTHIA event generator_ , _Eur. Phys. J. C_ 76 (2016) 202 [1510.07436].
  * [43] B. Blok and P. Gunnellini, _Dynamical approach to MPI four-jet production in Pythia_ , _Eur. Phys. J. C_ 75 (2015) 282 [1503.08246].
  * [44] B. Cabouat, J. R. Gaunt and K. Ostrolenk, _A Monte-Carlo Simulation of Double Parton Scattering_ , _JHEP_ 11 (2019) 061 [1906.04669].
  * [45] B. Cabouat and J. R. Gaunt, _Combining single and double parton scatterings in a parton shower_ , _JHEP_ 10 (2020) 012 [2008.01442].
  * [46] A. Cafarella, C. Coriano and M. Guzzi, _Precision Studies of the NNLO DGLAP Evolution at the LHC with CANDIA_ , _Comput. Phys. Commun._ 179 (2008) 665 [0803.0462].
  * [47] G. P. Salam and J. Rojo, _A Higher Order Perturbative Parton Evolution Toolkit (HOPPET)_ , _Comput. Phys. Commun._ 180 (2009) 120 [0804.3755].
  * [48] M. Botje, _QCDNUM: Fast QCD Evolution and Convolution_ , _Comput. Phys. Commun._ 182 (2011) 490 [1005.1481].
  * [49] V. Bertone, S. Carrazza and J. Rojo, _APFEL: A PDF Evolution Library with QED corrections_ , _Comput. Phys. Commun._ 185 (2014) 1647 [1310.1394].
  * [50] V. Bertone, _APFEL++: A new PDF evolution library in C++_ , _PoS_ DIS2017 (2018) 201 [1708.00911].
  * [51] S. Weinzierl, _Fast evolution of parton distributions_ , _Comput. Phys. Commun._ 148 (2002) 314 [hep-ph/0203112].
  * [52] A. Vogt, _Efficient evolution of unpolarized and polarized parton distributions with QCD-PEGASUS_ , _Comput. Phys. Commun._ 170 (2005) 65 [hep-ph/0408244].
  * [53] A. Candido, F. Hekhorn and G. Magni, _EKO: evolution kernel operators_ , _Eur. Phys. J. C_ 82 (2022) 976 [2202.02338].
  * [54] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _Eur. Phys. J. C_ 75 (2015) 132 [1412.7420].
  * [55] M. Diehl, R. Nagar and F. J. Tackmann, _ChiliPDF: Chebyshev interpolation for parton distributions_ , _Eur. Phys. J. C_ 82 (2022) 257 [2112.09703].
  * [56] M. Diehl and T. Kasemets, _Positivity bounds on double parton distributions_ , _JHEP_ 05 (2013) 150 [1303.0842].
  * [57] R. Kirschner, _Generalized Lipatov-Altarelli-Parisi Equations and Jet Calculus Rules_ , _Phys. Lett. B_ 84 (1979) 266.
  * [58] V. P. Shelest, A. M. Snigirev and G. M. Zinovev, _The Multiparton Distribution Equations in QCD_ , _Phys. Lett. B_ 113 (1982) 325.
  * [59] A. M. Snigirev, _Double parton distributions in the leading logarithm approximation of perturbative QCD_ , _Phys. Rev. D_ 68 (2003) 114012 [hep-ph/0304172].
  * [60] F. A. Ceccopieri, _An update on the evolution of double parton distributions_ , _Phys. Lett. B_ 697 (2011) 482 [1011.6586].
  * [61] M. Buza, Y. Matiounine, J. Smith and W. L. van Neerven, _Charm electroproduction viewed in the variable flavor number scheme versus fixed order perturbation theory_ , _Eur. Phys. J. C_ 1 (1998) 301 [hep-ph/9612398].
  * [62] S. Moch, J. A. M. Vermaseren and A. Vogt, _The Three loop splitting functions in QCD: The Nonsinglet case_ , _Nucl. Phys._ B688 (2004) 101 [hep-ph/0403192].
  * [63] A. Vogt, S. Moch and J. A. M. Vermaseren, _The Three-loop splitting functions in QCD: The Singlet case_ , _Nucl. Phys._ B691 (2004) 129 [hep-ph/0404111].
  * [64] S. Moch, J. A. M. Vermaseren and A. Vogt, _The Three-Loop Splitting Functions in QCD: The Helicity-Dependent Case_ , _Nucl. Phys. B_ 889 (2014) 351 [1409.5131].
  * [65] S. Moch, J. A. M. Vermaseren and A. Vogt, _On ${\gamma}_{5}$ in higher-order QCD calculations and the NNLO evolution of the polarized valence distribution_, _Phys. Lett. B_ 748 (2015) 432 [1506.04517].
  * [66] J. Blümlein, P. Marquard, C. Schneider and K. Schönwald, _The three-loop unpolarized and polarized non-singlet anomalous dimensions from off shell operator matrix elements_ , _Nucl. Phys. B_ 971 (2021) 115542 [2107.06267].
  * [67] J. Blümlein, P. Marquard, C. Schneider and K. Schönwald, _The three-loop polarized singlet anomalous dimensions from off-shell operator matrix elements_ , _JHEP_ 01 (2022) 193 [2111.12401].
  * [68] W. Vogelsang, _Next-to-leading order evolution of transversity distributions and Soffer’s inequality_ , _Phys. Rev. D_ 57 (1998) 1886 [hep-ph/9706511].
  * [69] W. Vogelsang, _$Q^{2}$ evolution of spin dependent parton densities_, _Acta Phys. Polon. B_ 29 (1998) 1189 [hep-ph/9805295].
  * [70] J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M. Round et al., _The Transition Matrix Element $A_{gq}(N)$ of the Variable Flavor Number Scheme at $O(\alpha_{s}^{3})$_, _Nucl. Phys. B_ 882 (2014) 263 [1402.0359].
  * [71] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel et al., _The 3-Loop Non-Singlet Heavy Flavor Contributions and Anomalous Dimensions for the Structure Function $F_{2}(x,Q^{2})$ and Transversity_, _Nucl. Phys. B_ 886 (2014) 733 [1406.4654].
  * [72] A. Behring, I. Bierenbaum, J. Blümlein, A. De Freitas, S. Klein and F. Wißbrock, _The logarithmic contributions to the $O(\alpha^{3}_{s})$ asymptotic massive Wilson coefficients and operator matrix elements in deeply inelastic scattering_, _Eur. Phys. J. C_ 74 (2014) 3033 [1403.6356].
  * [73] M. Diehl, J. R. Gaunt, P. Plößl and A. Schäfer, _Two-loop splitting in double parton distributions_ , _SciPost Phys._ 7 (2019) 017 [1902.08019].
  * [74] M. Diehl, J. R. Gaunt and P. Plößl, _Two-loop splitting in double parton distributions: the colour non-singlet case_ , _JHEP_ 08 (2021) 040 [2105.08425].
  * [75] S. Bailey, T. Cridge, L. A. Harland-Lang, A. D. Martin and R. S. Thorne, _Parton distributions from LHC, HERA, Tevatron and fixed target data: MSHT20 PDFs_ , _Eur. Phys. J. C_ 81 (2021) 341 [2012.04684].
  * [76] ATLAS collaboration, _Measurement of the production cross section of prompt $J/\psi$ mesons in association with a $W^{\pm}$ boson in $pp$ collisions at $\sqrt{s}=$ 7 TeV with the ATLAS detector_, _JHEP_ 04 (2014) 172 [1401.2831].
  * [77] L. N. Trefethen, _Approximation Theory and Approximation Practice_. Society for Industrial and Applied Mathematics, 2012.
  * [78] P. Prince and J. Dormand, _High order embedded Runge-Kutta formulae_ , _Journal of Computational and Applied Mathematics_ 7 (1981) 67 .
  * [79] W. Giele et al., _The QCD / SM working group: Summary report_ , in _Physics at TeV colliders. Proceedings, Euro Summer School, Les Houches, France, May 21-June 1, 2001_ , hep-ph/0204316.
  * [80] M. Dittmar et al., _Working Group I: Parton distributions: Summary report for the HERA LHC Workshop Proceedings_ , hep-ph/0511119.
  * [81] M. Diehl, T. Kasemets and S. Keane, _Correlations in double parton distributions: effects of evolution_ , _JHEP_ 05 (2014) 118 [1401.1233].
  * [82] M. Diehl, R. Nagar and P. Plößl, _Quark mass effects in double parton distributions_ , 2212.07736.
  * [83] M. Diehl, P. Plößl and A. Schäfer, _Proof of sum rules for double parton distributions in QCD_ , _Eur. Phys. J. C_ 79 (2019) 253 [1811.00289].
  * [84] M. Diehl, F. Fabry and A. Vladimirov, _Two-loop evolution kernels for colour dependent double parton distributions_ , _JHEP_ 05 (2023) 067 [2212.11843].