Heavy-quark mass effects in Higgs plus jets production

We study the production of a Standard Model Higgs boson in the gluon-fusion channel at the 13 TeV LHC. Our results are accurate to the next-to-leading order in QCD, bar for the lack of some two-loop amplitudes, for up to two extra jets and are matched to the Pythia8 Monte Carlo. We address the impact, at the level of inclusive rates and of differential distributions, of the merging of samples characterised by different final-state multiplicities, and of the effects induced by top and bottom masses through heavy-quark loop diagrams. We find that both the merging and the heavy-quark masses must be included in the calculation in order to realistically predict observables of experimental interest.


Introduction
The properties of the narrow-width particle discovered in 2012 by ATLAS [1] and CMS [2] at the LHC have so far been found to be fully compatible with those of the Higgs boson in the Standard Model, where the main production mechanism proceeds through heavyquark (denoted by Q henceforth) loops in gluon fusion. The latter is an example of a so-called loop-induced process, in which squared one-loop amplitudes enter already at the Born level. The Higgs is attached to the loops by means of an HQQ vertex, hence the requirement that the quark be heavy, in order to have a sizable Yukawa coupling y Q . In the SM, contributions where Q is the top quark are thus largely dominant; however, those where bottom quarks circulate in the loops are not entirely negligible either, predominantly because they can interfere with the Q = t ones. Incidentally, this mechanism is therefore responsible for a much more involved situation in BSM theories that feature a larger y b /y t ratio w.r.t. that of the SM.
Ignoring for the moment these issues, and considering only the Q = t contributions, a very significant simplification may occur in the computation of the Higgs fully-inclusive cross section. The condition m h < 2m t implies that the top degrees of freedom can be safely integrated out. By doing so, the three-point one-loop ggH function is identified with an elementary vertex (and likewise for the four-and five-point functions), whose corresponding interaction Lagrangian is added to that of the SM (minus the top) to define an effective field theory (EFT). From a more physical viewpoint, this is equivalent to saying that the virtuality of the gg system, equal to the Higgs mass, is not large enough to resolve the top-quark loop, whose kinematical complications can therefore be ignored. In this way, at any given order (in perturbative QCD) EFT computations feature one loop less than those in the full theory: for example, the Born gg → H process is associated with a treelevel (as opposed to one-loop) diagram. EFT results at the NLO in QCD have proved (somehow heuristically) to be an excellent approximation of those obtained in the SM [3][4][5][6], once top-mass effects are accounted for at the LO. This has given a strong motivation JHEP08(2016)006 to pursue calculations at yet higher orders, and predictions have thus been obtained at the NNLO [7][8][9] and even at the N 3 LO [10].
While such theoretical results are extremely significant, it is important to bear in mind that a deeper understanding of the m t → ∞ regime would be desirable; moreover, the knowledge of differential distributions is especially crucial. Firstly, these must be used in the computation of the acceptance corrections that allow one to relate the visible and the total cross sections. Secondly, they can be employed either in the direct comparison to measured differential data (see e.g. refs. [11][12][13][14][15][16] for recent LHC results), which will become more numerous and precise in the near future, or in the context of multivariate analyses. Thirdly, and connected to the previous item, experimental data may have to be characterised in terms of fixed jet multiplicities, and theory predictions must match such a characterisation. At the differential level, in those phase-space regions where a Higgs and/or jets with sizable energy can be found (the rule of thumb being e.g. p T (H) m H ), the EFT is an increasingly poor approximation, since the virtuality of the incoming-parton pair becomes sufficiently large to resolve even top-quark loops. It is therefore mandatory to obtain predictions for all the relevant observables by retaining the exact heavy-quark mass dependences in the loops; in this way, one is also able to treat bottom-quark contributions in a sensible manner.
In the course of the past few years, full-SM predictions that lift the EFT approximation have started to appear. For inclusive cross sections, the NNLO results of refs. [17][18][19] (obtained by means of an asymptotic expansion in 1/m t ) have confirmed earlier findings at the NLO, that top-mass effects have a very small impact. Differential observables are much more laborious to compute, and only a handful of results exist that go beyond the m t → ∞ approximation (conversely, there has been a significant activity within the EFT, too vast to be covered here; for recent reviews, see e.g. refs. [20][21][22]). At fixed NLO, subleading 1/m t predictions are available for the Higgs and hardest-jet p T [23,24]. Mass effects are also accounted for in analytically-resummed transverse momentum spectra [25][26][27]. As far as perturbative computations matched to parton showers are concerned, multi-jet merged predictions at the leading order have been first presented in ref. [28] (see also ref. [29]). NLO+PS simulations have also been available for a while for inclusive H + 0j production, in the POWHEG [30] and MC@NLO [31] (see also ref. [32]) frameworks, with the latter subsequently employed [29] within the MEPS@NLO [33] merging scheme. Finally, NLOaccurate mass effects have been included in an NNLOPS implementation [34].
The aim of this work is to assess the impact of heavy-quark masses on observables of current experimental interest, as well as to give realistic, hadron-level, differential predictions for both the Higgs boson and several accompanying jets. This is achieved in the context of merged simulations, where MC@NLO [35] samples of a given H + nj multiplicity are consistently combined by means of the FxFx method [36] for n = 0, 1, 2; the Pythia8 [37] parton shower is used throughout. Born and real-emission amplitudes are computed exactly at one loop for all samples (therefore, the one-loop amplitudes with the highest mutiplicity considered feature three final-state partons, and enter the real corrections of the n = 2 sample). Conversely, the virtual amplitudes are computed exactly at two loops [6,38,39] for n = 0; when n = 1, 2, the one-loop EFT virtuals are used in-JHEP08(2016)006 stead, rescaled by the ratio of the full over the EFT Born amplitude squared. While we note that our approach thus employs all of the most accurate matrix elements available to date for n ≤ 2, it would be interesting to re-assess our findings should the exact two-loop amplitudes with n = 1, 2 become available. 1 In summary, the present work has a scope similar to that of refs. [28,29]; it differs from the former in that it features NLO effects, and from the latter in the merging scheme employed, as well as in the facts that it includes (some of the) two-loop and bottom-quark contributions, and that two-jet observables are also NLO-accurate.
The paper is organised as follows. In section 2 we describe the nature of the simulations we have performed, as well as their detailed setups, including input parameters. In section 3 we present predictions for inclusive rates and differential observables. We draw our conclusions in section 4.

