text:
# Matching NLO QCD Computations and Parton Shower Simulations11footnote 1Work supported in part by the UK Particle Physics and Astronomy Research Council and by the EU Fourth Framework Programme ‘Training and Mobility of Researchers’, Network ‘Quantum Chromodynamics and the Deep Structure of Elementary Particles’, contract FMRX-CT98-0194 (DG 12 - MIHT).

## 1 Introduction

The reliable prediction of cross sections and final-state distributions for QCD processes is important not only as a test of QCD but also for the design of collider experiments and new particle searches. All systematic approaches to this problem are based on perturbation theory, usually truncated at next-to-leading order (NLO), or in a few cases next-next-to-leading order (NNLO). Such calculations yield the best available results for sufficiently inclusive observables, in kinematic regions where higher-order corrections are not enhanced. However, in many cases a more exclusive description of final states and/or a wider kinematic coverage are needed.

For the description of exclusive hadronic final states, perturbative calculations have to be combined with a model for the conversion of partonic final states into hadrons (hadronization). Existing hadronization models are in remarkably good agreement with a wide range of data, after tuning of model parameters . However, these models operate on partonic states with high multiplicity and low relative transverse momenta, which are obtained from a parton shower  or dipole cascade  approximation to QCD dynamics and not from fixed-order calculations.

To extend the coverage of perturbative calculations into regions where higher-order contributions are enhanced, one needs to identify and resum the enhanced terms to all orders. The resummed expression can then be combined with the fixed-order result, after subtracting resummed terms of the same order to avoid double-counting. This has been done explicitly for many important observables and processes, but a general procedure has not yet been developed. In principle, the parton shower provides such a procedure, if the problem of double-counting can be overcome. Combining NLO and parton-shower computations also offers the possibility of interfacing to hadronization models, thus overcoming the main deficiencies of existing calculations.

Methods for combining fixed-order and parton-shower or dipole-cascade computations in QCD have been discussed in refs. -. Many of those papers concern the rather different problem of combining leading-order (tree level) matrix elements for hard multi-jet production with parton showers or dipole cascades for jet fragmentation. The only fully-implemented schemes for matching NLO computations and parton showers that have been used to produce phenomenological results - have been based on the so-called phase-space slicing method for regularising the infrared divergences encountered at NLO. The slicing method is simple but has the disadvantage of introducing a slicing parameter, which should in principle be sent to zero at the end of the calculation but in practice cannot be chosen too small. This can introduce uncontrolled errors and/or large negative event weights. The slicing method has therefore tended to be replaced where possible in NLO calculations by subtraction methods, in which no approximations need to be made and event weights can be controlled by improving the choice of subtraction terms.

The subtraction method has been considered in refs. . The main concern of those papers was the correct matching of NLO matrix element calculations to the scale and scheme for next-to-leading logarithmic (NLL) parton showers. Up to now, no satisfactory algorithms for NLL showering have been developed. When such algorithms are available, the methods proposed there will make it possible for event generators to generate the NLO evolution of parton distribution functions (PDFs) and hard process cross sections, without the explicit input of NLO PDFs. Existing parton shower generators, however, correspond only to LL evolution, and therefore the matching criteria of refs.  operate beyond the level of precision that can be achieved at present.

In the present paper, we present a method for matching NLO matrix elements and existing parton shower generators based on the subtraction method for NLO calculations. First, in sect. 2, we outline the general objectives and limitations of our approach. The method is then illustrated in sect. 3 in the context of a simplified toy model. The case of QCD is dealt with in sects. 4 and 5. In sect. 4 results are presented for the process of W+W- production in hadronic collisions. In sect. 5 we discuss briefly the application of the method to other, more complex, processes. Concluding comments are given in sect. 6. All the technical details are collected in the appendices, where we also comment on the matching between NLO computations and parton shower approaches in the framework of the slicing method.

## 2 General framework

### 2.1 Objectives

It is undeniable that parton shower Monte Carlos (MC’s) provide more general and flexible tools compared to NLO computations. By running an MC we get events, which we believe to be a faithful description of the real events detected by real experiments; this is one of the reasons why the MC’s are so popular with experimentalists. On the other hand, NLO computations emphasize the role of (infrared-safe) observables; the presence of negative contributions, and the necessity of considering inclusive quantities in order to get rid of infinities, make it very problematic to even talk about events. Still, NLO computations have their virtues; they can handle hard emissions, and they can estimate total rates with better accuracy than MC’s.

In this paper, we aim to construct an MC that works like a regular MC, but on top of that knows how to treat hard emissions, and can compute rates to NLO accuracy. For brevity, we call this tool an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. Clearly, there is a lot of freedom in constructing an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. We shall follow a minimal approach. Namely, we want to keep the original MC as unmodified as possible; all major changes have to be carried out on the NLO codes. Amongst other things, this implies that our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ will not improve the logarithmic accuracy of the original MC, in those regions of the phase space where resummation is needed. As far as the NLO calculations are concerned, we shall not perform any approximations, neither at the level of matrix elements, nor at the level of phase space; thus, we shall adopt the subtraction method in order to cancel the infrared divergences that arise in the intermediate steps of the computation.

To be more definite, we require $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ to have the following characteristics:

* The output is a set of events, which are fully exclusive.
* Total rates are accurate to NLO.
* NLO results for all observables are recovered upon expansion of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results in $`\alpha _\mathrm{s}`$.
* Hard emissions are treated as in NLO computations.
* Soft/collinear emissions are treated as in MC.
* The matching between hard- and soft-emission regions is smooth.
* MC hadronization models are adopted.

In sect. 4, we shall construct an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ fulfilling these conditions. However, we find that the technical complications of QCD tend to hide the basic ideas upon which our approach is based. Thus, we prefer to start from an oversimplified (and unrealistic) case, where the structure of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ is apparent; this is done in sect. 3. The reason for doing this is that we shall see that the approach that leads to the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ in a simple case can actually be adopted in the QCD case as well, which will only be technically more involved.

### 2.2 Negative weights

We do not require the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ to be positive definite; events may have negative weights. However, it should be pointed out that these negative-weight events are of a completely different nature from the negative contributions that appear in NLO computations. In particular, the distributions of the positive- and negative-weight events are separately finite in our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ (the same is not true in an NLO computation). This implies the possibility of unweighted event generation,222In this paper, we adopt a slightly different notation from that commonly used by experimentalists. Namely, we call unweighted events any set of events whose weights are all identical, up to a sign. In other words, the weights of our unweighted events are of the form $`\pm w`$, where $`w`$ is a constant. as is customary in ordinary MC’s. The situation is actually completely analogous to that occurring in polarized collisions, where physical observables are obtained as the difference of two positive-definite and finite quantities. Furthermore, we shall see that in our case the number of negative-weight events is reasonably small, and thus the number of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ events necessary to get a smooth distribution is comparable to that in ordinary MC’s.

### 2.3 Double counting

When the MC generates events with real parton emission, it generates kinematical configurations that are also taken into account by the NLO computation: the possibility of having the same kinematical configuration from the MC and from the NLO may lead to double counting. We shall adopt the following quantitative definition of double counting:

An $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ is affected by double counting if its prediction for any observable, at the first order beyond the Born approximation333By the Born approximation we mean the lowest order at which the hard process contributes to any observable. in the expansion in the coupling constant, is not equal to the NLO prediction.

According to this definition, double counting may correspond to either an excess or a deficit in the prediction, at any point in phase space. This includes contributions from real emission and virtual corrections. Generally speaking, the MC fills the phase space efficiently only in the soft and collinear regions. In these regions, the leading behaviour of the kinematical distribution of one-particle emission is given by the Born cross section, times a kernel that describes the soft or collinear emission. This holds both for the MC and for the NLO. Thus, in the soft and collinear regions of the one-particle-emission phase space, the MC and NLO results coincide when only the leading terms are considered. Elsewhere, the MC is not reliable and the NLO provides the best estimate. The problem is to merge these two descriptions.

## 3 Toy model studies

In this section we study a toy model that allows a simple discussion of the key features of the NLO computation, of the MC approach, and of the matching between the two. We assume that a system can radiate massless particles (which we call photons), whose energy we denote by $`x`$, with $`0xx_s1`$, $`x_s`$ being the energy of the system before the radiation. After the radiation, the energy of the system is $`x_s^{}=x_sx`$. The system can undergo several further emissions; on the other hand, one photon cannot split further.

In a perturbative computation, the Born term corresponds to no emissions. The first non-trivial order in perturbation theory gets contribution from those diagrams with one and only one emission, with either virtual or real photons. We write the corresponding contributions to the cross section as follows:

$`\left({\displaystyle \frac{d\sigma }{dx}}\right)_\mathrm{B}`$ $`=`$ $`B\delta (x),`$ (1)
$`\left({\displaystyle \frac{d\sigma }{dx}}\right)_\mathrm{V}`$ $`=`$ $`a\left({\displaystyle \frac{B}{2ϵ}}+V\right)\delta (x),`$ (2)
$`\left({\displaystyle \frac{d\sigma }{dx}}\right)_\mathrm{R}`$ $`=`$ $`a{\displaystyle \frac{R(x)}{x}},`$ (3)
for the Born, virtual, and real contributions respectively; $`a`$ is the coupling constant, possibly times a colour factor; $`B`$ and $`V`$ are constant with respect to $`x`$, and

$$\underset{x0}{lim}R(x)=B.$$
(4)
The constant $`B`$ appears in eqs. (2) and (4) since we expect the residue of the leading singularity of the virtual and real contributions to be given by the Born term, times a suitable kernel. We take this kernel equal to 1, since this simplifies the computations, and it is not restrictive. Finally, $`ϵ`$ is the parameter entering dimensional regularization in $`42ϵ`$ dimensions.

The task of predicting an infrared-safe observable $`O`$ to NLO accuracy amounts to computing the following quantity

$$O=\underset{ϵ0}{lim}_0^1𝑑xx^{2ϵ}O(x)\left[\left(\frac{d\sigma }{dx}\right)_\mathrm{B}+\left(\frac{d\sigma }{dx}\right)_\mathrm{V}+\left(\frac{d\sigma }{dx}\right)_\mathrm{R}\right],$$
(5)
where $`O(x)`$ is the observable as a function of $`x`$, possibly times a set of $`\mathrm{\Theta }`$ functions defining a histogram bin. By “infrared-safe” here we simply mean that the integral in eq. (5) exists. The main technical problem in eq. (5) is due to the presence of the regularising parameter $`ϵ`$; in order to have an efficient numerical procedure, it is mandatory to extract the pole in $`ϵ`$ from the real contribution, in this way cancelling analytically the pole explicitly present in the virtual contribution.

### 3.1 NLO: slicing versus subtraction

Within the context of this toy model, the usual methods of cancelling the pole in $`ϵ`$ are as follows. In the slicing method, a small parameter $`\delta `$ is introduced into the real contribution in the following way:

$$O_\mathrm{R}=_0^\delta 𝑑xx^{2ϵ}O(x)\left(\frac{d\sigma }{dx}\right)_\mathrm{R}+_\delta ^1𝑑xx^{2ϵ}O(x)\left(\frac{d\sigma }{dx}\right)_\mathrm{R}.$$
(6)
In the first term on the r.h.s. of this equation we expand $`O(x)`$ and $`R(x)`$ in Taylor series around 0, and keep only the first term; the smaller $`\delta `$, the better the approximation. On the other hand, the second term in eq. (6) does not contain any singularity, and we can just set $`ϵ=0`$. We obtain

$`O_\mathrm{R}`$ $`=`$ $`aBO(0){\displaystyle _0^\delta }𝑑x{\displaystyle \frac{x^{2ϵ}}{x}}+{\displaystyle _\delta ^1}𝑑xO(x)\left({\displaystyle \frac{d\sigma }{dx}}\right)_\mathrm{R}+𝒪(\delta )`$ (7)
$`=`$ $`a\left({\displaystyle \frac{1}{2ϵ}}+\mathrm{log}\delta \right)BO(0)+a{\displaystyle _\delta ^1}𝑑x{\displaystyle \frac{O(x)R(x)}{x}}+𝒪(\delta ,ϵ).`$ (8)
Using this result in eq. (5), we get the NLO prediction for $`O`$ as given in the slicing method:

$$O_{\mathrm{slice}}=BO(0)+a\left[\left(B\mathrm{log}\delta +V\right)O(0)+_\delta ^1𝑑x\frac{O(x)R(x)}{x}\right]+𝒪(\delta ).$$
(9)
The terms $`𝒪(\delta )`$ cannot be computed, and are neglected; in general, one must check that these neglected terms are numerically small, by plotting $`O_{\mathrm{slice}}`$ versus $`\delta `$ in a suitable range in $`\delta `$.

In the subtraction method, no approximation is performed. One writes the real contribution as follows:

$$O_\mathrm{R}=aBO(0)_0^1𝑑x\frac{x^{2ϵ}}{x}+a_0^1𝑑x\frac{O(x)R(x)BO(0)}{x^{1+2ϵ}}.$$
(10)
The second term on the r.h.s. does not contain singularities, and we can set $`ϵ=0`$:

$$O_\mathrm{R}=a\frac{B}{2ϵ}O(0)+a_0^1𝑑x\frac{O(x)R(x)BO(0)}{x}.$$
(11)
Therefore, the NLO prediction as given in the subtraction method is:

$$O_{\mathrm{sub}}=BO(0)+a\left[VO(0)+_0^1𝑑x\frac{O(x)R(x)BO(0)}{x}\right].$$
(12)
We rewrite this in a slightly different form, which will be more suited to discussion of the matching with the MC:

$$O_{\mathrm{sub}}=_0^1𝑑x\left[O(x)\frac{aR(x)}{x}+O(0)\left(B+aV\frac{aB}{x}\right)\right].$$
(13)
### 3.2 Toy Monte Carlo

In a treatment based on Monte Carlo methods, the system can undergo an arbitrary number of emissions (branchings), with probability controlled by the Sudakov form factor, defined for our toy model as follows:

$$\mathrm{\Delta }(x_1,x_2)=\mathrm{exp}\left[a_{x_1}^{x_2}𝑑z\frac{Q(z)}{z}\right],$$
(14)
where $`Q(z)`$ is a monotonic function with the following general properties:

$$0Q(z)1,\underset{z0}{lim}Q(z)=1,\underset{z1}{lim}Q(z)=0.$$
(15)
Note that the MC is also well defined without the last condition in eq. (15). We discuss later the specific functional form assumed for $`Q`$. If $`x_s`$ is the energy of the system before the first branching occurs, then $`\mathrm{\Delta }(x,x_s)`$ is the probability that no photon be emitted with energy $`z`$ such that $`xzx_s`$. At fixed $`x`$, this probability tends to 1 when $`a0`$; that is, no emission is possible in the zero-coupling limit. On the other hand, when $`a\mathrm{}`$, we have $`\mathrm{\Delta }(x,x_s)0`$, which means that it is impossible that the system does not emit. The toy MC code we have written implements ordered emissions in the photon energies444This may violate energy-momentum conservation, as is usual in a leading logarithm implementation. and proceeds through the following steps:

* Define $`x_0`$ as the lower bound on the photon energy after a branching; thus $`x_0`$ plays the same role as the cutoffs used to define the Sudakov form factor in QCD MC’s.
* Define $`x_\mathrm{M}`$ to be the maximum energy available to the photon; before the first branching, we have $`x_\mathrm{M}=x_s`$, where $`x_s`$ is the energy of the system.
* Pick a random number $`R_1[0,1]`$, and solve for $`z`$ the equation

$$\mathrm{exp}\left[a_z^{x_\mathrm{M}}\frac{dz^{}}{z^{}}\right]\left(\frac{z}{x_\mathrm{M}}\right)^a=R_1z=x_\mathrm{M}R_1^{1/a}.$$
(16)
* Pick a random number $`R_2[0,1]`$. If $`Q(z)<R_2`$, let $`x_\mathrm{M}=z`$, reject the branching, and return to step 2).
* If $`z<x_0`$, reject the branching and exit. If $`z>x_0`$, a photon with energy $`z`$ has been emitted. We let $`x_\mathrm{M}=z`$, and iterate the procedure going back to 2).

This procedure ensures that the photon emissions have the distribution prescribed by the Sudakov form factor in eq. (14). As can be seen from that equation, the probability of having soft emission (i.e., with $`z<x_0`$), and thus of terminating the shower, is $`\mathrm{\Delta }(x_0,x_s)`$, which tends to one when $`a0`$ (weak coupling), and when $`x_sx_0`$ (no energy left for non-soft branching). In this toy model, the hardest emission is always the first one.

### 3.3 Matching NLO and MC

In order to introduce some useful notation, we start from the LO perturbative result, which does not pose any problems, since it is in fact what is implemented in ordinary MC’s. At the LO, the system has not yet lost any energy to perturbative photon emission. Thus, the maximum energy available to the photon in the first branching is 1. For each event, we set $`x_\mathrm{M}=1`$, and proceed as described in sect. 3.2. The total exclusive rate will have to be equal to $`B`$, according to eq. (1); thus, if we generate $`N`$ events, each event has a weight equal to $`B/N`$.

We can formally read this procedure from eq. (13), by setting $`a=0`$ and performing the formal substitution

$$O(x)I_{\mathrm{MC}}(O,x_\mathrm{M}(x)),$$
(17)
where $`I_{\mathrm{MC}}(O,x_\mathrm{M})`$ stands for interface-to-MC, and indicates symbolically the distribution in the observable $`O`$ as obtained by running the MC starting from a given $`x_\mathrm{M}`$, and giving each event the weight $`1/N`$; its precise definition is given in app. B.1. For $`x=0`$ we need $`x_\mathrm{M}=1`$; we do not need to specify $`x_\mathrm{M}(x)`$ for other values of $`x`$ at present. Interfacing LO with MC, we thus obtain

$$\left(\frac{d\sigma }{dO}\right)_{\mathrm{MC}\mathrm{@}\mathrm{LO}}=BI_{\mathrm{MC}}(O,1).$$
(18)
This equation reminds us of the total rate ($`B`$), and the maximum energy available for photons at the first branching ($`x_\mathrm{M}=1`$).

#### 3.3.1 Naive subtraction

The first attempt at extending eq. (18) to NLO amounts simply to substituting eq. (17) into eq. (13). We get

$$\left(\frac{d\sigma }{dO}\right)_{\mathrm{naive}}=_0^1𝑑x\left[I_{\mathrm{MC}}(O,x_\mathrm{M}(x))\frac{aR(x)}{x}+I_{\mathrm{MC}}(O,1)\left(B+aV\frac{aB}{x}\right)\right].$$
(19)
This equation suggests the following procedure:

* Pick at random $`0x1`$.
* Generate an MC event with $`x_\mathrm{M}(x)`$ as maximum energy available to the photon in the first branching; attach to this event the weight $`w_{\mathrm{EV}}=aR(x)/x`$.
* Generate another MC event (a “counter-event”) with $`x_\mathrm{M}=1`$; attach to this event the weight $`w_{\mathrm{CT}}=B+aVaB/x`$.
* Repeat the first three steps $`N`$ times, and normalize with $`1/N`$.

Unfortunately, this procedure is bound to fail, since the weights $`w_{\mathrm{EV}}`$ and $`w_{\mathrm{CT}}`$ diverge as $`x0`$. A small cutoff $`\rho `$ can then be introduced; as $`x\rho `$, we would have $`w_{\mathrm{EV}}aB/\rho `$ and $`w_{\mathrm{CT}}aB/\rho `$. Both these weights are very large in absolute value; this means that the eventual unweighting procedure, performed by the MC, will be highly inefficient. Furthermore, as we shall see in sect. 3.4, eq. (19) has problems of double counting, defined according to the criterion in sect. 2.3.

#### 3.3.2 Modified subtraction

To avoid the problems outlined above, we define the modified subtraction method by

