NNLOPS simulation of Higgs boson production

We detail a simulation of Higgs boson production via gluon fusion, accurate at next-to-next-to-leading order in the strong coupling, including matching to a parton shower, yielding a fully exclusive, hadron-level description of the final-state. The approach relies on the POWHEG method for merging the NLO Higgs plus jet cross-section with the parton shower, and on the MINLO method to simultaneously achieve NLO accuracy for inclusive Higgs boson production. The NNLO accuracy is reached by a reweighting procedure making use of the HNNLO program.


Introduction
Presently, the discovery of a new spin-zero particle in the search for the Standard Model (SM) Higgs boson at the Large Hadron Collider (LHC), by the ATLAS and CMS collaborations [1,2], has focused the physics agenda on studying its properties, to search for possible departures from the SM predictions. The accurate measurement of the couplings and quantum numbers of the new particle forms the basis of this endeavour. Thus far, with present data, the newly discovered boson shows no significant deviations from the SM expectations.
One of the fundamental limitations in our ability to study the nature of the new particle is the precision afforded by our theoretical predictions in describing its production and decay rates, and the associated experimental acceptances. Now, following the initial discovery, with the mass of the would-be Higgs boson known, it is reasonable to expect that the associated data will grow at a considerable rate when the LHC restarts in 2014/15. It is foreseeable then, that as the statistical errors on the experimental measurements diminish, the accuracy of the theoretical predictions will become an even more pressing issue, as work continues towards constraining deviations from the apparent SM behaviour.
The anticipation of new and exciting discoveries at the LHC has catalysed significant progress in addressing concerns regarding the precision of theoretical predictions (for a recent review see e.g. Ref. [3] and references therein). Notably, the last ten years have witnessed remarkable developments in precision fixed order QCD calculations for the main (gluon fusion) Higgs boson production channel, beginning with next-to-next-to-leading order (NNLO) predictions for the total inclusive cross section [4][5][6], followed by analogous computations of differential quantities [7,8]. Most recently, landmark calculations have been carried out in deriving the full or partial NNLO QCD corrections to a number of 2 → 2 partonic scattering processes [9][10][11][12], including gg → H + jet [13].
While the emergence of fully differential fixed order NNLO calculations represents an important theoretical breakthrough, the practical value of which is underlined by their use in the Higgs boson discovery analyses [1,2], their description of the final state is comprised of at most two QCD partons. Multiple parton emission, resummation effects, hadronisation and the underlying event are not accounted for. All of the latter class of contributions represent corrections formally beyond NNLO for inclusive quantities. They are however very relevant for the exclusive description of the final state. Furthermore, depending on the nature of the observables under consideration, they may sum to give corrections as significant as the NNLO ones, or supersede them altogether, in kinematic regions where fixed order perturbation theory becomes unreliable.
In parallel with advances in fixed order computations, significant progress has been made on the description of exclusive final states, most notably through the development of methods for the consistent inclusion of next-to-leading order (NLO) matrix elements in parton shower Monte Carlo event generators (Nlops) [14,15]. In less than ten years since their inception, these matching schemes have evolved to the point where the construction of Nlops simulations, for even complex multi-leg processes, proceeds with a high level of automation [16][17][18][19], with little or nothing required in the form of Nlops expertise. Like the fixed order NNLO computations, the value of Nlops simulations is well established and they too have been used in the Higgs boson discovery analyses.
The rapid development in Nlops techniques in the last three years has led to the creation of a number of public codes for most processes of interest at the LHC, among which those of the form pp → X +n jets, with X being e.g. a W/Z/H boson [19][20][21][22][23][24]. In view of the positive experience with matrix-elements and parton shower merging schemes (Meps) [25][26][27][28], several efforts have appeared in the literature towards the developments of methods that allow the merging of NLO samples with different associated jet multiplicities [29][30][31][32][33][34][35][36][37][38][39][40]. We now focus upon the example of inclusive Higgs and Higgs plus one jet generator (henceforth H and Hj). Ideally, the merging procedure should yield a generator that is NLO accurate for both fully inclusive observables, and for observables requiring the presence of one associated jet. Recall that the H generator is NLO accurate for inclusive observables but only LO accurate in the description of the associated jet, while the Hj generator is NLO accurate in the description of the associated jet, but it cannot be used to compute fully inclusive observables. It is clear that the construction of a 'merged' generator of this kind is a first step in the development of a NNLO accurate generator matched to a parton shower. In fact, the latter is better than the former only for inclusive quantities, where the former is only NLO accurate.
Among the new multi-jet merging schemes the Minlo approach, as it is set out in [39], is unique in that it does not resort to the introduction of an unphysical merging scale cut to partition the phase space into a 0-jet region, to be populated by the H Nlops simulation, and a ≥ 1-jet region, to be filled by the Hj one. On the contrary, the Minlo recipe acts to extend the reach of the underlying Hj NLO computation so as to return NLO accurate predictions for observables inclusive with respect to all QCD radiation, without the introduction of any additional unphysical parameters. As a consequence, the Hj-Minlo event generator provides all of the same fixed order accuracy as a fully differential NNLO computation, except fully inclusive observables, for which it is only NLO accurate. In ref. [39] it was argued that, by a simple reweighting procedure, full NNLO accuracy can actually be achieved with such a generator. In the present work we implement this reweighting, yielding a first true Nnlops simulation for Higgs boson production. The fact that the standalone Hj-Minlo simulation already achieves NLO accuracy for inclusive Higgs boson production observables is at the heart of this development. Were it not for this feature, the magnitude of the weights required to correct the inclusive distributions to NNLO accuracy would be such that the NLO precision for Hj observables is then lost. The remainder of this paper is structured as follows. In section 2 we present the theoretical framework underlying the Nnlops method, starting with the basic formulation in section 2.1, followed by simple refinements in section 2.2. Section 3 describes our prescription for the determination of theoretical uncertainties in the new method. In section 4 we present a selection of results obtained having implemented our Nnlops method, with the aim of probing and validating the procedure. Lastly, in section 5, we present our conclusions and comment on further developments.

