NLO Higgs boson production plus one and two jets using the POWHEG BOX, MadGraph4 and MCFM

We present a next-to-leading order calculation of Higgs boson production plus one and two jets via gluon fusion interfaced to shower Monte Carlo programs, implemented according to the POWHEG method. For this implementation we have used a new interface of the POWHEG BOX with MadGraph4, that generates the codes for generic Born and real processes automatically. The virtual corrections have been taken from the MCFM code. We carry out a simple phenomenological study of our generators, comparing them among each other and with fixed next-to-leading order results.


Introduction
The search for the Standard Model Higgs boson is entering the endgame phase. Since the start of the LHC running, the allowed mass range for the Standard Model Higgs boson has been greatly curtailed. Furthermore there are tantalizing, but inconclusive, hints in the remaining low mass region [1,2]. Should these hints be confirmed by further data, it will be a matter of some urgency to examine the properties of the new state (or states). In the low mass region, the Standard Model Higgs boson is predicted to have about ten decay modes with branching fractions greater than one per mille, so there will be a number of channels to be studied.
In order to extract the Higgs boson couplings to fermions and gauge bosons from these different channels, detailed information about the production rates will be required [3,4,5]. In this context, the Higgs boson + 2 jet process enters principally as an irreducible background for the vector boson fusion (VBF) process [6,7,8]. This process is of particular importance in the extraction of the HW W and HZZ couplings and in confirming if the detected scalar particle is the Higgs boson responsible for electroweak symmetry breaking. In order to distinguish VBF Higgs boson signal from backgrounds, stringent cuts are required on the Higgs boson decay products as well as on the two forward quark jets which are characteristic for VBF. Information on the CP properties of the Higgs boson can be extracted by studying the azimuthal distributions of the two hardest jets [4,9], in events with a large ( 3 units) rapidity separation between them. In fact, the signature of the VBF Higgs boson production are two jets well separated in rapidity. In addition, we expect very low hadronic activity between the two hardest jets, due to the exchange of colourless vector bosons in the t channel, in contrast to what is expected for the gluon-fusion production mechanism. For the gluon-fusion processes, the higher-order corrections are significant [10,11,12,13], and a detailed understanding of the structure of the radiation pattern is needed in order to access the efficiency of the central-jet veto [14,15].
The Higgs boson search channels are often subdivided according to the number of associated jets, for a number of reasons. First, a differing number of associated jets can imply a different source of background. For example, in the search for a Higgs boson decaying to W + W − , the backgrounds have different origins in the differing jet bins, so it proves advantageous to analyze the different jet multiplicities separately. Second, in the two-jet bin, the new VBF production channel will come into play. For these reasons it is pressing to provide NLO predictions for gluon-fusion initiated Higgs boson + jets processes in the POWHEG formalism, so that experimentalists can incorporate the best theoretical information into their analyses.
Phenomenological studies for the production of a Higgs boson in association with two jets, at the parton level and at fixed order, have been available for a while. In refs. [14,15], the calculation was performed at leading order, with the gluons coupling to the Higgs boson via a top-quark loop, retaining the exact top-quark mass dependence (m t ) in the whole calculation. In that paper it was shown that the large-m t limit provides an excellent approximation to the full m t dependence when the Higgs boson mass, m H , is small compared to the top-pair threshold. The large-m t limit was found to break down for m H > m t and when jet transverse momenta become large (p j T m t ). However, large dijet invariant masses do not invalidate the m t → ∞ limit, as long as the Higgs boson mass and the jet transverse momenta are small enough, less than the top-quark mass in practice. This observation opened the possibility to compute the NLO corrections to Higgs boson + 2 jet production in gluon fusion, in the large-m t limit, starting with an effective coupling of the Higgs boson to the gluon field [12,13] In this work we have implemented Higgs boson plus one-and two-jet production in gluon fusion, in the large top-quark mass limit, using the POWHEG method. While for the first process a matched NLO+shower calculation has already appeared in the literature [16], the second one has never been performed before. We have built our generators using the POWHEG BOX framework [17], with the virtual corrections taken from the MCFM program [18], and the Born , colour-correlated Born, spin-correlated Born and real contributions computed using a new POWHEG BOX interface to the MadGraph4 [19,20] program, which is fully generic and can be applied to any process. Thus, the aim of this paper is twofold: • To present a new generic interface of MadGraph4 to the POWHEG BOX that allows for the automatic generation of the code for the Born amplitude, the Born colour-and spincorrelated contributions, for the real-radiation amplitude and for the Born colour flow in the leading-colour approximation. These ingredients are all needed in the implementation of a new process in the POWHEG BOX. This interface can be used for any process that can be generated with MadGraph4. Thanks to this interface, in order to construct a POWHEG BOX generator, one needs only to provide the Born phase space and the virtual corrections.
• To illustrate the Higgs plus one and two jet generators, by comparing their outputs to the corresponding NLO results and amongst themselves.
The paper is organized as follows: in sec. 2, we present the new interface of the POWHEG BOX to MadGraph4. In sec. 3, we discuss the approximations used in our calculation and how we extract the virtual corrections from the MCFM code. In sec. 4 we give more details about the phase-space generator for Higgs boson plus one and two jets, and how we deal with the divergences present at the Born level. After briefly recalling a few crucial POWHEG BOX parameter settings, we discuss some phenomenology in sec. 5. In particular we compare the NLO differential cross sections with the results obtained with the POWHEG BOX, at several levels of approximation. Finally, in sec. 6, we summarize our findings. A few technical details of the MadGraph4 interface and of the PYTHIA setup are collected in the two appendices.

