aainstitutetext:  Università di Milano-Bicocca and INFN, Sezione di Milano-
Bicocca, Piazza della Scienza 3, 20126 Milano, Italybbinstitutetext: INFN,
Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano,
Italyccinstitutetext: Rudolf Peierls Centre for Theoretical Physics,
University of Oxford, Clarendon Laboratory, Parks Road, Oxford OX1 3PU, United
Kingdomddinstitutetext:  All Souls College, Oxford OX1 4AL, United Kingdom

# Multiplicative-Accumulative matching of NLO calculations with parton showers

Paolo Nason c,d Gavin P. Salam

(Received: date / Accepted: date)

###### Abstract

We propose a new approach for combining next-to-leading order (NLO) and parton
shower (PS) calculations so as to obtain three core features: (a)
applicability to general showers, as with the MC@NLO and POWHEG methods; (b)
positive-weight events, as with the KrkNLO and POWHEG methods; and (c) all
showering attributed to the parton shower code, as with the MC@NLO and KrkNLO
methods. This is achieved by using multiplicative matching in phase space
regions where the shower overestimates the matrix element and accumulative
(additive) matching in regions where the shower underestimates the matrix
element, an approach that can be viewed as a combination of the MC@NLO and
KrkNLO methods.

††preprint: OUTP-21-26P

##  1 Introduction

Next-to-leading (NLO) order event generators interfaced to Parton Showers
(NLOPS from now on) have, in the past decades, become the state of the art for
the simulation of hard collider processes. The MC@NLO algorithm was the first
one to be proposed Frixione:2002ik . Soon after, the POWHEG method was
proposed Nason:2004rx  in order to overcome the problem of negative weights
inherent to MC@NLO. Novel methods were proposed later Platzer:2012bs ;
Lonnblad:2012ix ; Jadach:2015mza . Among them, the KrkNLO method, like the
POWHEG method, has the characteristic of being free from negative weights,
although its applicability is at the moment restricted to relatively simple
processes.

Both the MC@NLO and the KrkNLO methods are intimately intertwined with the
parton shower generator that is adopted. In fact, the two methods provide the
corrections one must apply to the parton shower result in order to achieve
next-to-leading order (NLO) accuracy. This is at variance with POWHEG, which
is largely independent from the shower generator. In fact, POWHEG takes care
of the generation of the hardest event in such a way that NLO accuracy is
preserved, while the PS takes care of the remaining (less hard) radiation.
This feature is convenient, since one can generate events in the Les Houches
format Boos:2001cv  and then shower them with any available PS generator that
complies with the Les Houches format specifications.

In view of recent developments aimed at the improvement of the logarithmic
accuracy of the shower (see Refs. Dasgupta:2018nvj ; Dasgupta:2020fwr ;
Hamilton:2020rcu ; Karlberg:2021kwr , Forshaw:2019ver ; Forshaw:2020wrq ;
Holguin:2020joq  and Nagy:2020rmk ; Nagy:2020dvz ) it is legitimate to ask
whether POWHEG will maintain the same advantage, i.e. whether minor
modifications to the POWHEG algorithm will be sufficient to guarantee that,
when interfaced to an already NLL accurate shower, NLL accuracy will be
maintained, or whether it will be more feasible to correct an already NLL
accurate shower to achieve also NLO accuracy, or rather whether both
alternatives will be feasible, and will productively compete among each other.
It is therefore important to explore the full range of options for matching
NLO calculations and parton showers.

In this note, we will address certain limitations of MC@NLO and KrkNLO, and
show how they can be overcome. The MC@NLO approach is widely applicable, and
widely used, but has the undesirable feature of negative weights.111Recent
proposals for the reduction of the negative-weight fraction include Refs.
Frederix:2020trv ; Andersen:2021mvw ; Danziger:2021xvr . On the other hand,
the KrkNLO method has positive weights, but it is difficult to extend it to
generic processes. We will show that combining the two methods one could
achieve both positive weights and unrestricted applicability.

The paper is organised as follows, in section 2 we will compare the POWHEG,
MC@NLO and KrkNLO methods by formulating them in a common language. We will do
this by extending the formulation of the POWHEG and MC@NLO methods of ref.
Nason:2012pr  also to the KrkNLO method. In section 3 we will present the
combined method, and in section 4 we will present some of its possible
variants. Finally, in section 5 we present our conclusions.

##  2 POWHEG, MC@NLO and KrkNLO