Details on the calculations
We consider four different types of simulations, all accurate to NLO 2 in QCD and matched to parton showers, which we list here roughly in order of increasing complexity: • inclusive H + 0j in the EFT, 3 denoted by inc EFT ; • inclusive H + 0j in the full SM, denoted by inc M ; • merged H + nj (n ≤ 2) in the EFT, denoted by FxFx EFT ; • merged H + nj (n ≤ 2) in the full SM, denoted by FxFx M .
It is by now common practice [26] in the context of matrix-element computations combined with resummation (be it analytical or achieved through parton showers) to treat separately the contributions due to squared top-loop amplitudes from those stemming from bottom loops (which thus include the interference between top-loop and bottom-loop amplitudes). This allows one to assign different resummation or shower scales to such contributions, with the scale relevant to the latter (of the order of a few bottom masses) smaller than that relevant to the former (of the order of m H ). One needs to keep in mind that this procedure is pragmatic, but does not address directly the three-scale problem posed by the bottom-loop contributions, a problem that must still be considered as open (for recent considerations, see e.g. ref. [41]). The presence of tagged final-state jets further complicates the picture, because it potentially introduces extra hard scales. This implies that the merging of samples stemming from bottom loops is necessarily quite heuristic, given the 1 Such amplitudes might induce some differences, typically in the kinematics regions characterised by large mass scales (e.g. above the top-quark threshold), w.r.t. the approximated ones we have used, if effects that do not factorise the Born would turn out to be important, as was recently observed [40] in a case that bears some analogy with the present one, namely di-Higgs hadroproduction.
2 Strictly speaking, in order to formally attain such an accuracy in the full SM we would need to compute the virtual amplitudes at two loops, which we do only in the case of the lowest multiplicity. 3 These are pure EFT results: there is no rescaling (e.g. overall) by the full-SM over EFT Born matrix elements.