$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{msub}}`$ $`=`$ $`{\displaystyle _0^1}dx[I_{\mathrm{MC}}(O,x_\mathrm{M}(x)){\displaystyle \frac{a[R(x)BQ(x)]}{x}}`$ (20)
$`+I_{\mathrm{MC}}(O,1)(B+aV+{\displaystyle \frac{aB[Q(x)1]}{x}})].`$
This amounts to saying that, in eq. (19), we subtract and add the quantity

$$I_{\mathrm{MC}}(O,x_\mathrm{M})\frac{aBQ(x)}{x},$$
(21)
using $`x_\mathrm{M}=x_\mathrm{M}(x)`$ in the first and $`x_\mathrm{M}=1`$ in the second term introduced in this way. We point out that $`Q`$ is the same function that appears in the definition of the Sudakov form factor, eq. (14). By virtue of the second condition in eq. (15), the weights $`w_{\mathrm{EV}}`$ and $`w_{\mathrm{CT}}`$ attached to events and counter-events are now separately convergent as $`x0`$.

The two terms involving $`Q(x)`$, that we inserted in eq. (20), are not identical; therefore, the procedure which leads from eq. (19) to eq. (20) is not a subtraction in the usual sense of an NLO computation. Still, these two terms do not contribute to the observable $`O`$ at $`𝒪(a)`$, because they are compensated by analogous terms due to the parton shower $`BI_{\mathrm{MC}}(O,1)`$; this will be shown with explicit examples in sect. 3.4, and in general in app. B.3. We point out that the function $`aQ(x)/x`$ that appears in eq. (21) does not depend upon the MC cutoff $`x_0`$. It follows that the dependence upon $`x_0`$ of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, eq. (20), is identical to that of the ordinary MC, being entirely due to $`I_{\mathrm{MC}}`$, and is therefore power-like.

It should be clear from the description of the procedure given above that the dependence upon $`O`$ in eqs. (19) and (20) is purely formal, serving only to allow us to write such equations in a compact form. We stress that, in the actual implementation of all the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$’s described in this paper, the computation of observables proceeds exactly as in standard MC’s: namely, events are generated, showered, and hadronized without reference to any specific observable.

Equation (20) is our master equation for the toy model. We shall show in what follows that the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ defined in this way meets the requirements given in sect. 2.1. In sect. 4, we shall use eq. (20) to construct a QCD $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ by analogy with the toy model.

### 3.4 Observables

#### 3.4.1 Exclusive observable

As an example of an “exclusive” observable in the toy model, we shall consider the quantity

$$y=\mathrm{max}(x_1,\mathrm{},x_n),$$
(22)
where $`n`$ denotes the number of emissions in one given event (in the case of no emission, the event is not used to fill the histograms). This observable is the energy of the hardest photon of the event. It is particularly convenient because in the MC approach it coincides with the energy of the first photon emitted in the shower; when $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ is considered, $`y`$ is the larger of the energy of the photon emitted perturbatively, and that of the first photon in the shower (if any). This fact allows us to obtain analytical results for MC, NLO and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$.

The NLO prediction for this observable is (here and in the following, we consider only $`y>x_0`$ in order to simplify the analytical results)

$$\left(\frac{d\sigma }{dy}\right)_{\mathrm{NLO}}=a\frac{R(y)}{y},$$
(23)
while the MC result, normalized to the Born cross section, is

$$\left(\frac{d\sigma }{dy}\right)_{\mathrm{MC}}=aB\frac{Q(y)}{y}\mathrm{\Delta }(y,1),$$
(24)
which is just the differential cross section at $`x=y`$, computed using eq. (14), times the probability that there is no harder emission at $`x>y`$.

In order to study the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results, we have to consider the quantities $`I_{\mathrm{MC}}(y,1)`$ and $`I_{\mathrm{MC}}(y,x_\mathrm{M}(x))`$ that will be substituted into eq. (19) or (20). In the case of $`I_{\mathrm{MC}}(y,1)`$, the hardest photon comes from the MC and we have

$$I_{\mathrm{MC}}(y,1)=a\frac{Q(y)}{y}\mathrm{\Delta }(y,1).$$
(25)
In the case of $`I_{\mathrm{MC}}(y,x_\mathrm{M}(x))`$, the hardest photon may come from the NLO contribution or the MC, and so we find

$`I_{\mathrm{MC}}(y,x_\mathrm{M}(x))`$ $`=`$ $`\mathrm{\Delta }(\mathrm{min}\{x,x_\mathrm{M}(x)\},x_\mathrm{M}(x))\delta (yx)`$ (26)
$`+a{\displaystyle \frac{Q(y)}{y}}\mathrm{\Delta }(y,x_\mathrm{M}(x))\mathrm{\Theta }(yx)\mathrm{\Theta }(x_\mathrm{M}(x)y).`$
As a shorthand notation, in the following we shall denote the contributions to $`d\sigma /dy`$ from $`I_{\mathrm{MC}}(y,1)`$ and $`I_{\mathrm{MC}}(y,x_\mathrm{M}(x))`$ by $`d\sigma _𝕊/dy`$ and $`d\sigma _{}/dy`$ respectively, since they correspond to a $`𝕊`$tandard MC evolution (i.e., with no prior NLO real emission), and to a $``$ard MC evolution (whose initial condition includes one NLO real emission). This notation will be extensively used from sect. 4.1 on. In the case of the naive subtraction method, eq. (19), these contributions are separately divergent, but their sum is not. We find

$$\left(\frac{d\sigma _𝕊}{dy}\right)_{\mathrm{naive}}=a\frac{Q(y)}{y}\mathrm{\Delta }(y,1)\left[B+aVaB_0^1\frac{dx}{x}\right].$$
(27)
The results for $`d\sigma _{}/dy`$ depend somewhat on the form of the function $`x_\mathrm{M}(x)`$ introduced in eq. (17). We define $`x_e`$ to be the solution of $`x=x_\mathrm{M}(x)`$, i.e., $`x_e=x_\mathrm{M}(x_e)`$, and we assume this solution to be unique. We further assume $`x_\mathrm{M}(x)`$ to be a monotonically decreasing function; thus, $`x<x_\mathrm{M}(x)`$ for $`x<x_e`$. Then:

$``$ for $`y<x_e:`$

$$\left(\frac{d\sigma _{}}{dy}\right)_{\mathrm{naive}}=a\frac{R(y)}{y}\mathrm{\Delta }(y,x_\mathrm{M}(y))+a^2\frac{Q(y)}{y}_0^y\frac{dx}{x}R(x)\mathrm{\Delta }(y,x_\mathrm{M}(x));$$
(28)
$``$ for $`y>x_e:`$

$$\left(\frac{d\sigma _{}}{dy}\right)_{\mathrm{naive}}=a\frac{R(y)}{y}+a^2\frac{Q(y)}{y}_0^{x_\mathrm{M}^1(y)}\frac{dx}{x}R(x)\mathrm{\Delta }(y,x_\mathrm{M}(x)).$$
(29)
Expanding in $`a`$, we find

$$\left(\frac{d\sigma }{dy}\right)_{\mathrm{naive}}\left(\frac{d\sigma _𝕊}{dy}\right)_{\mathrm{naive}}+\left(\frac{d\sigma _{}}{dy}\right)_{\mathrm{naive}}=\frac{a}{y}[R(y)+BQ(y)]+𝒪(a^2).$$
(30)
Therefore the naive subtraction method suffers from double counting, according to the definition in sect. 2.3, especially at small $`y`$, where $`Q(y)1`$.

For the modified subtraction method, eq. (20), we find instead that both contributions are finite:

$$\left(\frac{d\sigma _𝕊}{dy}\right)_{\mathrm{msub}}=a\frac{Q(y)}{y}\mathrm{\Delta }(y,1)\left[B+aV+aB_0^1𝑑x\frac{Q(x)1}{x}\right],$$
(31)
while

$``$ for $`y<x_e:`$

$`\left({\displaystyle \frac{d\sigma _{}}{dy}}\right)_{\mathrm{msub}}`$ $`=`$ $`{\displaystyle \frac{a}{y}}[R(y)BQ(y)]\mathrm{\Delta }(y,x_\mathrm{M}(y))`$ (32)
$`+`$ $`a^2{\displaystyle \frac{Q(y)}{y}}{\displaystyle _0^y}{\displaystyle \frac{dx}{x}}[R(x)BQ(x)]\mathrm{\Delta }(y,x_\mathrm{M}(x));`$
$``$ for $`y>x_e:`$

$$\left(\frac{d\sigma _{}}{dy}\right)_{\mathrm{msub}}=\frac{a}{y}[R(y)BQ(y)]+a^2\frac{Q(y)}{y}_0^{x_\mathrm{M}^1(y)}\frac{dx}{x}[R(x)BQ(x)]\mathrm{\Delta }(y,x_\mathrm{M}(x)).$$
(33)
Expanding in $`a`$, we now get

$$\left(\frac{d\sigma }{dy}\right)_{\mathrm{msub}}=\frac{aR(y)}{y}+𝒪(a^2),$$
(34)
in agreement with the NLO prediction for all values of $`y`$. Furthermore, for $`y0`$ we have the expected resummation of leading logarithms,

$$\left(\frac{d\sigma }{dy}\right)_{\mathrm{msub}}\left(\frac{aB}{y}+𝒪(a^2)\right)\mathrm{\Delta }(y,1)+\mathrm{PST}\mathrm{for}y0,$$
(35)
where PST stands for power-suppressed terms, i.e., terms that are not singular for $`y0`$. We stress that the behaviour of $`d\sigma /dy`$ in eq. (35) is dictated by eq. (31), i.e., the contribution of $``$ events is power suppressed for $`y0`$, since $`(R(y)BQ(y))/y=𝒪(1)`$ in this limit. We also notice that the functional form in eq. (35) is the same as that in eq. (24). Thus, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ resums large logarithms of $`y`$ in the same way as does the MC. A general derivation of this result for a wide class of observables is given in app. B.4.

#### 3.4.2 Inclusive observable

As an “inclusive” observable, we shall consider the fully inclusive distribution of the photon energies; thus, each photon emitted, either at the NLO level or by the MC shower, will contribute to the physical observable, which we shall denote by $`z`$. The NLO and MC results (for $`z>x_0`$) are:

$`\left({\displaystyle \frac{d\sigma }{dz}}\right)_{\mathrm{NLO}}`$ $`=`$ $`a{\displaystyle \frac{R(z)}{z}},`$ (36)
$`\left({\displaystyle \frac{d\sigma }{dz}}\right)_{\mathrm{MC}}`$ $`=`$ $`aB{\displaystyle \frac{Q(z)}{z}}.`$ (37)
Eq. (36) coincides with eq. (23), which is a trivial consequence of the fact that only one emission occurs at NLO. On the other hand, the MC result in eq. (37) differs from the exclusive case, eq. (24), since all emissions now contribute, not only the first. As a consequence, eq. (37) has a logarithmically-divergent integral at $`z0`$, whereas this does not happen in eq. (24).

In order to obtain the predictions of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, we need to find the analogues of eqs. (25) and (26), in the case of the variable $`z`$. We obtain

$`I_{\mathrm{MC}}(z,1)`$ $`=`$ $`a{\displaystyle \frac{Q(z)}{z}},`$ (38)
$`I_{\mathrm{MC}}(z,x_\mathrm{M}(x))`$ $`=`$ $`\delta (zx)+a{\displaystyle \frac{Q(z)}{z}}\mathrm{\Theta }(x_\mathrm{M}(x)z).`$ (39)
The term with the $`\delta `$-function in eq. (39) accounts for real emission at the NLO level, which always contributes to the inclusive spectrum when it occurs; the other term is the contribution of the MC shower.

To compute the $`z`$ distribution for the naive subtraction method, we use eq. (38) and eq. (39) in eq. (19), to obtain

$$\left(\frac{d\sigma }{dz}\right)_{\mathrm{naive}}=\frac{a}{z}\left[R(z)+Q(z)\left(\sigma _{tot}a_{x_\mathrm{M}^1(z)}^1\frac{dx}{x}R(x)\right)\right],$$
(40)
where $`\sigma _{tot}`$ is the NLO total cross section

$$\sigma _{tot}=B+aV+a_0^1\frac{dx}{x}[R(x)B].$$
(41)
As before, we may expand in $`a`$ to find

$$\left(\frac{d\sigma }{dz}\right)_{\mathrm{naive}}=\frac{a}{z}[R(z)+BQ(z)]+𝒪(a^2),$$
(42)
which shows again that the naive subtraction method suffers from double counting.

For the modified subtraction method we use eq. (38) and eq. (39) in eq. (20) to find

$`\left({\displaystyle \frac{d\sigma }{dz}}\right)_{\mathrm{msub}}`$ $`=`$ $`{\displaystyle \frac{a}{z}}[R(z)BQ(z)`$ (43)
$`+Q(z)(\sigma _{tot}a{\displaystyle _{x_\mathrm{M}^1(z)}^1}dx{\displaystyle \frac{R(x)BQ(x)}{x}})].`$
Expanding in $`a`$ we now find

$$\left(\frac{d\sigma }{dz}\right)_{\mathrm{msub}}=\frac{aR(z)}{z}+𝒪(a^2);$$
(44)
thus, as in the case of the observable $`y`$, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ prediction at $`𝒪(a)`$ coincides with the NLO result for all $`z`$. The small-$`z`$ limit of eq. (43) is

$$\underset{z0}{lim}z\left(\frac{d\sigma }{dz}\right)_{\mathrm{msub}}=a\sigma _{tot}.$$
(45)
From eq. (37) we obtain instead

$$\underset{z0}{lim}z\left(\frac{d\sigma }{dz}\right)_{\mathrm{MC}}=aB.$$
(46)
Therefore, at $`z0`$ the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ result is equal to that of the standard MC, times $`\sigma _{tot}/B`$, which is the K-factor for the total rate.

### 3.5 Results

In this section we present toy model results for NLO, MC, and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. The following parameters have been used in eqs. (1), (2) and (3):

$$a=0.3,B=2,V=1.$$
(47)
For the real emission distribution in eq. (3) we have taken

$$R(x)=B+x(1+x/2+20x^2),$$
(48)
so that the total cross section is

$$\sigma _{tot}=B+a\left(V+\frac{95}{12}\right).$$
(49)
Note that this implies a large K-factor, $`\sigma _{tot}/B=2.34`$. Finally, for the MC cutoff parameter we used

$$x_0=0.02,$$
(50)
and for the maximum energy function in eq. (17) we chose the form $`x_\mathrm{M}(x)=1x`$. In the left panel of fig. 1 we show the NLO result obtained with eq. (13) for the observable $`y`$ defined in eq. (22) (solid histogram), compared to the result obtained by running the MC (dashed histogram), with $`Q(x)1`$; the integral of the latter result has been normalized to the NLO rate (i.e., each event is given a weight equal to $`\sigma _{tot}/N`$). As a cross check for the numerical implementation against analytical results, we also plot as solid lines the functions given in eqs. (23) and (24) (for the latter, $`B`$ has been substituted with $`\sigma _{tot}`$, in order to be consistent with the normalization of the dashed histogram), which in fact appear to be in perfect agreement with the corresponding histograms. The shapes of the NLO and MC results are completely different. In fact, the function $`R(x)`$ has been chosen to introduce a clear difference between the NLO and MC results, for a better understanding of the effect of the matching between NLO and MC.

To study the matching of NLO and MC, we consider an MC with a dead zone. In our toy model, this corresponds to an MC that cannot emit photons with energy $`x>x_{dead}`$ (in the numerical simulations we set $`x_{dead}=0.6`$). We can easily get a toy MC with a dead zone by suitably choosing the function $`Q`$ in eq. (14). We first introduce a smoothing function of the form

$$G(x)=\frac{c^2(1x)^{2\beta }}{x^{2\alpha }+c^2(1x)^{2\beta }},$$
(51)
where $`\alpha `$, $`\beta `$ and $`c`$ are free parameters, and then consider three different implementations of the dead zone:

$`1.Q(x)=\mathrm{\Theta }(x_{dead}x);`$ (52)
$`2.Q(x)=\mathrm{\Theta }(x_{dead}x)G(x/x_{dead}),\mathrm{with}\alpha =1,\beta =1,c=1;`$ (53)
$`3.Q(x)=\mathrm{\Theta }(x_{dead}x)G(x/x_{dead}),\mathrm{with}\alpha =2,\beta =1,c=8.`$ (54)
The results for the pure MC evolution are presented in the right panel of fig. 1; the dependence upon $`Q`$ is actually fairly big, since we made extreme choices in eqs. (52)–(54). It is thus striking that the corresponding results obtained with the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, shown in the left panel of fig. 2, do not display big differences from each other. All three $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ curves are quite close to the pure NLO result (and coincide with it for $`y>x_{dead}`$), and differ sizably from it only in the region of small $`x`$, where the effect of multiple emission is expected to be important.555Since $`y`$ is an exclusive variable, and the total rate is that of the NLO result, the differences between the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ and NLO curves at $`y>x_0`$ are in fact compensated by the region $`y<x_0`$ in the total rate. In this region, however, the three $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ curves are close to each other, notwithstanding the differences shown in the right panel of fig. 1. We are here in a region where the leading logarithmic behaviour of the MC is highly dominant, and thus the impact of $`Q`$ (which is subleading in logarithmic accuracy) is moderate. We see that the curve corresponding to eq. (52) displays a less regular behaviour than the others at $`yx_{dead}`$; however, the effect is very moderate, being at worst of $`𝒪(a^2)`$.

The results for the inclusive $`z`$ distribution are presented in the right panel of fig. 2. Again, we can see that the three $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ curves are remarkably close to each other, in spite of the differences in the implementations of the dead zone. The ratio of the cross section in the leftmost bin of the solid histogram over the NLO result for the same quantity is 2.23. This value is fairly close to the K-factor $`\sigma _{tot}/B=2.34`$, that we should find according to eqs. (45) and (46); the small difference is due to fact that this bin is not quite at the limit $`z0`$.

We can thus conclude that $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ reduces the numerical dependence upon subleading logarithms, which is seen when running the MC alone; this is because those regions where the subleading terms are becoming numerically important get sizable contributions from the hard NLO emissions in the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, which are not present in the MC, and which eventually become dominant. On the other hand, in the soft-dominated regions, the shapes of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ curves just reproduce those of the MC curves, as we expect.

## 4 The QCD case: W+W- production

We now face the problem of constructing an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ in QCD. In order to simplify as much as possible all the technical details related to the colour structure of the hard process, we consider here the production of W+W- pairs in hadronic collisions. This process also has a fairly simple structure of soft and collinear singularities. However, we argue that the method we propose is actually valid for any kind of hard production process.

We shall use HERWIG  as our reference parton shower MC. Other choices are of course possible, and our approach is not tailored to work with HERWIG, although some of the technical details, reported in the appendices, do depend on the specific implementation of the hard process and parton showering in that program.

In order to define our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ for W pair production, we rely on the results previously presented for the toy model. First we emphasize the main points of the toy $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. Then we make a parallel between the contributions to the toy model cross section, and the contributions to the W+W- cross section. Using these results, we proceed to construct the QCD $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, by analogy with the toy $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$.

### 4.1 Toy model revisited

Let us consider again our master equation for the toy model $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, eq. (20). The $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ approach consists in running the ordinary MC, with initial conditions and weights determined by the NLO calculation. Depending upon the initial conditions, we distinguish two classes of events, which correspond to $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$ and $`I_{\mathrm{MC}}(O,1)`$. We denote the class corresponding to $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$ by $``$, which stands for “hard”, and reminds us that one emission already occurred, due to NLO, prior to the MC run. The class corresponding to $`I_{\mathrm{MC}}(O,1)`$ is denoted by $`𝕊`$, which stands for “standard”; this reminds us that this initial condition is the same as for the ordinary MC runs.

It is important to realize that for $`𝕊`$ events the finite part of the virtual contribution, $`aV`$, the real counterterm, $`aB/x`$, and the term $`aBQ(x)/x`$, serve exclusively to fix the normalization of this class of events, relative to that of the class $``$. In other words, it is irrelevant that the latter two terms explicitly depend upon $`x`$: in the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, the spectrum of the observable $`O`$ resulting from $`𝕊`$ events is entirely due to MC evolution. In contrast, the NLO hard emissions do contribute to the kinematics of $``$ events, and thus help to determine $`O`$.

In order to obtain in QCD a structure analogous to that outlined above, we therefore have to collect all the terms of the NLO cross section that do not contribute to hard emissions, and use them to compute the normalization of class $`𝕊`$ events; these events result from MC evolutions whose initial conditions are identical to those of an ordinary MC. On the other hand, the kinematics of NLO hard emissions determine the initial conditions for MC evolutions for $``$ events.

Finally, we have already seen in the case of the observables $`y`$ and $`z`$ (see also app. B.3 for a more formal proof) that at $`𝒪(a)`$ the dependence upon $`Q(x)`$ disappears. Roughly speaking, this happens because the contribution due to the Born term $`B`$, in class $`𝕊`$ events, cancels at $`𝒪(a)`$ with the contribution of the terms $`\pm aBQ(x)/x`$. We observe that, in the toy model, $`aBQ(x)/x`$ is the $`𝒪(a)`$ term in the expansion of the MC result (see eqs. (24) and (37)). We therefore guess that the role of $`aBQ(x)/x`$ in a QCD $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ is played by the $`𝒪(\alpha _\mathrm{s})`$ term in the expansion of the result of the ordinary MC, when only the shower is taken into account (i.e., the hadronization step is skipped). We shall eventually prove (see app. B.3) that this is in fact the case.

### 4.2 NLO cross section

We now deal explicitly with W+W- production in hadronic collisions, and we start by writing the NLO cross section. In this paper, we use the results of ref. , which we rewrite in a form more suited to our present goal; it should be clear that any W+W- cross section (or at least, one based upon the subtraction method) can be cast in the same form. Owing to the factorization theorem, we have

$$d\sigma (P_1,P_2)=\underset{ab}{}𝑑x_1𝑑x_2f_a^{(H_1)}(x_1)f_b^{(H_2)}(x_2)𝑑\widehat{\sigma }_{ab}(x_1P_1,x_2P_2),$$
(55)
for the collision of hadrons $`H_1`$ and $`H_2`$ with momenta $`P_1`$ and $`P_2`$ respectively. The sum runs over those parton flavours that give partonic subprocesses contributing to W pair production at NLO. We symbolically denote these processes as follows:

$`q\overline{q}`$ $``$ $`W^+W^{},`$ (56)
$`q\overline{q}`$ $``$ $`W^+W^{}g,`$ (57)
$`qg`$ $``$ $`W^+W^{}q,`$ (58)
where $`q=u,d,\mathrm{}`$ and eq. (58) includes the corresponding antiquark processes. Equation (56) contributes to Born and virtual emission terms, whereas eqs. (57) and (58) contribute to real emission terms.

The information on the hard process in eq. (55) is contained in the subtracted partonic cross sections, $`d\widehat{\sigma }_{ab}`$ (the hat means that counterterms are implicit in the expression). We write these cross sections as follows (eq. (3.45) of ref. ):

$$d\widehat{\sigma }_{ab}=d\sigma _{ab}^{(b)}+d\sigma _{ab}^{(sv)}+d\widehat{\sigma }_{ab}^{(c+)}+d\widehat{\sigma }_{ab}^{(c)}+d\widehat{\sigma }_{ab}^{(f)},$$
(59)
where:

* $`d\sigma _{ab}^{(b)}`$ is the Born contribution. It has $`22`$ kinematics, and corresponds to eq. (1) of the toy model.
* $`d\sigma _{ab}^{(sv)}`$ is the contribution of the virtual diagrams, and of all those terms, resulting from various subtractions, that have $`22`$ kinematics and are of $`𝒪(\alpha _\mathrm{s})`$. This term corresponds to the finite part of eq. (2) of the toy model.
* $`d\widehat{\sigma }_{ab}^{(c\pm )}`$ are the contributions of the finite remainders resulting from the subtraction of the initial-state collinear divergences ($`+`$ and $``$ indicate radiation from partons $`a`$ and $`b`$ respectively). These terms do not have an analogue in the toy model, since in the toy model there are no parton densities. They have $`23`$ kinematics, with the light outgoing parton always collinear to the beam-line; we refer to these configurations as quasi-$`22`$ kinematics. We stress that $`d\widehat{\sigma }_{ab}^{(c\pm )}`$ are subtracted quantities, and that the counterterms have $`22`$ kinematics.
* $`d\widehat{\sigma }_{ab}^{(f)}`$ is the contribution of the real emission diagrams, with its divergences suitably subtracted; it has $`23`$ kinematics (the counterterms have $`22`$ or quasi-$`22`$ kinematics). The unsubtracted part of this term corresponds to eq. (3) of the toy model.

We point out that all the terms on the r.h.s. of eq. (59) are separately finite. Also, they do not contain any large logarithms, and so no cancellation of large numbers occurs among them. This property is a feature of the subtraction method, and does not hold in the slicing method.

Items a)-d) above achieve the identification between terms contributing to the QCD cross section and to the toy model cross section. We now go into more detail, by writing the quantities appearing in eq. (59) more explicitly .

We denote by $`d\varphi _2(s)`$ and $`d\varphi _3(s)`$ the phase spaces of the W+W- pair and of the W+W- pair plus a light parton, respectively; they depend upon the kinematical variables of these particles, and upon $`s=(p_1+p_2)^2`$, the c.m. energy squared of the partons that initiate the hard scattering, whose momenta we denote by $`p_1=x_1P_1`$ and $`p_2=x_2P_2`$. The Born term is

$$d\sigma _{ab}^{(b)}(p_1,p_2)=_{ab}^{(b)}(p_1,p_2)d\varphi _2(s),$$
(60)
where $`_{ab}^{(b)}`$ is the Born matrix element squared, summed over spin and colour degrees of freedom (the average factor for initial-state ones is understood), times the flux factor $`1/(2s)`$. The soft-virtual contribution is

$$d\sigma _{ab}^{(sv)}(p_1,p_2)=\left[𝒬_{ab}_{ab}^{(b)}(p_1,p_2)+_{ab}^{(v,reg)}(p_1,p_2)\right]d\varphi _2(s),$$
(61)
where $`𝒬_{ab}`$ is a constant with respect to the kinematics of the partonic process, and $`_{ab}^{(v,reg)}`$ is the finite part of the virtual contribution. As with any other finite part resulting from the cancellation of singularities, its definition is to a certain extent arbitrary; however, any redefinition is compensated by suitable changes in the other terms contributing to the cross section. The initial-state collinear terms are:

$$d\widehat{\sigma }_{ab}^{(c\pm )}=d\sigma _{ab}^{(c\pm )}|_{\mathrm{ev}}d\sigma _{ab}^{(c\pm )}|_{\mathrm{ct}},$$
(62)
where

$`d\sigma _{ab}^{(c+)}(p_1,p_2)|_{\mathrm{ev}}`$ $`=`$ $`{\displaystyle \frac{\alpha _\mathrm{s}}{2\pi }}\left[\mathrm{log}{\displaystyle \frac{s(1x)^2}{\mu ^2}}P_{ca}^{(0)}(x)P_{ca}^{(1)}(x)\right]_{cb}^{(b)}(xp_1,p_2)d\varphi _2(xs)dx,`$
$`d\sigma _{ab}^{(c)}(p_1,p_2)|_{\mathrm{ev}}`$ $`=`$ $`{\displaystyle \frac{\alpha _\mathrm{s}}{2\pi }}\left[\mathrm{log}{\displaystyle \frac{s(1x)^2}{\mu ^2}}P_{cb}^{(0)}(x)P_{cb}^{(1)}(x)\right]_{ac}^{(b)}(p_1,xp_2)d\varphi _2(xs)dx.`$
In eq. (62), the term $`d\sigma _{ab}^{(c\pm )}|_{\mathrm{ct}}`$ (the “ct” stands for counterterm) is constructed in order to cancel the divergences present in the collinear “event” contribution $`d\sigma _{ab}^{(c\pm )}|_{\mathrm{ev}}`$; these divergences are soft ($`x1`$), as can be seen explicitly in eqs. (LABEL:cpev) and (LABEL:cmev). The counterterm has a form that is fully specified in a given subtraction scheme (see appendix A). In eqs. (LABEL:cpev) and (LABEL:cmev), $`\mu `$ is the factorization mass scale (which also enters in the parton densities), and $`P_{ab}^{(i)}`$ is the $`i^{th}`$ term in the $`ϵ`$-expansion of the Altarelli-Parisi kernel in $`42ϵ`$ dimensions, taken for $`x<1`$ (i.e., no $`\delta (1x)`$ terms are included). Finally, we write the real emission term as

$`d\widehat{\sigma }_{ab}^{(f)}=d\sigma _{ab}^{(f)}|_{\mathrm{ev}}d\sigma _{ab}^{(f)}|_{\mathrm{ct}},`$ (65)
$`d\sigma _{ab}^{(f)}(p_1,p_2)|_{\mathrm{ev}}=_{ab}^{(r)}(p_1,p_2)d\varphi _3(s),`$ (66)
where $`_{ab}^{(r)}`$ is the matrix element squared resulting from the real emission diagrams (summed over all spins and colours, and averaged over those relevant to incoming partons), times the flux factor. This quantity has soft and initial-state collinear singularities, and thus $`d\sigma _{ab}^{(f)}|_{\mathrm{ev}}`$ has divergences, which are cancelled by the counterterms $`d\sigma _{ab}^{(f)}|_{\mathrm{ct}}`$ in eq. (65). As in the case of $`d\sigma _{ab}^{(c\pm )}|_{\mathrm{ct}}`$, these counterterms are fully specified only in a given subtraction scheme. However, one should keep in mind that the physical cross section is always strictly independent of the choice of the subtraction scheme.

### 4.3 Observables at NLO

As in the case of the toy model, a preliminary step in the construction of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ is the prediction of an observable at the NLO. We start by defining an unintegrated, fully-differential cross section

$$d\widehat{\mathrm{\Sigma }}_{ab}(x_1,x_2)=f_a^{(H_1)}(x_1)f_b^{(H_2)}(x_2)d\widehat{\sigma }_{ab}(x_1P_1,x_2P_2),$$
(67)
and analogous quantities, for all the contributions to the partonic cross section that appear on the r.h.s. of eq. (59). Any observable $`O`$ can depend, at the NLO, upon the kinematical variables of the W+, of the W-, and of the light outgoing parton. We denote by $`O(\mathrm{𝟑})`$ the observable computed in non-soft and non-collinear kinematics; this corresponds to $`O(x)`$ in the toy model. We also denote by $`O(\mathrm{𝟐})`$ and $`O(\stackrel{~}{\mathrm{𝟐}})`$ the observable computed in $`22`$ (i.e., soft) and in quasi-$`22`$ (i.e., collinear non-soft) kinematics respectively; these correspond to $`O(0)`$ in the toy model. In all cases, we understand the variables of the outgoing particles to be given in the laboratory frame. We then have

$`O`$ $`=`$ $`{\displaystyle \underset{ab}{}}{\displaystyle }dx_1dx_2d\varphi _3[O(\mathrm{𝟑}){\displaystyle \frac{d\mathrm{\Sigma }_{ab}^{(f)}}{d\varphi _3}}|_{\mathrm{ev}}+O(\mathrm{𝟐}){\displaystyle \frac{1}{_2}}({\displaystyle \frac{d\mathrm{\Sigma }_{ab}^{(b)}}{d\varphi _2}}+{\displaystyle \frac{d\mathrm{\Sigma }_{ab}^{(sv)}}{d\varphi _2}})`$ (68)
$`+O(\stackrel{~}{\mathrm{𝟐}}){\displaystyle \frac{1}{_{\stackrel{~}{2}}}}\left({\displaystyle \frac{d\mathrm{\Sigma }_{ab}^{(c+)}}{d\varphi _2dx}}|_{\mathrm{ev}}+{\displaystyle \frac{d\mathrm{\Sigma }_{ab}^{(c)}}{d\varphi _2dx}}|_{\mathrm{ev}}\right)O(\mathrm{𝟐}){\displaystyle \frac{1}{_{\stackrel{~}{2}}}}\left({\displaystyle \frac{d\mathrm{\Sigma }_{ab}^{(c+)}}{d\varphi _2dx}}|_{\mathrm{ct}}+{\displaystyle \frac{d\mathrm{\Sigma }_{ab}^{(c)}}{d\varphi _2dx}}|_{\mathrm{ct}}\right)`$
$`\{O(\mathrm{𝟐}),O(\stackrel{~}{\mathrm{𝟐}})\}{\displaystyle \frac{d\mathrm{\Sigma }_{ab}^{(f)}}{d\varphi _3}}|_{\mathrm{ct}}],`$
where in the last term $`\{O(\mathrm{𝟐}),O(\stackrel{~}{\mathrm{𝟐}})\}`$ indicates that the counterterms of the real contribution can have both $`22`$ and quasi-$`22`$ kinematics. The factors $`_2`$ and $`_{\stackrel{~}{2}}`$ guarantee that the $`22`$ and quasi-$`22`$ terms are correctly normalized upon integration over the full three-body phase space; the following equations must therefore hold

$`{\displaystyle 𝑑\varphi _3(s)}`$ $`=`$ $`_2{\displaystyle 𝑑\varphi _2(s)},`$ (69)
$`{\displaystyle 𝑑\varphi _3(s)}`$ $`=`$ $`_{\stackrel{~}{2}}{\displaystyle 𝑑\varphi _2(xs)𝑑x},`$ (70)
the integration being performed over the whole phase space. Eq. (68) should be compared with eq. (13). The factors $`1/_2`$ and $`1/_{\stackrel{~}{2}}`$ are not present in the toy model, simply because the integral over the phase space $`dx`$ in that case is equal to 1; if this were not the case, a factor analogous to $`1/_2`$ would have to be inserted in eq. (13). The term $`d\mathrm{\Sigma }_{ab}^{(f)}/d\varphi _3|_{\mathrm{ev}}=`$$`f_af_b_{ab}^{(r)}`$ corresponds to $`aR(x)/x`$ in eq. (13). The term $`d\mathrm{\Sigma }_{ab}^{(b)}/d\varphi _2=`$$`f_af_b_{ab}^{(b)}`$ corresponds to $`B`$, the Born contribution, while the term $`d\mathrm{\Sigma }_{ab}^{(sv)}/d\varphi _2=`$ $`f_af_b(𝒬_{ab}_{ab}^{(b)}+_{ab}^{(v,reg)})`$ corresponds to $`V`$, the finite part of the virtual contribution. As we said already, the contributions $`d\mathrm{\Sigma }_{ab}^{(c\pm )}`$ have no analogues in the toy model, simply because there are no parton densities there. Finally, $`d\mathrm{\Sigma }_{ab}^{(f)}/d\varphi _3|_{\mathrm{ct}}`$ represents the counterterms necessary to cancel the divergences of the real matrix element; it is thus equivalent to the last term in eq. (13), $`aB/x`$.

### 4.4 Matching NLO and MC: the QCD $`\mathrm{𝐌𝐂}\mathbf{@}\mathrm{𝐍𝐋𝐎}`$

Eq. (68) is very similar to eq. (13). We could try to follow the same steps that brought us from eq. (13) to eq. (20): replace the observables in eq. (68) with suitable interface-to-MC’s and then modify the equation with a quantity analogous to that introduced in eq. (21). There is, however, a technical problem, related to the initial conditions of class $`𝕊`$ events. In the toy model, these initial conditions are unambiguously determined by letting $`x0`$ in the corresponding system-plus-one-photon configurations. This procedure cannot be automatically generalized to QCD, where matrix elements have soft and collinear singularities (corresponding to $`\mathrm{𝟐}`$ and $`\stackrel{~}{\mathrm{𝟐}}`$ configurations respectively). This may suggest the necessity of considering two $`𝕊`$ classes, defined as the sets of those events resulting from MC evolutions that start from $`\mathrm{𝟐}`$ and $`\stackrel{~}{\mathrm{𝟐}}`$ initial conditions respectively. However, the MC cannot handle $`\stackrel{~}{\mathrm{𝟐}}`$ initial conditions, since final-state light partons are strictly collinear to the beam-line in these cases.

The solution can be found by observing that, at least as far as the variables of the W+ and W- are concerned, $`\mathrm{𝟐}`$ and $`\stackrel{~}{\mathrm{𝟐}}`$ only differ by a longitudinal boost. The variables of the outgoing light parton are of no concern here, since they cannot contribute to the kinematics of any observables in the configurations $`\mathrm{𝟐}`$ and $`\stackrel{~}{\mathrm{𝟐}}`$, the parton being either soft of collinear to the beam line. The integration over the momentum fractions $`x_1`$ and $`x_2`$ in eq. (68) implicitly contains such a boost. Roughly speaking, for fixed $`x_1`$ and $`x_2`$ we have $`\mathrm{𝟐}\stackrel{~}{\mathrm{𝟐}}`$, but given $`x_1`$ and $`x_2`$, there exist $`\stackrel{~}{x}_1`$ and $`\stackrel{~}{x}_2`$ such that $`\mathrm{𝟐}(x_i)=\stackrel{~}{\mathrm{𝟐}}(\stackrel{~}{x}_i)`$. In practice, we adopt the following procedure.666To the best of our knowledge, this method was invented in order to improve the numerical stability of the results of ref. . We treat each term under the integral sign in eq. (68) independently (we assume some regularization procedure has been adopted in the intermediate steps of the computation), performing a change of integration variables that is in general different for each term:

$$(x_1,x_2)(z_1,z_2),x_1=x_1^{(\alpha )}(z_1,z_2,\varphi _3),x_2=x_2^{(\alpha )}(z_1,z_2,\varphi _3),$$
(71)
where $`\varphi _3`$ denotes here the set of independent three-body phase-space variables. The functional forms of these changes of variables will be given in app. A.4. Here we are only interested in the fact that eq. (71) must be such that $`\mathrm{𝟐}(x_i)=\stackrel{~}{\mathrm{𝟐}}(\stackrel{~}{x}_i)`$. We define, for each of the functions $`\mathrm{\Sigma }_{ab}`$ that contribute to eq. (68),

$$d\overline{\mathrm{\Sigma }}_{ab}^{(\alpha )}(z_1,z_2,\varphi _3)=\frac{(x_1^{(\alpha )},x_2^{(\alpha )})}{(z_1,z_2)}d\mathrm{\Sigma }_{ab}^{(\alpha )}(x_1^{(\alpha )}(z_1,z_2,\varphi _3),x_2^{(\alpha )}(z_1,z_2,\varphi _3),\varphi _3),$$
(72)
where $`(x_1^{(\alpha )},x_2^{(\alpha )})/(z_1,z_2)`$ is the Jacobian corresponding to eq. (71). Then, as we shall show in app. A.4, eq. (68) can be recast in the following form:

$`O`$ $`=`$ $`{\displaystyle \underset{ab}{}}{\displaystyle }dz_1dz_2d\varphi _3\{O(\mathrm{𝟑}){\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(f)}}{d\varphi _3}}|_{\mathrm{ev}}+O(\mathrm{𝟐})[{\displaystyle \frac{1}{_2}}({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(b)}}{d\varphi _2}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(sv)}}{d\varphi _2}})`$ (73)
$`+{\displaystyle \frac{1}{_{\stackrel{~}{2}}}}\left({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c+)}}{d\varphi _2dx}}|_{\mathrm{ev}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c)}}{d\varphi _2dx}}|_{\mathrm{ev}}\right){\displaystyle \frac{1}{_{\stackrel{~}{2}}}}\left({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c+)}}{d\varphi _2dx}}|_{\mathrm{ct}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c)}}{d\varphi _2dx}}|_{\mathrm{ct}}\right)`$
$`{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(f)}}{d\varphi _3}}|_{\mathrm{ct}}]\}.`$
We stress that the only difference between eq. (68) and eq. (73) is that the latter does not contain $`O(\stackrel{~}{\mathrm{𝟐}})`$ any longer, thanks to the definition of the change of variables, which ensures that $`\mathrm{𝟐}(x_i)=\stackrel{~}{\mathrm{𝟐}}(\stackrel{~}{x}_i)`$. We also stress that no approximation has been made in eq. (73); it is just a different way of writing the NLO cross section.

It should be clear that the necessity for the change of variables in eq. (71) is due to the simultaneous presence in eq. (68) of $`O(\mathrm{𝟐})`$ and $`O(\stackrel{~}{\mathrm{𝟐}})`$, which is in turn due to the form of the hadronic cross section in eq. (55), where the luminosity is factorized. One could directly start with an NLO cross section that does not contain $`O(\stackrel{~}{\mathrm{𝟐}})`$; clearly, our formulae hold in this case as well, since one just needs to set $`(z_1,z_2)(x_1,x_2)`$.

It should also be clear that our goal is to obtain an expression for the NLO cross section that only contains $`\mathrm{𝟑}`$ and $`\mathrm{𝟐}`$ configurations. Although we find that we can easily achieve this through eq. (71), we stress that other methods could be devised, and that our final formula for the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ will hold, regardless of the procedure adopted to eliminate $`\stackrel{~}{\mathrm{𝟐}}`$ configurations.

The quantities that contribute to $`O(\mathrm{𝟑})`$ and $`O(\mathrm{𝟐})`$ in eq. (73) correspond to those that contribute to $`O(x)`$ and $`O(0)`$ in the toy model, respectively, according to the analogies discussed in sect. 4.2. The only exceptions are the terms $`d\widehat{\sigma }_{ab}^{(c\pm )}`$, which do not have analogues in the toy model. However, as shown in app. A.4, these terms must be treated in the same way as the real matrix element collinear counterterms. This is in fact to be expected: $`d\widehat{\sigma }_{ab}^{(c\pm )}`$ are strictly related to the collinear counterterms, being finite remainders resulting from the subtraction of the collinear singularities of the real matrix elements.

Equation (73) is therefore strictly analogous to eq. (13). This allows us to extend to the QCD case the concept of the event classes $``$ and $`𝕊`$. In turn, this implies that we need two interface-to-MC’s, as in the toy model, which we denote by $`I_{\mathrm{MC}}(O,\mathrm{𝟑})`$ and $`I_{\mathrm{MC}}(O,\mathrm{𝟐})`$. These are the analogues of $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$ and $`I_{\mathrm{MC}}(O,1)`$ in the toy model. We can thus turn eq. (73) into a prescription for a $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, by replacing $`O(\mathrm{𝟑})`$ by $`I_{\mathrm{MC}}(O,\mathrm{𝟑})`$, and $`O(\mathrm{𝟐})`$ by $`I_{\mathrm{MC}}(O,\mathrm{𝟐})`$ . In this way, we get the analogue of eq. (19), that is, a “naive” $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ that suffers from double counting and difficulties in the generation of unweighted events. However, precisely as in the case of the toy model, we can modify the naive prescription in order to get a QCD $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ that implements a modified subtraction, and fulfills the requirements given in sect. 2.1.

To accomplish this task, we need to find the analogue of the quantity $`aBQ(x)/x`$ of the toy model. As discussed in sect. 4.1, what we are looking for is the $`𝒪(\alpha _\mathrm{s})`$ term in the expansion of the result obtained by running the ordinary MC; we denote this quantity by

$$\frac{d\overline{\mathrm{\Sigma }}_{ab}}{d\varphi _3}|_{\mathrm{MC}}.$$
(74)
The specific form of the quantity in eq. (74) is dependent upon the implementation; different MC’s will induce different forms of $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$. On the other hand, $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ is independent of the procedure through which we eliminated $`\stackrel{~}{\mathrm{𝟐}}`$ configurations. The form we use in the present paper, and all the technical details concerning its determination, are reported in app. A.5. Here, we merely need to note that $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ must have the same leading singular behaviour as the matrix element contribution $`d\overline{\mathrm{\Sigma }}_{ab}^{(f)}`$ when approaching soft or collinear configurations, in order for the former quantity to be used as a local counterterm to the latter. As we shall discuss in app. A.5, this fact is not trivial, and we shall need to modify the $`𝒪(\alpha _\mathrm{s})`$ result obtained from the MC. For the time being, we treat $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ in the same way as $`aBQ(x)/x`$ is treated in the toy model (see eq. (21)): we subtract and add the quantity

$$I_{\mathrm{MC}}(O,𝐧)\frac{d\overline{\mathrm{\Sigma }}_{ab}}{d\varphi _3}|_{\mathrm{MC}}$$
(75)
under the integral sign, setting $`𝐧=\mathrm{𝟑}`$ in the negative term and $`𝐧=\mathrm{𝟐}`$ in the positive one. This gives a modified subtraction formula analogous to that in eq. (20):

$`{\displaystyle \frac{d\sigma }{dO}}`$ $`=`$ $`{\displaystyle \underset{ab}{}}{\displaystyle }dz_1dz_2d\varphi _3\{I_{\mathrm{MC}}(O,\mathrm{𝟑})({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(f)}}{d\varphi _3}}|_{\mathrm{ev}}{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}}{d\varphi _3}}|_{\mathrm{MC}})`$ (76)
$`+I_{\mathrm{MC}}(O,\mathrm{𝟐})[{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(f)}}{d\varphi _3}}|_{\mathrm{ct}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}}{d\varphi _3}}|_{\mathrm{MC}}+{\displaystyle \frac{1}{_2}}({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(b)}}{d\varphi _2}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(sv)}}{d\varphi _2}})`$
$`+{\displaystyle \frac{1}{_{\stackrel{~}{2}}}}({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c+)}}{d\varphi _2dx}}|_{\mathrm{ev}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c)}}{d\varphi _2dx}}|_{\mathrm{ev}}){\displaystyle \frac{1}{_{\stackrel{~}{2}}}}({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c+)}}{d\varphi _2dx}}|_{\mathrm{ct}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c)}}{d\varphi _2dx}}|_{\mathrm{ct}})]\}.`$
Eq. (76) is our proposal for a QCD $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. In order to use this prescription, we need to write it in a form more suited to implementation in a computer code, which will be done in the next section. We stress that, as in the case of the toy model, the dependence upon $`O`$ in eq. (76) is only formal: the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ generates events without reference to any observable.

### 4.5 QCD $`\mathrm{𝐌𝐂}\mathbf{@}\mathrm{𝐍𝐋𝐎}`$: practical implementation

The straightforward implementation of eq. (76) results in a computer code which is not optimally efficient, in particular in the case of unweighted event generation. We therefore implement our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ through the following procedure.

We start by computing the following integrals:

$`I_{}`$ $`=`$ $`{\displaystyle \underset{ab}{}}{\displaystyle 𝑑z_1𝑑z_2𝑑\varphi _3\left(\frac{d\overline{\mathrm{\Sigma }}_{ab}^{(f)}}{d\varphi _3}|_{\mathrm{ev}}\frac{d\overline{\mathrm{\Sigma }}_{ab}}{d\varphi _3}|_{\mathrm{MC}}\right)}`$ (77)
$`I_𝕊`$ $`=`$ $`{\displaystyle \underset{ab}{}}{\displaystyle }dz_1dz_2d\varphi _3[{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(f)}}{d\varphi _3}}|_{\mathrm{ct}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}}{d\varphi _3}}|_{\mathrm{MC}}`$ (78)
$`+{\displaystyle \frac{1}{_2}}\left({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(b)}}{d\varphi _2}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(sv)}}{d\varphi _2}}\right)+{\displaystyle \frac{1}{_{\stackrel{~}{2}}}}\left({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c+)}}{d\varphi _2dx}}|_{\mathrm{ev}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c)}}{d\varphi _2dx}}|_{\mathrm{ev}}\right)`$
$`{\displaystyle \frac{1}{_{\stackrel{~}{2}}}}({\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c+)}}{d\varphi _2dx}}|_{\mathrm{ct}}+{\displaystyle \frac{d\overline{\mathrm{\Sigma }}_{ab}^{(c)}}{d\varphi _2dx}}|_{\mathrm{ct}})].`$
Note that these two integrals are finite, owing to the properties of $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$. Clearly,

$$\sigma _{tot}=I_𝕊+I_{},$$
(79)
where $`\sigma _{tot}`$ is the total rate at the NLO level. We also compute the integrals $`J_{}`$ and $`J_𝕊`$, defined as in eqs. (77) and (78) respectively, except that the integrands are taken in absolute value (so that, if the integrands were positive-definite, we would have $`I_i=J_i`$). We use a modified version of the adaptive integration routine BASES , which allows one to compute at the same time the integral of a given function, and the integral of the absolute value of that function. BASES is an evolution of the original VEGAS package .

We then make use of the routine SPRING, which is part of the BASES package. Roughly speaking, for any positive-definite function $`f(\stackrel{}{x})`$, SPRING returns a set of values $`\{\stackrel{}{x}_i\}_{i=1}^n`$ distributed according to the function $`f`$; the number $`n`$ must be given as an input. Thus, SPRING effectively achieves the generation of $`n`$ unweighted events (with weights equal to 1). In our case, the absolute values of the integrands in eqs. (77) and (78) are given as input functions to SPRING. The numbers of events generated in the two cases are (possibly with a roundoff to the next integer)

$$N_{}=N_{tot}\frac{J_{}}{J_𝕊+J_{}},N_𝕊=N_{tot}\frac{J_𝕊}{J_𝕊+J_{}},$$
(80)
where $`N_{tot}`$ is the total number of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ events that we eventually want to generate. It is easy, within SPRING, to keep track of the sign of the integrand function, before its absolute value is taken. Thus, SPRING achieves unweighted event generation, with the weights equal to +1 or –1.

We symbolically denote the results of the SPRING runs as follows:

$`I_{}`$ $`\stackrel{\mathrm{SPRING}}{}`$ $`\{\mathrm{𝟑}_i,𝐟_i,w_i^{()}\}_{i=1}^N_{},`$ (81)
$`I_𝕊`$ $`\stackrel{\mathrm{SPRING}}{}`$ $`\{\mathrm{𝟐}_i,𝐟_i,w_i^{(𝕊)}\}_{i=1}^{N_𝕊};`$ (82)
thus, eqs. (81) and (82) correspond to events belonging to classes $``$ and $`𝕊`$ respectively. For a given $`i`$, $`\{\mathrm{𝟑}_i,𝐟_i,w_i^{()}\}`$ gives the complete description of a three-body event; $`\mathrm{𝟑}_i`$ represents the kinematics (in the lab frame), $`𝐟_i=\{a_i,b_i,c_i\}`$ are the flavours of the initial- $`(a_i,b_i)`$ and final-state $`(c_i)`$ partons, and $`w_i^{()}=\pm 1`$ is the weight of the event returned by SPRING. The event $`\{\mathrm{𝟐}_i,𝐟_i,w_i^{(𝕊)}\}`$ is analogous, except that it has $`22`$ kinematics (and thus $`𝐟_i=\{a_i,b_i\}`$).

The events in eqs. (81) and (82) are the initial conditions for MC evolution. In the case of $``$ events, eq. (81), the kinematics $`\mathrm{𝟑}_i`$ is determined by $`(z_1,z_2,\varphi _3)`$. The choice of the flavours $`𝐟_i`$ is more involved. One possibility is to fix the pair of flavours $`(a,b)`$ that appears in the sum in eq. (77) (and thus to fix $`𝐟`$, since $`c`$ in this case is unambiguously determined by $`a`$ and $`b`$), to compute the integral, and to repeat the computation for all possible flavour pairs. This procedure is however inconvenient from the point of view of computing time. A faster approach is simply to compute three integrals of the type in eq. (77), for the cases $`ab=q\overline{q}`$, $`qg`$ and $`\overline{q}g`$ summed over quark flavours. For each of these integrals, the quark flavours are then determined by statistical methods; this can give a wrong assignment of flavours on event-by-event basis, but guarantees that the “distribution” in the flavours is returned correctly. The reason for not computing a single integral is that subprocesses with $`q\overline{q}`$, $`qg`$ and $`\overline{q}g`$ initial states often give contributions of different signs to the cross section, which is a source of errors when statistical methods are used.

In the case of $`𝕊`$ events, eq. (82), the kinematic configuration $`\mathrm{𝟐}_i`$ is unambiguously defined as the soft limit of the corresponding $`23`$ kinematics, determined by $`(z_1,z_2,\varphi _3)`$. For the flavour assignment, we proceed similarly to what was done in the case of $``$ events. However, here the relevant flavours are those of the partons entering the $`22`$ hard process which factorizes in the case of soft or collinear emissions. Thus, we use the Born contribution $`d\overline{\mathrm{\Sigma }}_{ab}^{(b)}`$ to assign the flavour. Notice that this can result in errors of $`𝒪(\alpha _\mathrm{s})`$ in flavour assignment. However, the effects induced in this way are $`𝒪(\alpha _\mathrm{s}^2)`$, and thus beyond control. On the other hand, with this prescription one achieves at $`𝒪(\alpha _\mathrm{s})`$ the cancellation of the MC terms $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ introduced in the modified subtraction, i.e., one avoids double counting. More details on this point are given in app. B.3. It should be noticed that, since the Born contribution is always positive, there is no need to compute the integral in eq. (78) separately for $`q\overline{q}`$, $`qg`$ and $`\overline{q}g`$ initial states.

The above statistical method has the virtue of being easy to generalize to other processes. We point out that eqs. (79) and (80), together with those in the remainder of this section, are trivially extended to the case in which more than one integral is computed for $``$ and/or $`𝕊`$ events.

The weights in eqs. (81) and (82) have the following properties:

$$\frac{J_{}}{N_{}}\underset{i=1}{\overset{N_{}}{}}w_i^{()}=I_{},\frac{J_𝕊}{N_𝕊}\underset{i=1}{\overset{N_𝕊}{}}w_i^{(𝕊)}=I_𝕊,$$
(83)
and thus, using eq. (80),

$$\sigma _{tot}=\frac{J_{}+J_𝕊}{N_{tot}}\left(\underset{i=1}{\overset{N_{}}{}}w_i^{()}+\underset{i=1}{\overset{N_𝕊}{}}w_i^{(𝕊)}\right).$$
(84)
We then define the set of events

$$\{𝐤_i,𝐟_i,w_i\}_{i=1}^{N_{tot}}=\{\mathrm{𝟑}_i,𝐟_i,w_i^{()}\}_{i=1}^N_{}\{\mathrm{𝟐}_i,𝐟_i,w_i^{(𝕊)}\}_{i=1}^{N_𝕊},$$
(85)
where $`𝐤_i`$ denotes either $`\mathrm{𝟑}_i`$ or $`\mathrm{𝟐}_i`$, and

$$w_i=w_i^{(𝕊,)}\frac{J_{}+J_𝕊}{N_{tot}}\underset{i=1}{\overset{N_{tot}}{}}w_i=\sigma _{tot}.$$
(86)
We constructed a Fortran code whose output is the set of events $`\{𝐤_i,𝐟_i,w_i\}_{i=1}^{N_{tot}}`$. We point out that the merging of the two event sets on the r.h.s. of eq. (85) is done randomly; therefore, the probability of finding a $`𝕊`$ ($``$) event at an arbitrary place in the event list is a constant, equal to $`N_𝕊/N_{tot}`$ ($`N_{}/N_{tot}`$). This set of events is passed to HERWIG; there, it replaces the routine that generates the hard process. We stress that HERWIG has not been modified; we just added a routine that reads the parton configurations coming from the NLO code, and translates them in the standard HERWIG form. Each configuration generates a shower, and all the other operations requested by the user. Notice that the unweighting procedure is entirely dealt with by the NLO code; HERWIG reads the weights provided by the NLO code, and does not discard any events (apart from the few killed in the ordinary shower/hadronization process). According to eq. (86), in order to get the total cross section right, exactly $`N_{tot}`$ events must be processed by HERWIG . However, should only $`N_{ev}<N_{tot}`$ events be processed, the correct rate will be recovered by multiplying all the physical results by a factor of $`N_{tot}/N_{ev}`$, since all weights have the same absolute value. In other words, any subset of the event set $`\{𝐤_i,𝐟_i,w_i\}_{i=1}^{N_{tot}}`$ can act as an event set itself, provided that the weights are multiplied by a common factor, in order for eq. (86) to hold.

The presence of negative-weight events is a distinctive feature of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$’s. Negative weights imply larger statistical errors than in the case of ordinary MC’s, when the same total number of events is generated. In the $``$ and $`𝕊`$ samples of eqs. (81) and (82), the fraction of negative-weight events is

$$f_i=\frac{1}{2}\left(1\frac{I_i}{J_i}\right),$$
(87)
for $`i=𝕊,`$; clearly, the closer $`J_i`$ is to $`I_i`$, the smaller is $`f_i`$. We can in fact reduce the number of negative-weight events in our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ by varying $`J_𝕊`$. This quantity depends upon the form of the counterevents, or, in a given subtraction scheme, upon the free parameters that define the subtraction777These are $`\stackrel{~}{\rho }`$ and $`\omega `$ in the standard subtraction, $`\zeta `$ in the $`\zeta `$-subtraction: see app. A.2 and app. A.3.. Notice that this is not in contradiction with the properties of the subtraction method, since $`J_𝕊`$ is not a physical quantity. Thus, we can adjust these free parameters in order to minimize the difference $`J_𝕊I_𝕊`$; the solution can only be found through numerical methods, because of the complexity of the integrand appearing in $`J_𝕊`$. Further details on this point are given in app. A.3.

We point out that the integral $`J_{}`$ does not depend upon any free parameters; thus, the number of negative-weight $``$ events is a constant for fixed input (physical) parameters in eq. (76). We remind the reader that these events will mainly populate the hard emission region.

We remark finally that, owing to the procedure described in sect. 4.4, the $`22`$ kinematic configurations resulting from the integral in eq. (78) can all be obtained as the soft limit of the $`23`$ kinematic configurations identified by $`(z_1,z_2,\varphi _3)`$. This suggests writing the integration measure in eq. (78) with the two-body phase space explicitly factorized, $`d\varphi _3=d\varphi _2d\stackrel{~}{\varphi }_3`$, and integrating over $`d\stackrel{~}{\varphi }_3`$ before passing the integrand function to SPRING. It should be clear that this intermediate integration step does not affect our capability of obtaining fully-exclusive final states; however, it does have an impact on the number of negative-weight events produced by the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. In this paper, we do not exploit this possibility, leaving the option open for future developments.

### 4.6 Generalities on matched computations

We now consider the case of an observable $`O`$ whose kinematically-allowed range is $`O_s<O<O_h`$. We assume that the region $`OO_h`$ receives only contributions due to hard emissions, and that large logarithms are present in any fixed-order cross section for $`OO_s`$, where a resummation is therefore necessary. In accordance with eq. (76), we write a differential cross section obtained from the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ as follows:

$$\frac{d\sigma }{dO}=\left(\frac{d\sigma }{dO}\right)_\mathrm{H}\left(\frac{d\sigma }{dO}\right)_{\mathrm{RS0}}+\left(\frac{d\sigma }{dO}\right)_{\mathrm{RS}}.$$
(88)
The quantity $`(d\sigma /dO)_{\mathrm{RS}}`$ is the contribution of all $`𝕊`$ events, i.e. of the MC evolutions $`I_{\mathrm{MC}}(O,\mathrm{𝟐})`$. The quantities $`(d\sigma /dO)_\mathrm{H}`$ and $`(d\sigma /dO)_{\mathrm{RS0}}`$ are the contributions of the first and second terms under the integral sign in eq. (76) respectively, and thus both are associated with $``$ events. We expect them to be divergent for $`OO_s`$, only their difference being finite. However, for $`OO_h`$, we can assume that $`(d\sigma /dO)_\mathrm{H}`$ and $`(d\sigma /dO)_{\mathrm{RS0}}`$ are separately finite. This is because this region is by hypothesis dominated by hard emission, which can only occur through real-emission diagrams; this prevents the emitted parton over whose variables we integrate in eq. (76) from being soft or collinear. The quantities introduced in eq. (88) have the following properties:

$`\left|\left({\displaystyle \frac{d\sigma }{dO}}\right)_\mathrm{H}\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{RS0}}\right|\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{RS}}\mathrm{for}OO_s;`$ (89)
$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_\mathrm{H}=\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{NLO}}+𝒪(\alpha _\mathrm{s}^2)\mathrm{for}OO_h;`$ (90)
$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{RS}}=\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{RS0}}+𝒪(\alpha _\mathrm{s}^2)\mathrm{for}OO_h.`$ (91)
Eqs. (90) and (91) imply that, for $`OO_h`$, we recover the expected result $`d\sigma /dO=(d\sigma /dO)_{\mathrm{NLO}}`$, which holds up to $`𝒪(\alpha _\mathrm{s}^2)`$ (NNLO) terms. Strictly speaking, eq. (91) is not necessary for this to happen; $`(d\sigma /dO)_{\mathrm{RS}}`$ and $`(d\sigma /dO)_{\mathrm{RS0}}`$ include the leading (and possibly next-to-leading) tower of logarithms $`\mathrm{log}O/O_h`$, or more precisely $`\mathrm{log}[(OO_s)/(O_hO_s)]`$, which become small for $`OO_h`$. However, eq. (91) implies that $`d\sigma /dO(d\sigma /dO)_{\mathrm{NLO}}`$ also for those values of $`O`$ which are not too close to $`O_h`$, but still outside the region where resummation is necessary. The point here is the size of the unknown $`𝒪(\alpha _\mathrm{s}^2)`$ coefficients. There is in fact no guarantee that these coefficients are small; if they are big, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ result can differ substantially from the NLO one, although the accuracy of the two computations is the formally of the same order in perturbation theory. This problem is always encountered when a fixed-order computation is matched with a resummed computation. A possible solution (see for example refs. ) would be to rewrite eq. (88) as follows

$$\frac{d\sigma }{dO}=\left(\frac{d\sigma }{dO}\right)_\mathrm{H}+𝒟(O)\left[\left(\frac{d\sigma }{dO}\right)_{\mathrm{RS}}\left(\frac{d\sigma }{dO}\right)_{\mathrm{RS0}}\right],$$
(92)
where $`𝒟(O)`$ is a regular function, such that $`𝒟(O_h)0`$ for $`OO_h`$, and $`𝒟(O)1`$ for $`OO_s`$, up to terms suppressed by powers of $`(OO_s)/(O_hO_s)`$. The perturbative and logarithmic accuracies of eq. (92) are the same as those of eq. (88); however, the function $`𝒟`$ can be suitably chosen in order to suppress the “MC” contribution in the region in which the NLO computation can be trusted. Notice that the presence of $`𝒟(O)`$ implies that the total rate predicted by eq. (92) differs from that predicted at the NLO level, even when no cuts are applied.

The problem with eq. (92) is that the function $`𝒟`$ is specific to the observable considered; thus, eq. (92) cannot be implemented directly in an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. However, one could try to devise a solution, similar to that of eq. (92), that works in the context of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$’s. Clearly, in the hard-emission region the function $`𝒟`$ suppresses the contribution of $`𝕊`$ events, and that part of $``$ events that corresponds to the $`𝒪(\alpha _\mathrm{s})`$ expansion of the MC result. In order to achieve this result without any reference to a specific observable, a modification of the showering mechanism in the MC is necessary. In particular, one could multiply the Altarelli-Parisi kernel appearing in the Sudakov by a function $`Q(\xi ,z)`$ of the showering variables $`\xi `$ and $`z`$, that tends to 1 in the soft and/or collinear limits, and to zero in the hard emission region. It should be noted that this function is completely analogous to the function $`Q(x)`$ that we used in the toy MC, eq. (14). Exactly as in the case of the toy model, in order to satisfy the requirement that the $`𝒪(\alpha _\mathrm{s})`$ expansion of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ result be identical to the NLO result, we would need to perform the formal substitution

$$d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}Q(\xi ,z)d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}},$$
(93)
in eq. (76). Since this procedure would modify the integrand in eqs. (77) and (78), it would also have an impact on the number of negative-weight events generated. In this paper, we shall not pursue this issue further, since we prefer to keep the parton shower MC unmodified.

### 4.7 Results: W+W- and jet observables

In this section we present results obtained with our W+W- $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, implemented as described in sect. 4.5. We do not seek to give a phenomenological description that corresponds to a specific experimental configuration; rather, we wish to show the differences between $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, standard HERWIG MC, and NLO results. For this reason, we also present distributions obtained by integrating over the whole phase space, i.e., in most cases we do not apply any acceptance cuts. In the case in which cuts are applied, we only consider the W+ and W- variables, and not the variables of their decay products. We also consider jet distributions, the jets being reconstructed by means of the $`k_\mathrm{T}`$-clustering algorithm . We show results for $`pp`$ collisions at $`\sqrt{S}=14`$ TeV. We adopt the low-$`\alpha _\mathrm{s}`$ set of MRST99 parton densities , since this set has a $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ value which is rather close to that used as HERWIG default; this value of $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ ($`\mathrm{\Lambda }_5^{\overline{\mathrm{MS}}}=164`$ MeV) is used for all our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, MC, and NLO runs. In the case of standard HERWIG MC runs, we give each event (we generate unweighted events) a weight equal to $`\sigma _{tot}/N_{tot}`$, with $`N_{tot}`$ the total number of events generated. For fully exclusive distributions with no cuts applied, this is equivalent to normalizing HERWIG results to the total NLO rate $`\sigma _{tot}`$. All the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ and MC results (but not, of course, the NLO ones) include the hadronization of the partons in the final state.

We present in fig. 3 the rapidity of the W+, with (right panel) or without (left panel) a cut on the transverse momentum $`p_\mathrm{T}^{\mathrm{W}^+}`$. The solid, dashed and dotted histograms show the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, MC and NLO results, respectively. This distribution is fairly inclusive, and we expect it to be predicted well by NLO QCD in a wide range. From the figure, we see that the NLO and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results are extremely close to each other in the whole range considered, whether or not there is a $`p_\mathrm{T}^{\mathrm{W}^+}`$ cut. The MC result is also similar, but clearly broader than the other two.

A similarly inclusive variable is presented in fig. 4. The transverse momentum of the W+ is shown with (right panel) or without (left panel) a cut on the rapidity $`y^{\mathrm{W}^+}`$. The previous comment applies in this case as well: NLO and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ basically coincide, whereas the MC result is softer than the other two. The same pattern can for example be found for the invariant mass of the W pair. We conclude that, for this kind of observable, the $`𝒪(\alpha _\mathrm{s}^2)`$ effects are very small, and NLO and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ are almost equivalent. This also implies that any possible reshuffling of the momenta, due to the hadronization phase in $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, has negligible impact on the colourless W’s, as we expect.

We now turn to the case of more exclusive quantities, such as correlations between W+ and W- variables. In fig. 5 we present $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$, the modulus of the vector sum of the transverse momenta of the two W’s. NLO computations cannot predict this observable in the region $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}0`$, because of a logarithmic divergence for $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}0`$; on the other hand, NLO is expected to give reliable predictions at large $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$. The MC behaves in the opposite way; thanks to the cascade emission of soft and collinear partons, it can effectively resum the distribution around $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}=0`$; however, its results are not reliable in the large-$`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ region, which is mainly populated by events in which a very hard parton recoils against the W pair. Apart from the unreliability of the MC shower approximation to the matrix element in the large-$`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ region, this region cannot really be filled efficiently by the shower; this is easily understood by observing that, at fixed partonic c.m. energy, the large-$`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ region coincides with the dead zone (the border of the dead zone in terms of $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ can be worked out by using the results of app. A.5).

This complementary behaviour of NLO and MC approaches can be seen clearly in fig. 5, independent of the cuts on the rapidities and transverse momenta of the bosons. In the tail, the NLO cross section is much larger than the MC, simply because hard emissions are correctly treated only in the former. For $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}0`$, the difference between the two histograms shows the effect of all-order resummation; clearly, no meaningful comparison between NLO and data can be attempted in this region. It is therefore reassuring that the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ result interpolates the MC and NLO results smoothly. In the small-$`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ region, the shape of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ curve is identical to that of the MC result. This is evidence of the fact that MC and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ resum large logarithms at the same level of accuracy, as argued in app. B.4. When $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ grows large, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ tends to the NLO result, as expected. Again, the hadronization has no significant impact on the W+W- system.

In fig. 6 we present the distribution of the difference between the azimuthal scattering angles (i.e., those in the plane transverse to the beam direction) of the W+ and W-. This distribution cannot be reliably predicted by fixed-order QCD computations in the region $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}\pi `$; in fact, the NLO prediction diverges logarithmically for $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}\pi `$. This can be seen in the insets of the plots, where the variable $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}`$ has been plotted versus $`(\pi \mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)})/\pi `$ on a logarithmic scale, in order to visually enhance the region $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}\pi `$. We can see that in this region the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ and MC results have identical shapes, as in the case of the observable $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ near zero discussed above. The other end of the spectrum, i.e. the tail $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}0`$, might be considered as analogous to the region of large $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$, in that this region is dominated, in an NLO computation, by hard single-parton emission. However, there are important differences from the point of view of MC treatment. As discussed earlier, the region of high $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ is not filled efficiently by the MC shower. On the other hand, the $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}0`$ tail for low $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ is more easily populated by the MC. This is confirmed in the left-hand panel of fig. 6, where we can see that the NLO and MC results are quite close to each other at $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}0`$ when no cuts are applied. Thus, the region $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}0`$ receives contributions both from hard parton emissions (which dominate the configurations in which the W’s have large transverse momenta), and from multiple soft or collinear parton emissions (when the W’s have small transverse momenta). NLO can treat the former but not the latter, and MC can treat the latter but not the former. $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, on the other hand, is expected to handle both emissions correctly, while still avoiding any double counting. When cuts on the transverse momenta of the W’s are applied, the contribution from multiple soft or collinear emissions becomes less important, as may be seen in the right-hand panel of fig. 6. Overall, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results are close to the NLO predictions in the $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}0`$ tail; as we shall show later, the differences are compatible with $`𝒪(\alpha _\mathrm{s}^2)`$ effects.