The interface to MadGraph4
The MadGraph4 package [19,20] generates squared tree-level matrix elements in an automatic way. It is based on Feynman diagrams, using HELAS routines [21] as building blocks which are subsequently combined into colour-ordered amplitudes and written in a FOR-TRAN code. Furthermore, the MadDipole [22,23,24] extension of MadGraph4 and the MadFKS tool [25] can generate the spin-and colour-correlated Born squared amplitudes. The MadGraph4 package can thus construct all the ingredients needed to interface a NLO calculation to the POWHEG BOX (see the complete list in the introduction of ref. [17]), with the exception of the Born phase space and the virtual matrix elements. This makes the MadGraph4 framework an ideal environment to complement the POWHEG BOX in the creation of all tree-level squared matrix elements and in this section we describe a newly developed interface between these two codes. For a given process, the routines and parameters that are generated by MadGraph4 for the POWHEG BOX are (see ref. [17] for more details on the notation and conventions used): • the multiplicity of the Born and real-emission processes, nlegborn and nlegreal, respectively, defined in the nlegborn.h include file; • the list of Born and real-emission flavour structures flst_born(k=1,nlegborn,j=1:flst_nborn), flst_real(k=1,nlegreal,j=1:flst_nreal), respectively, to be initialised by a call to the init_processes routine; • the routine setborn(p(0:3,1:nlegborn),bflav(1:nlegborn),born, bornjk(1:nlegborn,1:nlegborn),bmunu(0:3,0:3,1:nlegborn)) which computes, for a given set of four-momenta p and flavour structure bflav, the Born squared matrix element born, the colour-correlated, bornjk, and the spincorrelated one, bmunu; • the routine that computes the real emission squared matrix elements amp2 (stripped of a factor α S /(2π)), for a given momentum configuration p and flavour structure rflav setreal(p(0:3,nlegreal),rflav(1:nlegreal),amp2); • a colour flow assignment to the Born squared matrix elements, in the leading colour approximation, needed by the showering programs. The colour is assigned on statistical grounds, based on information that is cached during the call to the squared matrix elements, by the routine borncolour_lh. In addition to this routine, the interface also provides a similar one for the real squared matrix elements, realcolour_lh. In the present version of the POWHEG BOX this routine is not used; instead a different method is employed to assign the colour to the real-radiation matrix elements, as discussed in sec. 8 of ref. [17]; • an interface to the Les Houches parameter input card to specify the physics-model parameters through the routine init_couplings.
In contrast to the default way of generating amplitudes with MadGraph4, all the squared matrix elements for the various flavour structures are written in the same directory, making sure that the files and routines have different names. This allows one to compile the code into a single library that contains all the matrix elements for all the flavour structures. An interface is provided that concatenates the flavour vectors bflav or rflav into strings, which are used as unique identifiers of each of the squared matrix elements. More technical details on the use of the MadGraph4 interface to the POWHEG BOX can be found in appendix A.