JHEP08(2016)006
delicate role that mass scales play in merging procedures. In view of this, and of the fact that (even with a more sophisticated treatment of the three-scale problem w.r.t. what is done at present) the resummation or shower scale naturally associated with bottom-loop contributions will result in relatively soft jets, we have adopted the simple solution of treating these contributions purely as an inclusive H + 0j sample. 4 Such a sample is then added to the samples stemming from top loops: to the inclusive one to give inc M , and to the FxFx-merged one to give FxFx M . By adopting a rather sketchy notation, we can summarise the previous discussion as follows: with each sub-sample (in round brackets) meant to be NLO-accurate, and the FxFx "operator" on the r.h.s. of eq. (2.2) understood to implement the FxFx merging. In a few occasions we shall also present full-SM results by retaining only top-loop contributions. By using the notation of eqs. (2.1) and (2.2), these are therefore: Our calculations are performed by means of MadGraph5 aMC@NLO [42], which automates all ingredients relevant to the simulation of LO and NLO cross sections, including the matching to parton showers. The FKS method [43,44] (automated in the module MadFKS [45]) for the subtraction of the infrared (IR) singularities of real-emission matrix elements underpins all NLO-accurate results. The computations of one-loop amplitudes are carried out by switching dynamically between two integral-reduction techniques, OPP [46] or Laurent-series expansion [47], and TIR [48][49][50]. These have been automated in the module MadLoop [51], which in turn exploits CutTools [52] or Ninja [53,54], together with an in-house implementation of the OpenLoops optimisation [55]; we point out that, starting from v2.3.0, MadLoop can handle loop-induced processes automatically [56].
The status of MadGraph5 aMC@NLO just described allows one to perform all simulations outlined before without any human intervention in the EFT case, and with a minimal amount of it when working in the full SM. However, the latter is very time consuming, owing to the presence of numerically-computed loop amplitudes in all the contributions to, and the stages of, the computations, which renders the generation of millions of events rather impractical on a medium-size cluster. We have therefore proceeded as follows: • All of the matrix elements relevant to the (H + 0j) α sub-samples, with α = t 2 , b 2 , bt, have been replaced, by means of the aMCSusHi script [32], with their SusHi [57] analytically-computed counterparts. Note, in particular, that this implies that twoloop virtuals are employed in such sub-samples, which are integrated and unweighted directly.

