aainstitutetext: INFN, Sezione di Genova, Via Dodecaneso 33, I-16146, Genoa,
Italybbinstitutetext: University of Cambridge, Cavendish Laboratory, J.J.
Thomson Avenue, Cambridge, UK

# Correcting for cutoff dependence in backward evolution of QCD parton showers

S. Frixione b B.R. Webber<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Monte Carlo event generators for hard hadronic collisions depend on the
evolution of parton showers backwards from a high-scale subprocess to the
hadronization scale. The evolution is treated as a branching process with a
sequence of resolvable parton emissions. The criterion of resolvability
involves cutoffs that determine the no-emission probability (NEP) for a given
range of the evolution scale. Existing event generators neglect cutoff-
dependent terms in the NEP that, although formally power-suppressed, can have
significant phenomenological effects. We compute such terms and study their
consequences. One important result is that it is not possible for the backward
shower to faithfully reproduce the cutoff-independent parton distribution
functions (PDFs) used to generate it. We show that the computed NEP
corrections mitigate but do not eliminate this problem. An alternative
approach is to use cutoff-dependent PDFs that are consistent with the
uncorrected NEP. Then one must apply cutoff-dependent corrections to hard
subprocess matrix elements. We compute those corrections to the first
nontrivial order for the Drell-Yan process and for Higgs production by gluon
fusion.

##  1 Introduction

Although the precision of predictions of short-distance cross sections using
QCD perturbation theory has greatly increased in recent years, it remains true
that their comparison with experimental data relies to a large extent on
parton-shower based Monte Carlo event generators111For a review see
Buckley:2019kjt . (MCEGs) for the estimation of non-perturbative and
approximate higher-order perturbative effects. The connection between a
measured cross section at a hadron collider and that stemming from the short-
distance subprocess can be described in an inclusive sense using perturbation
theory, factorization theorems and parton distribution functions (PDFs).
However, for the more exclusive description required for the estimation of
experimental effects a reliable MCEG is essential.

A key component of any MCEG for hadronic collisions is a backward parton
shower that links the short-distance subprocess to the incoming hadrons via an
iterative parton branching procedure. For reasons of Monte Carlo efficiency,
the shower starts at the high virtuality scale of the subprocess, with
appropriate parton flavours and momentum fractions, and ends at the lower
scale of hadron formation. For example, in the production of a $Z^{0}$ boson
at leading order, the parton showers should be initiated by a quark-antiquark
pair of equal flavour with invariant mass within the $Z^{0}$ line width. If
the showers were generated forwards from the hadron scale, as is normally done
in PDF evolution, then the efficiency for finding a pair of the same flavour
with an appropriate invariant mass would be unacceptably low.

A special feature of the backward parton shower Sjostrand:1985xi ;
Marchesini:1987cf  is that it must be “guided” by input PDFs, which are
supposed to ensure that the ensemble of parton flavours and momentum fractions
in the shower at any intermediate scale remains consistent with those PDFs.
Compared to forward evolution, this implies modifications to both the
probability of branching as a function of scale, and the distribution of
momentum fractions within each branching. However, the branching process
necessarily involves a sequence of resolvable parton emissions, defined by
some cutoffs, while the PDFs are normally taken from global fits that satisfy
evolution equations Gribov:1972ri ; Altarelli:1977zs ; Dokshitzer:1977sg  that
contain no such cutoffs. This could give rise to systematic biases that, as
far as we are aware, have not been studied so far and are the focus of the
present paper.

In sect. 2 we present a general analysis of PDF evolution equations, not
limited to QCD or any particular perturbative order but suited to the
discussion of issues related to the resolvability of emissions. We pay
particular attention to the ambiguities in the treatment of unresolved and
virtual contributions, and the choices inherent in their resolution. Section 3
examines the backward MC showering process in this framework, in particular
the key concept of the non-emission probability (NEP), which governs the
evolution of the shower in a way supposedly consistent with a given set of
guiding PDFs. We show that neither of the NEP expressions in current use is
formally correct in the presence of cutoffs. However, we find that there is no
fully satisfactory formulation of the NEP as long as the guiding PDFs satisfy
the normal cutoff-independent evolution equations.

Section 4 applies the general approach of the preceding sections to the case
relevant to the most widely-used MCEGs (before any matching or merging),
namely that of leading-order QCD. We show results on the NEP formulations in
current use, the improved expression derived in sect. 3, and their effects on
MC backward evolution. The general conclusion is that, while the improved
expression performs best, all formulations fail to achieve satisfactory
consistency with the guiding cutoff-independent PDFs.

We therefore turn in sect. 5 to an alternative approach, in which the guiding
PDFs obey evolution equations that incorporate the same cutoffs as the
backward parton shower. We show that such PDFs can be made exactly consistent
with the constraints of flavour and momentum conservation, and verify that the
corresponding NEP ensures consistency between the guiding PDFs and the MC
results. Of course, if this approach were implemented, the cutoff-dependent
PDFs would be specific to the cutoffs employed in a particular MCEG, and would
need to be extracted from dedicated global fits. Furthermore, in those fits
the factorization of PDFs and short-distance cross sections implies that the
latter will also be modified by cutoff-dependent terms. In sect. 6 we derive a
general expression for these cutoff corrections to the first nontrivial order
in QCD, and illustrate its application to the processes of lepton pair and
Higgs boson production. Finally in sect. 7 we summarize our main results and
conclusions.

Appendix A contains a more detailed discussion of the relation between the
evolution equations and the backward MC process. Appendix B presents a toy
model in which all emissions are unresolvable, designed to illuminate the
ambiguities and difficulties in defining the NEP.

##  2 PDF evolution equations

In this section, we recast the evolution equations for the PDFs in a form
which is suited to a parton-shower Monte Carlo (MC)-like approach. The
evolution variable $\mu^{2}$ has canonical dimension of mass squared; its
specific nature is not relevant here.

The starting point is the evolution equations Gribov:1972ri ; Altarelli:1977zs
; Dokshitzer:1977sg , which we write as follows:

$\frac{\partial F(x)}{\partial\log\mu^{2}}={\mathbb{O}}\otimes_{x}F\,,$ (1)

where

${\mathbb{O}}\otimes_{x}F=\int_{0}^{1}\frac{dz}{z}\,{\mathbb{O}}(z)\,F(x/z)\,,$
(2)

with the understanding that $F(x/z)=0$ for $z<x$. We assume to be working in a
$d\equiv 1+2(N_{u}+N_{d})$-dimensional flavour space, where $F$ is a column
vector whose $d$ individual components $(F)_{i}$ are the PDFs $f_{i}$ of the
various partons, and ${\mathbb{O}}$ is a $d\times d$ matrix, whose elements in
the $\overline{\rm MS}$ factorisation scheme are the splitting kernels;
additional terms are present in a non-$\overline{\rm MS}$ factorisation
scheme. Both $F$ and ${\mathbb{O}}$ are $x$-space objects, that depend on
$\mu^{2}$ as well; in the notation, either or both dependences may be included
explicitly or understood. It is safe to assume, at least up to the NLO and in
any factorisation scheme, that the most general form of ${\mathbb{O}}$ is:

${\mathbb{O}}(z)=\left[{\mathbb{A}}(z)\right]_{+}+{\mathbb{B}}\,\delta(1-z)+{\mathbb{C}}(z)\,,$
(3)

where

$\left({\mathbb{A}}(z)\right)_{ij}=\delta_{ij}A_{i}(z)\,,\;\;\;\;\;\;\;\;\left({\mathbb{B}}\right)_{ij}=\delta_{ij}B_{i}\,,\;\;\;\;\;\;\;\;\left({\mathbb{C}}(z)\right)_{ij}=C_{ij}(z)\,,$
(4)

with $1\leq i,j\leq d$ the parton indices; note that ${\mathbb{C}}$ is in
general non-diagonal. $A_{i}(z)$ and $C_{ij}(z)$ are regular functions of $z$,
and $B_{i}$ are constants in $z$; all of them depend on $\mu^{2}$. Typically,
$A_{i}(z)$ diverges when $z\to 1$, and the plus prescription in eq. (3)
regularises that divergence; any divergence at $z\to 0$ is not regularised.
Equation (3) encompasses one of the forms in which the NLO splitting kernels
are usually written, namely:

$\sum_{k=0}^{1}\left(\frac{\alpha(\mu^{2})}{2\pi}\right)^{k+1}{\mathbb{P}}^{[k]}(z)=\widetilde{{\mathbb{A}}}(z)\left[\frac{1}{1-z}\right]_{+}+\widetilde{{\mathbb{B}}}\,\delta(1-z)+{\mathbb{C}}(z)\,,$
(5)

with $\widetilde{{\mathbb{A}}}(z)$ finite at $z=1$. Indeed, it is a matter of
applying the definition of the plus distribution to show that, when the
following relationships

$\displaystyle{\mathbb{A}}(z)=\frac{\widetilde{{\mathbb{A}}}(z)}{1-z}\,,\;\;\;\;\;\;\;\;{\mathbb{B}}=\widetilde{{\mathbb{B}}}+\int_{0}^{1}dz\,\frac{\widetilde{{\mathbb{A}}}(z)-\widetilde{{\mathbb{A}}}(1)}{1-z}$
(6)

hold, then the r.h.s.’s of eqs. (3) and (5) are identical to one another.

In order to proceed, we introduce the following symbols:

$\Theta_{ij,z}^{\scriptscriptstyle\rm
IN}=\Theta(\epsilon_{ij}^{{\scriptscriptstyle\rm
L}}<z<1-\epsilon_{ij}^{{\scriptscriptstyle\rm
U}})\,,\;\;\;\;\;\;\;\;\Theta_{ij,z}^{\scriptscriptstyle\rm OUT}\equiv
1-\Theta_{ij,z}^{\scriptscriptstyle\rm
IN}=\Theta(z<\epsilon_{ij}^{{\scriptscriptstyle\rm
L}})+\Theta(z>1-\epsilon_{ij}^{{\scriptscriptstyle\rm U}})\,,$ (7)

with $1-\epsilon_{ij}^{{\scriptscriptstyle\rm
L}}-\epsilon_{ij}^{{\scriptscriptstyle\rm U}}>0$,
$\epsilon_{ij}^{{\scriptscriptstyle\rm L}}>0$, and
$\epsilon_{ij}^{{\scriptscriptstyle\rm U}}>0$. The parameters
$\epsilon_{ij}^{{\scriptscriptstyle\rm L}}$ and
$\epsilon_{ij}^{{\scriptscriptstyle\rm U}}$ are flavour (and possibly scale)
dependent cutoffs, which help to define an inner
($\Theta_{ij,z}^{\scriptscriptstyle\rm IN}$) and an outer
($\Theta_{ij,z}^{\scriptscriptstyle\rm OUT}$) region; the former (latter) will
be associated with resolved (unresolved) emissions in the $z$ space for the
branching:

$j(1)\;\longrightarrow\;i(z)+k(1-z)\;\;\;\;\Longleftrightarrow\;\;\;\;P_{ij}(z)\,.$
(8)

Equation (7) then implies that $\epsilon_{ij}^{{\scriptscriptstyle\rm L}}$ and
$\epsilon_{ij}^{{\scriptscriptstyle\rm U}}$ limit from below the fractional
energy of parton $i$ and recoil system $k$, respectively. At the LO, the
recoil system is a parton itself, unambiguously determined by $i$ and $j$, so
that it may be denoted by $k=j\ominus i$. This suggests introducing the
$1+N_{u}+N_{d}$ parameters:

$\epsilon_{g}\,,\;\epsilon_{u}\,,\;\epsilon_{d}\,,\ldots\,,$ (9)

and setting:

$\epsilon_{ij}^{{\scriptscriptstyle\rm
L}}=\epsilon_{i}\,,\;\;\;\;\;\;\;\;\epsilon_{ij}^{{\scriptscriptstyle\rm
U}}=\epsilon_{j\ominus i}\,.$ (10)

This implies that the lower bounds on the fractional energies depend solely on
the individual parton identities, rather than on the splitting types. In
keeping with what has been done so far, the quantities defined in eq. (7) can
be arranged compactly in two matrices,
${\mathbb{T}}_{z}^{\scriptscriptstyle\rm IN}$ and
${\mathbb{T}}_{z}^{\scriptscriptstyle\rm OUT}$, whose elements are:

$\left({\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}\right)_{ij}=\Theta_{ij,z}^{\scriptscriptstyle\rm
IN}\,,\;\;\;\;\;\;\;\;\left({\mathbb{T}}_{z}^{\scriptscriptstyle\rm
OUT}\right)_{ij}=\Theta_{ij,z}^{\scriptscriptstyle\rm OUT}\,.$ (11)

For any function $g(z)$ and pair of parton indices $(i,j)$, we can exploit the
following identity:

$\displaystyle\left[g(z)\right]_{+}$ $\displaystyle=$
$\displaystyle\left[g(z)\,\Theta_{ij,z}^{\scriptscriptstyle\rm
OUT}\right]_{+}+\left[g(z)\,\Theta_{ij,z}^{\scriptscriptstyle\rm
IN}\right]_{+}$ (12) $\displaystyle=$
$\displaystyle\left[g(z)\,\Theta_{ij,z}^{\scriptscriptstyle\rm
OUT}\right]_{+}+\left(g(z)\,\Theta_{ij,z}^{\scriptscriptstyle\rm
IN}\right)+\left(-\int_{0}^{1}d\omega\,g(\omega)\Theta_{ij,\omega}^{\scriptscriptstyle\rm
IN}\right)\delta(1-z)\,,\phantom{aaa}$

and rewrite eq. (3) as follows:

${\mathbb{O}}(z)=\left[{\mathbb{A}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
OUT}\right]_{+}+{\mathbb{A}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}+\overline{{\mathbb{B}}}\,\delta(1-z)+{\mathbb{C}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
OUT}+{\mathbb{C}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm IN}\,,$ (13)

where by $\circ$ we have denoted the element-by-element matrix multiplication,
e.g.:

$\left({\mathbb{A}}\circ{\mathbb{T}}\right)_{ij}=\left({\mathbb{A}}\right)_{ij}\left({\mathbb{T}}\right)_{ij}$
(14)

and:

$\overline{{\mathbb{B}}}={\mathbb{B}}-\int_{0}^{1}dz\,{\mathbb{A}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}\,.$ (15)

For those operators for which eq. (6) holds, eq. (15) can be written in the
equivalent form222We point out that eq. (13) is unchanged. This implies, in
particular, that $\widetilde{{\mathbb{A}}}(z)/(1-z)$ is inside the plus
prescription in the first term on the r.h.s..:

$\overline{{\mathbb{B}}}=\widetilde{{\mathbb{B}}}+\int_{0}^{1}dz\,\frac{\widetilde{{\mathbb{A}}}(z)-\widetilde{{\mathbb{A}}}(1)}{1-z}\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
OUT}-\int_{0}^{1}dz\,\frac{\widetilde{{\mathbb{A}}}(1)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}}{1-z}\,.$ (16)

With eq. (13), we can write the evolution equations as follows:

$\frac{\partial
F(x)}{\partial\log\mu^{2}}={\mathbb{O}}\otimes_{x}F={\mathbb{W}}\left[F\right](x)+{\mathbb{Z}}\left[F\right](x)+\overline{{\mathbb{B}}}F(x)\,,$
(17)

where:

$\displaystyle{\mathbb{W}}\left[F\right](x)$ $\displaystyle=$
$\displaystyle\Big{(}\big{[}{\mathbb{A}}\circ{\mathbb{T}}^{\scriptscriptstyle\rm
OUT}\big{]}_{+}+{\mathbb{C}}\circ{\mathbb{T}}^{\scriptscriptstyle\rm
OUT}\Big{)}\otimes_{x}F\,,$ (18) $\displaystyle{\mathbb{Z}}\left[F\right](x)$
$\displaystyle=$
$\displaystyle\Big{(}\big{[}{\mathbb{A}}+{\mathbb{C}}\big{]}\circ{\mathbb{T}}^{\scriptscriptstyle\rm
IN}\Big{)}\otimes_{x}F\,.$ (19)

With the scalar functions introduced in eq. (4), these are ($(F)_{i}=f_{i}$):

