Practical improvements and merging of POWHEG simulations for vector boson production

In this article we generalise POWHEG next-to-leading order parton shower (NLOPS) simulations of vector boson production and vector boson production in association with a single jet, to give matrix element corrected MENLOPS simulations. In so doing we extend and provide, for the first time, an exact and faithful implementation of the MENLOPS formalism in hadronic collisions. We also consider merging the resulting event samples according to a phase space partition defined in terms of an effective jet clustering scale. The merging scale is restricted such that the component generated by the associated production simulation does not impact on the NLO accuracy of inclusive vector boson production observables. The dependence of the predictions on the unphysical merging scale is demonstrated. Comparisons with Tevatron and LHC data are presented.


Introduction
In order to discover and distinguish the possible production of new particles at the Tevatron and LHC, a thorough understanding of the backgrounds is fundamental. In this respect vector boson production is an all too common obstacle in searches for physics beyond the standard model. Furthermore, given the high energies afforded by these colliders, in particular the LHC, and the high-p T regions of phase space to which new physics searches are sensitive, the emission of associated QCD radiation, manifest as hard and soft jets, is ubiquitous.
The role played by W and Z boson production at hadron colliders is of course not simply negative. Production of these particles in such large quantities facilitates precision measurements of e.g. the W boson mass, compensating their obscuring of direct new physics signals by probing indirect ones. W and Z production are also an important tool in furthering our knowledge of other fundamentals such as parton distributions, their uncertainties and also detector / jet energy calibration. From a purely practical point of view there is therefore a clear need for accurate and detailed simulations of these processes.
Of course vector boson production is not the only standard model process for which precision event generators are well motivated and so it comes as little surprise, that in recent years the nearing of LHC data-taking has fueled a resurgence in the research and development of Monte Carlo simulations. Ground-breaking work carried out by various groups in 2000-2004 demolished significant long standing problems, radically improving the fully exclusive, hadron-level, description of parton shower simulations to consistently include multi-leg tree-order matrix elements [1][2][3][4][5] and, separately, exact next-to-leading order corrections [6][7][8]; these are referred to as Meps and Nlops simulations respectively.
In the following years many simulation programs have been constructed, implementing these new methods and techniques, for a wide variety of collider physics processes, and their worth proven in various experimental applications. Progress in this direction has rapidly evolved to the point where there presently exist three public computer packages automating, or partially automating, the construction of Nlops programs, such that little or no expert knowledge is required on the part of the user to produce a simulation for a given process [9][10][11][12]. Moreover, these packages and the methods underlying them can now be considered to have withstood the test of complexity, with the recent arrival of public Nlops programs for processes involving rich colour structures and infrared divergences at leading order [13][14][15], as well as high multiplicity processes [11,12,16,17].
Having reached this degree of maturity, with systematic differences of Nlops schemes now well understood and code production reaching industrial scales, it becomes relevant to discuss whether further improvements are possible and how the output of the new automated Nlops packages may be best used to best effect. In Refs. [18][19][20][21][22][23] a number of novel and sophisticated methods for matching parton shower simulations to multiple higher order one-loop and tree level matrix elements have already been put forward. Applications of these techniques to hadron collider processes are eagerly anticipated.
More recently, two somewhat less ambitious Menlops methods have been introduced, in Refs. [24] and [25], for simple processes, combining a single Nlops event generator with a Meps simulation, enhancing the description of multi-jet events in the former. What these methods lack in accuracy with respect to their forerunners appears to be offset by their relative simplicity, making use of, and requiring little modification to, existing Nlops and Meps programs, with a number of applications to hadron collider processes having been performed in the referenced articles.
We take the view that there is still more to be gained by a continued bottom-up approach to improving on the established crop of Nlops and Meps simulations, as found in the case of Menlops simulations, building on existing programs rather than starting from scratch.
Our first aim in this work has been to see if we can improve on the Menlops implementations in Refs. [24] and [25], to realize the method exactly, simply and without approximation. In other words, we seek to implement the Menlops method for vector boson production without the introduction of an unphysical merging scale, in a way which manifestly respects the unitarity of the Powheg method and thus NLO accuracy.
Our second objective has been to extend the Menlops method to the case of processes involving a final-state jet at the leading order. Since the leading order matrix elements in such cases typically contain infrared divergences, these simulations require unphysical transverse momentum cut-offs to make event generation possible. This has two important practical consequences. Firstly, the description of the region in which the vector boson has low transverse momentum, where the great bulk of vector bosons are produced, is also unphysical, being highly sensitive to the cut and lacking all-orders resummation of large logarithms in that region. Secondly, when studying jet-associated production as a signal process, the users of these programs must always check that the predictions from the program do not depend on the cut. Typically this mandates using a generation cut far below those defining all observables, resulting in a very substantial fraction of events that never make it into any part of the analysis. Thus, here, by extending the Menlops method, we mean effectively including the relevant lower order, sub-leading matrix element, and the all-orders resummation required for a physical description of the low p T region. Note that the methodology here borrows much from the Sudakov reweighting / vetoing procedures in Meps merging schemes.
Lastly, we have considered merging event samples from the two Menlops simulations, to give one NLO accurate in the description of both fully inclusive and inclusive vector boson plus jet production observables. To this end we apply similar reasoning to Ref. [24], populating the region in which the vector boson has a low transverse momentum with events generated by the single vector boson Menlops simulation and elsewhere with events from the vector boson plus jet Menlops program. While requiring that both classes of observables be described with NLO accuracy puts constraints on the values that the unphysical merging scale can take, in all other respects the merging scale dependence is never worse than in the CKKW(-L) / MLM approaches. However, discounting the fact that renormalization and factorization scales are perhaps not always chosen optimally in fixed order calculations, concerns about the scale dependence may be more a matter of theoretical correctness, rather than practical importance. To minimize the dependence on the merging scale beyond that here, in particular, to retain NLO accuracy for both classes of observable independently of the merging scale, is to our understanding, tantamount to achieving genuine NNLO-parton shower matching; a goal far in advance of this simple study, beyond our humbler aims of systematically getting the most from existing and forthcoming simulations by minimal and modular interventions.
In section 2 we describe the extension of the Powheg Nlops simulations to Menlops ones, as well as the business of merging the resulting event samples. Here, when possible, we shall aim to be conceptual, proving that all relevant technical details have been sufficiently understood and controlled by reference to Monte Carlo validation plots. We kindly ask the reader to bear in mind that the Powheg formalism has been derived in considerable detail [7][8][9]26], so a fully self-contained presentation here is not feasible. In section 3 we demonstrate the improvements to be gained through comparisons to relevant Tevatron and LHC measurements. Our findings and conclusions are summarized in section 4.

Method
In this section we elaborate on the steps involved in formulating event samples NLO accurate in their description of fully inclusive and single-jet inclusive vector boson production observables. As we have described in the introduction, this proceeds by merging the output of two Menlops simulations. Thus, our presentation is organized here as follows. First, in section 2.1, we introduce some basic ideas and notation for Powheg Nlops matching, we then go on to describe how these simulations of W and Z production may be enhanced to produce genuine Menlops simulations in section 2.2, this is followed by the case of the associated production channel in section 2.3. Finally, the combination of the Menlops samples is described in section 2.4.
Since we will already show in this section results from the programs we have used and created, for illustrative and validation purposes, we first declare some of the main technical parameters used to run them. All Powheg and Menlops events in this paper have been generated using the CTEQ6M [27] parton distribution functions, interfaced via Lhapdf [28], with the corresponding value of Λ QCD . In the case of vector boson production, the renormalization and factorization scales assumed in evaluating theB functions, distributing the underlying Born kinematics (Sect. 2.1), are given by the vector boson invariant mass, while for associated production the transverse momentum of the initial hh → V j state is used. The Powheg Box defaults for the W and Z boson masses and widths are: m W = 80.398 GeV, Γ W = 2.141 GeV, m Z = 91.188 GeV, Γ Z = 2.486 GeV. The value of the QED coupling used is given by α (m Z ) = 1/128.93. Throughout the article we have used the Rivet analysis framework [29], including the FastJet package [30], to study our event samples. The Monte Carlo validation analyses in this section then include a set of minimal cuts on the pseudorapidities and transverse momenta of the charged leptons emanating from the boson decays, |η| < 3.5, p T > 25 GeV and, in the case of W production, a missing transverse energy cut, / E T ≥ 25 GeV. Lastly, in order to develop the bare Powheg and Menlops events to include effects of parton showering, hadronization and multiple interactions we have used the Pythia 8.150 [31][32][33] program with its default tune and PDF set.