We assume that the phase space with radiation $\Phi$ can be written in terms
of an underlying Born phase space $\Phi_{\rm B}$ and three radiation
variables, indicated collectively as $\Phi_{\rm rad}$. We also assume that the
mapping from $\Phi$ to $\Phi_{\rm B}$ is such that in the singular (collinear
or soft) limit, the Born configuration matches the full phase space with the
collinear pair merged into a single parton or with the soft particle removed.
Such mappings are easy to realise for processes with a single singular region,
while for more complex processes one must separate the real cross section into
contributions having a single singular region. In the following illustrative
discussion we ignore these complications.

We define the following quantity, a function of the underlying Born
kinematics,

$\bar{B}_{\text{\sc s}}(\Phi_{\rm B})=B_{0}(\Phi_{\rm B})+V(\Phi_{\rm B})+\int
R_{\text{\sc s}}(\Phi)\mathrm{d}\Phi_{\rm rad},$ (1)

where the full phase space $\Phi$ is defined as function of $\Phi_{\rm B}$ and
$\Phi_{\rm rad}$, $B_{0}$ is the Born cross section, $V$ comprises the virtual
corrections and, for hadron initiated processes, the collinear counterterms
integrated at fixed underlying Born, and $R_{\text{\sc s}}{}$ is a part of the
real cross section that includes all soft and collinear singularities. In
other words, $R-R_{\text{\sc s}}{}$ is non-singular. Notice that if
$R_{\text{\sc s}}{}$ is taken equal to $R$, $\bar{B}_{\text{\sc s}}{}$ is the
inclusive cross section at fixed underlying Born. We imagine that
renormalisation has been carried out, and that the infrared divergences
arising in $V$ and in the $\mathrm{d}\Phi_{\rm rad}$ integral of $R_{s}$ are
regularised in some way. Notice also that $\bar{B}_{\text{\sc s}}{}$ is
finite, since the singularities present in $V$ cancel those arising
integrating the $R_{\text{\sc s}}$ term.

We also define the Sudakov form factor associated with $R_{\text{\sc s}}{}$

$\displaystyle S(t,\Phi_{\rm B})$
$\displaystyle=\exp\left[-\int_{t_{\Phi}>t}\frac{R_{\text{\sc
s}}(\Phi)}{B_{0}(\Phi_{\rm B})}\mathrm{d}\Phi_{\rm rad}\right],$ (2)

where $\Phi$ is defined in terms of the variables $\Phi_{\rm B}$ and
$\Phi_{\rm rad}$, and $t$ is some definition of hardness, depending upon the
full phase space with radiation. One may think, for example, that $t$ is the
relative transverse momentum of the splitting pair. If $R_{\text{\sc s}}$ was
taken equal to $R$, the Sudakov form factor would represent the probability
for not radiating anything harder than $t$.

The hardest event cross section can be represented in both POWHEG and MC@NLO
as

$\mathrm{d}\sigma=\bar{B}_{\text{\sc s}}{}(\Phi_{\rm B})S(t_{\rm
cut},\Phi_{\rm B})\mathrm{d}\Phi_{\rm B}+\bar{B}_{\text{\sc s}}{}(\Phi_{\rm
B})\,S(t_{\Phi},\Phi_{\rm B})\times\frac{R_{\text{\sc
s}}(\Phi)}{B_{0}(\Phi_{\rm B})}\,\theta(t_{\Phi}-t_{\rm
cut})\mathrm{d}\Phi\,+\\\ +\left[R(\Phi)-R_{\text{\sc
s}}{}(\Phi)\right]\mathrm{d}\Phi,$ (3)

where $\mathrm{d}\Phi=\mathrm{d}\Phi_{\rm B}\,\mathrm{d}\Phi_{\rm rad}$, and
$t_{\rm cut}$ represents a lower limit for radiation, needed to avoid the
Landau-pole singularities. Events are generated with a probability
proportional to each term of the cross section formula. The first two terms
are generated with a Monte Carlo technique. In fact, they satisfy the shower
unitarity equation

$S(t_{\rm cut},\Phi_{\rm B})+\int_{t_{\Phi}>t_{\rm cut}}S(t_{\Phi},\Phi_{\rm
B})\times\frac{R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm B})}\mathrm{d}\Phi_{\rm
rad}=1,$ (4)

which follows from the fact that the expression under the integral sign is an
exact differential. In POWHEG, the generation of the hardness $t$ is uniform
in the Sudakov form factor, and one can generate events with the standard
shower algorithm by equating a random number with the Sudakov form factor, and
solving for $t$. If the $t$ value so obtained is above $t_{\rm cut}$, the
radiation kinematics is generated, and the event with the hardest radiation is
fed to a shower generator, which takes care of adding subsequent (less hard)
radiation.