$\displaystyle\Big{(}\big{[}{\mathbb{A}}\circ{\mathbb{T}}^{\scriptscriptstyle\rm
OUT}\big{]}_{+}\otimes_{x}F\Big{)}_{i}$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\int_{0}^{1}dz\,A_{i}(z)\Theta_{ii,z}^{\scriptscriptstyle\rm
OUT}\left[\frac{\Theta(z\geq
x)}{z}f_{i}\left(\frac{x}{z}\right)-f_{i}(x)\right],\phantom{aa}$ (20)
$\displaystyle\Big{(}\big{(}{\mathbb{A}}\circ{\mathbb{T}}^{\scriptscriptstyle\rm
IN}\big{)}\otimes_{x}F\Big{)}_{i}$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\int_{0}^{1}dz\,A_{i}(z)\Theta_{ii,z}^{\scriptscriptstyle\rm
IN}\,\frac{\Theta(z\geq x)}{z}f_{i}\left(\frac{x}{z}\right),$ (21)
$\displaystyle\Big{(}\big{(}{\mathbb{C}}\circ{\mathbb{T}}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}\big{)}\otimes_{x}F\Big{)}_{i}$
$\displaystyle\\!\\!=\\!\\!$
$\displaystyle\sum_{j}\int_{0}^{1}dz\,C_{ij}(z)\Theta_{ij,z}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}\,\frac{\Theta(z\geq
x)}{z}f_{j}\left(\frac{x}{z}\right).$ (26)

By construction, the r.h.s. of eq. (17) is cutoff-independent. One can show
that the contribution to $({\mathbb{W}}\left[F\right](x))_{i}$ from the
splitting $j\to ik$ is power-suppressed when $x<1-\epsilon_{j\ominus i}$.
Overall, $({\mathbb{W}}\left[F\right](x))_{i}$ cannot be power-suppressed when
$x>1-\epsilon$, with $\epsilon=\min_{j}\epsilon_{j\ominus i}$, because in that
region $({\mathbb{Z}}\left[F\right](x))_{i}=0$, since for such $x$ values one
has $\Theta_{ij,z}^{\scriptscriptstyle\rm IN}\Theta(z\geq x)=0$ for any $z$
and $j$. Therefore, the cutoff dependence of ${\mathbb{W}}\left[F\right](x)$
must cancel that of the $\overline{{\mathbb{B}}}F(x)$ term, which is in
general logarithmic (see e.g. eq. (15)). From a physical viewpoint,
${\mathbb{Z}}\left[F\right]$ describes resolved (owing to
${\mathbb{T}}^{\scriptscriptstyle\rm IN}$) real emissions with
$\max(x,\epsilon)\leq z\leq 1-\epsilon$, while $\overline{{\mathbb{B}}}F(x)$
describes virtual emissions (being proportional to $F(x)$). The term
${\mathbb{W}}\left[F\right]$ is a remainder333This distinction between
${\mathbb{Z}}$ and ${\mathbb{W}}$ is not entirely precise, owing to the
possible flavour dependence of the cutoffs, which implies that certain
kinematical configurations are resolvable only for certain types of
branchings. The underpinning physical picture is nevertheless correct. that
arises from the fact that the kernels of the evolution equations are not
ordinary functions, but distributions that involve subtractions; from a
physics viewpoint, it may be associated with branchings resolvable in the
scale but not in Bjorken $x$.

With MC applications in mind, eq. (17) can be further manipulated by writing:

$\overline{{\mathbb{B}}}(\mu^{2})=\frac{\mu^{2}}{{\mathbb{S}}(\mu^{2})}\,\frac{\partial{\mathbb{S}}(\mu^{2})}{\partial\mu^{2}}\,,$
(27)

with

$\left({\mathbb{S}}(\mu^{2})\right)_{ij}=\delta_{ij}S_{i}(\mu^{2})\,,\;\;\;\;\;\;\;\;\frac{1}{{\mathbb{S}}}\equiv\left({\mathbb{S}}\right)^{-1}\;\;\;\Longrightarrow\;\;\;\left(\frac{1}{{\mathbb{S}}}\right)_{ij}=\delta_{ij}\frac{1}{S_{i}}\,,$
(28)

and

$S_{i}(\mu^{2})=\exp\left[\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\overline{B}_{i}(\kappa^{2})\right]\;\;\;\;\Longleftrightarrow\;\;\;\;{\mathbb{S}}(\mu^{2})=\exp\left[\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\overline{{\mathbb{B}}}(\kappa^{2})\right].$
(29)

In other words, ${\mathbb{S}}$ is a diagonal matrix that collects the Sudakov
form factors. As such, it may seem that the sign in the exponent in eq. (29)
is the opposite w.r.t. the standard one, but in fact this is not the case, as
can be understood from eq. (15). With eq. (27), eq. (17) can be cast as
follows:

$\frac{\partial}{\partial\mu^{2}}\left(\frac{1}{{\mathbb{S}}(\mu^{2})}F(\mu^{2})\right)=\frac{1}{\mu^{2}\,{\mathbb{S}}(\mu^{2})}\,\Big{(}{\mathbb{W}}\left[F\right](\mu^{2})+{\mathbb{Z}}\left[F\right](\mu^{2})\Big{)}\,,$
(30)

which can be put in an integrated form, thus:

$F(\mu^{2})=\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,F(\mu_{0}^{2})+\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\kappa^{2})}\,\Big{(}{\mathbb{W}}\left[F\right](\kappa^{2})+{\mathbb{Z}}\left[F\right](\kappa^{2})\Big{)}\,,$
(31)

or, alternatively, thus444The r.h.s. of eq. (32) features the multiplication
of two column vectors, which is meant as an element-by-element multiplication.
Since no confusion is possible with the multiplications that feature the
transpose of a column vector, no special symbol has been introduced here.:

$\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,\frac{F(\mu_{0}^{2})}{F(\mu^{2})}=\exp\left[-\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{1}{F(\kappa^{2})}\,\Big{(}{\mathbb{W}}\left[F\right](\kappa^{2})+{\mathbb{Z}}\left[F\right](\kappa^{2})\Big{)}\right]\,.$
(32)

We stress again that eqs. (31) and (32) are fully equivalent to eq. (17) but,
being in an integrated form, they also include the information on the initial
conditions ($F(\mu_{0}^{2})$). In turn, they are all equivalent to the
original evolution equation, eq. (1). Thus, in spite of the fact that they
feature cutoff-dependent kernels (${\mathbb{B}}$, ${\mathbb{Z}}$, and
${\mathbb{W}}$), the PDFs that solve them are cutoff-independent. In fact, if
one were interested only in determining the PDFs, the solution of eq. (1)
(best obtained in Mellin space) would be much more straightforward than that
of eqs. (31) or (32). The primary interest of the latter equations is in the
fact that they are expressed in terms of the same quantities that are used in
initial-state parton showers; as such, they can be regarded as giving
consistency conditions among these quantities that initial-state parton
showers (which assume knowledge of the PDFs) must respect. We shall show later
that, in the context of the current approaches used in MCs, this is not quite
the case.

###  2.1 Ambiguities and choices

One of the ingredients of the manipulation of the PDF evolution equations is
the definition of the Sudakov form factors. We point out that, by means of eq.
(27), we have defined them by exponentiating the entire virtual term that
appears in eq. (17). This is not mandatory, and in fact it may lead to
problems (see e.g. sect. 5). In the context of a more flexible approach, we
start by writing the rightmost term on the r.h.s. of eq. (17) as follows:

$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}F(x)+\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm IN}F(x)\,,$ (33)

for any two quantities $\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm OUT}$
and $\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm IN}$ such that:

$\overline{{\mathbb{B}}}=\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}+\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm IN}\,.$ (34)

Then, we define the Sudakov factors thus

${\mathbb{S}}(\mu^{2})=\exp\left[\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
IN}(\kappa^{2})\right],$ (35)

rather than with eq. (29), and we include the contribution from
$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm OUT}$ in the ${\mathbb{W}}[F]$
functional. In order to do that, and also in view of future use (see also
appendix A), it turns out to be convenient to introduce the two evolution
operators:

$\displaystyle{\mathbb{O}}^{\scriptscriptstyle\rm OUT}(z)$ $\displaystyle=$
$\displaystyle\left[{\mathbb{A}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
OUT}\right]_{+}+\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}\,\delta(1-z)+{\mathbb{C}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
OUT}\,,$ (36) $\displaystyle{\mathbb{O}}^{\scriptscriptstyle\rm IN}(z)$
$\displaystyle=$
$\displaystyle{\mathbb{A}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}+\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
IN}\,\delta(1-z)+{\mathbb{C}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}\,,$ (37)

which, loosely speaking, account for emissions in the outer (unresolved) and
inner (resolved) regions, respectively. By construction (see eq. (13)):

${\mathbb{O}}(z)={\mathbb{O}}^{\scriptscriptstyle\rm
OUT}(z)+{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\;\;\;\;\Longrightarrow\;\;\;\;\frac{\partial
F(x)}{\partial\log\mu^{2}}={\mathbb{O}}^{\scriptscriptstyle\rm
OUT}\otimes_{x}F+{\mathbb{O}}^{\scriptscriptstyle\rm IN}\otimes_{x}F\,,$ (38)

and

$\displaystyle{\mathbb{W}}\left[F\right](x)$ $\displaystyle=$
$\displaystyle{\mathbb{O}}^{\scriptscriptstyle\rm OUT}\otimes_{x}F$ (39)
$\displaystyle=$
$\displaystyle\Big{(}\big{[}{\mathbb{A}}\circ{\mathbb{T}}^{\scriptscriptstyle\rm
OUT}\big{]}_{+}+{\mathbb{C}}\circ{\mathbb{T}}^{\scriptscriptstyle\rm
OUT}\Big{)}\otimes_{x}F+\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}F(x)\,,$ (40)
$\displaystyle{\mathbb{Z}}\left[F\right](x)+\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
IN}F(x)$ $\displaystyle=$ $\displaystyle{\mathbb{O}}^{\scriptscriptstyle\rm
IN}\otimes_{x}F\,.$ (41)

As was anticipated, owing to eq. (34) the expression of ${\mathbb{W}}[F]$ in
eq. (40) is in general not the same as that in eq. (18), while that of
${\mathbb{Z}}[F]$ is still given by eq. (19). The crucial thing is that, by
taking into account the redefinition of the Sudakov factor and of the
${\mathbb{W}}[F]$ functional, the integrated form of the evolution equation is
still given by eq. (31) or eq. (32).

While eq. (33) is so far largely arbitrary, given the interpretation of
${\mathbb{W}}$ it is wise to require that:

$\lim_{\epsilon\to 0}\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}=0\,,\;\;\;\;\;\;\;\;\epsilon=\\{\epsilon_{ij}^{{\scriptscriptstyle\rm
L}},\epsilon_{ij}^{{\scriptscriptstyle\rm U}}\\}_{ij}\,.$ (42)

This constraint implies that the Sudakov form factors that one would obtain by
choosing different
$\overline{{\mathbb{B}}}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}$ would differ from one another by
terms suppressed by powers of the cutoffs. In terms of the quantities that
appear in eq. (3), the above can be rewritten by exploiting eq. (15), thus:

$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}={\mathbb{B}}^{\scriptscriptstyle\rm
OUT}\,,\;\;\;\;\;\;\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
IN}={\mathbb{B}}^{\scriptscriptstyle\rm
IN}-\int_{0}^{1}dz\,{\mathbb{A}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}\,,$ (43)

with

${\mathbb{B}}={\mathbb{B}}^{\scriptscriptstyle\rm
OUT}+{\mathbb{B}}^{\scriptscriptstyle\rm
IN}\,,\;\;\;\;\;\;\;\;\lim_{\epsilon\to 0}{\mathbb{B}}^{\scriptscriptstyle\rm
OUT}=0\,.$ (44)

We stress that associating the entire second term on the r.h.s. of eq. (15)
with $\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm IN}$ is merely a sensible
choice, but a choice nevertheless. For example, we could have associated a
$(1-\epsilon)$ fraction of it with
$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm IN}$, and the remaining
$\epsilon$ fraction with $\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}$. In the following, we shall not exploit this option, and always employ
eqs. (43) and (44), so that the flexibility in choosing
$\overline{{\mathbb{B}}}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}$ will be entirely controlled by
the choice of
${\mathbb{B}}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}$.

In the cases where eq. (6) holds, eqs. (36) and (37) become:

$\displaystyle{\mathbb{O}}^{\scriptscriptstyle\rm OUT}(z)$ $\displaystyle=$
$\displaystyle\left[\frac{\widetilde{{\mathbb{A}}}(z)}{1-z}\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
OUT}\right]_{+}+\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}\,\delta(1-z)+{\mathbb{C}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
OUT}\,,$ (45) $\displaystyle{\mathbb{O}}^{\scriptscriptstyle\rm IN}(z)$
$\displaystyle=$
$\displaystyle\frac{\widetilde{{\mathbb{A}}}(z)}{1-z}\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}+\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
IN}\,\delta(1-z)+{\mathbb{C}}(z)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}\,,$ (46)

where, by taking eq. (16) into account:

$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}=\widetilde{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}+\int_{0}^{1}dz\,\frac{\widetilde{{\mathbb{A}}}(z)-\widetilde{{\mathbb{A}}}(1)}{1-z}\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
OUT}\,,\;\;\;\;\;\;\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
IN}=\widetilde{{\mathbb{B}}}^{\scriptscriptstyle\rm
IN}-\int_{0}^{1}dz\,\frac{\widetilde{{\mathbb{A}}}(1)\circ{\mathbb{T}}_{z}^{\scriptscriptstyle\rm
IN}}{1-z}\,,$ (47)

with

$\widetilde{{\mathbb{B}}}=\widetilde{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}+\widetilde{{\mathbb{B}}}^{\scriptscriptstyle\rm
IN}\,,\;\;\;\;\;\;\;\;\lim_{\epsilon\to
0}\widetilde{{\mathbb{B}}}^{\scriptscriptstyle\rm OUT}=0\,.$ (48)

Here, the same remark made after eq. (44) applies: namely, the association of
the two rightmost terms of eq. (16) with
$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm OUT}$ and
$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm IN}$, respectively, as is done
in eq. (47) is a choice we shall always adhere to, and for the operators of
this form the flexibility in choosing
$\overline{{\mathbb{B}}}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}$ will be controlled by the choice
of
$\widetilde{{\mathbb{B}}}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}$.

There is an easy way to enforce the conditions in eqs. (44) and (48) that is,
once again, quite arbitrary, but that allows an easy interpretation from a
physical viewpoint, and leads to the Sudakov form factors which are typically
adopted at the LO in QCD (see sect. 4). Namely, one finds functions
$b_{ij}(\omega)$ and $\tilde{b}_{ij}(\omega)$ which are bounded from above and
below, and are such that:

$B_{j}=\int_{0}^{1}d\omega\sum_{i}b_{ij}(\omega)\,,\;\;\;\;\;\;\;\;\widetilde{B}_{j}=\int_{0}^{1}d\omega\sum_{i}\tilde{b}_{ij}(\omega)\,,$
(49)

and defines:

$B^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm
OUT}$\cr}}}}_{j}=\int_{0}^{1}d\omega\sum_{i}b_{ij}(\omega)\,\Theta_{ij,\omega}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm
OUT}$\cr}}}}\,,\;\;\;\;\;\;\;\;\widetilde{B}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm
OUT}$\cr}}}}_{j}=\int_{0}^{1}d\omega\sum_{i}\tilde{b}_{ij}(\omega)\,\Theta_{ij,\omega}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}\,.$ (50)

Note that in the case where eq. (6) holds, ${\mathbb{B}}$ (and therefore the
functions $b_{ij}(\omega)$) need not necessarily be introduced. If one still
finds it convenient to do so (e.g. to use both the form of eq. (3) and that of
eq. (5)), eqs. (6) and (49) imply:

$\int_{0}^{1}d\omega\,b_{ii}(\omega)=\int_{0}^{1}d\omega\,\tilde{b}_{ii}(\omega)+\int_{0}^{1}dz\,\frac{\widetilde{A}_{i}(z)-\widetilde{A}_{i}(1)}{1-z}\,.$
(51)

Clearly, the easiest way to achieve this and be consistent with eqs. (43)–(48)
is to work with a local version of eq. (51), namely:

$b_{ii}(z)=\tilde{b}_{ii}(z)+\frac{\widetilde{A}_{i}(z)-\widetilde{A}_{i}(1)}{1-z}\,.$
(52)

We anticipate that at the LO in QCD the choice of the $b_{ij}(\omega)$ and
$\tilde{b}_{ij}(\omega)$ functions along the lines presented above leads to
Sudakov form factors expressed as integrals of the LO splitting kernels (see
sect. 4 for more details). However, this discussion should render it clear
that, even at the LO in QCD, this is a choice that is not dictated by any
fundamental principle, but by convenience and ease of interpretation.