Higgs boson production in gluon fusion
The amplitudes for a Higgs boson in association with three, four or five partons, needed for the Higgs boson plus one and two jet cross sections at NLO, are calculated using an effective Lagrangian to express the coupling of the gluons to the Higgs field [26] where the trace is over the colour degrees of freedom. At the order required in this paper, the coefficient C is given in the MS scheme by [27,28] where v is the vacuum expectation value of the Higgs field, v = 246 GeV. We have used the automatic interface described in sec. 2 to generate the code for the Born amplitude, the Born colour-and spin-correlated contributions and for the real cross section. In appendix A we give some more details on the adopted procedure. The one-loop amplitudes for the Higgs boson plus three and four parton processes are extracted from MCFM as described below.
The interface between MCFM and the POWHEG BOX is fairly straightforward. Since in the POWHEG BOX implementation the Born, real and subtraction terms are available independently, MCFM needs only to return the pure virtual contributions. The interface to MCFM transfers the electroweak parameters, scale and scheme choices. Once these are established, calls to the virtual routines of MCFM can be made on an event-by-event basis.
Since the routines that fill the values of the electroweak parameters (and the other information) are generic, this interface could be fairly easily extended to include other processes currently implemented in MCFM. However the normal MCFM routines are designed to return matrix elements that are summed over the flavour of the final-state partons. Therefore, in order to correctly interface to the POWHEG BOX, one must make small modifications to the MCFM code such that it returns the matrix elements for individual final-state flavour combinations.

Implementation of the Hj and Hjj generators
In the following, we use the notation H, Hj and Hjj to refer to Higgs boson generators that, at the Born level, describe the production of a Higgs boson plus zero, one and two partons.
Using the automated MadGraph4 interface described in sec. 2, we have built the routines for the Born amplitude, the Born colour-and spin-correlated contributions and the real squared amplitudes, directly in the format required by the POWHEG BOX. The generation of the amplitudes is fast, taking only a few minutes for the Hjj process. The virtual corrections were extracted from the MCFM code, as described in sec. 3.1. At this point, the only missing ingredient for completing our generators is the Born phase space.

Phase space for the Hj generator
The phase space routine for the Hj Born process is trivial, since it is just a two-body phase space. It has been implemented with the possibility to activate an optional cut on the transverse momentum of the Higgs boson, and with an effective importance sampling of the small transverse-momentum region. The cut is necessary if one wants to generated unweighted events. Alternatively the cut can be set to zero, and the cross section is weighted with a suppression factor, that is a function of the underlying-Born kinematics. This factor is equal to where p T is the Higgs boson transverse momentum in the underlying-Born configuration and p 0 is set by the user 1 . Events are then generated uniformly in the cross section times F , and with a weight proportional to 1/F . The normalization of the weight is such that the cross section for events passing a set of cuts is given by the sum of all weights for the events passing the cut divided by the total number of generated events.

Phase space for the Hjj generator
The phase space for the Hjj generator was built using the same factorized phase space that POWHEG uses for the real kinematics. In other words, we treat the Hjj Born phase space as the real phase space for the Hj process with one extra emission. We have allowed the possibility of performing importance sampling in the regions where the Born cross section becomes singular. In more detail, labelling 1 and 2 as the two final-state light partons in the Hjj process, we write the Hjj phase space using the following identity dΦ H12 = dΦ H12 N d 12 where E i is the energy of the i-th parton in the center-of-mass frame, and where d 1 , d 2 and d 12 are phase-space functions that vanish respectively when parton 1 is collinear to the initial-state partons, parton 2 is collinear to the initial-state partons, and partons 1 and 2 are collinear to each other. Their precise form is given in eqs. (4.23)-(4.26) of ref. [17]. We then factorize each phase-space factor in (4.2) according to the formula The subscripts H1 or H2 characterize the underlying Born, while the superscript in Φ rad specifies the radiation process. Thus, for example, dΦ H2 is the underlying-Born phase space of the Higgs boson together with parton 2, and dΦ (21) rad is the radiation phase space corresponding to parton 2 emitting parton 1. The notation dΦ  partons 1 or 2 are emitted by the initial-state partons. The factorization of the phase space into underlying Born and radiation phase space for both initial-and final-state radiation is the default one used in POWHEG, and is described in detail in secs. 5.1 and 5.2 of ref. [36]. The decomposition in eq. (4.4) is such that appropriate importance sampling is performed in all singular regions. In fact, the factors in each term damp all but one singular region, and the corresponding factorized phase space performs importance sampling precisely in that region. Ideally, the phase-space integration should be performed with an integrator that can sum over a discrete variable. The POWHEG BOX integrator [37] does not have this feature at the moment. Thus, we divided the range of one extra integration variable into four segments, mapping each segment to one of the phase space components.
The phase space of eq. (4.4) is not the default one used in the Hjj generator. It is activated by setting the variable fullphsp in the powheg input file. The default phase space is simply with no importance sampling at all. We have in fact observed that the loss of efficiency due to the increased number of calls to the matrix-element routines overwhelms any benefit arising from the improved importance sampling.
As for the case of the Hj generator, the Hjj generator also includes the implementation of a Born suppression factor for the suppression of the singularities of the underlying-Born amplitude, in the born suppression routine. It has the form where p 0 is a parameter that characterizes the minimum jet energy where some accuracy is required, p T 1 and p T 2 are the transverse momenta (with respect to the beam axis) of the two final-state partons, and that can be interpreted as the transverse momentum of parton 1 with respect to parton 2 or vice-versa, depending upon which is the softest. This suppression factor can also play the role of a generation cut, if p 0 is chosen small enough.