In the case of MC@NLO, the implementation of the first term is generally more
involved. If the shower is ordered in transverse momentum, the hardest
emission is the first, and the radiation mechanism is the same as in POWHEG,
except that it is implemented within the shower generator, rather than in the
NLO program. In case of an angular ordered shower, large angle soft radiation
is generated first, and the hardest radiation occurs somewhere down the
shower.222Taking $t$ to be the angular scale of the shower means that $t$ does
not represent a hardness. As a result, in the absence of an infrared cutoff,
the first emission is dominated by the infinitely soft region (where
$R_{\text{\sc s}}=R$) rather than by the hard region. It was shown in ref.
Nason:2004rx  that, in this case, by suitable rearrangement of the Sudakov
factors for each emission, one reconstructs the transverse momentum Sudakov
form factor in formula (3).

The last term in Eq. (3) is non-singular, and thus is dominated by hard
radiation. It is handled essentially in the same way in MC@NLO and POWHEG.

The NLO accuracy of Eq. (3) can be demonstrated by computing the expectation
value of a generic infrared safe observable $O(\Phi)$ as follows

$\displaystyle\langle O\rangle$ $\displaystyle=\int\mathrm{d}\Phi_{\rm
B}\,\bar{B}_{\text{\sc s}}{}(\Phi_{\rm B})\Bigg{\\{}S(t_{\rm cut},\Phi_{\rm
B})\,O(\Phi_{\rm B})+\int_{t_{\Phi}>t_{\rm cut}}S(t_{\Phi},\Phi_{\rm
B})\times\frac{R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm
B})}O(\Phi)\mathrm{d}\Phi_{\rm rad}\Bigg{\\}}\,+$ $\displaystyle\hskip
160.00024pt\qquad\qquad+\int\left[R(\Phi)-R_{\text{\sc
s}}{}(\Phi)\right]O(\Phi)\mathrm{d}\Phi,$ (5a) or equivalently,
$\displaystyle\langle O\rangle$ $\displaystyle=\int\mathrm{d}\Phi_{\rm
B}\bar{B}_{\text{\sc s}}{}(\Phi_{\rm B})\Bigg{\\{}S(t_{\rm cut},\Phi_{\rm
B})\,O(\Phi_{\rm B})+\int_{t_{\Phi}>t_{\rm cut}}S(t_{\Phi},\Phi_{\rm
B})\times\frac{R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm B})}O(\Phi_{\rm
B})\mathrm{d}\Phi_{\rm rad}\,+$ $\displaystyle\hskip
120.00018pt+\int_{t_{\Phi}>t_{\rm cut}}S(t_{\Phi},\Phi_{\rm
B})\times\frac{R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm
B})}[O(\Phi)-O(\Phi_{\rm B})]\mathrm{d}\Phi_{\rm rad}\Bigg{\\}}\,+$
$\displaystyle\hskip 160.00024pt\qquad\qquad+\int\left[R(\Phi)-R_{\text{\sc
s}}{}(\Phi)\right]O(\Phi)\mathrm{d}\Phi\,.$ (5b)

In view of the unitarity equation (4), the first two terms in the curly
bracket of equation (5) collapse into $O(\Phi_{\rm B})$. Furthermore, the
factor $[O(\Phi)-O(\Phi_{\rm B})]$ in the third term in the curly bracket
suppresses the singular region, so that up to higher order terms the Sudakov
form factor can be omitted, and the lower integration limit can be set to zero
up to a power-suppressed correction. We thus get

$\displaystyle\langle O\rangle$ $\displaystyle=\int\mathrm{d}\Phi_{\rm
B}\bar{B}_{\text{\sc s}}{}(\Phi_{\rm B})\Bigg{\\{}O(\Phi_{\rm
B})+\int\frac{R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm B})}[O(\Phi)-O(\Phi_{\rm
B})]\mathrm{d}\Phi_{\rm rad}\Bigg{\\}}\,+$ $\displaystyle\hskip
120.00018pt+\int\left[R(\Phi)-R_{\text{\sc
s}}{}(\Phi)\right]O(\Phi)\mathrm{d}\Phi+{\cal O}({\rm NNLO})\,,$ (6a)
$\displaystyle=\int\mathrm{d}\Phi_{\rm B}\left\\{\bar{B}_{\text{\sc
s}}{}(\Phi_{\rm B})O(\Phi_{\rm B})+\int R_{\text{\sc
s}}(\Phi)[O(\Phi)-O(\Phi_{\rm B})]\mathrm{d}\Phi_{\rm rad}\right\\}\,+$
$\displaystyle\hskip 120.00018pt+\int\left[R(\Phi)-R_{\text{\sc
s}}{}(\Phi)\right]O(\Phi)\mathrm{d}\Phi+{\cal O}({\rm NNLO})\,.$ (6b)