The separation of the virtual terms in eq. (33) stemming from eq. (34)
encompasses the case where such a separation is not considered. Even after
making a definite choice for the cutoffs, one can continuously pass from one
scenario to the other by means of the replacements555Although in general the
parameter $\lambda$ can be flavour-dependent, for our purposes such a
dependence can be neglected.

$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}\;\longrightarrow\;\lambda\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}\,,\;\;\;\;\;\;\;\;\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
IN}\;\longrightarrow\;\overline{{\mathbb{B}}}-\lambda\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}\,,$ (53)

with $0\leq\lambda\leq 1$ in all quantities that feature a dependence on
$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm OUT}$ and/or
$\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm IN}$.

##  3 Monte Carlo backward evolution

When an MC generates initial-state parton showers, the PDFs are thought to be
given. They are employed to “guide” the backward evolution, and thus increase
the efficiency of the latter. One usually assumes that consistency demands
that the longitudinal momentum left after all branchings have occurred be
distributed according to the given PDFs (this identification holds in a
statistical sense; it is exact only after an infinite number of showers have
been carried out). However, since only resolved branchings (i.e., those with
$\epsilon<z<1-\epsilon$) can be generated, the identification above can be
true only in the resolved region. In fact, as we shall show, even in the
resolved region MCs are generally not able to reconstruct the PDFs.
Ultimately, this arises from the fact that the PDF evolution equations are
expressed as convolution integrals, and thus the derivative w.r.t. the scale
of the PDF at a given $x$ receives contributions from all $z$’s, with $x\leq
z\leq 1$. In other words, the unresolved region feeds into the resolved region
as well as itself. This is unavoidable: for PDF evolution, the separation
between the resolved and unresolved regions is totally arbitrary, and has no
bearing on the final form of the PDFs.

Conversely, MCs cannot function without a clear separation between resolved
and unresolved regions, i.e. without the introduction of cutoffs. As is well
known, this leads to the possibility of generating showers by means of an
iterative Markovian process, one of whose key ingredients is the inversion of
the so-called non-emission probability (NEP henceforth), which gives one the
scale at which the next parton branching occurs. The usual argument adopted
for deriving the NEP associated with initial-state emissions exploits a
partonic picture of the PDFs, whereby these “count” the number of partons at
any given values of the Bjorken $x$ and scale. As was said before, one
identifies the ${\mathbb{Z}}$ and ${\mathbb{W}}$ contributions to the PDF
evolution as associated with resolvable branchings and with branchings
resolvable in $\mu$ but not in $z$, respectively. Therefore, the number of
partons of flavour $i$ that do not undergo branchings of any type in the range
$(\mu_{0}^{2},\mu^{2})$ is equal to:

$\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\,f_{i}(x,\mu_{0}^{2})\,,$ (54)

while that of partons that either do not branch, or branch in a manner
unresolvable in $x$, is equal to:

$\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\,f_{i}(x,\mu_{0}^{2})+\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{S_{i}(\mu^{2})}{S_{i}(\kappa^{2})}\,\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(x,\kappa^{2})\,.$
(55)

The NEP is defined as the fraction of partons that do not branch in a
resolvable manner between any two scales. The difference between forward and
backward evolution is simply the reference relative to which that fraction is
measured, because the elementary branching mechanism must not be affected by
the direction of the evolution. For an evolution in the range
$(\mu_{0}^{2},\mu^{2})$, if the evolution is forwards (backwards) the
reference is $f_{i}(x,\mu_{0}^{2})$ ($f_{i}(x,\mu^{2})$). Thus, eq. (55) leads
to the non-emission probabilities for the forward and backward evolution of a
parton of type $i$,

Forward: $\displaystyle\phantom{aaa}{\rm
NEP}_{i}=\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}+\frac{1}{f_{i}(x,\mu_{0}^{2})}\,\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{S_{i}(\mu^{2})}{S_{i}(\kappa^{2})}\,\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(x,\kappa^{2}),$
(56) Backward: $\displaystyle\phantom{aaa}{\rm
NEP}_{i}=\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\,\frac{f_{i}(x,\mu_{0}^{2})}{f_{i}(x,\mu^{2})}+\frac{1}{f_{i}(x,\mu^{2})}\,\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{S_{i}(\mu^{2})}{S_{i}(\kappa^{2})}\,\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(x,\kappa^{2}).\phantom{aaaa}$
(57)

Here we are concerned with the backward case, which will henceforth always be
implied. Then from eqs. (57) and (31) one also obtains:

${\rm
NEP}_{i}=1-\frac{1}{f_{i}(x,\mu^{2})}\,\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{S_{i}(\mu^{2})}{S_{i}(\kappa^{2})}\,\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i}(x,\kappa^{2})\,,$
(58)

consistently with the meaning of the ${\mathbb{Z}}[F]$ functional. In a
backward evolution, which proceeds from larger to smaller scales, starting
from a given $\mu^{2}$ one obtains the “next” scale $\mu_{0}^{2}<\mu^{2}$ by
solving for $\mu_{0}^{2}$ the equation

$r={\rm NEP}_{i}\,,$ (59)

with $0<r<1$ a uniform random number, and $i$ given. After selecting the
branching channel and its momentum fraction, the procedure is iterated until a
$\mu_{0}^{2}$ value is obtained that is smaller than some pre-defined
threshold (the so-called hadronization scale). However, MCEGs do not literally
solve eq. (59), but either Marchesini:1987cf

$r={\rm NEP}_{i}^{({\rm
R})}\equiv\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\,\frac{f_{i}(x,\mu_{0}^{2})}{f_{i}(x,\mu^{2})}\,,$
(60)

or Sjostrand:1985xi

$r={\rm NEP}_{i}^{({\rm
E})}\equiv\exp\left[-\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{1}{f_{i}(x,\kappa^{2})}\,\big{(}{\mathbb{Z}}\left[F\right]\big{)}_{i}(x,\kappa^{2})\right]\,,$
(61)

where the superscript R or E indicates that a ratio or an exponential
approximation for the NEP has been used, respectively. The crucial thing,
which follows directly from the evolution equations as given in eq. (31) or
eq. (32), is that:

$\mbox{if}\phantom{aaaa}{\mathbb{W}}\left[F\right]=0\phantom{aaaa}\mbox{then}\phantom{aaaa}{\rm
NEP}_{i}={\rm NEP}_{i}^{({\rm R})}={\rm NEP}_{i}^{({\rm E})}\,.$ (62)

Thus, in the resolved region the solution of eq. (59) coincides with that of
eqs. (60) or (61) up to terms suppressed by powers of the cutoffs (since in
that region ${\mathbb{W}}$ vanishes with the cutoffs). This is the reason why,
in standard current practice, eq. (60) and eq. (61) are considered equivalent
to one another – effects suppressed by powers of the cutoffs are
systematically neglected. This is in fact a dangerous position to take, given
that MC cutoffs are often not particularly small, their effects can accumulate
over the course of evolution, and there is no other mechanism that forces
${\mathbb{W}}$ to vanish bar the vanishing of the cutoffs.

It is therefore instructive to see how the three NEP expressions considered
above differ from each other when power-suppressed cutoff effects are not
neglected. We start by considering the probability distribution of the scale
of the next branching which, according to eqs. (59), (60), and (61), is given
by the derivative w.r.t. $\mu_{0}^{2}$ of the respective NEP. By direct
computation, and by employing the evolution equations (31) and (32), we
obtain:

$\displaystyle\frac{\partial}{\partial\log\mu_{0}^{2}}\,{\rm
NEP}_{i}=\frac{1}{f_{i}(x,\mu^{2})}\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\,\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i}(x,\mu_{0}^{2})\,,$
(63) $\displaystyle\frac{\partial}{\partial\log\mu_{0}^{2}}\,{\rm
NEP}_{i}^{({\rm
R})}=\frac{1}{f_{i}(x,\mu^{2})}\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\Big{[}\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(x,\mu_{0}^{2})+\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i}(x,\mu_{0}^{2})\Big{]},$
(64) $\displaystyle\frac{\partial}{\partial\log\mu_{0}^{2}}\,{\rm
NEP}_{i}^{({\rm E})}=$ (65)
$\displaystyle\phantom{aaaaaa}\frac{1}{f_{i}(x,\mu^{2})}\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\,\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i}(x,\mu_{0}^{2})\,\exp\\!\left[\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{1}{f_{i}(x,\kappa^{2})}\,\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(x,\kappa^{2})\right].$

Equation (63), since it factors out ${\mathbb{Z}}$ that is non-null only for
resolved emissions, shows that ${\rm NEP}_{i}$ is consistent with the
requirement that the NEP be associated with the fraction of partons that do
not branch in a resolvable manner. This may seem to be the case also for ${\rm
NEP}_{i}^{({\rm E})}$, but in fact the exponentiated ${\mathbb{W}}$ term in
eq. (65) introduces a spurious extra cutoff dependence w.r.t. the evolution
generated by means of ${\rm NEP}_{i}$. Finally, in eq. (64) the ${\mathbb{W}}$
and ${\mathbb{Z}}$ contributions are on the same footing: this is because, as
the comparison between eqs. (54) and (55) shows, ${\rm NEP}_{i}^{({\rm R})}$
is actually the NEP for no branchings, regardless whether they are resolved or
unresolved in $x$.

In appendix A we discuss in detail the implications of eqs. (63)–(65) for the
requirement that MC backward evolution allows one to reconstruct the PDFs
given in input to the parton shower. The bottom line is that, in practice,
such a reconstruction always fails. It can be made to formally succeed with
${\rm NEP}_{i}^{({\rm R})}$, while if ${\rm NEP}_{i}$ is adopted one can
reconstruct PDFs where all non-resolved contributions are consistently
neglected; the same is true for ${\rm NEP}_{i}^{({\rm E})}$ if a branching-by-
branching reweighting is applied.

The above suggests that ${\rm NEP}_{i}$ and ${\rm NEP}_{i}^{({\rm E})}$ are
closer to each other than either is to ${\rm NEP}_{i}^{({\rm R})}$. This can
be also seen in another way, by considering the differences between any two of
these quantities. From eqs. (57) and (60) we obtain:

${\rm NEP}_{i}^{({\rm R})}-{\rm
NEP}_{i}=-\frac{S_{i}(\mu^{2})}{f_{i}(\mu^{2})}\,\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{1}{S_{i}(\kappa^{2})}\,\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(\kappa^{2})\equiv{\cal
O}(\alpha_{\scriptscriptstyle S})\,,$ (66)

whereas from eqs. (57) and (61):

$\displaystyle{\rm NEP}_{i}^{({\rm E})}-{\rm NEP}_{i}$
$\displaystyle\phantom{aa}=\frac{S_{i}(\mu^{2})}{f_{i}(\mu^{2})}\,\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}S_{i}(\kappa^{2})}\left(\frac{f_{i}(\mu_{0}^{2})}{f_{i}(\kappa^{2})}\,\frac{S_{i}(\kappa^{2})}{S_{i}(\mu_{0}^{2})}-1\right)\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(\kappa^{2})$
$\displaystyle\phantom{aa=}+\frac{1}{2}\,\frac{f_{i}(\mu_{0}^{2})}{f_{i}(\mu^{2})}\,\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\left(\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{1}{f_{i}(\kappa^{2})}\,\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(\kappa^{2})\right)^{2}+\ldots$
$\displaystyle\phantom{aa}=-\frac{S_{i}(\mu^{2})}{f_{i}(\mu^{2})}\,\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}f_{i}(\kappa^{2})}\int_{\mu_{0}^{2}}^{\kappa^{2}}\frac{d\rho^{2}}{\rho^{2}S_{i}(\rho^{2})}\Big{(}\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(\rho^{2})+\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i}(\rho^{2})\Big{)}\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(\kappa^{2})$
$\displaystyle\phantom{aa=}+\frac{1}{2}\,\frac{f_{i}(\mu_{0}^{2})}{f_{i}(\mu^{2})}\,\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\left(\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{1}{f_{i}(\kappa^{2})}\,\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(\kappa^{2})\right)^{2}+\ldots$
(67) $\displaystyle\phantom{aa}\equiv{\cal O}(\alpha_{\scriptscriptstyle
S}^{2})\,,$ (68)

where the ellipsis represents terms with three or more ${\mathbb{W}}$ terms,
and we have used eqs. (31) and (32). The powers of $\alpha_{\scriptscriptstyle
S}$ in eqs. (66) and (68) stem from having regarded both the PDFs and the
Sudakovs as quantities of perturbative ${\cal O}(1)$, while both
${\mathbb{W}}$ and ${\mathbb{Z}}$ are of ${\cal O}(\alpha_{\scriptscriptstyle
S})$ (see eqs. (18) and (19)).

We finally note that, when power-suppressed effects are not neglected, the
simple probabilistic interpretation upon which MCs rely to perform initial-
state backward evolution may lose validity. In all cases, this can be seen to
come from the fact that ${\mathbb{W}}[F]$ has no definite sign. Thus, from eq.
(32) one sees that ${\rm NEP}_{i}^{({\rm R})}$ is not necessarily monotonic,
and from eq. (31) that ${\rm NEP}_{i}$ is not necessarily positive. In both
cases, this implies that, in some regions of the phase space (typically, at
large Bjorken $x$), these NEPs are actually not cumulative probability
distributions, and therefore that the solution of eq. (60) or eq. (59) may not
exist, or may not be unique. As far as ${\rm NEP}_{i}^{({\rm E})}$ is
concerned, it is positive definite, monotonic, and bounded by one; however, as
was discussed in relation to eq. (65), its physical interpretation is unclear.

We conclude by remarking that, while ${\rm NEP}_{i}$ may turn out to be
negative, it generally is positive. One can in fact turn the requirement that
it be positive into a tool to determine the cutoffs in a physically-meaningful
manner, in the sense of limiting the impact of non-resolvable emissions to an
extent that allows one to recover a probabilistic interpretation.

Another way to approach the problem is to acknowledge the fact that PDFs and
initial-state parton showers are inherently incompatible at some level, and to
construct MC-specific PDFs by means of which all issues are removed ab initio.
This option will be discussed in sect. 5.

##  4 The LO QCD case

The general approach of sects. 2 and 3 can be applied to the case which is
currently the most relevant to MC simulations, namely that where only the LO
evolution kernels are considered. In order to simplify our discussion, we
assume all quarks to be massless, and ignore complications due to the presence
of mass thresholds; thus, we shall not need to specify the individual
flavours.

We denote the LO kernels, i.e. the elements of ${\mathbb{P}}^{[0]}$ in eq.
(5), as follows666Bearing in mind that at the LO there are no $q\bar{q}$
kernels, in the notation we need not distinguish quarks and antiquarks.:

$\displaystyle P_{qq}(z)$ $\displaystyle=$ $\displaystyle
C_{\scriptscriptstyle F}\left(\frac{1+z^{2}}{1-z}\right)_{+}\,,$ (69)
$\displaystyle P_{gq}(z)$ $\displaystyle=$ $\displaystyle
C_{\scriptscriptstyle F}\,\frac{1+(1-z)^{2}}{z}\,,$ (70) $\displaystyle
P_{qg}(z)$ $\displaystyle=$ $\displaystyle T_{\scriptscriptstyle
F}\left(z^{2}+(1-z)^{2}\right)\,,$ (71) $\displaystyle P_{gg}(z)$
$\displaystyle=$ $\displaystyle 2C_{\scriptscriptstyle
A}\left(\frac{z}{(1-z)_{+}}+\frac{1-z}{z}+z(1-z)\right)+\gamma(g)\delta(1-z)\,,$
(72)

with ($N_{\scriptscriptstyle F}=N_{u}+N_{d}$):

$\gamma(g)=\frac{11C_{\scriptscriptstyle A}-4T_{\scriptscriptstyle
F}N_{\scriptscriptstyle F}}{6}\,.$ (73)

We also denote by $\hat{P}_{ij}$ the ordinary function obtained from the
kernels $P_{ij}$ above by discarding the endpoint contributions (i.e. by
turning plus distributions into ordinary functions, and by ignoring
contributions proportional to $\delta(1-z)$). From eqs. (69)–(72) one can read
off the quantities introduced in eqs. (4) and (5), since at this order:

$\big{(}{\mathbb{O}}(z)\big{)}_{ij}=\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\Big{(}{\mathbb{P}}^{[0]}\Big{)}_{ij}\equiv\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\,P_{ij}(z)\,.$ (74)

Thus:

$\displaystyle\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\,A_{q}(z)=C_{\scriptscriptstyle
F}\,\frac{1+z^{2}}{1-z}\,,\;\;\;\;B_{q}=0\,,\;\;\;\;C_{qq}(z)=0\,,$ (75)
$\displaystyle\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\,C_{gq}(z)=C_{\scriptscriptstyle F}\,\frac{1+(1-z)^{2}}{z}\,,$ (76)
$\displaystyle\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\,C_{qg}(z)=T_{\scriptscriptstyle F}\left(z^{2}+(1-z)^{2}\right)\,,$ (77)
$\displaystyle\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\,\widetilde{A}_{g}(z)=2C_{\scriptscriptstyle
A}\,z\,,\;\;\;\;\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\,\widetilde{B}_{g}=\gamma(g)\,,\;\;\;\;\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\,C_{gg}(z)=2C_{\scriptscriptstyle
A}\left(\frac{1-z}{z}+z(1-z)\right).\phantom{aaaa}$ (78)

The case of a quark is straightforward: in view of eq. (75), by making the
simplest choice777At this stage, this is not mandatory. We shall comment
further on this point (see sect. 5.1).
$B_{q}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}=0$, eq. (43) leads to:

$\displaystyle\overline{B}_{q}^{\scriptscriptstyle\rm OUT}$ $\displaystyle=$
$\displaystyle 0\,,$ (79)
$\displaystyle\overline{B}_{q}^{\scriptscriptstyle\rm IN}$ $\displaystyle=$
$\displaystyle-\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\hat{P}_{qq}(z)\Theta_{qq,z}^{\scriptscriptstyle\rm
IN}$ (80) $\displaystyle\equiv$
$\displaystyle-\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\,\frac{1}{2}\left(\hat{P}_{qq}(z)\Theta_{qq,z}^{\scriptscriptstyle\rm
IN}+\hat{P}_{gq}(z)\Theta_{gq,z}^{\scriptscriptstyle\rm IN}\right)\,,$ (81)

with the form in eq. (81) identical to that in eq. (80) thanks to the
$z\leftrightarrow 1-z$ symmetry of both the splitting kernels and their
respective integration limits (owing to eq. (10)). The corresponding quark
Sudakov factor is obtained by inserting
$\overline{B}_{q}^{\scriptscriptstyle\rm IN}$ into eq. (35). The case of the
gluon is slightly more involved. We use:

$\displaystyle\frac{2\pi}{\alpha_{\scriptscriptstyle S}}\,\tilde{b}_{gg}(z)$
$\displaystyle=$ $\displaystyle\frac{C_{\scriptscriptstyle
A}}{1-z}+\frac{C_{\scriptscriptstyle A}}{z}-\frac{1}{2}\,\hat{P}_{gg}(z)\equiv
C_{\scriptscriptstyle A}\Big{(}2-z+z^{2}\Big{)}\,,$ (82)
$\displaystyle\frac{2\pi}{\alpha_{\scriptscriptstyle S}}\,\tilde{b}_{qg}(z)$
$\displaystyle=$ $\displaystyle-\frac{1}{2}\hat{P}_{qg}(z)\,,$ (83)

which indeed satisfies eq. (49) given $\gamma(g)$ of eq. (73) (note that the
sum over flavours includes both quarks and antiquarks). With this, eqs. (47)
and (78) lead to:

$\displaystyle\overline{B}_{g}^{\scriptscriptstyle\rm OUT}$ $\displaystyle=$
$\displaystyle-\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\,C_{\scriptscriptstyle
A}\int_{0}^{1}dz\,z(1-z)\Theta_{gg,z}^{\scriptscriptstyle\rm
OUT}-\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\,\frac{1}{2}\sum_{q,\bar{q}}\int_{0}^{1}dz\,\hat{P}_{qg}(z)\Theta_{qg,z}^{\scriptscriptstyle\rm
OUT}\,,$ (84) $\displaystyle\overline{B}_{g}^{\scriptscriptstyle\rm IN}$
$\displaystyle=$ $\displaystyle-\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\,\frac{1}{2}\int_{0}^{1}dz\left(\hat{P}_{gg}(z)\Theta_{gg,z}^{\scriptscriptstyle\rm
IN}+\sum_{q,\bar{q}}\hat{P}_{qg}(z)\Theta_{qg,z}^{\scriptscriptstyle\rm
IN}\right)\,.$ (85)

The form of eq. (85) stems from exploiting:

$\int_{0}^{1}dz\,\Theta_{gg,z}^{\scriptscriptstyle\rm
IN}\left(-\frac{2C_{\scriptscriptstyle A}}{1-z}+\frac{C_{\scriptscriptstyle
A}}{1-z}+\frac{C_{\scriptscriptstyle A}}{z}\right)=0\,,$ (86)

which is due to the fact that (see eq. (10)):

$\epsilon_{gg}^{{\scriptscriptstyle\rm
L}}=\epsilon_{g}\,,\;\;\;\;\epsilon_{gg}^{{\scriptscriptstyle\rm
U}}=\epsilon_{g}\;\;\;\;\Longrightarrow\;\;\;\;\Theta_{gg,z}^{\scriptscriptstyle\rm
IN}=\Theta(\epsilon_{g}<z<1-\epsilon_{g})\,.$ (87)

If the range in $z$ defined by $\Theta_{gg,z}^{\scriptscriptstyle\rm IN}$ were
not symmetric under $z\leftrightarrow 1-z$,
$\overline{B}_{g}^{\scriptscriptstyle\rm IN}$ would still be well defined, but
eqs. (47), (50), and (82) would not lead to a result solely expressed in terms
of $\hat{P}_{gg}$ for the part proportional to $C_{\scriptscriptstyle A}$.
Finally, we point out that eqs. (81) and (85), which enter the quark and gluon
Sudakov form factors (in the latter case, only when $\lambda=1$),
respectively, have the usual form of the integrals of the splitting kernels
over the inner region. The reader is encouraged to bear in mind that this is a
consequence of several arbitrary choices, which we have outlined in sect. 2.1.

By using the results above, those of eqs. (19) and (40), and the replacement
of eq. (53), we obtain the following after some trivial algebra:

$\displaystyle\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{q}(x)$
$\displaystyle=$ $\displaystyle\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}\frac{dz}{z}\,\Theta(z\geq
x)\left[\Theta_{qq,z}^{\scriptscriptstyle\rm
IN}\hat{P}_{qq}(z)f_{q}\left(\frac{x}{z}\right)+\Theta_{qg,z}^{\scriptscriptstyle\rm
IN}\hat{P}_{qg}(z)f_{g}\left(\frac{x}{z}\right)\right],\phantom{aa}$ (88)
$\displaystyle\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{q}(x)$
$\displaystyle=$ $\displaystyle\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\,\Bigg{\\{}\Theta_{qq,z}^{\scriptscriptstyle\rm
OUT}\hat{P}_{qq}(z)\left[\frac{1}{z}\,f_{q}\left(\frac{x}{z}\right)\Theta(z\geq
x)-f_{q}(x)\right]$ (89)
$\displaystyle\phantom{\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\,\Theta_{\epsilon}}+\frac{1}{z}\,\Theta_{qg,z}^{\scriptscriptstyle\rm
OUT}\hat{P}_{qg}(z)f_{g}\left(\frac{x}{z}\right)\Theta(z\geq
x)\Bigg{\\}}+\lambda\overline{B}_{q}^{\scriptscriptstyle\rm
OUT}f_{q}(x)\,,\phantom{aaa}$

and:

$\displaystyle\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{g}(x)$
$\displaystyle=$ $\displaystyle\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}\frac{dz}{z}\,\Theta(z\geq
x)\left[\Theta_{gg,z}^{\scriptscriptstyle\rm
IN}\hat{P}_{gg}(z)f_{g}\left(\frac{x}{z}\right)+\sum_{q,\bar{q}}\Theta_{gq,z}^{\scriptscriptstyle\rm
IN}\hat{P}_{gq}(z)f_{q}\left(\frac{x}{z}\right)\right],\phantom{aaaaa}$ (90)
$\displaystyle\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{g}(x)$
$\displaystyle=$ $\displaystyle\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\,\Bigg{\\{}\Theta_{gg,z}^{\scriptscriptstyle\rm
OUT}\left[\frac{\hat{P}_{gg}(z)}{z}f_{g}\left(\frac{x}{z}\right)\Theta(z\geq
x)-\frac{2C_{\scriptscriptstyle A}z}{1-z}\,f_{g}(x)\right]$ (91)
$\displaystyle\phantom{\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz}+\frac{\Theta(z\geq
x)}{z}\sum_{q,\bar{q}}\Theta_{gq,z}^{\scriptscriptstyle\rm
OUT}\hat{P}_{gq}(z)f_{q}\left(\frac{x}{z}\right)\Bigg{\\}}+\lambda\overline{B}_{g}^{\scriptscriptstyle\rm
OUT}f_{g}(x)\,.\phantom{aaa}$

We note that, owing to eq. (79), the last term on the r.h.s. of eq. (89) is
null, independent of the value of $\lambda$; the reader must bear in mind that
this is a choice (see footnote 7). If one chooses $\lambda=1$, eqs. (82) and
(83) allow one to rewrite eq. (91) in the seemingly more familiar form:

$\displaystyle\\!\\!\\!\\!\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{g}(x)$
$\displaystyle=$ $\displaystyle\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\,\Bigg{\\{}\Theta_{gg,z}^{\scriptscriptstyle\rm
OUT}\hat{P}_{gg}(z)\left[\frac{1}{z}f_{g}\left(\frac{x}{z}\right)\Theta(z\geq
x)-\frac{1}{2}f_{g}(x)\right]$ (92)
$\displaystyle\phantom{a}+\sum_{q,\bar{q}}\left(\frac{1}{z}\,\Theta_{gq,z}^{\scriptscriptstyle\rm
OUT}\hat{P}_{gq}(z)f_{q}\left(\frac{x}{z}\right)\Theta(z\geq
x)-\frac{1}{2}\,\Theta_{qg,z}^{\scriptscriptstyle\rm
OUT}\hat{P}_{qg}(z)f_{g}\left(x\right)\right)\\!\Bigg{\\}}.\phantom{aaa}$

Again, here a simplification has been made thanks to the fact that the
analogue of eq. (86) holds with $\Theta_{gg,z}^{\scriptscriptstyle\rm
IN}\to\Theta_{gg,z}^{\scriptscriptstyle\rm OUT}$ there, given eq. (87).
Moreover, we observe that this is also a direct consequence of the fact that
the subtraction term in eq. (91) is proportional to $z/(1-z)$, as opposed to
$1/(1-z)$ – the definition of $\overline{{\mathbb{B}}}$ respects the
convention for the plus prescription mentioned in footnote 2.

Equation (92) does not offer any specific advantages w.r.t. eq. (91). In
addition to being valid only when $\lambda=1$, it may seem to feature
uncancelled divergences stemming from the second term in the integrand. In
fact, this is not the case, as one can easily see by regularising the
integral. However, such a regularisation is not practical in the context of
numerical computations. A better alternative is to exploit the
$z\leftrightarrow 1-z$ symmetry of the $\hat{P}_{gg}(z)$ and $\hat{P}_{qg}(z)$
kernels and eq. (87) (as well as its analogue for the $g\to q\bar{q}$
branching), and to obtain a manifestly-finite integral by means of either of
the formal replacements:

$\frac{1}{2}f_{g}\left(x\right)\;\longrightarrow\;\Theta\\!\left(z\geq\frac{1}{2}\right)f_{g}\left(x\right)\,,\;\;\;\;\;\;\;\;\frac{1}{2}f_{g}\left(x\right)\;\longrightarrow\;z\,f_{g}\left(x\right)\,,$
(93)

in the second and fourth terms of the integrand.

###  4.1 Results on backward evolution

We present here some results obtained within the leading-order framework
outlined above. Since our objective is to illustrate issues raised in previous
sections, rather than to perform realistic phenomenology, we consider two
cases of a universal, flavour-independent cutoff
$\epsilon_{ij}^{{\scriptscriptstyle\rm
L}}=\epsilon_{ij}^{{\scriptscriptstyle\rm U}}=\epsilon$. The first is
relatively large and scale-independent, $\epsilon=0.1$, while the second is
slightly more realistic from a parton-shower MC viewpoint, being scale
dependent and defined as $\epsilon=\;\mbox{(2 GeV)}/q$, where $q$ is the mass
scale relevant to the current computation (e.g. in the Sudakov factor of eq.
(35), $q=\sqrt{\kappa^{2}}$). For the leading-order PDFs we adopt the CT18LO
set of ref. Yan:2022pzl ; the argument of $\alpha_{\scriptscriptstyle S}$ is
taken to be a mass-scale squared and, in keeping with ref. Yan:2022pzl , we
have $\alpha_{\scriptscriptstyle S}(m_{\scriptscriptstyle Z}^{2})=0.135$.

Figure 1:  Non-emission probability (NEP) for backward evolution of up quarks
and gluons with $\mu=100$ GeV, according to NEP (57) (black, solid), ${\rm
NEP}^{({\rm R})}$ (60) (blue, dashed) and ${\rm NEP}^{({\rm E})}$ (61) (red,
dotted). The three sets of curves correspond to $x=0.01$ (lowest), 0.1, and
0.5 (highest). The open circles (green) show the NEP computed with cutoff-
dependent PDFs, to be discussed in sect. 5.

Figure 1 shows the resulting true NEP (57) (black, solid) and the
approximations ${\rm NEP}^{({\rm R})}$ (60) (blue, dashed) and ${\rm
NEP}^{({\rm E})}$ (61) (red, dotted), for up-quarks and gluons as a function
of $\mu_{0}$, with $\mu=100$ GeV. In each panel, the three sets of curves are
for $x=0.01$, 0.1, and 0.5, from lowest to highest, respectively. One sees
that, for the range of $\mu_{0}$ shown, ${\rm NEP}^{({\rm E})}$ is closer to
the true NEP than ${\rm NEP}^{({\rm R})}$, as could be anticipated from eqs.
(66) and (68). ${\rm NEP}^{({\rm R})}$ becomes a poorer approximation with
increasing $x$ (because there ${\mathbb{W}}\left[F\right]$ tends to be large
and negative), eventually possibly becoming non-monotonic and/or greater than
unity at high $x$ (the latter e.g. in the case of the down quark, which is not
shown in the figure).

Figures 2 and 3 show results of MC backward evolution from 1 TeV to 10 GeV
using the NEP (57) (black crosses), and the approximation ${\rm NEP}^{({\rm
R})}$ (60) (blue vertical crosses) or ${\rm NEP}^{({\rm E})}$ (61) (red
boxes). Here $10^{7}$ unweighted MC events were generated starting at $\mu=1$
TeV, with a probability distribution of momentum fraction $x$ and flavour $i$

$\frac{dP_{i}}{dx}=xf_{i}(x,\mu^{2})\,,$ (94)

using the momentum sum rule

$\sum_{j}\int_{0}^{1}dx\,xf_{j}(x,\mu^{2})=1$ (95)

as normalization. Following the selection of the next branching scale
$\mu_{0}$ according to the relevant NEP, the momentum fraction $x^{\prime}$
and flavour $j$ of the branching parent was chosen according to the
distribution

$\frac{dP_{j}}{dx^{\prime}}=\frac{1}{x^{\prime}}\Theta_{ij,x/x^{\prime}}^{\scriptscriptstyle\rm
IN}\hat{P}_{ij}(x/x^{\prime})f_{j}(x^{\prime},\mu_{0}^{2})\Bigg{/}\sum_{k}\int_{x}^{1}\frac{dz}{z}\,\Theta_{ik,z}^{\scriptscriptstyle\rm
IN}\hat{P}_{ik}(z)f_{k}(x/z,\mu_{0}^{2})\,.$ (96)

Note that this implies that only resolvable emissions were generated; although
this is the standard practice, it is not necessarily what the various NEPs
employed here would dictate – more details on this point are given in app. A
(see in particular eqs. (151) and (152)). For computational speed, evolution
was discretized on a $(500,70)$-node grid in $(x,\mu)$. The procedure was
iterated until the next branching scale fell below 10 GeV.

In the case of ${\rm NEP}^{({\rm R})}$, we have seen that it may be non-
monotonic or larger than one, in which case the solution of eq. (60) was
chosen larger than or equal to the value of $\mu_{0}$ for which ${\rm
NEP}^{({\rm R})}$ has its minimum. This may account for a part of the large
discrepancies between the results of using ${\rm NEP}^{({\rm R})}$ and the
true NEP or ${\rm NEP}^{({\rm E})}$ at high $x$.