We now turn to the discussion of jet observables. In fig. 7 we present the transverse energy distribution of the hardest jet of each event (left), and of the inclusive jets (right). In this case, we only compare $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ and MC results; NLO jet results are trivial for this process, since there is only one “jet”, which coincides with a parton. In the case of the hardest jet, the same discussion as in the case of $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ applies; the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ resums large logarithms at small $`E_\mathrm{T}`$ in the same way as the MC does, and can also treat the large-$`E_\mathrm{T}`$ region, where the MC fails. In the case of inclusive jets, the spectrum diverges for $`E_\mathrm{T}0`$, because of the increasing number of jets with smaller and smaller $`E_\mathrm{T}`$. In this region, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ and MC results (rescaled by the K-factor) coincide. We stress that this is actually the same pattern that we have already discussed in the toy model: see eqs. (45) and (46).

The NLO results presented so far have been obtained by setting the renormalization ($`\mu _\mathrm{R}`$) and factorization ($`\mu _\mathrm{F}`$) scales, that we collectively refer to as NLO scales ($`\mu _{\mathrm{NLO}}`$), equal to a reference scale $`\mu _0`$, defined as follows

$$\mu _0^2=\frac{1}{2}\left[\left(M_\mathrm{T}^{\mathrm{W}^+}\right)^2+\left(M_\mathrm{T}^\mathrm{W}^{}\right)^2\right],$$
(94)
where $`M_\mathrm{T}^{\mathrm{W}^+}`$ and $`M_\mathrm{T}^\mathrm{W}^{}`$ are the transverse masses of the W+ and W- respectively. The same scale choice has been adopted to compute the integrals given in eqs. (77) and (78), that is, to obtain the sets of $``$ and $`𝕊`$ events eventually given to the MC as initial conditions.