The POWHEG BOX parameter setting
We have turned on the bornzerodamp flag by default in both the Hj and the Hjj generators. The purpose of this flag is explained in ref. [17]. This results in a considerable speedup of the Hjj code. However, no appreciable differences were observed in the results obtained without the bornzerodamp option. The separation of the singular regions is controlled in the POWHEG BOX by the parameters par diexp and par dijexp that are set to 1 by default (see ref. [17], section 4.7). For the results presented here we have chosen to set them to 2, which leads to slightly better stability for the set of distributions that we have considered. We note that these parameters were also set to 2 in the POWHEG BOX dijet generator [38]. Higher values of these parameters correspond to sharper separation of singular regions.

Phenomenology
In this section, we present a phenomenological study of our generators. This study does not aim to produce results with realistic experimental cuts, but rather to compare the predictions of the generators for some key observables. In particular, we compare the generator results to the fixed NLO result and with the POWHEG output at the level of the generation of the hardest radiation (i.e. before interfacing the result to a parton shower program), after the shower with no hadronization effects, and after hadronization.
We present results for the LHC running at 7 TeV, computed using the CTEQ6M parton distribution function (pdf) set [39]. The same calculations can easily be performed with any other available set [40,41], but a study of pdf effects is beyond the scope of the present paper. Jets are reconstructed using the anti-k T jet algorithm [42], with R = 0.5 and default recombination scheme. Cuts of 20, 50 and 100 GeV on the final-state jets are considered.
The factorization-and renormalization-scale choice deserves a more detailed discussion. In view of the large NLO corrections for these processes, the scale dependence is quite large. This is also a consequence of the fact that Higgs boson production in gluon fusion starts at order α 2 S , so that the Hj and the Hjj processes are of order α 3 S and α 4 S , respectively. While for H production the natural scale choice is of the order of the Higgs boson mass, in the case of Hj production one may have a two-scale problem, if the transverse momentum of the jet is much smaller or much larger than the Higgs boson mass. In addition, besides the Higgs boson mass, one may also consider the Higgs boson transverse momentum, which may be more appropriate for very small or very large Higgs boson transverse momentum p H T , or the Higgs boson transverse mass m H T = m 2 H + (p H T ) 2 , which may be more appropriate for large Higgs boson transverse momentum. For the Hjj case, there are even more possibilities. Although a full study of scale dependence would be very valuable, we will not perform it in the present work. In the Hj case we will limit ourselves to the fixed Higgs boson mass scale choice and to transverse momentum of the Higgs boson in the underlying-Born configuration. In the Hjj case, we will consider the fixed Higgs boson mass scale choice and to theĤ T scale choice at the underlying-Born level, wherê and p T i are the final-state parton transverse momenta in the underlying-Born kinematics. All results shown in the next sections have been computed with m H = 120 GeV and Γ H = 0.00575 GeV. The Higgs boson momentum has been generated distributed according to a Breit-Wigner function, with fixed width. In the following, we will label "LHE" the results obtained at the level of the POWHEG first emission, "PY" the results obtained with the POWHEG+PYTHIA combination with the hadronization and underlying event switched off, and with "NLO" the fixed NLO results.
For some observables, there is overlap among the various generators. For example, the Higgs boson transverse momentum, as well as the one-jet multiplicity and the transverse momentum of the leading jet, are described by both the H and the Hj generators.
The H generator gives a plausible description for these quantities also at small transverse momenta, since it provides a resummation of transverse-momentum logarithms that the Hj program does not provide. On the other hand, at large transverse momenta, the Hj generator has NLO accuracy, something that the H generator does not have. For the same quantities, the Hjj generator cannot be used, since it requires the presence of a second jet. The two-jet multiplicity, as well as the second hardest jet transverse-momentum distribution, are provided by both the Hj and the Hjj programs, again with a different level of accuracy. A comparison of the generators in the regions where they overlap will be carried out in sec. 5.4.