Theoretical framework
In ref. [39] it was proven that the Hj-Minlo computation is NLO accurate for Higgs plus one jet (Hj) and Higgs (H) inclusive observables. As noted in the introduction and outlined in ref. [39], this property is crucial for the promotion of the Hj-Minlo simulation from Nlops to Nnlops. In this section we describe in detail the procedure for reaching NNLO accuracy. For details regarding how the Hj-Minlo simulation first attains NLO accuracy for inclusive Higgs production observables we refer the reader back to ref. [39].

Method
Let us abbreviate by dσ MINLO the cross section obtained from the Hj-Minlo event generator, fully differential in the final state phase space, Φ, at the level of the hardest emission events, i.e. prior to showering. On integration, this distribution reproduces the leading order, O α 2 S , and next-to-leading order contributions, O α 3 S , in the perturbative expansion of the Higgs boson rapidity distribution. In addition, unenhanced spurious terms entering at O α 4 S and higher are present, since the fully differential cross section includes all order contributions in the Sudakov form factors, and contributions of order O α 4 S to the Higgs plus jet cross section. By analogy we denote the conventional fixed order, next-to-next-toleading order differential distribution by dσ NNLO . Since the rapidity distributions derived from dσ MINLO and dσ NNLO are formally identical up to and including O α 3 S terms, it follows that their ratio is equal to one up to O α 2 S terms: where the c i are simply constant O (1) coefficients.
While it is obvious that by this reweighting the inclusive rapidity distribution acquires NNLO accuracy, the crucial point here is that the NLO (i.e. O(α 4 S )) accuracy of the cross section in the presence of jets (that starts at order α 3 S ) is maintained, since the reweighting factor combined with this cross section yields spurious terms of order α 5 S and higher. We stress again that, were it not for the special property that the Hj-Minlo generator reproduces the conventional fixed order result up to and including NLO terms, W (y) would yield relative corrections of O(α S ), thus spoiling the NLO accuracy of Higgs plus one jet distributions.
We shall now demonstrate that the Hj-Minlo generator reweighted with the procedure outlined above achieves O(α 4 S ) accuracy for all observables. To begin with we must prove the following theorem: A parton level Higgs boson production generator that is accurate at O(α 4 S ) for all IR safe observables that vanish with the transverse momenta of all light partons, and that also reaches O(α 4 S ) accuracy for the inclusive Higgs rapidity distribution, achieves the same level of precision for all IR safe observables, i.e. it is fully NNLO accurate.
To this end, we consider a generic observable F that is an infrared safe function of the final state kinematics. Its value will be given by with a sum over final state multiplicities being implicit in the phase space integral. Infrared safety ensures that F has a smooth limit when the transverse momenta of the light partons vanish. Such a limit may only depend upon the Higgs boson's rapidity, y, since it is the only observable left when no other partons are resolved. We generically denote such a limit by F y . The value of F can be considered as the sum of two terms: F − F y + F y . Since, F − F y tends to zero with the transverse momenta of all the light partons, by hypothesis its value is given with O(α 4 S ) accuracy by the parton level generator. On the other hand, which is also exact at the O(α 4 S ) level by hypothesis. Thus, F = F − F y + F y is accurate at the O(α 4 S ) level, proving our theorem.
This theorem is easily generalized to arbitrary processes. It is enough to replace F y in the above proof with F Φ B , where Φ B are infrared safe quantities parametrising the associated Born phase space. The Hj-Minlo parton level generator is, by itself, one that fulfills the first condition required by our theorem; predicting IR safe observables that vanish when the transverse momentum of the light partons vanishes with O(α 4 S ) accuracy. In fact, this would be the case even without the Minlo improvement. The second hypothesised statement entering the theorem, regarding NNLO accuracy of the Higgs boson's rapidity spectrum, is realised by augmenting the Hj-Minlo generator by the reweighting procedure described above. The proof of O(α 4 S ) accuracy for these observables thus corresponds to the general proof of NLO accuracy of the Powheg procedure, given in refs. [41,42].
Observe, also, that for observables of the type F − F y , adding the full shower development does not alter the O(α 4 S ) accuracy of the algorithm, for the same reasons as in the case of the regular Powheg method. The only remaining worry one can then have concerns the possibility that the inclusive Higgs boson rapidity distribution is modified by the parton shower evolution at the level of O(α 4 S ) terms. However, our algorithm already controls the two hardest emissions with the required α 4 S accuracy. A further emission from the shower is thus bound to lead to corrections of higher order in α S . 1 We notice that the Nnlops event generator described here is NNLO accurate in the same sense in which the current MC@NLO or Powheg type generators are NLO, i.e. infrared safe observables are NLO or NNLO accurate. The degree of logarithmic accuracy, leading, next-to-leading or next-to-next-to-leading, is contingent upon the way in which the associated Sudakov form factors are implemented.
It should be clear from our discussion that reweighting can be performed at the partonic level either before or after the shower. It cannot be carried out at the level of the generation of the underlying Born configuration of the Hj-Minlo generator, since the Higgs rapidity changes after radiation. 2