Inserting into the above equation the expression for $\bar{B}_{\text{\sc s}}$,
Eq. (1), we get

$\displaystyle\langle O\rangle=\int\mathrm{d}\Phi_{\rm B}[B(\Phi_{\rm
B})+V(\Phi_{\rm B})]O(\Phi_{\rm B})+\int R(\Phi)O(\Phi)\mathrm{d}\Phi+{\cal
O}({\rm NNLO}),$ (7)

which is the correct NLO expression for the observable.

POWHEG implements Eq. (3) directly. $R_{\text{\sc s}}{}$ is defined as

$R_{\text{\sc s}}{}(\Phi)=F(\Phi)R\,,$ (8)

where $F(\Phi)\leq 1$, and approaches 1 as $\Phi$ approaches the singular
region, so that $R_{\text{\sc s}}{}$ carries all the singularity structure of
$R$.

In MC@NLO, $R_{\text{\sc s}}{}$ is the shower approximation to $R$, typically
given by the Born term times a DGLAP splitting function. Thus, in MC@NLO one
computes a Born configuration using the $\bar{B}_{\text{\sc s}}{}$ function,
and passes it directly to the Shower Monte Carlo. The $R-R_{\text{\sc s}}{}$
terms are generated separately, and fed directly to the shower. Negative
weights can arise at this stage, since there is no guarantee that
$R-R_{\text{\sc s}}{}$ is positive.

Using the same language adopted here for POWHEG and MC@NLO, the key formula
for the KrkNLO method can be written as follows

$\mathrm{d}\sigma=\bar{B}_{\text{\sc s}}{}(\Phi_{\rm B})S(t_{\rm
cut},\Phi_{\rm B})\mathrm{d}\Phi_{\rm B}+\bar{B}_{\text{\sc s}}{}(\Phi_{\rm
B})\left\\{S(t_{\Phi},\Phi_{\rm B})\times\frac{R_{\text{\sc
s}}(\Phi)}{B_{0}(\Phi_{\rm
B})}\right\\}\times\left[\frac{R(\Phi)}{R_{\text{\sc
s}}(\Phi)}\right]\mathrm{d}\Phi\,,$ (9)

(with an implicit $\theta(t-t_{\rm cut})$ in $d\Phi$). In the literature
dealing with the KrkNLO method Jadach:2017ujd ; Jadach:2016qti ;
Jadach:2016viv ; Jadach:2016acv ; Jadach:2015mza ; Jadach:2015zsq ;
Jadach:2020xfl , particular emphasis has been put on the use of a specific
scheme for the parton densities, which considerably simplifies the expression
for the $\bar{B}_{\text{\sc s}}{}$ function. Here we are instead interested in
a simpler aspect of the method, which is that to generate NLO accurate
radiation, it uses a multiplicative correction, rather than the additive
correction of the MC@NLO method. The NLO accuracy of the KrkNLO formula can be
simply demonstrated by showing that KrkNLO is equivalent to MC@NLO at the NLO
level. In fact, we can rewrite formula (9) as

$\mathrm{d}\sigma=\bar{B}_{\text{\sc s}}{}(\Phi_{\rm B})S(t_{\rm
cut},\Phi_{\rm B})\mathrm{d}\Phi_{\rm B}+\bar{B}_{\text{\sc s}}{}(\Phi_{\rm
B})S(t_{\Phi},\Phi_{\rm B})\times\frac{R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm
B})}\times\left[\frac{R(\Phi)}{R_{\text{\sc
s}}(\Phi)}-1\right]\mathrm{d}\Phi+\\\ +\bar{B}_{\text{\sc s}}{}(\Phi_{\rm
B})S(t_{\Phi},\Phi_{\rm B})\frac{R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm
B})}\mathrm{d}\Phi\,.$ (10)

The middle term is now insensitive to the soft region, because the factor in
square brackets vanishes there, so we can drop the Sudakov form factor and,
neglecting terms of NNLO size, set $\bar{B}_{\text{\sc s}}{}/B_{0}=1$. By
doing this we recover exactly Eq. (3).

The KrkNLO method leads to positive weighted events. On the other hand, unlike
MC@NLO and POWHEG, its cross section at fixed underlying Born does not exactly
match the corresponding fixed order result, but differs from it by NNLO
terms.333Although some authors consider this to be an undesirable feature, it
does not constitute a real problem, since choices for uncontrolled NNLO terms
are made, for example, when choosing the scales, and there is no preferred way
to define an NLO result as far as the neglected NNLO terms are concerned. The
KrkNLO method, however, generates weighted events, so, the unweighting
efficiency may constitute a problem if one wants to apply the method to
generic processes without having to do process-by-process adjustments. As a
related problem, the shower generator may not cover the full radiation phase
space. This is the same as saying that $R_{\text{\sc s}}{}$ can become zero in
certain regions, in which case the method is not applicable. For these
reasons, it seems difficult to apply the method to generic processes in
automated framework, something that has been available in POWHEG and MC@NLO
since a long time.