In general, we refer to the scales used in the computations of eqs. (77) and (78) as $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ scales ($`\mu _{\mathrm{MC}\mathrm{@}\mathrm{NLO}}`$). We point out that the NLO and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ scales need not coincide. It actually seems more natural to choose different forms for them. The $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ scale should in some way match the scale used in the parton shower; although any mismatch between the two is formally of higher order than we are considering, a choice of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ scale motivated by MC considerations should result in smaller coefficients for the $`𝒪(\alpha _\mathrm{s}^2)`$ terms in the regions dominated by hard emissions. Notice that, when the NLO and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ scales do not coincide, the total rates predicted by NLO and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ computations are different.

We shall not pursue the issue of an appropriate $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ scale choice any further in the present paper. We shall limit ourselves to considering the variation of NLO and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ predictions obtained by setting $`\mu _\mathrm{R}=\mu _\mathrm{F}=\mu _0/2`$ and $`\mu _\mathrm{R}=\mu _\mathrm{F}=2\mu _0`$. We stress that, while this procedure is a commonly accepted way to obtain a rough estimate of the unknown higher-order terms in an NLO computation, it cannot be used to guess the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ scale that best matches the parton shower scale; for this purpose, functional forms other than that shown in eq. (94) would need to be considered. Still, the procedure can give an idea of the size of the neglected higher-order terms in the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ formalism. The results are presented in fig. 8 for the azimuthal correlation $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}`$; we have chosen this observable since it is the one in which the scale effects are largest. The solid and dotted bands show the variation in the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ and NLO predictions respectively. We observe a sizable reduction of the scale dependence when going from NLO to $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results. Also, in the $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}0`$ tail, the two bands are rather close to each other, and overlap in the case that cuts are applied. This is an indication of the fact that NLO and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results indeed differ by $`𝒪(\alpha _\mathrm{s}^2)`$ terms, and that the difference is less significant when hard cuts are applied, as we should expect.

An issue related to the choice of scale for the parton showers is the setting of the HERWIG parameter EMSCA, which enforces an upper limit on the transverse momenta of partons emitted in the showers. The coherence of gluon emission is ensured in HERWIG by using an angular variable to perform the shower evolution. In order to avoid double-counting of hard processes, it is then necessary to veto emissions with transverse momenta greater than the scale of the primary hard process, which is specified by EMSCA. If such emissions were allowed, they would constitute the true primary hard process. When the veto is applied, the harder emission does not occur but the evolution variable is reset to the value selected for the vetoed emission. In this way the constraint that the true hard scale is EMSCA is correctly taken into account together with QCD coherence in the shower.

In our case, when running the standard HERWIG MC for W+W- production, we set EMSCA=$`\sqrt{s}`$. In the case of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, there is no reason to have the same EMSCA for events belonging to the classes $``$ and $`𝕊`$. We have already pointed out that the class $`𝕊`$ must behave in the same way as the ordinary MC, in order to reproduce (shape-wise) its results in those regions where a resummation is needed. Thus, for $`𝕊`$ events, we set should also set EMSCA=$`\sqrt{s}`$. As far as $``$ events are concerned, we notice that the initial condition for the MC in the toy model can also be interpreted as a veto scale; in the case of toy $``$ events, the initial condition is $`x_\mathrm{M}(x)`$, which is imposed to be a monotonic function, with $`x_\mathrm{M}(0)=1`$ (see sect. 3.3 and sect. 3.4). A possible choice that reproduces this condition in the real case is EMSCA=$`\sqrt{s}2p_\mathrm{T}`$, where $`p_\mathrm{T}`$ is the transverse momentum of the final-state light parton, which is always present in $``$ events. Then, in the soft or collinear limits, this scale tends to that adopted in the case of $`𝕊`$ events. However, such a choice implies that if a very hard emission (i.e., $`p_\mathrm{T}\sqrt{s}/2`$) occurs at the NLO, only a very limited possibility is left for the MC to emit further, which can lead to problems in the hadronization process. Other choices are of course possible. In the toy model, it can be proven that different $`x_\mathrm{M}`$’s only affect logarithms beyond the accuracy of the MC. In the real case, we have verified that other reasonable choices of EMSCA do not induce visible differences in any of the observables studied.

Our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results are obtained without soft matrix element corrections , i.e. we set the HERWIG parameter SOFTME=.FALSE.. By doing so, we obtain the $`𝒪(\alpha _\mathrm{s})`$ HERWIG result reported in app. A.5, which is crucial in our implementation of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ to avoid double counting. Furthermore, the purpose of soft matrix element corrections is to ensure smooth matching on the boundary between the regions dominated by parton showers and by hard emissions, whereas in our formalism there is no such hard boundary. In order to have meaningful comparisons with standard MC results, we set SOFTME=.FALSE. in that case as well. We verified that the choice of SOFTME has in general only a small impact on the MC results for the observables we studied; the only visible effects are in the shapes of distributions in those regions in which a resummation is needed: $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}0`$, $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}\pi `$, and $`p_\mathrm{T}^{(jet)}0`$ for the hardest jet of the event. A non-negligible effect ($`10`$%) is also seen at $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}0`$, which is partly related to the behaviour at $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}\pi `$, since the total rate is independent of the value of SOFTME.

We should mention finally that we performed many studies of the stability of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results. We studied the dependences of the physical results upon the free parameters $`\alpha `$ and $`\beta `$ introduced in eqs. (LABEL:Gfundef) and (195), and found them completely negligible. This is a rather direct confirmation of the claim in app. A.5, that the angular distribution of soft partons is irrelevant in the determination of physical observables. We performed test $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ runs with increased values of the showering cutoffs VQCUT and VGCUT, and found no significant variations with respect to the results presented here; this supports the claim that MC and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ have the same dependences upon MC showering cutoffs. Also, we find that the parton densities used for MC backward evolution do not strongly influence the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results for this production process. The results obtained using low-$`\alpha _\mathrm{s}`$ MRST99 in the MC evolution are statistically compatible with those obtained using MRST-LO (low-$`\alpha _\mathrm{s}`$ MRST99 being used in both cases for the NLO part).

## 5 Other processes

The extension of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ to other processes, at least those with unambiguous colour flow, seems straightforward in principle. The key equation, (76), is formally valid for any kind of singularities. For final-state collinear singularities, the choice of the Bjorken $`x`$’s and the flavour assignment are trivial: one just has to give the MC the lower-order process that preceded the collinear splitting.

In processes where there are several possible colour flows, the procedure adopted in the HERWIG MC is first to select the kinematics according to the full differential cross section, including interference between colour flows, and then to assign a colour flow according to the relative contributions to the cross section neglecting interference (i.e. in the large-$`N_c`$ limit) . This procedure can easily be extended to the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ case. We point out that the colour flows can be determined independently for $`𝕊`$ and $``$ events. The correct NLO results will still be recovered upon expansion in $`\alpha _\mathrm{s}`$, since the cancellation of the MC terms (75) inserted in the modified subtraction formula, eq. (76), occurs on statistical basis, and not event by event. In order to construct the MC-like cross sections (74), the technique presented in this paper can be adopted: namely, the pure MC result is smoothly matched onto that predicted by perturbative QCD in the soft limit. In fact, independently of the choice of the showering variables, the former cross section is in general locally different from the latter, since it is only in the former that the radiation function is integrated over an azimuthal angle, in order to implement colour coherence in a MC-friendly manner. For processes involving emission from gluons, there are also azimuthal correlations in the collinear limit, which integrate to zero and are therefore not included in MC’s but cannot be neglected locally.

