aainstitutetext: Dipartimento di Fisica, Università di Genova and INFN,
Sezione di Genova,
Via Dodecaneso 33, 16146, Genoa, Italybbinstitutetext: Department of Physics
and Astronomy, University of Manchester, Manchester, M13 9PL,United Kingdom

# On sum rules for double and triple parton distribution functions and P
YTHIA ’s model of multiple parton interactions

Oleh Fedkevych b Jonathan R. Gaunt<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Multi-parton distributions in a proton, the nonperturbative quantities needed
to make predictions for multiple scattering rates, are poorly constrained from
theory and data and must be modelled. All Monte Carlo event generators that
simulate multiple parton interactions (e.g. P  YTHIA ) contain such a model of
multi-parton PDFs. One important theoretical constraint for the case of double
parton distributions is provided by the so-called number and momentum sum
rules. In this paper we investigate to what extent the double parton
distribution functions used in the P  YTHIA event generator obey these sum
rules. We also derive the number and momentum sum rules for the triple parton
distribution functions and discuss how one can use the P  YTHIA code to
construct triple parton distribution functions which approximately satisfy
these sum rules.

LU-TP-20-17

MCNET-22-15

##  1 Introduction

Double parton scattering (DPS) is a process where two hard interactions occur
in one hadron-hadron collision. Such processes have been measured in various
experiments for proton-proton (pp), proton-antiproton (p$\bar{\rm p}$) and
proton-nucleus (pA) collisions AxialFieldSpectrometer:1986dfj ; UA2:1991apc ;
CDF:1993sbj ; CDF:1997lmq ; CDF:1997yfa ; D0:2009apj ; ATLAS:2016ydt ;
LHCb:2012aiv ; ATLAS:2013aph ; CMS:2013huw ; D0:2014vql ; D0:2014owy ;
ATLAS:2014ofp ; CMS:2013slh ; LHCb:2015wvu ; D0:2015dyx ; D0:2015rpo ;
ATLAS:2016rnd ; CMS:2017han ; CMS:2021ijt ; LHCb:2020jse ; CMS:2021lxi ;
CMS:2022pio . DPS processes are sensitive to partonic correlations of various
types e.g. Korotkikh:2004bz ; Cattaruzza:2005nu ; Gaunt:2009re ;
Calucci:1997ii ; Frankfurt:2003td ; Calucci:1999yz ; Calucci:2009sv ;
Calucci:2010wg ; Rogers:2009ke ; Domdey:2009bg ; Flensburg:2011kj ;
Blok:2012mw ; Seymour:2013sya ; Blok:2013bpa ; Rinaldi:2018slz  and their
study can provide a deeper insight into hadron structure.

Double parton distribution functions (dPDFs) are crucial ingredients of
theoretical predictions for DPS, encoding the correlations between the two
partons arising from a given hadron111 Here, we shall use the term “dPDF” to
refer to the double parton density integrated over the partonic transverse
separation $\bm{y}$. We will discuss the definition in detail in Section 2. .
Unfortunately, due to the absence of experimental fits, the dPDFs remain
unknown. Therefore, in phenomenological studies of the DPS phenomenon one
still has to rely on certain model dependent assumptions about the shape of
dPDFs and their relation to standard single parton distribution functions
(sPDFs). There are different approaches to modelling dPDFs. To name some we
mention the Light-Front approach Rinaldi:2014ddl ; Rinaldi:2016mlk ;
Ceccopieri:2017oqe ; Rinaldi:2018bsf , the “bag model” of a proton
Chang:2012nw , the AdS/QCD method Traini:2016jru ; Rinaldi:2020ybv , chiral
quark models Broniowski:2019rmu  and the models based on the solution of so-
called “double” DGLAP evolution equations Gaunt:2009re ; Nagar:2019gij ;
Snigirev:2003cq ; Korotkikh:2004bz ; Diehl:2017kgu . First-principles
calculations of the lowest Mellin moments of double parton distributions have
also now been performed in lattice QCD Bali:2018nde ; Bali:2020mij ;
Bali:2021gel .

Apart from the aforementioned approaches to the dPDFs there exist a broad
class of Monte Carlo (MC) models of multiple partonic interactions (MPI) as
implemented in general purpose event generators, e.g. P  YTHIA
Sjostrand:2006za ; Sjostrand:2007gs ; Bierlich:2022pfr , H  ERWIG 7
Bellm:2015jjp  and S  HERPA Gleisberg:2008ta . The detailed description of MPI
models being used in P  YTHIA and H  ERWIG 7 can be found in corresponding
publications Sjostrand:1985vv ; Sjostrand:1987su ; Sjostrand:2004pf ;
Corke:2009tk ; Corke:2011yy ; Bahr:2008dy ; Gieseke:2016fpz . The MPI model of
S  HERPA is based upon the old MPI model of P  YTHIA Sjostrand:1987su  with
some modifications that are described in Ref. Gleisberg:2008ta 222 There is
also a new MPI model under development, see Ref. (Martin:2012nm, ). . These
MPI models are also widely used to simulate DPS processes.

Despite the variety of the models of dPDFs available on the market, there are
common criteria dictated by first principles QCD which, to some extent, have
to be preserved by each model of dPDFs. One such criteria is that the dPDFs
should satisfy the sum rules originally proposed by Gaunt and Stirling (GS sum
rules) Gaunt:2009re ; Blok:2013bpa ; Ceccopieri:2014ufa ; Diehl:2018kgr . In
this paper we demonstrate how one can use the MPI model of P  YTHIA to
construct sets of dPDFs which approximately obey the GS sum rules, and
quantify the extent to which these sum rules are satisfied. We start by
constructing asymmetric dPDFs and find that they satisfy the GS sum rules when
integrating over one $x$ (momentum fraction) argument at the level of a few
percent when the other $x$ argument is in between $10^{-6}$ and $0.8$. Then we
present our study of symmetrized dPDFs, obtained by symmetrizing the
asymmetric dPDFs with respect to the first and second parton arguments. We
find that these dPDFs obey the GS sum rules (when integrating over either $x$
argument) at about $20\%$ accuracy level except for $x$ fractions bigger than
0.4. As a benchmark we use the GS09 (Gaunt:2009re, ) set of dPDFs which
represents solutions of the inhomogeneous “double” DGLAP evolution equations
for an input which was constructed to approximately respect the GS sum rules.

Whereas the dPDFs and the corresponding sum rules have been extensively
studied in the literature, the case of triple parton distribution functions
(tPDFs) has received rather less attention. The evolution of the triple parton
densities was studied in Ref. Snigirev:2016uaq  and theoretical aspects of
correlations in MPI have been studied in Refs. Calucci:2009sv ; Calucci:2009ea
; Treleani:2012zi , however, the available phenomenological studies of triple
parton scattering (TPS) Maina:2009sj ; Maciula:2017wpe ; dEnterria:2017yhd ;
Maciula:2017meb ; dEnterria:2016yhy ; dEnterria:2016ids  rely on factorization
of tPDFs into products of three sPDFs which essentially neglects correlations
between partons333The recent CMS measurements CMS:2021qsn  suggest a non-
negligible TPS contribution to triple $J/\psi$ production. . Moreover, to the
best of our knowledge, no solutions of the tPDF evolution equations exist in
the literature. Therefore, in the absence of “first principles” sets of tPDFs,
it is useful to have a method to build tPDFs that account for correlations in
$x$-space related to momentum and valence number constraints in an approximate
way. In this work we demonstrate how one can use the P  YTHIA code to
accomplish this task.

This paper is organized as follows: in Section 2 we provide a short
description of the GS09 model and the model being used in the P  YTHIA event
generator and in Sections 3–4 we provide our findings. Namely, in Section 3 we
demonstrate how one can use P  YTHIA code to construct asymmetric and
symmetric dPDFs and study how well they satisfy the GS sum rules as well as
perform some toy phenomenological studies for the double Drell-Yan process.
Then, in Section 4 we derive number and momentum sum rules for the case of
TPS, and demonstrate how one can use the P  YTHIA code to create tPDFs which
approximately obey such sum rules. Finally, in Section 5 we give our
conclusions and some closing remarks.

##  2 GS09 and P  YTHIA models of dPDFs

###  2.1 Double parton distribution functions and the GS09 dPDFs

In this subsection we introduce the theoretical definition of dPDFs, and
discuss their key properties including the GS sum rules. We will also briefly
describe the GS09 dPDF set Gaunt:2009re , and discuss simplifying assumptions
that have been used in the past in order to utilise dPDFs directly in DPS
cross section predictions.

First of all, consider the total DPS cross section for the production of final
states $A$ and $B$ in the collision of hadrons $h_{A}$ and $h_{B}$:444 The
parton model derivation of the expression for the total DPS cross section was
given in Refs. (Paver:1982yp, ; Mekhfi:1983az, ; Diehl:2011tt, ; Diehl:2011yj,
) and the full QCD derivation was provided in Refs. (Diehl:2015bca, ;
Diehl:2017kgu, ; Diehl:2018wfy, ; Vladimirov:2017ksc, ).

$\displaystyle\sigma^{\rm DPS}_{AB}$
$\displaystyle=\frac{1}{1+\delta_{AB}}\sum\limits_{j_{1},\,j_{2},\,j_{3},j_{4}}\int_{|\bm{y}|>\max(Q^{-1}_{1},Q^{-1}_{2})}\,\prod\limits^{4}_{i=1}dx_{i}\,d^{2}y\,\Gamma_{j_{1}\,j_{2}/h_{A}}(x_{1},x_{2},\bm{y},Q_{1},Q_{2})$
$\displaystyle\mathrel{\phantom{=}}\times\,\Gamma_{j_{3}\,j_{4}/h_{B}}(x_{3},x_{4},\bm{y},Q_{1},Q_{2})\,\hat{\sigma}_{j_{1}\,j_{3}\rightarrow
A}\,\hat{\sigma}_{j_{2}\,j_{4}\rightarrow B}.$ (1)

Following Diehl:2011tt ; Diehl:2011yj ; Diehl:2017kgu , we shall refer to the
objects $\Gamma_{j_{1}\,j_{2}/h}(x_{1},x_{2},\bm{y},Q_{1},Q_{2})$ as the
double parton distributions (DPDs) (they are also referred to as “two parton
generalized parton distributions” or 2pGPDs in the literature Blok:2011bu ;
Blok:2012mw ; Gaunt:2012dd ). Loosely speaking, these can be interpreted as a
probability to find two partons $j_{1}$ and $j_{2}$ with longitudinal momentum
fractions $x_{1}$ and $x_{2}$ separated by transverse distance $\bm{y}$ in a
hadron $h$. A lower cut-off on the transverse separation $|\bm{y}|$ is
included since the region of very small $|\bm{y}|$ values is more
appropriately described by single parton scattering (for more details see Ref.
Diehl:2017kgu ).

Let us define the dPDFs as the integral of the DPD functions over $\bm{y}$, as
follows:

$D_{j_{1}\,j_{2}/h}(x_{1},x_{2},Q_{1},Q_{2})\equiv\int_{|\bm{y}|>\max(Q^{-1}_{1},Q^{-1}_{2})}d^{2}y\,\Gamma_{j_{1}\,j_{2}/h}(x_{1},x_{2},\bm{y},Q_{1},Q_{2}).$
(2)

The dPDFs in Eq. (2) taken at equal factorization scales
($Q_{1}=Q_{2}=Q\gg\Lambda_{QCD}$) obey so called inhomogeneous double DGLAP
(dDGLAP) evolution equations Kirschner:1979im ; Shelest:1982dg ;
Zinovev:1982be  which, at the leading logarithmic level, have the form

$\displaystyle\frac{dD_{j_{1}j_{2}}(x_{1},x_{2},t)}{dt}$
$\displaystyle=\frac{\alpha_{s}(t)}{2\pi}\sum\limits_{j^{\prime}_{1}}\,\int\limits^{1-x_{2}}_{x_{1}}\,\frac{dx^{\prime}_{1}}{x^{\prime}_{1}}\,P_{j^{\prime}_{1}\rightarrow
j_{1}}\left(\frac{x_{1}}{x^{\prime}_{1}}\right)\,D_{j^{\prime}_{1}j_{2}}(x^{\prime}_{1},x_{2},t)$
$\displaystyle\mathrel{\phantom{=}}+\,\frac{\alpha_{s}(t)}{2\pi}\sum\limits_{j^{\prime}_{2}}\,\int\limits^{1-x_{1}}_{x_{2}}\,\frac{dx^{\prime}_{2}}{x^{\prime}_{2}}\,P_{j^{\prime}_{2}\rightarrow
j_{2}}\left(\frac{x_{2}}{x^{\prime}_{2}}\right)\,D_{j_{2}j^{\prime}_{2}}(x_{1},x^{\prime}_{2},t)$
$\displaystyle\mathrel{\phantom{=}}+\,\frac{\alpha_{s}(t)}{2\pi}\,\sum\limits_{j^{\prime}}\,P_{j^{\prime}\rightarrow
j_{1}j_{2}}\left(\frac{x_{1}}{x_{1}+x_{2}}\right)\,f_{j^{\prime}}(x_{1}+x_{2},t)\,\frac{1}{x_{1}+x_{2}},$
(3)

where $f_{j^{\prime}}(x,t)$ are “standard” collinear sPDFs and we have dropped
the subscript $h_{A}$ in $D_{j_{1}j_{2}/h_{A}}$ in order not to overload the
notation555 Similar evolution equations for quark fragmentation functions were
derived in Refs. Puhala:1980ms ; Sukhatme:1980vs . . Here we use a
dimensionless evolution parameter $t=\log(Q^{2}/Q^{2}_{0})$ where $Q_{0}$ is
an arbitrary reference scale. The first two terms on the right hand side of
Eq. (3) involve standard DGLAP leading order splitting functions
$P_{j^{\prime}_{1}\rightarrow j_{1}}(x)$ and, despite the presence of two
Bjorken-$x$’es, essentially have the same structure as the standard DGLAP
splitting kernels. The last term, however, is new and is linked to the
reduction of the cut-off in Eq. (2) when $Q$ is increased. It is an
inhomogeneous term that involves the standard sPDFs and depends on
“$1\rightarrow 2$” splitting functions $P_{j^{\prime}\rightarrow
j_{1}j_{2}}(x)$ which are defined at the leading order in $\alpha_{s}$ via

$\displaystyle P_{i\rightarrow jk}(x)=P^{R}_{i\rightarrow
j}(x)\,\delta_{k\kappa(i,j)}$ (4)

where $\kappa(i,j)$ is the only parton that can be emitted at leading order
when parton $i$ splits to $j$ (e.g. if $i=g$ and $j=q$ then
$\kappa(i,j)=\bar{q}$) and $P^{R}_{i\rightarrow j}(x)$ are given by the real-
emission parts of the DGLAP splitting functions, see Ref. (Gaunt:2009re, ).
Therefore, in the following, we refer to the third term on the right hand side
of Eq. (3) as the “splitting” or $1\rightarrow 2$ term. As we have mentioned
before, the dPDFs in Eq. (3) are taken at the same values of the factorization
scales. The unequal scale dPDFs can be obtained from the equal scale dPDFs via
a single DGLAP evolution step as described in Ref. Gaunt:2009re .

The dPDFs obey the following sum rules:

$\displaystyle\sum\limits_{j_{2}}\int\limits^{1-x_{1}}_{0}dx_{2}\,x_{2}\,D_{j_{1}j_{2}}(x_{1},x_{2},Q)=(1-x_{1})\,f_{j_{1}}(x_{1},Q),$
(5)
$\displaystyle\int\limits^{1-x_{1}}_{0}dx_{2}\,D_{j_{1}j_{2v}}(x_{1},x_{2},Q)=\left(N_{j_{2v}}-\delta_{j_{1}j_{2}}+\delta_{j_{1}\bar{j}_{2}}\right)f_{j_{1}}(x_{1},Q).$
(6)

These dPDF sum rules were first introduced by Gaunt and Stirling in Ref.
Gaunt:2009re , who also established that Eqs. (5), (6) are preserved by dDGLAP
evolution provided that the dPDFs obey the sum rules at the initial evolution
scale $Q_{0}$ (see also Ref. Blok:2013bpa ). Thus, we refer to Eq. (5) and Eq.
(6) as the GS sum rules. The fundamental character of the GS sum rules was
recently demonstrated by Diehl, Plößl and Schäfer in Ref. Diehl:2018kgr  where
these relations were proved in all orders of perturbation theory for bare
dPDFs and for the dPDFs renormalized according to the $\overline{\rm MS}$
scheme (see also Ref. Gaunt:2012tfk ). The $\overline{\rm MS}$-renormalized
dPDFs coincide with the dPDFs in Eq. (2) up to corrections suppressed by
$\mathcal{O}(\alpha_{s})$ or $\mathcal{O}(\Lambda_{\rm QCD}^{2}/Q^{2})$
Diehl:2018kgr , such that the dPDFs in Eq. (2) also satisfy the GS sum rules
up to these corrections. In Ref. Gaunt:2009re , a set of dPDFs was constructed
which approximately satisfies the GS sum rules at the initial scale $Q_{0}=1$
GeV. These starting conditions were then evolved according to Eq. (3) to
higher scales, up to $Q^{2}=10^{9}$ GeV2 and for $x_{1,2}\geq 10^{-6}$. The
dPDFs generated by this procedure are referred to as the GS09 set666For some
phenomenological applications of the GS09 dPDFs see, e.g. Refs. Gaunt:2010pi ;
Maina:2010vh ; Fedkevych:2020cmd ..

If one assumes that $\Gamma_{j_{1}\,j_{2}/h}(x_{1},x_{2},\bm{y},Q_{1},Q_{2})$
can be written as a product of a longitudinal and transverse-dependent pieces,
with the transverse piece being some smoothly varying function with a width of
approximately the proton radius, then one can write the DPD in terms of the
dPDF as follows:

$\displaystyle\Gamma_{j_{1}\,j_{2}/h}(x_{1},x_{2},\bm{y},Q_{1},Q_{2})\approx
D_{j_{1}\,j_{2}/h}(x_{1},x_{2},Q_{1},Q_{2})\,F(\bm{y}),$ (7)

provided that $F(\bm{y})$ is normalised as $\int d^{2}yF(\bm{y})=1$. Under Eq.
(7), the DPS cross section in Eq. (1) may be simplified as follows:

$\displaystyle\sigma^{\rm DPS}_{AB}$ $\displaystyle=\frac{1}{\sigma_{\rm
eff}}\frac{1}{1+\delta_{AB}}\sum\limits_{j_{1},\,j_{2},\,j_{3},j_{4}}\int\,\prod\limits^{4}_{i=1}dx_{i}\,D_{j_{1}\,j_{2}/h_{A}}(x_{1},x_{2},Q_{1},Q_{2})\,D_{j_{3}\,j_{4}/h_{B}}(x_{3},x_{4},Q_{1},Q_{2})$
$\displaystyle\mathrel{\phantom{=}}\times\,\hat{\sigma}_{j_{1}\,j_{3}\rightarrow
A}\,\hat{\sigma}_{j_{2}\,j_{4}\rightarrow B},$ (8)

where $\sigma_{\rm eff}=\left[\int d^{2}y\,F^{2}(\bm{y})\right]^{-1}$ can be
seen as an effective DPS interaction area.

Nowadays, it is known that the factorization of the DPDs into the dPDFs and a
transverse-dependent piece cannot hold Blok:2011bu ; Gaunt:2011xd ;
Blok:2013bpa ; Gaunt:2012dd ; Ryskin:2011kk ; Ryskin:2012qx ; Manohar:2012pe ;
Diehl:2017kgu . Actually, the discussion above already indicates this – if Eq.
(7) held, then the cut-off in Eq. (2) would only have a power-suppressed
effect, and there would be no inhomogeneous term in Eq. (3). Nevertheless,
when using the GS09 dPDF set in cross section calculations below we will use
Eq. (8) – this is because we would like to make the most direct comparison
against the predictions of P  YTHIA , which uses a framework for DPS/MPI that
is most closely aligned with Eq. (8).