##  3 The new method

By comparing the MC@NLO and KrkNLO method using a common language, it becomes
clear that the two methods have much in common, and can in fact be merged in
such a way that the KrkNLO positivity is maintained, unweighted events can be
generated on the fly, and no issues arise from the limited coverage of the
phase space by the parton shower code. The merged method is defined by the
following formula

$\mathrm{d}\sigma=\bar{B}_{\text{\sc s}}{}(\Phi_{\rm B})\,S(t_{\Phi},\Phi_{\rm
B})\times\frac{R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm
B})}\times\left\\{1+\frac{R-R_{\text{\sc s}}{}}{R_{\text{\sc
s}}}\theta(R_{\text{\sc s}}-R)\right\\}\mathrm{d}\Phi+\theta(R-R_{\text{\sc
s}})\left[R-R_{\text{\sc s}}\right]\mathrm{d}\Phi\,.$ (11)

It is easy to see by inspection that this formula has the same NLO accuracy as
MC@NLO. In fact, the $\theta$ function in the first term of eq. (11) is
regular, and as such the Sudakov form factor multiplying it and the
$\bar{B}_{\text{\sc s}}{}/B_{0}$ ratio can both be set to 1 up to NNLO
corrections. Having done this, the contribution proportional to the $\theta$
function in the first term becomes identical to the last term, except for the
theta function, which has the opposite sign in the argument. The two terms
thus combine, and the theta function disappears, yielding the MC@NLO formula.

The method of Eq. (11) is immediately seen to have positive weights, since the
second term in eq. (11) is forced to be positive by the theta function, while
the factor in the curly bracket of the first term is forced by the theta
function to be less than 1. Notice that regions of phase space not populated
by the Monte Carlo only contribute to the additive term in the square bracket,
since they have $R_{\text{\sc s}}=0$.

Let us consider how the method of Eq. (11) might be implemented in a scenario
where it is the parton shower that is responsible for generating $R_{\text{\sc
s}}$.444Alternatively, one might also implement the approach in a POWHEG-style
scheme, where the NLO program takes responsibility for the first emission, and
$R_{\text{\sc s}}$ is not required to be smaller than $R$. We first consider
the case of a parton shower ordered in a genuine hardness variable (for
example, transverse momentum). As in the MC@NLO approach, the NLO program is
responsible for generating a sample of Born events, with a distribution
corresponding to the correct $\bar{B}_{\text{\sc s}}{}(\Phi_{\rm
B})\mathrm{d}\Phi_{\rm B}$ weight. The event is communicated to the parton
shower, which generates a first emission, i.e. accounting for a remaining
factor $S(t_{\Phi},\Phi_{\rm B})\frac{R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm
B})}\mathrm{d}\Phi/\mathrm{d}\Phi_{\rm B}$ in Eq. (11). Next, the parton
shower communicates the event, with its first emission, back to the NLO
program, which evaluates the contents of the curly bracket in Eq. (11). The
result of that evaluation, which is bounded to be between zero and one, is to
be used as an acceptance probability for the event. If the event is accepted,
the parton shower then continues from that point, generating the remaining
emissions. In addition, the NLO program is responsible for generating a second
(positive definite) sample of events, corresponding to the
$\theta(R-R_{\text{\sc s}}$) term of Eq. (11), covering the part of the real
phase space where the shower $R_{\text{\sc s}}$ underestimates the true real
matrix element. These events are passed to the parton shower code for normal
showering.

The above scheme requires a more connected workflow between the parton shower
and the NLO code than either of the MC@NLO or POWHEG approaches. Recall that
those approaches pass a Les Houches Event file to the shower code, and then
let the shower take over from there. In contrast, the scheme outlined above
requires additional action from the NLO code _after_ the first parton shower
emission. However, this is simply a technical consideration, which in our view
is a small price to pay for an approach that, like the MC@NLO scheme, leaves
responsibility for the first emission with the parton shower, while
eliminating negative weights.

A further remark concerns angular-ordered parton showers. For such showers,
the application of the acceptance probability (factor in curly brackets in Eq.
(11)) after the first emission would not be correct. Instead, one might
envisage an approach in which one runs the complete shower, uses a jet
algorithm to map the full event to the real phase space and then applies the
acceptance probability.

##  4 Variants of the method

Eq. (11) can be viewed as a special case of a family of approaches,
parameterised by a constant $c\geq 1$,