For some processes, such as QCD jet production, the total cross section $`\sigma _{tot}`$ is not defined without a cut; for example one may require the total amount of transverse energy to be bigger than a given quantity. The formulation of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ described in this paper should remain valid in the presence of such a cut.

In processes involving incoming photons, pointlike and hadronic components needs to be defined in a given subtraction scheme. Within the NLO approximation, a change of scheme gives rise to terms that are of higher order. However, for consistency the same scheme should be adopted for the generation of the NLO cross section and the backward evolution of the parton shower in the MC. We also point out that, in order for the procedure described in this paper to be applicable, the photons should have a continuous energy spectrum, such as that from an incoming electron beam, in which the photon energies are distributed according to the Weizsäcker-Williams function.

## 6 Conclusions and outlook

We have presented a practical proposal for constructing an event generator which, we believe, displays the main advantages of the NLO and MC approaches. For brevity we call this $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. A key ingredient is the use of a modified subtraction method for dealing with infrared and collinear singularities. The method has been illustrated using a toy model and implemented in full for the process of hadroproduction of W boson pairs. The results there look encouraging (although of course comparisons with data are for the future): all the observables investigated show the expected behaviour in regions where the NLO or MC results should be most reliable, with a smooth transition between hard and soft/collinear regions.

The method generates a fraction of events with negative weights, but in the cases we have studied the fraction is small and unweighted (signed) event generation remains rather efficient. One can reduce the number of negative weights for a given process by tuning the parameters specific to that process (these parameters have nothing to do with the global MC parameters, which are tuned once and for all by comparison to data). We shall consider in future upgrades of the implementation of our method the strategies outlined in sect. 4.5 and in sect. 4.6, which should reduce the number of $`𝕊`$ and $``$ negative-weight events respectively, using an intermediate integration step, and a modification of the showering algorithm. One could also consider more exotic possibilities, such as defining “superevents” by combining configurations that are close together in three-body phase space.

The implementation presented here is based on the HERWIG MC generator, but there appears no reason why other parton shower generators such as PYTHIA should not be used. The main points to be reconsidered when changing the MC generator will be the dead zones for hard emission, the matching of local counterterms in the soft region, and the relation between the 2- and 3-body configurations connected by the backward evolution in the showering initiated by initial-state partons (corresponding to the event projection that we used in this paper).

###### Acknowledgments.

We wish to thank CERN TH division for hospitality during the preparation of this work. Many informative discussions with S. Catani, G. Corcella, M. Dobbs, J.-Ph. Guillet, Z. Kunszt, M. Mangano, P. Nason, R. Pittau, G. Ridolfi, and M. Seymour are gratefully acknowledged. We also thank S. Catani, J. Collins, and M. Dobbs for comments on the manuscript.

## Appendix A W+W- cross sections

In this appendix, we give all the details relevant to cross sections at $`𝒪(\alpha _\mathrm{s}^0)`$ and $`𝒪(\alpha _\mathrm{s})`$. Although we study explicitly the case of W+W- production, some of the results are based on fairly general ideas, which can be applied to any other hard production process.

### A.1 Kinematics

We start by defining some kinematic variables for W+W- production. The momenta in the $`23`$ partonic processes are assigned as follows:

$$a(p_1)+b(p_2)c(k)+\mathrm{W}^+\mathrm{W}^{}(k_{\mathrm{WW}}).$$
(95)
Following ref. , we introduce the quantities

$$x=M_{\mathrm{WW}}^2/s,y=\mathrm{cos}\theta ,$$
(96)
where $`s=(p_1+p_2)^2`$, $`M_{\mathrm{WW}}`$ is the invariant mass of the W pair, and $`\theta `$ is the scattering angle of the outgoing light parton in the c.m. frame of the colliding partons. In this frame, $`1x`$ is the energy fraction of the outgoing light parton (in units of $`\frac{1}{2}\sqrt{s}`$). These quantities unambiguously identify the soft and initial-state collinear regions of the three-body phase space. In particular

$`x=1\mathrm{soft}(k^0=0),`$ (97)
$`y=1\mathrm{collinear}(\stackrel{}{k}\stackrel{}{p}_1),`$ (98)
$`y=1\mathrm{collinear}(\stackrel{}{k}\stackrel{}{p}_2).`$ (99)
The ranges allowed for $`x`$ and $`y`$ are

$$\rho x1,1y1,$$
(100)
where $`\rho =4m_\mathrm{W}^2/s`$. Defining

$$\beta =\sqrt{1\rho },\beta _x=\sqrt{1\rho /x},$$
(101)
the two-body and three-body phase spaces are

$`d\varphi _2(s)`$ $`=`$ $`{\displaystyle \frac{\beta }{16\pi }}d\mathrm{cos}\theta _\mathrm{W},`$ (102)
$`d\varphi _2(xs)`$ $`=`$ $`{\displaystyle \frac{\beta _x}{16\pi }}d\mathrm{cos}\theta _\mathrm{W},`$ (103)
$`d\varphi _3(s)`$ $`=`$ $`{\displaystyle \frac{s\beta _x}{512\pi ^4}}(1x)dxdyd\mathrm{cos}\theta _\mathrm{W}d\phi ,`$ (104)
where $`\theta _\mathrm{W}`$ and $`\phi `$ are the polar and azimuthal angles of the W+ in the W+W- c.m. frame. The normalization in eq. (104) is such that both $`\theta _\mathrm{W}`$ and $`\phi `$ range between 0 and $`\pi `$.

### A.2 Standard subtraction

Here we give the subtraction terms using the definition of ref. . The quantity $`𝒬_{ab}`$ appearing in eq. (61) reads

$$𝒬_{q\overline{q}}=\frac{C_\mathrm{F}\alpha _\mathrm{s}}{4\pi }\left(\mathrm{log}\frac{s}{\mu ^2}(6+16\mathrm{log}\stackrel{~}{\beta })+32\mathrm{log}^2\stackrel{~}{\beta }\frac{4}{3}\pi ^2\right),$$
(105)
and $`𝒬_{ab}=0`$ for $`\{a,b\}\{q,\overline{q}\}`$. We have introduced a free parameter $`\stackrel{~}{\rho }`$ such that

$$\rho \stackrel{~}{\rho }<1,\stackrel{~}{\beta }=\sqrt{1\stackrel{~}{\rho }}.$$
(106)
The collinear terms of eq. (62) read

$`d\widehat{\sigma }_{ab}^{(c+)}(p_1,p_2)`$ $`=`$ $`{\displaystyle \frac{\alpha _\mathrm{s}}{2\pi }}\{[\mathrm{log}{\displaystyle \frac{\omega s}{2\mu ^2}}\left({\displaystyle \frac{1}{1x}}\right)_{\stackrel{~}{\rho }}+2\left({\displaystyle \frac{\mathrm{log}(1x)}{1x}}\right)_{\stackrel{~}{\rho }}](1x)P_{ca}^{(0)}(x)`$ (107)
$`P_{ca}^{(1)}(x)\}^{(b)}_{cb}(xp_1,p_2)d\varphi _2(xs)dx,`$
and analogously for $`d\widehat{\sigma }_{ab}^{(c)}`$. The real emission contribution of eq. (65) reads

$`d\widehat{\sigma }_{ab}^{(f)}(p_1,p_2)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{1x}}\right)_{\stackrel{~}{\rho }}\left[\left({\displaystyle \frac{1}{1y}}\right)_\omega +\left({\displaystyle \frac{1}{1+y}}\right)_\omega \right]`$ (108)
$`\times (1x)(1y^2)^{(r)}(p_1,p_2)d\varphi _3(s).`$
The distributions introduced in eqs. (107) and (108) define the subtraction terms as follows:

$`{\displaystyle _\rho ^1}𝑑xh(x)\left({\displaystyle \frac{1}{1x}}\right)_{\stackrel{~}{\rho }}`$ $`=`$ $`{\displaystyle _\rho ^1}𝑑x{\displaystyle \frac{h(x)\mathrm{\Theta }(x\stackrel{~}{\rho })h(1)}{1x}},`$ (109)
$`{\displaystyle _\rho ^1}𝑑xh(x)\left({\displaystyle \frac{\mathrm{log}(1x)}{1x}}\right)_{\stackrel{~}{\rho }}`$ $`=`$ $`{\displaystyle _\rho ^1}𝑑x{\displaystyle \frac{h(x)\mathrm{\Theta }(x\stackrel{~}{\rho })h(1)}{1x}}\mathrm{log}(1x),`$ (110)
$`{\displaystyle _1^1}𝑑yg(y)\left({\displaystyle \frac{1}{1\pm y}}\right)_\omega `$ $`=`$ $`{\displaystyle _1^1}𝑑y{\displaystyle \frac{g(y)\mathrm{\Theta }(y1+\omega )g(1)}{1\pm y}},`$ (111)
where $`0<\omega 2`$, and $`h`$ and $`g`$ are generic regular functions. It is clear that a variation of the (arbitrary) parameters $`\stackrel{~}{\rho }`$ and $`\omega `$ amounts to a finite redefinition of the soft and collinear counterterms respectively. The cross section is strictly independent of the values of $`\stackrel{~}{\rho }`$ and $`\omega `$ in their whole range (i.e., these parameters do not need to be small, as opposed to the slicing method).

### A.3 $`𝜻`$–subtraction

We give here an alternative form of the subtraction terms, which has been used in our implementation of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. Although the results of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ do not depend on the particular form adopted for the subtraction terms, we find that the form introduced in this section gives a superior numerical stability, in a sense that, in the context of our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, allows to reduce the number of negative-weight events. See sect. 4.5 for more details. We define

$$𝒫(x,y)=(1x)^2(1y^2),$$
(112)
where $`x`$ and $`y`$ are defined in eq. (96). The function $`𝒫`$ has a simple physical meaning; in fact

$$p_\mathrm{T}^{\left(\mathrm{WW}\right)}=\frac{\sqrt{s}}{2}\sqrt{𝒫(x,y)},$$
(113)
where $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$ is the transverse momentum of the W pair, which coincides with that of the outgoing light parton at NLO. The $`\zeta `$–subtraction basically amounts to subtracting the counterterms when the condition

$$𝒫(x,y)<\zeta $$
(114)
is met, $`\zeta `$ being a numerical parameter, $`0<\zeta 1`$. More precisely, we define

$`{\displaystyle _\rho ^1}𝑑xh(x)\left({\displaystyle \frac{1}{1x}}\right)_𝒫`$ $`=`$ $`{\displaystyle _\rho ^1}𝑑x{\displaystyle \frac{h(x)\mathrm{\Theta }(\zeta 𝒫(x,y))h(1)}{1x}},`$ (115)
$`{\displaystyle _\rho ^1}𝑑xh(x)\left({\displaystyle \frac{F(x)}{1x}}\right)_𝒫`$ $`=`$ $`{\displaystyle _\rho ^1}𝑑x{\displaystyle \frac{h(x)\mathrm{\Theta }(\zeta 𝒫(x,y))h(1)}{1x}}F(x),`$ (116)
$`{\displaystyle _1^1}𝑑yg(y)\left({\displaystyle \frac{1}{1\pm y}}\right)_𝒫`$ $`=`$ $`{\displaystyle _1^1}𝑑y{\displaystyle \frac{g(y)\mathrm{\Theta }(\zeta 𝒫(x,y))\mathrm{\Theta }(y)g(1)}{1\pm y}},`$ (117)
where $`y`$ in eqs. (115), (116), and $`x`$ in eq. (117), are held fixed. We have

$`𝒬_{q\overline{q}}={\displaystyle \frac{C_\mathrm{F}\alpha _\mathrm{s}}{4\pi }}\left(\mathrm{log}{\displaystyle \frac{s}{\mu ^2}}(6+16\mathrm{log}\beta )+32\mathrm{log}^2\beta {\displaystyle \frac{4}{3}}\pi ^2+_s\mathrm{\Theta }(\overline{y})\right),`$ (118)
$`d\widehat{\sigma }_{ab}^{(c+)}(p_1,p_2)={\displaystyle \frac{\alpha _\mathrm{s}}{2\pi }}\{[\mathrm{log}{\displaystyle \frac{s}{2\mu ^2}}\left({\displaystyle \frac{1}{1x}}\right)_𝒫+\left({\displaystyle \frac{\mathrm{log}(1_c(x))}{1x}}\right)_𝒫`$
$`+2\left({\displaystyle \frac{\mathrm{log}(1x)}{1x}}\right)_𝒫](1x)P_{ca}^{(0)}(x)P_{ca}^{(1)}(x)\}^{(b)}_{cb}(xp_1,p_2)d\varphi _2(xs)dx,`$ (119)
$`d\widehat{\sigma }_{ab}^{(f)}(p_1,p_2)={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{1x}}\right)_𝒫\left[\left({\displaystyle \frac{1}{1y}}\right)_𝒫+\left({\displaystyle \frac{1}{1+y}}\right)_𝒫\right]`$
$`\times (1x)(1y^2)^{(r)}(p_1,p_2)d\varphi _3(s),`$ (120)
where

$`\overline{y}=\{\begin{array}{cc}\sqrt{1\zeta /\beta ^4}\hfill & \sqrt{\zeta }<1\rho ,\hfill \\ 0\hfill & \sqrt{\zeta }1\rho ,\hfill \end{array}`$ (123)
and

$`_s`$ $`=`$ $`4\mathrm{log}{\displaystyle \frac{\zeta }{\beta ^4}}\mathrm{log}{\displaystyle \frac{1+\overline{y}}{1\overline{y}}}2\mathrm{log}{\displaystyle \frac{1+\overline{y}}{1\overline{y}}}\mathrm{log}(4(1\overline{y}^2))`$ (124)
$`+4\mathrm{L}\mathrm{i}_2\left({\displaystyle \frac{1+\overline{y}}{2}}\right)4\mathrm{L}\mathrm{i}_2\left({\displaystyle \frac{1\overline{y}}{2}}\right),`$
$`_c(x)`$ $`=`$ $`\{\begin{array}{cc}\sqrt{1\zeta /(1x)^2}\hfill & x<1\sqrt{\zeta },\hfill \\ 0\hfill & x1\sqrt{\zeta }.\hfill \end{array}`$ (127)
We point out that the subtraction scheme defined in this appendix is rather general. In fact, the condition that determines the subtraction, eq. (114), applies to the variables of the parton that recoils against the system, regardless of the nature of the system itself.

We can now go back to the point raised in sect. 4.5, concerning the number of negative-weight events in our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. As mention there, this number is a (very complicated) function of the free parameters used to define the subtractions terms: $`\stackrel{~}{\rho }`$ and $`\omega `$ for standard subtraction, $`\zeta `$ for $`\zeta `$–subtraction. However, for a crude estimate of their best values, we can assume that the partonic cross sections do not depend upon them, and thus that the parameter dependence of $`J_𝕊`$ is only due to the size of the phase-space regions in which the subtraction of the counterterms is performed (this amounts to saying that the only parameter dependence taken into account is that due to the $`\mathrm{\Theta }`$ functions in eqs. (109)–(111) for the standard subtraction, and in eqs. (115)–(117) for the $`\zeta `$–subtraction). In general, there is a cancellation between the first two terms in eq. (78); to a certain extent, the term $`d\overline{\mathrm{\Sigma }}_{ab}^{(f)}/d\varphi _3|_{\mathrm{ct}}`$ can be tuned in order to cancel with $`d\overline{\mathrm{\Sigma }}_{ab}/d\varphi _3|_{\mathrm{MC}}`$. This explains why $`\zeta `$-subtraction is a better choice than standard subtraction in the present case: the region in which the counterterms are subtracted, eq. (114), is a better match to the dead zone, eq. (174), compared with the corresponding region for the standard subtraction.

In this paper we do not try to find an optimal value for the parameter $`\zeta `$. The crude estimate mentioned above gives $`\zeta =0.275`$; the fraction of negative-weight events that we get by using this value is 12%. This fraction appears to be remarkably stable against variations of $`\zeta `$ in the region $`0.05\zeta 0.5`$.

### A.4 Choice of Bjorken $`𝒙`$’s

Here we give the form of the change of variables in eq. (71). We stress that this is a convenient procedure for getting rid of the $`\stackrel{~}{\mathrm{𝟐}}`$ configurations that appear in eq. (68), which also has the virtue of being useful for defining the MC-like terms of eq. (75), as we shall show in app. A.5. These two problems could in principle be solved independently, and by different means, as it may be necessary in other MC implementations; however, our master formula, eq. (76), would still hold. We denote by

$$\mathrm{𝟑}(x_1,x_2,\varphi _3)$$
(128)
the three-body configuration that is obtained by using the three-body phase-space variables $`\varphi _3`$, and the momentum fractions $`x_1`$ and $`x_2`$; the momenta in eq. (128) are understood in the hadronic c.m. frame. Furthermore, we denote by

$$\mathrm{𝟐}_s(x_1,x_2,\varphi _3^{(s)})\mathrm{𝟑}(x_1,x_2,\varphi _3^{(s)})$$
(129)
the two-body configuration that is obtained by setting to zero the energy of the final-state light parton; $`\varphi _3^{(s)}`$ means “soft limit” of $`\varphi _3`$. In the parametrization given in eq. (104), this limit can be obtained by letting $`x1`$ (see eq. (97)). It has to be pointed out that any $`\mathrm{𝟐}`$ configuration appearing in eq. (68) can be viewed as the soft limit of a $`\mathrm{𝟑}`$ configuration; thus, eq. (129) is fully general. Analogously, we denote by

$`\stackrel{~}{\mathrm{𝟐}}_{c+}(x_1,x_2,\varphi _3^{(c+)})`$ $``$ $`\mathrm{𝟑}(x_1,x_2,\varphi _3^{(c+)})`$ (130)
$`\stackrel{~}{\mathrm{𝟐}}_c(x_1,x_2,\varphi _3^{(c)})`$ $``$ $`\mathrm{𝟑}(x_1,x_2,\varphi _3^{(c)})`$ (131)
the three-body configurations that are obtained when the final-state light parton becomes collinear to the incoming parton coming from the left ($`p_1k`$) and from the right ($`p_2k`$) respectively; these limits correspond to $`y1`$ and $`y1`$ respectively (see eqs. (98) and (99)).

We now look for changes of variables of the type given in eq. (71). There, the possibility exists to have a different form for each of the terms contributing to eq. (68). However, since we are looking for a solution to the equation $`\mathrm{𝟐}(x_i)=\stackrel{~}{\mathrm{𝟐}}(\stackrel{~}{x}_i)`$, we shall actually have to find a change of variables (71) for each kinematical configuration $`\mathrm{𝟐}`$ and $`\stackrel{~}{\mathrm{𝟐}}`$. Since only kinematics matters here, it follows that the collinear counterterms to $`d\mathrm{\Sigma }_{ab}^{(f)}|_{\mathrm{ev}}`$, and the terms $`d\mathrm{\Sigma }_{ab}^{(c\pm )}|_{\mathrm{ev}}`$, must be treated exactly in the same way in eq. (68), since they all have $`\stackrel{~}{\mathrm{𝟐}}`$ kinematics. Similarly, the Born term $`d\overline{\mathrm{\Sigma }}_{ab}^{(b)}`$, the soft-virtual term $`d\overline{\mathrm{\Sigma }}_{ab}^{(sv)}`$, the soft counterterm to $`d\mathrm{\Sigma }_{ab}^{(f)}|_{\mathrm{ev}}`$, and the counterterms $`d\mathrm{\Sigma }_{ab}^{(c\pm )}|_{\mathrm{ct}}`$, are treated in the same way, since they have $`\mathrm{𝟐}`$ kinematics. In summary, our task amounts to finding the three changes of variables

$$x_1=x_1^{(\alpha )}(z_1,z_2,\varphi _3),x_2=x_2^{(\alpha )}(z_1,z_2,\varphi _3),\alpha =s,c+,c,$$
(132)
such that, at fixed $`z_1`$, $`z_2`$ and $`\varphi _3`$,

$`\mathrm{𝟐}_s(x_1^{(s)}(z_1,z_2,\varphi _3),x_2^{(s)}(z_1,z_2,\varphi _3),\varphi _3^{(s)})`$
$`=\stackrel{~}{\mathrm{𝟐}}_{c+}(x_1^{(c+)}(z_1,z_2,\varphi _3),x_2^{(c+)}(z_1,z_2,\varphi _3),\varphi _3^{(c+)})`$
$`=\stackrel{~}{\mathrm{𝟐}}_c(x_1^{(c)}(z_1,z_2,\varphi _3),x_2^{(c)}(z_1,z_2,\varphi _3),\varphi _3^{(c)}),`$ (133)
where the equality holds for the W+ and W- variables only. It is then apparent that eqs. (68) and (73) are identical; in fact, if $`\mathrm{𝟐}=\stackrel{~}{\mathrm{𝟐}}`$ for the W+ and W- variables, then $`O(\mathrm{𝟐})=O(\stackrel{~}{\mathrm{𝟐}})`$, since the outgoing light parton cannot contribute to the kinematics of any infrared-safe observables, being either soft or collinear to the beam line.

The requirement $`\mathrm{𝟐}=\stackrel{~}{\mathrm{𝟐}}`$ amounts to two conditions, since in both configurations the transverse momentum of the W pair is equal to zero, and only two components of the four-momentum of the pair are left to be constrained. Notice that it is sufficient to constrain the pair momentum, since the internal variables of the pair are held fixed in eq. (133). Instead of imposing eq. (133) directly, it turns out to be convenient to follow the strategy developed in ref.  (where the procedure described here is called “event projection”), where $`\mathrm{𝟐}`$ and $`\stackrel{~}{\mathrm{𝟐}}`$ are related to $`\mathrm{𝟑}`$. This is done by imposing

$`O_1(\mathrm{𝟑}(z_1,z_2,\varphi _3))`$ $`=`$ $`O_1(\mathrm{𝟐}_s(x_1^{(s)},x_2^{(s)},\varphi _3^{(s)})),`$ (134)
$`O_2(\mathrm{𝟑}(z_1,z_2,\varphi _3))`$ $`=`$ $`O_2(\mathrm{𝟐}_s(x_1^{(s)},x_2^{(s)},\varphi _3^{(s)})),`$ (135)
and analogous equations for $`c+`$ and $`c`$; $`O_1`$ and $`O_2`$ are two arbitrary observables, defined with the variables of the W+W- pair; therefore, the conditions in eq. (133) are effectively implemented by relating $`\mathrm{𝟐}`$ and $`\stackrel{~}{\mathrm{𝟐}}`$ to $`\mathrm{𝟑}`$. It must be stressed that eqs. (134) and (135) are imposed in the hadronic c.m. frame for those observables that are not invariant under longitudinal boosts.

Notice that this procedure can be applied to any production process; in the general case, the observables $`O_1`$ and $`O_2`$ are defined with the variables of the system that recoils against the hard parton emitted by real diagrams at the NLO level. In the present case, one obvious choice is

$$O_1=M_{\mathrm{WW}}^2,O_2=y_{\mathrm{WW}},$$
(136)
where $`M_{\mathrm{WW}}^2`$ and $`y_{\mathrm{WW}}`$ are the invariant mass squared and the rapidity of the pair respectively. These two observables are easily expressed in terms of the momentum fractions $`x_{1,2}`$, and of the phase space variables in the parametrization given in eq. (104):

$`M_{\mathrm{WW}}^2(\mathrm{𝟑}(z_1,z_2,\varphi _3))M_{\mathrm{WW}}^2(z_1,z_2,x)=xz_1z_2S,`$ (137)
$`y_{\mathrm{WW}}(\mathrm{𝟑}(z_1,z_2,\varphi _3))y_{\mathrm{WW}}(z_1,z_2,x,y)=\mathrm{log}\xi +{\displaystyle \frac{1}{2}}\mathrm{log}{\displaystyle \frac{z_1}{z_2}},`$ (138)
where $`S`$ is the hadronic c.m. energy squared, and

$$\xi =\sqrt{\frac{2(1x)(1+y)}{2(1x)(1y)}}.$$
(139)
By imposing the conditions in eqs. (134) and (135), and using the explicit forms for $`M_{\mathrm{WW}}^2`$ and $`y_{\mathrm{WW}}`$ given in eqs. (137) and (138), we get the following:

$`x_1^{(s)}=z_1\xi \sqrt{x},x_2^{(s)}=z_2{\displaystyle \frac{\sqrt{x}}{\xi }},`$ (140)
$`x_1^{(c+)}=x_1^{(s)}/x,x_2^{(c+)}=x_2^{(s)},`$ (141)
$`x_1^{(c)}=x_1^{(s)},x_2^{(c)}=x_2^{(s)}/x.`$ (142)
It should be stressed that $`x_1^{(\alpha )}`$ and $`x_2^{(\alpha )}`$ as given in eqs. (140), (141), and (142) do not need to be less than 1 any longer; thus, the parton densities contained in $`\overline{\mathrm{\Sigma }}_{ab}`$ have to be defined as equal to zero if their arguments are larger than 1.