Preliminaries
At the heart of the Powheg Box simulation of vector boson production, as with all Powheg simulations, is the so-called hardest emission cross section [7] i.e. the differential cross section governing the distribution of the hardest (highest p T ) parton branching in each event: In Eq. 2.1 B(Φ V ) stands for the leading order cross section, dependent on the underlying Born kinematics Φ V . Analogously, R(Φ V j ) is the real emission cross section depending on the real kinematics, Φ V j . The real kinematics are defined by a mapping which takes as arguments the leading order kinematics, Φ V , and the radiative phase space variables, Φ j 1 , parametrizing the hardest parton branching with respect to the configuration Φ V , the R(Φ V j ) term can be understood to have absorbed the Jacobian associated to the factorisation dΦ V j = dΦ V dΦ j 1 . For the Powheg Box programs used in this study such mappings are explicitly given in Ref. [8]. 1B (Φ V ) is the next-to-leading order distribution of the same (infrared-safe) underlying Born kinematics where the virtual contribution V (Φ V ) here implicitly includes soft and collinear divergences, canceling those present in the real term; to avoid digressing we simply state that in general some regularization / subtraction scheme is adopted rendering the square bracket term in Eq. 2.3 finite and regular. Lastly, the Powheg Sudakov form factor is defined as where k T (Φ V j ) tends to the transverse momentum of the radiated parton in the collinear and soft limits. As in conventional parton shower simulations k T always has an implicit lower cut-off, p min T , of order Λ QCD . Note that from the point of view of the next-to-leading order corrections, the details concerning the internal parametrization of the Born kinematics are completely irrelevant, the key point there being the mappings defining the relationship between them and the real kinematics. As an heuristic aid, neglecting, for a moment, the leptonic decay products, Φ V , could be taken to comprise of the vector boson mass and rapidity, while the radiation phase space in the Powheg Box is mapped by the polar and azimuthal angles of the emitted parton with respect to the singular direction (here the beam axis) together with its energy.
We now quickly remind the reader of the two main properties of the hardest emission cross section, Eq. 2.1, through which it gives rise to NLO predictions. Firstly, as the transverse momentum of the hardest branching increases beyond the Sudakov peak region, the form factor, ∆, tends to one and, neglecting higher order terms of relative order α 2 S , the kinematics are distributed according to the real emission cross section R(Φ V j ): More subtly and, from the point of view of this study, more significantly, the second key feature of Eq. 2.1 through which NLO accuracy attained is its unitarity with respect to the Born kinematics, specifically, in the present context, that the term in square brackets integrates to one for any Φ V . This can be readily deduced by making use of the identity where the full kinematic dependence of ∆ has been suppressed for ease of notation.
1 As we will discuss in the next paragraphs, more generally the real emission cross section needs to be separated according to the number of counterterms required to make it finite on numerical integration, each one having its own unique phase space parametrization.
For the case of jet-associated vector boson production the basic form of the hardest emission cross section is unchanged. The principal difference is in the treatment of the real phase space and real cross section. In general the real cross section is separated into pieces, R = α R α , one for each infrared counterterm required to render the cross section regular and numerically integrable, moreover, for each contribution the real phase space is parametrized differently, such that the counterterms lend themselves to being integrated analytically over the radiative phase space i.e.
Thus, the Powheg hardest emission cross section in the case of vector boson plus jet production has the form and Note that in order to ease the proliferation of indices, the functions B, V , R, ∆ andB, for the direct and jet-associated vector boson production processes, are implicitly distinguished by their phase space arguments. Also, although such details are somewhat tangential to our course, the α-dependent Jacobians associated with the phase space factorisation dΦ α V jj → dΦ V j dΦ j 2 can be understood to have been absorbed in the R α terms, unless otherwise stated.
What is important to realize is that, fundamentally, the arguments concerning how NLO accuracy manifests in this more general case are unchanged with respect to that of vector boson production. In the high p T regime the Sudakov form factor tends to one and combines with the theB prefactor to adjust the real cross section by genuinely NNLO terms only. Likewise, as with Eq. 2.1, by virtue of the fact that the exponent of the Sudakov form factor, Eq. 2.10, is precisely equal to the coefficient multiplying it in Eq. 2.8, integration over the radiative phase space reduces the square bracket term to one, leaving a next-to-leading order distribution of the leading order vector boson plus jet kinematics. Of course, the fact that the integral over the radiative phase space returns justB, does not prove that NLO accuracy is obtained for general inclusive and semi-inclusive observables but hopefully it is sufficiently conducive that we can proceed here and refer the interested reader to section. 4.3 of Ref. [8], where a rigorous proof based on this property can be found. For the purposes of our discussion it is sufficient to appreciate that this unitarity is a necessary condition for the Powheg method to yield NLO predictions, deviations from unitarity amount to deviations from NLO accuracy in equal measure. With preliminary concepts and notation in hand we now go on to discuss how we augment the Powheg simulations of vector boson production and vector boson plus jet production.