Variant schemes
One can readily construct simple variants of the method discussed in sect. 2.1. In particular, rather than performing a global reweighting of all Hj-Minlo events using the fully inclusive Higgs boson rapidity distribution, one can instead consider splitting the cross section according to Recall that Powheg simulations limit the transverse momenta of branchings in the subsequent parton shower simulation to be less than the hardest emission generated with respect to the underlying Born (passed through the so called scalup variable), which here is the second emitted parton. 2 The Powheg Box implementation guarantees that the rapidity of the system comprising the Higgs plus the hardest parton remains the same after radiation, but not the rapidity of the Higgs itself.
where p T here represents some overall measure of the hardness of radiation in the event, with h a monotonic profile function such that lim p T →0 h(p T ) = 1, lim p T m H h(p T ) = 0, and simply reweight the dσ A component rather than the full cross section. A suitable form for the profile function is where β and γ are constant parameters. We reweight the Hj-Minlo events with the factor Multiplying the above equation by dσ MINLO δ(y − y(Φ)), using equations (2.7) and (2.8), and integrating over the full phase space we obtain identically Eq. (2.11) differs from the fixed order NNLO cross-section only by terms of order α 5 S . In this work we will adopt a further modification of the reweighting factor, that has the advantage of yielding a Higgs rapidity distribution that coincides exactly with the NNLO result: which leads precisely to dσ dy The purpose of the h profile function is quite similar to what is done sometimes in Powheg, when the real emission cross section is separated into a singular and a finite part [15,42,43]. The only difference, in this case, is that, rather than an inclusive LO-to-NLO correction, here we include an NLO-to-NNLO correction. This correction, in the fixed order calculation, is concentrated in the region of zero transverse momenta of the radiated partons, while in a resummed calculation like Nlops or Nnlops, this is no longer the case, the zero transverse momentum region being suppressed by a Sudakov form factor. Thus, the correction must be spread over a range of non-zero transverse momentum.
To facilitate a more intuitive understanding, we point out that in the limit γ → ∞, . Thus, taking for example the leading jet transverse momentum to define the argument p T of the h function, we see that dσ A and dσ B in eqs. (2.6-2.8) are nothing more than the usual 0-and ≥ 1-jet cross sections. Hence, in this limit, the reweighting procedure merely amounts to rescaling the weights of the 0-jet events by the ratio of their respective NNLO-to-NLO cross sections, albeit differentially in the Higgs boson's rapidity. More generally, moving away from γ = ∞ towards finite values, the effect on h (p T ) is to smear the step in the θ (β m H − p T ) function. The h profile function is therefore most easily thought of as a smeared step function. In this work we have only performed studies with γ = 2. Nevertheless we consider the results we find with this parameter choice to be wholly satisfactory.
Turning to the β parameter in the profile function, one sees that increasing β increases the dσ A component of the cross section relative to dσ B . In fact, in the limit β → ∞, dσ B = 0 (eq. 2.8) and one recovers the simple global event reweighting of eqs. (2.1-2.3). If we choose β ≈ 1, the NNLO correction factor in W (y, p T ) is applied in a region where radiation is not much harder than m H . It is also clear that we cannot take β 1; if we do so, the NNLO correction will be concentrated in a small region of the radiative phase space, p T β m H m H . It thus becomes a delta function as β → 0, spoiling the accuracy of the resummation. Effectively, β must be of order one to avoid such a pathology.
The β parameter shares some features with the ratio of the resummation scale to the heavy boson mass in matched NNLO analytic resummation calculations (see e.g. [44]). In conventional resummation calculations the resummation scale affects both the logarithms which are resummed and also how far both the hard NLO and NNLO virtual corrections are to be distributed along the transverse momentum spectrum. Here β plays an analogous role but it only affects the distribution of the hard NNLO virtual corrections; the argument of the logarithms being resummed and the distribution of NLO virtual corrections is unaffected by it. By the same analogy one should consider the 'sensible' range in which to vary β as being limited to the same range in which the resummation scale (divided by m H ) is varied in conventional analytic resummation calculations.
Before ending our discussion on the β parameter, we wish to emphasise that while its precise value is a source of systematic uncertainty, it is fundamentally different to the merging scales encountered in all other recent attempts to merge Nlops simulations for multi-jet processes; even if a NNLO reweighting of such simulations were to be admissible, e.g. by somehow having the equivalent NNLO reweighting function of the form 1 + O α 2 S , the merging scale would remain as an additional source of systematic uncertainty. Moreover, while the dependence on β is formally O(α 5 S ) or even zero in the case of inclusive quantities, in all other recent Nlops merging attempts, the dependence on the respective merging scales is O(α 4 S ) or worse. Lastly, stepping back from the technicalities of the profile function h, we point out that in these variant reweighting schemes the extension of the proof of NNLO accuracy to variables other than the inclusive Higgs boson rapidity spectrum follows that given earlier with only trivial adjustments.