Although the choice made in eq. (136) is arbitrary, we shall show later that, in the context of our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, it is actually dictated by the procedure used by the MC to generate the Bjorken $`x`$’s. In this procedure, the MC typically constrains either the invariant mass and the rapidity of the pair, or the invariant mass and the longitudinal momentum of the pair. In the current version of HERWIG the latter choice is adopted. We therefore also consider

$$O_1=M_{\mathrm{WW}}^2,O_2=k_{\mathrm{WW}},$$
(143)
where $`k_{\mathrm{WW}}`$ is the longitudinal momentum of the pair. With the parametrization of eq. (104), we get

$$k_{\mathrm{WW}}(\mathrm{𝟑}(z_1,z_2,\varphi _3))k_{\mathrm{WW}}(z_1,z_2,x,y)=\frac{\sqrt{S}}{4}\left[(1x)y(z_1+z_2)+(1+x)(z_2z_1)\right].$$
(144)
By imposing again eqs. (134) and (135), we now obtain

$`x_1^{(s)}={\displaystyle \frac{1}{2}}\left(\sqrt{B^2+4A}B\right),x_2^{(s)}={\displaystyle \frac{1}{2}}\left(\sqrt{B^2+4A}+B\right),`$ (145)
$`x_1^{(c+)}=x_1^{(s)}/x,x_2^{(c+)}=x_2^{(s)},`$ (146)
$`x_1^{(c)}=x_1^{(s)},x_2^{(c)}=x_2^{(s)}/x,`$ (147)
where