Generally speaking, all versions of the NEP perform poorly in reproducing the
backward evolution of the PDFs, especially outside the intermediate region
$0.01<x<0.1$; the true NEP performs best at higher $x$. However, the
fundamental problem remains that the PDF evolution generated by the MC results
from the accumulation of recoils against resolved emissions, whereas the
actual evolution results from both resolved and unresolved emissions888For an
early prescription to account for unresolved emissions in an average manner,
see ref. Bengtsson:1986gz . .

One possible approach that avoids this problem, while introducing others, is
to guide the backward MC with cutoff-dependent PDFs that are generated by
resolved emissions alone. We consider this approach in detail in the following
section.

Figure 2: Data points show up-quark and gluon PDFs at 10 GeV after MC backward
evolution from 1 TeV, guided by the CT18LO PDFs using the NEP (57) (black
crosses), and the approximation ${\rm NEP}^{({\rm R})}$ (60) (blue vertical
crosses) or ${\rm NEP}^{({\rm E})}$ (61) (red boxes). Solid curves show the
cutoff-independent PDFs at 10 GeV. Also shown (dashed) are the cutoff-
independent PDFs at the starting scale of 1 TeV, to illustrate the amount of
evolution. As in fig. 1, open circles (green) show results obtained with
cutoff-dependent PDFs, to be discussed in sect. 5.

Figure 3: As in fig. 2, but with a logarithmic $x$ scale.

##  5 PDF evolution with cutoff

In sect. 2 we have shown that the correct form of the backward initial-state
radiation NEP is that of eq. (57). However, since the NEP is associated with
the resolution-dependent (non-) emission of resolvable partons, there is an
inconsistency in using it to generate backward shower evolution guided by PDFs
that obey the resolution-independent evolution equations (1). Barring some ad
hoc reweighting procedure, this results in discrepancies between the guiding
PDFs and those generated by backward evolution, which we have illustrated in
the previous section.

An alternative approach is suggested by eq. (62). Namely, one defines a new
type of PDFs, which we denote by $F^{(\epsilon)}$, that obey the following
evolution equation999This possibility was pointed out but not explored in ref.
Marchesini:1987cf .:

$\frac{\partial
F^{(\epsilon)}(x)}{\partial\log\mu^{2}}={\mathbb{O}}^{\scriptscriptstyle\rm
IN}\otimes_{x}F^{(\epsilon)}\,.$ (97)

As the notation suggests, such PDFs depend on the cutoffs
$\epsilon=\\{\epsilon_{ij}^{{\scriptscriptstyle\rm
L}},\epsilon_{ij}^{{\scriptscriptstyle\rm U}}\\}_{ij}$. However, in view of
eqs. (38) and (39), and of the characteristics of ${\mathbb{W}}[F]$ (see in
particular eqs. (40) and (42)), we expect that $F^{(\epsilon)}$ and $F$ will
differ, in the resolved region, only by terms suppressed by some powers of the
cutoffs; conversely, in the unresolved region the differences between the two
are in general logarithmic in the cutoffs.

###  5.1 Flavour and momentum conservation

One interesting question that immediately emerges in the case of cutoff-
dependent evolution, eq. (97), is whether flavour and momentum are conserved.
By working again at the LO, we obtain for the integrated non-singlet
contribution of quark flavour $q$:

$\frac{\partial}{\partial\log\mu^{2}}\int_{0}^{1}dx\Big{(}f_{q}^{(\epsilon)}(x)-f_{\bar{q}}^{(\epsilon)}(x)\Big{)}=\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\,\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{qq}\int_{0}^{1}dy\Big{(}f_{q}^{(\epsilon)}(y)-f_{\bar{q}}^{(\epsilon)}(y)\Big{)}\,,$
(98)

having used the identity:

$\int_{0}^{1}dx\,g\otimes_{x}h=\int_{0}^{1}dz\,g(z)\int_{0}^{1}dy\,h(y)\,.$
(99)

With eqs. (37) and (75) we obtain:

$\displaystyle\int_{0}^{1}dz\,\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{qq}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}dz\,A_{q}(z)\Theta_{qq,z}^{\scriptscriptstyle\rm
IN}+\overline{B}_{q}^{\scriptscriptstyle\rm IN}$ (100) $\displaystyle=$
$\displaystyle 0$ (101) $\displaystyle=$ $\displaystyle
B_{q}^{\scriptscriptstyle\rm IN}\,.$ (102)

The result of eq. (101) stems from eq. (80), whereas that of eq. (102) is what
we would have obtained if we had not chosen $B_{q}^{\scriptscriptstyle\rm
IN}=0$ (see eq. (43)).

This gives us the opportunity to discuss a general property of the Sudakov
definition in the context of MC-compatible PDF evolution equations. In
particular, one observes that this definition is to a certain extent always
arbitrary. In the standard case, such an arbitrariness is associated with the
choice of the resolved region – in practice, with the choices of the cutoffs
and of the functional dependence upon them of the borders of the resolved
region. In the case of cutoff-dependent evolution, in addition to the above
there is the freedom associated with the choice of the parameters
${\mathbb{B}}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}$ (or
$\widetilde{{\mathbb{B}}}^{\mathrel{\vbox{\offinterlineskip\halign{#\cr${\scriptscriptstyle\rm~{}IN}$\cr\kern
0.5pt\cr${\scriptscriptstyle\rm OUT}$\cr}}}}$, the two being related to each
other by relationships such as eq. (51)), given ${\mathbb{B}}$ and the
constraints of eq. (44) (or given $\widetilde{{\mathbb{B}}}$ and the
constraints of eq. (48)). Once these choices have been made, what is
exponentiated ($\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm IN}$) is
determined unambiguously by eq. (43) or eq. (47). In the case of a quark,
$B_{q}=0$, and eq. (44) implies that $B_{q}^{\scriptscriptstyle\rm IN}$ can be
set equal to any function of the cutoffs that vanishes with them. Thus, while
eq. (101) shows that the cutoff-dependent evolution conserves flavour given
the choice of eq. (80), eq. (102) can be seen as constraining
$B_{q}^{\scriptscriptstyle\rm IN}$ by imposing that flavour be conserved. In
other words: the additional freedom of the cutoff-dependent evolution w.r.t.
the standard one discussed above, namely that associated with terms suppressed
by powers of the cutoffs, can be exploited to impose physical conditions (such
as flavour conservation) that, at variance with the case of standard
evolution, may not necessarily emerge in a natural manner.

Turning to the case of momentum conservation101010The issue of momentum
conservation in PDF evolution with a cutoff, in that case due to a modified
argument of $\alpha_{\scriptscriptstyle S}$, was considered in ref.
Skrzypek:2007zz ., we write:

$\displaystyle\frac{\partial}{\partial\log\mu^{2}}\int_{0}^{1}dx\,x\sum_{i}f_{i}^{(\epsilon)}(x)=$
(103) $\displaystyle\phantom{aaaa}\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\,z\Big{(}\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{gg}+\sum_{q,\bar{q}}\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{qg}\Big{)}\int_{0}^{1}dy\,y\,f_{g}^{(\epsilon)}(y)$
$\displaystyle\phantom{aa}\\!+\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\sum_{q}\int_{0}^{1}dz\,z\Big{(}\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{qq}+\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{gq}\Big{)}\int_{0}^{1}dy\,y\Big{(}f_{q}^{(\epsilon)}(y)+f_{\bar{q}}^{(\epsilon)}(y)\Big{)},$

having used the identity:

$\int_{0}^{1}dx\,x\,g\otimes_{x}h=\int_{0}^{1}dz\,z\,g(z)\int_{0}^{1}dy\,y\,h(y)\,.$
(104)

Let us start by considering the integral over $x$ in the third line of eq.
(103). By proceeding as was done for the manipulations of the flavour-
conservation case, we obtain:

$\displaystyle\int_{0}^{1}dz\,z\Big{(}\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{qq}+\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{gq}\Big{)}=$ (105)
$\displaystyle\phantom{aaaa}\int_{0}^{1}dz\,z\left[A_{q}(z)\Theta_{qq,z}^{\scriptscriptstyle\rm
IN}+C_{gq}(z)\Theta_{gq,z}^{\scriptscriptstyle\rm
IN}\right]+B_{q}^{\scriptscriptstyle\rm
IN}-\int_{0}^{1}dz\,A_{q}(z)\Theta_{qq,z}^{\scriptscriptstyle\rm
IN}=B_{q}^{\scriptscriptstyle\rm IN}\,,$

where the rightmost side follows from the direct computation of the integrals
that appear in the central expression. Since $B_{q}^{\scriptscriptstyle\rm
IN}=0$ as we have discussed above, the integral on the l.h.s. of eq. (105) is
thus equal to zero. Turning to the integral over $z$ in the second line of eq.
(103), we have:

$\displaystyle\int_{0}^{1}dz\,z\Big{(}\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{gg}+\sum_{q,\bar{q}}\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)\big{)}_{qg}(z)\Big{)}=$ (106)
$\displaystyle\phantom{aa}\int_{0}^{1}dz\,z\left[\left(\frac{\widetilde{A}_{g}(z)}{1-z}+C_{gg}(z)\right)\Theta_{gg,z}^{\scriptscriptstyle\rm
IN}+\sum_{q,\bar{q}}C_{qg}(z)\Theta_{qg,z}^{\scriptscriptstyle\rm
IN}\right]+\widetilde{B}_{g}^{\scriptscriptstyle\rm
IN}-\int_{0}^{1}dz\,\frac{\widetilde{A}_{g}(1)}{1-z}\Theta_{gg,z}^{\scriptscriptstyle\rm
IN}\,.$

By using eqs. (78), (82), and (83) (the latter two in eq. (50)), one sees that
the integral on the l.h.s. of eq. (106) is equal to zero. Combined with the
null result in eq. (105), this is the analogue of eq. (101), and concludes the
proof that the cutoff-dependent evolution conserves the momentum. Conversely,
we can proceed by analogy with eq. (102), and impose the l.h.s. of eq. (106)
to be equal to zero in order to determine
$\widetilde{B}_{g}^{\scriptscriptstyle\rm IN}$. If we then adopt the local
form (50):

$\widetilde{B}^{{\scriptscriptstyle\rm
IN}}_{g}=\int_{0}^{1}dz\left[\tilde{b}_{gg}(z)\Theta_{gg,z}^{\scriptscriptstyle\rm
IN}+\sum_{q,\bar{q}}\tilde{b}_{qg}(z)\Theta_{qg,z}^{\scriptscriptstyle\rm
IN}\right],$ (107)

we obtain

$\displaystyle\frac{2\pi}{\alpha_{\scriptscriptstyle S}}\,\tilde{b}_{gg}(z)$
$\displaystyle=$ $\displaystyle 2C_{\scriptscriptstyle
A}\,z\Big{(}2-z+z^{2}\Big{)}\,,$ (108)
$\displaystyle\frac{2\pi}{\alpha_{\scriptscriptstyle S}}\,\tilde{b}_{qg}(z)$
$\displaystyle=$ $\displaystyle-z\hat{P}_{qg}(z)\,,$ (109)

leading via eq. (35) to the following expression for the integrand in the
exponent of the gluon Sudakov form factor:

$\overline{B}_{g}^{\scriptscriptstyle\rm IN}=-\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\,z\left(\hat{P}_{gg}(z)\Theta_{gg,z}^{\scriptscriptstyle\rm
IN}+\sum_{q,\bar{q}}\hat{P}_{qg}(z)\Theta_{qg,z}^{\scriptscriptstyle\rm
IN}\right)\,.$ (110)

Although different from eq. (85), this expression is equally valid, since the
two integrands differ by a function that integrates to zero. In an analogous
manner, from eq. (105) we would obtain for the quark

$\overline{B}_{q}^{\scriptscriptstyle\rm IN}=-\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\int_{0}^{1}dz\,z\left(\hat{P}_{qq}(z)\Theta_{qq,z}^{\scriptscriptstyle\rm
IN}+\hat{P}_{gq}(z)\Theta_{gq,z}^{\scriptscriptstyle\rm IN}\right)\,,$ (111)

which again coincides with eq. (81), in spite of having a different integrand.
We point out that the strict equality of the results for
$\overline{B}_{g}^{\scriptscriptstyle\rm IN}$ stemming from eqs. (110) and
(85), and of those for $\overline{B}_{q}^{\scriptscriptstyle\rm IN}$ from eqs.
(111) and (81), relies among other things on the symmetry properties of the
$\Theta_{ij,z}^{\scriptscriptstyle\rm IN}$ functions. On the other hand, the
integrands of eqs. (110) and (111), at variance with those of eqs. (85) and
(81), do not have a $z\to 0$ singularity when $\epsilon\to 0$. This implies
that they lead to finite quantities also when completely removing the
constraints enforced by the lower cutoffs; such quantities can then be
employed to define Sudakov factors that differ from those used thus far by
terms suppressed by powers of the cutoff111111For examples of Sudakov form
factors in a different context, whose definitions do differ from one another
by cutoff-suppressed terms, see e.g. app. A of ref. Frederix:2020trv ..

###  5.2 Results on cutoff-dependent PDFs

Figure 4: Cutoff-dependent PDFs evolved from the CT18LO set at 100 GeV
backwards to 10 GeV (red, dashed) and forwards to 1 TeV (blue, dot-dashed),
compared to cutoff-independent PDFs (solid) at 10 GeV (red), 100 GeV (black)
and 1 TeV (blue). The ratio plots show the cutoff-dependent PDFs relative to
the cutoff-independent ones at the same scale.

Figure 4 shows examples of cutoff-dependent PDFs corresponding to the two
cutoff choices discussed in sect. 4.1. Starting from the CT18LO set at scale
$\mu=100$ GeV, the PDFs were evolved forwards to 1 TeV and backwards to 10 GeV
using eq. (97) in place of (1), with the flavour- and momentum-conserving
formulation described above.

As expected, the cutoff-dependent PDFs generally evolve more slowly with
increasing scale than the true PDFs, thus being generally softer below the
starting scale and harder above it. The relative differences grow with
increasing $x$ and are largest in the unresolved region $x>1-\epsilon$, where
the PDFs are however very small. The scale-dependent cutoff naturally leads to
PDFs much closer to the cutoff-independent ones at high scales.

Non-emission probabilities for backward evolution guided by the cutoff-
dependent PDFs, chosen to coincide with the CT18LO set at 10 GeV, are shown by
green circles in fig. 1 as a function of $\mu_{0}$ with $\mu=100$ GeV. Since
by construction ${\mathbb{W}}\left[F^{(\epsilon)}\right]\equiv 0$, all of the
expressions for the NEP are now equivalent, by virtue of eq. (62).
Empirically, the NEP for cutoff-dependent PDFs appears closest to ${\rm
NEP}^{({\rm E})}$ computed from cutoff-independent PDFs.

Results of backward MC evolution guided by the cutoff-dependent PDFs are also
shown by green circles in figs. 2 and 3. There, in contrast to fig. 4 but
analogously to fig. 1, the cutoff-dependent PDFs were chosen to coincide with
the CT18LO set at 10 GeV and evolved upwards to 1 TeV, where they were used as
the starting distributions for the backward MC. In this way, the backward-
generated MC distributions at 10 GeV should agree with the CT18LO set.
Compared to the results using cutoff-independent PDFs, agreement is indeed
greatly improved at all $x$ values. The small residual systematic
discrepancies are most likely due to accumulated errors from our
discretization of the backward evolution.

##  6 Cutoff-dependent cross sections

The cutoff-dependent PDFs emerging from eq. (97) imply that short-distance
cross sections must be cutoff-dependent too, in order for the l.h.s. of the
factorisation formula to be cutoff independent121212Up to terms one
perturbative order higher than those included in the computation of the short
distance cross sections.. We shall assume in what follows that the cutoff
dependence of the PDFs is solely due to their evolution. This implies that, at
the scale chosen as the starting point for PDF evolution, the initial
conditions must be cutoff independent; this is not mandatory, but doing
otherwise would require some modeling assumptions for the initial conditions.
In order to determine the cutoff-dependent terms of the cross section, we
consider the generic factorisation formula for one incoming leg, starting from
the cutoff-independent case:

$\sigma=F^{\rm
T}\star\hat{\Sigma}\equiv\sum_{i}\int_{0}^{1}dx\,f_{i}(x)\,\hat{\sigma}_{i}(x)\,.$
(112)

Here, we have denoted by $\hat{\Sigma}$ the column vector that collects all of
the (subtracted) short-distance cross sections
$\hat{\sigma}_{i}\equiv(\hat{\Sigma})_{i}$, whereas $\sigma$ is the hadron-
level cross section that results from the sum over all of the partonic
processes in eq. (112). The RGE invariance of $\sigma$ under factorisation-
scale variation is:

$0=\frac{\partial\sigma}{\partial\log\mu^{2}}=\frac{\partial F^{\rm
T}}{\partial\log\mu^{2}}\star\hat{\Sigma}+F^{\rm
T}\star\frac{\partial\hat{\Sigma}}{\partial\log\mu^{2}}=\left({\mathbb{O}}\otimes
F\right)^{\rm T}\star\hat{\Sigma}+F^{\rm
T}\star\frac{\partial\hat{\Sigma}}{\partial\log\mu^{2}}\,.$ (113)

It is a matter of algebra to show that, for any functions $g$, $h$, and $l$,
the following identity holds:

$\big{(}g\otimes h\big{)}\star l=g\star\big{(}h\star l)=h\star\big{(}g\star
l)\,.$ (114)

Equation (113) then implies:

$0=F^{\rm T}\star\frac{\partial\hat{\Sigma}}{\partial\log\mu^{2}}+F^{\rm
T}\star\left({\mathbb{O}}^{\rm T}\star\hat{\Sigma}\right).$ (115)

This equation must be true for any PDFs, and therefore:

$\frac{\partial\hat{\Sigma}}{\partial\log\mu^{2}}=-{\mathbb{O}}^{\rm
T}\star\hat{\Sigma}\;\;\;\;\Longleftrightarrow\;\;\;\;\frac{\partial\hat{\sigma}_{i}(x)}{\partial\log\mu^{2}}=-\sum_{j}\int_{0}^{1}dy\,\big{(}{\mathbb{O}}(y)\big{)}_{ji}\,\hat{\sigma}_{j}(xy)\,.$
(116)

By writing the perturbative expansion of the short-distance cross sections as
follows:

$\hat{\Sigma}=\hat{\Sigma}^{[0]}+\frac{\alpha_{\scriptscriptstyle
S}}{2\pi}\,\hat{\Sigma}^{[1]}+\ldots$ (117)

where all terms $\hat{\Sigma}^{[i]}$ include a factor
$\alpha_{\scriptscriptstyle S}^{b}$, with $b$ a process-dependent constant
(e.g. $b=0$ and $b=2$ for dilepton and Higgs production, respectively), and by
working with LO kernels where eq. (74) holds, eq. (116) implies:

$\hat{\Sigma}^{[1]}=-\log\frac{\mu^{2}}{q_{0}^{2}}\,\left(\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}{\mathbb{O}}^{{\rm
T}}\star\hat{\Sigma}^{[0]}\right)+C\equiv-\log\frac{\mu^{2}}{q_{0}^{2}}\,\left({\mathbb{P}}^{[0]{\rm
T}}\star\hat{\Sigma}^{[0]}\right)+C\,,$ (118)

with $q_{0}$ an arbitrary reference scale, and $C$ a column vector of
$\mu$-independent integration constants. The determination of $C$ can be done
by means of an explicit cross section calculation. For example, it can be read
from the FKS formalism Frixione:1995ms ; Frixione:1997np , where a term in the
same form as the leftmost one on the r.h.s. of eq. (118) is contained in the
so-called $(n+1)$-body degenerate contributions.

The derivation above can be repeated verbatim for the cutoff-dependent PDFs
and short-distance cross sections. Denoting the latter by
$\hat{\Sigma}^{(\epsilon)}$, owing to eq. (97) the analogue of eq. (118) reads
as follows:

$\hat{\Sigma}^{(\epsilon)[1]}=-\log\frac{\mu^{2}}{q_{0}^{2}}\,\left(\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm IN}\big{)}^{\rm
T}\star\hat{\Sigma}^{(\epsilon)[0]}\right)+C^{(\epsilon)}\,.$ (119)

Under our assumptions concerning the cutoff dependence discussed at the
beginning of this section, we may now set:

$\hat{\Sigma}^{(\epsilon)[0]}=\hat{\Sigma}^{[0]}\,.$ (120)

Furthermore, by choosing $q_{0}$ to coincide with the starting scale of the
PDF evolution, at $\mu=q_{0}$ we must have:

$\hat{\Sigma}^{(\epsilon)[1]}=\hat{\Sigma}^{[1]}\;\;\;\;\Longrightarrow\;\;\;\;C^{(\epsilon)}=C\,,$
(121)

and therefore, for a generic scale value:

$\hat{\Sigma}^{(\epsilon)[1]}=\hat{\Sigma}^{[1]}+\log\frac{\mu^{2}}{q_{0}^{2}}\,\left[\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\left({\mathbb{O}}-{\mathbb{O}}^{\scriptscriptstyle\rm IN}\right)^{{\rm
T}}\star\hat{\Sigma}^{[0]}\right]\equiv\hat{\Sigma}^{[1]}+\log\frac{\mu^{2}}{q_{0}^{2}}\,\left(\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\big{(}{\mathbb{O}}^{\scriptscriptstyle\rm OUT}\big{)}^{\rm
T}\star\hat{\Sigma}^{[0]}\right)\,.$ (122)

Equation (122) allows one to obtain the sought cutoff-dependent short-distance
cross sections given the cutoff-independent ones. The rightmost term on the
r.h.s. of eq. (122) is, as expected, suppressed by powers of the cutoff; we
shall call it the cutoff correction.

We note that by iteration of this procedure one can obtain the cutoff
correction to any perturbative order, in terms of contributions of lower
orders to the short-distance cross section and the cutoff-dependent and
cutoff-independent evolution kernels.

###  6.1 Results for cross sections

####  6.1.1 Drell-Yan process

As a first illustration of the use of cutoff-dependent PDFs with a cutoff-
corrected short-distance cross section, we consider the photon-induced ${\cal
O}(\alpha^{2})$ and ${\cal O}(\alpha^{2}\alpha_{\scriptscriptstyle S})$
contributions to the cross section for lepton pair production as a function of
pair invariant mass $M_{ll}$ at fixed hadronic collision energy $\sqrt{s}$.

Some results for $pp$ collisions at $\sqrt{s}=13$ TeV are shown in fig. 5. As
in sect. 4.1, we use the CT18LO leading-order PDFs, and we do so with both LO
(${\cal O}(\alpha^{2})$) and NLO (${\cal
O}(\alpha^{2}+\alpha^{2}\alpha_{\scriptscriptstyle S})$) short-distance cross
sections, for the latter of which we employ the $\overline{\rm MS}$
factorisation scheme. We again consider two cases of a universal, flavour-
independent cutoff $\epsilon_{ij}^{{\scriptscriptstyle\rm
L}}=\epsilon_{ij}^{{\scriptscriptstyle\rm U}}=\epsilon$: one relatively large
and scale-independent, $\epsilon=0.1$, the other scale-dependent,
$\epsilon=\;\mbox{(2 GeV)}/q$, with $q$ the relevant mass scale. The reference
scale $q_{0}$, at which the cutoff-dependent and cutoff-independent PDFs are
identical, is set equal to 10 GeV in the upper plots and to 100 GeV in the
lower ones. The scale for evaluation of the PDFs, $\alpha_{\scriptscriptstyle
S}$, and NLO corrections is taken to be $\mu^{2}=M_{ll}^{2}$ throughout.

At leading order there is no cutoff correction to the short-distance cross
section and so the discrepancies between the cutoff-dependent (blue, dashed)
and cutoff-independent (black, solid) LO results simply reflect those between
the corresponding PDFs. Since the cutoff-dependent PDFs evolve more slowly,
the resulting hadronic cross section initially falls below the true (i.e.
obtained with cutoff-independent PDFs) LO value for $M_{ll}>q_{0}$ but
eventually rises above it at higher values of $M_{ll}$ (higher $x$).
Correspondingly, for $q_{0}=100$ GeV, it lies above the true LO value when
$M_{ll}<q_{0}$.

At next-to-leading order, the cutoff correction comes into play and reduces
the discrepancy between the cutoff-dependent and true NLO results. The
reduction is strong around the reference scale $M_{ll}\sim q_{0}$, but only
modest above and very rapidly deteriorating below $q_{0}$. For the scale-
dependent cutoff with a low reference scale (the upper right plot), the effect
of the cutoff correction vanishes much more rapidly than that of the
difference in PDFs at high $M_{ll}$. This is because the relevant scale in the
cutoff correction is the local value $q=M_{ll}$, whereas the difference in
PDFs results from the accumulation of cutoff effects over the whole range from
$q_{0}$ to $M_{ll}$.

In summary, the comparison of the LO and NLO results shows that the NLO cutoff
correction of eq. (122) partly compensates for the differences between the
cutoff-independent predictions and those one would have obtained by employing
cutoff-dependent PDFs without the inclusion of such a correction in the short-
distance cross sections. In general, the use of cutoff-dependent PDFs together
with the correction (122) gives results for the Drell-Yan cross section that
are relatively close to the cutoff-independent NLO predictions, provided the
reference scale $q_{0}$ is close to, or not too far below, the dilepton mass.
We point out that, at this level of accuracy, a more systematic assessment of
the compensation mechanism just mentioned would require the definition of a
proper NLO cutoff-dependent PDF set.

Figure 5: Drell-Yan cross section at $\sqrt{s}=13$ TeV (photon-induced
contribution only), calculated using cutoff-dependent PDFs at leading order
(blue, dashed) and next-to-leading order (red, dot-dashed), compared to
corresponding results using cutoff-independent PDFs (solid).

####  6.1.2 Higgs boson production

Since the Drell-Yan process is quark dominated, we consider as a second
example the gluon fusion contribution to Higgs boson hadroproduction as a
function of the hadronic collision energy $\sqrt{s}$.

Figure 6: Higgs cross section at $\sqrt{s}=1-100$ TeV (gluon fusion
contribution only), calculated using cutoff-dependent PDFs at leading order
(blue, dashed) and next-to-leading order (red, dot-dashed), compared to the
corresponding results using cutoff-independent PDFs (solid).

Figure 6 shows results (in the infinite top mass approximation) for $pp$
collisions at $\sqrt{s}=1-100$ TeV. The PDFs, $\alpha_{\scriptscriptstyle S}$,
cutoffs and factorisation scheme are as in sect. 6.1.1, but now the scale used
in their evaluation is fixed at $\mu^{2}=m^{2}_{h}$. Thus the differences
between the cutoff-dependent (blue, dashed) and cutoff-independent (black,
solid) LO results simply reflect the different $x$ dependences of the
corresponding gluon PDFs. Since by construction the cutoff-dependent and
-independent PDFs coincide at the reference scale $q_{0}$, and the former
evolve more slowly, the corresponding LO result falls below the true (i.e.
obtained with cutoff-independent PDFs) value at high $\sqrt{s}$ (small $x$) as
long as $q_{0}<m_{h}$. This discrepancy is naturally much smaller for
$q_{0}=100\;\mbox{GeV}\;\sim m_{h}$ than for $q_{0}=10$ GeV.

The next-to-leading order corrections to Higgs production are very large,
comparable to the leading order. The relative differences between the cutoff-
dependent (red, dot-dashed) and true cutoff-independent (black, solid) NLO
results are reduced compared to leading order. Again they are naturally much
smaller for $q_{0}=100$ GeV than for $q_{0}=10$ GeV.

##  7 Conclusions

Our aim in this paper has been to study the extent to which the use of PDFs to
guide backward MC parton showering can be a consistent procedure. We have
shown that it cannot be fully so if normal cutoff-independent PDFs are used,
even if the non-emission probabilities (NEPs) currently in use in Monte Carlo
event generators (MCEGs) are corrected to account for the fact that they
generate only resolved parton emissions. The cutoffs inherent in the
resolution criteria lead to inconsistencies that are formally power-suppressed
in the resolved region. Nevertheless, these can accumulate to have large
effects when showers evolve over a wide range of scales, and increase with
$x$.

As an alternative, formally more consistent approach, we have considered the
use of cutoff-dependent PDFs, together with short-distance cross sections that
include compensating cutoff corrections. We have illustrated the extent to
which this compensation works at NLO in lepton pair and Higgs boson
production.

Obviously, if the use of cutoff-dependent PDFs for event generation at hadron
colliders is to be pursued, global PDF fits tailored to the sets of cutoffs in
the widely used MCEGs would need to be performed. In principle this seems a
straightforward matter of using the cutoff-dependent PDF evolution kernels and
corresponding subprocess cutoff corrections; at leading order, this could be a
worthwhile improvement on the current practice. Beyond leading order, however,
the whole concept of guided backward parton showering needs further
clarification.

## Acknowledgements

We are grateful to Torbjörn Sjöstrand and Mike Seymour for valuable comments
on the manuscript. This work has been partially supported by UK STFC HEP
Theory Consolidated grant ST/T000694/1. SF thanks the CERN TH division for the
kind hospitality during the course of this work.

##  Appendix A PDF reconstruction with MC backward evolution

In order to see how MC initial-state showers reconstruct PDFs, we first need
to find a solution of the evolution equations that renders the comparison with
MC-derived results as easy as possible. To this end, we employ eqs. (36) and
(37), and rewrite the latter as follows:

${\mathbb{O}}^{\scriptscriptstyle\rm
IN}(z)={\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN}(z)+\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm IN}\delta(1-z)\,,$ (123)

so that:

${\mathbb{Z}}\left[F\right](x)={\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN}\otimes_{x}F\,.$ (124)

In other words, ${\mathbb{O}}_{R}^{\scriptscriptstyle\rm IN}$ is the
contribution to the inner-region evolution operator
${\mathbb{O}}^{\scriptscriptstyle\rm IN}$ due to real (as opposed to virtual)
emissions. With eqs. (36) and (124) one obtains:

$\displaystyle M\Big{[}{\mathbb{W}}\left[F\right]\Big{]}$ $\displaystyle=$
$\displaystyle M\Big{[}{\mathbb{O}}^{\scriptscriptstyle\rm
OUT}\Big{]}M\Big{[}F\Big{]}\equiv{\mathbb{O}}_{N}^{\scriptscriptstyle\rm
OUT}\,F_{N}\,,$ (125) $\displaystyle
M\Big{[}{\mathbb{Z}}\left[F\right]\Big{]}$ $\displaystyle=$ $\displaystyle
M\Big{[}{\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN}\Big{]}M\Big{[}F\Big{]}\equiv{\mathbb{O}}_{R,N}^{\scriptscriptstyle\rm
IN}\,F_{N}\,,$ (126)

where by $M[g]\equiv g_{N}$ we have denoted the Mellin transform of a function
$g(x)$. Thus, the Mellin transform of eq. (31) reads as follows:

$F_{N}(\mu^{2})=\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,F_{N}(\mu_{0}^{2})+\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\kappa^{2})}\,\Big{(}{\mathbb{O}}_{N}^{\scriptscriptstyle\rm
OUT}(\kappa^{2})+{\mathbb{O}}_{R,N}^{\scriptscriptstyle\rm
IN}(\kappa^{2})\Big{)}\,F_{N}(\kappa^{2})\,.$ (127)

Equation (127) is a Volterra equation of the second kind131313Its kernel is
separable in the two relevant variables ($\mu^{2}$ and $\kappa^{2}$), which
leads to (at least in a one-dimensional flavour space) a closed-form solution;
this, however, is not of particular interest here, and will not be
considered., which is formally solved by a Neumann series:

$F_{N}(\mu^{2})=\sum_{k=0}^{\infty}F_{N}^{(k)}(\mu^{2})\,,$ (128)

with:

$\displaystyle F_{N}^{(0)}(\mu^{2})$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,F_{N}(\mu_{0}^{2})\,,$
(129) $\displaystyle F_{N}^{(k)}(\mu^{2})$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\int_{\mu_{0}^{2}}^{\mu^{2}}\Bigg{[}\prod_{p=1}^{k}\frac{d\kappa_{p}^{2}}{\kappa_{p}^{2}}\,\Theta\left(\kappa_{p+1}^{2}\leq\kappa_{p}^{2}\leq\kappa_{p-1}^{2}\right)\frac{{\mathbb{S}}(\kappa_{p-1}^{2})}{{\mathbb{S}}(\kappa_{p}^{2})}$
(130)
$\displaystyle\phantom{\int_{\mu_{0}^{2}}^{\infty}\Bigg{[}a}\times\Big{(}{\mathbb{O}}_{N}^{\scriptscriptstyle\rm
OUT}(\kappa_{p}^{2})+{\mathbb{O}}_{R,N}^{\scriptscriptstyle\rm
IN}(\kappa_{p}^{2})\Big{)}\Bigg{]}\frac{{\mathbb{S}}(\kappa_{k}^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,F_{N}(\mu_{0}^{2})\,,\phantom{aa}$

where the matrix product in eq. (130) has a left-to-right order, i.e. the
elements corresponding to $p=1$ ($p=k$) are the leftmost (rightmost) ones, and
we have defined:

$\kappa_{0}^{2}=\mu^{2}\,,\;\;\;\;\;\;\;\;\kappa_{k+1}^{2}=\mu_{0}^{2}\,.$
(131)

Equations (129) and (130) can be easily transformed back to the $x$ space:

$\displaystyle F^{(0)}(x,\mu^{2})$ $\displaystyle\\!=\\!$
$\displaystyle\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,F(x,\mu_{0}^{2})\,,$
(132) $\displaystyle F^{(k)}(x,\mu^{2})$ $\displaystyle\\!=\\!$
$\displaystyle\int_{0}^{1}\left[\prod_{q=1}^{k+1}dz_{q}\right]\delta\\!\left(x-\prod_{q=1}^{k+1}z_{q}\right)$
$\displaystyle\\!\\!\\!\\!\times\int_{\mu_{0}^{2}}^{\mu^{2}}\Bigg{[}\prod_{p=1}^{k}\frac{d\kappa_{p}^{2}}{\kappa_{p}^{2}}\,\Theta\left(\kappa_{p+1}^{2}\leq\kappa_{p}^{2}\leq\kappa_{p-1}^{2}\right)\frac{{\mathbb{S}}(\kappa_{p-1}^{2})}{{\mathbb{S}}(\kappa_{p}^{2})}$
$\displaystyle\phantom{\int_{\mu_{0}^{2}}^{\infty}\Bigg{[}a}\times\Big{(}{\mathbb{O}}^{\scriptscriptstyle\rm
OUT}(z_{p},\kappa_{p}^{2})+{\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN}(z_{p},\kappa_{p}^{2})\Big{)}\Bigg{]}\frac{{\mathbb{S}}(\kappa_{k}^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,F(z_{k+1},\mu_{0}^{2})\,.\phantom{aa}$

The presence of a Dirac $\delta$ in eq. (A) complicates its manipulation. We
therefore introduce the $k$ independent variables:

$y_{i}=\prod_{j=i+1}^{k+1}z_{j}\,,\;\;\;\;\;\;\;\;1\leq i\leq k\,,$ (134)

and the dummy variable

$y_{0}\equiv x\,,$ (135)

so that

$y_{i}=z_{i+1}y_{i+1}\;\;\;({\rm for}~{}0\leq i\leq
k-1)\,,\;\;\;\;\;\;\;\;y_{k}=z_{k+1}\,.$ (136)

In this way, by using the identity:

$1=\int_{0}^{1}\left(\prod_{i=1}^{k-1}dy_{i}\,\delta\big{(}y_{i}-z_{i+1}y_{i+1}\big{)}\right)dy_{k}\,\delta\big{(}y_{k}-z_{k+1}\big{)}\,,$
(137)

eq. (A) becomes:

$\displaystyle F^{(k)}(z,\mu^{2})$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\int_{0}^{1}\left[\prod_{q=1}^{k}\frac{dy_{q}}{y_{q}}\right]\int_{\mu_{0}^{2}}^{\mu^{2}}\Bigg{[}\prod_{p=1}^{k}\frac{d\kappa_{p}^{2}}{\kappa_{p}^{2}}\,\Theta\left(\kappa_{p+1}^{2}\leq\kappa_{p}^{2}\leq\kappa_{p-1}^{2}\right)\frac{{\mathbb{S}}(\kappa_{p-1}^{2})}{{\mathbb{S}}(\kappa_{p}^{2})}$
$\displaystyle\phantom{\int_{\mu_{0}^{2}}^{\infty}\Bigg{[}a}\times\left({\mathbb{O}}^{\scriptscriptstyle\rm
OUT}\left(\frac{y_{p-1}}{y_{p}},\kappa_{p}^{2}\right)+{\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN}\left(\frac{y_{p-1}}{y_{p}},\kappa_{p}^{2}\right)\right)\Bigg{]}\frac{{\mathbb{S}}(\kappa_{k}^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,F(y_{k},\mu_{0}^{2})\,.\phantom{aa}$

Note that from eq. (136):

$x\equiv y_{0}\leq y_{1}\leq\ldots y_{k-1}\leq y_{k}\,,$ (139)

which are automatically enforced by the requirement that the first arguments
of ${\mathbb{O}}^{\scriptscriptstyle\rm OUT}$ and
${\mathbb{O}}_{R}^{\scriptscriptstyle\rm IN}$ in eq. (A) be less than one. For
a given parton identity $i_{0}$, eq. (A) gives:

$\displaystyle f_{i_{0}}^{(1)}(x,\mu^{2})$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\sum_{i_{1}}\int_{0}^{1}\frac{dy_{1}}{y_{1}}\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa_{1}^{2}}{\kappa_{1}^{2}}\,\frac{S_{i_{0}}(\mu^{2})}{S_{i_{0}}(\kappa_{1}^{2})}$
$\displaystyle\phantom{aaaaa}\times\left(({\mathbb{O}}^{\scriptscriptstyle\rm
OUT})_{i_{0}i_{1}}\\!\left(\frac{x}{y_{1}},\kappa_{1}^{2}\right)+({\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN})_{i_{0}i_{1}}\\!\left(\frac{x}{y_{1}},\kappa_{1}^{2}\right)\right)\frac{S_{i_{1}}(\kappa_{1}^{2})}{S_{i_{1}}(\mu_{0}^{2})}\,f_{i_{1}}(y_{1},\mu_{0}^{2})\,,\phantom{aa}$
$\displaystyle f_{i_{0}}^{(2)}(x,\mu^{2})$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\sum_{i_{1}i_{2}}\int_{0}^{1}\frac{dy_{1}}{y_{1}}\frac{dy_{2}}{y_{2}}\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa_{1}^{2}}{\kappa_{1}^{2}}\frac{d\kappa_{2}^{2}}{\kappa_{2}^{2}}\,\Theta\left(\kappa_{2}^{2}\leq\kappa_{1}^{2}\right)\frac{S_{i_{0}}(\mu^{2})}{S_{i_{0}}(\kappa_{1}^{2})}$
$\displaystyle\phantom{aaaaa}\times\left(({\mathbb{O}}^{\scriptscriptstyle\rm
OUT})_{i_{0}i_{1}}\\!\left(\frac{x}{y_{1}},\kappa_{1}^{2}\right)+({\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN})_{i_{0}i_{1}}\\!\left(\frac{x}{y_{1}},\kappa_{1}^{2}\right)\right)\frac{S_{i_{1}}(\kappa_{1}^{2})}{S_{i_{1}}(\kappa_{2}^{2})}$
$\displaystyle\phantom{aaaaa}\times\left(({\mathbb{O}}^{\scriptscriptstyle\rm
OUT})_{i_{1}i_{2}}\\!\left(\frac{y_{1}}{y_{2}},\kappa_{2}^{2}\right)+({\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN})_{i_{1}i_{2}}\\!\left(\frac{y_{1}}{y_{2}},\kappa_{2}^{2}\right)\right)\frac{S_{i_{2}}(\kappa_{2}^{2})}{S_{i_{2}}(\mu_{0}^{2})}\,f_{i_{2}}(y_{2},\mu_{0}^{2})\,,\phantom{aa}$

and so forth.

We now assume that the parton type $i_{0}$, momentum fraction $x$, and two
scales $\mu^{2}$ and $\mu_{0}^{2}$ (with $\mu_{0}^{2}\leq\mu^{2}$) are given,
and we want to compute the probability that, starting the evolution at
$(x,\mu^{2})$, one eventually (i.e. after an arbitrary number of backward
emissions, including none) ends up emitting at a scale lower than
$\mu_{0}^{2}$. Such a probability is the sum of the probabilities $p_{k}$
associated with $k$ emissions, with $0\leq k\leq\infty$. As far as $k=0$ is
concerned, $p_{0}$ is equal to one minus the probability of emitting at scales
larger than $\mu_{0}^{2}$, in turn equal to the non-emission probability in
$(\mu_{0}^{2},\mu^{2})$ at $x$. Thus, one needs to start with a definite
choice for the latter; we begin by considering ${\rm NEP}^{({\rm R})}$ of eq.
(60). Hence:

$\displaystyle p_{0}$ $\displaystyle=$ $\displaystyle{\rm NEP}_{i_{0}}^{({\rm
R})}(x,\mu_{0}^{2},\mu^{2})=\frac{S_{i_{0}}(\mu^{2})}{S_{i_{0}}(\mu_{0}^{2})}\,\frac{f_{i_{0}}(x,\mu_{0}^{2})}{f_{i_{0}}(x,\mu^{2})}\,.$
(142) $\displaystyle=$
$\displaystyle\frac{f_{i_{0}}^{(0)}(x,\mu^{2})}{f_{i_{0}}(x,\mu^{2})}\,,$
(143)

having used eq. (132) in the second line. The case $k=1$ corresponds to one
emission at a scale $\kappa_{1}^{2}\in(\mu_{0}^{2},\mu^{2})$, followed by one
emission below $\mu_{0}^{2}$. As far as the probability associated with the
former is concerned, one first determines the scale at which it occurs by
solving

$r={\rm NEP}_{i_{0}}^{({\rm R})}(x,\kappa_{1}^{2},\mu^{2})$ (144)

for $\kappa_{1}^{2}$, given a uniform random number $r$. The solution is
discarded if $\kappa_{1}^{2}<\mu_{0}^{2}$, which gives the correct
normalisation. Indeed, the distribution in $\log\kappa_{1}^{2}$ induced by eq.
(144) is:

$\frac{\partial}{\partial\log\kappa_{1}^{2}}\,{\rm NEP}_{i_{0}}^{({\rm
R})}(x,\kappa_{1}^{2},\mu^{2})\,,\;\;\;\;\;\;\;\;\mu_{0}^{2}\leq\kappa_{1}^{2}\leq\mu^{2}\,,$
(145)

so that the total probability for such an emission is:

$\displaystyle\int_{\mu_{0}^{2}}^{\mu^{2}}d\log\kappa_{1}^{2}\frac{\partial}{\partial\log\kappa_{1}^{2}}\,{\rm
NEP}_{i_{0}}^{({\rm R})}(x,\kappa_{1}^{2},\mu^{2})$ $\displaystyle=$
$\displaystyle{\rm NEP}_{i_{0}}^{({\rm R})}(x,\mu^{2},\mu^{2})-{\rm
NEP}_{i_{0}}^{({\rm R})}(x,\mu_{0}^{2},\mu^{2})$ (146) $\displaystyle=$
$\displaystyle 1-{\rm NEP}_{i_{0}}^{({\rm R})}(x,\mu_{0}^{2},\mu^{2})\,.$

As was already reported in eq. (64), by means of an explicit computation we
obtain:

$\frac{\partial}{\partial\log\kappa_{1}^{2}}\,{\rm NEP}_{i_{0}}^{({\rm
R})}(x,\kappa_{1}^{2},\mu^{2})=\frac{1}{f_{i_{0}}(x,\mu^{2})}\frac{S_{i_{0}}(\mu^{2})}{S_{i_{0}}(\kappa_{1}^{2})}\Big{[}\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})+\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})\Big{]}.$
(147)

Having computed the probability density for an emission at
$\kappa_{1}^{2}>\mu_{0}^{2}$, we need to multiply it by the probability of the
next emission occurring below $\mu_{0}^{2}$. This is functionally the same
quantity as was computed for the $k=0$ case in eq. (142); presently, we must
use that result by replacing $\mu^{2}\to\kappa_{1}^{2}$, and $z$ and $i_{0}$
with the momentum fraction $y_{1}$ and parton type $i_{1}$ that have resulted
from the branching at $\kappa_{1}^{2}$, which may have changed w.r.t. the
original $z$ and $i_{0}$. In order to determine these, and the probabilities
associated with their choices, we introduce the functions:

${\cal
P}_{ij}(y;x,\kappa^{2})=\int_{y}^{1}\frac{d\omega}{\omega}\left[({\mathbb{O}}^{\scriptscriptstyle\rm
OUT})_{ij}\left(\frac{x}{\omega},\kappa^{2}\right)+({\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN})_{ij}\left(\frac{x}{\omega},\kappa^{2}\right)\right]f_{j}(\omega,\kappa^{2})\,,$
(148)

with $y\geq x$; these are such that:

$\sum_{j}{\cal
P}_{ij}(x;x,\kappa^{2})=\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i}(x,\kappa^{2})+\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i}(x,\kappa^{2})\,.$
(149)

We first define $i_{1}$ to be the smallest index that fulfills the following
inequality:

$r\leq\frac{\sum_{j}^{j\leq i_{1}}{\cal
P}_{i_{0}j}(x;x,\kappa_{1}^{2})}{\sum_{j}{\cal
P}_{i_{0}j}(x;x,\kappa_{1}^{2})}\,,$ (150)

with $r$ a uniform random number; in this way, the probability of obtaining a
given $i_{1}$ is equal to:

$\frac{{\cal
P}_{i_{0}i_{1}}(x;x,\kappa_{1}^{2})}{\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})+\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})}\,,$
(151)

owing to eq. (149). Next, we determine $y_{1}$ by solving for it the equation:

$1-r=\frac{{\cal P}_{i_{0}i_{1}}(y_{1};x,\kappa_{1}^{2})}{{\cal
P}_{i_{0}i_{1}}(x;x,\kappa_{1}^{2})}\,,$ (152)

with $r$ a uniform random number. Thus, the probability distribution
associated with $y_{1}$ is:

$\displaystyle-\frac{\partial}{\partial y_{1}}\,\frac{{\cal
P}_{i_{0}i_{1}}(y_{1};x,\kappa_{1}^{2})}{{\cal
P}_{i_{0}i_{1}}(x;x,\kappa_{1}^{2})}=$ (153)
$\displaystyle\phantom{aaaaaa}\frac{1}{{\cal
P}_{i_{0}i_{1}}(x;x,\kappa_{1}^{2})}\frac{1}{y_{1}}\left[({\mathbb{O}}^{\scriptscriptstyle\rm
OUT})_{i_{0}i_{1}}\left(\frac{x}{y_{1}},\kappa_{1}^{2}\right)+({\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN})_{i_{0}i_{1}}\left(\frac{x}{y_{1}},\kappa_{1}^{2}\right)\right]f_{i_{1}}(y_{1},\kappa_{1}^{2})\,.\phantom{aa}$

We finally obtain the sought probability by multiplying the results of eqs.
(147), (151), (153), and (142) (with $\mu^{2}\to\kappa_{1}^{2}$, $i_{0}\to
i_{1}$, and $x\to y_{1}$ in the latter), by integrating over all possible
intermediate scales $\kappa_{1}^{2}$ and momentum fractions $y_{1}$, and by
summing over all possible parton types $i_{1}$:

$\displaystyle p_{1}$ $\displaystyle=$
$\displaystyle\sum_{i_{1}}\int_{0}^{1}dy_{1}\int_{\log\mu_{0}^{2}}^{\log\mu^{2}}d\log\kappa_{1}^{2}$
(154)
$\displaystyle\phantom{aaa}\times\frac{1}{f_{i_{0}}(x,\mu^{2})}\frac{S_{i_{0}}(\mu^{2})}{S_{i_{0}}(\kappa_{1}^{2})}\Big{[}\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})+\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})\Big{]}$
$\displaystyle\phantom{aaa}\times\frac{{\cal
P}_{i_{0}i_{1}}(x;x,\kappa_{1}^{2})}{\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})+\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})}$
$\displaystyle\phantom{aaa}\times\frac{1}{{\cal
P}_{i_{0}i_{1}}(x;x,\kappa_{1}^{2})}\frac{1}{y_{1}}\left[({\mathbb{O}}^{\scriptscriptstyle\rm
OUT})_{i_{0}i_{1}}\left(\frac{x}{y_{1}},\kappa_{1}^{2}\right)+({\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN})_{i_{0}i_{1}}\left(\frac{x}{y_{1}},\kappa_{1}^{2}\right)\right]f_{i_{1}}(y_{1},\kappa_{1}^{2})$
$\displaystyle\phantom{aaa}\times\frac{S_{i_{1}}(\kappa_{1}^{2})}{S_{i_{1}}(\mu_{0}^{2})}\,\frac{f_{i_{1}}(y_{1},\mu_{0}^{2})}{f_{i_{1}}(y_{1},\kappa_{1}^{2})}\,.$