It is important to remark that although dPDFs satisfying the GS sum rules do
not directly appear in the DPS cross section Eq. (1), they are linked to the
DPDs that do, via Eq. (2). Any procedure that yields dPDFs that
(approximately) satisfy the sum rules is thus of interest and utility even
under the theoretically rigorous approach associated with Eq. (1). It is also
worth remarking that it is possible to adapt techniques used to construct
dPDFs satisfying the sum rules to construct DPDs which have both the
appropriate transverse dependence and which integrate to dPDFs satisfying the
sum rules – in Ref. Diehl:2020xyg  the techniques used to construct the GS09
dPDFs were adapted, extended and improved to yield such a set of DPDs.

###  2.2 The PYTHIA model of multiple parton distribution functions

Now let us briefly sketch the way the P  YTHIA event generator models
$n$-parton distribution functions (nPDFs) and, more specifically, dPDFs. We
remind the reader that these dPDFs (and $n$PDFs) are applied to compute cross
sections according to a framework that resembles Eq. (8); the transverse
dependence is assumed to essentially factor out. The generation of the first
hard interaction is performed in the usual way, for example, for the case of
massless $2\rightarrow 2$ scattering we have

$\displaystyle\frac{d\sigma}{dp^{2}_{\perp}}=\sum_{j_{1}\,j_{2}\,j_{3}\,j_{4}}\int
dx_{1}\,dx_{2}\,d\hat{t}\,f_{j_{1}}(x_{1},Q)\,f_{j_{2}}(x_{2},Q)\,\frac{d\sigma_{j_{1}\,j_{2}\rightarrow
j_{3}\,j_{4}}}{d\hat{t}}\,\delta\left(p^{2}_{\perp}-\frac{\hat{t}\hat{u}}{\hat{s}}\right).$
(9)

After the successful generation of the first hard interaction, one has to
generate other subsequent interactions by taking into account MPIs already
generated in a given pp collision. In the general MPI model of P  YTHIA these
interactions are generated at successively lower scales. There is also an
option to generate DPS events specifically (controlled by the option
“SecondHard”) – here both orderings of the two processes are simulated (see
the beginning of Section 3.3 for more details). To generate the additional
interaction(s) P  YTHIA dynamically modifies sPDFs in Eq. (9) according to a
history of previous interactions. In the following, while describing P  YTHIA
’s framework, we use subscripts “$r$” and “$m$” to distinguish between
unmodified (“raw”) and modified sPDFs.

In the first model of MPI (Sjostrand:1985vv, ; Sjostrand:1987su, ) P  YTHIA
could only preserve overall energy and momentum by reducing the value of
Bjorken-$x$ after each interaction according to

$\displaystyle x^{\prime}_{i}=\frac{x_{i}}{X}\,<\,1,$ (10)

where $x^{\prime}_{i}$ is the value of Bjorken-$x$ used for the $i$’th
interaction and $X\equiv 1-\sum^{i-1}_{j=1}x_{j}$ is a total longitudinal
momentum already taken by previous interactions. For example, if $i=2$ we have

$\displaystyle x^{\prime}_{2}=\frac{x_{2}}{1-x_{1}}\,<\,1,$ (11)

which implies $x_{1}+x_{2}<1$; thus, this procedure corresponds to having a
dPDF given by

$\displaystyle
D_{j_{1}\,j_{2}}(x_{1},x_{2},Q)=f_{j_{1}}^{r}(x_{1},Q)\,f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q),$ (12)

where $f_{j_{1}}^{r}(x_{1},Q)$ is the unmodified sPDF and

$\displaystyle f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)=\frac{1}{1-x_{1}}f_{j_{2}}^{r}\left(\frac{x_{2}}{1-x_{1}},Q\right)$
(13)

is the sPDF modified according to the first interaction involving parton
$j_{1}$ with momentum fraction $x_{1}$. Therefore,

$\displaystyle\sum\limits_{j_{2}}\,\int\limits^{1}_{0}dx_{2}\,\theta(1-x_{1}-x_{2})\,x_{2}\,f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)=1-x_{1}.$ (14)

That is, the momentum of all partons after the first interaction is $1-x_{1}$,
as desired. In general, after $n$ interactions we will have:

$\displaystyle\int\limits^{X}_{0}dx\,x\left(\sum\limits_{f}\,q_{fn}(x,Q)+g_{n}(x,Q)\right)=X.$
(15)

This approach, however, does not take into account changes in the quark
content of a proton. For example, if one probes a “valence” $d$-quark in the
first hard interaction then for the second interaction the number of “valence”
$d$-quarks in a beam remnant should be reduced to zero. Furthermore, if one
probes a $d$-quark sPDF at low-$x$ it means that with a high probability the
$d$-quark originates from a perturbative $1\rightarrow 2$ splitting of a gluon
into a quark-antiquark pair, which, in turn, implies that there is a left-over
“companion” antiquark which should be included in a beam remnant, see Fig. 1.
These effects, obviously, apply to all quark flavours.

a)

b)

Figure 1: Two different ways to generate $\left(d\,g\rightarrow
d\,g\right)\otimes\left(\bar{d}\,g\rightarrow\bar{d}\,g\right)$ DPS process in
the P

YTHIA

event generator: a) the initial state $d$-quark comes from a perturbative
$g\rightarrow d\,\bar{d}$ splitting process, b) the initial state $d$-quark
comes from a raw $d$-quark sPDF.

In 2004 the MPI model of P  YTHIA was significantly improved: $n$-parton PDFs
of arbitrarily high order $n$ were constructed such that changes in quark
content due to higher-scale interactions were taken into account. For the
first interaction a standard sPDF is assumed to be valid, as fitted to data.
The sPDF for the second interaction thereafter depends on the properties of
the first interaction, and so on Sjostrand:2004pf . The starting point of this
approach is to split sPDFs into valence and sea parts as

$\displaystyle f(x,Q)=f_{\mathrm{val}}(x,Q)+f_{\mathrm{sea}}(x,Q),$ (16)

where the flavour-antiflavour symmetry in the sea is assumed. If we denote the
raw valence sPDF for flavour $f$ by $q_{fv0}(x,Q)$, then the valence sPDF
after $n$ interactions is given by

$\displaystyle
q_{fvn}(x,Q)=\frac{N_{fvn}}{N_{fv0}}\,\frac{1}{X}\,q_{fv0}(x/X,Q),$ (17)

where $N_{fv0}$ is an original number of valence quarks of a flavour $f$,
$q_{fvn}$ and $N_{fvn}$ are distribution and number of valence quarks of the
flavour $f$ after $n$ interactions correspondingly. One can see that
modification in Eq. (17) changes the standard number rule

$\displaystyle\int\limits^{1}_{0}dx\,q_{fv0}(x,Q)=N_{fv0}$ (18)

to

$\displaystyle\int\limits^{1}_{0}dx\,q_{fvn}(x,Q)\,\theta(X-x)=\frac{N_{fvn}}{N_{fv0}}\,\int\limits^{X}_{0}\,\frac{dx}{X}\,q_{fv0}(x/X,Q)=N_{fvn}.$
(19)

Equation (17) allows to account for changes in the number of valence quarks.
However, it cannot be used on its own since a subtraction of a valence quark
also removes an average momentum fraction carried by $q_{fvn}(x,Q)$

$\displaystyle\langle
x_{fvn}(Q)\rangle\equiv\frac{\int\limits^{X}_{0}dx\,x\,q_{fvn}(x,Q)}{\int\limits^{X}_{0}dx\,q_{fvn}(x,Q)}=X\langle
x_{fv0}(Q)\rangle,$ (20)

and this will cause the algorithm to violate the momentum rule Eq. (15) (which
would be preserved if we only rescaled the momentum fraction in Eq. (17)). We
will describe how this problem is resolved shortly; however, let us first
describe how the $1\to 2$ splitting effects described earlier are accounted
for in the P  YTHIA code. Let us say that in a particular interaction we find
a sea quark $q_{s}$ in the proton. This must have been produced in conjunction
with a companion (anti-)quark, $q_{c}$:

$\displaystyle g\rightarrow q_{s}+q_{c}.$ (21)

Let us denote the momentum fraction of the parent gluon by $x_{g}$, that of
the sea quark by $x_{s}$ and that of the companion by $x=x_{g}-x_{s}$. Then,
according to Ref. Sjostrand:2004pf , we introduce an sPDF for the companion
quark, which is given by777 Note the resemblance between the structure of the
right hand side of Eq. (22) and the splitting term on the right hand side of
dDGLAP evolution equations, Eq. (3) – this is of course no coincidence, as
both refer to a perturbative $1\to 2$ splitting.

$\displaystyle q_{c}(x,x_{s})=C\,P_{g\rightarrow
q\bar{q}}\left(\frac{x_{s}}{x_{s}+x}\right)\frac{g(x_{s}+x)}{x_{s}+x},$ (22)

where $P_{g\rightarrow
q\bar{q}}\left(z\right)=\frac{1}{2}\left[z^{2}+(1-z)^{2}\right]$ and $C$ is a
normalization constant which is fixed by

$\displaystyle\int\limits^{1-x_{s}}_{0}dx\,q_{c}(x,x_{s})=1.$ (23)

Note that Eq. (22) is defined assuming no momentum loss due to the previous
hard interactions and hence has to be rescaled accordingly.

Using Eq. (22) one can compute an average momentum fraction introduced by a
companion quark to a beam remnant $X\langle x_{fc0}(x_{s})\rangle$ defined in
analogy with Eq. (20) as

$\displaystyle X\langle
x_{fc0}(x_{s})\rangle\equiv\frac{\int\limits^{X}_{0}dx\,x\,q_{fcn}(x,x_{s})}{\int\limits^{X}_{0}dx\,q_{fcn}(x,x_{s})}=X\int\limits^{1}_{0}dx\,xq_{fc0}(x,x_{s}),$
(24)

where index $f$ is used to distinguish between different flavours and
$q_{fcn}$ is linked to $q_{fc0}$ via the momentum rescaling:
$q_{fcn}(x,x_{s})\equiv q_{fc0}(x/X,x_{s})/X$. The average momentum introduced
by a companion quark together with the average momenta subtracted due to the
removal of valence quarks $X\langle x_{fv0}(Q)$ leads to the violation of Eq.
(15). The way to circumvent this issue is to let the sea and gluon
distributions fluctuate in such a way that both contributions are compensated
and the momentum rule Eq. (15) is preserved together with the number rules

$\displaystyle\int\limits^{X}_{0}dx\,q_{fvn}(x,Q)=N_{fvn},$ (25)
$\displaystyle\int\limits^{X}_{0}dx\,q_{fc_{j}n}(x,x_{j})=1,\forall j,$ (26)

where index $j$ is used to distinguish between different companion quarks of
the flavour $f$. The corresponding modifications of the sea and gluon PDFs are
given by

$\displaystyle q_{fs}(x,Q)\rightarrow aq_{fs}(x,Q),$ (27) $\displaystyle
g(x,Q)\rightarrow ag(x,Q),$ (28)

which allows one to rewrite the momentum rule as

$\displaystyle 1$ $\displaystyle=$
$\displaystyle\frac{1}{X}\int\limits^{X}_{0}dx\,x\left(\sum\limits_{f}\left[q_{fvn}(x,Q)+\sum\limits_{j}q_{fc_{j}n}(x,x_{j})+aq_{fs}(x,Q)\right]+ag_{n}(x,Q)\right)=$
(29) $\displaystyle=$ $\displaystyle a\left(1-\sum\limits_{f}N_{fv0}\langle
x_{fv0}\rangle\right)+\sum\limits_{f}N_{fvn}\langle
x_{fv0}\rangle+\sum\limits_{f,j}\langle x_{fc_{j}0}\rangle,$

where $\langle x_{fv0}\rangle$ and $\langle x_{fc_{j}0}\rangle$ are defined
according to Eqs. (20) and (24) respectively. Equation (29) fixes the value of
the constant $a$ in Eqs. (27) and (28) which allows to build a closed system
of equations valid for all sPDFs used to simulate an $n$’th interaction.

##  3 Comparison between the GS09 and P  YTHIA models of dPDFs

###  3.1 Asymmetric double parton distribution functions

Now let us study in detail how well dPDFs constructed out of the sPDFs being
used in the MPI model of P  YTHIA satisfy the GS sum rules. We shall note that
one can use raw and modified sPDFs to construct asymmetric P  YTHIA dPDFs as
in Eq. (12), or one can construct symmetric dPDFs according to

$\displaystyle D^{\rm
sym}_{j_{1}j_{2}}(x_{1},x_{2},Q)\equiv\frac{D_{j_{1}j_{2}}(x_{1},x_{2},Q)+D_{j_{2}j_{1}}(x_{2},x_{1},Q)}{2},$
(30)

where $D_{j_{1}j_{2}}(x_{1},x_{2},Q)$ and $D_{j_{2}j_{1}}(x_{2},x_{1},Q)$ are
given by Eq. (12). The asymmetric dPDFs have a strict ordering of the
arguments since the second sPDF $f^{m\leftarrow j_{1},x_{1}}_{j_{2}}(x_{2},Q)$
is always modified according to the first hard interaction involving a parton
$j_{1}$ carrying a momentum fraction $x_{1}$. Such ordering appears because in
the MPI model of P  YTHIA the MPIs are considered as semi-hard corrections to
the first hard process. However, one should keep in mind that the dPDFs should
be symmetric under simultaneous interchange of flavours, Bjorken-$x$’es and
factorization scales, and hence the GS sum rules should be satisfied when
either $x$ fraction is integrated over for the equal scale dPDFs. Nonetheless,
it is interesting to see to what extent the asymmetric dPDFs satisfy the sum
rules when integrating over the momentum fraction of the second “modified”
parton, as this is the scenario where, by design, the P  YTHIA dPDFs should
best satisfy the sum rules. Therefore, we start the discussion with the
asymmetric dPDFs as in Eq. (12) and then consider effects of symmetrization.
Since the sum rules apply to the case where the scales in the dPDF are equal,
we set $Q_{1}=Q_{2}$ and set this common scale to $91$ GeV.

First of all, let us check how the ansatz given by Eq. (12) satisfies the
momentum sum rule given by Eq. (5) (when integrating over the second parton
only). By substituting Eq. (12) into Eq. (5) we get

$\displaystyle\sum\limits_{j_{2}}\int\limits^{1-x_{1}}_{0}dx_{2}\,x_{2}\,f_{j_{1}}^{r}(x_{1},Q)f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)\overset{?}{=}(1-x_{1})\,f_{j_{1}}^{r}(x_{1},Q),$ (31)

or

$\displaystyle\sum\limits_{j_{2}}\int\limits^{1-x_{1}}_{0}dx_{2}\,x_{2}\,f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)\overset{?}{=}1-x_{1}.$ (32)

where we use the symbol $\overset{?}{=}$ to denote that this relation should
be true for the sum rule to be satisfied. Since $j_{2}$ in $\sum_{j_{2}}$ runs
over all parton species one can write Eq. (32) as

$\displaystyle\int\limits^{1-x_{1}}_{0}dx_{2}\,x_{2}\,\left(\sum\limits_{j=q,\bar{q}}f_{j}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)+f_{g}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)\right)\overset{?}{=}1-x_{1}.$ (33)

Looking at Eq. (15) with $n=2$, $X=1-x_{1}$, we see that this relation does
indeed hold. That is, the asymmetric P  YTHIA dPDFs satisfy the momentum sum
rule, when integrating over the second parton, by design. Later in this
section we will investigate how well Eq. (5) is satisfied by the asymmetric P
YTHIA dPDFs in practice. However, before doing that, let us first consider the
number sum rule given by Eq. (6). It is convenient to express Eq. (6) in terms
of a response function $R_{j_{1}j_{2}}$ defined as

$\displaystyle R_{j_{1}j_{2}}(x_{1},x_{2},Q)\equiv
x_{2}\,\frac{D_{j_{1}j_{2}}(x_{1},x_{2},Q)-D_{j_{1}\bar{j_{2}}}(x_{1},x_{2},Q)}{f_{j_{1}}^{r}(x_{1},Q)}.$
(34)

The meaning of $R_{j_{1}j_{2}}(x_{1},x_{2},Q)$ becomes clear if one neglects
momentum and number conservation and substitutes
$D_{j_{1}j_{2}}(x_{1},x_{2},Q)=f_{j_{1}}^{r}(x_{1},Q)f_{j_{2}}^{r}(x_{2},Q)$
in Eq. (34). The function $R_{j_{1}j_{2}}(x_{1},x_{2},Q)$ then turns simply
into a raw valence quark sPDF multiplied by a corresponding Bjorken-$x$:

$\displaystyle
R_{j_{1}j_{2}}(x_{1},x_{2},Q)=x_{2}\left[f_{j_{2}}^{r}(x_{2},Q)-f_{\bar{j_{2}}}^{r}(x_{2},Q)\right]=x_{2}f^{r}_{j_{2v}}(x_{2},Q).$
(35)

Using $R_{j_{1}j_{2}}(x_{1},x_{2},Q)$ one can write Eq. (6) as

$\displaystyle\int\limits^{1-x_{1}}_{0}\frac{dx_{2}}{x_{2}}\,R_{j_{1}j_{2}}(x_{1},x_{2},Q)\overset{?}{=}N_{j_{2v}}-\delta_{j_{1}j_{2}}+\delta_{j_{1}\bar{j}_{2}}.$
(36)

Since in the P  YTHIA model
$D_{j_{1}j_{2}}(x_{1},x_{2},Q)=f_{j_{1}}^{r}(x_{1},Q)f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)$ one can write $R_{j_{1}j_{2}}(x_{1},x_{2},Q)$ as

$\displaystyle
R_{j_{1}j_{2}}(x_{1},x_{2},Q)=x_{2}\,\left[f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)-f_{\bar{j_{2}}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)\right],$ (37)

which turns Eq. (36) into

$\displaystyle\int\limits^{1-x_{1}}_{0}dx_{2}\,\left[f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)-f_{\bar{j_{2}}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)\right]\overset{?}{=}N_{j_{2v}}-\delta_{j_{1}j_{2}}+\delta_{j_{1}\bar{j}_{2}},$
(38)

where one can consider the difference $f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)-f^{m\leftarrow j_{1},x_{1}}_{\bar{j_{2}}}(x_{2},Q)$ as
the response of the sPDF of a valence quark of a flavour $j_{2}$ to the first
hard interaction involving a parton of the flavour $f_{1}$ probed at the
Bjorken-$x$ equal to $x_{1}$.

In order to be able to perform numerical integration in Eq. (38) one has to be
able to access sPDFs modified by P  YTHIA to simulate a second hard
interaction. Normally, it is not possible since the corresponding part of the
MPI model is hidden from a standard P  YTHIA user. However, it can be achieved
by instantiation of objects which are members of the `BeamParticle` class as
explained in Appendix A. Each call of the routine `second_xPDF` from Appendix
A returns the value of the modified sPDF being used for the second interaction
which allows us to perform numerical integration in Eq. (38). One important
comment has to be made: the aforementioned method to get the value of modified
sPDF is essentially probabilistic which implies that the modified sPDF being
plotted as functions of $x_{2}$ will have strong oscillations which make
numerical integration in Eq. (38) unstable. Such oscillations appear due to
the different treatment of valence and sea sPDFs in the DPS (MPI) machinery of
P  YTHIA . Therefore, depending on the value of $x_{1}$ P  YTHIA decides in a
probabilistic manner if $f^{m\leftarrow j_{1},x_{1}}_{j_{2}}(x_{2},Q)$ should
be evaluated assuming $f^{r}_{j_{1}}(x_{1},Q)$ to be either valence or sea
sPDF. The value of $f^{m\leftarrow j_{1},x_{1}}_{j_{2}}(x_{2},Q)$ will differ
for each case which leads to the oscillations of $f^{m\leftarrow
j_{1},x_{1}}_{j_{2}}(x_{2},Q)$ as soon as $x_{1}$ approaches the region where
the valence and sea sPDFs overlap. However, the problem can be solved by
considering the values of $f^{m\leftarrow j_{1},x_{1}}_{j_{2}}(x_{2},Q)$
averaged over a sufficiently high number of calls of `second_xPDF` routine at
given values of $x_{1}$ and $x_{2}$. Therefore, in the following, we always
average $f^{m\leftarrow j_{1},x_{1}}_{j_{2}}(x_{2},Q)$ before performing
numerical integration (we use $10^{6}$ and $10^{4}$ averaging calls for the
checks of the number and momentum rules correspondingly).