Estimating uncertainties
We now examine the source of uncertainties in our NNLO generator, and set up a method for determining its theoretical errors.
The uncertainties in the Hj-Minlo generator are explored according to the prescription we gave in ref. [31]. There we have considered a 7-point scale variation, where all renormalization scales appearing in the Minlo procedure are multiplied by a scale factor K R , and the factorization scale is multiplied by a factor K F , with (K R , K F ) = (0.5, 0.5), (1, 0.5), (0.5, 1), (1, 1), (2, 1), (1,2), (2,2) . (3.1) We will consider the variation in our results induced by the above procedure. We compute dσ NNLO /dy using the Hnnlo program of ref. [8,45]. The theoretical uncertainty in this calculation can be estimated by performing a factorization and renormalization scale variation in the usual way. The choice of the central scale has been subject of some debate. The value m H has been used for a long time, but recently the value m H /2 seems to be preferred, on the grounds that it yields smaller NLO and NNLO corrections. We will thus take m H /2 as the central scale choice for the NNLO calculation. Also in the case of the NNLO calculation, we consider the variations of the renormalization and factorization scales of eq. (3.1), this time applied at the central value m H /2. This yields 49 variations in the Nnlops result. On the other hand, we found that by limiting ourself to K R = K F in the NNLO result no appreciable reduction of the scale variation envelope is observed. We therefore restrict ourselves to this case, thus ending up with 21 scale variation points. On top of this, we have freedom in the choice of the p T variable and of the constant β of eq. (2.9). In the implementation we have chosen to define the argument of the profile function, p T , as the transverse momentum of the hardest jet, computed according to the inclusive k T -algorithm with R = 0.7. This choice has the advantage that the region p T → 0 is approached only when all radiated partons have vanishing transverse momenta. We have verified that other choices, such as the transverse momentum of the Higgs, do not lead to significant differences. We assume as default β = 0.5, and consider variations between 0.5 and ∞.
In order to perform our study, we have generated a single sample of Hj-Minlo events. The seven scale variation combinations have been obtained by using the reweighting feature of the Powheg Box. 3 The integrals dσ MiNLO A/B /dy, needed in eq. (2.12), were performed and tabulated for each scale variation combination using the Hj-Minlo generated sample. Similarly, dσ HNNLO /dy was tabulated for each of the three scale variation points. The analysis is then performed by generating the Minlo event with given values of (K R , K F ), and multiplying its weight with the factor The central value is obtained by setting (K R , K F ) and (K R , K F ) equal to one, while to obtain the uncertainty band we apply this formula for all the seven (K R , K F ) and three The conservative rationale/ansatz adopted here, in estimating errors by varying the scales in the NNLO and NLO inputs in a fully independent way, is essentially that we regard the uncertainties in the normalizations of distributions, e.g. the transverse momentum spectrum of the Higgs boson, as being independent of the respective uncertainties in the shapes -at least in the region covered by the profile function, h(p T ), i.e. that which includes the low p T domain. The former are determined by the Hnnlo program, while the latter are due to the Hj-Minlo input. Outside of the low p T region, in the part corresponding to the 1−h(p T ) term in eq. (3.2), the uncertainty is given by the standard Hj-Minlo computation (which there corresponds to that of conventional NLO with µ R = µ F = p T for the central scale choice). Thus, in the low p T region the absolute uncertainty at a given point in a distribution which is not inclusive, e.g. a p T spectrum, is essentially given by the uncertainty in the shape at that point times the uncertainty in the normalization. Readers familiar with such matters may recognise the approach taken here as being analogous to the so-called efficiency method [46], used for estimating errors on cross sections in the presence of cuts; where uncertainties on the theoretical predictions for the total cross section and the associated efficiencies are assummed to be uncorrelated.

Phenomenological analysis
In this section we present our phenomenological results. We consider the production of a 125.5 GeV Higgs boson at the 8 TeV LHC. Throughout, we use the MSTW8NNLO [47] parton distribution functions for all our results, including the Minlo ones. We remark that also the CT [48], and NNPDF [49] collaborations have produced NNLO fits and could have been used in this context. However, here we are not interested in PDF comparisons, and will stick to a single set for simplicity. The Nnlops events include parton showering, as determined by Pythia 6 [50], with Perugia 0 tune [51] (PYTUNE(320)). We switched off hadronization and multi-parton interactions, in order to carry out a more sensible comparison with other parton level generators.
Throughout this paper, to define jets we used the anti-k T algorithm [52] as implemented in FastJet [53,54].
As the reader may have noticed, there is a slight tension between the choice of scale in the Hnnlo and Nnlops calculations, for observables which are not fully inclusive with respect to all QCD radiation. In particular, in Hnnlo, due to the nature of the calculation, it is not possible to use a dynamical scale, as is done in Hj-Minlo. 4 In this work we elect to use m H /2 as the central scale in Hnnlo, as input to our reweighting procedure, and in comparing to its predictions for inclusive observables: since this setting is favoured by the community of Higgs NNLO experts in determining the total inclusive cross section. On the other hand, for jet cross sections at moderate and large transverse momenta, that scale is generally considered to be too low. Thus, when comparing to jet observables, we have instead used µ F = µ R = m H as the central scale choice in the Hnnlo predictions (i.e. still maintaining µ F = µ R = 1 2 m H for the Hnnlo input to the Nnlops reweighting procedure). In all cases, to obtain the Hnnlo uncertainty band we perform a standard 7-point scale variation around the central value.