Therefore, from eq. (A):

$p_{1}=\frac{f_{i_{0}}^{(1)}(x,\mu^{2})}{f_{i_{0}}(x,\mu^{2})}\,.$ (155)

This procedure can manifestly be iterated, to obtain:

$p_{k}=\frac{f_{i_{0}}^{(k)}(x,\mu^{2})}{f_{i_{0}}(x,\mu^{2})}\;\;\;\;\Longrightarrow\;\;\;\;\sum_{k=0}^{\infty}p_{k}=1\,,$
(156)

having exploited the inverse Mellin transform of eq. (128).

Equation (156) shows that an evolution generated by means of ${\rm NEP}^{({\rm
R})}$ of eq. (60) and of the functions of eq. (148) for the backward steps in
the scale and $x$ spaces, respectively, allows one to recover the PDF used to
guide the evolution. However, it should be clear that this conclusion is
affected by a number of fallacies, that have to do with the NEP possibly being
non-monotonic and not bounded from above by one, as well as the probabilities
of eq. (151) possibly being negative. Both aspects have ultimately to do with
the fact that ${\rm NEP}^{({\rm R})}$ of eq. (60) also accounts for non-
resolvable contributions; as such, it is consistent that it be in agreement
(although only formally) in the sense of eq. (156) with the PDF, whose form is
determined by both resolved and non-resolved contributions. Moreover, note
that eq. (148) is not what is current employed in practical MC
implementations, which rather corresponds to that form with the
${\mathbb{O}}^{\scriptscriptstyle\rm OUT}$ contribution removed (see e.g. eq.
(96)): the proof above shows that, by doing so, one does not recover the PDF
after the evolution.

One can repeat this procedure by adopting the true NEP of either eq. (57) or
eq. (58) (the two coincide). The analogue of eq. (147) is (see eq. (63)):

$\frac{\partial}{\partial\log\kappa_{1}^{2}}\,{\rm
NEP}_{i_{0}}(x,\kappa_{1}^{2},\mu^{2})=\frac{1}{f_{i_{0}}(x,\mu^{2})}\frac{S_{i_{0}}(\mu^{2})}{S_{i_{0}}(\kappa_{1}^{2})}\,\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})\,.$
(157)

Because of this result, the analogues of the functions ${\cal P}_{ij}$ to be
employed in this case are obtained from those in eq. (148) by removing the
${\mathbb{O}}^{\scriptscriptstyle\rm OUT}$ contribution there. By doing so,
one arrives at the analogue of eq. (156), which reads:

$p_{k}=\frac{1}{f_{i_{0}}(x,\mu^{2})}\left(\left.f_{i_{0}}^{(k)}(x,\mu^{2})\right|_{{\mathbb{O}}^{\scriptscriptstyle\rm
OUT}\to 0}\right)\;\;\;\;\Longrightarrow\;\;\;\;\sum_{k=0}^{\infty}p_{k}\neq
1\,.$ (158)

This result need not be surprising: ${\rm NEP}_{i}$ correctly accounts for
resolved emissions only, while as was already said an actual PDF includes non-
resolved contributions. It may appear that the PDF could be recovered in the
context of this evolution by including a branching-by-branching correction
factor equal to:

$\left[({\mathbb{O}}^{\scriptscriptstyle\rm
OUT})_{i_{p-1}i_{p}}\left(\frac{y_{p-1}}{y_{p}},\kappa_{p}^{2}\right)+({\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN})_{i_{p-1}i_{p}}\left(\frac{y_{p-1}}{y_{p}},\kappa_{p}^{2}\right)\right]\Bigg{/}({\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN})_{i_{p-1}i_{p}}\left(\frac{y_{p-1}}{y_{p}},\kappa_{p}^{2}\right)\,,$ (159)

as well as correcting by means of the ratio ${\rm NEP}_{i_{0}}^{({\rm
R})}/{\rm NEP}_{i_{0}}$ the zero-emission contribution. Unfortunately, eq.
(159) does not work: ${\mathbb{O}}^{\scriptscriptstyle\rm OUT}$ and
${\mathbb{O}}_{R}^{\scriptscriptstyle\rm IN}$ have non-overlapping supports in
the $x$ space, and here $y_{p}$ has been generated by using only the latter
operator; it follows that the ratio of eq. (159) is in practice always equal
to one.

Finally, when adopting ${\rm NEP}^{({\rm E})}$ (61) for the NEP, the analogue
of eq. (147) reads as follows (see eq. (65)):

$\displaystyle\frac{\partial}{\partial\log\kappa_{1}^{2}}\,{\rm
NEP}_{i_{0}}^{({\rm E})}(x,\kappa_{1}^{2},\mu^{2})=$ (160)
$\displaystyle\phantom{aaaaaa}\frac{1}{f_{i_{0}}(x,\mu^{2})}\frac{S_{i_{0}}(\mu^{2})}{S_{i_{0}}(\kappa_{1}^{2})}\,\big{(}{\mathbb{Z}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa_{1}^{2})\,\exp\\!\left[\int_{\kappa_{1}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{1}{f_{i_{0}}(x,\kappa^{2})}\,\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i_{0}}(x,\kappa^{2})\right],$

having employed eq. (32). The similarity of this result with that of eq. (157)
suggests that also in this case one can obtain a PDF that stems from keeping
only resolved emissions by including a branching-by-branching correction
factor equal to:

$\exp\left[\int_{\kappa_{p}^{2}}^{\kappa_{p-1}^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{1}{f_{i_{p-1}}(y_{p-1},\kappa^{2})}\,\big{(}{\mathbb{W}}\left[F\right]\\!\big{)}_{i_{p-1}}(y_{p-1},\kappa^{2})\right].$
(161)

Since ${\mathbb{W}}[F]$ has no definite sign, this factor can be larger or
smaller than one. This is connected with the fact that the ${\cal
O}(\alpha_{\scriptscriptstyle S}^{2})$ coefficient in eq. (68) has no definite
sign, and its cumulative effect over successive backward branchings is
therefore typically smaller than naive coupling-constant power counting would
suggest, a point borne out in practice by the results shown in fig. 1.

##  Appendix B An academic model

A different perspective on the three forms of NEP considered in this paper can
be obtained in the context of an academic model, defined so that the only
branchings are of virtual origin. This can be achieved by setting:

${\mathbb{A}}={\mathbb{C}}=0\,.$ (162)

Before proceeding, we stress that eq. (162) defines the virtual contribution
in an unique manner only because one understands eq. (3) so that, from eq.
(15)

$\overline{{\mathbb{B}}}={\mathbb{B}}\,.$ (163)

A different model which still sets the real splitting kernels equal to zero
can be defined as follows:

$\widetilde{{\mathbb{A}}}={\mathbb{C}}=0$ (164)

which understands the form of eq. (5), so that from eq. (16):

$\overline{{\mathbb{B}}}=\widetilde{{\mathbb{B}}}\,.$ (165)

Thus: ${\mathbb{A}}=0$ implies eq. (163), while $\widetilde{{\mathbb{A}}}=0$
implies eq. (165), with ${\mathbb{B}}$ and $\widetilde{{\mathbb{B}}}$ still
related to one another by eq. (6), with $\widetilde{{\mathbb{A}}}\neq 0$
there. For example, in the case of the gluon at the LO in QCD:

$\displaystyle
A_{g}(z)=0\;\;\;\;\Longrightarrow\;\;\;\;\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\,\overline{B}_{g}=-\frac{C_{\scriptscriptstyle A}+4T_{\scriptscriptstyle
F}N_{\scriptscriptstyle F}}{6}\,,$ (166)
$\displaystyle\widetilde{A}_{g}(z)=0\;\;\;\;\Longrightarrow\;\;\;\;\frac{2\pi}{\alpha_{\scriptscriptstyle
S}}\,\overline{B}_{g}=\gamma(g)\,.$ (167)

Having clarified this point, we can solve the PDF evolution equations and
consider the MC-generated backward evolution for the model of eq. (162). We do
so by employing the parameter $\lambda$ introduced in eq. (53). We can solve
directly the PDF evolution equations, and obtain:

$F(\mu^{2})=\frac{{\mathbb{S}}_{\lambda=0}(\mu^{2})}{{\mathbb{S}}_{\lambda=0}(\mu_{0}^{2})}\,F(\mu_{0}^{2})\,.$
(168)

With eq. (31) we have instead

$F(\mu^{2})=\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,F(\mu_{0}^{2})+\lambda\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\kappa^{2})}\,\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}(\kappa^{2})F(\kappa^{2})\,.$ (169)

These two solutions coincide (as they should), since by using eq. (168) in the
second term on the r.h.s. of eq. (169) one obtains:

$\lambda\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\kappa^{2})}\,\overline{{\mathbb{B}}}^{\scriptscriptstyle\rm
OUT}(\kappa^{2})F(\kappa^{2})=\frac{{\mathbb{S}}_{\lambda=0}(\mu^{2})}{{\mathbb{S}}_{\lambda=0}(\mu_{0}^{2})}\,F(\mu_{0}^{2})-\frac{{\mathbb{S}}(\mu^{2})}{{\mathbb{S}}(\mu_{0}^{2})}\,F(\mu_{0}^{2})\,,$
(170)

which shows explicitly, in this simple case, the cancellation of dependence on
$\lambda$ on the r.h.s. of eq. (31).

We now consider an MC backward evolution. We start by noting that:

${\rm NEP}_{i}=1\,,\;\;\;\;\;\;\;\;{\rm NEP}_{i}^{({\rm
E})}=1\,,\;\;\;\;\;\;\;\;{\rm NEP}_{i}^{({\rm
R})}=\frac{S_{i,\lambda=0}(\mu_{0}^{2})}{S_{i,\lambda=0}(\mu^{2})}\,\frac{S_{i}(\mu^{2})}{S_{i}(\mu_{0}^{2})}\,.$
(171)

The leftmost result in eq. (171) is what we expect in view of the physical
interpretation of ${\rm NEP}_{i}$: only strictly resolved emissions may
contribute to it, and in this model no resolved emissions can be generated –
thus, the NEP must be equal to one. The middle result in eq. (171) shows that
this model is too simple to allow one to distinguish the behaviour of ${\rm
NEP}_{i}^{({\rm E})}$ from that of ${\rm NEP}_{i}$: the spurious terms of
${\mathbb{W}}[F]$ origin potentially present in the former case according to
eq. (65) are all identically equal to zero, being proportional to
${\mathbb{Z}}[F]=0$. Therefore, also in this case the NEP is equal to one.
Finally, the rightmost result in eq. (171) shows that for any $\lambda\neq 0$
the NEP is not equal to one since it receives, independently from one another,
both resolved and unresolved contributions, and the latter are different from
zero owing to the virtual contribution to ${\mathbb{W}}[F]$ (see eq. (40)).
Having said that, we point out that, in view of eq. (168), a sound
probabilistic interpretation of ${\rm NEP}_{i}^{({\rm R})}$ requires that:

$S_{i}(\mu^{2})<S_{i}(\mu_{0}^{2})\;\;\;\;\Longleftrightarrow\;\;\;\;\overline{B}_{i}^{\scriptscriptstyle\rm
IN}\equiv\overline{B}_{i}-\lambda\overline{B}_{i}^{\scriptscriptstyle\rm
OUT}<0\,,$ (172)

for any $\mu^{2}>\mu_{0}^{2}$. Therefore, in the simplest scenario
$\lambda=0$, this happens for the gluon (see eq. (166)) in the context of the
model of eq. (162) we work with; however, this would not happen had we chosen
the different model of eq. (164) (see eq. (167)). This simple example confirms
a general fact that has already been inferred before, namely that the
interpretation of ${\rm NEP}_{i}^{({\rm R})}$ as a NEP might become
problematic.

If we want to obtain the Neumann summands of eq. (A) relevant to the model of
eq. (162) we need to use the fact that:

$({\mathbb{O}}^{\scriptscriptstyle\rm
OUT})_{ij}\left(z,\kappa^{2}\right)=\lambda\delta_{ij}\overline{B}_{i}^{\scriptscriptstyle\rm
OUT}(\kappa^{2})\delta(1-z)\,,\;\;\;\;\;\;\;\;({\mathbb{O}}_{R}^{\scriptscriptstyle\rm
IN})_{ij}\left(z,\kappa^{2}\right)=0\,.$ (173)

The absence of off-diagonal terms leads to an immediate dramatic
simplification of eq. (A), which becomes:

$f_{i_{0}}^{(k)}(x,\mu^{2})=\frac{S_{i_{0}}(\mu^{2})}{S_{i_{0}}(\mu_{0}^{2})}\,f_{i_{0}}(x,\mu_{0}^{2})\,\frac{\lambda^{k}}{k!}\int_{\mu_{0}^{2}}^{\mu^{2}}\frac{d\kappa^{2}}{\kappa^{2}}\overline{B}_{i_{0}}^{\scriptscriptstyle\rm
OUT}(\kappa^{2})\,,$ (174)

a result that is also valid for $k=0$. By summing over $k$ one finds again the
solution of eq. (168), as one must by construction. Moreover, as we have
previous learned, eq. (174) can be seen as the MC contribution to the PDFs due
to showers that feature $k$ emissions. However, for this to be true in the
model defined by eq. (162), one would have to have chosen ${\rm
NEP}_{i}^{({\rm R})}$ as the NEP, since that is the only nontrivial NEP in
this context (see eq. (171)). Therefore, this simple example confirms the
previous general findings, namely that while ${\rm NEP}_{i}^{({\rm R})}$ is
not the correct non-emission probability, it nevertheless formally allows one
to recover the PDF given in input. Conversely, if ${\rm NEP}_{i}$ (or ${\rm
NEP}_{i}^{({\rm E})}$) had been adopted, both matrix elements in the analogue
of eq. (173) would be equal to zero, leading to a Neumann series whose terms
would be all equal to zero bar the first. The latter then coincides with the
reconstructed PDF, and reads as follows:

$f_{i_{0}}^{(0)}(x,\mu^{2})=\frac{S_{i_{0}}(\mu^{2})}{S_{i_{0}}(\mu_{0}^{2})}\,f_{i_{0}}(x,\mu_{0}^{2})\,.$
(175)

This is in general different from the exact solution of eq. (168). Clearly,
the model of eq. (162) is maximally perverse, since all emissions are
unresolved; it is therefore not particularly surprising that evolutions based
on NEPs that can only account for resolved emissions fail.

## References

  * (1) S. Alioli et al., _Monte Carlo event generators for high energy particle physics event simulation_ , 1902.01674.
  * (2) G. Marchesini and B. R. Webber, _Monte Carlo Simulation of General Hard Processes with Coherent QCD Radiation_ , _Nucl. Phys. B_ 310 (1988) 461–526.
  * (3) T. Sjostrand, _A Model for Initial State Parton Showers_ , _Phys. Lett. B_ 157 (1985) 321–325.
  * (4) V. N. Gribov and L. N. Lipatov, _Deep inelastic e p scattering in perturbation theory_ , _Sov. J. Nucl. Phys._ 15 (1972) 438–450.
  * (5) G. Altarelli and G. Parisi, _Asymptotic Freedom in Parton Language_ , _Nucl. Phys._ B126 (1977) 298–318.
  * (6) Y. L. Dokshitzer, _Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics._ , _Sov. Phys. JETP_ 46 (1977) 641–653.
  * (7) M. Yan, T.-J. Hou, P. Nadolsky and C. P. Yuan, _CT18 global PDF fit at leading order in QCD_ , _Phys. Rev. D_ 107 (2023) 116001, [2205.00137].
  * (8) M. Bengtsson, T. Sjostrand and M. van Zijl, _Initial State Radiation Effects on $W$ and Jet Production_, _Z. Phys. C_ 32 (1986) 67.
  * (9) M. Skrzypek, S. Jadach, K. Golec-Biernat and W. Placzek, _QCD ISR Monte Carlo with alpha(k(T))_ , _Acta Phys. Polon. B_ 38 (2007) 2369–2378.
  * (10) R. Frederix, S. Frixione, S. Prestel and P. Torrielli, _On the reduction of negative weights in MC@NLO-type matching procedures_ , _JHEP_ 07 (2020) 238, [2002.12716].
  * (11) S. Frixione, Z. Kunszt and A. Signer, _Three jet cross-sections to next-to-leading order_ , _Nucl. Phys. B_ 467 (1996) 399–442, [hep-ph/9512328].
  * (12) S. Frixione, _A General approach to jet cross-sections in QCD_ , _Nucl. Phys. B_ 507 (1997) 295–314, [hep-ph/9706545].