Apart from the dPDFs constructed out of the P  YTHIA sPDFs we also consider
the GS09 dPDFs. In order to check how well GS09 dPDFs satisfy the number sum
rule we use the publicly available GS09 grids with 600 grid points GS09web .
Since the GS09 dPDFs are based upon the MSTW2008 LO PDF set Martin:2009iq  we
use it in all our computations involving sPDFs.

Additionally, as a baseline, we use so-called “naive” model of dPDFs where one
replaces a dPDF by a product of two sPDFs and a $\theta$-function to enforce
the kinematic constraint $x_{1}+x_{2}\leq 1$:

$\displaystyle
D_{j_{1}j_{2}}(x_{1},x_{2},Q)=f^{r}_{j_{1}}(x_{1},Q)\,f^{r}_{j_{2}}(x_{2},Q)\,\theta(1-x_{1}-x_{2}).$
(39)

This approach neglects correlations in $x$-space and violates the GS sum
rules. Moreover, as was shown in Refs. (Kirschner:1979im, ; Shelest:1982dg, ;
Zinovev:1982be, ), this ansatz does not satisfy the dDGLAP evolution
equations.

Now let us check how well P  YTHIA and GS09 models satisfy the number rule
given by Eq. (6) taking as an example the $R_{uu}(x_{1},x_{2},Q)$ and
$R_{dd}(x_{1},x_{2},Q)$ response functions. In Fig. 2 we compare response
functions evaluated with the P  YTHIA and GS09 dPDFs against raw valence $u$\-
and $d$-quark sPDFs. For our checks we consider six different fixed values of
$x_{1}$, namely $10^{-6},\,10^{-3},\,10^{-1},\,0.2,\,0.4$ and $0.8$. One can
see how the $1\rightarrow 2$ splitting, simulated according to Eq. (22), leads
to negative values of the response functions at small values of $x_{2}$. For
example, if in the first interaction we probe a $u$-quark at $x_{1}=10^{-6}$
it is most likely a sea quark which was created due to a perturbative
splitting of a gluon into a $u\bar{u}$ pair. Therefore, the remaining
companion $\bar{u}$-quark has to be included into a $\bar{u}$-quark sPDF which
means that $f_{\bar{u}}^{r}(x_{2},Q)\,<\,f_{\bar{u}}^{m\leftarrow
u,x_{1}}(x_{2},Q)$ for $x_{2}\sim 10^{-6}$. It implies that
$R_{uu}(x_{1},x_{2},Q)$ becomes negative if $x_{2}\sim 10^{-6}$. In Fig. 2 one
can see how the minimum of $R_{uu}(x_{1},x_{2},Q)$ propagates with the value
of $x_{1}$ (see $x_{1}=10^{-6}$ and $x_{1}=10^{-3}$ cases). The same
considerations, obviously, apply to the $g\rightarrow d\bar{d}$ splitting. By
contrast, when $x_{1}\geq 10^{-1}$ it becomes more likely to pick a valence
quark and, therefore, P  YTHIA starts to use Eq. (17) more frequently than Eq.
(22), which results in a decrease of the peak of the valence distributions. We
also would like to note that in the case of the $R_{uu}$ response function we
observe almost a factor of two reduction in the height of the peak of the
modified valence quark distribution comparing to the raw valence distribution
if $x_{1}=0.1$ and that in case of the $R_{dd}$ response function the peak of
the modified valence quark distribution almost disappears completely (as one
would naively expect).

a)

b)

Figure 2: The responses of the valence quarks as functions of $x_{2}$ at fixed
$x_{1}$. The factorization scale $Q=91$ GeV. Comparison between the GS09 and P

YTHIA

dPDFs: a) $R_{uu}(x_{1},x_{2},Q)$ response function, b)
$R_{dd}(x_{1},x_{2},Q)$ response function.

The results given in Fig. 2 demonstrate that P  YTHIA and GS09 models give
similar predictions for the response of the valence $u$\- and $d$-quarks. Now
let us check numerically how well both models satisfy the number sum rule. In
Table 1 and Table 2 we provide a numerical check of Eq. (6) for the $u$\- and
$d$-quarks correspondingly. We see that both models satisfy Eq. (6) at the
accuracy level of a few percent for most values of $x_{1}$. The violation of
Eq. (6) by GS09 dPDFs for the case $x_{1}=10^{-6}$ is a boundary effect due to
the lower limit of the GS09 grids we use. One can see in Fig. 2 that if
$x_{1}=10^{-6}$, a significant part of the response function “leaks” outside
of the plotting range (compare $x_{1}=10^{-6}$ and $x_{1}=10^{-3}$ cases).
When working with the P  YTHIA dPDFs we set the lower integration limit to
$10^{-8}$ in order to avoid the aforementioned boundary effects (the lowest
accessible $x$ value in P  YTHIA is equal to $10^{-10}$). Comparing results of
numerical integration obtained with P  YTHIA and GS09 dPDFs for the
$x_{1}=10^{-6}$ case shown in Tables 1 and 2 we see that by setting a lower
integration limit to $10^{-8}$ in the P  YTHIA code one can avoid violation of
the number sum rule due to the boundary effects.

$x_{1}$ |  $N_{u_{v1}}$ P  YTHIA |  $N_{u_{v1}}$ GS09 |  $N_{u_{v1}}$ Naive
---|---|---|---
$10^{-6}$ | 1.022 | 1.526 | 2.006
$10^{-3}$ | 1.009 | 1.000 | 2.006
$10^{-1}$ | 1.003 | 1.036 | 2.005
0.2 | 1.002 | 1.016 | 2.005
0.4 | 1.001 | 0.949 | 1.997
0.8 | 0.999 | 0.862 | 1.708
Table 1: Numerical integration with respect to $x_{2}$ over $R_{uu}(x_{1},x_{2},Q)$ response function at fixed $x_{1}$. Here $N_{u_{v1}}$ stands for remaining number of valence $u$-quarks after the first hard interaction. In the ideal situation when the GS sum rules are perfectly satisfied $N_{u_{v1}}=1$.  $x_{1}$ |  $N_{d_{v1}}$ P  YTHIA |  $N_{d_{v1}}$ GS09 |  $N_{d_{v1}}$ Naive
---|---|---|---
$10^{-6}$ | 0.015 | 0.633 | 0.999
$10^{-3}$ | 0.001 | 0.029 | 0.999
$10^{-1}$ | -0.002 | 0.004 | 0.999
0.2 | -0.001 | 0.006 | 0.999
0.4 | -0.001 | 0.002 | 0.998
0.8 | -0.001 | 0.000 | 0.897
Table 2: Same as in Table 1 but for the $R_{dd}(x_{1},x_{2},Q)$ response
function. In the ideal situation when the GS sum rules are perfectly satisfied
$N_{d_{v1}}=0$.

Now let us check numerically how well the GS09 and asymmetric P  YTHIA dPDFs
satisfy the momentum sum rule. In order to do that we write Eq. (5) and Eq.
(32) as

$\displaystyle\frac{1}{1-x_{1}}\int\limits^{1-x_{1}}_{0}dx_{2}\,x_{2}\,\sum\limits_{j_{2}}\,\frac{D_{j_{1}j_{2}}(x_{1},x_{2},Q)}{f^{r}_{j_{1}}(x_{1},Q)}\overset{?}{=}1,$
(40)
$\displaystyle\frac{1}{1-x_{1}}\int\limits^{1-x_{1}}_{0}dx_{2}\,x_{2}\,\sum\limits_{j}f_{j}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)\overset{?}{=}1,$ (41)

and check how well the LHS of Eq. (40) and Eq. (41) match unity when we insert
the GS09 dPDFs and the modified P  YTHIA sPDFs respectively.

The results of the numerical checks are given in Table 3. We see that, as
expected from Eq. (33), P  YTHIA dPDFs obey the momentum rule at the permille
accuracy level and that GS09 dPDFs satisfy the momentum rule at a few percent
accuracy level. We also note that the naive dPDFs obey the momentum sum rule
at the few permille accuracy level if $x_{1}$ is in the range
$10^{-6}-10^{-3}$, but show a strong disagreement if $x_{1}\geq 0.1$. The
reason for this becomes clear if one substitutes Eq. (39) in Eq. (5) which
leads to

$\displaystyle\frac{1}{1-x_{1}}\int\limits^{1-x_{1}}_{0}dx_{2}\,x_{2}\,\sum\limits_{j_{2}}f^{r}_{j_{2}}(x_{2},Q)\overset{?}{=}1,$
(42)

which holds only in the $x_{1}\rightarrow 0$ limit.

$x_{1}$ | $j_{1}$ | P  YTHIA dPDFs. | GS09 dPDFs. | Naive dPDFs.
---|---|---|---|---
$10^{-6}$ | $d$ | 0.996 | 0.931 | 0.996
$10^{-3}$ | $d$ | 0.996 | 1.009 | 0.997
$10^{-1}$ | $d$ | 1.003 | 1.026 | 1.106
0.2 | $d$ | 1.007 | 1.006 | 1.244
0.4 | $d$ | 1.009 | 0.980 | 1.649
0.8 | $d$ | 1.009 | 0.854 | 4.131

Table 3: Check of the momentum sum rule for the P

YTHIA

, GS09 and naive models of dPDFs. In the ideal situation when the momentum sum
rule is perfectly satisfied the number in each cell should be equal to unity.

Similarly to the case of the number sum rule, it is interesting to visualize
the behaviour of the integrands in Eqs. (40), (41). In Fig. 3 we plot the
response function for the momentum sum rule:

$\displaystyle\mathcal{R}_{j_{1}}\equiv\frac{1}{1-x_{1}}\sum_{j_{2}}\frac{D_{j_{1}j_{2}}(x_{1},x_{2},Q)}{f^{r}_{j_{1}}(x_{1},Q)}\overset{\text{P\scalebox{0.8}{YTHIA}}}{=}\frac{1}{1-x_{1}}\,\sum\limits_{j_{2}}f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)$ (43)

as a function of $x_{2}$ for six fixed values of $x_{1}$, and for the case in
which $j_{1}=d$. In Fig. 3 a) we compare results obtained with P  YTHIA dPDFs
(solid lines) and GS09 dPDFs (dotted lines) whereas in Fig. 3 b) we compare P
YTHIA predictions against results obtained with naive dPDFs (dashed lines).
The blue lines in Fig. 3 correspond to the sum over all flavours in Eq. (43)
whereas red, magenta and orange lines stand for partial contributions from
gluon, $d$\- and $\bar{d}$-quarks correspondingly. Finally, the green lines
are given by contribution from all quarks except $d$\- and $\bar{d}$-flavours.

a)

b)

Figure 3: Behaviour of the momentum sum rule response function, as defined in
Eq. (43), for the dPDF models considered. Solid, dotted and dashed lines
correspond to the P

YTHIA

dPDFs, GS09 and naive dPDFs correspondingly.

By comparing results obtained with P  YTHIA and GS09 dPDFs shown in Fig. 3 a)
we see that both models give similar predictions for all partial contributions
except in the $x_{1}=10^{-6}$ and $x_{1}=0.8$ cases where the GS09 curves lie
consistently below the P  YTHIA ones. It is also noteworthy that, for all
$x_{1}$ values probed, the bulk of area under the curves lies above
$x_{2}=10^{-6}$ (in contrast to the corresponding plots for the number sum
rule in Fig. 2, where we find that for $x_{1}=10^{-6}$ an important part of
the curve lies below $x_{2}=10^{-6}$). This means that the value of 0.931 at
$x_{1}=10^{-6}$ in Table 3 cannot be ascribed to the fact that the GS09 grids
are limited to $x_{i}>10^{-6}$, but is instead linked to the GS09 curve being
somewhat “too low” at larger $x_{2}$ values $\sim 10^{-2}-10^{-1}$ for
$x_{1}=10^{-6}$. This deficiency in the extent to which the GS09 set satisfies
the momentum sum rule at the lowest $x_{1}$ values only appears after
evolution – when one constructs Table 3 at $Q=1$ GeV, one finds a value of
0.997 at $x_{1}=10^{-6}$, and the sum rule integrand curve for GS09 lies
almost exactly on top of the P  YTHIA one for this $x_{1}$ value. We were not
able to isolate a simple reason why this deficiency appears after evolution –
it is presumably due to deficiencies at larger $x_{i}$ values migrating
downwards, but the region $x_{1}=10^{-6},x_{2}=10^{-2}-10^{-1}$ is fed by
various $x$ values and flavours during evolution, making it very difficult to
identify what features of the initial conditions are driving this effect.
Since the main goal of this paper is in any case to study dPDFs and tPDFs
generated using the P  YTHIA model, and the GS09 results are used as a
baseline, we defer further investigation of this phenomenon to future studies.

Now let us discuss the momentum response functions for the naive dPDFs, as
shown in Fig. 3 b). We see that if $x_{1}$ is equal to $10^{-6}$ or $10^{-3}$
there is no visible difference between the naive and P  YTHIA response
functions (solid and dashed lines practically coincide) which is in agreement
with the results of the numerical integration shown in Table 3. We also see
that starting from $x_{1}=10^{-1}$ the naive model starts to consistently
overestimate the sum rule integrand. One would expect such an overestimate
from the fact that the naive momentum sum rule curves do not go to zero at
$x_{2}=1-x_{1}$ as they should, but their trajectories are such that they
would go to zero at $x_{2}=1$, if there were not the $\theta$-function cutting
them off. This effect is hardly noticeable for very small $x_{1}\lesssim
10^{-1}$, but becomes more and more significant at large $x_{1}$ – see e.g.
the $x_{1}=0.8$ curve of Fig. 3 b). This increasing overshoot results in the
naive model momentum sum rule integral deviating more and more from $1$ as
$x_{1}$ increases, as we observed in Table 3.

Before moving on, let us briefly comment on the case of unequal scale dPDFs in
P  YTHIA . Using the machinery discussed in Appendix A, we can access
asymmetric unequal-scale dPDFs, with either $Q_{1}<Q_{2}$ or $Q_{1}>Q_{2}$.
The sPDF of the “second” parton is always modified as in Section 2.2. For such
dPDFs it is clear that, when integrating over the second parton, the momentum
and number and sum rules will hold for the unequal scale P  YTHIA dPDFs to a
very similar extent as they held for the equal scale P  YTHIA dPDFs just
discussed888We also checked this directly for the $Q_{1}=45$, $Q_{2}=91$ GeV
and $Q_{2}=45$, $Q_{1}=91$ GeV cases.. It is worth remarking that it is not
necessarily expected theoretically that the dPDFs should satisfy the GS09 sum
rules for the case of $Q_{1}\neq Q_{2}$. If one follows the procedure
discussed in Section 2.1 to obtain the unequal scale dPDFs from the equal
scale dPDFs (and it is the case that those equal scale dPDFs satisfy the sum
rules) then the unequal scale dPDFs will satisfy the sum rules when
integrating over the parton with higher $Q$, but not when integrating over the
parton with lower $Q$.

###  3.2 Symmetrized double parton distribution functions

In Section 3.1 we demonstrated that the asymmetric dPDFs constructed as in Eq.
(12) satisfy the GS sum rules in the second argument at a few percent accuracy
level in the wide range of Bjorken-$x$’s from $10^{-6}$ up to $0.8$. We also
found that both P  YTHIA and GS09 dPDFs have similar behaviour of the response
functions as shown in Fig. 2 and Fig. 3. However, as we have mentioned before,
QCD requires dPDFs to be symmetric under simultaneous interchange of parton
flavours, Bjorken-$x$’es and factorization scales. Obviously, the dPDFs
defined by Eq. (12) do not satisfy this requirement; in this section we check
how well the symmetrized dPDFs defined in as Eq. (30) satisfy the GS sum
rules. We start with the momentum sum rule. By substituting Eq. (30) in Eq.
(5) we get

$\displaystyle\int\limits^{1-x_{1}}_{0}dx_{2}\,x_{2}\,\mathcal{R}^{\ast}_{j_{1}}\overset{?}{=}1,$
(44)

where, similarly to Eq. (43), we defined

$\displaystyle\mathcal{R}^{\ast}_{j_{1}}\equiv\frac{1}{1-x_{1}}\sum_{j_{2}}\frac{D^{\rm
sym}_{j_{1}j_{2}}(x_{1},x_{2},Q)}{f^{r}_{j_{1}}(x_{1},Q)}.$ (45)

In Table 4 we check Eq. (44) for the symmetrized P  YTHIA dPDFs for the case
where $j_{1}=u$. As a reference we also add the check of Eq. (32) for
asymmetric P  YTHIA dPDFs as well as for the GS09 and naive dPDFs. We see that
if $x_{1}\ =10^{-6}-10^{-1}$ the symmetric, asymmetric, and naive dPDFs
satisfy the GS momentum sum rule at the 2$\%$ accuracy level; the GS09 dPDFs
satisfy the sum rules to a similar extent except at $x_{1}=10^{-6}$, as
discussed in the previous section. We also note that the difference between
the considered models increases with the value of $x_{1}$. For example, if
$x_{1}=0.8$ the asymmetric P  YTHIA dPDFs and symmetric GS09 dPDFs satisfy the
momentum sum rule at the 1$\%$ and 13$\%$ accuracy level correspondingly,
whereas the symmetrized P  YTHIA dPDFs and naive dPDFs are demonstrating
strong violations of momentum conservation.

$x_{1}$ | $j_{1}$ | P  YTHIA dPDFs. | P  YTHIA sym. dPDFs. | GS09 dPDFs. | Naive dPDFs.
---|---|---|---|---|---
$10^{-6}$ | $u$ | 0.996 | 0.979 | 0.931 | 0.996
$10^{-3}$ | $u$ | 0.996 | 0.980 | 1.007 | 0.997
$10^{-1}$ | $u$ | 1.008 | 1.014 | 0.994 | 1.106
0.2 | $u$ | 1.010 | 1.047 | 0.978 | 1.244
0.4 | $u$ | 1.011 | 1.133 | 0.948 | 1.649
0.8 | $u$ | 1.011 | 1.679 | 0.882 | 4.131

Table 4: Test of the momentum sum rule for the symmetrized P

YTHIA

dPDFs.

The additional terms introduced by the symmetrization procedure also lead to
violation of the number sum rule given by Eq. (6). Let us illustrate this
statement by considering $u\bar{u}$ and $\bar{u}u$ flavour combinations as the
most prominent examples. Similarly to Eq. (34) we define a response function
for the symmetrized P  YTHIA dPDFs as

$\displaystyle R^{\ast}_{j_{1}j_{2}}(x_{1},x_{2},Q)\equiv x_{2}\,\frac{D^{\rm
sym}_{j_{1}j_{2}}(x_{1},x_{2},Q)-D^{\rm
sym}_{j_{1}\bar{j_{2}}}(x_{1},x_{2},Q)}{f_{j_{1}}^{r}(x_{1},Q)},$ (46)

which can be written as

$\displaystyle
R^{\ast}_{j_{1}j_{2}}(x_{1},x_{2},Q)=\frac{x_{2}}{2\,f^{r}_{j_{1}}(x_{1})}\left[\tilde{R}_{j_{1}j_{2}}(x_{1},x_{2})-\tilde{R}_{j_{1}\bar{j}_{2}}(x_{1},x_{2})\right],$
(47)