Results for Hj production
We have generated a sample with 2M events at the fixed scale and 7.5M at the running scale. The event generation time is approximately 25 seconds for 1000 events on a typical CPU.
We have simulated Hj production with two scale choices: µ F = µ R = m H and the p T of the underlying-Born parton µ F = µ R = p UB T . All the following plots will come in pairs, with the left plot referring to µ F = µ R = m H , and the right one to µ F = µ R = p UB T . We compare the fixed NLO results, the POWHEG hardest-emission results (LHE) and POWHEG+PYTHIA (PY) ones, where the hadronization and underlying-event effects in PYTHIA have been turned off (see appendix B). We show results for a few physical observables with three different cuts applied to the final-state jets: p T cuts of 20, 50 and 100 GeV. In fig. 1 we compare the jet multiplicity for one and two jets, at the NLO, LHE and PY levels, for jet cuts of 20, 50 and 100 GeV. We notice that the one-jet multiplicity is similar in the three results. The PYTHIA multiplicity is smaller than the LHE, which can be understood as a result of showering off the first jet. More marked differences can be seen in the cross sections with two jets: here the LHE result is clearly larger than the NLO one, showing however a very different pattern as a function of the jet cut, depending upon the scale choice. We will clarify the origin of these patterns later, when we discuss the transverse-momentum spectrum of the jets. We point out, however, that the NLO result for the two-jet multiplicity is of order α 4 S , and thus has a marked scale dependence. In the left plot the scale is equal to the Higgs boson mass, whilst in the plot on the right the scale is of the order of the minimum jet p T . For a jet cut of 20 GeV these scales are widely separated, which explains the large difference in the NLO results. In fig. 2 we show the transverse momentum of the Higgs boson. The LHE, the PY and the NLO results are in remarkable agreement for this quantity, which is expected, due its inclusiveness. The Higgs boson rapidity distribution is plotted in fig. 3 while the  hardest jet p T in fig. 4. Both these quantities display a behaviour similar to the Higgs boson transverse momentum, again expected due to their inclusiveness. We now turn to more exclusive quantities, beginning with the transverse momentum Thus, for small p T , the NLO corrections are smaller than for large p T . In POWHEG, one first generates the underlying-Born configuration, according to the cross sectionB, which is the total inclusive cross section at a fixed underlying-Born configuration. The radiation kinematics is then generated using a shower technique. As a result the whole distribution for the radiation is amplified by a K-factor equal tō B/B [43,44,45]. In this case, due to the growth of the NLO corrections as a function of p T , the amplification of the radiated-jet distribution due to theB/B K-factor increases as a function of the transverse momentum, a trend that is visible in the right figure. Conversely, no such effect is present for fixed scales. However, in this last case, one should recall that, in the LHE events, one power of α S is effectively evaluated at the transverse momentum of the jet rather than at the fixed scale, thus yielding a decrease in the cross section which is also visible in the left plot. The larger value of the LHE cross section with respect to the NLO one is related to the large K-factor, i.e. is due to the fact that the hardest radiation is amplified by a factorB/B.  Figure 6 is similar to fig. 5 except that the LHE result is obtained setting the bornonly flag to true in the POWHEG BOX input file. When this flag is true, the POWHEG BOX generator uses the Born cross section, rather than the NLOB function, to generate the underlying-Born configuration. We can see that, in this case, the LHE result is much closer to the NLO result, except for small transverse momenta, where the Sudakov damping becomes manifest. This results confirms that the enhancement of the transverse-momentum distribution of the second hardest jet is indeed due to theB/B K-factor.
In fig. 7 we show the second hardest jet transverse momentum for different cuts on the first jet p T . From the figure it is clear that, when the second jet has transverse momentum above the first jet p T cut, the leading jet is forced to have larger transverse momentum, and the cross section falls more rapidly. The comparison among the LHE, PY and NLO curves shows the typical pattern, with the NLO diverging at small transverse momenta, the LHE being instead suppressed in that region, but raising above the NLO result because of theB/B K-factor.  In fig. 8 we plot the relative transverse momenta of the partons within the hardest jet p rel,j 1 T . This quantity is defined as follows: we perform a longitudinal boost to a frame where the rapidity of the first jet j 1 is zero. In this frame we compute where k j is the momentum of the first jet, and p i the momenta of the partons clustered within the first jet. This quantity, when computed at the LHE level, displays a marked difference from the corresponding NLO result. As discussed in ref. [38], this can be easily understood if we remember that the LHE result is suppressed by a Sudakov form factor, which requires that no harder radiation has been emitted from either the initial-or the final-state partons. The large bin at p rel,j 1 T = 0 in the LHE result is due to events in which there is only one parton in the jet, most likely initial-state radiation events. On the other hand, the hardest radiation provided by PYTHIA should restore a shape closer to the NLO result, although amplified by theB/B factor. Further amplification of the showered result is due to the fact that the shower uses a running coupling evaluated at the scale of the radiation, while the NLO result uses higher scales, and to the presence of multiple emissions in the shower. It is clear that this distribution, being determined mostly by the shower program, is quite sensitive to the shower model and tuning, and to the interface between the shower and POWHEG BOX.