Menlops improved simulation of vector boson production
In Ref. [24] it was described how one could instill NLO accuracy in Meps-type simulations for simple processes. In that case the key consideration was that the effective hardest emission cross section of Meps simulations are not in general unitary. The square bracket terms in the formulae analogous to Eq. 2.1 do not contain precisely the real emission cross section R(Φ V j ) but something equivalent to it up to relative corrections O (α S ); the modifications being due to several factors, the presence of higher order tree-level matrix elements in the simulation being one of the more obvious ones. Furthermore, the Sudakov form factor multiplying the approximate real cross section aims to be that of the parton shower simulation rather than one based on exponentiating R(Φ V j ) or its Meps approximate version. In the Meps, therefore, the integral over the radiative phase space then rather gives a function N (Φ V ) ≈ 1 + O(α S ). Based on these assertions it was noted in Ref. [24] that Meps event samples could be made to respect unitarity and reproduce NLO predictions in the same way as Powheg by reweighting their events according toB(Φ V )/N (Φ V ).
Here we take a simpler different approach, albeit one adhering to the same principles, viewing the Powheg program for the jet-associated vector boson production channel as a 'perfect' (unitary) Meps apparatus, merging vector boson plus jet and vector boson plus dijet final states.
To explain how we formulate our Menlops simulation in this section first note that, in practice, the hardest emission cross section in a Powheg simulation is implemented by initially generating an unweighted set of leading order kinematics according to the nextto-leading orderB function, then, with this configuration in hand, the radiative variables are generated with respect to it, according to the distribution in square brackets in Eq. 2.1 / 2.8. This latter radiation generation mechanism is implemented in essentially the same way as the long established matrix element corrections technique used in parton shower simulations [34][35][36][37], with the relevant distributions, ∼ ∆ R/B, being sampled according to the same veto algorithm. Thus, here we proceed by first generating a parton level configuration according to the hardest emission cross section in the vector boson production Powheg code, Eq. 2.1, with no alterations, saving the events in the Les Houches format. These events are then passed as input to a slightly modified version of the vector boson plus jet simulation, which overrides the generation of leading order configurations according tō B(Φ V j ), taking the kinematics and flavour structure instead from the Les Houches event file input. By propagating the bare Powheg events output from the vector boson production program to the radiation generating apparatus of the associated production simulation we obtain doubly radiative events distributed according to (2.11) this follows straightforwardly from the modular form of the hardest emission cross sections and Monte Carlo implementation.
In this way we exploit the unitarity of the radiation machinery in the vector boson plus jet program to best effect. As described in section. 2.1, the final term in braces is an exact differential with the integral over the radiative phase space Φ j 2 being identically equal to one for all Φ V j , thus the distributions of Born and real kinematics of the vector boson production Nlops should be completely untouched in this way. At a technical level, that this is true for the real kinematics is trivial, since those events are fed directly to the radiation generation routines of the vector boson plus jet simulations, on the other hand a much more subtle and, in fact, obligatory, verification of this procedure is to check that the predictions for fully inclusive quantities are completely unchanged by the generation of the secondary radiation in this way.
Before proceeding to the validation let us also discuss the resummation properties of Eq. 2.11. Since the Powheg Box programs implement the hardest emission cross sections using the FKS subtraction formalism, the real cross section R is separated into pieces singular in only one collinear direction. Moreover, the damping factors giving rise to the separation are such that each R α is relatively suppressed by a factor 1/p 2 T as the direction of the emitted parton deviates from the collinear one associated to α, regardless of the positioning of other partons in the event. In fact, as other collinear singular regions are approached R α is defined to vanish [8].
At least from the point of view of FKS subtraction, the hardest emission generation in the Powheg programs and Eq. 2.11 is not greatly different from a conventional parton shower simulation: each leg (α) effectively having associated to it its own Sudakov form factor, with R α /B in the exponent (equal to an Altarelli-Parisi splitting function in the collinear limit), and its own local definition of p T . 2 For each emission the radiating leg is, with high probability, the one with the smallest relative transverse momentum separation and the distribution carries essentially the same Sudakov suppression as in a parton shower. Furthermore, the separation of the real cross section into unique collinear singular regions, R α , and the associated damping factors, effectively produce a strong p T ordering p T,1 > p T,2 in Eq. 2.11. 3 We stress, however, that while R α is heavily damped away from the collinear region associated to α, and vanishes when other collinear configurations are approached, there is no hard cut-off, p T,1 > p T,2 , as in a parton shower, so disordered configurations can occur, albeit at a reduced rate -indeed they must occur if one is to reproduce the radiation pattern of the double emission matrix element.
While in our opinion these points, and the formulae they are based on, are very compelling evidence that the logarithmic accuracy of the Menlops simulation, defined by Eq. 2.11, should be no worse than a conventional parton shower simulation -in particular, here, we mean that the Sudakov form factor is NLL accurate -we make no allusions about the fact that we have not presented a rigorous proof of this. Hence, in what follows, we must check that this is the case and that the quality of the vector boson production Sudakov form factor is not degraded. The yellow band corresponds to the projection of the Monte Carlo statistical errors on the reference data (the first entry in the legend) onto the ratio plot.
Since we shall introduce a second, different, Menlops simulation in the next section, and in Sect. 2.4 a merging of the two, we distinguish the one in this section with a suffix ∞, referring to it throughout as Menlops ∞ , with the double-hard emission cross section dσ ∞ . The ∞ labeling is indicative of the fact that the merged samples, to be discussed in Sect. 2.4, are exclusively comprised of events from this simulation in the limit that the p T scale associated with the phase space partitioning there tends to infinity.
In Fig. 1 we present a comparison of predictions for some inclusive observables obtained with the idealized Menlops ∞ simulation described above and also the standard single emission Powheg Box vector boson production simulation, at the level of the Les Houches file events i.e. the level at which events are input to shower Monte Carlo programs. All of the distributions show a remarkable level of agreement between the Menlops ∞ predictions (cyan triangles) and those obtained with the regular Powheg Box vector production programs underlying them (red lines). This is strong evidence that the double-hard emission cross section, Eq. 2.11, has been realized in practice, substantiating our earlier claim that, by unitarity, the fully inclusive Menlops ∞ predictions, Eq. 2.11, should be basically identical to those of the underlying Powheg program.  Figure 2 shows distributions of the Z and W boson p T spectra at Tevatron and LHC energies respectively. In both cases we have focused on the low p T Sudakov region, with a mind to testing our earlier assertions regarding the quality of the resummation in the Menlops ∞ prescription (Eq. 2.11). Here, since the distribution is sensitive to multiple emission effects, and because the Powheg vector boson production cross section contains just one emission, while Eq. 2.11 contains two, to facilitate a meaningful comparison we have evolved all events with Pythia 8 to include parton shower effects (but not hadronization or multiple interactions). Again, the predictions here vindicate our claims in regards to how the Menlops ∞ procedure should preserve the quality of the NLL resummation in the underlying Powheg vector boson simulation, with the differences in the peak region being 10 %.
From a purely technical point of view, considering how the Powheg Box program distributes radiation and constructs events [8,9], one expects some shuffling of the transverse momentum spectrum associated with the single emission Powheg-V events by the secondary emission which this Menlops ∞ procedure dresses on them. 4 More to the point it is natural that this agitation acts to smear out the peak rather than sharpen it, as is seen to be the case here.
Note that the same trend is seen to basically the same extent (a depression of the peak region of 10 %) when the basic cuts on the leptons in the final state in the Rivet analysis are removed. We also note a very similar sized downward shift in the Sudakov peak is seen in the matched NNLL+NLO p T spectra of Ref. [38] with respect to the NLL+LO one (again, with no cuts applied). Although a great deal more study would be required before drawing firm conclusions in this regard, and although we in no way claim to have approached a similar level of accuracy, the changes seen in that case are intriguingly similar to those shown here in going from the pure Powheg description to this Menlops one. In fact, we believe the explanation for the effects given there applies here too. In Ref. [38] it is stated that additive higher order corrections in the intermediate-and high-p T regions, coupled with the fact that the NNLO total cross section is only 3% greater than the NLO one, serves to redistribute radiation such that the Sudakov peak decreases. 5 In our case, by unitarity, the total cross section is in fact completely unchanged by the Menlops ∞ addition, while the same higher order real corrections enter as in Ref. [38], thus one finds a similar redistribution of radiation.
From our point of view, however, the most significant point to take away from Ref. [38] is that the theoretical uncertainties due to scale variations alone, for the NNLL+NLO calculation, already bracket the 10 % effect we see here, while for the NLL+LO calculation (which we would expect to have equivalent accuracy to [8,26]) they can change the peak height by ±25 %. We also remark that when the full range of non-perturbative effects are activated in the simulation the Sudakov peak is eroded further, with the differences between the Powheg and Menlops ∞ predictions becoming negligible.
Lastly, to end the discussion on the p T spectrum, we draw attention to Fig. 3, displaying a related but more inclusive observable with which to expect agreement, namely, the integrated 0-jet rate, as a function of the k T -jet clustering scale. This quantity is closely related to the integral of the vector boson p T spectrum. Here one can see very small differences at the low end of the distribution, in keeping with those seen in the p T spectrum, however, unlike that case, the more inclusive nature of the observable (the fact that it is cumulative) allows for an averaging effect, and sure enough, the two distributions rapidly become indistinguishable like those in Fig. 1. Again, this is a strong check of the unitarity and NLO accuracy of the double-hard emission cross section and its implementation.
Turning our attention to slightly less inclusive observables in Fig. 4 we show the leading jet transverse momentum and rapidity spectra. Somewhat reassuringly the rapidity spectrum exhibits no discernible differences between the two different simulations; since it is still a rather inclusive observable and the distributions governing its generation are essentially the same in either case, the good agreement seen in that case is to be expected.
The predictions for the leading jet p T spectrum also agree very well in the region where one expects them to i.e. the low-/ intermediate-p T region. On the contrary, the distributions come apart at high transverse momentum, the Menlops simulation proving to be harder there, due to the inclusion of the doubly radiative events there [39]. Transverse momentum of leading jet dσ/dp ⊥ (jet
In figure 5 we shift focus to quantities directly sensitive to the emission of two or more jets. Since the Powheg vector boson program only generates single radiation events while the Menlops ∞ enhancement generates a further emission, in order to have a comparison here we have evolved the events with Pythia to include parton shower effects (only). As with Fig. 4 the plots shown here lend themselves to a fairly intuitive understanding. In particular, we see that since both Powheg-V and Menlops simulations include the first order real emission corrections to vector boson production, the distributions of the 0-jet to 1-jet transition rate (top left plot) given by the two approaches are in good agreement, as are their predictions for the 0-and 1-jet cross sections (bottom left plot). On the other hand, once the presence of a second hard jet is probed, as in the 1-jet to 2-jet transition rate, in the upper-right corner of the figure, and the next-to-leading jet p T spectrum, in the lower-right corner, the softer spectrum of the Powheg-V simulation becomes apparent, with the Menlops result tending to that of the Powheg-Vj programs at high p T . The deficit of hard radiation from the Powheg-Z and Powheg-W simulations in these distributions is simply due to the fact that they are reliant on the soft-collinear parton shower approximation to generate the second jet, whereas the Menlops simulation defers this dependency to the next-to-next-to-leading jet. The Menlops generator distributes the second jet according to the vector boson plus dijet tree-level matrix element, in much the same way as the Powheg-Vj programs. As we shall now go on to discuss, and as is clear from these plots, while the Powheg-Vj simulations describe multi-jet events well, they lack the resummation of Sudakov-logs required to describe the low-p T region and hence inclusive observables.  Note that the Menlops simulation arrived at here is ideal in the sense that no unphysical scales have been introduced in its formulation, in keeping with the theoretical proposal of Ref. [24], improving on the approximate implementations carried out there and in Ref. [25].
Finally, before leaving this section, we wish to comment on the assignment of the colour structure and interfacing to the parton shower. While the flavour structure and Φ V j are read in from the Les Houches event file, the colour flow is assigned to the event in the same way as in the unmodified Powheg-Vj program. This is a two stage process. Firstly, the large-N C planar cross sections associated to the hh → V j state are evaluated at Φ V j and a flow assigned with a probability proportional to its cross section. Secondly, the colour structure for the subsequently emitted radiation, parametrized by Φ j 2 and α, is added as if it were a 1 → 2 parton branching from the leg of the hh → V j state associated to α. This is also the method of colour assignment used in Sect. 2.3, the only technical difference there being that the flavour structure and Φ V j kinematics are generated by the program itself, as opposed to being read in from a Les Houches file.
Note that a more careful treatment of the colour structure would include truncated showering effects, the absence of which has been suggested as a potential cause of spurious high-mass clusters in hadronization [5]. While we have not investigated such effects here, we have checked that at the hadron level we do not observe any of the distortions in differential jet rates seen in Ref. [5]. We suggest that this may be due to the fact that, as well as being a very much suppressed effect, the matrix element-parton shower interface occurs at a fixed scale in Meps merging schemes like that studied in Ref. [5], whereas in Powheg the analogous scale is variable, being set by the transverse momentum of the emitted radiation. We stress that the issue of truncated showers does not arise through any of the Menlops themes in this article but stems from the Powheg method in general, moreover, all studies to date in this context suggest that its effects are negligible [40][41][42]. In any case we note that such effects could be included for the processes studied here by modest extensions of the methods in Refs. [40,43].

Menlops improved simulation of vector boson plus jet production
In this section we aim to promote the Powheg simulation of vector boson plus jet production to a Menlops one. Whereas, up to now in the literature, the Menlops acronym has stood for Nlops simulations enhanced by higher multiplicity tree level matrix elements, now we extend it to cover also the inclusion of lower multiplicity ones; in this case by resummation and soft-collinear factorisation. By way of distinguishing the simulation here from that of Sect. 2.2 we append the suffix 0, referring to it as Menlops 0 with a double-hard emission cross section denoted dσ 0 . The naming derives from the fact that the distributions resulting from this simulation correspond to the limit where the merging scale to be introduced in Sect. 2.4 is set to zero.
The cross section for vector boson production develops large logarithms, ∼ log Q 2 /p 2 T , at all orders in perturbation theory, as the boson's transverse momentum becomes small, Q being the off-shell boson mass. In this region fixed order perturbation theory then becomes unstable, with leading order and next-to-leading order predictions for the vector boson p T spectrum diverging to plus and minus infinity respectively, as p T → 0. While the Powheg-Vj programs exhibit precisely the latter unphysical behavior, the same large logarithms are summed to all orders in α S in the Powheg-V programs, leading to a stable and physical description of the Sudakov region and, hence, also inclusive observables.
Thus, here we begin by clarifying the differences between the Powheg-V and Powheg-Vj programs in the low transverse momentum region, in order to better understand when the NLO computation of vector boson plus jet production can be expected to begin to break down and with a view to 'mapping' the behavior of the former onto that of the latter when this occurs. We point out that while the role of the renormalization and factorization scales may appear prosaic, it is central to the discussion here, hence, it is often referred to in the following. In Ref. [26] the authors established the equivalence of the resummed vector boson transverse momentum spectrum of Ellis and Veseli [44] (based on the DDT formalism [45]) and the Powheg hardest emission cross section Eq. 2.1. Specifically, it was shown how the Powheg Sudakov form factor for divector boson production relates to the next-to-leading log transverse momentum space Sudakov form factor in Ref. [44]. The same proof trivially holds here for vector boson production since the final-state is still colorless. In particular, neglecting terms beyond NLL accuracy it was proven that wherex and Q are evaluated in terms of the underlying Born kinematics, Φ V , and T NLL is the exponent of the next-to-leading log DDT Sudakov form factor with C F and C A the usual Casimir factors. The function A (α S ) is defined as (2.14) Substituting the expression for the Sudakov form factor in Eq. 2.12 into the Powheg-V hardest emission cross section, Eq. 2.1, and noting that the factorization scale there in B is Q, while in R/B it is p T , all PDF factors in the combinationB/B × ∆ cancel. Hence, to NLL accuracy one has where we-reiterate that the renormalization and factorization scales in R here are set to the p T of the W / Z boson. 7 We also remark that in order to achieve NLL accuracy in the Powheg Sudakov form factors α S is replaced by A (α S ) in the real cross section therein and, hence, in Eq. 2.15 [8,26]. Lastly, we note that in writing Eq. 2.15 we have taken the ratioB/B as being equal to that of the associated PDFs, which is a legitimate replacement at this level of approximation [26].
Using unitarity to integrate out the Φ j 2 phase space, the analogous expression to Eq. 2.15 in the vector boson plus jet Nlops is simply 14 GeV. The plots of the vector boson p T spectra following this discussion are consistent with this estimate. Let us make clear that this bound is lower than that which would be obtained were different renormalization / factorization scale choices made in Eq. 2.16 (see e.g. Fig. 1 of Ref. [38]) i.e. logarithmically enhanced corrections are already in effect to some extent above this threshold but they are offset by the resummation implicit in the scale choices.
Explicitly, we have takenp T = 14 GeV in Z production andp T = 12.5 GeV in W production throughout. While the values inferred from the aforementioned inequality should give a conservative estimate forp T , to put things on a more solid footing, one should ideally check thatp T is chosen above, yet close to, the point where the fixed order and resummed spectra just start to diverge. The values quoted here fulfill the latter criterion (see Fig. 7); in practice they should be adequate for all vector boson production applications.
Belowp T missing logarithmically enhanced terms in the DDT Sudakov form factor begin to render the NLO corrections to vector boson plus jet production invalid by themselves. Just abovep T the combination of the NLO corrections and renormalization and factorization scale choices serve to well model the onset of large logarithmic Sudakov effects -we take the fact that the Powheg-V and Powheg-Vj transverse momentum spectra in Fig. 7 agree to within 5% in this region to be good evidence of that.
Having noted the various scale dependencies and the projected domain of validity of the vector boson plus jet NLO computation, our Menlops strategy for enhancing the Powheg-Vj simulations reduces to cutting off the NLO contributions toB(Φ V j ) where they start to become compromised by the neglect of large logarithms at higher orders, p T,1 p T , and instead resumming those same Sudakov logs, in precisely the same way as is done in the Powheg-V programs. Thus, we write the modified Menlops 0 cross section as where P(p T,1 ) is the following switching function, 8 dσ V j is the unmodified NLO hardest emission cross section of the Powheg-Vj programs (Eq. 2.8) and dσ NLL V j a reweighting of dσ V j , at leading order (B(Φ V j ) → R(Φ V j )), intended to replicate the Φ V j dependence of the Powheg-V programs: This being so, when evaluating the leading order single emission cross section, R(Φ V j ), in Eq. 2.19, we replace α S → A (α S ), as when generating radiation in the Powheg-V programs. Furthermore, as in that case, R(Φ V j ) is multiplied by a ratio of PDFs, in the form Note that the latter ratio comprises a tower of large logarithms, the same large logarithmic content as the ratioB/B in the conventional Powheg-W / -Z cross section formulae. While we do not claim that the description of fully inclusive observables is better than leading order here, given the guiding principle amounts to trying to match, as far as practically possible, the Φ V j dependence of dσ V , we can at least correct the overall normalization using a K -factor, K =¯ B / B , for the vector boson production process (determined with µ R = µ F = m V throughout). Lastly, we include the Powheg Sudakov form factor ∆(Φ V , p T,1 ) leading to a double-hard emission cross section much like that of the Menlops ∞ simulation (Eq. 2.11), with the vector boson production hardest emission cross section recovered up to the replacementB → K B on integrating out Φ j 2 . All of the modifications and reweightings above are carried out by straightforward modifications to theB(Φ V j ) function, except for the addition of the Sudakov form factor. With the Monte Carlo generation cut on p T,1 set to that defining unresolvable radiation in the Powheg Box programs (p min T = 0.8 GeV), the modified event generation process including Sudakov form factor suppression proceeds as follows: 1. A vector boson plus jet configuration, Φ V j , is generated according to Eq. 2.17 but with the term ∆ (Φ V , p T,1 ) in dσ NLL V j replaced by one. 8 The width of the switching threshold is determined by , which we have taken to be 2.5 GeV throughout -the fact that the switching threshold has a finite width allows for small O (1%) differences between dσ NLL V j and dσV j , in the vicinity of p merge T , to be smoothed out.

2.
A random number R is generated and if R > P(p T,1 ) the event generation continues as in the unmodified Powheg-Vj program at step 7.
4. Kinematics for vector boson production are assembled from Φ V .

The
Powheg-V radiation generating code is used to generate a 'ghost' emission, Φ j 1 , with respect to Φ V .
, the vector boson plus jet configuration, Φ V j , is rejected and the event generation is restarted at step 1. 9 7. Radiation is generated from the Φ V j state exactly as in the Powheg-Vj program, completing the Powheg event.
Although generating the vector boson plus jet configurations at such low values of p T and then vetoing them in large numbers may seem computationally expensive, the actual efficiency associated with the rejection procedure above (with the generation cut set to 0.8 GeV) is not so costly, the acceptance rate being a little under 20%. When all modifications are considered together, the full event generation procedure actually proceeds much faster than normal. The switching off of the NLO corrections, when the p T of the underlying Born configuration dips below ∼ 14 GeV, means one only needs to evaluate relatively simple matrix elements there, moreover, the fact that the distribution is then manifestly positive in this region, greatly reduces the number of negative weight events and hence the number of folded-integrations [14,46] required to remove them. We have used a single folding of the real emission phase space Φ j 1 and Φ j 2 (a so-called 2-1-1 folding in the notation of Ref. [14]), simply out of a desire to test the programming apparatus under general conditions, rather than any inclination to remove negative weights. In the event samples created with the unmodified Powheg-Vj programs used in this study, the fraction of negative weight events was found to be ∼ 20 %, while in all other Powheg-V and Menlops samples it was below ∼ 0.5%. Although the computational speed here proved to be more than adequate for the purposes of this study, it is clear that the efficiency of the Sudakov suppression procedure can be improved without too much effort.
We remark that once thep T threshold is crossed and the NLO corrections inB(Φ V j ) are turned off, the basic idea behind the procedure here is really no different to the Sudakov suppression applied to matrix elements in the CKKW(-L) and MLM tree-level merging methods. On a basic level, the whole approach here is a direct application of those same principles.
Unlike the general tree-level merging techniques, the Powheg formalism, since it aims at NLO accuracy, is defined and realized with a high degree of technical precision, affording 9 Given the interpretation of the Sudakov form factor as a no-emission probability, events with kinematics ΦV , pT,1 then pass this veto with probability ∆ (ΦV , pT,1), thus generating the Sudakov suppression term that was set to one in the 'crude' distribution in step 1. more control over the implementation here. Note, in particular, that the inverse mapping of the phase space, Φ V j → {Φ V , Φ j 1 } , is clearly and unambiguously defined; this follows from the fact that the mapping covers the full single emission phase space, without dead zones, without the need for any kind of emitter-spectator assignment, and because it is, in this case, uniquely attributable to initial-state radiation. The clear specification of the formalism, combined with this phase space inversion, underpins all of the modifications described here, ultimately allowing us to reconstruct, up to the replacementB → K B, exactly the same distribution for Φ V and Φ j 1 as in the Powheg-W and Powheg-Z programs. This exactness of the implementation should reveal itself in a number of ways, as in the case of Sect. 2.2. In particular, predictions for inclusive W / Z production should closely trace their genuine NLO Powheg-V counterparts, the only source of differences beingB → K B, also the large-log resummation in regions sensitive to Sudakov effects should maintain its quality.
For completeness we note that in the Powheg-V simulations some events occur with non-radiative kinematics, Φ V only, being physically interpreted as ones in which the emitted radiation is unresolvable and in a region dominated by non-perturbative dynamics. These events are neglected in Eq. 2.19, comprising just a fraction ∆(p min T ) < 1 % of the total sample in the Powheg-V programs -a negligible contribution given that the accuracy aimed at for fully inclusive observables in this section is just leading order.
As in section 2.2, here we begin to demonstrate that the Menlops 0 method has been faithfully implemented according to Eq. 2.19 by first considering fully inclusive observables in Fig. 6. Immediately one can see, as expected, that the vector boson plus jet production programs give markedly different predictions to the Powheg-V and Menlops 0 simulations. This should not come as any surprise, since the behavior of the NLO calculation underlying the former becomes unphysical when the observables include contributions from the Sudakov region. Recall that the vector boson plus jet programs include an unphysical cut-off on the p T of the underlying Born kinematics, to screen the associated collinear singularity, without which event generation is not possible. We have taken this Monte Carlo generation cut to be 5 GeV here. Although we have not shown it explicitly, taking a different value of the cut for this type of observable will produce significant variations in the results here, in fact by making the cut sufficiently small (2 − 3 GeV) all of the predictions become negative [14]. On the contrary, since the large logarithms associated with the same Sudakov region are resummed in Eq. 2.19, the Menlops 0 formulation exhibits no such unphysical behavior. In fact, in this case, the technical cut has been taken down to that used in generating radiation in the vector boson production programs, 0.9 GeV, moreover, the dependence of the predictions upon the cut is physical. Pleasingly, we see that the level of agreement in the predictions of the Powheg vector boson production and the Menlops 0 simulations is by-and-large very good, with the two rarely deviating by more than 10%; recall that the overall normalization is to a large extent fixed to the NLO total cross section by the constant factor, K, in the low-p T region. Continuing as in section 2.2, we next look to check the implementation of resummation in the Sudakov region, where the transverse momentum of the vector boson becomes small in Eq. 2.19. As noted earlier, the similarities of Eq. 2.19 with the double-hard emission cross section in Eq. 2.11 of Sect. 2.1 are apparent, moreover, the surrounding discussion on resummation there applies here too, without modification. Here again, in figure 7, one can see good agreement between the Menlops 0 and Powheg-V predictions, with deviations being limited to about 10%, which we take as confirmation that the accuracy of the Sudakov resummation in either case is the same.
The Powheg-Vj prediction is shown in blue in Fig. 7, the fact that it fails to describe the low p T end of spectrum is no surprise, as we have already cited it as a motivating factor for the modifications proposed here. We have chosen to show it since we believe it substantiates our earlier statements regarding when the NLO calculation can be expected to diverge from the resummed predictions, ∼ 14 GeV. Note that this value for when one can expect the resummed and fixed order predictions to exhibit deviations is noticeably lower than that which one might have inferred from Fig. 1 of Ref. [38]. Note, however, that in the case of Ref. [38] the renormalization and factorization scales are set to the Z boson mass, whereas the Powheg-Vj program uses the p T of the Z boson in the underlying Born configuration and is thus in closer correspondence with the DDT formalism.  Figure 8: Menlops 0 (magenta) and Menlops ∞ (cyan triangles) predictions for the low transverse momentum region in Z production at the Tevatron (left) and W production at the LHC (right). The results shown here include the effects of multiple emissions through parton showering; the simulation of hadronization and multiple parton interactions has been omitted. These plots serve to prove that the non-trivial implementation of Sudakov suppression effects has been carried out properly in the Menlops 0 case.
In order to prevent the distributions obscuring one another and due to its importance in the context of the modifications being discussed in this section, in figure 8 we show the same plots again but this time comparing the Menlops 0 and Menlops ∞ output. If the hardest emission cross section formulae have been implemented precisely according to Eqs. 2.11 and 2.19 the two sets of predictions should be essentially identical in the peak region, as is seen to be the case.  Figure 9: Powheg-W (red), Powheg-Wj (blue) and Menlops 0 (magenta squares) predictions for the leading jet transverse momentum and rapidity spectra for current LHC pp beam energies.
In figure 9 we compare the predictions of the leading jet transverse momentum and rapidity distributions among the Powheg-W, Powheg-Wj and Menlops ∞ approaches. As expected, since these distributions are inclusive with respect to the jet-associated vector boson production process and since they do not receive contributions from regions of phase space in which the W boson has a low p T , the Powheg-Wj and Menlops 0 are indistinguishable from one another. In the latter two simulations a small 10 − 20 % NLO enhancement of the total vector boson plus jet cross section is discernible from the lowp T end of the transverse momentum spectrum and again, more clearly, in the rapidity spectrum. The jet transverse momentum spectrum clearly shows the NLO enhancement increasing with the p T [39]. Recall that a similar but slightly less robust, enhancement was visible in the leading jet transverse momentum spectrum in the Menlops ∞ case (Fig. 4). This is not unexpected since in that case the simulation should offer the same leading order description of two-jet events (whose rate is underestimated by the single emission Powheg-V programs) which serve to harden the leading jet p T spectrum. Figure 10 focuses more on observables sensitive to the emission of multiple hard jets. The trends shown here are readily understandable in terms of the basic features and frailties of the various generators. In the upper right-hand plot the differential 0 → 1-jet rate is presented, physically representing the exclusive k T -jet clustering scale at which a 0jet event becomes resolved as a 1-jet one. Here one can see good agreement among all three approaches except in the low transverse momentum / Sudakov region where the Monte Carlo generation cut and the lack of any large-log resummation in the Powheg-Zj program becomes evident with respect to the other two resummed predictions. The same feature is apparent in the differential 1 → 2-jet rate, for the same reason. In this case, however, at the high p T end of the spectrum, the NLO and partial NLO corrections to the jet-associated vector boson production cross section reveal themselves, leading to an excess of the Powheg-Zj and Menlops 0 results over, what is for this observable, in this region of phase space, only a LO prediction from the Powheg-Z program. The lower plots, displaying the exclusive jet multiplicities and transverse momentum spectrum of the next-to-leading jet, show the now-familiar picture of the Powheg-Vj and Menlops 0 simulations predicting more radiation in the high transverse momentum regions of phase space, where, by construction, they agree exactly with one another. Before moving on though we draw attention to the one place where these two predictions do not agree in the lower two plots viz the 0-jet multiplicity bin; there the Menlops 0 prediction correctly interpolates to that of the vector boson production Nlops.   Figure 10: Observables sensitive to the emission of at least one hard jet. Powheg-V, Powheg-Vj and Menlops 0 results are shown in red, blue and magenta respectively. As in all earlier plots, the sub-plot displays the ratio of the various predictions with respect to that given by the first entry in the legend; the inset yellow band corresponds to the statistical errors on the reference data.

Merging Menlops samples
The Menlops enhancement described in Sect. 2.2 augments the Powheg Nlops simulation of vector boson production freely with next-to-leading order real emission matrix elements, while that of Sect. 2.3 includes, through the factorisation theorem and Sudakov resummation, lower order matrix elements in the vector boson plus jet simulation. Neither simulation therefore includes a full set of NLO corrections to both vector boson production or jet-associated vector boson production.
A practical, if seemingly crude, way to improve on the predictions of either Menlops simulation, is to simply combine the event samples they produce appropriately, such that results for both fully inclusive observables and also observables inclusive with respect to the vector boson plus jet final-state, are accurate at the NLO level. To this end we follow a similar line of reasoning to Ref. [24], in particular, the discussion in Sect. 3.
In Ref. [24] an approximation to the exact Menlops method is implemented where events from Nlops and Meps samples are filtered according to their jet multiplicity, us-ing the exclusive k T -jet clustering algorithm, at the so-called Menlops merging scale, and arranged in 0-jet, 1-jet and 2-jet sub-samples. The Nlops sub-samples are generally considered to offer a better description of the 0-jet and 1-jet events than their Meps counterparts since both simulations possess the same tree-level matrix elements, yet, at least for the 0-jet class, the Powheg simulation also includes virtual corrections. Conversely, the Meps simulations contain higher multiplicity matrix elements that the Nlops simulations do not, so the 2-jet sample is regarded as having a better description there. Thus, in the pragmatic implementation of Ref. [24], the final Menlops samples are constructed by complementing events from the 0-and 1-jet Powheg sub-samples with those in the 2-jet Meps one, in carefully specified proportions; the key consideration being that the unitary violating Meps component should be restricted so as not to diminish the NLO accuracy of inclusive observables.
Here the situation is somewhat analogous. Rather than implement a multi-particle phase space partition according to a conventional jet measure, we use the transverse momentum of the vector boson plus jet kinematics, at the level of the bare (pre-shower) events. As noted in Ref. [13], in principle one can assemble an infrared safe jet algorithm from the radiation phase space mappings underlying the Powheg Box implementation [9], hence the partitioning here is in effect much like that of a Durham-type jet algorithm. From a theoretical point of view though, it is more appealing to base the partition on the underlying Born kinematics, since these are more directly interpretable in terms of the hardest emission cross section formulae, Eqs. 2.11, 2.19. To create the doubly NLO accurate event sample we then populate the region p T,1 < p merge T using the Menlops ∞ program and the region above p merge T using the Menlops 0 simulation -we adjust the number of events taken from either sample appropriately, to relieve (small) differences in their event weights, giving uniform ±1 weighted events in the final sample.
Plainly, with no modifications to the component Nlops programs, the dependence on the unphysical phase space partition is such that as it tends to high values one recovers the single vector boson Nlops results, having, we note, only a leading-log, soft-collinear, approximation for the next-to-leading jet. Conversely, as the p T defining the merging scale tends to low values, the associated vector boson production Nlops is recovered, for which the description of fully inclusive observables is (negatively) divergent and unphysical. Hence by merging instead the output of the Menlops simulations we mitigate the unwanted scale dependence to the same extent as in the CKKW method.
In order to preserve the NLO accuracy of inclusive observables in Ref. [24] it was important that the merging scale be bounded from below, to prevent too much of the tree level Meps component entering and spoiling it, specifically, it was stipulated that the Meps events form no more than a fraction O(α S ) of the total. Equally, in the present case, while the merging scale dependence is lessened by combining Menlops improved samples, if we require the sample to be doubly NLO accurate, in the sense indicated by the second paragraph in this section, it must be subject to restrictions. Here, as the merging scale is decreased towards zero, the Menlops 0 sample is recovered, which is strictly only LO accurate in the description of fully inclusive observables (modulo the NLO normalization K-factor), hence, the key point is again to bound the merging scale from below. Following the same logic as Ref. [24], it is clear that NLO accuracy for inclusive observables will be preserved, up to fluctuations of NNLO significance, if one restricts the fraction of Menlops 0 events in the sample to be less than O(α S ).
In practice, it is almost certainly the case that very much bolder merging scale choices, taking in a much greater Menlops 0 fraction would in fact not cause deviations in the predictions for inclusive observables. This point is clearly made in Fig. 6, where there is a remarkable level of agreement between the latter and the genuinely NLO Powheg-V predictions. This is not something to be deceived by though. Whereas the Menlops 0 kinematics, Φ V , there are distributed according to a K-factor enhanced leading order distribution, B(Φ V ), the Powheg-V predictions include full NLO correlations throughB(Φ V ). Therefore, while it may be the case that the combination of the K-factor and B(Φ V ) work well, a more complete study than this ought to expose the leading order sensitivity to changes in the renormalization and factorization scales. In order to have a truly doubly NLO accurate sample the lower bound on the merging scale should be respected.
We have enforced this criterion in building our merged samples: in the case of W production at the LHC, merging the two Menlops samples at a scale p merge T = 35 GeV, leads to a fraction of 11 − 12 % of Menlops 0 events in the final sample 10 and in the case of Z production at the Tevatron p merge T = 25 GeV leads to the same amount. Note that in combining the samples, obviously the Menlops ∞ and Menlops 0 cross sections are preserved in the regions either side of the cut, thus the final total cross section can differ from the NLO fully inclusive Powheg-V / Menlops ∞ results. In all cases we find the difference between the Menlops ∞ result and those of the merged samples -which we refer to as Menlops 25 and Menlops 35 on account of the merging scales used -is not greater than 1 %. This latter observation is unremarkable given the fraction of Menlops ∞ events allowed and the agreement shown in Fig. 6 for inclusive quantities, we merely note it for technical reasons.
Furthermore, in practical applications, provided that the vector boson p T is used as the renormalization and factorization scale in the fixed order components, the latter constraint renders the added resummation in the Menlops 0 simulations pointless; so long as the merging scale is not too small one can obviously use the jet-associated vector boson production program, unmodified, in the high p T region. We therefore employ the Menlops 0 simulation in merging more as a matter of theoretical correctness. In this regard, when considering the dependency of the predictions on the unphysical merging scale, the Sudakov suppression factors take on their practical as well as theoretical significance, in the same way as in the CKKW method. Actually the merging scale dependence here is substantially reduced with respect to the CKKW case since, as we shall see, the two component sub-samples are exactly unitary, moreover, they both include precisely all the same tree-level matrix elements and exact same resummation. In particular, from the point of view of the region below p merge T , we expect that using the Powheg-V simulation instead of its Menlops enhanced version would be notable as a marked dependence on the merging scale in multi-jet observ- 10 Setting p merge T = 20 GeV in LHC W boson production yields a sample 25 % of which is Menlops 0 events; note that the figures in section 3 indicate that inclusive predictions are rather insensitive to maximal merging scale variations (as does Fig. 6).
ables. As noted in the introduction, to lessen the dependence on the merging scale beyond this, would seem to be equivalent to the task of achieving genuine NNLO-parton shower matching, something which is beyond the more modest, practical, aims of this work.
A final, related, technical point concerns the presence of discontinuities in the differential jet rates and p T spectra. We claim a level of matching at least as good as in the CKKW method, in particular, that all logarithms pertaining to p T,1 and hence p merge T are resummed. On the other hand, as is by now clear, the same is not true for the finite, non-logarithmic, NLO (virtual) corrections. This alone defines the level of ambiguity in the present merging. Below p merge T we have full NLO corrections to vector boson production observables only, while above it there are only full NLO corrections to vector boson plus jet production. Of course, one expects and finds that the NLO corrections to the latter enhance jet rates and p T spectra with respect to the former by O(α S ) at intermediate / high p merge T -these are, after all, the differences we hope to take into account by merging -leading to the possibility that small discontinuities may occur around p merge T on merging the two samples. Note, however, that the exact value of the merging scale should not be taken too seriously, it is wrong, for instance, to imagine that the quality of either component sample instantly degrades either side of p merge T ; there is certainly some blurring in that regard. This being the case, to smoothen any small NLO differences occurring between the different jet rates in the vicinity of the merging partition, we do not implement an immediate step-function cut-off at p merge T , but rather a smooth sinusoidal damping function of the form Eq. Transverse momentum of leading jet dσ/dp ⊥ (jet  Figure 11: The 0-to 1-jet differential jet rate (left) and leading jet transverse momentum spectrum (right), in W boson production at the LHC. In each figure we show, as green crosses, the prediction from the merging of the component Menlops 0 and Menlops ∞ samples, with the merging scale set to p merge T = 35 GeV. The magenta lines and cyan triangles, respectively, show the corresponding predictions from the component Menlops 0 and Menlops ∞ samples. This collider setup and these distributions were selected on the basis that they show the matching ambiguity between the two samples in its worst light, specifically, the difference between the magenta and cyan distributions at p merge T .
By way of showing the worst that can be expected in terms of such discontinuities at the merging scale, we show in figure 11 the 0-to 1-jet differential jet rate and the leading jet p T spectrum, which we found to be most sensitive. Here one can see the Menlops 0 and Menlops ∞ predictions disagree at the level of ∼ 5 % in the vicinity of the 35 GeV merging scale, with the Menlops 35 merged prediction interpolating between them. Due to the size of the differences probed here the Monte Carlo statistical errors make the interpolation difficult to see, we suffice to say that in this regard we have used here event samples comprising of well over one million events. Other distributions such as the vector boson p T spectrum exhibit less sensitivity / ambiguity in the region of the cut, moreover, the magnitude of the NLO corrections at the Tevatron is smaller, hence, there Menlops 0 and Menlops ∞ results are more-or-less indistinguishable in the low p T region near p merge T . Many other distributions from all of the Menlops samples will now be studied in comparison to data in Sect. 3.

Tevatron and LHC predictions
Here we confront the native Powheg Box programs, their Menlops enhancements and the merged Menlops samples, with Tevatron and early LHC data. The analyses presented in this section are therefore carried out having developed the events output by the Powheg and Menlops simulations to the hadron level and included multiple interaction effects. To this end we have used the Pythia 8.150 program, following closely the recommendations set out in Ref. [31] but otherwise maintaining the default settings. 11 We have interfaced Pythia to the HepMC package [47] and hence on to the Rivet framework [29] to perform the experimental analyses. The selection of results is intended to be representative.
We have endeavored to keep the colouring of the predictions the same here as in earlier plots in the article. As in all preceding histograms, each sub-plot shows the ratio of each type of prediction to that placed first in the legend, here the experimental data. Also, as before, in each case the shaded yellow band depicts the projection of the errors on the reference data into the same ratio distribution.

Comparison to Tevatron data
In figure 12 we show the Z boson rapidity spectrum as measured by the CDF collaboration at the Tevatron [48]. On the left Powheg-Z, Powheg-Zj and merged Menlops predictions are compared to data, while on the right the Menlops 0 and Menlops ∞ component samples are shown along with the same merged Menlops 25 sample (p merge T = 25 GeV). Given the inclusivity of the measurement here (Fig. 12) the Powheg-Zj prediction unsurprisingly fails to describe the data; we present it for completeness and to illustrate the improvement rendered to it by the Menlops 0 enhancement described in Sect. 2.3, as shown in magenta on the right-hand plot. The Powheg-Z and merged Menlops sample predictions (red line and green crosses) offer a good description of the shape of the distribution while their overall normalization is negatively offset from the data by about 10 %, an acceptable disagreement for an NLO computation and one echoed by the other, independent, theoretical predictions given in Ref. [48]. From the point of view of validating the Menlops procedure, we draw attention to the fact that the result of merging the Menlops samples is in complete agreement with the Powheg-Z one for this fully inclusive observable. From the same perspective, the stability of the predictions against extremal variations in the merging scale, in the right-hand plot, is remarkable. While the level of agreement between the three Menlops predictions is striking we do not wish to over emphasize it; by virtue of the rules applied to the composition of the merged sample, the Menlops ∞ and Menlops 25 results should be in near perfect agreement, whereas, in principle, one could have expected some marginal deformations of the shape of the Menlops 0 result with respect to them, similar in size to those seen between NLO and K-factor adjusted LO predictions. We cannot claim that the excellent level of agreement between the Menlops 0 results and the other two Menlops ones is not, to some extent, coincidental -as noted in Sect. 2.3, the true leading order nature of the shape of fully inclusive Menlops 0 predictions would reveal itself as a sensitivity to renormalization and factorization scale variations, while the genuine NLO quality of the other two Menlops samples would show them to be suitably insensitive.
The two key points to take away here are: firstly, that the Menlops 0 enhancement of the Powheg-Zj program means that it now gives a physical prediction for this inclusive observable where it did not before, moreover, one that has an uncanny resemblance to the fully NLO predictions; the second important point is that the Menlops ∞ and Menlops 25 results are in complete agreement with that of the NLO Powheg-Z program, from which they should have inherited this prediction. Figure 13 shows the same set of comparisons with respect to a recent DØ measurement [49] of an inclusive quantity, φ * η , said to be closely related to the Z boson transverse momentum: φ * η = tan (φ acop /2) sin θ * η , with φ acop defined as π minus the azimuthal separation of the muonic decay products and θ * η a measure of the scattering angle of the dimuons with respect to the proton beam in their rest frame. From the definition one can see that when the final-state leptons are back-to-back in azimuth, as would be the case if there were no additional radiation in the final state, φ * η tends to zero. Hence the region φ * η = 0 is highly correlated with the Sudakov region in the Z boson transverse momentum spectrum.
It follows that, despite the plots being normalized to unit area, the Powheg-Zj prediction undershoots the data for φ * η → 0, while the others offer reasonable agreement there. The fact that this area of the plots is associated with the Sudakov peak region also accounts for the small experimental errors there. By-and-large all of the conclusions drawn in regard to the Z boson rapidity spectrum, in particular the Menlops predictions, apply here unchanged. What is worth noting in addition, is the binning of the prediction according to the rapidity of the Z boson. The Sudakov form factor, amongst others resumming large-logs in the Menlops 0 improved Powheg-Zj simulation, depends on this, hence, it is interesting and reassuring to see that the low φ * η region is described equally well by this program in the |y Z | < 1 bin as in the 1 < |y Z | < 2 one.  MC/data Figure 15: On the left we compare Powheg simulations of vector boson production and vector boson plus jet production, as well as a Menlops prediction, to a recent DØ measurement of the Z boson p T spectrum, in the dimuon channel [51]. The Menlops event sample was formed by merging Menlops 0 and Menlops ∞ samples at an underlying Born p T of 25 GeV. The right-hand plot shows the maximum variation in Menlops predictions associated to changing the merging scale. The yellow band signifies the projection of the errors on the reference data -the first entry in the legend -into the ratio MC / data.
In figures 14 and 15 we show the transverse momentum spectra of W and Z bosons as measured by the Tevatron DØ collaboration in Refs. [50] and [51] respectively. In all cases the agreement of the various Monte Carlo predictions is quite pleasing, the only exception being the Powheg-Wj and Powheg-Zj programs' failure to describe the low p T region. As can be seen on the left-hand side of each of the figures, this problem is remedied well in the Menlops 0 enhanced versions, which include the appropriate resummation of large Sudakov logarithms there (magenta lines). Once again the agreement between the merged Menlops prediction (green crosses) and those of its two component samples is very good.
As with the case of the Menlops 0 result for the Z boson rapidity spectrum, the fact that the Menlops ∞ simulation replicates the high-p T behavior of the Powheg-Vj and other Menlops predictions in this region is noteworthy, since its description there is only partially next-to-leading order -through the inclusion of the double real emission matrix elements and associated Sudakov suppression -the others include full NLO corrections to this observable throughout the high p T region. While it is the case that the high-p T description of the Menlops ∞ program appears in good agreement with its fully NLO counterparts, one expects that it exhibits a little more renormalization and factorization scale sensitivity there.
Here the first main point is that the Menlops results give an improved description over the full p T range with respect to their unenhanced versions, in particular the description of the low-p T region is always physical and the predictions exhibit no more variation among the improved programs than in the originals (which, for this observable, was already quite mild). The second important point is that the merged Menlops 25 sample can be considered to give a genuine NLO+NLL accurate prediction all through the spectrum and that here the sensitivity of the prediction to the exact choice of merging scale in a 'real-life' application should be practically negligible. Figure 16 displays CDF measurements of the 1-, 2-and 3-jet cross sections, as well as inclusive jet p T spectra underneath. The various Monte Carlo simulation predictions are relatively straightforward to understand in terms of their constituent matrix elements. Since all of the observables demand the presence of at least one jet, the Powheg-Z predictions obviously suffer with respect to all of the others, since it includes only a tree-level single emission matrix element for Z production, with additional radiation being generated according to the parton shower approximation. In particular, the Powheg-Z predictions have a characteristically softer radiation pattern associated with them, manifested here as diminished multi-jet cross sections and inclusive jet p T spectra. It is then pleasing to see that the Menlops ∞ improvement (Sect. 2.2), which generates a secondary emission, according to the exact double emission matrix element, from the Powheg-Z single emission configurations, gives a set of predictions in better agreement with experimental data and the other approaches. On a related note, one can also see the general trend persisting from the previous figures, whereby the Menlops-improved results are in closer agreement with one another than the underlying Powheg programs.  While the general agreement of the improved approaches with the data in Fig. 16 is good, clearly, in regards to the analysis carried out there, all of the Menlops approaches offer no practical improvements over the description afforded by the Powheg-Zj program; the low p T region is excluded by the nature of the observables, hence, there is no need for Sudakov resummation (provided optimal scale choices are made) or NLO corrections to the underlying inclusive Z production process. Thus, with this unique application in mind, the advantages of Menlops event samples are rendered moot, only becoming useful when viewed as a tool to perform many different analyses.Alternatively, more importantly, when Z production is being considered as a background rather than a signal process, the distinction between the inclusive and associated Z plus jet production processes is not often clear and generally applicable tools are required. It follows that the ultimate practical significance of these plots (and in fact all plots in this section) is in verifying the ability of the Menlops approach to capture the best description from the Nlops simulations throughout phase space and to show that the sensitivity of merged Menlops predictions to the merging scale, p merge T , is remarkably mild. In particular, comparing the left-and right-hand columns, that the merging scale sensitivity is minimized when the component sub-samples are made from Menlops improved programs as opposed to unimproved ones.  Displayed at the top of figure 17 are distributions of the azimuthal separation of the Z boson and leading jet (above) in the muonic decay channel at the Tevatron [53], as well as the corresponding rapidity separation (below). Here the description of the data is byand-large fair in all of the approaches, the most prominent source of discrepancy being in the region corresponding to a large azimuthal separation of the leading jet and the Z boson in the Powheg-Z case. Again, as in the discussion surrounding Fig. 16, this is somewhat expected due to the fact that this simulation resorts to a soft-collinear parton shower approximation to generate events there, while the other methods make use of the tree-level double emission matrix element. In connection with this, one can again see a tendency for the Menlops predictions (right) to be in closer correspondence with one another than the unimproved Powheg ones (left), re-asserting that the p merge T dependence is minimized by building the merged sample from Menlops improved Powheg simulations rather than unimproved ones.

Comparison with LHC data
In figures 18 and 19 we compare Powheg-W, Powheg-Wj and Menlops predictions to Atlas collaboration measurements of the N -jet cross sections (0 N 4), as well as the leading and next-to-leading jet transverse momentum spectra, in the electron and muon decay channels respectively [54].
Regarding the leading and next-to-leading jet p T spectra, all of the approaches, corrected and uncorrected, offer a reasonable description of the data. The soft-collinear parton shower approximation for the second jet in the Powheg-W simulation does not reveal itself to the same extent as in Fig. 16 -this should be exposed by greater statistics in future experimental data sets -nevertheless, one can see a systematic tendency for the Powheg-W program to undershoot the experimental measurements in the higher p T bins. On the other hand, the relative softness of the radiation pattern in the Powheg-W program is much more in evidence in the N -jet cross sections where it can be seen to greatly underestimate that of higher multiplicity events. As in Figs. 16, 17 this deficiency of the Powheg-W simulation is cured in the Menlops 0 enhanced version (magenta lines).
From the point of view of the advantages to be had by the Menlops approaches here, considering just the jet p T spectra, there is not much to be gained in practical terms with respect to the Powheg-Wj simulation; the discussion given earlier surrounding figure 16 largely applying here too. On the other hand, looking also to the N -jet cross sections one can see, as expected, the Powheg-Wj program fails due to its unphysical description in the region where the vector boson p T is small, thus here the Menlops approach can start to pay dividends. The same advantage will of course be present again in other measurements also covering both the low-and high-p T domains, most obviously, the vector boson p T spectrum, as seen already in Figs. 14 and 15.
Having said that, let us re-emphasize the point made earlier in regards to the jet p T spectra in figure 16, namely that we foresee the main advantage in the Menlops method being in the modeling of W and Z production as background processes, where the significance of the various jet multiplicity bins is not generally clear. Hence, again, the real practical value of these plots is to demonstrate that the Menlops approach always embodies the best of the Nlops simulations and to show that the sensitivity of merged Menlops pre-dictions to the merging scale is marginalized (as illustrated by the plots on the right-hand side of the figures).  MC/data Figure 19: Here, on the left, we confront Powheg-W (red), Powheg-Wj (blue) and Menlops 35 (green crosses) predictions with Atlas collaboration measurements [54] of jet multiplicity and jet p T spectra in the muon decay channel. On the right-hand side we show the results from the corresponding Menlops 0 (magenta) and Menlops ∞ (cyan dashes) improved versions of these simulations, along with that obtained by merging their output ( Menlops 35 ).

Conclusions
Firstly in this article, we set out to investigate whether Powheg simulations of simple processes could be promoted to Menlops ones without approximation [24]. This was described and demonstrated in Sect. 2.2, where only a fairly modest effort was required; effectively recycling much of the radiation generating machinery of Powheg vector boson plus jet simulations so as to generate next-to-next-to-leading order emissions according to double emission matrix elements, instead of the parton shower approximation. The method was shown to preserve the NLO accuracy of inclusive observables and the logarithmic accuracy of those sensitive to Sudakov effects. In all cases the Menlops enhancement, labeled Menlops ∞ , provided a much improved description of Tevatron and LHC data with respect to that of the underlying Powheg-V program, specifically, for observables sensitive to the emission of more than one jet. Given the relative ease of the implementation and the findings of the validation exercise, we believe that the same technique can be readily and successfully applied to other processes, in particular, that it would be straightforward to include the vector boson plus three partons matrix elements in the same way. Our second goal was to seek improvement and similarly extend the reach of vector boson plus jet simulations (the Powheg-Vj programs). Here we have again worked in such a way as to minimize the interference with existing programming, directly utilizing the radiation generating apparatus of the Powheg-V programs to induce Sudakov suppression effects when the vector boson has low transverse momentum. The method was shown to preserve the NLO accuracy of vector boson plus jet observables and the all-orders resummation of Sudakov logarithms, while at the same time giving predictions for fully inclusive observables in remarkably good agreement with genuine NLO ones (better than one has right to expect). The Menlops 0 enhancements lead to a much improved description of Tevatron and LHC data with respect to that of the underlying Powheg-Vj program, for the case where the definition of the observables had some overlap with the Sudakov region, e.g. vector boson p T and rapidity distributions, as well as the 0-jet cross section (for which the latter program exhibits unphysical behavior). We also remark that the benefits here were not restricted to the program's output but also extended to its run-time, which was found to be greatly reduced due to it evaluating less complex NLO cross sections in the Sudakov region.
Following these modifications we have considered the combination of the resulting event samples in such a way as to yield a merged sample improving on both of the components: inheriting both the NLO accuracy of the enhanced vector boson production simulation Menlops ∞ and that of the corresponding vector boson plus jet simulation Menlops 0 . We have discussed in Sect. 2.4 how the merging scale must be bounded from below, to enforce that the fraction of Menlops 0 events in the final sample not exceed O (α S ) of the total, to maintain the NLO accuracy of inclusive observables.
As noted in Sect. 2.4, if optimal renormalization and factorization scale choices are used, the latter constraint should effectively render the resummation in the vector boson plus jet simulations null in the context of merging samples, with the Menlops 0 sub-sample being open to replacement by one from the normal Powheg-Vj program. We therefore use the Menlops improved vector boson plus jet simulation in merging out of theoretical correctness. On the contrary, the same is not true for the Menlops ∞ simulation, which if replaced by its unenhanced Powheg-V program would lead to two-and three-jet events being underestimated by the parton-shower soft collinear approximation (exactly how much depends on the value of the merging scale). In any case, by merging Menlops samples rather than Nlops ones we minimize the merging scale dependence: even for the most extreme variations we seldom find the Menlops 0 , Menlops ∞ and merged Menlops predictions out of agreement by more than 10% (typically this occurs in regions were the predictions are all no better than leading order anyway). The same cannot be said of the unimproved Powheg results. With this last point duly noted, it seems likely that a true matching of the Nlops simulations, without any merging scale dependence, would offer little in the way of practical improvements beyond this pragmatic approach.
Good while this may sound, in regards to actually describing W and Z (plus jet) production as well defined signal processes, none of these enhancements seem to offer an improved description over what could be done before, by taking care to use the right tool for the right job (in the right way). Inclusive vector boson production measurements are described as well by Powheg-V programs as by any Menlops ones and the same is true in regards to the Powheg-Vj programs and vector boson plus jet measurements. In real terms, the main potential improvement with regard to testing QCD is one of convenience, through being able to perform inclusive and exclusive analyses with the one simulation / event sample. On the other hand, when viewed as less well defined background processes, one does not in general make any clear distinction between vector boson and vector boson plus jet production. It is in these applications that we believe the approach taken here can be expected to offer real gains. Also, it is from this perspective that the trouble taken to minimize the merging scale dependence, the validation exercises and the favorable comparisons with real data, take on their true significance.
While the latter simple approach (Sect. 2.4) is not an exact theoretical solution to the Nlops-Nlops matching problem, based on the comparison with data in Sect. 3, it would appear to be an effective, convenient and extensible one. Comparisons with a truly scale-free method will be interesting when one is realized in the near future.
In conclusion we wish to recommend the Menlops procedures documented here, in particular the merging of Menlops-improved Powheg simulation output, as an efficient and clear means to get the most from existing and future Nlops simulations.

Acknowledgments
KH is very grateful to Paolo Nason and Torbjörn Sjöstrand for patiently answering many questions regarding this work, and particularly to Bryan Webber for reading and suggesting improvements to all versions of the manuscript. KH thanks Italian Grid Infrastructure, TheoPhys VO and INFN Milano-Bicocca computing support, without whose help and resources this work would not have been possible. The work of SA is supported by the Helmholtz Gemeinschaft under contract VH-HA-101 (Alliance Physics at the Terascale) and by Deutsche Forschungsgemeinschaft in Sonderforschungsbereich/Transregio 9. ER and SA also acknowledge financial support from the LHCPhenoNet network under the Grant Agreement PITN-GA-2010-264564.