where, omitting factorization scale labels, we defined

$\displaystyle\tilde{R}_{j_{1}j_{2}}(x_{1},x_{2})$ $\displaystyle\equiv$
$\displaystyle f^{r}_{j_{1}}(x_{1})f^{m\leftarrow
j_{1},x_{1}}_{j_{2}}(x_{2})+f^{r}_{j_{2}}(x_{2})f^{m\leftarrow
j_{2},x_{2}}_{j_{1}}(x_{1}).$ (48)

The second term on the right hand side of Eq. (48) is new and leads to the
violation of the number sum rule by the symmetrized P  YTHIA dPDFs. Let us
consider first the $u\bar{u}$ flavour combination. The
$R^{\ast}_{u\bar{u}}(x_{1},x_{2},Q)$ response functions are given in Fig. 4 a)
and the corresponding checks of the number sum rule are provided in Table 5.
The checks for the $\bar{u}u$ flavour combinations are given in Fig. 4 b) and
Table 6. We see that in the $u\bar{u}$ case the deviation of the number sum
rule integral from the expected value of $N_{\bar{u}_{v1}}\equiv-
N_{u_{v1}}=-1$ is about $10-25\%$ depending on the value of $x_{1}$. However,
in the $\bar{u}u$ case the deviations become stronger. In particular, the sum
rule integral shows a deviation from the expected value of $3$ by around
16$\%$ for moderate ($[10^{-3},10^{-1}]$) values of $x_{1}$ and a large
deviation for the valence region (about 130$\%$ deviation if $x_{1}=0.8$).

$x_{1}$ |  $N_{\bar{u}_{v1}}$ P  YTHIA |  $N_{\bar{u}_{v1}}$ P  YTHIA sym. |  $N_{\bar{u}_{v1}}$ GS09 |  $N_{u_{v1}}$ Naive
---|---|---|---|---
$10^{-6}$ | -1.022 | -1.227 | -1.526 | -2.006
$10^{-3}$ | -1.007 | -0.847 | -1.000 | -2.006
$10^{-1}$ | -1.004 | -0.925 | -1.036 | -2.005
0.2 | -1.003 | -0.928 | -1.016 | -2.005
0.4 | -1.002 | -0.884 | -0.949 | -1.997
0.8 | -0.998 | -0.740 | -0.862 | -1.708

Table 5: Numerical integration with respect to $x_{2}$ over
$R^{\ast}_{u\bar{u}}(x_{1},x_{2},Q)$ response function at fixed $x_{1}$. In
the ideal situation when the GS sum rules are perfectly satisfied
$N_{\bar{u}_{v1}}=-1$.

$x_{1}$ |  $N_{u_{v1}}$ P  YTHIA |  $N_{u_{v1}}$ P  YTHIA sym. |  $N_{u_{v1}}$ GS09 |  $N_{u_{v1}}$ Naive
---|---|---|---|---
$10^{-6}$ | 2.989 | 2.961 | 2.163 | 2.006
$10^{-3}$ | 3.003 | 3.351 | 2.982 | 2.006
$10^{-1}$ | 3.002 | 3.491 | 2.968 | 2.005
0.2 | 3.002 | 3.580 | 2.912 | 2.005
0.4 | 3.000 | 3.858 | 2.833 | 1.997
0.8 | 2.994 | 7.048 | 2.667 | 1.708

Table 6: Numerical integration with respect to $x_{2}$ over
$R_{\bar{u}u}(x_{1},x_{2},Q)$ response function at fixed $x_{1}$. In the ideal
situation when the GS sum rules are perfectly satisfied $N_{u_{v1}}=3$.

a)

b)

Figure 4: Same as in Fig. 2 but for the $u\bar{u}$ and $\bar{u}u$ flavour
combinations: a) $R^{\ast}_{u\bar{u}}$ response function, b)
$R^{\ast}_{\bar{u}u}$ response function.

In order to find out the origin of the peak of $R^{\ast}_{\bar{u}u}$ at
$x_{1}=0.8$ that is responsible for this large deviation, we study in more
detail the terms on the right hand side of Eq. (47). Let us start with the
$u\bar{u}$ flavour combination. The corresponding terms are given by

$\displaystyle\frac{x_{2}}{2\,f^{r}_{u}(x_{1})}\,\tilde{R}_{u\bar{u}}(x_{1},x_{2})$
$\displaystyle=$ $\displaystyle\frac{x_{2}}{2}\left[f^{m\leftarrow
u,x_{1}}_{\bar{u}}(x_{2})+\frac{f^{r}_{\bar{u}}(x_{2})f^{m\leftarrow\bar{u},x_{2}}_{u}(x_{1})}{f^{r}_{u}(x_{1})}\right],$
(49)
$\displaystyle\frac{x_{2}}{2\,f^{r}_{u}(x_{1})}\,\tilde{R}_{uu}(x_{1},x_{2})$
$\displaystyle=$ $\displaystyle\frac{x_{2}}{2}\left[f^{m\leftarrow
u,x_{1}}_{u}(x_{2})+\frac{f^{r}_{u}(x_{2})f^{m\leftarrow
u,x_{2}}_{u}(x_{1})}{f^{r}_{u}(x_{1})}\right].$ (50)

In Fig. 5 a) we plot the four different terms given in Eqs. (49) – (50) (red
and blue curves). The difference between Eq. (49) and Eq. (50) (which gives us
$R^{\ast}_{u\bar{u}}$) is given by the green dash-dotted line. We see that all
four terms have rather similar behaviour. Now let us consider the $\bar{u}u$
case which yields

$\displaystyle\frac{x_{2}}{2\,f^{r}_{\bar{u}}(x_{1})}\,\tilde{R}_{\bar{u}u}(x_{1},x_{2})$
$\displaystyle=$
$\displaystyle\frac{x_{2}}{2}\left[f^{m\leftarrow\bar{u},x_{1}}_{u}(x_{2})+\frac{f^{r}_{u}(x_{2})f^{m\leftarrow
u,x_{2}}_{\bar{u}}(x_{1})}{f^{r}_{\bar{u}}(x_{1})}\right],$ (51)
$\displaystyle\frac{x_{2}}{2\,f^{r}_{\bar{u}}(x_{1})}\,\tilde{R}_{\bar{u}\bar{u}}(x_{1},x_{2})$
$\displaystyle=$
$\displaystyle\frac{x_{2}}{2}\left[f^{m\leftarrow\bar{u},x_{1}}_{\bar{u}}(x_{2})+\frac{f^{r}_{\bar{u}}(x_{2})f^{m\leftarrow\bar{u},x_{2}}_{\bar{u}}(x_{1})}{f^{r}_{\bar{u}}(x_{1})}\right].$
(52)

a)

b)

Figure 5: Different terms contributing to a)
$R^{\ast}_{u\bar{u}}(x_{1},x_{2})$ and b) $R^{\ast}_{\bar{u}u}(x_{1},x_{2})$.

a)

b)

Figure 6: Decomposition into factors of a)
$f^{r}_{\bar{u}}(x_{2})f^{m\leftarrow\bar{u},x_{2}}_{u}(x_{1})/f^{r}_{u}(x_{1})$
and b) $f^{r}_{u}(x_{2})f^{m\leftarrow
u,x_{2}}_{\bar{u}}(x_{1})/f^{r}_{\bar{u}}(x_{1})$ terms.

In Fig. 5 b) we plot the four different terms given in Eqs. (51)–(52). Unlike
the $u\bar{u}$ case we see that now different terms have a different
behaviour. Namely, we see that the term proportional to
$f^{r}_{u}(x_{2})f^{m\leftarrow
u,x_{2}}_{\bar{u}}(x_{1})/f^{r}_{\bar{u}}(x_{1})$ (red dotted line) has the
biggest contribution and is largely responsible for the large peak given by
the green dash-dotted line. Let us study the aforementioned ratio in more
details. In Fig. 6 a) we plot a decomposition of the
$f^{r}_{\bar{u}}(x_{2})f^{m\leftarrow\bar{u},x_{2}}_{u}(x_{1})/f^{r}_{u}(x_{1})$
term and in Fig. 6 b) we plot its analogue for the $\bar{u}u$ flavour
combination $f^{r}_{u}(x_{2})f^{m\leftarrow
u,x_{2}}_{\bar{u}}(x_{1})/f^{r}_{\bar{u}}(x_{1})$. By comparing the middle
plots in Fig. 6 we see that $f^{m\leftarrow\bar{u},x_{2}}_{u}(x_{1})$ rapidly
drops down in the region between $10^{-2}$ and $10^{-1}$ whereas
$f^{m\leftarrow u,x_{2}}_{\bar{u}}(x_{1})$ becomes peaked. The reason for such
a different behaviour is as follows. In the case of
$f^{m\leftarrow\bar{u},x_{2}}_{u}(x_{1}=0.8)$ we are looking at how finding a
$\bar{u}$ quark with momentum $x_{2}$ affects the distribution of $u$ quarks
with momentum $x_{1}=0.8$. Since there is a sizeable distribution of $u$
quarks at large $x$ (arising from the valence contribution), the boost to this
distribution from the companion quark effects is negligible, and the only
visible effect from the $\bar{u}$ is just the “squeezing” caused by momentum
conservation. That is, according to Section 2.2, for the middle plot in Fig. 6
a) at a given factorization scale we have

$\displaystyle
f^{m\leftarrow\bar{u},x_{2}}_{u}(x_{1})\approx\frac{1}{1-x_{2}}f^{r}_{u}(\frac{x_{1}}{1-x_{2}}).$
(53)

However, in case of the middle plot in Fig. 6 b) the situation changes. The
$\bar{u}$-quark sPDF, having no valence component, is much smaller at large
$x$ than the $u$-quark sPDF999At $Q=91$ GeV, we have $f^{r}_{u}(x_{1}=0.8)\sim
10^{-3}$ whereas $f^{r}_{\bar{u}}(x_{1}=0.8)\sim 10^{-5}$. . In this case, the
companion quark effects are noticeable, particularly when $x_{2}$ is close in
magnitude with $x_{1}$. This results in a peak in the middle plot in Fig. 6
b). Additionally, the comparison between upper plots in Fig. 6 shows us that
the $x_{2}f^{r}_{u}(x_{2})$ term stays roughly the same in the region of
interest (changes from 0.75 to 0.5) whereas $x_{2}f^{r}_{\bar{u}}(x_{2})$
drops from 0.6 down to zero. Combining upper and middle plots in Fig. 6
together and dividing by $2\,f^{r}_{u}(x_{1}=0.8)$ or
$2\,f^{r}_{\bar{u}}(x_{1}=0.8)$ we get a smoothly decreasing function in the
lower plot in Fig. 6 a) and a peaked function in the lower plot in Fig. 6 b)
correspondingly. Therefore, we conclude that the large peak in the response
functions for the symmetrized dPDFs shown in Fig. 4 b) is due to the companion
quark contribution to the $\bar{u}$ sPDF $f^{m\leftarrow
u,x_{2}}_{\bar{u}}(x_{1})$. We also shall note that similar large peaks occur
for other flavour combinations where a sea quark sPDF receives a large
contribution from a $g\rightarrow q\bar{q}$ splitting, e.g. for the $s\bar{s}$
flavour combination as shown in Fig. 7.

Figure 7: Same as in Fig. 4 but for the $R^{\ast}_{s\bar{s}}$ response
function.

Thus, we see that the symmetric dPDFs constructed using the MPI model of P
YTHIA do not satisfy the sum rules perfectly. Nevertheless, we shall note that
the largest violation of the GS sum rules by dPDFs constructed as in Eq. (30)
occur at very large values of Bjorken-$x$’es (if $x\geq 0.4$) where the
absolute values of dPDFs are strongly suppressed. For the rest of the
considered Bjorken-$x$’es the symmetric dPDFs obey the GS sum rules at
$10-25\%$ accuracy level. It would be interesting to modify the P  YTHIA
algorithm to suppress the large peaks in Figs. 4 and 7 such that we obtain
symmetric dPDFs which obey the GS sum rules at a better accuracy level (for
example, by suppressing the probability to find a companion quark accompanying
a quark, when both partons have large $x$ fractions). However, such a study is
beyond the scope of this paper.

###  3.3 Double Drell-Yan production

In Section 3.1 and Section 3.2 we studied how well asymmetric and symmetric P
YTHIA dPDFs satisfy the GS sum rules, and compared these to the GS09 dPDFs.
Notably, we saw that the response functions defined by Eqs. (34) and (43)
match rather closely between the GS09 and P  YTHIA dPDFs.

Given this close match, an interesting question is to what extent the P  YTHIA
and GS09 dPDFs are similar overall. Here we study this question via
phenomenological MC simulations using both sets of dPDFs. We also include
predictions obtained using the naive dPDFs defined in Eq. (39), as a
benchmark. It is most interesting to study a DPS process where the effects of
the GS sum rules can potentially play an important role - i.e. one which
receives contributions associated with two quarks of the same (anti)flavour
belonging to the same hadron. Therefore, we concentrate on a four-lepton
production through the double Drell-Yan (dDY) process Goebel:1979mi ;
Mekhfi:1983az ; Halzen:1986ue ; Kom:2011nu ; Myska:2011ji ; Kasemets:2012pr ;
Krasny:2013aca , defined as
$\left(q\bar{q}\rightarrow\gamma^{\ast}/Z^{\ast}\rightarrow
2l\right)\otimes\left(q\bar{q}\rightarrow\gamma^{\ast}/Z^{\ast}\rightarrow
2l\right)$, where $l$ labels either an electron or muon final state.

Before presenting our results, let us note that the older versions of P  YTHIA
(version $\leq$ 8.235) were using asymmetric DPS luminosities

$\displaystyle\mathcal{L}_{j_{1},j_{2},j_{3},j_{4}}(x_{1},x_{2},x_{3},x_{4})=D_{j_{1},j_{2}}(x_{1},x_{2})\,D_{j_{3},j_{4}}(x_{3},x_{4})$
(54)

to simulate DPS production, whereas starting from version 8.240 P  YTHIA uses
the symmetric ones101010More precisely, staring from P  YTHIA version 8.240
Eq. (55) is used when the option ``SecondHard— is switched on. However, if the
event is produced through the standard MPI machinery then a simple product of
sPDFs as in Eq. (54) is used.

$\displaystyle\mathcal{L}^{\rm
sym}_{j_{1},j_{2},j_{3},j_{4}}(x_{1},x_{2},x_{3},x_{4})$
$\displaystyle=\frac{1}{2}\,\mathcal{L}_{j_{1},j_{2},j_{3},j_{4}}(x_{1},x_{2},x_{3},x_{4})$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{2}\,\mathcal{L}_{j_{2},j_{1},j_{4},j_{3}}(x_{2},x_{1},x_{4},x_{3}),$
(55)

where for the sake of simplicity we skipped factorization scale labels.
However, if the number of generated events is high enough then MC simulations
with Eq. (54) and Eq. (55) lead to the same results, since both equivalent
ways to produce the $A\,B$ final state
($j_{1}\,j_{3}\,\rightarrow\,A\,\otimes\,j_{2}\,j_{4}\,\rightarrow\,B$ and
$j_{2}\,j_{4}\,\rightarrow\,A\,\otimes\,j_{1}\,j_{3}\,\rightarrow\,B$) will be
used approximately equal number of times. To check this statement we simulated
the DPS distributions presented in this section with P  YTHIA 8.235 and P
YTHIA 8.240, and found basically no difference between both programs.
Therefore, throughout this section we use the symmetrized DPS luminosities Eq.
(55) only and omit the “$\rm sym$” labels to keep our notation
simple111111Bear in mind that the usage of the symmetic DPS luminosities as in
Eq. (55) does not correspond to working with symmetric dPDFs. The luminosities
obtained using symmetric dPDFs are given by:  $\displaystyle D^{\rm
sym}_{j_{1},j_{2}}(x_{1},x_{2})\,D^{\rm sym}_{j_{3},j_{4}}(x_{3},x_{4})$
$\displaystyle=\frac{1}{2}\,\mathcal{L}^{\rm
sym}_{j_{1},j_{2},j_{3},j_{4}}(x_{1},x_{2},x_{3},x_{4})$
$\displaystyle\mathrel{\phantom{=}}+\frac{1}{2}\,\mathcal{L}^{\rm
sym}_{j_{1},j_{2},j_{4},j_{3}}(x_{1},x_{2},x_{4},x_{3}).$ (56) By comparing
Eq. (56) with Eq. (55) we see that the symmetrization procedure as in P  YTHIA
8.240 does not take into account the second term on the right hand side of Eq.
(56)..

To obtain the predictions with the GS09 and naive dPDFs, we change the weight
of each P  YTHIA event by multiplying it by a corresponding ratio of DPS
luminosities:121212To perform the reweighting procedure one has to either use
the method we describe in Appendix A to compute $\mathcal{L}^{\rm Pythia}$ or
to be able to access DPS luminosities used by P  YTHIA for a given DPS event.
The latter can be achieved, for example, by adding new members to the ``Event—
class to store the information about the sPDFs being used to simulate the
second interaction. This additional information has to be added to the event
record in the ``ProcessLevel::nextOne(Event& process)— method of the
``ProcessLevel— class after the successful generation of the second hard
interaction. The information assigned in such a way is then accessible
together with the other standard information stored in the event record and
can be used to evaluate $\omega^{\rm GS09}$ and $\omega^{\rm Naive}$ during
the generation procedure.

$\displaystyle\omega^{\rm GS09}$ $\displaystyle=$ $\displaystyle\omega^{\rm
Pythia}\,\frac{\mathcal{L}^{\rm GS09}}{\mathcal{L}^{\rm Pythia}},$ (57)
$\displaystyle\omega^{\rm Naive}$ $\displaystyle=$ $\displaystyle\omega^{\rm
Pythia}\,\frac{\mathcal{L}^{\rm Naive}}{\mathcal{L}^{\rm Pythia}}.$ (58)

Before discussing our results we would also like to note that according to the
P  YTHIA model the total DPS cross section is given by

$\displaystyle\sigma_{\rm DPS}=\frac{\sigma_{\rm
ND}}{2}\left(\frac{\sigma_{1}}{\sigma_{\rm ND}}\right)^{2},$ (59)

where $\sigma_{\rm ND}$ is a total non-diffractive cross section
Sjostrand:2004pf . This differs by a factor of $\sigma_{\rm eff}/\sigma_{\rm
ND}$ from the usual “pocket formula” result: $\sigma_{\rm
DPS}=\sigma^{2}_{1}/2\sigma_{\rm eff}$. The transition from Eq. (59) to the
pocket formula expression in P  YTHIA , however, is used only in the MPI
machinery and is not performed if the flag `SecondHard` is turned on. This
implies that the DPS cross sections in our simulations are evaluated according
to Eq. (59) with $\sigma_{\rm ND}$ instead of $\sigma_{\rm eff}$. However,
since we are interested only in the relative difference between distributions
of DPS events produced with different models of dPDFs, the aforementioned
difference in the value of the DPS cross section does not affect our analysis.

The differential distributions generated with the naive, P  YTHIA and GS09
dPDFs are given in Fig. 8. We perform binning in terms of the maximal rapidity
difference ${\rm\Delta Y}\equiv{\rm max}|y_{i}-y_{j}|$ of the produced leptons
Kom:2011bd . The absolute values of rapidities of all produced leptons were
requested to be $\leq 5$. The complete P  YTHIA setup we use is given in
Appendix B. Since we are essentially interested in comparison between the
different models of dPDFs we do not impose any additional cuts on the final
state leptons. Instead, in order to have better control of the probed values
of the Bjorken-$x$’es, we consider different values of the collision energy,
namely 3, 7, 13 and 20 TeV. We fix the renormalization scale equal to the
factorization scale, and set the factorization scale of both interactions to
91 GeV ( $\sim M_{Z}$).