Higgs boson rapidity spectrum
We begin by showing in fig. 1 the fully inclusive Higgs rapidity distribution. This, by construction, should be identical in the Hnnlo and Nnlops calculations.
On the left, in the red shaded area, one can see the scale uncertainty band predicted by the Nnlops simulation, with the conventional fixed order Hnnlo result superimposed as green points. The lower panel shows the ratio with respect to the Nnlops prediction obtained with its central scale choice. On the right we have made the same plots as on the left but with the Hnnlo predictions replacing those of the Nnlops and vice versa; the scale uncertainty bands are formed as described in Sec. 3. In the following we will compare the Nnlops to other results with plots of the same kind. As expected, for this observable the two calculations are in full agreement, both for their central values and scale uncertainty envelopes; the latter being approximately ±10% in size.

Higgs boson transverse momentum
Here, to begin with, we wish to discuss the evolution of the Nnlops program's prediction, at each of the main stages of the simulation process, as part of its validation and in order to provide relevant background, before comparing it to state-of-the-art resummed calculations. In figure 2 we show how the Higgs boson transverse momentum spectrum is affected at the various phases of the event generation process in the underlying Hj-Minlo simulation (as described in sects. accuracy if the small transverse momentum region is excluded.
The first most obvious feature is the difference between the various Hj-Minlo predictions and those of the conventional fixed order program Hnnlo in the low p T region, with the latter exhibiting unphysical divergent behaviour, and the former displaying, instead, the anticipated, physical, Sudakov peak. In the high p T tail region all of the predictions are in  good agreement; we have verified that they are within each other's scale uncertainties. Nevertheless, the lower renormalization and factorization scales at high transverse momentum in the Hnnlo result leads to a less steeply falling cross-section in that region, as evident from the lower panel where the ratio to the Hj-Minlo showered result is displayed. We remind the reader that for p T > m H the Hj-Minlo program reverts to µ R = µ F = p T with the Minlo Sudakov form factor also being equal to one in this region.
Looking among the three Hj-Minlo predictions, at high and moderate p H T the three distributions are almost indistinguishable from one another, agreeing at the level of ∼2-3%, very much within the Hj-Minlo scale uncertainty envelope -a nearly flat band, with a width of approximately ±20% in the ratio subplot (see fig. 3 of ref. [39]). This remarkable agreement is easily understood as being due to the inclusivity of the observable considered here (for more details see Appendix A).
Differences among the Minlo predictions only become apparent in the low transverse momentum region below 50 GeV. We attribute these more prominent deviations as being due to the amplification of NNLO sized differences (see Appendix A), between the hardest emission cross section (blue) and enhanced fixed order (red), by large logarithms of p H T /m H . We also note the expected vanishing of the Hj-Minlo fixed order prediction (red) at low p H T , by virtue of the fact that the NLO computation includes a Sudakov form factor which is a function of the Higgs boson's transverse momentum. The Hj-Minlo hardest emission cross section and parton shower level predictions smear out this region, by transforming underlying Born configurations with low p H T into configurations where the Higgs boson has zero transverse momentum, through generating the hardest emission and also subsequent, shower emissions in the latter case.
The differences in the region below 50 GeV can also be understood from a different point of view, by noting that for the case that the observable O is the total inclusive cross section, it is clear that the Hj-Minlo enhanced fixed order prediction and its derivatives, namely, the hardest emission cross section and subsequent parton showered predictions, all agree identically; thus, a relatively low contribution in the vicinity of the p H T ∼ 0 GeV in the case of the enhanced fixed order prediction (red) must be compensated by it having a relatively high contribution elsewhere, in this case the region ∼ 15 − 50 GeV.
In summary, we have seen that the Hj-Minlo predictions at the NLO, Les Houches, and showered level are in close agreement, the largest discrepancy (near 10%) occurring in the Sudakov region, where effects beyond NLO are numerically significant, for reasons which are well understood.
The final ingredient to reach the NNLO accuracy is the inclusion of the reweighting procedure discussed in sections 2.1 and 2.2. In figure 3 we display the effect of the inclusion  of the NNLO reweighting with respect to the Hj-Minlo result, for β = ∞ and β = 1 2 . In the β = ∞ case, the NNLO reweighting can be well modeled by an overall K -factor, that does not modify the shape of the transverse momentum distribution at all. This is easily understood, since in practice the reweighting factor has a fairly mild dependence upon the rapidity. By introducing a finite β we do instead alter the shape of the transverse momentum distribution, since, in this case, the K -factor is only applied to the lower portion of the p T spectrum. We observe that the NNLO correction factor is quite large, around 1.5, in the small transverse momentum region, where the bulk of the cross section lies. We remind the reader that in carrying out the reweighting here, we have set µ F = µ R = 1 2 m H in the Hnnlo program and used the default Hj-Minlo settings (which correspond well, in the case of inclusive quantities, to conventional NLO predictions with µ F = µ R = m H ). Had we chosen µ F = µ R = m H in determining the Hnnlo input to the reweighting procedure, the correction factor would be near 1.3.   In figures 4 and 5, we compare the Nnlops (see eq. (3.2)) with the HqT [55, 56] result for two choices of the β parameter in the profile function. The uncertainty band is the envelope of the 21-point scale variation illustrated in section 3. We used the 'switched' output of HqT, forming the related uncertainty band from the envelope of the seven results obtained by independent variations of µ R and µ F , by a factor of two, symmetrically, about µ R = µ F = 1 2 m H , while keeping the resummation scale always fixed to 1 2 m H . Pleasingly, we see that the Nnlops and HqT results are almost completely contained within each other's uncertainty band in the region of moderate transverse momenta. We have verified that at high transverse momentum the HqT prediction agrees identically with that of Hnnlo, since the 'switched' output in the former uses the fixed order result in this region. It follows that here we see the HqT spectrum falling less rapidly than that of the Nnlops simulation at large p H T . As was seen in fig. 3 and remarked upon in the related discussion, in the case of β = ∞, the Nnlops result is very well approximated by that of Hj-Minlo multiplied by a uniform NNLO-to-NLO K -factor of 1.5, leaving the slope of the distribution unchanged. On the other hand, for β = 1 2 ( fig. 5) the K -factor enhancement is predominantly concentrated in the region p T 1 2 m H , yielding a modest but marked change in the shape of the Nnlops distribution. In this case, the discrepancy observed earlier, in fig. 2, between the fixed order calculation (Hnnlo) and the Hj-Minlo result at high transverse momentum remains unaltered, both in terms of the shape and normalization of the spectrum in that region. Once the partial NNLO calculation of Higgs plus one jet [13] will be complete, a comparison to it will certainly provide further insight.
Notice that the scale uncertainty band in the Nnlops calculation, for high Higgs p T , is larger in the β = 1/2 than in the β = ∞ case. This is easily understood. By reweighting we reduce the scale dependence of the Hj-Minlo result for the inclusive cross section, so that at the end we have a scale variation that is appropriate to the NNLO calculation. In other words, in the β = ∞ case, reweighting also partially compensates the scale variation in the large transverse momentum tail. This leads to an underestimate of the theoretical error in this region, since the high transverse momentum tail, that is only computed with NLO accuracy, gets a scale variation of relative NNLO order. This is not the case for β = 1/2, for which the associated scale uncertainty is characteristic of the NLO accuracy at high p T .
We remark that the choice of the β parameter is related, to some extent, to the choice of the resummation scale in conventional resummed calculations. In the latter, the hard function, i.e. the factor 1 + C 1 α S + . . . in front of the resummed expression, 5 which embodies hard virtual corrections to the leading order process, is present in the resummed component in the region below the resummation scale p T < Q. Whereas, in conventional resummation calculations Q governs the p T range for the full NNLO and NLO hard virtual corrections, as well as affecting the argument of the logarithms being resummed, in the Nnlops case the scale βm H only determines the extent of the second order hard virtual corrections ∼ C 2 α 2 S . Motivated by this correspondence, in this work we favour values of β not larger than one, the value 1 2 corresponding to the preferred choice of resummation scale in HqT (Q = 1 2 m H ). In the following we will thus stick to the β = 1 2 choice as our default. We remind the reader, however, that in the β = ∞ case, the modification of the spectrum at high transverse momentum is an effect of order α 5 S , i.e. beyond our intended accuracy. In finishing this discussion we note that in the β = ∞ case the level of agreement between HqT and the Nnlops is better than between HqT and the Nlops Powheg Higgs production program in fig. 22 of ref. [58], while for the optimal setting β = 1 2 the agreement between HqT and the Nnlops is quite satisfactory.