$`A`$ $`=`$ $`xz_1z_2,`$ (148)
$`B`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[(1x)y(z_1+z_2)+(1+x)(z_2z_1)\right].`$ (149)
Note that eq. (146) is identical to eq. (141), and eq. (147) is identical to eq. (142). This gives a consistent picture when the results given above are used in the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, eq. (76). In fact, event projection has been introduced in the context of an NLO computation, and the constraints have been imposed on final-state quantities only. However, in an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ the initial-state configuration of the hard process is relevant as well, since it defines the initial conditions for the MC evolutions implicit in $`I_{\mathrm{MC}}(O,\mathrm{𝟑})`$ and $`I_{\mathrm{MC}}(O,\mathrm{𝟐})`$. Eqs. (140), (141), and (142) \[and eqs. (145), (146), and (147)\] all enter the definition of $`𝕊`$ events, and since in general $`x_i^{(s)}x_i^{(c+)}x_i^{(c)}`$, there is an ambiguity for the prescription of the momenta of the initial state partons that initiate the backward evolution. However, as explained in sect. 4.5, the initial conditions for class $`𝕊`$ events are the $`22`$ hard processes; the initial states of these processes coincide, by definition, with those of the $`22`$ contributions to $`𝕊`$ events, but do not coincide in general with those of the quasi-$`22`$ contributions. If we denote the momentum fractions of the initial-state partons entering the $`22`$ processes as $`x_1^{}`$ and $`x_2^{}`$, owing to the fact that $`1x`$ is the energy fraction of the outgoing light parton (see app. A.1), we have

$`x_1^{(s)^{}}=x_1^{(s)},x_2^{(s)^{}}=x_2^{(s)},`$ (150)
$`x_1^{(c+)^{}}xx_1^{(c+)}=x_1^{(s)},x_2^{(c+)^{}}x_2^{(c+)}=x_2^{(s)},`$ (151)
$`x_1^{(c)^{}}x_1^{(c)}=x_1^{(s)},x_2^{(c)^{}}xx_2^{(c)}=x_2^{(s)}.`$ (152)
Thus, $`x_i^{(s)^{}}=x_i^{(c+)^{}}=x_i^{(c)^{}}`$, and there is no ambiguity in the momenta of the initial state partons. It should be noted that an ambiguity here would in any case affect our results at $`𝒪(\alpha _\mathrm{s}^2)`$, since it would only be due to quasi-$`22`$ contributions; however, the absence of any ambiguity fits nicely in a picture where MC and NLO give different descriptions of the same physics.

In appendix A.5, we shall need to consider the soft and collinear limits of $`x_i^{(s)}`$ and $`x_i^{(c\pm )}`$. By direct computation, we obtain the following

$`x1`$ $`\{\begin{array}{cc}x_1^{(s)}z_1,\hfill & x_2^{(s)}z_2,\hfill \\ x_1^{(c+)}z_1,\hfill & x_2^{(c+)}z_2,\hfill \\ x_1^{(c)}z_1,\hfill & x_2^{(c)}z_2,\hfill \end{array}`$ (156)
$`y1`$ $`\{\begin{array}{cc}x_1^{(s)}xz_1,\hfill & x_2^{(s)}z_2,\hfill \\ x_1^{(c+)}z_1,\hfill & x_2^{(c+)}z_2,\hfill \end{array}`$ (159)
$`y1`$ $`\{\begin{array}{cc}x_1^{(s)}z_1,\hfill & x_2^{(s)}xz_2,\hfill \\ x_1^{(c)}z_1,\hfill & x_2^{(c)}z_2,\hfill \end{array}`$ (162)
where we did not consider the $`y1`$ limit of $`x_i^{(c)}`$ and the $`y1`$ limit of $`x_i^{(c+)}`$, since these limits are not relevant in the computation of singular cross sections. We stress that eqs. (156)–(162) hold regardless of whether eq. (136) or eq. (143) has been used. This is a manifestation of the factorization properties of QCD, as will become clear in app. A.5.

### A.5 Dead zone and subtraction term

Here we define the quantity $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$, introduced in eq. (74). We make use of some of the results of ref. ; although ref.  concerns only single vector boson production, the kinematics carry over to the present case if we simply replace the momentum of the vector boson by that of the pair. According to the discussion of sect. 4.1, we have to consider the cross section that we would obtain by keeping the $`𝒪(\alpha _\mathrm{s})`$ term in the expansion of the result of the MC, which we denote by $`d\sigma |_{\mathrm{MC}}`$. Roughly speaking, $`d\sigma |_{\mathrm{MC}}`$ is equal to the Born cross section, times the $`𝒪(\alpha _\mathrm{s})`$ term in the expansion of the Sudakov form factor. More precisely

$`d\sigma |_{\mathrm{MC}}`$ $`=`$ $`{\displaystyle \underset{abc}{}}dx_1^{(\mathrm{MC})}dx_2^{(\mathrm{MC})}{\displaystyle \frac{\alpha _\mathrm{s}}{2\pi }}d\sigma _{ab}^{(b)}(x_1^{(\mathrm{MC})}P_1,x_2^{(\mathrm{MC})}P_2)`$ (163)
$`\times ({\displaystyle \frac{d\xi _+}{\xi _+}}{\displaystyle \frac{dz_+}{z_+}}\mathrm{\Theta }(z_+^2\xi _+)P_{ac}^{(0)}(z_+)f_c^{(H_1)}(x_1^{(\mathrm{MC})}/z_+)f_b^{(H_2)}(x_2^{(\mathrm{MC})})`$
$`+{\displaystyle \frac{d\xi _{}}{\xi _{}}}{\displaystyle \frac{dz_{}}{z_{}}}\mathrm{\Theta }(z_{}^2\xi _{})P_{bc}^{(0)}(z_{})f_a^{(H_1)}(x_1^{(\mathrm{MC})})f_c^{(H_2)}(x_2^{(\mathrm{MC})}/z_{})).`$
Equation (163) has been obtained directly from the QCD Sudakov form factor; in doing so, we have eliminated any cutoffs, precisely as in the case of the toy model. This implies that eq. (163) has soft and collinear divergences, which is what we want, since from it we derive $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$, namely a local counterterm to the real emission matrix element. The two terms on the r.h.s. account for the emissions from the incoming partons. If strongly-interacting partons were present in the final state at the Born level, corresponding terms should be added. The variables $`\xi _\pm `$ and $`z_\pm `$ are those used by the HERWIG MC  in the parton showering process. The $`\mathrm{\Theta }`$ functions define the region where an emission can occur during the shower; therefore, from eq. (163) we can easily read the form of the dead zone, which is by definition the region where no emission can occur. We have

$$\mathrm{Dead}\mathrm{zone}:\xi _+>z_+^2,\xi _{}>z_{}^2,$$
(164)
where both conditions must be fulfilled. Finally, the $`1/z_\pm `$ factor that appears in the arguments of the parton densities in eq. (163) is due to the backward evolution in the shower from initial-state partons. In the case of emission from final-state partons, the arguments of the parton densities would simply read $`x_1^{(\mathrm{MC})}`$ and $`x_2^{(\mathrm{MC})}`$.

The values of $`x_1^{(\mathrm{MC})}`$ and $`x_2^{(\mathrm{MC})}`$ are determined by the MC (through energy-momentum conservation) at the Born level, after the generation of the $`\mathrm{𝟐}`$ configuration which is used as initial condition for the shower. At the end of the shower, the momenta of the final-state particles are boosted, the boost being determined by imposing a relation between the variables of the W pair (in general, of the system produced at the Born level) before and after the shower. It follows that, after the shower, $`x_1^{(\mathrm{MC})}`$ and $`x_2^{(\mathrm{MC})}`$ have still the values determined by the $`\mathrm{𝟐}`$ configuration; however, after the shower $`x_1^{(\mathrm{MC})}`$ and $`x_2^{(\mathrm{MC})}`$ cannot be interpreted any longer as Bjorken $`x`$’s associated with the incoming partons.

This procedure is identical (when we consider the $`𝒪(\alpha _\mathrm{s})`$ result, that is, the emission of a single parton) to the one we defined for event projection in app. A.4. In other words, the $`\mathrm{𝟑}`$ configuration implicit in the l.h.s. of eq. (163) is related to the $`\mathrm{𝟐}`$ configuration that initiates the shower by means of constraints such as those in eqs. (134) and (135).888This is the reason for introducing event projection in app. A.4, instead of just imposing eq. (133). In a completely general case, one should map the three-body configurations obtained from the MC onto the $`\mathrm{𝟑}`$ configurations of the NLO cross section, thus defining a $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ that acts as a local counterterm in eq. (76). It follows that, if we label the $`\mathrm{𝟑}`$ configurations according to the conventions of this paper, $`\mathrm{𝟑}\mathrm{𝟑}(z_1,z_2,\varphi _3)`$, we have

$$x_1^{(\mathrm{MC})}=x_1^{(s)},x_2^{(\mathrm{MC})}=x_2^{(s)}.$$
(165)
The explicit form of $`x_i^{(s)}`$ depends on the observables $`O_i`$ used by the MC to implement eqs. (134) and (135). A priori, these need not coincide with those given in eq. (136) or eq. (143). In practice, HERWIG uses the variables in eq. (143), and therefore the results in eq. (145) give the explicit form of $`x_1^{(\mathrm{MC})}`$ and $`x_2^{(\mathrm{MC})}`$. The relations between $`x_1^{(\mathrm{MC})}`$, $`x_2^{(\mathrm{MC})}`$ and $`z_1`$, $`z_2`$ allow us to recast eq. (163) in the following form

$$d\sigma |_{\mathrm{MC}}=\underset{ab}{}dz_1dz_2d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}},$$
(166)
where

$`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ $`=`$ $`{\displaystyle \frac{1}{z_+}}{\displaystyle \frac{(x_1^{(s)},x_2^{(s)})}{(z_1,z_2)}}f_a^{(H_1)}(x_1^{(s)}/z_+)f_b^{(H_2)}(x_2^{(s)})d\sigma _{ab}^+|_{\mathrm{MC}}`$ (167)
$`+`$ $`{\displaystyle \frac{1}{z_{}}}{\displaystyle \frac{(x_1^{(s)},x_2^{(s)})}{(z_1,z_2)}}f_a^{(H_1)}(x_1^{(s)})f_b^{(H_2)}(x_2^{(s)}/z_{})d\sigma _{ab}^{}|_{\mathrm{MC}}`$
is the quantity that the appears in eq. (76), and

$`d\sigma _{ab}^+|_{\mathrm{MC}}`$ $`=`$ $`{\displaystyle \frac{\alpha _\mathrm{s}}{2\pi }}{\displaystyle \frac{d\xi _+}{\xi _+}}dz_+\mathrm{\Theta }(z_+^2\xi _+)P_{ca}^{(0)}(z_+)d\sigma _{cb}^{(b)}(x_1^{(s)}P_1,x_2^{(s)}P_2),`$ (168)
$`d\sigma _{ab}^{}|_{\mathrm{MC}}`$ $`=`$ $`{\displaystyle \frac{\alpha _\mathrm{s}}{2\pi }}{\displaystyle \frac{d\xi _{}}{\xi _{}}}dz_{}\mathrm{\Theta }(z_{}^2\xi _{})P_{cb}^{(0)}(z_{})d\sigma _{ac}^{(b)}(x_1^{(s)}P_1,x_2^{(s)}P_2),`$ (169)
where the sum over the parton label $`c`$ is now understood.

In order to prove that $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ in eq. (167) can act as a local counterterm in eq. (76), we have to express it in terms of the variables we use to describe the NLO cross section. This can be easily done as explained in ref. : we relate the variables $`\xi _\pm `$ and $`z_\pm `$ used in the HERWIG parton showering to Lorentz invariant quantities, which we then express in terms of the phase-space variables of eq. (104). Following ref. , we define

$`t_k=(p_1k)^2={\displaystyle \frac{s}{2}}(1x)(1y),`$ (170)
$`u_k=(p_2k)^2={\displaystyle \frac{s}{2}}(1x)(1+y),`$ (171)
where the momenta are assigned as in eq. (95), and $`x`$ and $`y`$ have been defined in eq. (96). We have

$`t_k`$ $`=`$ $`{\displaystyle \frac{1z_+}{z_+^2}}\xi _+M_{\mathrm{WW}}^2,`$ (172)
$`{\displaystyle \frac{t_ku_k}{s}}`$ $`=`$ $`{\displaystyle \frac{(1z_+)^2}{2z_+^2}}\xi _+(2\xi _+)M_{\mathrm{WW}}^2,`$ (173)
for a shower initiated by parton $`a`$. A shower initiated by the other incoming parton ($`b`$) is described in terms of $`\xi _{}`$ and $`z_{}`$, which are in turn related to $`t_k`$ and $`u_k`$ by the same functional forms given in eqs. (172) and (173), with the formal exchange $`t_ku_k`$. We can now substitute eqs. (96), (170), and (171) into eqs. (172) and (173), to find $`\xi _\pm (x,y)`$ and $`z_\pm (x,y)`$, whose explicit form we do not report here, since it is not needed in the following.

We can now find the dead zone boundary in the $`x,y`$ plane; from eq. (164) we obtain

$$\mathrm{Dead}\mathrm{zone}:\rho xx_{\mathrm{DZ}},\left|y\right|<Y_{\mathrm{DZ}}(x),$$
(174)
where both conditions are met simultaneously, and

$`x_{\mathrm{DZ}}`$ $`=`$ $`{\displaystyle \frac{7+\sqrt{17}}{16}},`$ (175)
$`Y_{\mathrm{DZ}}(x)`$ $`=`$ $`1{\displaystyle \frac{x}{1x}}\left(3\sqrt{1+8x}\right).`$ (176)
We now consider eqs. (168) and (169). Using eq. (104), we obtain

$`{\displaystyle \frac{d\sigma _{ab}^+(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}`$ $`=`$ $`{\displaystyle \frac{16\pi \alpha _\mathrm{s}}{s}}{\displaystyle \frac{P_{ca}^{(0)}(z_+)}{(1x)\xi _+}}{\displaystyle \frac{(z_+,\xi _+)}{(x,y)}}_{cb}^{(b)}(x_1^{(s)}P_1,x_2^{(s)}P_2)\mathrm{\Theta }(yY_{\mathrm{DZ}}(x)),`$
$`{\displaystyle \frac{d\sigma _{ab}^{}(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}`$ $`=`$ $`{\displaystyle \frac{16\pi \alpha _\mathrm{s}}{s}}{\displaystyle \frac{P_{cb}^{(0)}(z_{})}{(1x)\xi _{}}}{\displaystyle \frac{(z_{},\xi _{})}{(x,y)}}_{ac}^{(b)}(x_1^{(s)}P_1,x_2^{(s)}P_2)\mathrm{\Theta }(yY_{\mathrm{DZ}}(x)),`$
where $`p_i=z_iP_i`$ (according to the definition given in eq. (67)); therefore, $`s=z_1z_2S`$.

We can now consider the soft and collinear limits of $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ given in eq. (167). We observe that, owing to the procedure adopted by the MC to determine $`x_1^{(\mathrm{MC})}`$ and $`x_2^{(\mathrm{MC})}`$, the luminosities that appear in eq. (167) have non-trivial soft and collinear limits. We make use of the fact, obtained by direct computation, that

$$\underset{x1}{lim}z_\pm =1,\underset{y1}{lim}z_+=x,\underset{y1}{lim}z_{}=x.$$
(179)
Therefore, using eqs. (156)–(162), we obtain

$`\underset{x1}{lim}f_a^{(H_1)}(x_1^{(s)}/z_+)f_b^{(H_2)}(x_2^{(s)})=\underset{x1}{lim}f_a^{(H_1)}(x_1^{(s)})f_b^{(H_2)}(x_2^{(s)}/z_{})`$ (180)
$`=`$ $`\underset{y1}{lim}f_a^{(H_1)}(x_1^{(s)}/z_+)f_b^{(H_2)}(x_2^{(s)})=\underset{y1}{lim}f_a^{(H_1)}(x_1^{(s)})f_b^{(H_2)}(x_2^{(s)}/z_{})`$
$`=`$ $`f_a^{(H_1)}(z_1)f_b^{(H_2)}(z_2),`$
where the collinear limits have been considered only for the luminosities attached to short-distance cross sections which are singular in those limits. In the soft and collinear limits, $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ will have to cancel the singularities of $`d\overline{\mathrm{\Sigma }}_{ab}^{(f)}`$ and $`d\overline{\mathrm{\Sigma }}_{ab}^{(c\pm )}`$. Owing to the event projection procedure described in app. A.4, the luminosities attached to the latter two terms behave as follows in the soft and collinear limits:

$`\underset{x1}{lim}f_a^{(H_1)}(x_1^{(s)})f_b^{(H_2)}(x_2^{(s)})=\underset{y1}{lim}f_a^{(H_1)}(x_1^{(c+)})f_b^{(H_2)}(x_2^{(c+)})`$ (181)
$`=`$ $`\underset{y1}{lim}f_a^{(H_1)}(x_1^{(c)})f_b^{(H_2)}(x_2^{(c)})=f_a^{(H_1)}(z_1)f_b^{(H_2)}(z_2),`$
where we used eqs. (156)–(162). Eqs. (180) exactly match eqs. (181). It is also straightforward to prove that

$$\frac{1}{z_\pm }\frac{(x_1^{(s)},x_2^{(s)})}{(z_1,z_2)}1$$
(182)
in the soft and collinear limits, and that the corresponding Jacobians in the NLO cross section have the same behaviour. This implies that the properties of $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ as a local counterterm are determined by the behaviour of its short-distance parts, given in eqs. (LABEL:HWxsecpl) and (LABEL:HWxsecmn).

When considering eqs. (LABEL:HWxsecpl) and (LABEL:HWxsecmn), we notice that these formulae are valid for any parton flavours $`a`$ and $`b`$. However, in our case the soft and collinear limits are all non-trivial only in the case $`a=q`$, $`b=\overline{q}`$, with $`c=g`$; we shall therefore explicitly discuss only this case. We expect the forms of the leading soft and collinear singularities of the cross sections given in eqs. (LABEL:HWxsecpl) and (LABEL:HWxsecmn) to be dictated by QCD factorization theorem. In other words, we expect the following equations to hold

$`{\displaystyle \frac{d\sigma _{q\overline{q}}^+(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}+{\displaystyle \frac{d\sigma _{q\overline{q}}^{}(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}\stackrel{x1}{}8\pi \alpha _\mathrm{s}C_\mathrm{F}{\displaystyle \frac{p_1p_2}{p_1kp_2k}}_{q\overline{q}}^{(b)}(p_1,p_2)`$
$`={\displaystyle \frac{64\pi \alpha _\mathrm{s}}{s}}{\displaystyle \frac{C_\mathrm{F}}{(1x)^2(1y^2)}}_{q\overline{q}}^{(b)}(p_1,p_2),`$ (183)
$`{\displaystyle \frac{d\sigma _{q\overline{q}}^+(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}\stackrel{y1}{}{\displaystyle \frac{4\pi \alpha _\mathrm{s}}{p_1k}}P_{qq}^{(0)}(x)_{q\overline{q}}^{(b)}(xp_1,p_2)`$
$`={\displaystyle \frac{16\pi \alpha _\mathrm{s}}{s}}{\displaystyle \frac{P_{qq}^{(0)}(x)}{(1x)(1y)}}_{q\overline{q}}^{(b)}(xp_1,p_2),`$ (184)
$`{\displaystyle \frac{d\sigma _{q\overline{q}}^{}(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}\stackrel{y1}{}{\displaystyle \frac{4\pi \alpha _\mathrm{s}}{p_2k}}P_{qq}^{(0)}(x)_{q\overline{q}}^{(b)}(p_1,xp_2)`$
$`={\displaystyle \frac{16\pi \alpha _\mathrm{s}}{s}}{\displaystyle \frac{P_{qq}^{(0)}(x)}{(1x)(1+y)}}_{q\overline{q}}^{(b)}(p_1,xp_2).`$ (185)
Note that in the case of the soft limit we sum the short distance cross sections, thanks to eq. (180). Explicit computations of the collinear limits $`y\pm 1`$ in eqs. (LABEL:HWxsecpl) and (LABEL:HWxsecmn) give eqs. (184) and (185). On the other hand, in the soft limit $`x1`$ we find the following:

$`{\displaystyle \frac{d\sigma _{q\overline{q}}^+(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}+{\displaystyle \frac{d\sigma _{q\overline{q}}^{}(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}\stackrel{x1}{}{\displaystyle \frac{64\pi \alpha _\mathrm{s}}{s}}C_\mathrm{F}_{q\overline{q}}^{(b)}(p_1,p_2)`$
$`\times \left({\displaystyle \frac{2\mathrm{\Theta }(y+1/3)}{(1x)^2(1y)(3+y)}}+{\displaystyle \frac{2\mathrm{\Theta }(y+1/3)}{(1x)^2(1+y)(3y)}}\right).`$ (186)
It appears therefore that the choice of the showering variables in HERWIG violates the factorization properties of QCD cross sections. However, the situation is not as bad as it may seem. In fact, the kernel in eq. (186) has the following property:

$$\underset{\epsilon 0}{lim}_{1+\epsilon }^{1\epsilon }𝑑y\left(\frac{2\mathrm{\Theta }(y+1/3)}{(1y)(3+y)}+\frac{2\mathrm{\Theta }(y+1/3)}{(1+y)(3y)}\frac{1}{1y^2}\right)=0.$$
(187)
Eq. (187) implies that, although (in the absence of soft matrix element corrections ) HERWIG generates an incorrect angular distribution of soft gluons, the total amount of radiated soft energy is identical to that computed in perturbative QCD. This is reassuring, and explains why HERWIG correctly describes physical observables, eq. (186) notwithstanding: any IR safe observable is not sensitive to the details of the angular distribution of soft gluons.

Nevertheless, eq. (187) prevents us from using $`d\overline{\mathrm{\Sigma }}|_{\mathrm{MC}}`$ in the construction of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. In fact, $`d\overline{\mathrm{\Sigma }}|_{\mathrm{MC}}`$ acts as a local counterterm in eq. (76), and its soft behaviour (eq. (186)) prevents the necessary cancellation.

This problem can be solved in two different ways. The most straightforward procedure would be to define different showering variables, that would give short-distance cross sections with the correct soft limit. However, this would imply a modification of the existing MC code in one of its crucial parts, the showering algorithm. As already stressed, in this paper we prefer not to modify any substantial part of HERWIG.

The second solution relies on the observation that the angular distribution of soft gluons is irrelevant to $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ physical results. We need a local counterterm to cancel the divergences of the real matrix element, but in the soft emission region the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ cross section is dominated by the pure MC contribution. On the other hand, for hard MC emission, we must have full control of the gluon emitted by the MC, since only in this way will the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ result match the NLO result. It follows that we can choose a form for $`d\mathrm{\Sigma }|_{\mathrm{MC}}`$ identical to that we get by using eqs. (LABEL:HWxsecpl) and (LABEL:HWxsecmn) for non-soft emission, and a form that can act as a local counterterm for real matrix elements in the soft region. One possibility is therefore the following

$`{\displaystyle \frac{d\sigma _{q\overline{q}}^+(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}`$ $`=`$ $`{\displaystyle \frac{16\pi \alpha _\mathrm{s}}{s}}({\displaystyle \frac{P_{qq}^{(0)}(z_+)}{(1x)\xi _+}}{\displaystyle \frac{(z_+,\xi _+)}{(x,y)}}\mathrm{\Theta }(yY_{\mathrm{DZ}}(x))𝒢(x)`$ (188)
$`+{\displaystyle \frac{P_{qq}^{(0)}(x)}{(1x)(1y)}}(1𝒢(x)))^{(b)}_{cb}(x_1^{(s)}P_1,x_2^{(s)}P_2),`$
$`{\displaystyle \frac{d\sigma _{q\overline{q}}^{}(p_1,p_2)}{d\varphi _3}}|_{\mathrm{MC}}`$ $`=`$ $`{\displaystyle \frac{16\pi \alpha _\mathrm{s}}{s}}({\displaystyle \frac{P_{qq}^{(0)}(z_{})}{(1x)\xi _{}}}{\displaystyle \frac{(z_{},\xi _{})}{(x,y)}}\mathrm{\Theta }(yY_{\mathrm{DZ}}(x))𝒢(x)`$ (189)
$`+{\displaystyle \frac{P_{qq}^{(0)}(x)}{(1x)(1+y)}}(1𝒢(x)))^{(b)}_{ac}(x_1^{(s)}P_1,x_2^{(s)}P_2),`$
where $`𝒢`$ is a largely arbitrary function, such that

$$\underset{x1}{lim}𝒢(x)=0,\underset{xx_{\mathrm{DZ}}}{lim}𝒢(x)=1.$$
(190)
A simple choice is

$`𝒢(x)={\displaystyle \frac{(1T(x))^{2\alpha }}{T(x)^{2\alpha }+(1T(x))^{2\alpha }}},T(x)`$ $`=`$ $`\{\begin{array}{cc}(x\stackrel{~}{x}_{\mathrm{DZ}})/(1\stackrel{~}{x}_{\mathrm{DZ}})\hfill & \stackrel{~}{x}_{\mathrm{DZ}}<x<1,\hfill \\ 0\hfill & x<\stackrel{~}{x}_{\mathrm{DZ}},\hfill \end{array}`$ (193)
where

$$\stackrel{~}{x}_{\mathrm{DZ}}=1(1x_{\mathrm{DZ}})\beta .$$
(195)
The parameters in eqs. (LABEL:Gfundef) and (195) must be chosen in the following ranges: $`\alpha 1`$, $`0<\beta 1`$. The dependence of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ results upon these parameters will show the sensitivity of physical observables to the details of soft gluon emission.

We point out that the procedure outlined above is not necessary when the matrix element is not singular in the soft limit. This happens for the partonic subprocesses $`qgW^+W^{}q`$. In such a case, we simply set $`𝒢1`$ in eqs. (188) and (189).

Finally, it is worth noting that the interesting features of eqs. (186) and (187) are of a general nature, since the hard production process is factorized into $`_{q\overline{q}}^{(b)}`$ in eq. (186). Thus, the problem raised here is not specific to W+W- production. We expect that the angular distributions of soft gluons in HERWIG do not agree with those prescribed by the factorization theorems of QCD. Therefore, the local counterterm $`d\overline{\mathrm{\Sigma }}|_{\mathrm{MC}}`$ that we need in order to construct the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ will not coincide in general with the $`𝒪(\alpha _\mathrm{s})`$ HERWIG result. As in the present case, we expect that physical observables will not display any dependence upon the form chosen for $`d\overline{\mathrm{\Sigma }}|_{\mathrm{MC}}`$ in the soft limit. In all cases, we must verify that the total amount of soft gluon energy given by the MC agrees with the QCD prediction, i.e., that the analogue of eq. (187) holds.

## Appendix B Technicalities of $`\mathrm{𝐌𝐂}\mathbf{@}\mathrm{𝐍𝐋𝐎}`$

In this appendix, we present formal proofs of some of the features of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. We do so in the case of the toy model. However, we believe that the corresponding properties hold in the case of QCD as well, as we shall point out explicitly; the toy model just allows us to simplify the notation.

### B.1 Construction of $`𝑰_{\mathrm{𝐌𝐂}}`$

Our first aim is to give a precise definition of $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$. We identify an MC event with the set of energies of the $`n`$ photons emitted in the shower, $`\{x_k\}_{k=1}^n`$; the dependence upon $`x_\mathrm{M}(x)`$ is understood. We define an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ event as the corresponding MC event plus the energy of the photon emitted at the NLO level: $`\{x,x_k\}_{k=1}^n`$. In the case of soft emission, $`x=0`$, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ event simply coincides with the MC event. Any observable $`O`$ will be in general a multi-valued function of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ events, $`O(\{x,x_k\}_{k=1}^n)=\{O_1,\mathrm{},O_p\}`$, with $`1pn+1`$, the $`O_j`$’s being values in the physical range of $`O`$. We then define

$$𝒩(\{x,x_k\}_{k=1}^n;\overline{O},\delta O)=\underset{j=1}{\overset{p}{}}\mathrm{\Theta }(\overline{O}+\delta O/2O_j)\mathrm{\Theta }(O_j\overline{O}+\delta O/2),$$
(196)
where $`\overline{O}`$ is any given value that the observable $`O`$ can assume in its physical range. $`𝒩`$ is then just the number of entries in the bin $`(\overline{O}\delta O/2,\overline{O}+\delta O/2)`$ resulting from the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ event $`\{x,x_k\}_{k=1}^n`$. Therefore $`0𝒩pn+1`$. Starting from a given $`x_\mathrm{M}(x)`$, we generate $`N`$ MC events, and we define

$$I_{\mathrm{MC}}(\overline{O},x_\mathrm{M}(x))=\underset{N\mathrm{}}{lim}\underset{\delta O0}{lim}\frac{1}{N}\frac{1}{\delta O}\underset{i=1}{\overset{N}{}}𝒩(\{x,x_k^{(i)}\}_{k=1}^{n_i};\overline{O},\delta O),$$
(197)
$`n_i`$ being the number of photons, with energies $`x_k^{(i)}`$, emitted in event number $`i`$. Thus in general $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$ is the differential cross section that one obtains by running the MC starting at $`x_\mathrm{M}(x)`$, defining $`O`$ through all photons emitted in the shower, plus that emitted at the NLO level, normalized by the total number of MC events generated.

If we assume that the MC is switched off, then from eq. (197) we obtain

$$I_{\mathrm{MC}}(O,x_\mathrm{M}(x))=\delta (OO(x)),$$
(198)
$`O(x)`$ being the single-valued function used at the NLO level to describe the observable $`O`$. Eq. (198) can be proven by using eq. (197) in the computation of $`𝑑OI_{\mathrm{MC}}(O)f(O)`$, for any arbitrary function $`f(O)`$. With the same technique, we can also find a closed form for $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$ in the general case. We denote by $`P_n(x_1,\mathrm{},x_n)`$ the probability that the toy MC emits $`n`$ photons with energies
$`x_1>x_2>\mathrm{}>x_n`$. This probability can be obtained from eq. (14), and reads as follows

$`P_0`$ $`=`$ $`\mathrm{\Delta }(x_0,x_\mathrm{M}(x)),`$ (199)
$`P_n(x_1,\mathrm{},x_n)`$ $`=`$ $`a^n\mathrm{\Delta }(x_0,x_\mathrm{M}(x))\mathrm{\Theta }(x_\mathrm{M}(x)x_1){\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{Q(x_i)}{x_i}}\mathrm{\Theta }(x_ix_{i+1}),`$ (200)
where in the case of $`n`$ emissions we define $`x_{n+1}=x_0`$. We also denote by $`_n`$ the probability of having $`n`$ emissions, regardless of the photon energies. Clearly, $`_0P_0`$. For $`n>0`$, we can compute this quantity by integrating eq. (200), to obtain

$$_n_0^1\left(\underset{i=1}{\overset{n}{}}dx_i\right)P_n(x_1,\mathrm{},x_n)=\frac{a^n}{n!}\left[_{x_0}^{x_\mathrm{M}(x)}𝑑z\frac{Q(z)}{z}\right]^n\mathrm{\Delta }(x_0,x_\mathrm{M}(x)).$$
(201)
With eqs. (199) and (201) we can check directly that the toy MC conserves probability:

$$_0+\underset{n=1}{\overset{\mathrm{}}{}}_n=1.$$
(202)
We now introduce the functions $`O_{n;j}(x;x_1,\mathrm{},x_n)`$ that return the multi-valued observable $`O`$ (for single-value observables, we simply have $`j1`$). Observing that the sum over all events in eq. (197) can be split into sums over events with a fixed number of emissions, we easily get the following:

$`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))=_0\delta (OO(x))`$
$`+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle _0^1}\left({\displaystyle \underset{i=1}{\overset{n}{}}}dx_i\right)P_n(x_1,\mathrm{},x_n){\displaystyle \underset{j=1}{\overset{p(x;x_1,\mathrm{},x_n)}{}}}\delta (OO_{n;j}(x;x_1,\mathrm{},x_n)),`$ (203)
where the integer $`p(x;x_1,\mathrm{},x_n)`$ in general depends on the kinematics (this happens with jet observables, for example). Eq. (203) can be used, together with the explicit forms of the emission probabilities $`P_n`$, to recover the results for $`I_{\mathrm{MC}}(y)`$ and $`I_{\mathrm{MC}}(z)`$ given in eqs. (25), (26), (38), and (39); these equations have been obtained by direct computation, which constitutes a cross check of eq. (203).

It should be noted that eq. (203) is obtained from eq. (197) without using the explicit form of the emission probabilities $`P_n`$. Thus, both equations hold in the QCD case as well, with obvious formal substitutions (the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ event $`\{x,x_k\}_{k=1}^n`$ gets replaced by an $`𝐧+\mathrm{𝟏}`$ kinematical configuration, and the information on parton flavours must be inserted).

### B.2 Total rate for exclusive observables

We define an observable to be exclusive if its values can be obtained through a single-valued function, regardless of the number of the photons emitted in the shower. The hardest jet energy $`y`$, considered in section 3.4, is an observable of this kind, whereas the fully inclusive photon energy $`z`$ is not. It is clear that in the exclusive case the function $`𝒩`$ defined in app. B.1 has the following properties:

$`𝒩(\{x,x_k\}_{k=1}^n;\overline{O},\delta O)=0,1\overline{O};`$ (204)
$`!\overline{O}\mathrm{for}\delta O0\mathrm{such}\mathrm{that}:𝒩(\{x,x_k\}_{k=1}^n;\overline{O},\delta O)=1.`$ (205)
Eq. (205) holds for any given $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ event. Then

$$𝑑OI_{\mathrm{MC}}(O,x_\mathrm{M}(x))=1.$$
(206)
We can also get this results directly from eq. (203), using eq. (202), and the fact that the observable $`O`$ is a single-valued function (and thus $`p(x;x_1,\mathrm{},x_n)=1`$). It is then apparent that eq. (206) also holds in QCD. We can use eq. (206) in eq. (20); we get

$`{\displaystyle 𝑑O\left(\frac{d\sigma }{dO}\right)_{\mathrm{msub}}}`$ $`=`$ $`{\displaystyle _0^1}dx[\left({\displaystyle }dOI_{\mathrm{MC}}(O,x_\mathrm{M}(x))\right){\displaystyle \frac{a[R(x)BQ(x)]}{x}}`$ (207)
$`+\left({\displaystyle }dOI_{\mathrm{MC}}(O,1)\right)(B+aV+{\displaystyle \frac{aB[Q(x)1]}{x}})]`$
$`=`$ $`{\displaystyle _0^1}𝑑x\left[B+aV+a{\displaystyle \frac{R(x)B}{x}}\right]\sigma _{tot}.`$
In the particular case of $`Oy`$, this result could be obtained from the results for $`d\sigma /dy`$ given in sec. 3.4, by including the full $`x_0`$ dependence that has been neglected there.

Since eq. (206) holds in QCD, we conclude that eq. (207) also holds in QCD. This is in fact confirmed numerically by the results for the exclusive observables (such as $`p_\mathrm{T}^{\mathrm{W}^+}`$, $`p_\mathrm{T}^{\left(\mathrm{WW}\right)}`$, $`\mathrm{\Delta }\varphi ^{\left(\mathrm{WW}\right)}`$) presented in sect. 4.7 (in the case that no cuts are applied).

### B.3 Expansion to $`𝓞\mathbf{(}𝜶_𝐬\mathbf{)}`$

We now consider the expansion to order $`a`$ and $`\alpha _\mathrm{s}`$ of our formulae for the toy model and QCD $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, eqs. (20) and (76) respectively. We start with the toy model. Since we are not interested in $`𝒪(a^2)`$ terms, only zero- and one-photon emissions from the MC are relevant. We have

$$I_{\mathrm{MC}}(O,x_\mathrm{M}(x))=\delta (OO(x))+𝒪(a).$$
(208)
In eq. (208) we do not need to know the $`𝒪(a)`$ term, since $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$ is already multiplied by $`a`$ in eq. (20). On the other hand

$`I_{\mathrm{MC}}(O,1)`$ $`=`$ $`\mathrm{\Delta }(x_0,1)\delta (OO(0))+a{\displaystyle _{x_0}^1}𝑑t{\displaystyle \frac{Q(t)}{t}}\mathrm{\Delta }(x_0,1)\delta (OO(t))+𝒪(a^2)`$
$`=`$ $`\left(1a{\displaystyle _{x_0}^1}𝑑t{\displaystyle \frac{Q(t)}{t}}\right)\delta (OO(0))+a{\displaystyle _{x_0}^1}𝑑t{\displaystyle \frac{Q(t)}{t}}\delta (OO(t))+𝒪(a^2).`$
The first term on the r.h.s. of eq. (B.3) is the contribution to the observable in the case that the MC does not emit any photons. The corresponding observable must be computed in the soft limit, with a coefficient which is the probability of no MC emission, that is, the Sudakov form factor for the range $`(x_0,1)`$. The second term corresponds to the one-emission contribution. We stress that eqs. (208) and (B.3) could have been obtained directly from eq. (203). We can now substitute eqs. (208) and (LABEL:IMCsoftexpt) into eq. (20). We get:

$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{msub}}`$ $`=`$ $`{\displaystyle _0^1}dx[\delta (OO(x)){\displaystyle \frac{a[R(x)BQ(x)]}{x}}`$ (211)
$`+\delta (OO(0))\left(B+aV{\displaystyle \frac{aB}{x}}\right)`$
$`+aB\delta (OO(0))\left({\displaystyle \frac{Q(x)}{x}}{\displaystyle _{x_0}^1}𝑑t{\displaystyle \frac{Q(t)}{t}}\right)`$
$`+aB{\displaystyle _{x_0}^1}dt\delta (OO(t)){\displaystyle \frac{Q(t)}{t}}]+𝒪(a^2).`$
We thus have

$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{msub}}`$ $`=`$ $`{\displaystyle _0^1}𝑑x\left[\delta (OO(x)){\displaystyle \frac{aR(x)}{x}}+\delta (OO(0))\left(B+aV{\displaystyle \frac{aB}{x}}\right)\right]`$ (212)
$`+`$ $`aB{\displaystyle _0^{x_0}}𝑑x{\displaystyle \frac{Q(x)}{x}}\left[\delta (OO(0))\delta (OO(x))\right]+𝒪(a^2).`$
If we use eq. (212) to compute $`O=𝑑OO(d\sigma /dO)`$, we get the NLO result of eq. (13), up to uncontrolled $`𝒪(a^2)`$ terms, and up to terms suppressed by powers of $`x_0`$, which are numerically very small, and can be eliminated by letting $`x_00`$. In the case of QCD, the role of $`x_0`$ is played by some soft and collinear cutoffs, which cannot be set to zero. Thus, in the real case we have power-suppressed corrections to the $`𝒪(\alpha _\mathrm{s})`$ result. Notice that MC’s do not in any case control terms of this kind, and the dependence upon the cutoff is in general power-like. Since, as discussed in sect. 3.3.2, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ and the MC behave the same with respect to the cutoff, the result in eq. (212) implies that the dependence upon the cutoff of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ result beyond $`𝒪(\alpha _\mathrm{s})`$ is also power-like.

Had we applied the same procedure using naive subtraction, eq. (19), we would have obtained the following (after letting $`x_00`$):

$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{naive}}`$ $`=`$ $`{\displaystyle _0^1}𝑑x\left[\delta (OO(x)){\displaystyle \frac{aR(x)}{x}}+\delta (OO(0))\left(B+aV{\displaystyle \frac{aB}{x}}\right)\right]`$ (213)
$`+aB{\displaystyle _0^1}𝑑x{\displaystyle \frac{Q(x)}{x}}\left[\delta (OO(x))\delta (OO(0))\right]+𝒪(a^2).`$
The second integral on the r.h.s. of this equation vanishes upon integration over the whole range of $`O`$, and the total rate (if meaningful) is then identical to that obtained from eq. (212). On the other hand, if we compute $`O`$ the coefficient of $`O(x)`$ reads $`a(R(x)+BQ(x))/x`$, instead of the $`aR(x)/x`$ that we get from eq. (212). This proves in general that the naive subtraction method suffers from double counting, as already shown in eq. (30) and (42).

We now turn to the case of the QCD $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. The situation is completely analogous to that of the toy model. The NLO result can be recovered only if the $`𝒪(\alpha _\mathrm{s})`$ terms of the expansion of $`I_{\mathrm{MC}}(O,\mathrm{𝟐})`$, multiplied by the Born term, cancel the contributions of the terms $`d\overline{\mathrm{\Sigma }}|_{\mathrm{MC}}`$ in eq. (76). The QCD case is however more involved, owing to the fact that, in the case of emission from initial-state partons, the Sudakov factor has a non-trivial dependence upon the parton densities. The analogue of eq. (LABEL:IMCsoftexpt) in the QCD case reads as follows

$`I_{\mathrm{MC}}(O,\mathrm{𝟐})=(1{\displaystyle \frac{\alpha _\mathrm{s}(\mu _{\mathrm{MC}})}{2\pi }}{\displaystyle _{q_c}}{\displaystyle \frac{d\xi _+}{\xi _+}}{\displaystyle \frac{dz_+}{z_+}}\mathrm{\Theta }(z_+^2\xi _+)P_{ac}^{(0)}(z_+){\displaystyle \frac{f_c^{(H_1)}(x_1^{(\mathrm{MC})}/z_+,\mu _{\mathrm{MC}})}{f_a^{(H_1)}(x_1^{(\mathrm{MC})},\mu _{\mathrm{MC}})}}`$
$`{\displaystyle \frac{\alpha _\mathrm{s}(\mu _{\mathrm{MC}})}{2\pi }}{\displaystyle _{q_c}}{\displaystyle \frac{d\xi _{}}{\xi _{}}}{\displaystyle \frac{dz_{}}{z_{}}}\mathrm{\Theta }(z_{}^2\xi _{})P_{bc}^{(0)}(z_{}){\displaystyle \frac{f_c^{(H_2)}(x_2^{(\mathrm{MC})}/z_{},\mu _{\mathrm{MC}})}{f_b^{(H_2)}(x_2^{(\mathrm{MC})},\mu _{\mathrm{MC}})}})\delta (OO(\mathrm{𝟐}))`$
$`+\mathrm{}`$ (214)
where the dots indicate two more terms of $`𝒪(\alpha _\mathrm{s})`$, identical to those appearing in eq. (214) except for the formal substitution $`\delta (OO(\mathrm{𝟐}))\delta (OO(\mathrm{𝟑}))`$. In eq. (214) we understand that the initial-state partons that initiate the shower have flavours $`a`$ and $`b`$. The lower bound on the integration range, $`q_c`$, collectively denotes the cutoffs that prevent MC emissions from being soft and/or collinear, and plays the same role as $`x_0`$ in eq. (LABEL:IMCsoftexpt). The arguments of the parton densities have already been discussed in app. A.5. We also indicate explicitly the dependence upon the hard scale $`\mu _{\mathrm{MC}}`$, used by the MC in the shower; in principle, this scale depends on the kinematics of the emission. However, this is not important in what follows, and the notation we use simplifies the discussion. The Born term that gets multiplied by eq. (214) reads as follows:

$$\frac{d\overline{\mathrm{\Sigma }}_{ab}^{(b)}}{d\varphi _2}=\frac{(x_1^{(s)},x_2^{(s)})}{(z_1,z_2)}f_a^{(H_1)}(x_1^{(s)},\mu _{\mathrm{MC}\mathrm{@}\mathrm{NLO}})f_b^{(H_2)}(x_2^{(s)},\mu _{\mathrm{MC}\mathrm{@}\mathrm{NLO}})\frac{d\sigma _{ab}^{(b)}}{d\varphi _2}(x_1^{(s)}P_1,x_2^{(s)}P_2).$$
(215)
The $`𝒪(\alpha _\mathrm{s})`$ terms in the product of eqs. (214) and (215) must cancel the contributions of $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$ in eq. (76). It is apparent that a necessary condition for this to happen is $`x_i^{(\mathrm{MC})}x_i^{(s)}`$; this is precisely eq. (165). However, the present proof adds one piece of information which could not be inferred from eq. (165): in fact, that equation merely states the fact that the functional form of $`x_i^{(\mathrm{MC})}`$ can be obtained through the procedure of event projection; however, it does not constrain the $`x_i^{(\mathrm{MC})}`$ variables to be equal to the $`x_i^{(s)}`$ variables used in the event projection at the NLO level. This logical step is made only here. Once the condition on the Bjorken $`x`$’s is fulfilled, we achieve the desired cancellation only if the short-distance cross section resulting from eqs. (214) and (215) match those in $`d\overline{\mathrm{\Sigma }}_{ab}|_{\mathrm{MC}}`$. The former have the form reported in eqs. (LABEL:HWxsecpl) and (LABEL:HWxsecmn), whereas we have chosen the latter as in eqs. (188) and (189); therefore, we seemingly have uncancelled $`𝒪(\alpha _\mathrm{s})`$ terms. However, as argued in app. A.5, the difference between the two forms cannot induce any effects on infrared-safe observables. This is in fact what we observe in practice (see sect. 4.7). In other words, we do not have event-by-event cancellation, but we do have cancellation on a statistical basis. We remind the reader that the uncancelled $`𝒪(\alpha _\mathrm{s})`$ terms would disappear event by event if the MC had the correct angular distribution for soft gluon emission. We point out that, owing to the presence of the cutoff $`q_c`$ in eq. (214), we still have a small uncancelled $`𝒪(\alpha _\mathrm{s})`$ contribution, which is the analogue of the last term on the r.h.s. of eq. (212). As in the case of the toy model, this contribution is suppressed by powers of the cutoff. Notice that, because of the simultaneous presence of soft and collinear singularities in QCD, these powers can in general be multiplied by powers of the logarithm of the cutoff, but they remain negligible.

We also mention the fact that the flavour structure of eq. (214) is non-trivial. For the cancellation at $`𝒪(\alpha _\mathrm{s})`$ to occur, we need to ensure that the flavour assignment in the initial conditions for $`𝕊`$ events matches the flavour structure of the Born term, eq. (215). This is in fact guaranteed by the prescription described in sect. 4.5.

We remark finally that the mismatch between $`\mu _{\mathrm{MC}}`$ and $`\mu _{\mathrm{MC}\mathrm{@}\mathrm{NLO}}`$ is formally of higher order. The two scales have fairly different physical meaning, and one cannot in general assume that $`\mu _{\mathrm{MC}}=\mu _{\mathrm{MC}\mathrm{@}\mathrm{NLO}}`$. However, since the evolution of the parton densities in the MC is only leading-logarithmic, in the present paper we have made that assumption. A limited study of scale dependence has been presented in sect. 4.7.

### B.4 Resummation of large logarithms

Here we give arguments to indicate that the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ resums large logarithms in the same way as the ordinary MC. As previously, we start with the toy model, by reconsidering the results of section 3.4 for the variable $`z`$. We assume that eq. (37) is replaced by the more general form

$$\left(\frac{d\sigma }{dz}\right)_{\mathrm{MC}}=BS(z,a).$$
(216)
The function $`S`$ is largely arbitrary, but it should at least have the following properties:

* $`S`$ has a power expansion in $`a`$, starting with an $`𝒪(a)`$ term;
* the leading behaviour at small $`z`$ is $`a/z`$.

The function in eq. (37) has these properties. Equation (216) is, however, more general; we can for example assume that $`S`$ correctly takes into account next-to-leading logarithmic terms. Eqs. (38) and (39) now become

$`I_{\mathrm{MC}}(z,1)`$ $`=`$ $`S(z,a),`$ (217)
$`I_{\mathrm{MC}}(z,x_\mathrm{M}(x))`$ $`=`$ $`\delta (zx)+S(z,a)\mathrm{\Theta }(x_\mathrm{M}(x)z).`$ (218)
We thus get

$`\left({\displaystyle \frac{d\sigma }{dz}}\right)_{\mathrm{msub}}`$ $`=`$ $`S(z,a)\left[\sigma _{tot}+a{\displaystyle \frac{R(z)BQ(z)}{zS(z,a)}}a{\displaystyle _{x_\mathrm{M}^1(z)}^1}𝑑x{\displaystyle \frac{R(x)BQ(x)}{x}}\right].`$ (219)
Thus, since $`S(z,a)`$ diverges as $`1/z`$ for $`z0`$, we have

$$\left(\frac{d\sigma }{dz}\right)_{\mathrm{msub}}\sigma _{tot}S(z,a)+\mathrm{}$$
(220)
where the neglected terms are either constant for $`z0`$, or vanish in this limit. Eq. (220) generalizes eq. (45). Formally, the $`𝒪(a)`$ terms in $`\sigma _{tot}`$ produce a new tower of NLL terms when multiplied by the LL tower in $`S(z,a)`$. But since $`\sigma _{tot}`$ is factorized, this new tower of logarithms does not produce any change in the $`z`$ spectrum of the emitted photons compared to the MC. This is in fact evident if we normalize the $`z`$ distribution to the total exclusive rate, which is $`\sigma _{tot}`$ in the case of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, and $`B`$ in the case of MC: in both cases, we get $`S(z,a)`$. Furthermore, the terms neglected in eq. (45) do not lead to any logarithmic-enhanced behaviour. Thus, the leading and whatever subleading logarithmic behaviour in $`z`$ of soft photon emission is the same in MC and $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ evolution.

Notice that eq. (220) implies that the ratio of the number of soft photons to the number of showers is the same as in the ordinary MC. This fact, and the fact that the energy spectrum of the emitted photons is also the same, suggests that, for an arbitrary observable $`O`$ that gets resummed through ordinary MC methods, the resummation performed by the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ has the same logarithmic accuracy as the MC.

However, the preceding arguments are fairly closely linked to the details of the toy model, which makes it difficult to extend them to the QCD case. We therefore now give another argument, which can easily be extended to QCD; this argument works only for a given class of observables, which however includes most of the common measurable quantities. Namely, we consider the case of an exclusive observable $`O`$ such that, when $`OO_s`$, any fixed-order cross section contains large logarithms $`\mathrm{log}(OO_s)`$. Furthermore, we assume that, regardless of the number of photons emitted, at NLO and/or by the MC, when $`OO_s`$ the energies of all photons emitted are constrained to be small; in other words, any hard photon emission moves $`O`$ away from $`O_s`$. This requirement is not very restrictive; still, there can be observables which do not meet it, and have large logarithms in the cross sections, but these logarithms will be multiplied by at least an extra factor of $`a`$ relative to the case we treat.

Our assumption implies that, when we consider $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$ for $`OO_s`$, the energy $`x`$ of the photon emitted at NLO is effectively a function of $`O`$:

$$xx(O)=\kappa (OO_s)+𝒪\left((OO_s)^2\right),$$
(221)
where we used the fact that $`x(O_s)=0`$. We can use eq. (221) to compute the NLO integrated cross section. We have

$`\overline{\sigma }_{\mathrm{NLO}}(O){\displaystyle _{O_s}^O}𝑑\overline{O}\left({\displaystyle \frac{d\sigma }{d\overline{O}}}\right)_{\mathrm{NLO}}`$
$`={\displaystyle _0^1}𝑑x\left[{\displaystyle _{O_s}^O}𝑑\overline{O}\delta (\overline{O}O(x)){\displaystyle \frac{aR(x)}{x}}+{\displaystyle _{O_s}^O}𝑑\overline{O}\delta (\overline{O}O(0))\left(B+aV{\displaystyle \frac{aB}{x}}\right)\right],`$
(222)
where we have used eq. (13). With eq. (221) we obtain

$$\overline{\sigma }_{\mathrm{NLO}}(O)=B+aV+_0^1𝑑x\left[\frac{aR(x)}{x}\mathrm{\Theta }(\kappa (OO_s)x)\frac{aB}{x}\right].$$
(223)
Consistently with eq. (4), we assume that

$$R(x)=B+\underset{n=1}{\overset{\mathrm{}}{}}R_nx^n,$$
(224)
and thus eq. (223) gives

$$\overline{\sigma }_{\mathrm{NLO}}(O)=B+aB\mathrm{log}\left(\kappa (OO_s)\right)+aV+a\underset{n=1}{\overset{\mathrm{}}{}}\frac{R_n}{n}\left(\kappa (OO_s)\right)^n.$$
(225)
The first two terms on the right-hand side are the first two terms in the LL tower that is resummed by MC means. The third and fourth terms are NLL, the latter being also power suppressed. We point out that since the only large logarithm in the cross section has to be $`\mathrm{log}(OO_s)`$, the coefficient $`\kappa `$ must be of order unity.

We now consider the case of the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. We denote the normalized pure MC result as follows:

$$\overline{\sigma }_{\mathrm{MC}}(O)_{O_s}^O𝑑\overline{O}I_{\mathrm{MC}}(\overline{O},1).$$
(226)
We also notice that, owing to the fact that $`x_\mathrm{M}(0)=1`$, and assuming that $`x_\mathrm{M}(x)`$ has a power expansion in $`x`$ around $`x=0`$, we have (for $`OO_s`$)

$$_{O_s}^O𝑑\overline{O}I_{\mathrm{MC}}(\overline{O},x_\mathrm{M}(x))=_{O_s}^O𝑑\overline{O}\left(I_{\mathrm{MC}}(\overline{O},1)+\mathrm{NLL}\right)\mathrm{\Theta }(\kappa (OO_s)x).$$
(227)
We can now repeat the same computation that led us to eq. (225), for the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, eq. (20). In the case of $`𝕊`$ events, analogously to what happens for the contribution of the 0-body events in eq. (222), the restricted range in $`O`$ does not constrain the integration range in $`x`$. Thus

$`\overline{\sigma }_𝕊(O)`$ $``$ $`{\displaystyle _{O_s}^O}𝑑\overline{O}{\displaystyle \frac{d\sigma _𝕊}{d\overline{O}}}`$ (228)
$`=`$ $`\overline{\sigma }_{\mathrm{MC}}(O)\left[B+aV+aB{\displaystyle _0^1}𝑑x{\displaystyle \frac{Q(x)1}{x}}\right];`$
we therefore recover the pure MC result, up to $`𝒪(a)`$ terms in the normalization. On the other hand, the restricted range in $`O`$ does constrain the integration range in $`x`$ in the case of $``$ events, as shown in eq. (227). We obtain

$`\overline{\sigma }_{}(O)`$ $``$ $`{\displaystyle _{O_s}^O}𝑑\overline{O}{\displaystyle \frac{d\sigma _{}}{d\overline{O}}}`$ (229)
$`=`$ $`a\left(\overline{\sigma }_{\mathrm{MC}}(O)+\mathrm{NLL}\right){\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{R_nBQ_n}{n}}\left(\kappa (OO_s)\right)^n,`$
where $`Q_n`$ are the coefficients of the expansion of $`Q(x)`$ around $`x=0`$. Therefore, eq. (229) does not contribute to the large-logarithm resummation for $`OO_s`$. Taking into account eq. (228), this allows us to conclude that the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ resums large logarithms in the same way as the ordinary MC does.

It is clear that the argument given above works in the same way in the QCD case. One has to pay attention to the fact that, in QCD, soft and collinear singularities exist. In general, the restricted range in $`O`$ constrains the emissions to be either soft or collinear, not necessarily soft and collinear. This implies that in QCD each term of the power series that appears in eq. (229) can be multiplied by a power of $`\mathrm{log}(OO_s)`$. However, this does not change our conclusion, that $``$ events do not contribute large logarithms to the physical cross section.

## Appendix C $`\mathrm{𝐌𝐂}\mathbf{@}\mathrm{𝐍𝐋𝐎}`$ based upon slicing method

In this appendix, we comment on the approach of refs. , where a formalism is presented to construct an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ using the slicing method at the NLO level. This formalism is called $`\mathrm{\Phi }`$-space Veto in ref. . It has the main objective and virtue of ensuring that only positive event weights are generated. However, we argue that this is at the expense of being able to recover the NLO results upon expansion in $`\alpha _\mathrm{s}`$ (this is because the resulting $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ does not have an expansion in $`\alpha _\mathrm{s}`$), and a smooth matching between different kinematic regions.

In the course of the discussion, we use the toy model, which is sufficient to highlight the main features of the method. We also propose a formula for an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ based upon the slicing method, which does not have any of the problems that affect the $`\mathrm{\Phi }`$-space Veto method, and fulfills the definition of $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ given in sect. 2.1 (but as a consequence is not guaranteed to generate positive event weights).

The $`\mathrm{\Phi }`$-space Veto relies on the result for the NLO cross section obtained with the slicing method, eq. (9). The first step consists in fixing the parameter $`\delta `$, by imposing the condition that eq. (9) does not contain any contribution proportional to $`O(0)`$. In other words, one seeks $`\delta _0`$ such that999This idea was first introduced, in a slightly different form, in ref. .

$$B+a\left(B\mathrm{log}\delta _0+V\right)=0.$$
(230)
Thus

$$\mathrm{log}\delta _0=\frac{B+aV}{aB}\delta _0=\mathrm{exp}\left(\frac{B+aV}{aB}\right).$$
(231)
With such a choice, eq. (9) becomes

$$O_{\mathrm{slice}}=a_{\delta _0}^1𝑑x\frac{O(x)R(x)}{x}+𝒪(\delta _0).$$
(232)
Given the fact that the parameter $`\delta _0`$ is not free any longer, being related to the Born and virtual contributions as shown in eq. (231), the uncontrolled $`𝒪(\delta _0)`$ terms in eq. (232) may be non-negligible. An improvement of eq. (232) is proposed in ref. , which uses a subtraction method in a restricted part of the phase space (this procedure is called “hybrid subtraction”): however, this improvement is of no importance in what follows, and will be neglected here (also, it is not used in current implementations of refs. ). As anticipated, the reason for fixing the slicing parameter $`\delta `$ can be seen in eq. (232); this equation does not contain any negative contributions, and this fact is used in refs.  to define an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ with no negative-weight events. However, the resulting $`\delta _0`$ is not analytic around $`a=0`$, which seems to us to cause problems, as will be discussed below.

The next step in the $`\mathrm{\Phi }`$-space Veto method is the introduction of a free parameter $`\delta _{\mathrm{PS}}`$, which fulfills the condition $`\delta _0<\delta _{\mathrm{PS}}<1`$. In ref. , it is suggested that one set $`\delta _{\mathrm{PS}}`$ equal to the boundary of the dead zone, in the case that HERWIG is used to construct the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. In general, it seems difficult to prove that this condition is consistent with the requirement $`\delta _0<\delta _{\mathrm{PS}}`$, especially since the boundary of the dead zone has a different (process dependent) shape from the region removed by slicing, and may depend on different variables.

In the language of the present paper, the introduction of $`\delta _{\mathrm{PS}}`$ implies the definition of a modified slicing formula:

$$O_{\mathrm{mslice}}=a_0^1𝑑x\left[O(x)\frac{R(x)}{x}\mathrm{\Theta }(x\delta _{\mathrm{PS}})+O(0)\frac{R(x)}{x}\mathrm{\Theta }(x\delta _0)\mathrm{\Theta }(\delta _{\mathrm{PS}}x)\right].$$
(233)
From this equation, we can construct a toy $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ as explained in sect. 3.3, paying attention to the fact that, in the $`\mathrm{\Phi }`$-space Veto method, the maximum energy available to the photon in the first branching is required to be $`\delta _{\mathrm{PS}}`$ (instead of 1) for class $`𝕊`$ events. Therefore

$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{mslice}}`$ $`=`$ $`a{\displaystyle _0^1}dx[I_{\mathrm{MC}}(O,x_\mathrm{M}(x)){\displaystyle \frac{R(x)}{x}}\mathrm{\Theta }(x\delta _{\mathrm{PS}})`$ (234)
$`+I_{\mathrm{MC}}(O,\delta _{\mathrm{PS}}){\displaystyle \frac{R(x)}{x}}\mathrm{\Theta }(x\delta _0)\mathrm{\Theta }(\delta _{\mathrm{PS}}x)].`$
Eq. (234) can now be used to compute the predictions of the $`\mathrm{\Phi }`$-space Veto method, in the same way that eq. (20) was used for our $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$. We consider here the exclusive observable $`y`$, introduced in sect. 3.4. We can substitute eq. (26) directly into eq. (234); on the other hand, we need to use

$$I_{\mathrm{MC}}(y,\delta _{\mathrm{PS}})=a\frac{Q(y)}{y}\mathrm{\Delta }(y,\delta _{\mathrm{PS}})\mathrm{\Theta }(\delta _{\mathrm{PS}}y)$$
(235)
instead of eq. (25). We find the following:

$$\left(\frac{d\sigma _𝕊}{dy}\right)_{\mathrm{mslice}}=a^2\frac{Q(y)}{y}\mathrm{\Delta }(y,\delta _{\mathrm{PS}})\mathrm{\Theta }(\delta _{\mathrm{PS}}y)_{\delta _0}^{\delta _{\mathrm{PS}}}𝑑x\frac{R(x)}{x}.$$
(236)
We also find:

$``$ for $`y<x_e:`$

$$\left(\frac{d\sigma _{}}{dy}\right)_{\mathrm{mslice}}=\left[a\frac{R(y)}{y}\mathrm{\Delta }(y,x_\mathrm{M}(y))+a^2\frac{Q(y)}{y}_{\delta _{\mathrm{PS}}}^y\frac{dx}{x}R(x)\mathrm{\Delta }(y,x_\mathrm{M}(x))\right]\mathrm{\Theta }(y\delta _{\mathrm{PS}});$$
(237)
$``$ for $`y>x_e:`$

$$\left(\frac{d\sigma _{}}{dy}\right)_{\mathrm{mslice}}=a\frac{R(y)}{y}\mathrm{\Theta }(y\delta _{\mathrm{PS}})+a^2\frac{Q(y)}{y}\mathrm{\Theta }(x_\mathrm{M}^1(y)\delta _{\mathrm{PS}})_{\delta _{\mathrm{PS}}}^{x_\mathrm{M}^1(y)}\frac{dx}{x}R(x)\mathrm{\Delta }(y,x_\mathrm{M}(x)).$$
(238)
At $`y0`$, we expect the result of the ordinary MC, eq. (24), to be reproduced, up to terms of higher logarithmic order, and up to multiplicative factors of higher order in $`a`$, i.e., we expect an equation similar to (35) to hold. Eq. (236) appears to fail to achieve this result, since it starts at $`𝒪(a^2)`$. However, $`\delta _0`$ also depends upon $`a`$; thus, in order to reach firm conclusions, the integral in eq. (236) needs to be evaluated explicitly. To do that, we assume again that $`R(x)`$ has a power series expansion in $`x`$, eq. (224). Eq. (236) then becomes

$$\left(\frac{d\sigma _𝕊}{dy}\right)_{\mathrm{mslice}}=a^2\frac{Q(y)}{y}\mathrm{\Delta }(y,\delta _{\mathrm{PS}})\mathrm{\Theta }(\delta _{\mathrm{PS}}y)\left[B\mathrm{log}\frac{\delta _{\mathrm{PS}}}{\delta _0}+\underset{n=1}{\overset{\mathrm{}}{}}\frac{R_n}{n}\left(\delta _{\mathrm{PS}}^n\delta _0^n\right)\right].$$
(239)
Using eq. (231) we now find

$`\left({\displaystyle \frac{d\sigma _𝕊}{dy}}\right)_{\mathrm{mslice}}`$ $`=`$ $`a{\displaystyle \frac{Q(y)}{y}}\mathrm{\Delta }(y,\delta _{\mathrm{PS}})\mathrm{\Theta }(\delta _{\mathrm{PS}}y)[B+aV+aB\mathrm{log}\delta _{\mathrm{PS}}`$ (240)
$`+a{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{R_n}{n}}(\delta _{\mathrm{PS}}^n\mathrm{exp}(n{\displaystyle \frac{B+aV}{aB}}))].`$
This equation has the same functional form in $`y`$ as eq. (24), up to uncontrolled NLL terms (due to the fact that the second argument of the Sudakov form factor is $`\delta _{\mathrm{PS}}`$ instead of 1), and formally starts at $`𝒪(a)`$. However, it has the rather unpleasant feature that it cannot be expanded as a power series in $`a`$, since it is not analytic at $`a=0`$. We therefore have the paradoxical situation that a quantity which is based on perturbation theory cannot be expanded in powers of the coupling constant. This is in fact due to the form of $`\delta _0`$, eq. (231). Still, if we assume that $`a^2\delta _0^n`$ is of $`𝒪(a^2)`$, then eq. (240) has the correct Sudakov behaviour for $`y0`$ (this is also formally consistent with the necessity of neglecting powers of $`\delta _0`$ in the slicing procedure; however, we believe that this is not justified here, since $`\delta _0`$ is not a free parameter any longer). However, with the same assumption we conclude that the $`\mathrm{\Phi }`$-space Veto method has double counting problems, since for $`y<\delta _{\mathrm{PS}}`$ the result of eq. (240) does not coincide with the pure NLO result. On the other hand, for $`y>\delta _{\mathrm{PS}}`$ the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ and NLO results agree up to $`𝒪(a^2)`$ terms (see eqs. (237) and (238)). Therefore, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ prediction has a discontinuity at $`y=\delta _{\mathrm{PS}}`$.

We now define an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ based upon the slicing method that does not have the problems outlined above. We can easily accomplish this task by noting that the slicing formulae can be obtained, at the NLO level, by making suitable approximations in the subtraction formulae. We start from eq. (9), which we re-write as follows

$$O_{\mathrm{slice}}=_\delta ^1𝑑x\left[O(x)\frac{aR(x)}{x}+O(0)\left(B+aV+aB\mathrm{log}\delta \right)\right].$$
(241)
From now on, we do not write explicitly the terms of order $`𝒪(\delta )`$ that are always neglected in slicing formulae. Notice that eqs. (9) and (241) in fact coincide up to $`𝒪(\delta )`$ terms; the two equations would be exactly identical had we multiplied $`O(0)`$ by $`1/(1\delta )`$. Using eq. (241) we now define an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$, following the procedure described in sect. 3.3. We obtain

$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{mslice}}`$ $`=`$ $`{\displaystyle _\delta ^1}dx[I_{\mathrm{MC}}(O,x_\mathrm{M}(x)){\displaystyle \frac{a[R(x)BQ(x)]}{x}}`$ (242)
$`+I_{\mathrm{MC}}(O,1)(B+aV+aB\mathrm{log}\delta +aB{\displaystyle \frac{Q(x)}{x}})].`$
We stress that this formula has the same features as a “standard” slicing formula, such as eq. (241): the parameter $`\delta `$ is free (and small), the integration over $`x`$ does not extend down to zero, and the counterterm of the subtraction method is replaced by a term proportional to $`\mathrm{log}\delta `$. There is, however, one important difference. In eq. (241), the weights attached to $`O(x)`$ and $`O(0)`$ are divergent when $`\delta 0`$ (and can now be negative). This does not happen in eq. (242) for the weights attached to $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$ and $`I_{\mathrm{MC}}(O,1)`$. In the latter term, however, this property holds only if the integration over the phase space is performed before the limit $`\delta 0`$ is carried out. This is inconvenient, because the unweighting procedure is performed in practice before the integration. In the spirit of the slicing method, we should integrate the term $`aBQ(x)/x`$ in the weight attached to $`I_{\mathrm{MC}}(O,1)`$ analytically, which may or may not be possible, depending upon the form of $`Q(x)`$. In order to simplify the discussion, we assume

$$Q(x)=\mathrm{\Theta }(x_{\mathrm{dead}}x).$$
(243)
Using this form in eq. (242) we obtain

$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{mslice}}`$ $`=`$ $`{\displaystyle _\delta ^1}𝑑xI_{\mathrm{MC}}(O,x_\mathrm{M}(x)){\displaystyle \frac{a[R(x)B\mathrm{\Theta }(x_{\mathrm{dead}}x)]}{x}}`$ (244)
$`+I_{\mathrm{MC}}(O,1)\left(B+aV+aB\mathrm{log}x_{\mathrm{dead}}\right).`$
In general, the weights attached to $`I_{\mathrm{MC}}(O,x_\mathrm{M}(x))`$ and $`I_{\mathrm{MC}}(O,1)`$ in this equation are not positive definite. However, it is clear that unweighting (in the sense of generating events with weights of equal magnitude) can easily be achieved, since all the cancellations between large numbers occur under the integral sign. In this sense, the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ based upon the slicing method is similar to the $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ based upon the subtraction method, eq. (20). However, while in eq. (20) the cancellation is achieved by numerical methods, the generalization of eq. (244) can be obtained from eq. (242) only through an analytical integration. This is similar to what happens at the NLO level; in the subtraction method, the residue of the soft/collinear singularity is integrated numerically, whereas in the slicing method an analytical integration is performed.