$\mathrm{d}\sigma=\bar{B}_{\text{\sc s}}{}(\Phi_{\rm B})\,S(t_{\Phi},\Phi_{\rm
B})\times\frac{c\,R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm
B})}\times\left\\{1+\frac{R-c\,R_{\text{\sc s}}{}}{c\,R_{\text{\sc
s}}}\theta(c\,R_{\text{\sc s}}-R)\right\\}\mathrm{d}\Phi\,+\\\
+\theta(R-c\,R_{\text{\sc s}})\left[R-c\,R_{\text{\sc
s}}\right]\mathrm{d}\Phi\,.$ (12)

This formula was obtained by replacing $R_{\text{\sc s}}\to c\,R_{\text{\sc
s}}$ in all occurrences where it appears in Eq. (11), except for the Sudakov
form factor. First of all, we should convince ourselves that this formula has
the correct behaviour near the soft limit, and that it is NLO accurate. This
is seen immediately by writing (12) as

$\mathrm{d}\sigma=\bar{B}_{\text{\sc s}}{}(\Phi_{\rm B})\,S(t_{\Phi},\Phi_{\rm
B})\times\frac{c\,R_{\text{\sc s}}(\Phi)}{B_{0}(\Phi_{\rm
B})}\times\left\\{1+\frac{R-c\,R_{\text{\sc s}}{}}{c\,R_{\text{\sc
s}}}\right\\}\mathrm{d}\Phi\,,\\\ +\theta(R-c\,R_{\text{\sc
s}})(R-c\,R_{\text{\sc s}})\left[1-\frac{\bar{B}_{\text{\sc s}}{}(\Phi_{\rm
B})}{B_{0}(\Phi_{\rm B})}\,S(t_{\Phi},\Phi_{\rm B})\right]\mathrm{d}\Phi\,.$
(13)

The first line is equal to the KrkNLO formula, eq. (9) (leaving aside the
$t_{\rm cut}$, for simplicity). In the second line, the factor in square
brackets is of order $\alpha_{\scriptscriptstyle S}$ and it multiplies a term
dominated by the hard region (and thus also of order
$\alpha_{\scriptscriptstyle S}$), since

$\theta(R-c\,R_{\text{\sc s}})(R-c\,R_{\text{\sc s}})\leq\theta(R-R_{\text{\sc
s}})(R-R_{\text{\sc s}})\,.$

Accordingly, the second line is of NNLO order, and we can conclude that
formula (12) is equivalent to the KrkNLO formula up to NNLO corrections. For a
shower where $R_{\text{\sc s}}/R$ is always larger than some value $r$, taking
$c\geq 1/r$ causes Eq. (12) to reduces to the KrkNLO approach over all of
phase space.

The implementation of formula (12) goes as follows. The soft events are
oversampled by a factor of $c$, and accepted with a probability proportional
to the expression in the curly bracket, while the hard events are generated in
a standard way.

There are several reasons why it may be of use to have such a family of
approaches parameterised by $c$. One is that it is useful to have a parameter
to help gauge systematic uncertainties associated with terms that are beyond
the accuracy of the method. In fact, the $c$ parameter also gauges the amount
of radiation that is multiplied by the inclusive $K$-factor, with respect to
the amount that is added in as hard radiation. In this sense, it would play a
similar role to the hdamp parameter in POWHEG Nason:2004rx ; Alioli:2008tz .
Another reason is that the $\theta$-functions in Eqs. (11), (12) can,
conceivably, introduce non-smoothness in phase-space coverage even when the
underlying $R_{\text{\sc s}}$ and $R$ functions are smooth. Sampling over a
range of $c$ values would allow one to address that issue.555It is not clear
to us that this would be a genuine problem. A worse issue is potentially
present in the MC@NLO approach (and also in Eqs. (11) and (12)) if a shower
has a discontinuity in the distribution of the hardness variable. We are not
aware of this having caused significant problems in practice, possibly because
subsequent showering smoothens any discontinuities. Were it to be a problem,
one solution could be to sample over a range of shower starting scales.

So far, the variants that we have discussed involve the rejection of events.
For showers with an ordering variable that involve a genuine hardness scale,
such as transverse momentum, it is also possible to envisage a variant where
instead of rejecting the event with probability $\frac{R_{\text{\sc
s}}-R}{R_{\text{\sc s}}}\theta(R_{\text{\sc s}}-R)$, one rejects the shower’s
first emission with that same probability. If that first emission is rejected,
at a value of the ordering variable that we label $t_{1}$, the shower then
continues from that scale $t_{1}$, based on an event without the first
emission. As the shower continues, each time the shower again attempts to
create a new first emission, that emission continues to be rejected with
probability $\frac{R_{\text{\sc s}}-R}{R_{\text{\sc s}}}\theta(R_{\text{\sc
s}}-R)$. Once the shower has generated a first emission that is accepted, the
shower continues as normal. This ensures that the shower (with first emission
rejection) remains unitary.666As a matching algorithm to ensure the correct
matrix element for the first emission, it bears strong similarities to the
algorithm of Refs. Bengtsson:1986et ; Seymour:1994df . That algorithm works
within the assumption that one can find some constant $c$ such $cR_{\text{\sc
s}}>R$ over the full phase space, evading the need for an additive term. The
Born normalisation factor that multiplies the shower generation then needs to
be modified to read