Leading jet transverse momentum
In this subsection we turn from the Higgs boson p T spectrum to that of the leading jet, p j 1 T . In figure 6 we show results for the leading jet transverse momentum distribution, all accurate up to and including O α 4 S contributions: the conventional NLO QCD prediction from Hnnlo with µ R = µ F = m H (black), the Hj-Minlo enhanced fixed order prediction as described in sects. 2 and 3 of ref. [39] (red), the Hj-Minlo hardest emission cross section (Hj-Minlo LH, blue), and the ensuing parton shower level prediction from Pythia (green). In the case of R = 0.7, in the left-hand plot, we see qualitatively the same pattern of results as for the Higgs boson transverse momentum spectrum. At least to leading order in perturbation theory, the jet and the Higgs boson recoil against each other with equal and opposite momenta in the transverse plane. At higher orders this picture is modified, with multiple parton emissions leaking energy outside the jet depleting its transverse momentum. We note that, for the cases we consider, all R dependence of the cross section originates from real emission contributions to the Hj process; more precisely, for the cross section to have any R dependence there must be at least two partons in the final-state. Since the radius R = 0.7 is quite large, one expects any radiation leakage to be small, indeed, in the limit R → ∞ the jet algorithm clusters all final state partons together into a single jet, which, by momentum conservation, must exactly recoil against the Higgs boson. This being the case we offer the same explanation for the pattern of results shown in the left of fig. 6 as in fig. 2.
For R = 0.7 jets we note just one small difference, namely, that in the case of the Higgs boson transverse momentum spectrum, in the vicinity of p T = 0 GeV, the enhanced fixed order prediction (red) was greatly suppressed relative to the hardest emission cross section and the showered result (blue and green) -a deficit which was compensated by a slight excess in the p T range 15-50 GeV -while here it is much more compatible with them, in fact, if anything there is a slight excess of the former over the latter. This feature in the Higgs boson transverse momentum spectrum was understood as being due to the fact that the Sudakov form factor in the Minlo formulation of ref. [39] is directly a function of the boson's p T , giving (formally) an infinite suppression at p T = 0 GeV. In the present case, instead, a small p T of the hardest jet does not imply the vanishing of the Higgs p T even at the NLO level, and such suppression is not present, leading to better agreement in the p T < 50 GeV region at each stage of the event generation process.
For radius R = 0.4, the leading jet transverse momentum spectrum is shown on the right of fig. 6. As in the case of R = 0.7 at high transverse momentum the Hnnlo result is less steeply falling than the Hj-Minlo ones, due to the different scale choice. We notice that for R = 0.4 the Hj-Minlo+Pythia result is about 10% smaller than the Hj-Minlo-LH one in the mid-to-high p T range. In fact, as pointed to in the preceding paragraph, multiple emission parton shower effects cause the progenitor partons in the Les Houches events to lose energy and hence reduce the cross sections for any jet which had previously been associated to them. This leakage of radiation outside the jet is clearly amplified for smaller jet radii.
In fig. 7 we have plotted Nnlops and NNLL+NNLO JetVHeto [59] predictions for the jet veto efficiency, ε (p T,veto ), defined as the cross section for Higgs boson production events containing no jets with transverse momentum greater than p T,veto , divided by the respective total inclusive cross section, The jets considered here are formed according to the anti-k T jet algorithm [52], for a variety of different jet radii; from top to bottom in fig. 7 we have, pairwise, R =0.4, 0.5, 1.0. In the left-hand column, in the red shaded area, we show the scale uncertainty band predicted by the Nnlops simulation, with the central NNLL+NNLO resummed prediction of JetVHeto superimposed in green (matching scheme-(a), µ R = µ F = µ Q = 1 2 m H , µ Q being the resummation scale). The lower panel shows the ratio with respect to the Nnlops prediction obtained with its central scale choice. On the right we have made the same plots as on the left but with the JetVHeto predictions replacing those of the Nnlops and vice versa.
The uncertainty band in the JetVHeto results is the envelope of a seven point variation of µ R and µ F by a factor of two. This is in contrast to the band associated with it in the predictions of ref. [59], where additionally resummation scale and matching scheme variations were included in the envelope. Thus the JetVHeto error band here is considerably smaller than that shown in ref. [59]. However, in order to have a like-for-like comparison to the Nnlops band we have restricted the JetVHeto uncertainty estimate to the same class of variations. Notwithstanding this more limited evaluation of the uncertainties, we see the two sets of predictions still lie within each other's error bands, except for a barely visible excursion of the Nnlops outside the low edge of the JetVHeto envelope in the R = 1.0 plot ( fig. 7, bottom-right). In all cases the central predictions of the Nnlops and JetVHeto programs are never out of agreement by more than 5-6%, thus it becomes possible to meaningfully use the former in conjunction with the more comprehensive, conservative, uncertainties of the latter in real analysis.