Figure 8: Comparison between ${\rm\Delta Y}={\rm max}|y_{i}-y_{j}|$ DPS
distributions generated with the GS09, P

YTHIA

and naive dPDFs. The collision energy $\sqrt{S}$ is equal to $3$, $7$, $13$
and $20$ TeV.

We see that the ratio $d\sigma_{\rm{Naive}}/d\sigma_{\rm{Pythia}}$ is fairly
flat in ${\rm\Delta Y}$, and remains rather close to $1$ for the whole range
of ${\rm\Delta Y}$ (except at the 3 TeV collision energy where we see some
increase at small and large values of ${\rm\Delta Y}$), whereas
$d\sigma_{\rm{GS09}}/d\sigma_{\rm{Pythia}}$ decreases as one approaches the
largest values of ${\rm\Delta Y}$. We see also that the ratio
$d\sigma_{\rm{GS09}}/d\sigma_{\rm{Pythia}}$ remains more or less constant at
small and moderate ${\rm\Delta Y}$ as $\sqrt{S}$ is varied, whilst there is a
noticeable decrease in this ratio at large ${\rm\Delta Y}$ as $\sqrt{S}$
decreases.

Let us now try to explain all of these features, starting with the very large
${\rm\Delta Y}$ region (and in particular the last bin with ${\rm\Delta
Y}\in[9.5,10]$). We expect this ${\rm\Delta Y}$ bin to probe a very asymmetric
configuration of $x$ fractions in the dPDF, with the values of $x$ increasing
as $\sqrt{S}$ is lowered. This is confirmed in the plots of Fig. 9, which also
shows the values of $x$ probed at each value of $\sqrt{S}$.

a)

b)

Figure 9: Distribution of Bjorken $x$-es belonging to the same hadron for the
bin with ${\rm\Delta Y}\in[9.5,10]$. Left plot: distribution of the minimal
and maximal Bjorken $x$-es at 3, 7, 13 and 20 TeV collision energies. Right
plot: distribution of the DPS events contributing to the ${\rm\Delta
Y}\in[9.5,10]$ bin in the $x_{1}$ \- $x_{2}$ plane for the $x_{2}<x_{1}$ case.

Why does the ratio $d\sigma_{\rm{GS09}}/d\sigma_{\rm{Pythia}}$ decrease as we
decrease $\sqrt{S}$ and probe larger $x$ values, whilst
$d\sigma_{\rm{Naive}}/d\sigma_{\rm{Pythia}}$ stays fairly close to $1$? To try
to understand this, we use a simplified setup, in which we consider only one
partonic contribution; the one in which a $u$ and a $\bar{u}$ from one proton
collide with a $\bar{u}$ and a $u$ in the other. This sort of contribution
should be dominant at large $\Delta{\rm Y}$, as we can have a right-moving
valence quark in the first subprocess, and a left-moving valence quark in the
second. We fix the $x$ values in the dPDFs to correspond to the production of
two heavy states of mass $M=\sqrt{x_{1}\,x_{2}\,S}=13$ GeV, with the first
being produced with rapidity ${\rm Y}=\log(x_{1}/x_{2})/2=4.5$, and the other
being produced with rapidity ${\rm Y}=-4.5$ (yielding ${\rm\Delta Y}=9$
overall), and for simplicity we consider the luminosity rather than the cross
section. The corresponding results are given in Fig. 10 where we plot ratios
of DPS luminosities as a function of collision energy $\sqrt{S}$. As we
mentioned before, Fig. 8 is generated using the particular choice of
factorization scale $Q=91$ GeV. However, in Fig. 10 we choose to generate the
ratios for a range of factorization scales in between $1$ and $180$ GeV, as
this helps us explain the observed behaviour at $Q=91$ GeV.

Figure 10: Ratios of the DPS luminosities at different values of factorization
scale $Q$ as functions of collision energy $\sqrt{S}$ for the $u$-$\bar{u}$
flavour combination.

Let us first study the $Q$-dependence of the $\mathcal{L}^{\rm
Pythia}/\mathcal{L}^{\rm Naive}$ ratio which is driven by the modification of
the sPDFs affecting the $\mathcal{L}^{\rm Pythia}$ luminosity. We see in the
upper plot of Fig. 10 that at $Q=91$ GeV the ratio $\mathcal{L}^{\rm
Pythia}/\mathcal{L}^{\rm Naive}$ changes if $\sqrt{S}\in[3,9]$ TeV and becomes
almost constant if $\sqrt{S}>9$ TeV, which is consistent with our results
shown in Fig. 8131313Note that here we consider a given flavour combination
and fixed $x$-values. Therefore, we expect only qualitative agreement between
Figs. 8 and 10.. We also see that the ratio $\mathcal{L}^{\rm
Pythia}/\mathcal{L}^{\rm Naive}$ demonstrates a considerable change as the
factorization scale increases from $1$ to $180$ GeV (compare orange and green
curves). The $Q$-dependence becomes milder as $Q$ increases, and is very
gradual for $Q\geq 45$ GeV. To investigate this behaviour let us note that for
a given configuration ($x_{1}=x_{3},\,x_{2}=x_{4}$,
$j_{1}=j_{3},\,j_{2}=j_{4}$) $\mathcal{L}^{\rm Pythia}/\mathcal{L}^{\rm
Naive}$ is given by

$\frac{D^{\rm Pythia}_{u\bar{u}}(x_{1},x_{2})\,D^{\rm
Pythia}_{\bar{u}u}(x_{2},x_{1})}{D^{\rm Naive}_{u\bar{u}}(x_{1},x_{2})\,D^{\rm
Naive}_{\bar{u}u}(x_{2},x_{1})}=\left[\frac{f^{m\leftarrow
u,x_{1}}_{\bar{u}}(x_{2},Q)}{f^{r}_{\bar{u}}(x_{2},Q)}\right]\left[\frac{f^{m\leftarrow\bar{u},x_{2}}_{u}(x_{1},Q)}{f^{r}_{u}(x_{1},Q)}\right].$
(60)

The kinematics we consider implies that $x_{1}\in[10^{-2},1]$ whereas $x_{2}$
is well between $10^{-4}$ and $10^{-6}$ (for instance, if $\sqrt{S}=3$ TeV
$x_{1}\approx 0.4$ and $x_{2}\approx 4.81\times 10^{-5}$). We found that for
the Bjorken-$x$’es corresponding to Fig. 10 the second factor on the right
hand side of Eq. (60) is very close to unity whereas the first one is
responsible for the observed differences. Since the first factor on the right
hand side of Eq. (60) implies that in the first interaction we probe a valence
$u$-quark and a sea $\bar{u}$-quark in the second one, then, according to the
P  YTHIA model, we deal with two modification mechanisms: first we squeeze all
sPDFs according to Eqs. (13), (17) and then rescale the sea sPDFs by
multiplying them by the $a$-factor as in Eqs. (27), (28) to preserve overall
momentum conservation. We found that the rescaling of the $\bar{u}$ sPDFs
according to Eq. (27) increases the value of $f^{m\leftarrow
u,x_{1}}_{\bar{u}}$ whereas the squeezing according to Eq. (13) tends to make
it smaller. The effect of the squeezing depends on the rate at which the
$\bar{u}$ PDF grows at small $x$ – if the $\bar{u}$ PDF at small $x$ is
proportional to $x^{-b}$, then taking only the squeezing effect into account
we would have

$\displaystyle\frac{f^{m\leftarrow
u,x_{1}}_{\bar{u}}(x_{2},Q)}{f^{r}_{\bar{u}}(x_{2},Q)}=(1-x_{1})^{b-1}.$ (61)

For the MSTW set we use, $b>1$ and grows rapidly at small $Q\sim 1-2$ GeV
before reaching an approximately stable value of $b\sim 1.5$ at larger $Q$
values. Thus, the squeezing effect corresponds to a suppression of the ratio
$f^{m\leftarrow u,x_{1}}_{\bar{u}}/f^{r}_{\bar{u}}$, which increases at low
$Q$ but stabilises at larger $Q$ values. This suppression also increases as
$x_{1}$ increases, or $\sqrt{S}$ decreases.

The “$a$-rescaling” effect is initially quite large – we find $a\sim 1.2-1.3$
at $Q=1$ GeV, which explains the (perhaps unexpected!) fact that
$\mathcal{L}^{\rm Pythia}/\mathcal{L}^{\rm Naive}$ significantly exceeds $1$
at this $Q$ value. This enhancement rapidly reduces to $\sim 1.1$ at $Q\sim 5$
GeV, and thereafter changes rather slowly. At $Q=1$ GeV, there is a noticeable
decrease in $a$ as $\sqrt{S}$ increases, but at higher $Q$ values the
behaviour of $a$ is rather flat in $\sqrt{S}$.

Thus, we see that the reduction of $\mathcal{L}^{\rm Pythia}/\mathcal{L}^{\rm
Naive}$ as $Q$ increases is caused both by the squeezing suppression getting
stronger and the $a$-rescaling enhancement getting weaker. At large $Q$ values
such as $Q=91$ GeV the downward trend in the ratio as $\sqrt{S}$ decreases is
mainly driven by the squeezing effect getting stronger, as the $a$-rescaling
enhancement is roughly constant with $\sqrt{S}$.

Now let us consider the ratio $\mathcal{L}^{\rm GS09}/\mathcal{L}^{\rm Naive}$
shown in the middle plot of Fig. 10141414Note that the oscillation of the
ratios involving the $\mathcal{L}^{\rm GS09}$ luminosities in Fig. 10 are
numerical artifacts of the linear interpolation used in the GS09 grids for the
Bjorken-$x$’es $\geq 0.1$. . At $Q=1$ GeV this ratio is very close to unity,
which is essentially by design – when the $x$ value of one parton is very
small (as we have here), the GS09 inputs are constructed to match closely with
the naive dPDFs, as this will ensure the momentum sum rule is well-satisfied
when integrating over the other parton (c.f. Eq. (42)). At $Q=1$ GeV the
differences between the GS09 and naive dPDFs predominantly manifest themselves
in the region where both $x$ values in the dPDF are large (see Ref.
Gaunt:2009re ). However, as $Q$ starts to increase we observe some important
deviations of the ratio $\mathcal{L}^{\rm GS09}/\mathcal{L}^{\rm Naive}$ from
1 in Fig. 10 caused by the dDGLAP evolution. We see that if $\sqrt{S}$ is
approximately in the $[3,10]$ TeV interval the ratio $\mathcal{L}^{\rm
GS09}/\mathcal{L}^{\rm Naive}$ is smaller than one, whereas at large
$\sqrt{S}$ we observe the opposite behaviour. This observation is consistent
with the picture of the dDGLAP evolution from Ref. Gaunt:2009re ; at the
starting value of the evolution scale $Q=1$ GeV the differences between the
GS09 and naive dPDFs are located in the region where both $x_{i}$ are large,
however, as $Q$ grows the differences propagate down to smaller values of
$x_{1},x_{2}$ leading to a decrease of the GS09 $D_{u\bar{u}}$ dPDFs relative
to the naive ones. Taking into account the fact that in our analysis the
decrease in $\sqrt{S}$ leads to an increase of both $x_{1}$ and $x_{2}$ (see
the distributions in Fig. 9) we conclude that if $\sqrt{S}\in[3,10]$ TeV the
$x_{i}$ values enter the region significantly affected by the aforementioned
propagation of the differences between two models down to the smaller values
of momentum fractions. However, as $\sqrt{S}$ increases both $x_{1}$ and
$x_{2}$ become smaller which lead us to the region of small $x$’es where the
$1\rightarrow 2$ splitting contributions due to the last term on the right
hand side of Eq. (3) become important. Since the $1\rightarrow 2$ feeding term
causes an increase of the parton densities we conclude that the change of the
behaviour of the $\mathcal{L}^{\rm GS09}/\mathcal{L}^{\rm Naive}$ ratio at
large values of $\sqrt{S}$ is due to the positive contribution from
$1\rightarrow 2$ splittings.

Having now studied and understood the behaviour of the $\mathcal{L}^{\rm
Pythia}/\mathcal{L}^{\rm Naive}$ and $\mathcal{L}^{\rm GS09}/\mathcal{L}^{\rm
Naive}$ ratios, we can now understand the behaviour of $\mathcal{L}^{\rm
GS09}/\mathcal{L}^{\rm Pythia}$ as this can simply be constructed by taking
the quotient of the ratios already studied. In particular, we note that at
very low $Q$’s the ratio $\mathcal{L}^{\rm GS09}/\mathcal{L}^{\rm Pythia}$
lies significantly below $1$ predominantly due to the $a$-rescaling
enhancement in P  YTHIA . At large $Q$, the trend in $\sqrt{S}$ is mainly
driven by momentum/number sum rule suppression in the GS09 dPDFs, which has
started to propagate downwards into the region where one $x_{i}$ is large and
the other small. As one can see the behaviour of the $\mathcal{L}^{\rm
GS09}/\mathcal{L}^{\rm Pythia}$ at $Q=91$ GeV is consistent with the
distributions shown in Fig. 8.

Now let us focus on the small $\Delta\rm Y$ region in Fig. 8, and initially
restrict our attention to the behaviour of
$d\sigma_{\rm{Naive}}/d\sigma_{\rm{Pythia}}$. We saw in Fig. 8 that this ratio
stays within 5% of $1$ for the whole range of $\Delta\rm Y$ at 7, 13 and 20
TeV collision energies. However, if $\sqrt{S}=3$ TeV we observe about 10%
deviation of the aforementioned ratio from $1$ at small $\Delta\rm
Y\in[0,0.5]$.

The $\Delta\rm Y\in[0,0.5]$ bin should probe dPDFs where the two $x_{i}$
values are of the same magnitude – this is confirmed in Fig. 11 where we plot
the distribution of Bjorken-$x$’es for the P  YTHIA events contributing to
this bin. We also see that the distribution of events in this common $x$ value
is approximately flat with the bulk of events having $x\in[10^{-5},10^{-1}]$.
Furthermore, unlike the $\Delta\rm Y\in[9.5,10]$ bin which predominantly
probes dPDFs with flavour structure $q_{i}\bar{q_{j}}$
($i,j=\\{u,d,s,c,b\\}$), the events in the $\Delta\rm Y\in[0,0.5]$ bin arise
from an approximately even mix of events where a $q_{i}q_{j}$ in one proton
collide with a $\bar{q}_{i}\bar{q}_{j}$ in the other, and events where we have
a $q_{i}\bar{q}_{j}$ vs. $\bar{q}_{i}q_{j}$ collision – see Table 7. Finally,
we found that amongst the five (anti-)quark flavours we consider, the
combinations of $\bar{u}$, $\bar{d}$, $u$ and $d$ flavours occurs most
frequently.

a)

b)

Figure 11: Distribution of Bjorken-$x$’es belonging to the proton with the $D_{j_{1}j_{2}}(x_{1},x_{2},Q)$ dPDF for the bin with ${\rm\Delta Y}\in[0,0.5]$. Left plot: distribution of the minimal and maximal Bjorken-$x$’es for 3, 7, 13 and 20 TeV collision energies. Right plot: distribution of the DPS events contributing to the ${\rm\Delta Y}\in[0,0.5]$ bin in the $x_{1}$ \- $x_{2}$ plane for the $x_{2}<x_{1}$ case. $\sqrt{S}$ | $\Delta\rm Y\in[0,0.5]$ | $\Delta\rm Y\in[9.5,10]$
---|---|---
$\sqrt{S}=3$ TeV, $q_{i}\,q_{j}$ | 53.3% | 19.2%
$\sqrt{S}=7$ TeV, $q_{i}\,q_{j}$ | 51.3% | 29.9%
$\sqrt{S}=13$ TeV, $q_{i}\,q_{j}$ | 50.7% | 38.5%
$\sqrt{S}=20$ TeV, $q_{i}\,q_{j}$ | 50.1% | 42.7%
Table 7: Flavour composition of the dPDFs belonging to the same proton and
contributing to the first and last $\Delta\rm Y$ bins
($q_{i}\,q_{j}\,+\,q_{i}\,\bar{q}_{j}=100\%$).

In order to get a better understanding of
$d\sigma_{\rm{Naive}}/d\sigma_{\rm{Pythia}}$ in the region $\Delta\rm
Y\in[0,0.5]$, we plot the ratio $\mathcal{L}^{\rm Naive}/\mathcal{L}^{\rm
Pythia}$ where we set all Bjorken-$x$’es to the same value in Fig. 12. In this
figure we consider combinations of $u$ and $d$ flavours having different
behaviour, namely $\mathcal{L}_{uu\bar{u}\bar{u}}$,
$\mathcal{L}_{u\bar{u}\bar{u}u}$, $\mathcal{L}_{ud\bar{u}\bar{d}}$ and
$\mathcal{L}_{u\bar{d}\bar{u}d}$ (here we skip ${dd\bar{d}\bar{d}}$ and
${d\bar{d}\bar{d}d}$ flavour combinations since they have behaviour similar to
the analogous combinations of $u$ and $\bar{u}$ quarks). We see that for all
considered flavour combinations, except for the $u\bar{u}\bar{u}u$ case, the
ratio is bigger than unity for the whole range of $x$. The decrease of the
ratio $\mathcal{L}^{\rm Naive}/\mathcal{L}^{\rm Pythia}$ for the
$u\bar{u}\bar{u}u$ flavour combination below $1$ is due to the companion quark
contribution from $g\rightarrow u\bar{u}$ splitting in the P  YTHIA case. We
also see that the ratio corresponding to the $uu\bar{u}\bar{u}$ flavour
combination increases faster compared to the $ud\bar{u}\bar{d}$ and
$u\bar{d}\bar{u}d$ cases. This behaviour can be explained by the fact that for
the $uu$ dPDF we invoke not only momentum squeezing mechanism in P  YTHIA , as
in Eq. (10), but also have a suppression effect due to the change in remaining
number of valence quarks, as in Eq. (17). It is also interesting to note that
the effects due to the momentum and number sum rule constraints can penetrate
down to rather small Bjorken-$x$ values $\sim 10^{-2}-10^{-1}$. Since at
$\sqrt{S}=3$ TeV collision energy we probe Bjorken-$x$’es at larger values
than at 7 - 20 TeV (see blue distributions in Fig. 11 a)), we conclude that
the increase in the ratio of naive to P  YTHIA distributions at small values
of $\Delta\rm Y$ at 3 TeV collision energy, as shown in Fig. 8, is due to the
number and momentum sum rule effects which manifest themselves at the larger
values of Bjorken-$x$’es accessible at the 3 TeV collision energy.

Figure 12: Ratio of the DPS luminosities
$\mathcal{L}^{\mathrm{Naive}}/\mathcal{L}^{\mathrm{Pythia}}$ as a function of
$x$ where $x_{1}=x_{2}=x_{3}=x_{4}=x$.

Let us finally turn to considering the ratio
$d\sigma_{\rm{GS09}}/d\sigma_{\rm{Pythia}}$ at small to moderate values of
$\rm\Delta Y\in[0,5]$. We see in Fig. 8 that this ratio is approximately
constant for $\rm\Delta Y\in[0,5]$, and is $\sim 1.1$ across all values of
$\sqrt{S}$. This enhancement of GS09 over P  YTHIA is linked to $1\to 2$
splitting contributions to the GS09 dPDFs during dDGLAP evolution, see the the
last term on the right hand side of Eq. (3). To illustrate this, in Fig. 13 we
plot the ratio of the $\Delta\rm Y$ distributions obtained by solving dDGLAP
evolution equation with and without the $1\rightarrow 2$ splitting
contribution. As one can see the evolution effects due to the $1\rightarrow 2$
splittings increase the $\Delta\rm Y$ distributions reasonably uniformly in
the whole range of binning but somewhat more at the small values of $\Delta\rm
Y$ than at the large values of $\Delta\rm Y$. One notes that the enhancement
factor when including the $1\to 2$ splitting effects is $\sim 1.30-1.35$ at
small $\rm\Delta Y$, which is rather larger than the factor $\sim 1.1$
enhancement of GS09 over P  YTHIA that we see in Fig. 8. This implies that
there is some suppression in GS09 from momentum/number sum rule suppression
effects migrating to lower $x$ values during evolution, which partially
compensates the $1\to 2$ splitting enhancement.