$\tilde{B}_{\text{\sc s}}(\Phi_{\rm B})=B_{0}(\Phi_{\rm B})+V(\Phi_{\rm
B})+\int\min[R(\Phi),R_{\text{\sc s}}(\Phi)]\mathrm{d}\Phi_{\rm rad}\,.$ (14)

With this variant, the distribution of the hardest radiation is given by

$\mathrm{d}\sigma=\tilde{B}_{\text{\sc s}}{}(\Phi_{\rm
B})\,\tilde{S}(t_{\Phi},\Phi_{\rm B})\times\frac{\min(R(\Phi),R_{\text{\sc
s}}(\Phi))}{B_{0}(\Phi_{\rm B})}\mathrm{d}\Phi+\theta(R-R_{\text{\sc
s}})\left[R-R_{\text{\sc s}}\right]\mathrm{d}\Phi\,,$ (15)

where $\tilde{S}(t_{\Phi},\Phi_{\rm B})$ is the Sudakov form factor that is
effectively obtained as a result of the emission rejection procedure,

$\displaystyle\tilde{S}(t,\Phi_{\rm
B})=\exp\left[-\int_{t_{\Phi}>t}\frac{\min[R(\Phi),R_{\text{\sc
s}}(\Phi)]}{B_{0}(\Phi_{\rm B})}\mathrm{d}\Phi_{\rm rad}\right].$ (16)

##  5 Conclusion

In this article we have proposed an algorithm for NLOPS generators that, like
the MC@NLO approach, relies on the parton shower code for driving the
showering steps, but avoids the negative-weight events that complicate the
practical use of the MC@NLO approach. The approach can be viewed as a hybrid
version of the MC@NLO and KrkNLO methods. The method should be straightforward
to implement, requiring no more information than is already used for MC@NLO
codes, though it should be noted that it will probably require a closer
integration of the showering and fixed-order codes. Following the tradition of
associating an acronym to NLOPS methods, a suitable one for this method could
be MAcNLOPS, reflecting the combination of Multiplicative and Accumulative
(additive) steps.

## Acknowledgements

P.N. wishes to thank Stanislaw Jadach, Andrzej Siódmok, Silvia Ferrario
Ravasio, James Whitehead and Giulia Zanderighi for helpful conversations.
G.P.S. wishes to thank his collaborators on the PanScales project (Melissa van
Beekveld, Mrinal Dasgupta, Frédéric Dreyer, Basem El-Menoufi, Silvia Ferrario
Ravasio, Keith Hamilton, Alexander Karlberg, Rok Medves, Pier Monni, Ludovic
Scyboz, Alba Soto Ontoso, Grégory Soyez and Rob Verheyen) for discussions on
these and related topics. P.N. acknowledges support from Fondazione Cariplo
and Regione Lombardia, grant 2017-2070, and from INFN. G.P.S. acknowledges
support from a Royal Society Research Professorship
(RP$\backslash$R1$\backslash$180112), the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No. 788223, PanScales) and the Science and Technology
Facilities Council (STFC) under grant ST/T000864/1.