Results for Hjj production
We have generated a sample with 2.5M events both at fixed and running scales. The event generation time is approximately 25 minutes for 1000 events on a typical CPU.
We have run the Hjj program for two scale choices, µ F = µ R = m H and µ F = µ R =Ĥ T . All the following plots will come in pairs, the left one referring to the first scale choice, and the right one referring to the second scale choice. We begin by showing in fig. 9 the two-and three-jet multiplicities in the Hjj process. We find considerable agreement of the LHE, PY and NLO output for both the two-and three-jet multiplicities for the left plot, which corresponds to the choice of scale µ F = µ R = m H . The right plot corresponds to theĤ T choice of scale. In this plot the two-jet multiplicity agrees for all cuts, while the three-jet multiplicity displays marked differences between the LHE results and the NLO ones. For the 100 GeV transverse-momentum cut, the fully showered result and the NLO also differ. The reason for these differences is the following: with the m H scale choice the K-factor is near one, and thus theB/B amplification that usually enhances the spectrum of the radiated jet (that in this case is the third jet) is not effective. TheĤ T scale is considerably larger than m H , especially with a large p T cut, which implies a reduction of the Born cross section and an increase of the K-factor. In fig. 10 we show the transverse-momentum distribution of the Higgs boson in a Hjj sample with a 20 GeV p T cut on the two hardest jets. Good agreement is found between the LHE, PY and NLO curves, except for small Higgs transverse momenta, where differences of the order of 30% are found. This is not surprising, in view of the required presence of two relatively soft jets. Radiation off these jets is Sudakov suppressed in the LHE result with respect to the NLO one. Since this radiation would deplete the jets, the LHE result suffers less depletion, and is thus larger. On the other hand, the shower degrades the energy of the jets lowering the cross section, an effect visible in the PY result. For large Higgs boson p T , at least one of the jets is forced to be hard, and thus these effects lose importance. In fig. 11 we show the Higgs boson rapidity distribution in Hjj for p T cuts of 20, 50 and 100 GeV on the two hardest jets. Again, since this is an inclusive distribution, it displays good agreement between the LHE, PY and NLO results. The ratio is only displayed for the 20 GeV cut. The three predictions become even more consistent for larger p T cuts, as one would expect.  In figs. 12 and 13 we display the transverse-momentum distribution for the first and second jet. The first distribution bears some similarity to the Higgs transverse momentum, for inclusiveness reasons. Observe that the small momentum disagreement between the LHE and the NLO result observed for the Higgs transverse momentum case is less evident in the hardest jet plot, and even less so in the second jet spectrum. This is explained by the fact that the same point in the abscissa corresponds to an increasing hardness of the event in the Higgs boson, first jet and second jet p T distributions.
Turning to less inclusive quantities, we show in fig. 14 the transverse-momentum distribution of the third hardest jet. We recognize here the typical behaviour of the LHE results, with the damping of the small-momentum growth found in the NLO result, and with the typical increase due to theB/B factor at higher momenta. As already anticipated in the discussion on the jet multiplicity, we notice that the increase is modest for the m H choice of scale, but is very large for theĤ T scale. In fig. 15 we display the third-jet transverse momentum using three different cuts on the p T of the first and second jet. The pattern is similar to what was already observed in the Hj case. In the left plot, where the scale is chosen equal to m H , we see a better concordance of the LHE and NLO results as the transverse momentum increases.
In figs. 16 and 17 we show the p rel,j T distribution for the hardest and second hardest jet respectively. Here too the pattern is similar to what already observed for the Hj case. The LHE distribution is strongly suppressed, due to the fact that small p rel,j T values also imply that no initial-state radiation has taken place above that scale, and the distribution is dominated by the shower effects, that are not obtained with the same scale choice as the NLO result. We note again that theB/B K-factor plays a role here, especially for the  right plots, that have a larger value because of the use of theĤ T scale.