JHEP08(2016)006
• The one-loop matrix elements for the other sub-samples, (H + 1j) t 2 and (H + 2j) t 2 , are stored in a library generated by MadLoop. The library includes a wrapper that allows one to call the individual matrix elements.
• The sub-samples of the previous item are first integrated and unweighted by working in the EFT. The resulting hard events are then reweighted by using the relevant 5 full-SM over EFT ratios, with the numerators computed on the fly using the MadLoop library constructed before. This is the same procedure used to compute di-Higgs production in gluon fusion in refs. [58,59].
A few comments are in order. Firstly, the reweighting procedure described above does not entail any approximation. Its only implication is that the final event samples, given in input to the MC, are weighted rather than unweighted (as it would be customary in Mad-Graph5 aMC@NLO). Although a secondary unweighting could be envisaged, we did not consider it in this work. Secondly, the Higgs-plus-three-partons top-loop amplitudes enter (although in regions that do not overlap kinematically) both the (H + 0j) t 2 contribution (as real corrections) where they are computed by SusHi, and the (H + 1j) t 2 contribution (as Borns) where they are computed by MadLoop. Given the strict equivalence between SusHi and MadLoop one-loop results, the (H + 1j) t 2 bit could also be computed with SusHi; this would be faster, but rather more awkward from a procedural point of view, and we have therefore refrained from doing it. Thirdly, the reweighting by the full-SM matrix elements adopted in this paper is carried out by using, in part, software routines which are not public. However, the hard events will be made publicly available. 6 Furthermore, future public versions of MadGraph5 aMC@NLO will contain NLO-accurate reweighting capabilities, equivalent to the procedure adopted here. The matrix elements have been computed by using the input parameters 7 reported in table 1. The central renormalisation and factorisation scales (µ 0 ) are set in merged simulations as dictated by the FxFx prescription [36], and for inclusive results equal to H T /2. In both cases, uncertainties are computed by varying independently the two scales in the range [µ 0 /2, 2µ 0 ]; this is done automatically by MadGraph5 aMC@NLO and does not entail re-running, thanks to the method of ref. [61]. The merging-scale dependence has been assessed by comparing the results obtained by setting µ Q = 20, 30, and 50 GeV, with 30 GeV being our default choice that corresponds to the central predictions to be shown in section 3. We have also considered the very extreme choice µ Q = 70 GeV, but since the corresponding results differ only marginally from those obtained with µ Q = 50 GeV we refrain from reporting them here.
In Pythia8, of which we have used both v8.210 and v8.212 (with identical results), the parameters have been set equal to their defaults, possibly as defined by the Mad-Graph5 aMC@NLO interface. The Higgs boson is kept stable, and underlying events are 5 When very close to the IR limits of the real-emission matrix elements, Born matrix elements are employed in the reweighting, in order to ensure a more stable numerical behaviour, and a marginally faster procedure. Furthermore, as is customary Borns are also used to reweight the virtuals. 6 The interested reader is advised to contact the authors.  not generated. We point out that all the FxFx-related operations are now fully automated in Pythia8; more details can be found in ref. [62], where a phenomenology validation of the FxFx method is performed; at variance with ref. [62], no event oversampling is carried out in the present paper. The jets that enter the definitions of physical quantities are reconstructed at the hadron level by means of the anti-k T algorithm [63] with R = 0.4 (as implemented in FastJet [64]), by keeping those with: We finally recall that MadGraph5 aMC@NLO gives the user some control on the eventby-event shower scale Q sh to be passed to the MC in the case of inclusive simulations, choosing it at random in the range where √ s 0 is the Born-level partonic c.m. energy, and α, f 1 , and f 2 are user-defined numerical constants -more details on this can be found e.g. in section 2.4.4 of ref. [42] and in section 3.2 of ref. [65]. This possibility is exploited in the present context to assign different shower scales to the three sub-samples that appear on the r.h.s. of eq. (2.1) (and thus also to the two rightmost terms on the r.h.s. of eq. (2.2)), in keeping with the procedure adopted in refs. [32,[66][67][68]. In particular, by choosing the settings suggested in ref. [67], the peak values of the shower-scale distributions for the three sub-samples are Q t 2 sh ∼ m H /2, Q b 2 sh ∼ 21 GeV, and Q bt sh ∼ 31 GeV.

Inclusive rates
We start by presenting inclusive rates, which we collect in table 2. Each of the four leftmost columns presents the results obtained with one of our standard merged or inclusive simulations. In the rightmost column we report the contributions of the bottom-loop induced processes to the full-SM merged and inclusive predictions. In other words, for any given row the entries in the FxFx M and inc M columns include the entry in the σ b column; in turn, σ b is the cross section that emerges from the sub-samples (H +0j) b 2 and (H +0j) btsee eqs. (2.1) and (2.2). All results are obtained with central hard and merging scales; by using an error notation, we also show the fractional hard-and merging-scale systematics, obtained as discussed in section 2.
The first row of table 2 features the results obtained without imposing any final-state cuts. The predictions in the next three rows are relevant to requiring a given number of JHEP08(2016)006  exclusive (N jet = 0, 1) or inclusive (N jet ≥ 2) jets among those that satisfy the cuts in eq. (2.5). The two rows at the bottom present cross sections in two phase-space regions defined by VBF-like cuts: where j 1 and j 2 denote the hardest and second-hardest jet of the event, respectively. We observe that the hard-scale uncertainties of the merged results are rather similar to those of the inclusive ones, if only slightly larger, consistently with the fact that they receive contributions from matrix elements that feature higher powers of α S w.r.t. those that enter the inclusive samples. Such uncertainties are significantly larger, for all of the acceptance regions considered, than the merging-scale ones, in spite of the fact that the latter are associated with a very conservative choice for the range of µ Q . This is particularly striking in the case of the two VBF-like regions, for which the small merging-scale dependence implies that the descriptions, given by the matrix elements and by Pythia8, of kinematic configurations where jets tend to be close to the beam line are mutually quite consistent with each other. This is clearly an MC-dependent statement, and is e.g. at variance with what has been found in the past with Herwig6 both at the LO and the NLO [69]: it is wise to always assess the merging-scale systematics associated with a given parton shower, especially in longitudinally-dominated regions.
The impact of the bottom-induced contributions, measured by the ratio of σ b over FxFx M or inc M , is non-negligible in the regions dominated by small jet multiplicities, but irrelevant elsewhere. This is a consequence of the choices made for the treatment of such contributions (see section 2), and consistent with the idea that bottom-loop amplitudes should describe the emission of mostly soft jets. As expected, the opposite behaviour is associated with merging, whose impact becomes more significant the larger the jet multiplicity. Even in the fully-inclusive case (first row of table 2), merged rates are larger by JHEP08(2016)006 about 5.5% w.r.t. the inclusive ones. The cross section thus increases (we recall that FxFx is a non-unitary prescription), and this increase, as we shall see from the differential results of section 3.2, is essentially driven by the contributions of the sub-samples (H + 1j) t 2 and (H + 2j) t 2 .
In order to assess in a more transparent manner the effects of merging, we report in table 3 the ratios of the cross sections for a given number of inclusive or exclusive jets, over the corresponding fully-inclusive cross sections. Note that these results are not independent from those of table 2, but are obtained from the latter (by using only the central values). For example, the first, second, and third entry from the left in the first row of table 3 are evaluated by computing the ratio of the entry in the second, third, and fourth row, respectively, of the first column of table 2, over the first entry in the same column of that table. As can be seen from table 3, the merging decreases the N jet = 0 cross sections by a relative 10%, while it increases the N jet = 1 and N jet ≥ 2 ones by a relative 15% and 25%, respectively. These fractions are essentially independent of whether heavy-quark mass effects are included or not. However, masses do have an impact for a given class of calculations (merged or inclusive); in particular, they decrease the N jet = 0 rates (owing to bottom-mass effects, as we shall discuss below), while they increase the N jet = 1 and N jet ≥ 2 ones.
A different way of looking at mass effects stems from the computation of the ratios of full-SM over EFT cross sections, whose results we show in table 4 for both the merged and the inclusive simulations. In the two rightmost columns of the table, the SM predictions are computed by excluding bottom-loop contributions (i.e. the numerators correspond to eqs. (2.3) and (2.4)). The most striking feature of table 4 is that such effects factorise almost exactly in the case of all the rates except the VBF ones; in other words, they are independent of whether merging is carried out or not, except for longitudinally-dominated configurations. In absolute value, mass effects are rather small when both top-and bottomloops are taken into account: almost completely negligible for fully-inclusive cross sections, 8 and at most ±3% for the other rates considered here. The effects are larger, for the fullyinclusive and lowest-jet multiplicity cases, when only massive top quarks are included in the computations.

Differential distributions
We now turn to presenting differential observables, constructed with the Higgs and the jets four-momenta. We begin by considering the transverse momentum of the Higgs boson; no cuts are imposed. We display this quantity in figure 1, whose layout is the following. The conclusions that can be drawn from figure 1 are the following. The effects of multi-jet merging start to be quite significant at p T (H) 50 GeV; for p T (H) m H , the inclusive results are outside of the hard-scale uncertainty band of FxFx M . However, with increasing transerve momentum inc EFT falls less rapidly than FxFx M , and as a consequence of that at about p T (H) ∼ 350 GeV it is again nearly compatible with the lower border of the uncertainty band of the latter. This is ultimately due to the fact that top-quark mass effects start to be visible for p T (H) 250 GeV, where they suppress the full-SM results w.r.t. their EFT counterparts. As can be seen from the middle inset, by comparing the histograms with the symbols, heavy-quark mass effects almost exactly factorise w.r.t. the merging procedure: they affect equally the merged and the inclusive predictions, which is quite consistent with what has been already observed for inclusive rates in section 3.1. We note that this applies both to the large-and to the small-p T (H) region. In the latter, for p T (H) 50 GeV, the bottom-loop contributions do have a non-negligible impact on the shape of the distribution, in keeping with what previously found [25,26,30]. Finally, the theoretical systematics that affect the FxFx M result also have a similar pattern as those relevant to inclusive rates: namely, on the whole transverse-momentum range considered, hard-scale uncertainties largely dominate over merging-scale ones. The latter are in fact at most ±10%, and typically much smaller than that. In summary, multi-jet merging is more relevant to obtaining a sensible prediction for the p T (H) spectrum than the exact JHEP08(2016)006 treatment of heavy-quark loops when statistics is limited, because this restricts one to relatively small transverse momenta and/or to use bins of large widths. However, as soon as one is able to access the high-p T region and to better resolve the details of the spectrum (including its low-p T end), then mass effects must mandatorily be taken into account. The Higgs transverse momentum can also be observed in a more differential manner, by requiring a given number of accompanying jets; this is potentially very relevant to experimental analyses which employ N jet categorisation. We present our predictions for this quantity in figure 2, by requiring N jet = 0, N jet = 1, and N jet ≥ 2. The layout of the figure is the following. In the main frame, the solid histograms show the FxFx M results for the three N jet cases, together with the respective hard-scale uncertainty bands. Each of the JHEP08(2016)006 three insets, which have an indentical layout, is relevant to one N jet case. All results shown there are ratios of different numerators over the same denominator, that coincides with the central FxFx M prediction. Dot-dashed, dashed, and dotted histograms are obtained by setting the numerator equal to FxFx EFT , inc M , and inc EFT , respectively. Open circles and crosses correspond to the µ Q = 20 GeV and µ Q = 50 GeV FxFx M cross sections, and thus give the fractional uncertainty associated with merging-scale choice. Finally, the shaded area is the fractional hard-scale systematics.

JHEP08(2016)006
Our findings can be summarised as follows. For p T (H) 40 GeV, the N jet = 0 result becomes smaller than the N jet = 1 one, and falls very rapidly. At p T (H) ∼ 100 GeV, the N jet = 1 and N jet ≥ 2 predictions are comparable, with the latter becoming dominant as the JHEP08(2016)006 transverse momentum increases; the importance of a proper multi-jet merging is therefore obvious. All types of simulation give very similar results in the N jet = 0 bin; in particular, they are all compatible with each other within the FxFx M hard-scale uncertainty (which largely dominates over the merging-scale one). Things do change when requiring non-null jet multiplicities. When N jet = 1, the impact of merging is visible for p T (H) 50 GeV, and that of heavy-quark masses for p T (H) 250 GeV. Although these values are roughly the same as those relevant to the inclusive Higgs p T , the pattern of the corresponding effects is different; note, in particular, that the inc EFT prediction is significantly harder than in the inclusive case. However, this feature should not be over-interpreted; in the presence of a tagged jet, inc EFT will depend to a certain extent on the MC adopted for the parton-shower phase. As for N jet = 0, hard-scale systematics dominates. The N jet ≥ 2 results follow roughly their N jet = 1 counterparts, with two notable exceptions. Firstly, merging effects are present also at small p T (H) (mostly owing to configurations where two jets harder than the Higgs are present in the transverse plane, and recoil against each other). Secondly, in that region the merging-scale dependence is not completely negligible. We see, in particular, that the µ Q = 50 GeV result tends to be closer to the inclusive ones, which is consistent with the fact that such merging-scale choice emphasises more the role of the parton shower over that of the matrix elements, thus diminishing the probability of finding two jets harder than the Higgs. Finally, we note that fairly small mass effects are also present at low p T (H) (although difficult to see in figure 2), especially in the N jet = 0 and N jet ≥ 2 bins, where they are predominantly induced by bottom-and top-loop contributions, respectively.
We next consider the transverse momenta of the four hardest jets, which we display in figure 3; note that this figure has the same layout as figure 2. These observables are correlated, to different extents, to the inclusive or jet-binned Higgs transverse momentum previously discussed, and they thus have some characteristics whose patterns are similar to those of the latter. In particular, the onset of merging or of heavy-quark mass effects occurs at transverse momenta of the same order as those relevant to the Higgs spectra. On the other hand, the impact of merging is generally larger than for the Higgs p T (except perhaps in the case of the hardest jet), and it dramatically increases when one considers jets that are increasingly subleading. This applies also to jets which are formally beyond the matrix-element accuracy of the simulations performed in this work, and it need not be surprising. The same feature has been discussed e.g. in ref. [62] in the context of Z/W+jets production (which is analogous to Higgs production as far as merging is concerned), and validated against data there. Essentially, the benefits of merging extend beyond those naively expected from the underlying matrix elements, because the latter provide the parton showers with much more sensible initial conditions than those available in inclusive simulations. In support of this fact, we also observe that the merging-scale dependence is again much smaller than the hard-scale one, and all jets feature roughly the same stability against µ Q variations (the oscillations in some of the p T tails are clearly not resulting from a pattern associated with µ Q , but merely reflect a lack of statistics).
We now move to consider a few correlations between the two hardest jets of the event, which we present in figures 4-7. All these figures have exactly the same layout as figure 1, to which we refer the reader for detailed explanations. We begin in figure 4 with the azimuthal JHEP08(2016)006 distance; this observable has a distinctive shape, and can be used e.g. to help assess the presence of a BSM signal in data. As is clear from the figure, the impact of merging is quite significant, especially shape-wise. Conversely, heavy-quark mass effects are barely visible in the central predictions, and amply within the hard-scale systematics. This has to be expected, since the present observable is dominated by configurations in which jets are produced at their p T thresholds. As in all previous cases, merging-scale uncertainties are much smaller than their hard-scale counterparts. Having said that, interestingly there is a trend at ∆φ(j 1 , j 2 ) → π, where the µ Q = 50 GeV result is less peaked than the others. In other words, for such a merging scale the two hardest jets are less back-to-back than JHEP08(2016)006 for smaller merging scales; this is the same phenomenon responsible for the small-p T (H) behaviour when N jet ≥ 2, previously discussed in figure 2. The invariant mass of the pair of the two hardest jets and their rapidity distance are shown in figure 5 and figure 6, respectively. These two observables are important in that they are used to define VBF regions; therefore, their shapes control the extent to which a possible VBF signal is "contaminated" by that due to gluon fusion. Apart from the obvious, observable-specific differences, what has been said about ∆φ(j 1 , j 2 ) applies to M (j 1 , j 2 ) and ∆y(j 1 , j 2 ) as well. In particular, the effects of merging are evident, and especially prominent in the case of the rapidity correlation. They harden the spectrum of the invariant mass distribution, most notably at intermediate masses (∼ 150 − 400 GeV), and decrease the rapidity separation between the two jets, so that the merged results have relatively more events than the inclusive ones towards small ∆y(j 1 , j 2 ). Heavy-quark mass effects are only visible in the invariant-mass spectrum (see the middle inset of figure 5), where size-wise they are roughly as important as the merging ones.  a softening of the distribution w.r.t. the EFT behaviour, which softening in the large-mass tail seems to affect only the merged predictions. While lack of statistics may be an issue in this region, this trend is related to the breakdown of the factorisation of mass effects remarked in section 3.1 for the case of VBF-like inclusive cross sections. This is most likely due to an increased difference between the high-multiplicity matrix elements computed with or without mass effects, since in this kinematic region multiple large scales become relevant in the computation of the loop amplitudes. Finally, also for the two observables discussed here the merging-scale systematics is much smaller than that due to hard-scale variations. We remind the reader that, in the case of ∆y(j 1 , j 2 ), this is a statement that depends rather strongly on the MC adopted; see section 3.1 for related observations on VBF inclusive rates. We finally present in figure 7 the azimuthal distance between the two hardest jets within VBF 1 cuts. We find patterns which are similar to those observed in the case of its inclusive companion of   Figure 6. Same as in figure 1, for the rapidity distance between the two hardest jets.
the two figures, one sees that in the present case the effects of merging are much larger in absolute value. The almost complete de-correlation between the two jets resulting from the inclusive simulations demands that merging be carried out, in order to achieve physicallysensible results. Heavy-quark mass effects are unimportant, and hard-scale uncertainty is much larger than the merging-scale one, in spite of the latter being more significant than in all the other cases considered so far.

Conclusions
In this paper we have studied the production of a Standard Model Higgs boson in association with jets through the gluon-fusion channel, presenting results for both inclusive rates and differential observables relevant to the 13 TeV LHC. We have done so by considering several types of simulations, in all of which matrix elements are matched to parton  Figure 7. Same as in figure 1, for the azimuthal distance between the two hardest jets, within the VBF 1 acceptance region.
showers at the NLO accuracy in QCD according to the MC@NLO formalism. We have systematically compared predictions that stem from inclusive samples with those based on the consistent merging, by means of the FxFx method, of sub-samples characterised by different parton-level multiplicities, for up to two extra jets at the NLO. Within each type of approach, inclusive or merged, we have evaluated the underlying matrix elements both in the EFT where the Higgs couples directly to gluons, and in the SM by computing exactly the relevant top-and bottom-loop amplitudes (except for the two-loop virtual contributions to the one-and two-jet cross sections, which we have approximated). All of our computations are carried out in the MadGraph5 aMC@NLO framework, and we have adopted the Pythia8 parton shower.

JHEP08(2016)006
The primary conclusion of this work is that both merging and heavy-quark-mass effects must be taken into proper account in order to predict realistically observables of experimental interest, that will further our understanding of the particle discovered by ATLAS and CMS. The impact of merging is especially prominent, while in order to be sensitive to quark masses it is important to be able to probe kinematic regions dominated by large-p T emissions, and to have a good resolution power to allow for a sufficiently fine binning, which helps uncover such effects at small transverse momenta as well. Both features are characteristic of datasets with substantial statistics.
For a sufficient number of events, merging and mass effects provide one with very distinctive features, that have the further advantage of being fairly observable-dependent, thus allowing for stringent comparisons between predictions and measurements. The theoretical systematics are dominated by the dependences on the factorisation and renormalisation scales. Conversely, variations of the FxFx merging scale, which we have carried out in a very large range in order to be conservative, are almost always negligible, with a very few exceptions in which they are still smaller than those induced by the hard scales. However, we point out that this is an observable-and MC-dependent conclusion; the assessement of the merging-scale systematics is an important task that must always be performed.