## References

  * (1) S. Frixione and B. R. Webber, _Matching NLO QCD computations and parton shower simulations_ , _JHEP_ 06 (2002) 029, [hep-ph/0204244].
  * (2) P. Nason, _A New method for combining NLO QCD with shower Monte Carlo algorithms_ , _JHEP_ 11 (2004) 040, [hep-ph/0409146].
  * (3) S. Plätzer, _Controlling inclusive cross sections in parton shower + matrix element merging_ , _JHEP_ 08 (2013) 114, [1211.5467].
  * (4) L. Lönnblad and S. Prestel, _Merging Multi-leg NLO Matrix Elements with Parton Showers_ , _JHEP_ 03 (2013) 166, [1211.7278].
  * (5) S. Jadach, W. Płaczek, S. Sapeta, A. Siódmok and M. Skrzypek, _Matching NLO QCD with parton shower in Monte Carlo scheme — the KrkNLO method_ , _JHEP_ 10 (2015) 052, [1503.06849].
  * (6) E. Boos et al., _Generic user process interface for event generators_ , in _Physics at TeV colliders. Proceedings, Euro Summer School, Les Houches, France, May 21-June 1, 2001_ , 2001.  hep-ph/0109068.
  * (7) M. Dasgupta, F. A. Dreyer, K. Hamilton, P. F. Monni and G. P. Salam, _Logarithmic accuracy of parton showers: a fixed-order study_ , _JHEP_ 09 (2018) 033, [1805.09327].
  * (8) M. Dasgupta, F. A. Dreyer, K. Hamilton, P. F. Monni, G. P. Salam and G. Soyez, _Parton showers beyond leading logarithmic accuracy_ , _Phys. Rev. Lett._ 125 (2020) 052002, [2002.11114].
  * (9) K. Hamilton, R. Medves, G. P. Salam, L. Scyboz and G. Soyez, _Colour and logarithmic accuracy in final-state parton showers_ , 2011.10054.
  * (10) A. Karlberg, G. P. Salam, L. Scyboz and R. Verheyen, _Spin correlations in final-state parton showers and jet observables_ , _Eur. Phys. J. C_ 81 (2021) 681, [2103.16526].
  * (11) J. R. Forshaw, J. Holguin and S. Plätzer, _Parton branching at amplitude level_ , _JHEP_ 08 (2019) 145, [1905.08686].
  * (12) J. R. Forshaw, J. Holguin and S. Plätzer, _Building a consistent parton shower_ , _JHEP_ 09 (2020) 014, [2003.06400].
  * (13) J. Holguin, J. R. Forshaw and S. Plätzer, _Improvements on dipole shower colour_ , _Eur. Phys. J. C_ 81 (2021) 364, [2011.15087].
  * (14) Z. Nagy and D. E. Soper, _Summations of large logarithms by parton showers_ , _Phys. Rev. D_ 104 (2021) 054049, [2011.04773].
  * (15) Z. Nagy and D. E. Soper, _Summations by parton showers of large logarithms in electron-positron annihilation_ , 2011.04777.
  * (16) 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].
  * (17) J. R. Andersen and A. Maier, _Unbiased Elimination of Negative Weights in Monte Carlo Samples_ , 2109.07851.
  * (18) K. Danziger, S. Höche and F. Siegert, _Reducing negative weights in Monte Carlo event generation with Sherpa_ , 2110.15211.
  * (19) P. Nason and B. Webber, _Next-to-Leading-Order Event Generators_ , _Ann. Rev. Nucl. Part. Sci._ 62 (2012) 187–213, [1202.1251].
  * (20) S. Jadach, G. Nail, W. Płaczek, S. Sapeta, A. Siódmok and M. Skrzypek, _KrkNLO in Herwig 7_ , _Acta Phys. Polon. B_ 48 (2017) 1121–1132.
  * (21) S. Jadach, G. Nail, W. Płaczek, S. Sapeta, A. Siodmok and M. Skrzypek, _Monte Carlo simulations of Higgs-boson production at the LHC with the KrkNLO method_ , _Eur. Phys. J. C_ 77 (2017) 164, [1607.06799].
  * (22) S. Jadach, W. Płaczek, S. Sapeta, A. K. Siodmok and M. Skrzypek, _New simpler methods of matching NLO corrections with parton shower Monte Carlo_ , _PoS_ LL2016 (2016) 020, [1607.00919].
  * (23) S. Jadach, W. Płaczek, S. Sapeta, A. Siodmok and M. Skrzypek, _Parton distribution functions in Monte Carlo factorisation scheme_ , _Eur. Phys. J. C_ 76 (2016) 649, [1606.00355].
  * (24) S. Jadach, W. Płaczek, S. Sapeta, A. Siódmok and M. Skrzypek, _NLO Corrections to Hard Process in Parton Shower MC - KrkNLO Method_ , _Acta Phys. Polon. B_ 46 (2015) 2089.
  * (25) S. Jadach, _On the universality of the KRK factorization scheme_ , _Acta Phys. Polon. B_ 51 (2020) 1363, [2004.04239].
  * (26) S. Alioli, P. Nason, C. Oleari and E. Re, _NLO Higgs boson production via gluon fusion matched with shower in POWHEG_ , _JHEP_ 04 (2009) 002, [0812.0578].
  * (27) M. Bengtsson and T. Sjostrand, _A Comparative Study of Coherent and Noncoherent Parton Shower Evolution_ , _Nucl. Phys. B_ 289 (1987) 810–846.
  * (28) M. H. Seymour, _Matrix element corrections to parton shower algorithms_ , _Comput. Phys. Commun._ 90 (1995) 95–101, [hep-ph/9410414].