Hadronization effects and matching ambiguities
In this section we focus upon two topics: the effects of hadronization and the ambiguities related to matching the LHE result to the shower. In figs. 18 and 19 we display the transverse-momentum distribution of the radiated jet for different cuts on the first jet p T , for the Hj and Hjj generators respectively.
The curves labelled "PYscalup" are obtained with a non-default determination of the scalup parameter to be set in the Les Houches interface, which limits the hardness of the radiation of the following shower. While normally in POWHEG scalup is set to the hardest momentum of the radiation that POWHEG generates, in the "PYscalup" results, we determine it by finding the smallest transverse momentum of the LHE, which can be either the transverse momentum of any parton relative to the beam axis, or of any parton relative  to any other parton, computed in the center-of-mass frame. Because of the way the realradiation contribution is separated into singular contributions, the two choices differ only for subleading configurations, and one expects only a minor effect due to this change. The plots in both figures confirm this expectation. Similarly, we see that hadronization and underlying-event effects (indicated with HAD in the figures) have a sizable impact on the distributions only for small transverse momenta, as expected.

Comparison between the H, Hj and Hjj generators
In this section, we compare a few distributions that are described by more than one available POWHEG BOX generator. This comparison can be considered as a first step in the direction of merging POWHEG BOX samples with increasing number of jets. We first consider a comparison between the Higgs boson transverse-momentum spectrum obtained using the H and the Hj generators. The H generator describes this dis- tribution in the whole p T range, with the Sudakov region obtained with next-to-leading logarithmic accuracy, and with leading-order accuracy for the high-momentum region. On the other hand, the Hj generator describes this distribution with NLO accuracy, but only for transverse momenta above the Sudakov region. In these comparisons, the H sample was obtained following the recommendation of ref. [46], i.e. the parameter hfact was set to m H /1.2 in the powheg.input file. With this setup, the resulting Higgs boson p T distribution is in remarkable agreement with the output of the HqT program [47, 48,49,50]. In addition, we apply a K-factor of 1. 32 fig. 20, which demonstrates the agreement between the results for the H sample and the output of the HqT program, for the Higgs boson p T distribution [46]. This comparison is performed using a scale equal to the Higgs boson mass, since HqT accepts only a fixed renormalization and factorization scale. We also show two predictions for the Hj generator, using both the fixed and running-scale choices. The Hj generator predictions with a fixed scale are also in excellent agreement with HqT for large transverse momenta. This is not surprising, since, in this kinematic region, it is using the same O(α 4 S ) result as HqT. At low p T , the Hj results begin to feel the lack of soft-gluon resummation effects that are included in HqT. The Hj result computed with the dynamical scale shows a substantially similar pattern, undershooting the HqT one by about 25%.
Turning to the comparison of the Hj and Hjj generators for the transverse-momentum distribution of the second hardest jet, we notice that the Hj generator computes this distribution down to small values of the transverse momentum, since it includes the appropriate Sudakov effects. At large transverse momenta, however, it only has leading-order accuracy.
The Hjj program has instead NLO accuracy at large transverse momenta, but does not include fully resummed Sudakov effects at small p T . We do not have any higher-accuracy calculation for this distribution, as in the previous case. The results are displayed in fig. 21. Notice that, for the m H scale choice, the matching is very good, something that we expect since the K-factor is close to one in this case. Matching with the running scales leads, instead, to non-negligible differences, although a stability plateau is present roughly above 60 GeV, provided that the transverse momentum is not too large.