Next-to-leading jet transverse momentum
We now briefly discuss the description of the second hardest jet given by our Nnlops generator. First of all, we remind the reader that our generator describes this distribution only at tree level accuracy, but certainly introducing incomplete higher order corrections. This is illustrated in fig. 8 where the second jet transverse momentum distribution obtained with the Hnnlo, the Hj-Minlo-NLO, the Hj-Minlo-LH and the fully showered Hj-Minlo . The jet veto efficiency, ε (p T,veto ), is defined as the cross section for Higgs boson production events containing no jets with transverse momentum greater than p T,veto , divided by the respective total inclusive cross section. In both plots shown above we display the jet veto efficiency as a function of the cut p T,veto . In the green shaded area, one can see the scale uncertainty band obtained from the Nnlops simulation (see Sect. 3 for details regarding this uncertainty estimate), with the NNLL+NNLO uncertainty band from the JetVHeto program [46,59] superimposed in red. The lower pane displays the same quantities as a ratio with respect to the central Nnlops prediction. The Nnlops predictions here were obtained with the default profile function (β = 1 2 ) used in determining the NNLO reweighting W (y, p T ).
generators. We see fair agreement between the Hnnlo and the Hj-Minlo-NLO predictions, as expected, with previously discussed differences in scale assignments accounting for discrepancies between the two in the high transverse momentum tail. On the other hand, the Hj-Minlo-LH result is markedly higher than the fixed order results. Since the second jet is certainly the Powheg hardest radiation in this case, we identify the cause of this increase as due to the fact that Powheg multiplies the hardest radiation spectrum by its NLO K -factor. Notice that in the case of R = 0.4, subsequent shower radiation leads to a transverse momentum spectrum that is in better agreement with the fixed order calculation. This is due to the fact that energy leakage outside the cone due to showering softens the spectrum. However, the fact that this effect competes with the K -factor effect, up to the point of nearly canceling it, has to be considered accidental. In fact, for R = 0.7, where the effect of energy leakage is negligible, we see no such compensation. We also remind the reader that the so called K -factor effect is formally of order α 5 S or higher, beyond our intended accuracy.
The scale uncertainty in the radiation of the hardest jet is usually underestimated in Powheg. In fact, the spectrum of the hardest radiation is generated with a scale independent procedure. This spectrum is multiplied by the underlying Born cross section, that has a next-to-leading order scale dependence. Thus, also the Powheg spectrum displays formally a next-to-leading order scale dependence. In the case of the Nnlops generator, the reweighting procedure constrains even more the underlying Born scale variation, forcing it to integrate up to the NNLO scale dependence. This is apparent from figs. 9 and 10. In fig. 9 we compare the Nnlops and Hnnlo results for the second hardest jet transverse momentum distribution. We use our default method to compute the uncertainty band in the Nnlops result, while in the Hnnlo case we present the standard 7-points scale variation taking µ F = µ R = m H for the central value. We see that the Hnnlo envelope includes the Nnlops one, the latter being considerably smaller. In fig. 10 we compare the Hj-Minlo showered result with our Nnlops output. We see that the main difference in the central value prediction is due to the fact that the NNLO K -factor is mainly applied when the transverse momentum of the hardest jet is below half  the Higgs mass, and that we expect the hardest jet transverse momentum to be just slightly above that of the second jet. The two results also tend to lie outside of each other's error band in the small transverse momentum region, a further indication that the scale variation uncertainty is too small for this tree level distribution. Furthermore, the Nnlops band is smaller than the Hj-Minlo one, as anticipated earlier.