Figure 13: Comparison between ${\rm\Delta Y}={\rm max}|y_{i}-y_{j}|$ DPS
distributions generated with the GS09 dPDFs produced by solving dDGLAP
evolution equation with and without $1\rightarrow 2$ splitting. The collision
energy $\sqrt{S}$ is equal to $3$, $7$, $13$ and $20$ TeV.

To finish this section, we would like to comment that despite the similarities
between the response functions computed with P  YTHIA and GS09 dPDFs, as shown
in Figs. 2 and 3, the $\Delta\rm Y$ distributions generated by these models of
dPDFs differ from one another, especially at small and large values of
rapidity separation. This implies some important differences between the dPDF
sets themselves. We also note that the $\Delta\rm Y$ distributions generated
with naive and P  YTHIA dPDFs also show some important differences not only at
large values of rapidity separation (something one would expect, since large
values of $\Delta\rm Y$ allow to probe dPDFs with one large $x$ value, where
momentum sum rule effects play a significant role) but also at small values of
$\Delta\rm Y$ at the 3 and 7 TeV collision energies.

##  4 Sum rules for triple PDFs

In the previous sections we studied how well the asymmetric and symmetric
dPDFs constructed with P  YTHIA code obey the GS sum rules. Now we would like
to extend the analysis performed in Sections 3.1–3.2 to the case of triple
PDFs. However, since the sum rules for tPDFs do not exist in the literature,
we first need to establish what these sum rules actually are; in Section 4.1
and Section 4.2 we derive the GS sum rules for the bare and renormalized tPDFs
correspondingly. Then, in Section 4.3, we investigate how well the P  YTHIA
tPDFs satisfy these sum rules.

###  4.1 Sum rule for bare tPDFs

We start by writing down the GS sum rules generalized to the case of bare
tPDFs. The momentum sum rule is given by

$\displaystyle\sum\limits_{j_{3}}\int\limits^{1-x_{1}-x_{2}}_{0}dx_{3}~{}x_{3}~{}T^{B}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})=(1-x_{1}-x_{2})D^{B}_{j_{1}j_{2}}(x_{1},x_{2}),$
(62)

where $T^{B}_{j_{1}j_{2}j_{3}}$ and $D^{B}_{j_{1}j_{2}}(x_{1},x_{2})$ are the
bare triple and double parton distribution functions respectively. To prove
this relation we apply the same method as was used to show that GS sum rules
hold for bare dPDFs in Ref. Gaunt:2012tfk . First of all let us write down the
light-cone representation Harindranath:1998pd ; Zhang:1993is ; Zhang:1993dd ;
Harindranath:1993de ; Harindranath:1998pc  of dPDFs and tPDFs:

$\displaystyle D^{B}_{j_{1}j_{2}}(x_{1},x_{2})$
$\displaystyle=\sum\limits_{N,\\{\beta_{i}\\}}\int[dz]_{N}~{}[d^{2}{\bm{k}}]_{N}~{}\left|\Phi_{N}\left(\\{\beta_{i},z_{i},{\bm{k}}_{i}\\}\right)\right|^{2}~{}\sum\limits^{N}_{k}~{}\delta(x_{1}-z_{k})~{}\delta_{j_{1}p_{k}}$
$\displaystyle\mathrel{\phantom{=}}\times\sum\limits^{N}_{l\neq
k}~{}\delta(x_{2}-z_{l})~{}\delta_{j_{2}p_{l}},$ (63) $\displaystyle
T^{B}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})$
$\displaystyle=\sum\limits_{N,\\{\beta_{i}\\}}\int[dz]_{N}~{}[d^{2}{\bm{k}}]_{N}~{}\left|\Phi_{N}\left(\\{\beta_{i},z_{i},{\bm{k}}_{i}\\}\right)\right|^{2}~{}\sum\limits^{N}_{k}~{}\delta(x_{1}-z_{k})~{}\delta_{j_{1}p_{k}}$
$\displaystyle\mathrel{\phantom{=}}\times\sum\limits^{N}_{l\neq
k}~{}\delta(x_{2}-z_{l})~{}\delta_{j_{2}p_{l}}~{}\sum\limits^{N}_{m\neq
k,l}~{}\delta(x_{3}-z_{m})~{}\delta_{j_{3}p_{m}},$ (64)

where

$\displaystyle[dz]_{N}$ $\displaystyle\equiv$
$\displaystyle\prod\limits^{N}_{i=1}~{}dz_{i}~{}\delta\left(1-\sum\limits^{N}_{i}~{}z_{i}\right),$
(65) $\displaystyle\left[d^{2}{\bm{k}}\right]_{N}$ $\displaystyle\equiv$
$\displaystyle\prod\limits^{N}_{i=1}~{}d^{2}{\bm{k}}_{i}~{}\delta^{2}\left(\sum\limits^{N}_{i}~{}{\bm{k}}_{i}\right),$
(66)

and $\Phi_{N}\left(\\{\beta_{i},z_{i},{\bm{k}}_{i}\\}\right)$ is the bare
light-cone wave function to find the hadron being composed from $N$ partons
with longitudinal momentum fractions $\\{z_{i}\\}$, transverse momenta
$\\{{\bm{k}}_{i}\\}$ and other quantum numbers $\\{\beta_{i}\\}$ including
flavour indices $\\{p_{i}\\}$ ($i=1...N$). The sum $\sum_{N,\\{\beta_{i}\\}}$
in Eq. (63) and Eq. (64) runs over all distinguishable Fock states which
reproduce appropriate hadronic quantum numbers.

By substituting Eq. (64) in Eq. (62) we get

$\displaystyle\sum\limits_{j_{3}}\int\limits^{1-x_{1}-x_{2}}_{0}dx_{3}~{}x_{3}~{}T^{B}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})=\sum\limits_{N,\\{\beta_{i}\\}}\int[dz]_{N}~{}[d^{2}{\bm{k}}]_{N}~{}\left|\Phi_{N}\left(\\{\beta_{i},z_{i},{\bm{k}}_{i}\\}\right)\right|^{2}$
$\displaystyle\mathrel{\phantom{=}}\times\sum\limits_{j_{3}}\int\limits^{1-x_{1}-x_{2}}_{0}dx_{3}~{}x_{3}~{}\sum\limits^{N}_{k}~{}\delta(x_{1}-z_{k})~{}\delta_{j_{1}p_{k}}~{}\sum\limits^{N}_{l\neq
k}~{}\delta(x_{2}-z_{l})~{}\delta_{j_{2}p_{l}}~{}\sum\limits^{N}_{m\neq
k,l}~{}\delta(x_{3}-z_{m})~{}\delta_{j_{3}p_{m}}$
$\displaystyle=\sum\limits_{N,\\{\beta_{i}\\}}\int[dz]_{N}~{}[d^{2}{\bm{k}}]_{N}~{}\left|\Phi_{N}\left(\\{\beta_{i},z_{i},{\bm{k}}_{i}\\}\right)\right|^{2}~{}\sum\limits^{N}_{k}~{}\delta(x_{1}-z_{k})~{}\delta_{j_{1}p_{k}}~{}\sum\limits^{N}_{l\neq
k}~{}\delta(x_{2}-z_{l})~{}\delta_{j_{2}p_{l}}$
$\displaystyle\mathrel{\phantom{=}}\times\sum\limits_{j_{3}}\sum\limits^{N}_{m\neq
k,l}z_{m}~{}\delta_{j_{3}p_{m}}.$ (67)

The last sum in Eq. (67) can be written as

$\displaystyle\sum\limits_{j_{3}}\sum\limits^{N}_{m\neq
k,l}z_{m}~{}\delta_{j_{3}p_{m}}=\sum\limits^{N}_{m}z_{m}-z_{k}-z_{l}=1-x_{1}-x_{2},$
(68)

where we used the fact that $\sum^{N}_{m}z_{m}=1$ and that $z_{k}=x_{1}$ and
$z_{l}=x_{2}$. Therefore, Eq. (67) turns into

$\displaystyle\sum\limits_{j_{3}}\int\limits^{1-x_{1}-x_{2}}_{0}$
$\displaystyle dx_{3}~{}x_{3}~{}T^{B}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})$
$\displaystyle=$
$\displaystyle\left(1-x_{1}-x_{2}\right)\sum\limits_{N,\\{\beta_{i}\\}}\int[dz]_{N}~{}[d^{2}{\bm{k}}]_{N}~{}\left|\Phi_{N}\left(\\{\beta_{i},z_{i},{\bm{k}}_{i}\\}\right)\right|^{2}$
$\displaystyle\times\sum\limits^{N}_{k}~{}\delta(x_{1}-z_{k})~{}\delta_{j_{1}p_{k}}~{}\sum\limits^{N}_{l\neq
k}~{}\delta(x_{2}-z_{l})~{}\delta_{j_{2}p_{l}}$ $\displaystyle=$
$\displaystyle\left(1-x_{1}-x_{2}\right)D^{B}_{j_{1}j_{2}}(x_{1},x_{2}).$ (69)

Now let us concentrate on the generalization of the number rule Eq. (6) for
the case of bare tPDFs. First of all let us note that

$\displaystyle\int\limits^{1-x_{1}-x_{2}}_{0}dx_{3}~{}T^{B}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})$
$\displaystyle=\sum\limits_{N,\\{\beta_{i}\\}}\int[dz]_{N}~{}[d^{2}{\bm{k}}]_{N}~{}\left|\Phi_{N}\left(\\{\beta_{i},z_{i},{\bm{k}}_{i}\\}\right)\right|^{2}$
(70)
$\displaystyle\mathrel{\phantom{=}}\times\sum\limits^{N}_{k}~{}\delta(x_{1}-z_{k})~{}\delta_{j_{1}p_{k}}~{}\sum\limits^{N}_{l\neq
k}~{}\delta(x_{2}-z_{l})~{}\delta_{j_{2}p_{l}}~{}\sum\limits^{N}_{m\neq
k,l}~{}\delta_{j_{3}p_{m}},$

and, therefore,

$\displaystyle\int\limits^{1-x_{1}-x_{2}}_{0}dx_{3}~{}T^{B}_{j_{1}j_{2}j_{3v}}(x_{1},x_{2},x_{3})\equiv\int\limits^{1-x_{1}-x_{2}}_{0}dx_{3}~{}\left[T^{B}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})-T^{B}_{j_{1}j_{2}\bar{j}_{3}}(x_{1},x_{2},x_{3})\right]$
$\displaystyle=\sum\limits_{N,\\{\beta_{i}\\}}\int[dz]_{N}~{}[d^{2}{\bm{k}}]_{N}~{}\left|\Phi_{N}\left(\\{\beta_{i},z_{i},\bm{k}_{i}\\}\right)\right|^{2}\sum\limits^{N}_{k}~{}\delta(x_{1}-z_{k})~{}\delta_{j_{1}p_{k}}~{}\sum\limits^{N}_{l\neq
k}~{}\delta(x_{2}-z_{l})~{}\delta_{j_{2}p_{l}}$
$\displaystyle\mathrel{\phantom{=}}\times\left(\sum\limits^{N}_{m\neq
k,l}~{}\delta_{j_{3}p_{m}}-\sum\limits^{N}_{m\neq
k,l}~{}\delta_{\bar{j}_{3}p_{m}}\right).$ (71)

The two sums in the parentheses in Eq. (71) can be written as

$\displaystyle\sum\limits^{N}_{m\neq k,l}~{}\delta_{j_{3}p_{m}}$
$\displaystyle=$
$\displaystyle\sum\limits^{N}_{m}~{}\delta_{j_{3}p_{m}}-\delta_{j_{3}j_{1}}-\delta_{j_{3}j_{2}}=N_{j_{3}|\\{\beta_{i}\\}}-\delta_{j_{3}j_{1}}-\delta_{j_{3}j_{2}},$
$\displaystyle\sum\limits^{N}_{m\neq k,l}~{}\delta_{\bar{j}_{3}p_{m}}$
$\displaystyle=$
$\displaystyle\sum\limits^{N}_{m}~{}\delta_{\bar{j}_{3}p_{m}}-\delta_{\bar{j}_{3}j_{1}}-\delta_{\bar{j}_{3}j_{2}}=N_{\bar{j}_{3}|\\{\beta_{i}\\}}-\delta_{\bar{j}_{3}j_{1}}-\delta_{\bar{j}_{3}j_{2}},$
(72)

where $N_{j_{3}|\\{\beta_{i}\\}}$ ($N_{\bar{j}_{3}|\\{\beta_{i}\\}}$) is the
number of $j_{3}$ ($\bar{j}_{3}$) partons in state $\\{\beta_{i}\\}$. The
difference between $N_{j_{3}|\\{\beta_{i}\\}}$ and
$N_{\bar{j}_{3}|\\{\beta_{i}\\}}$ is the number of valence quarks, which is
independent of $\\{\beta_{i}\\}$:

$\displaystyle
N_{j_{3v}}=N_{j_{3}|\\{\beta_{i}\\}}-N_{\bar{j}_{3}|\\{\beta_{i}\\}}.$ (73)

Inserting Eqs. (72) and (73) into Eq. (71), we obtain

$\displaystyle\int\limits^{1-x_{1}-x_{2}}_{0}dx_{3}~{}T^{B}_{j_{1}j_{2}j_{3v}}(x_{1},x_{2},x_{3})$
$\displaystyle=\left(N_{j_{3v}}-\delta_{j_{3}j_{1}}-\delta_{j_{3}j_{2}}+\delta_{\bar{j}_{3}j_{1}}+\delta_{\bar{j}_{3}j_{2}}\right)$
$\displaystyle\mathrel{\phantom{=}}\times D^{B}_{j_{1}j_{2}}(x_{1},x_{2}),$
(74)

which is Eq. (6) generalized for the case of bare tPDFs.

###  4.2 Sum rules for renormalized tPDFs

We have shown that Eq. (62), (74), hold for the bare distributions. Now let us
sketch the all order proof for the renormalized tPDFs following the method
similar to the one used in Ref. Diehl:2018kgr . Using notation of Refs.
Diehl:2011yj ; Diehl:2018kgr , the bare sPDF, dPDF and tPDF are defined as

$\displaystyle f^{B}_{j_{1}}(x_{1})$
$\displaystyle=(x_{1}p^{+})^{-n_{1}}\int\,\frac{dz^{-}_{1}}{2\pi}\,e^{ix_{1}z^{-}_{1}p^{+}}\,\langle
p|\mathcal{O}_{j_{1}}(0,z_{1})|p\rangle\biggr{\rvert}_{z^{+}_{1}=0,\,\bm{z}_{1}=\bm{0}},$
(75) $\displaystyle D^{B}_{j_{1}j_{2}}(x_{1},x_{2})$
$\displaystyle=(x_{1}p^{+})^{-n_{1}}(x_{2}p^{+})^{-n_{2}}2p^{+}\int\,\frac{dz^{-}_{1}}{2\pi}\,\frac{dz^{-}_{2}}{2\pi}\,dy^{-}\,d^{2-2\varepsilon}\bm{y}\,e^{i(\sum_{i=1}^{2}x_{i}z^{-}_{i})p^{+}}$
$\displaystyle\mathrel{\phantom{=}}\times\langle
p|\mathcal{O}_{j_{1}}(y,z_{1})\mathcal{O}_{j_{2}}(0,z_{2})|p\rangle\biggr{\rvert}_{z^{+}_{i}=y^{+}=0,\,\bm{z}_{i}=\bm{0}},$
(76) $\displaystyle T^{B}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})$
$\displaystyle=(x_{1}p^{+})^{-n_{1}}(x_{2}p^{+})^{-n_{2}}(x_{3}p^{+})^{-n_{3}}(2p^{+})^{2}$
$\displaystyle\mathrel{\phantom{=}}\times\int\,\frac{dz^{-}_{1}}{2\pi}\,\frac{dz^{-}_{2}}{2\pi}\,\frac{dz^{-}_{3}}{2\pi}\,dy_{1}^{-}dy_{2}^{-}\,d^{2-2\varepsilon}\bm{y}_{1}\,d^{2-2\varepsilon}\bm{y}_{2}\,e^{i(\sum_{i=1}^{3}x_{i}z^{-}_{i})p^{+}}$
$\displaystyle\mathrel{\phantom{=}}\times\langle
p|\mathcal{O}_{j_{1}}(y_{1},z_{1})\mathcal{O}_{j_{2}}(y_{2},z_{2})\mathcal{O}_{j_{3}}(0,z_{3})|p\rangle\biggr{\rvert}_{z^{+}_{i}=y_{i}^{+}=0,\,\bm{z}_{i}=\bm{0}},$
(77)

where we use the light-cone representation for four-vectors:
$v^{\pm}=(v^{0}\pm v^{3})/\sqrt{2}$, $\bm{v}=(v^{1},v^{2})$. The quantity
$n_{j}$ is $0$ if $j$ is a quark or antiquark, and $1$ if $j$ is a gluon. The
twist-two operators $\mathcal{O}$ in Eqs. (75)-(77) are bilinear in partonic
fields and are defined as

$\displaystyle\mathcal{O}_{q}(y,z)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,\bar{q}\left(y-\frac{z}{2}\right)\,\gamma^{+}\,q\left(y+\frac{z}{2}\right),$
(78) $\displaystyle\mathcal{O}_{\bar{q}}(y,z)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\,\bar{q}\left(y+\frac{z}{2}\right)\,\gamma^{+}\,q\left(y-\frac{z}{2}\right),$
(79) $\displaystyle\mathcal{O}_{g}(y,z)$ $\displaystyle=$ $\displaystyle
G^{+i}\left(y-\frac{z}{2}\right)\,G^{+i}\left(y+\frac{z}{2}\right),$ (80)

where $q,$ and $\bar{q}$ are the quark and antiquark fields, $G$ is the gluon
field strength, and $i$ is a transverse component151515Note that these
definitions are strictly speaking only appropriate when we are using light-
cone gauge $A^{+}=0$. In other gauges Wilson line factors appear in the
definitions of the operators $\mathcal{O}_{j}$..

As was noted in Ref. Diehl:2011yj , there is a UV singularity in Eq. (76) when
$\bm{y}=0$ and in Eq. (77) when $\bm{y}_{i}=0$ or $\bm{y}_{1}=\bm{y}_{2}$. In
Section 2 a cut-off on $\bm{y}$ was used to regulate the integral for the
dPDF, but as was also mentioned there, the dPDF sum rules only hold exactly
when we use dimensional regularization and $\overline{\rm MS}$ renormalization
– thus, here we will follow the latter approach. For the bare dPDF and tPDF in
Eqs. (76) and (77) we should change the integrations with respect to
$d^{2}\bm{y}_{(i)}$ to $d^{2-2\varepsilon}\bm{y}_{(i)}$ (in fact in the
previous section, integrations over transverse momenta $\bm{k}_{i}$ should be
similarly altered to regularize the UV divergences, although this does not
affect the overall argument presented there).

Generalizing the renormalization equation for the dPDFs, we write down the
following equation for the renormalization of the tPDFs:

$\displaystyle T_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3},Q)$
$\displaystyle=\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}Z_{j_{1}j^{\prime}_{1}}(Q)\otimes_{1}Z_{j_{2}j^{\prime}_{2}}(Q)\otimes_{2}Z_{j_{3}j^{\prime}_{3}}(Q)\otimes_{3}T^{B}_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}(z_{1},z_{2},z_{3})$
$\displaystyle\mathrel{\phantom{=}}+\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j^{\prime}_{1}}(Q)\otimes_{1}Z_{j_{2}j_{3},j^{\prime}_{2}}(Q)\otimes_{23}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}(z_{1},z_{2})+\mathrm{(permutations)}$
$\displaystyle\mathrel{\phantom{=}}+\sum\limits_{j^{\prime}_{1}}Z_{j_{1}j_{2}j_{3},j^{\prime}_{1}}(Q)\otimes_{123}f^{B}_{j^{\prime}_{1}}(z_{1}),$
(81)