In the general case in which $`Q(x)`$ does not give an easy analytical integral, we can use a mixture of analytical and numerical methods. Namely, we use the identity

$$Q(x)=\mathrm{\Theta }(x_{th}x)+Q(x)\mathrm{\Theta }(x_{th}x)$$
(245)
in the last term of eq. (242), where $`x_{th}`$ is an arbitrary constant. We readily obtain

$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{mslice}}`$ $`=`$ $`{\displaystyle _\delta ^1}dx[I_{\mathrm{MC}}(O,x_\mathrm{M}(x)){\displaystyle \frac{a[R(x)BQ(x)]}{x}}`$ (246)
$`+I_{\mathrm{MC}}(O,1)(B+aV+aB\mathrm{log}x_{th}+aB{\displaystyle \frac{Q(x)\mathrm{\Theta }(x_{th}x)}{x}})].`$
We can now prove that eq. (242) has no double counting, and resums soft logarithms in the same way as standard MC does. We follow the procedures introduced in app. B.3 and B.4. Using eqs. (208) and (LABEL:IMCsoftexpt), we get

$`\left({\displaystyle \frac{d\sigma }{dO}}\right)_{\mathrm{mslice}}`$ $`=`$ $`{\displaystyle _\delta ^1}𝑑x\left[\delta (OO(x)){\displaystyle \frac{aR(x)}{x}}+\delta (OO(0))\left(B+aV+aB\mathrm{log}\delta \right)\right]`$ (247)
$`+`$ $`aB{\displaystyle _0^1}𝑑x{\displaystyle \frac{Q(x)}{x}}[\delta (OO(0))\delta (OO(x))]\left[\mathrm{\Theta }(x\delta )\mathrm{\Theta }(xx_0)\right]`$
$`+`$ $`𝒪(a^2),`$
which has to be compared to eq. (212). As in that case, we recover the NLO result by considering $`𝑑OO(d\sigma /dO)`$, up to $`𝒪(a^2)`$ and power-suppressed terms. Notice that we also have terms of $`𝒪(\delta )`$ in eq. (247), consistently with the approximation upon which the slicing method is based. The distribution of the photon energies can be obtained as explained in app. B.4. We get

$`\left({\displaystyle \frac{d\sigma }{dz}}\right)_{\mathrm{mslice}}`$ $`=`$ $`S(z,a)\left[\sigma _{tot}+a{\displaystyle \frac{R(z)BQ(z)}{zS(z,a)}}+aB{\displaystyle _{\mathrm{max}(\delta ,x_\mathrm{M}^1(z))}^1}𝑑x{\displaystyle \frac{Q(x)}{x}}\right],`$ (248)
which is the analogue of eq. (219). It is apparent that the same conclusions as in app. B.4 apply here for $`z0`$.

Finally, notice that, by using eq. (206), we obtain, in the case of an exclusive variable

$$𝑑O\left(\frac{d\sigma }{dO}\right)_{\mathrm{mslice}}=_\delta ^1𝑑x\left[B+aV+aB\mathrm{log}\delta +\frac{aR(x)}{x}\right]\sigma _{tot},$$
(249)
which is the result that we also find in app. B.2. We conclude that eq. (242) defines an $`\mathrm{MC}\mathrm{@}\mathrm{NLO}`$ that has the same formal properties as that defined in eq. (20).