Conclusions
In the present work we have developed generators for the gluon-fusion production of a Higgs boson in association with one or two jets, in the large top-quark mass limit. We have examined the output of the two generators, comparing them with fixed NLO calculations, among each other, and with the fully inclusive Higgs boson POWHEG generator. The features of the distributions we have examined are all well understood and reflect our expectations for a typical POWHEG generator. In general, we find that quantities inclusive in the radiated jet (i.e. the second jet in H + 1 jet and the third jet in H + 2 jets) are in good agreement with the NLO result, while distributions in the radiated jet reflect the features of Sudakov suppression and NLO enhancement that are typically found in NLO generators matched with a shower. We see no indication of problems related to the increase in complexity when going from the inclusive Higgs boson production to the associated production with one and two jets, other than an increase in the amount of computer time required for the calculation.
The development of the Higgs boson production code has been achieved using a new interface to the MadGraph4 code, that has also been presented in this work. The use of this interface has considerably simplified the construction of the generator. The interface is fully generic, and so we expect that the implementation of new processes will be greatly simplified with its use. In addition, the code takes advantage of compact expressions for Higgs boson plus four parton virtual amplitudes that had previously been collected in MCFM.
The code of our new Higgs boson production generators can be accessed via the POWHEG BOX svn repository (see the POWHEG BOX web page http://powhegbox.mib.infn.it for instructions). The new MadGraph4 interface is also available there, so that people willing to develop their own POWHEG code may benefit from its use.

A. The interface with MadGraph4: technical details
The MadGraph4 interface is available under the POWHEG-BOX/MadGraphStuff directory. To use the interface, create a process directory in the POWHEG-BOX directory, and copy the whole content of the MadGraphStuff directory to this new process directory. To specify the process, edit the file Cards/proc_card.dat and set the process and the physics model as in MadGraph4. Always enter the real emission process, i.e. the Born process plus an extra jet j. For example, to generate Higgs boson plus two jets at a proton-proton collider, we entered p p > H j j j QCD=3 QED=0 HIG=1 It is recommended to set the parameters QCD and QED to be exactly equal to the number of strong and electroweak interactions in the real-emission process (excluding the couplings present in the effective vertex, if any). In the case of the heft model (see below), it is also needed to set the parameter HIG=1 to allow for the effective Higgs boson to gluon coupling to be included in the Feynman diagrams.
The ordering of the particles in the process should follow the POWHEG BOX conventions: 1. first particle: incoming particle with positive rapidity 2. second particle: incoming particle with negative rapidity 3. from the third particle onward: final-state particles ordered as follows • colourless particles first, • massive coloured particles, • massless coloured particles.
The default MadGraph4 interface has been validated to work with the following three physics models: • sm: this is the default MadGraph4 model with a massive top and bottom quark, diagonal CKM matrix and massless electrons and muons.
• smckm: this is the same as the sm model, however, the full CKM matrix can be specified in the parameter card.
• heft: this is the same as the default sm model with the inclusion of the effective coupling between gluons and the Higgs boson in the large top-quark mass limit. The Higgs boson Yukawa interaction with bottom quarks is neglected in this model.
Although the interface has not been tested with other physics models, it is straightforward to extend the interface to work with any simple extension of the Standard Model.
To generate the process, execute the NewProcess.sh script. This will compile the MadGraph4 code, generate the (correlated) Born and real-emission squared amplitudes and create the libraries that can be linked to the POWHEG BOX. In particular, libmadgraph.a contains the matrix elements, libdhelas3.a the HELAS routines and libmodel.a the physics model. Furthermore, the following files are written in the current directory: • Born.f: it contains the routines setborn to compute the (correlated) Born squared matrix elements and borncolour_lh that assigns the colour flow to a Born process.
• real.f: it contains the routines setreal to compute the real-emission squared matrix elements.
• Cards/param_card.dat: this is the input card where all the model parameters need to be specified.
• init_couplings.f: this contains the init_couplings routine that sets all the model parameters specified in the param_card.dat. The coupling constants that depend on the strong coupling st_alpha or from the event kinematics should be specified in the routine set_ebe_couplings, that is updated event-by-event.
• coupl.inc: it contains the common blocks for all the couplings used by MadGraph4.
To complete the implementation of a process in the POWHEG BOX, the user must provide the Born phase-space in the file Born_phsp.f and the virtual squared matrix elements in the file virtual.f. Also, no information on possible intermediate resonances in the matrix elements is kept. The user needs to specify explicitly in the routine finalize_lh (in the Born.f file) which resonances should be written in the LHE file, so that the shower Monte Carlo program can deal with them correctly. Finally, the resulting code can be compiled by executing the command $ make pwhg_main

B. PYTHIA setup
The sequence of PYTHIA calls we have used in the calculation of the results presented in sec. 5 is the following: • without hadronization