where we use the same notation as in Ref. Diehl:2018kgr . The definitions of
the convolution symbols in this expression are given by:

$\displaystyle A\otimes_{1}B$ $\displaystyle=$
$\displaystyle\int\frac{dz}{z}A\left(\frac{x_{1}}{z}\right)B(z),$ (82)
$\displaystyle A\otimes_{12}B$ $\displaystyle=$
$\displaystyle\int\frac{dz}{z^{2}}A\left(\frac{x_{1}}{z},\frac{x_{2}}{z}\right)B(z),$
(83) $\displaystyle A\otimes_{123}B$ $\displaystyle=$
$\displaystyle\int\frac{dz}{z^{3}}A\left(\frac{x_{1}}{z},\frac{x_{2}}{z},\frac{x_{3}}{z}\right)B(z).$
(84)

The terms on the first two lines of Eq. (81) are analogous to the two terms in
the renormalization of the dPDFs. The quantity $Z_{ij}$ is the renormalization
factor that appears in the renormalization of the sPDFs, and $Z_{ij,k}$ is a
renormalization factor linked with the $1\to 2$ splitting singularity. The
term “(permutations)” in Eq. 81 refers to two additional terms which have the
same structure as the first term on the second line of that equation, but
which have $\\{j_{1},j_{2},j_{3}\\}$ replaced by $\\{j_{2},j_{1},j_{3}\\}$ or
$\\{j_{3},j_{1},j_{2}\\}$. The final term in Eq. (81) is new and is associated
with a $1\to 3$ splitting process.

Let us use the particular implementation of the $\overline{\rm MS}$ scheme
discussed in Eqs. (79) and (80) of Ref. Diehl:2018kgr . Then, the structures
of the renormalization factors are as follows:

$\displaystyle Z_{ij}$
$\displaystyle=\delta_{ij}\delta(1-x)+\sum_{n=1}^{\infty}\alpha_{s}^{n}\sum_{m\geq
1}\dfrac{Z_{nm;ij}}{\varepsilon^{m}},$ (85) $\displaystyle Z_{jk,i}$
$\displaystyle=\sum_{n=1}^{\infty}\alpha_{s}^{n}\sum_{m\geq
1}\dfrac{Z_{nm;jk,i}}{\varepsilon^{m}},$ (86) $\displaystyle Z_{jkl,i}$
$\displaystyle=\sum_{n=2}^{\infty}\alpha_{s}^{n}\sum_{m\geq
1}\dfrac{Z_{nm;jk,i}}{\varepsilon^{m}}.$ (87)

Note that $Z_{jk,i}$ begins at order $\alpha_{s}$, whilst $Z_{jk,i}$ begins at
order $\alpha_{s}^{2}$ (since they are related to the $1\to 2$ and $1\to 3$
processes respectively). This then means that, in the $\overline{\rm MS}$
scheme, they contain no finite pieces, only pole terms.

#### Momentum sum rule.

To prove the momentum sum rule for the renormalized tPDFs, let us consider:

$\displaystyle\Delta^{j_{1}j_{2}}_{\text{mtm}}\equiv\sum_{j_{3}}\int_{3}X_{3}T_{j_{1}j_{2}j_{3}}-(1-X_{1}-X_{2})D_{j_{1}j_{2}},$
(88)

where here (and in what follows) we use the same notation as is introduced in
Section 6.1 of Diehl:2018kgr . For the momentum sum rule to hold,
$\Delta^{j_{1}j_{2}}_{\text{mtm}}$ must be zero. The first term on the right
hand side of Eq. (88), after inserting the expression for
$T_{j_{1}j_{2}j_{3}}$ from Eq. (81), becomes:

$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}\sum_{j_{3}}\int_{3}X_{3}Z_{j_{3}j^{\prime}_{3}}\otimes_{3}T^{B}_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}$
$\displaystyle+$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j_{2},j^{\prime}_{1}}\otimes_{12}\sum_{j_{3}}\int_{3}X_{3}Z_{j_{3}j^{\prime}_{2}}\otimes_{3}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}+\mathrm{(permutations)}$
$\displaystyle+$
$\displaystyle\sum\limits_{j^{\prime}_{1}}\sum_{j_{3}}\int_{3}X_{3}Z_{j_{1}j_{2}j_{3},j^{\prime}_{1}}\otimes_{123}f^{B}_{j^{\prime}_{1}}.$
(89)

We now simplify each of the terms in this expression. For the first term in
Eq. (4.2), we have:

$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}\sum_{j_{3}}\int_{3}X_{3}Z_{j_{3}j^{\prime}_{3}}\otimes_{3}T^{B}_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}\left(\sum_{j_{3}}\int
XZ_{j_{3}j^{\prime}_{3}}\right)\int_{3}X_{3}T^{B}_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}(\hbox{\pagecolor{red!50}{$\displaystyle
1$}}-\hbox{\pagecolor{yellow!50}{$\displaystyle
X_{1}$}}-\hbox{\pagecolor{green!50}{$\displaystyle
X_{2}$}})D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}.$ (90)

Going from the first line to the second line, we have used Eq. (64) from Ref.
Diehl:2018kgr , and going from the second to the third line, we have used the
momentum sum rule for the bare tPDFs, plus the momentum sum rule for the
$Z_{ij}$: $\sum_{j_{3}}\int XZ_{j_{3}j^{\prime}_{3}}=1$ (see Eq. (97) in Ref.
Diehl:2018kgr ).

For the second term in Eq. (4.2), we have:

$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j_{2},j^{\prime}_{1}}\otimes_{12}\sum_{j_{3}}\int_{3}X_{3}Z_{j_{3}j^{\prime}_{2}}\otimes_{3}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j_{2},j^{\prime}_{1}}\otimes_{12}\left(\sum_{j_{3}}\int
XZ_{j_{3}j^{\prime}_{2}}\right)\int_{2}X_{2}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j_{1}^{\prime}}Z_{j_{1}j_{2},j_{1}^{\prime}}\otimes_{12}(1-X)f^{B}_{j_{1}^{\prime}}.$
(91)

Once again we use Eq. (64) from Ref. Diehl:2018kgr  to go from the first to
the second line, and then we use the momentum sum rule for the $Z_{ij}$ along
with the momentum sum rule for the bare dPDFs (as proven in Ref. Diehl:2018kgr
).

In one of the terms denoted by “(permutations)”, we have:

$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1},j^{\prime}_{1}}\otimes_{1}\sum_{j_{3}}\int_{3}X_{3}Z_{j_{2}j_{3},j^{\prime}_{2}}\otimes_{23}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1},j^{\prime}_{1}}\otimes_{1}\sum_{j_{3}}\left(\int_{3}X_{3}Z_{j_{2}j_{3},j^{\prime}_{2}}\right)\otimes_{2}\left(X_{2}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}\right)$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1},j^{\prime}_{1}}\otimes_{1}\left[(\hbox{\pagecolor{green!50}{$\displaystyle
1$}}-\hbox{\pagecolor{blue!50}{$\displaystyle
X_{2}$}})Z_{j_{2}j_{2}^{\prime}}\right]\otimes_{2}\left(X_{2}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}\right).$
(92)

Here we use Eq. (66) from Ref. Diehl:2018kgr  to go from the first to the
second line, followed by the momentum sum rule for the $Z_{ij,k}$ factors as
given in Eq. (103) of Ref. Diehl:2018kgr . The other term in the
“(permutations)” has the same structure but with $1\leftrightarrow 2$.

Let us now consider the second term in Eq. (88). Expressing $D_{j_{1}j_{2}}$
in terms of bare quantities, using Eq. (70) from Ref. Diehl:2018kgr , we have:

$\displaystyle-$ $\displaystyle(1-X_{1}-X_{2})D_{j_{1}j_{2}}$ $\displaystyle=$
$\displaystyle(\hbox{\pagecolor{yellow!50}{$\displaystyle
X_{1}$}}+\hbox{\pagecolor{blue!50}{$\displaystyle
X_{2}$}}-\hbox{\pagecolor{red!50}1})\sum_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}+(X_{1}+X_{2}-1)\sum\limits_{j^{\prime}}Z_{j_{1}j_{2},j^{\prime}}\otimes_{12}f^{B}_{j^{\prime}}.$
(93)

Now many of the terms on the right hand side of the expression for
$\Delta_{\text{mtm}}$ cancel. The terms in the  ,  ,  boxes cancel against
each other directly, whilst the terms in the  boxes cancel against the other
“(permutations)” term that we have not written out explicitly. All the terms
containing $D^{B}$ on the right hand side then cancel out, and we are left
with:

$\displaystyle\Delta^{j_{1}j_{2}}_{\text{mtm}}=$
$\displaystyle\sum\limits_{j^{\prime}}\sum_{j_{3}}\int_{3}X_{3}Z_{j_{1}j_{2}j_{3},j^{\prime}}\otimes_{123}f^{B}_{j^{\prime}}+\sum\limits_{j^{\prime}}Z_{j_{1}j_{2},j^{\prime}}\otimes_{12}(1-X)f^{B}_{j^{\prime}}$
$\displaystyle+(X_{1}+X_{2}-1)\sum\limits_{j^{\prime}}Z_{j_{1}j_{2},j^{\prime}}\otimes_{12}f^{B}_{j^{\prime}}$
$\displaystyle=$
$\displaystyle\sum_{j^{\prime}k^{\prime}}\left(\sum_{j_{3}}\int_{3}X_{3}Z_{j_{1}j_{2}j_{3},j^{\prime}}-(1-X_{1}-X_{2})Z_{j_{1}j_{2},j^{\prime}}\right)\otimes_{12}\left(XZ^{-1}_{j^{\prime}k^{\prime}}\right)\otimes\left(Xf_{k^{\prime}}\right).$
(94)

To reach the last line we have used:

$\displaystyle\int_{3}X_{3}^{n}T\otimes_{123}A=\left(\int_{3}X_{3}^{n}T\right)\otimes_{12}X^{n}A$
(95)

and

$\displaystyle(X_{1}+X_{2})\left(D\otimes_{12}A\right)=\left([X_{1}+X_{2}]D\right)\otimes_{12}\left(XA\right).$
(96)

The quantity $\Delta^{j_{1}j_{2}}_{\text{mtm}}$ must be finite, as it is
expressed in terms of renormalized quantities. However, in Eq. (4.2) we see
that it can be expressed in terms of a pure pole factor convolved with a
finite quantity ($f_{k^{\prime}}$). The only way this can be true is if both
left and right hand side are zero - thus $\Delta^{j_{1}j_{2}}_{\text{mtm}}=0$
and the momentum sum rule for tPDFs holds at the renormalized level.

Given that $\Delta^{j_{1}j_{2}}_{\text{mtm}}=0$, we obtain from Eq. (4.2) a
sum rule for the $Z_{ijk,l}$ factors:

$\displaystyle\sum_{j_{3}}\int_{3}X_{3}Z_{j_{1}j_{2}j_{3},j^{\prime}}=(1-X_{1}-X_{2})Z_{j_{1}j_{2},j^{\prime}}.$
(97)

#### Number sum rule.

The proof of the number sum rule for renormalized tPDFs proceeds in an
analogous fashion. Consider:

$\displaystyle\Delta^{j_{1}j_{2}}_{\text{num}}\equiv\int_{3}T_{j_{1}j_{2}j_{3v}}-(N_{j_{3}v}+\delta_{j_{3}\bar{j}_{1}}+\delta_{j_{3}\bar{j}_{2}}-\delta_{j_{3}j_{1}}+\delta_{j_{3}j_{2}})D_{j_{1}j_{2}}.$
(98)

If this is zero, the number sum rule holds. We consider
$\int_{3}T_{j_{1}j_{2}j_{3v}}$, and express it in terms of bare quantities
using Eq. (81):

$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}\int_{3}Z_{j_{3v}j^{\prime}_{3}}\otimes_{3}T^{B}_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}$
$\displaystyle+$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j_{2},j^{\prime}_{1}}\otimes_{12}\int_{3}Z_{j_{3v}j^{\prime}_{2}}\otimes_{3}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}+\mathrm{(permutations)}$
$\displaystyle+$
$\displaystyle\sum\limits_{j^{\prime}_{1}}\int_{3}Z_{j_{1}j_{2}j_{3v},j^{\prime}_{1}}\otimes_{123}f^{B}_{j^{\prime}_{1}}.$
(99)

Using Eq. (64) from Ref. Diehl:2018kgr , the number sum rule for the $Z_{ij}$
(see Eq. (90) from Ref. Diehl:2018kgr ), and the number sum rule for the bare
tPDFs, we may express the first term of Eq. (4.2) as:

$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}\int_{3}Z_{j_{3v}j^{\prime}_{3}}\otimes_{3}T^{B}_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}\left(\delta_{j_{3}j^{\prime}_{3}}-\delta_{j_{3}\bar{j}^{\prime}_{3}}\right)\int_{3}T^{B}_{j^{\prime}_{1}j^{\prime}_{2}j^{\prime}_{3}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}(\hbox{\pagecolor{red!50}{$\displaystyle
N_{j_{3}v}$}}+\hbox{\pagecolor{yellow!50}{$\displaystyle\delta_{j_{3}\bar{j}^{\prime}_{1}}$}}+\hbox{\pagecolor{green!50}{$\displaystyle\delta_{j_{3}\bar{j}^{\prime}_{2}}$}}-\hbox{\pagecolor{yellow!50}{$\displaystyle\delta_{j_{3}j^{\prime}_{1}}$}}-\hbox{\pagecolor{orange!50}{$\displaystyle\delta_{j_{3}j_{2}^{\prime}}$}})D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}.$
(100)

For the second term, we use Eq. (64) from Ref. Diehl:2018kgr , the number sum
rule for the $Z_{ij}$, and the number sum rule for the bare dPDFs, to obtain:

$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j_{2},j^{\prime}_{1}}\otimes_{12}\int_{3}Z_{j_{3v}j^{\prime}_{2}}\otimes_{3}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j_{2},j^{\prime}_{1}}\otimes_{12}\left(\delta_{j_{3}j_{2}^{\prime}}-\delta_{j_{3}\bar{j}_{2}^{\prime}}\right)\int_{2}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j_{1}^{\prime}}Z_{j_{1}j_{2},j_{1}^{\prime}}\otimes_{12}(N_{j_{3}v}+\delta_{j_{3}\bar{j}^{\prime}_{1}}-\delta_{j_{3}j^{\prime}_{1}})f^{B}_{j_{1}^{\prime}}.$
(101)

In one of the terms in the “(permutations)”, we can use Eq. (66) from Ref.
Diehl:2018kgr  and the number sum rule for the $Z_{ij,k}$ (as given in Eq.
(94) of Ref. Diehl:2018kgr ):

$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1},j^{\prime}_{1}}\otimes_{1}\int_{3}Z_{j_{2}j_{3v},j^{\prime}_{2}}\otimes_{23}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1},j^{\prime}_{1}}\otimes_{1}\left(\int_{3}Z_{j_{2}j_{3v},j^{\prime}_{2}}\right)\otimes_{2}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}$
$\displaystyle=$
$\displaystyle\sum\limits_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1},j^{\prime}_{1}}\otimes_{1}(\hbox{\pagecolor{blue!50}{$\displaystyle\delta_{j_{3}\bar{j}_{2}}$}}-\hbox{\pagecolor{purple!50}{$\displaystyle\delta_{j_{3}j_{2}}$}}-\hbox{\pagecolor{green!50}{$\displaystyle\delta_{j_{3}\bar{j}^{\prime}_{2}}$}}+\hbox{\pagecolor{orange!50}{$\displaystyle\delta_{j_{3}j_{2}^{\prime}}$}})Z_{j_{2}j_{2}^{\prime}}\otimes_{2}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}.$
(102)

The other term in the “(permutations)” has the same structure, but with
$1\leftrightarrow 2$.

We express the second term in Eq. (98) in terms of bare quantities using Eq.
(70) from Ref. Diehl:2018kgr :

$\displaystyle-$
$\displaystyle(N_{j_{3}v}+\delta_{j_{3}\bar{j}_{1}}+\delta_{j_{3}\bar{j}_{2}}-\delta_{j_{3}j_{1}}-\delta_{j_{3}j_{2}})D_{j_{1}j_{2}}$
$\displaystyle=$ $\displaystyle(-\hbox{\pagecolor{red!50}{$\displaystyle
N_{j_{3}v}$}}-\hbox{\pagecolor{yellow!50}{$\displaystyle\delta_{j_{3}\bar{j}_{1}}$}}-\hbox{\pagecolor{blue!50}{$\displaystyle\delta_{j_{3}\bar{j}_{2}}$}}+\hbox{\pagecolor{yellow!50}{$\displaystyle\delta_{j_{3}j_{1}}$}}+\hbox{\pagecolor{purple!50}{$\displaystyle\delta_{j_{3}j_{2}}$}})\sum_{j^{\prime}_{1}j^{\prime}_{2}}Z_{j_{1}j^{\prime}_{1}}\otimes_{1}Z_{j_{2}j^{\prime}_{2}}\otimes_{2}D^{B}_{j^{\prime}_{1}j^{\prime}_{2}}$
$\displaystyle+$
$\displaystyle(-N_{j_{3}v}-\delta_{j_{3}\bar{j}_{1}}-\delta_{j_{3}\bar{j}_{2}}+\delta_{j_{3}j_{1}}+\delta_{j_{3}j_{2}})\sum\limits_{j^{\prime}}Z_{j_{1}j_{2},j^{\prime}}\otimes_{12}f^{B}_{j^{\prime}}.$
(103)

As for the momentum sum rule, a number of cancellations now occur. The terms
in the  ,  ,  ,  ,  boxes cancel against each other directly, whilst the terms
in the  boxes cancel against the other “(permutations)” term that we have not
written out explicitly. Once again, the terms with $D^{B}$ cancel out, and we
obtain (using Eq. (95)):

$\displaystyle\Delta^{j_{1}j_{2}}_{\text{num}}$
$\displaystyle=(-N_{j_{3}v}-\delta_{j_{3}\bar{j}_{1}}-\delta_{j_{3}\bar{j}_{2}}+\delta_{j_{3}j_{1}}+\delta_{j_{3}j_{2}})\sum\limits_{j^{\prime}}Z_{j_{1}j_{2},j^{\prime}}\otimes_{12}f^{B}_{j^{\prime}}$
$\displaystyle\mathrel{\phantom{=}}+\sum\limits_{j_{1}^{\prime}}Z_{j_{1}j_{2},j_{1}^{\prime}}\otimes_{12}(N_{j_{3}v}+\delta_{j_{3}\bar{j}^{\prime}_{1}}-\delta_{j_{3}j^{\prime}_{1}})f^{B}_{j_{1}^{\prime}}$
$\displaystyle\mathrel{\phantom{=}}+\sum\limits_{j^{\prime}_{1}}\int_{3}Z_{j_{1}j_{2}j_{3v},j^{\prime}_{1}}\otimes_{123}f^{B}_{j^{\prime}_{1}}$
$\displaystyle=\sum_{j^{\prime}_{1},k^{\prime}}\left(\int_{3}Z_{j_{1}j_{2}j_{3v},j^{\prime}_{1}}+\left[-\delta_{j_{3}\bar{j}_{1}}-\delta_{j_{3}\bar{j}_{2}}+\delta_{j_{3}j_{1}}+\delta_{j_{3}j_{2}}+\delta_{j_{3}\bar{j}^{\prime}_{1}}-\delta_{j_{3}j^{\prime}_{1}}\right]Z_{j_{1}j_{2},j_{1}^{\prime}}\right)$
$\displaystyle\qquad\quad\otimes_{12}Z^{-1}_{j_{1}^{\prime}k^{\prime}}\otimes
f_{k^{\prime}}.$ (104)