Conclusions
In this paper we have presented a Nnlops implementation of Higgs boson production in hadronic collisions, thus yielding the first example of a NNLO calculation matched to a parton shower. We observe that our method is the NNLO extension of the current Nlops methods. In fact, in both MC@NLO and Powheg, NLO accuracy is achieved without the need of unphysical separation scales; NLO accuracy is granted for infrared safe observables, and logarithmic accuracy is maintained at the level of the shower Monte Carlo. These same features are achieved at the NNLO level by our procedure.
Presently the method has been applied to Higgs production in the large m t limit. The full inclusion of finite mass effects would require a calculation of Higgs plus jet production at NLO including finite mass effects, which is currently not available. It is possible, however, to include such mass effects at least at O(α 3 S ), if we are allowed to treat them as small corrections. One could incorporate such effects by reweighting the Hj-Minlo events by the ratio of the full m b , m t mass dependent Higgs+jet cross section at O(α 3 S ) over its large m t limit, computed at the underlying Born level. Then one could proceed as in the present work, reweighting the events using the recent HNNLO-V2 calculation [60], which includes mass effects up to order O(α 3 S ). We leave this possibility to future studies. Our procedure relies upon the results of ref. [39], in which a method was presented for extending an Nlops simulation of Higgs plus jet production, so as to simultaneously deliver NLO accuracy for inclusive quantities. At present, the method of ref. [39] has been applied, besides the Higgs case, to the Drell-Yan processes and the Higgsstrahlung process [61]. It can be generalized easily to all reactions involving the production of a heavy colourless system. For more general processes, we see no obstacles to the implementation of the procedure of ref. [39], except for the increased complexity in the calculation of the needed resummation formulae. Thus, the present method for building Nnlops generator may be extended to more complex processes, subject to the availability of the corresponding NNLO calculation, and of the suitable extension of the procedure of ref. [39].

Acknowledgements
We gratefully acknowledge Carlo Oleari for important contributions made earlier in the development of the Minlo method. It is a pleasure to thank Andrea Banfi, Massimiliano Grazzini and Gavin Salam for helpful comments and discussions. GZ is supported by Science and Technology Facility Council. We acknowledge ESI (ER and GZ) and KITP (GZ) for hospitality while part of this work was carried out.
A Inclusivity of the Higgs p T spectrum Following sect 4.3 of ref. [41], suppressing indices for simplicity, the difference between conventional NLO and Powheg predictions for an observable O, O PWG = O NLO +δ O , is given by where the real phase space Φ, is factorised into that of the underlying Born kinematics, Φ B , and its complement parametrising those of the emitted parton, Φ rad . B and B refer to the leading order cross section differential in the Born variables and, respectively, its NLO equivalent. R denotes the real cross section and ∆ the Sudakov form factor associated to emission from the underlying Born. From eq. (A.1) it is clear that Powheg is NLO accurate for arbitrary infrared safe observables, up to NNLO sized contributions; the factor containing the difference of the observable functions in the rightmost bracket allowing one to effectively replace the Sudakov form factor by 1. In addition, it is also clear from eq. (A.1), that the more inclusive the observable is with respect to the underlying Born kinematics the greater will be the effect of the rightmost bracket in nullifying δ O . Indeed, if the observable is fully inclusive with respect to the Born kinematics δ O vanishes identically. The Higgs p T is not simply a function of the Born kinematics, Φ B , as implemented in the Hj-Minlo code, but it does have the property that it is rather inclusive, allowing the rightmost bracket in eq. (A.1) to significantly diminish the difference between the Minlo enhanced fixed order prediction and that of the hardest emission cross section, as seen in fig. 2.