According to a similar argument as was used for $\Delta_{\text{mtm}}$, we find
that $\Delta_{\text{num}}=0$ – that is, the number sum rule holds for the
renormalized tPDFs. We also obtain a number sum rule for the $Z_{ijk,l}$
factors:

$\displaystyle\int_{3}Z_{j_{1}j_{2}j_{3v},j^{\prime}_{1}}=\left(\delta_{j_{3}\bar{j}_{1}}+\delta_{j_{3}\bar{j}_{2}}-\delta_{j_{3}j_{1}}-\delta_{j_{3}j_{2}}-\delta_{j_{3}\bar{j}^{\prime}_{1}}+\delta_{j_{3}j^{\prime}_{1}}\right)Z_{j_{1}j_{2},j_{1}^{\prime}}.$
(105)

###  4.3 Check of the GS sum rules for tPDFs constructed with P  YTHIA code

In Sections 4.1 \- 4.2 we proved the GS sum rules generalized for the case of
tPDFs. Now let us demonstrate how one can build tPDFs which approximately obey
Eq. (62) and Eq. (74). As in Section 3.1 we begin the discussion with the case
of the asymmetric tPDFs which are given by

$\displaystyle T_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3},Q)$ $\displaystyle=$
$\displaystyle f^{r}_{j_{1}}(x_{1},Q)~{}f^{m\leftarrow
j_{1},x_{1}}_{j_{2}}(x_{2},Q)~{}f^{m\leftarrow
j_{1},x_{1};j_{2},x_{2}}_{j_{3}}(x_{3},Q).$ (106)

As for the dPDFs, we only check the sum rules when integrating over the final
parton (parton 3 here) – this is the case in which we expect the sum rules to
be satisfied the best. To construct tPDFs we use the P  YTHIA code as it is
explained in Appendix A. In this section we always set all scales in the tPDF
to be equal (this is the case for which the tPDF sum rules should hold), and
set the value of this common scale to $91$ GeV.

First of all let us check the momentum sum rule given by Eq. (62). Inserting
Eq. (106) into Eq. (62), we get:

$\displaystyle\frac{1}{1-x_{1}-x_{2}}\,\sum\limits_{j_{3}}\int\limits^{1-x_{1}-x_{2}}_{0}dx_{3}~{}x_{3}~{}f^{m\leftarrow
j_{1},x_{1};j_{2},x_{2}}_{j_{3}}(x_{3},Q)\overset{?}{=}1.$ (107)

The results of the numerical check of Eq. (107) are given in Table 8. The
problem of oscillations for the P  YTHIA dPDFs mentioned in Section 3 becomes
stronger for the case of tPDFs. Therefore, we increase a number of averaging
calls by one order of magnitude comparing to the case of dPDFs. For our check
we choose $j_{1}=j_{2}=u$, and vary $x_{1}$ between $10^{-6}$ and $0.8$
intervals while keeping $x_{2}$ equal to $10^{-4}$. As a baseline for our
computation we use naive tPDFs defined in analogy with Eq. (39) as

$\displaystyle
T_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3},Q)=f^{r}_{j_{1}}(x_{1},Q)~{}f^{r}_{j_{2}}(x_{2},Q)~{}f^{r}_{j_{3}}(x_{3},Q)~{}\theta(1-x_{1}-x_{2}-x_{3}).$
(108)

$x_{1}$ | $x_{2}$ | $j_{1}$ | $j_{2}$ | P  YTHIA tPDFs | Naive tPDFs
---|---|---|---|---|---
$10^{-6}$ | $10^{-4}$ | $u$ | $u$ | 0.996 | 0.996
$10^{-3}$ | $10^{-4}$ | $u$ | $u$ | 0.996 | 0.996
$10^{-1}$ | $10^{-4}$ | $u$ | $u$ | 1.008 | 1.106
0.2 | $10^{-4}$ | $u$ | $u$ | 1.010 | 1.244
0.4 | $10^{-4}$ | $u$ | $u$ | 1.011 | 1.649
0.8 | $10^{-4}$ | $u$ | $u$ | 1.011 | 4.119
Table 8: Test of the momentum sum rule for the tPDFs. In the ideal situation
when the momentum sum rule is perfectly satisfied the number in each cell
should be equal to unity.

Now let us check the number rule given by Eq. (74). In order to do that,
similarly to Eq. (34), we define

$\displaystyle R_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3},Q)\equiv
x_{3}\,\frac{T_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3},Q)-T_{j_{1}j_{2}\bar{j_{3}}}(x_{1},x_{2},x_{3},Q)}{f_{j_{1}}^{r}(x_{1},Q)f_{j_{2}}^{m\leftarrow
j_{1},x_{1}}(x_{2},Q)},$ (109)

which can be seen as a response of the valence sPDF $f_{j_{3v}}(x_{3},Q)$ to
the two interactions involving partons $j_{1}$ and $j_{2}$ with momentum
fractions $x_{1}$ and $x_{2}$ scales respectively. Due to the large number of
different flavour combinations we restrict ourselves to a few illustrative
cases. Let us first look at the $uuu$ and $u\bar{u}u$ combinations. The
results of our numerical checks are given in Tables 9 \- 10 and the
corresponding response functions are given by blue lines in Fig. 14. We see
that the aforementioned $1\rightarrow 2$ splitting mechanism now give rise to
a non-trivial modification of the valence $u$-quark sPDF after the first two
interactions. As a consequence, the modified valence $u$-quark sPDF can
possess several minima and maxima.

$x_{1}$ | $x_{2}$ |  $N_{u_{v2}}$ P  YTHIA |  $N_{u_{v2}}$ Naive
---|---|---|---
$10^{-6}$ | $10^{-4}$ | 0.025 | 2.006
$10^{-3}$ | $10^{-4}$ | 0.011 | 2.006
$10^{-1}$ | $10^{-4}$ | 0.007 | 2.005
0.2 | $10^{-4}$ | 0.006 | 2.005
0.4 | $10^{-4}$ | 0.005 | 1.997
0.8 | $10^{-4}$ | 0.002 | 1.708
Table 9: Integration over $R_{uuu}$ response function with respect to $x_{3}$ at fixed $x_{1}$, $x_{2}$. In the ideal situation when the GS sum rules are perfectly satisfied $N_{u_{v2}}=0$. $x_{1}$ | $x_{2}$ |  $N_{u_{v2}}$ P  YTHIA |  $N_{u_{v2}}$ Naive
---|---|---|---
$10^{-6}$ | $10^{-4}$ | 2.019 | 2.006
$10^{-3}$ | $10^{-4}$ | 2.005 | 2.006
$10^{-1}$ | $10^{-4}$ | 2.001 | 2.005
0.2 | $10^{-4}$ | 2.000 | 2.005
0.4 | $10^{-4}$ | 1.999 | 1.997
0.8 | $10^{-4}$ | 1.995 | 1.708
Table 10: Integration over $R_{u\bar{u}u}$ response function with respect to
$x_{3}$ at fixed $x_{1}$, $x_{2}$. In the ideal situation when the GS sum
rules are perfectly satisfied $N_{u_{v2}}=2$.

a)

b)

Figure 14: The responses of the valence $u$-quark sPDF computed with
asymmetric (blue) and symmetric (green) tPDFs: a) $uuu$ flavour combination
and b) $u\bar{u}u$ flavour combination.

According to Eq. (74) the integrations over the $uuu$ and $u\bar{u}u$ response
functions should yield

$\displaystyle\int\limits^{1-x_{1}-x_{2}}_{0}\frac{dx_{3}}{x_{3}}~{}R_{uuu}(x_{1},x_{2},x_{3},Q)$
$\displaystyle=$ $\displaystyle
N_{u_{v}}-\delta_{uu}-\delta_{uu}+\delta_{\bar{u}u}+\delta_{\bar{u}u}=0,$
(110)
$\displaystyle\int\limits^{1-x_{1}-x_{2}}_{0}\frac{dx_{3}}{x_{3}}~{}R_{u\bar{u}u}(x_{1},x_{2},x_{3},Q)$
$\displaystyle=$ $\displaystyle
N_{u_{v}}-\delta_{uu}-\delta_{u\bar{u}}+\delta_{\bar{u}u}+\delta_{\bar{u}\bar{u}}=2.$
(111)

We see that for most of the cases the results of the numerical integration
given in Table 9 and Table 10 agree with Eq. (110) and Eq. (111) up to a few
per mille.

Very similar results are obtained for the other flavour combinations. For
example, for the $s\bar{s}s$ case the number sum rule reads as follows:

$\displaystyle\int\limits^{1-x_{1}-x_{2}}_{0}\frac{dx_{3}}{x_{3}}~{}R_{s\bar{s}s}(x_{1},x_{2},x_{3},Q)$
$\displaystyle=$ $\displaystyle
N_{s_{v}}-\delta_{ss}-\delta_{s\bar{s}}+\delta_{\bar{s}s}+\delta_{\bar{s}\bar{s}}=0.$
(112)

The response functions for this case are plotted in Fig. 15, and numerically
we find that the integral on the left hand side of Eq. (112) does give zero up
to deviations of order 0.001, for the same $x$ values as were probed in the
$uuu$ and $u\bar{u}u$ cases. In Fig. 15 we see once again that the companion
quarks produced via the $g\rightarrow q\,\bar{q}$ splitting introduce
additional strange and anti-strange quarks to the beam remnant, leading to
peaks and troughs appearing in the response function.

Figure 15: Same as in Fig. 14 but for the $s\bar{s}s$ flavour combination.

Now let us move to discuss symmetrized tPDFs – we recall that the tPDFs should
be symmetric when the flavour, $x$ fraction and scale arguments for parton $i$
are exchanged with those of parton $j$ (where $i,j\in\left\\{1,3\right\\}$),
and that the (equal-scale) tPDFs should satisfy the sum rules when integrating
over any of the partons, not just the final one. Similarly to Eq. (30) we
define symmetrized P  YTHIA tPDFs as

$\displaystyle T^{\rm
sym}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3},Q)=\frac{1}{3!}\sum\limits_{\\{1,2,3\\}}T_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3},Q),$
(113)

where the sum on the right hand side of Eq. (113) runs over all partonic
permutations. For the sake of simplicity we omit the factorization scale label
in what follows (recall we take $Q=91$ GeV in the numerics). The momentum rule
Eq. (62), therefore, becomes

$\displaystyle\frac{1}{1-x_{1}-x_{2}}\sum_{j_{3}}\int\limits^{1-x_{1}-x_{2}}_{0}dx_{3}x_{3}~{}\frac{T^{\rm
sym}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})}{D^{\rm
sym}_{j_{1}j_{2}}(x_{1},x_{2})}$ $\displaystyle\overset{?}{=}1.$ (114)

The numerical check of Eq. (114) for the $j_{1}=j_{2}=u$ case is given in
Table 11. We see that when both $x_{1}$ and $x_{2}$ are small the sum rule is
pretty well satisfied (with deviations of the order of only $4\%$), but the
violations become larger when one $x$ fraction approaches $1$, reaching the
$70\%$ level when $x_{1}=0.8,x_{2}=10^{-4}$.

$x_{1}$ | $x_{2}$ | $j_{1}$ | $j_{2}$ | P  YTHIA tPDFs | P  YTHIA tPDFs sym. | Naive tPDFs
---|---|---|---|---|---|---
$10^{-6}$ | $10^{-4}$ | $u$ | $u$ | 0.996 | 0.965 | 0.996
$10^{-3}$ | $10^{-4}$ | $u$ | $u$ | 0.996 | 0.967 | 0.996
$10^{-1}$ | $10^{-4}$ | $u$ | $u$ | 1.008 | 0.998 | 1.106
0.2 | $10^{-4}$ | $u$ | $u$ | 1.010 | 1.029 | 1.244
0.4 | $10^{-4}$ | $u$ | $u$ | 1.011 | 1.117 | 1.649
0.8 | $10^{-4}$ | $u$ | $u$ | 1.011 | 1.719 | 4.119
Table 11: Test of the momentum sum rule for the tPDFs. In the ideal situation
when the momentum sum rule is perfectly satisfied the number in each cell
should be equal to unity.

Similarly to the case of the symmetrized dPDFs we define the symmetrized
response functions as

$\displaystyle R^{\ast}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})\equiv
x_{3}\,\frac{T^{\rm sym}_{j_{1}j_{2}j_{3}}(x_{1},x_{2},x_{3})-T^{\rm
sym}_{j_{1}j_{2}\bar{j_{3}}}(x_{1},x_{2},x_{3})}{D^{\rm
sym}_{j_{1}j_{2}}(x_{1},x_{2})},$ (115)

where $D^{\rm sym}_{j_{1}j_{2}}$ are symmetrized P  YTHIA dPDFs defined by Eq.
(30). The numerical checks for $R^{\ast}_{uuu}$ and $R^{\ast}_{u\bar{u}u}$ are
given in Tables 12 and 13 and in Fig.14. We see that for $R^{\ast}_{uuu}$ and
$R^{\ast}_{u\bar{u}u}$, the number sum rules are preserved reasonably well –
at most of the $x_{1},x_{2}$ points tested, we have deviations of order
$0.2-0.3$. At certain points we see larger deviations up to $\sim 0.5-0.6$ –
for the $uuu$ case, this occurs when one of the $x$ values gets close to $1$,
but interestingly, for the $u\bar{u}u$ case, this occurs when the $x$
fractions are at the smallest values tested ($x_{1}=10^{-6},x_{2}=10^{-4}$).

$x_{1}$ | $x_{2}$ |  $N_{u_{v2}}$ P  YTHIA |  $N_{u_{v2}}$ P  YTHIA sym. |  $N_{u_{v2}}$ Naive
---|---|---|---|---
$10^{-6}$ | $10^{-4}$ | 0.025 | 0.108 | 2.006
$10^{-3}$ | $10^{-4}$ | 0.011 | -0.276 | 2.006
$10^{-1}$ | $10^{-4}$ | 0.007 | -0.232 | 2.005
0.2 | $10^{-4}$ | 0.006 | -0.242 | 2.005
0.4 | $10^{-4}$ | 0.005 | -0.317 | 1.997
0.8 | $10^{-4}$ | 0.002 | -0.589 | 1.708
Table 12: Same as in Table 9 but for the symmetrized $R^{\ast}_{uuu}$ response function. In the ideal situation when the GS sum rules are perfectly satisfied $N_{u_{v2}}=0$. $x_{1}$ | $x_{2}$ |  $N_{u_{v}}$ P  YTHIA |  $N_{u_{v}}$ P  YTHIA sym. |  $N_{u_{v}}$ Naive
---|---|---|---|---
$10^{-6}$ | $10^{-4}$ | 2.019 | 2.542 | 2.006
$10^{-3}$ | $10^{-4}$ | 2.005 | 2.154 | 2.006
$10^{-1}$ | $10^{-4}$ | 2.001 | 2.188 | 2.005
0.2 | $10^{-4}$ | 2.000 | 2.189 | 2.005
0.4 | $10^{-4}$ | 1.999 | 2.161 | 1.997
0.8 | $10^{-4}$ | 1.995 | 2.079 | 1.708
Table 13: Same as in Table 10 but for symmetrized $R^{\ast}_{u\bar{u}u}$
response function. In the ideal situation when the GS sum rules are perfectly
satisfied $N_{u_{v2}}=2$.

Whereas the symmetrized tPDFs in Fig. 14 show relatively small violations of
the GS sum rules the problem of the large peaks discussed in Section 3 also
shows for certain tPDFs, for instance for the $s\bar{s}s$ combination in Fig.
15. Here we see a large violation of the number sum rule when $x_{1}$ is very
large, $x_{1}=0.8$, similar to what we saw for the dPDF case in Fig. 7.

##  5 Summary and conclusions

We have demonstrated how one can use the MPI model of the P  YTHIA event
generator to construct asymmetric and symmetric dPDFs which approximately obey
the Gaunt-Stirling (GS) number and momentum sum rules. The asymmetric dPDFs
constructed in this way obey the GS sum rules at a few percent accuracy level
in one of the arguments. We found that the symmetrized dPDFs obey the GS sum
rules at about 20$\%$ accuracy level for $x$ fractions smaller than 0.4. The
concept of a “response function” was introduced for both the momentum and
number sum rules, a function of $x_{1}$ and $x_{2}$ constructed out of dPDFs
and sPDFs that, when integrated over $x_{2}$, should give a constant value
according the sum rules. The precise definitions of these response functions
are given in Eqs. (34) and (43). An interesting observation we made is that
even though only the “areas” under the response function curves in the $x_{2}$
direction are fixed by the sum rules, the overall response functions computed
with P  YTHIA and GS09 dPDFs show very similar behaviour.

We also compared the GS09 and P  YTHIA dPDFs in terms of their predictions for
the double Drell-Yan cross section, focussing in particular on the
distribution in the maximal rapidity separation between the leptons
($\Delta\rm Y$). We observed differences in shape between the GS09 and P
YTHIA $\Delta\rm Y$ distributions, with discrepancies up to the $10\%$ level.
This indicates that although the response function curves are rather similar
between the two models, there are notable differences between the dPDFs
themselves.

Additionally, we generalized the GS sum rules to the case of tPDFs. We
demonstrated that these hold at the bare level using the light-cone wave
function approach, and showed that they also hold at the renormalized level
following the approach of Ref. Diehl:2018kgr . We demonstrated how one can use
the P  YTHIA code to construct the asymmetric and symmetric tPDFs which
approximately obey the corresponding version of the GS sum rules. We argue
that in the absence of available sets of tPDFs one can use P  YTHIA tPDFs as a
good first approximation for phenomenological studies of the TPS phenomena.

Finally, we shall note that both GS09 and P  YTHIA approaches do not take into
account transverse dependence ($\bm{y}$-dependence) between partons, and that
the cross section formulae for DPS or TPS actually involve the quantities that
include this transverse dependence, the double parton distributions (DPDs) or
triple parton distributions (TPDs). Theoretical studies of DPS phenomena
(Blok:2011bu, ; Blok:2013bpa, ; Ryskin:2011kk, ; Ryskin:2012qx, ;
Manohar:2012pe, ; Diehl:2017kgu, ; Diehl:2021gvs, ) have demonstrated that
proper account of the $\bm{y}$-dependence results in contributions with
$1\rightarrow 2$ partonic splittings being numerically more important than in
approaches neglecting transverse dependence of dPDFs, e.g. GS09. Even though
dPDFs and tPDFs do not appear directly in the rigorous treatment of the DPS
and TPS cross section formulae, they are nonetheless linked to the DPDs and
TPDs which do – the former can be obtained from the latter via an appropriate
integral over transverse separation parameter $\bm{y}$. One could, presumably,
extend the P  YTHIA approach to construct DPDs and TPDs approximately
satisfying sum rules, following a similar approach to Ref. (Diehl:2020xyg, ).
The DPDs and TPDs constructed in a such way would be a useful low-scale input
for new tools (e.g. dShower (Cabouat:2019gtm, ; Cabouat:2020ssr, )) that
follow the theoretical rigorous approach to multiple scattering.

## Acknowledgments

The work of OF has received funding from the European Union’s Horizon 2020
research and innovation programme as part of the Marie Skłodowska-Curie
Innovative Training Network MCnetITN3 (grant agreement no. 722104), the
Deutsche Forschungsgemeinschaft (DFG) through the Research Training Group “GRK
2149: Strong and Weak Interactions - from Hadrons to Dark Matter” and the
curiosity-driven grant “Using jets to challenge the Standard Model of particle
physics” from Università di Genova. The work of JRG is supported by the Royal
Society through Grant URF\R1\201500.