Higgs boson pair production merged to one jet

We develop a Monte Carlo event generator for Higgs Boson pair production merged to exact one-jet matrix elements. The matrix elements are generated with OpenLoops and event generation is performed with the HERWIG++ general-purpose event generator. This allows us to simulate fully-exclusive hadronic final states with accurate description of the kinematics of the leading jet in conjunction with a parton shower. We use the implementation to examine in detail the systematic uncertainties which result from the merging procedure. We assess the magnitude of the impact of the merging on experimental searches of Standard Model di-Higgs production that aim to constrain the Higgs boson self-coupling. We find that the use of a merged sample can reduce theoretical systematic uncertainties in the efficiencies of cuts on certain observables. This constitutes the most accurate simulation of the process available to date. The Monte Carlo event generator developed for this project is available as an add-on to the HERWIG++ event generator at http://www.physik.uzh.ch/data/openloops/download/projects/hhmerge/.


Context and introductory remarks
The Large Hadron Collider ATLAS and CMS experiments confirmed the existence of a scalar particle consistent with the Standard Model (SM) Higgs Boson [1][2][3][4][5]. However the quest for understanding the mechanism behind the electroweak symmetry breaking (EWSB) does not end with the discovery of the Higgs Boson; measuring its couplings to the SM fields is an important and long-term task that the ATLAS and CMS experiments, as well as future collider experiments, are expected to undertake. It is crucial, moreover, to determine whether the realisation of the mechanism of the EWSB is indeed SM-like. This can be investigated by examining the Higgs potential which, after EWSB in the minimal prescription, can be written as Within the SM we have λ SM hhh = λ SM hhhh = m 2 h /(2v 2 ) 0.13 for a Higgs boson mass of m h 125 GeV. The discovery of the Higgs boson only indicates the size of the curvature of the potential around the local minimum, coming from the quadratic term. To confirm the form of the potential, the measurement of higher-order terms is necessary. At the LHC, these terms can be probed directly via double or triple Higgs Boson production. The tiny cross section for triple Higgs production makes it impossible to perform any meaningful measurement in the foreseeable future, even during the full life-time of the LHC [6,7]. Higgs Boson pair production, on the other hand, is certainly challenging but not impossible to observe at the LHC. Interesting phenomenological studies were performed more than 10 years ago [8][9][10][11][12][13] and more recently, owing to the discovery of the Higgs Boson as well as the development of boosted jet techniques, the subject has undergone a lively rejuvenation [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28].
Despite the fact that several interesting and in-depth phenomenological studies of inclusive Higgs Boson pair production at the LHC (pp → hh + X) have been performed, the Monte Carlo event simulation of the process has relied so far only on leading-order matrix elements with the addition of parton showers to simulate the extra QCD radiation. 1 Exceptions to this are two recent studies which examined the exclusive one-and two-jet channels in the full theory, with the full top mass dependence, (i.e. pp → hhj + X and pp → hhjj + X) and contrasted these to results obtained in the effective theory [16,30]. It is important to stress, however, that the kinematical properties of inclusive final states can be substantially altered by the inclusion of higher-order matrix elements. This is especially true in the inclusive hh+X process, which is predominantly gluon-gluon initiated, and hence is inevitably accompanied by a copious amount of QCD radiation. Thus, the accuracy, and hence reliability, of the kinematics of inclusive di-Higgs searches will certainly benefit from the inclusion of the exact real-emission higher-order matrix elements.

Di-Higgs production at higher orders
The pp → hh + X process at hadron colliders is loop-induced at leading order, proceeding via a heavy quark loop. The leading-order gluon fusion diagrams are shown in figure 1. Evidently, a next-to-leading order calculation would involve, among others, diagrams with two loops that involve heavy fermions, and hence two mass scales (the fermion mass and the Higgs Boson mass). Such diagrams currently lie at the frontier of higher-order loop calculations. Consequently, this impedes the implementation of a matched next-to-leading order (NLO) plus shower simulation.
The effective theory approximation in the heavy top mass limit that has been employed in single Higgs boson production has been shown to be insufficient to describe the kinematics of the hh process, both at leading order [22,31], and at higher orders [16,30]. Nevertheless, inclusive NLO [9,32] and NNLO cross section calculations [20,21] have JHEP03(2014)126 been performed in the effective theory approximation, giving an estimate of the size of the higher-order corrections in the full theory.
Thus, in the absence of the full NLO calculation, the best one can do is merge samples of different jet multiplicities in a consistent way, carefully avoiding any issues that may arise due to double-counting or phase space region mismatch. Such simulations have been shown to reliably describe the kinematical properties of experimental data (see for example the relevant experimental CMS [33] and ATLAS [34] analyses), modulo the correct normalisation taken from higher-order cross section calculations. Here we perform such a merging of samples, including the full top and bottom mass dependence in the fermion loops of figure 1, as well as the higher-order real emission diagrams which we examine below.
This paper is organised in the following way: in section 2 we briefly describe the OpenLoops generator for one-loop matrix elements and provide cross sections for the various contributing exclusive channels. In section 3 we present results and examine the systematic uncertainties associated with the merging prescription, and in section 4 we investigate the phenomenological implications of including the merged higher-order matrix elements. We present our conclusions in section 5.

The OpenLoops matrix element generator
The OpenLoops generator is based on the open-loops algorithm [35] for the efficient evaluation of one-loop matrix elements. The algorithm employs a numerical recursion to construct the loop momentum dependence of the numerator of loop amplitudes combined with tensor integral reduction. The tensor integrals are computed by the Collier library, which implements the Denner-Dittmaier reduction procedure for the numerically stable evaluation of tensor integrals [36,37] and the scalar integrals of ref. [38]. Incidentally, using tensor integrals allows for a high degree of optimisation through caching, since the integrals can be shared across different Feynman diagrams for all helicity and colour configurations. For on-shell reduction approaches this is only possible when the loop amplitude is interfered with a tree amplitude [35] and therefore not in calculations of loop induced processes, like the one presented here.

Higgs pair production matrix elements
Like in the case of single Higgs boson production, the hh production cross section at hadron colliders is dominated by the gluon fusion channels. For a more detailed dissection of the leading order cross section, we refer the reader to section 2 of ref. [18].
The classes of higher-order real emission diagrams that we include in our calculation are shown in figure 2. The real emission process also contains diagrams with qg,qg and qq initial states which are subdominant but non-negligible and must be included for consistent merging since the parton shower introduces g → qq splittings on the initial state gluons of the 0-jet matrix elements [39]. It is important to note that the diagrams that involve radiation from the heavy quark loop, are not included in any limit in the parton shower Monte Carlos with which we merge with. Hence, an intrinsic assumption of the merging procedure is that these diagrams are sub-dominant with respect to the initial state radiation in the parton shower-dominated regime. The OpenLoops process libraries to compute matrix elements have been interfaced with HERWIG++. These can be used stand-alone, i.e. without the merging, to perform studies of leading-order hh production, or hh + j production. In table 1 we present the cross sections for the different sub-processes contributing to pp → hhj + X, where j is an associated parton. 2 Here, and throughout this paper, we use the 4-flavour MSTW2008nlo 68% confidence level parton density functions [41][42][43]. Obviously, even with the relatively high p ⊥ cut of 60 GeV, the real emission sub-processes possess a cross section that is comparable to the leading order gluon fusion process. This is an indication that they are indeed significant and have to be considered for an accurate description of the kinematics, even in an inclusive hh + X analysis.

Merging methods
In order to obtain a realistic simulation of processes involving associated high-p T jet production, e.g. W/Z/Higgs+jets, the parton shower approximation for the generation of soft and collinear QCD radiation must be supplemented by high multiplicity leading-order matrix elements. Matrix element-parton shower merging schemes, such as the so-called MLM [44][45][46] and CKKW [44,45,[47][48][49][50] methods, have been developed for this purpose. These methods work by partitioning phase space, by means of a jet algorithm, such that 2 The HERWIG++ implementation has been crossed-check against the SHERPA event generator [40].   (7) 60.6(9) 27.1(4) 0.79(1) Table 1. Cross sections for the partonic pp → hh + X and for the sub-processes contributing to pp → hhj + X at 14, 33 and 100 TeV. For the case of real emission, a cut of p ⊥ > 60 GeV was placed on the associated parton. The factorisation/renormalisation scales were both fixed to µ = m h + p ⊥,j , where p ⊥,j is the transverse momentum of the associated parton in the centre of mass frame.
the distribution of jets corresponds to that of the partons in the matrix elements, while the distribution of radiation inside the jets is appropriately developed by the shower. In addition, both the MLM and CKKW algorithms augment the distribution of radiation in the matrix element region with Sudakov suppression effects, not present in the matrix elements themselves, thus smoothing the transition from one radiation pattern to another at the phase space partition. 34 [55][56][57][58] includes an implementation of the MLM merging scheme. The current version of the merging algorithm has been validated against its FORTRAN HERWIG [59] counterpart for several processes. For the purposes of this project, the implementation has undergone minor modifications, to accommodate the use of internally-generated matrix elements. We use this algorithm in conjunction with the parton shower in order to merge the two. We fix the factorisation and renormalization scales to be equal, , where ν is a parameter which we vary, m h and p hh ⊥ are the Higgs boson mass and the transverse momentum (as defined in the centre-of-mass frame of the hard process) of the Higgs boson pair respectively. Note that for the LO hh process, p hh ⊥ = 0 and hence this implies that µ = νm h for all, even showered, LO samples. We call the merging scale E Tclus , inspired by the way the MLM method is implemented in the HERWIG++ generator. We call the lowest-order sample '0-jet' and the sample including one real emission '1-jet'. Broadly speaking, after showering is performed in HERWIG++, the MLM method will effectively veto all events in the '0-jet' sample that contain a jet with transverse momentum larger than E Tclus . This will result in what we will call the '0-jet exclusive' sample. In the showered '1-jet' sample the MLM algorithm will effectively veto any events with jets that have not 'matched' 5 the given extra parton produced in association with the Higgs boson pair, as well as events that contain jets harder than the 'matched' jet. The resulting sample is called '1-jet inclusive', meaning it contains no 0-jet contributions JHEP03(2014)126 but contains jets coming from the shower, of lower p ⊥ than the matrix element parton. For more details on the algorithm see, for example, [46].
Recently it has been shown [60] that the merging of samples of different multiplicities may result in kinks due to the presence of a significant mismatch in the description of extra emissions between the parton shower and the matrix element calculations in the region chosen for merging. In ref. [60], it was suggested that one can obtain smooth matching by the use of a smoothing 'D-function' that contains two scales instead of a single scale. Sudakov reweighting was also used to achieve an even smoother matching. Here we perform a variant of the former method: we generate a merging scale randomly in a given interval according to a given distribution. This merging 'range' is then characterised by two scales, an 'average' merging scaleĒ Tclus and a 'variation' scale clus . The merging scale is then randomly chosen on an event-by-event basis using two different 'schemes'. The first scheme uses a Sine function ('Sinusoidal'), and the second a linear function ('Uniform'), where x ∈ [0, 1] is a uniform random number. The effect of these schemes is to smooth out the unphysical discontinuities, resulting in a continuous merging of the shower and matrix element descriptions.
To further improve the merging between the two samples, we perform 'α S -reweighting' of the 1-jet matrix elements according to (schematically) where p hh ⊥ is the transverse momentum of the Higgs boson pair as defined in the centreof-mass frame of the hard process (or the transverse momentum of the associated extra parton) and µ is the renormalization scale. This is to accommodate the difference between the scale that the shower uses in the calculation of the strong coupling constant α S wrt. the scale used in the matrix elements. In practice the effect of the reweighting is small, especially in comparison to the uncertainties arising from variations of µ and the merging scale parameters. All of the results in the rest of the paper include α S -reweighting.

Systematic uncertainties
In what follows we present results obtained at parton level, using the Rivet analysis framework version 1.8.3 [61] and the anti-k T algorithm [62] with R = 0.4. In all the calculations, the Higgs boson mass was chosen to be m h = 125 GeV and the top quark and bottom quark masses to be 174.2 GeV and 4.7 GeV respectively, with all the widths set to zero. The renormalization and factorisation scales were both set to equal µ: µ F = µ R = µ.
We first examine the effect of the two schemes suggested to facilitate smooth merging, 'Sinusoidal' and 'Uniform', for E T Clus = 60 GeV and clus = 30 GeV. We compare to the  Transverse momentum of leading jet dσ/dp ⊥ (jet 1) [pb/GeV] purely showered LO sample ('un-merged') with µ = m h . In figure 3 we show the transverse momentum of the di-Higgs system and the transverse momentum of a Higgs boson, p hh ⊥ and p h ⊥ respectively, the distance between the two Higgs bosons, ∆R(h, h) and the p ⊥ of the leading jet. It is evident by examining the plots that the 'Uniform' scheme provides stronger smoothing than the 'Sinusoidal' scheme. Considering that the disagreement between the un-merged sample and the merged samples in the merging regions is large, we suggest the use the 'Uniform' scheme for merging in the hh process and we employ this in the rest of this study.
In figure 4 we examine the effect of the variation of the scale µ on various distributions for the merged samples. We vary the scale µ between µ = m h + p hh ⊥ and µ = 4(m h + p hh ⊥ ), while fixingĒ Tclus = 60 GeV and clus = 30 GeV. For comparison, we also show the equivalent un-merged scale variation between µ = m h and µ = 4m h . It is evident that, modulo normalisation differences originating from the scale variation, the general shapes of the distributions exhibit reasonable stability over the range of the chosen scales for  the merged sample. More importantly, the scale variation in the observables p hh ⊥ , p ⊥ of the leading jet and ∆R(h, h), is substantially reduced with respect to the leading-order showered samples. This is particularly true in the regions where the parton shower is not expected to provide a good description of the additional radiation, namely in the ∆R(h, h) π region and the high-p ⊥ regions. These improvements should not come as a surprise, as the considered observables are leading-order accurate for the merged sample in those regions, versus leading-logarithmic for the showered leading-order sample. The distribution of transverse momentum of a single Higgs boson, p h ⊥ , is not dominated by the extra radiation and thus the improvement is only marginal. Transverse momentum of leading jet dσ/dp ⊥ (jet 1) [pb/GeV] We again compare to the un-merged sample with µ = m h . The average merging scale was set toĒ Tclus = 60 GeV and the scale µ was set to µ = m h + p hh ⊥ . Evidently, smoother merging of the samples can be achieved using higher values of clus .
In figure 6 we vary the average clustering scale,Ē Tclus , while keeping µ = m h + p hh ⊥ and clus = 30 GeV. We again compare to the un-merged sample with µ = m h . Evidently, a relatively large systematic uncertainty comes from varyingĒ Tclus . This is due to the fact that changingĒ Tclus alters the regions that the parton shower and matrix element calculation contribute in. For a lowerĒ Tclus , the transition is smoother and we will usē E Tclus = 50 GeV for the phenomenological studies of the next section. We will contrast this toĒ Tclus = 70 GeV where appropriate.  Transverse momentum of leading jet dσ/dp ⊥ (jet 1) [pb/GeV] Finally, we show the envelope of the variation of both the merging scaleĒ Tclus (in [50,70] GeV) and the variation scale clus (in [10,30] GeV) in figure 7. In all of the samples, the scale was set to µ = m h + p hh ⊥ . The ratio sub-plot in the figure is taken with respect to the case where E Tclus = 60 GeV and clus = 20 GeV and the error bars represent the Monte Carlo statistical uncertainty on that sample. For the p hh ⊥ distribution, the uncertainty due to the variation of E Tclus and clus is O(20%) up to p hh ⊥ ∼ 250 GeV and grows to over ∼ 40% at higher values. The  the region where the merging scale becomes significant, [50,70] GeV, and are then reduced to O(10 − 20%) variations up to p ⊥ ∼ 300 GeV.

Phenomenological implications
It is important to examine the implications of the merging on realistic phenomenological analyses of Higgs boson pair production at the LHC. We do this by focussing on an example of a decay channel with a relatively large branching ratio, hh → (bb)(τ + τ − ). This has been examined in detail in [14,15,30]. We do not attempt here to perform a detailed signal versus background study; instead, we wish to show the magnitude of the effect of using the merged sample in a realistic analysis. We only focus on the top-anti-top background We consider the case of a 14 TeV LHC, and normalise all hh inclusive cross sections to the NNLO cross section obtained within the effective theory in [20], σ N N LO hh = 40.2 fb. We consider four different samples, un-merged with scales set to µ = m h and µ = 2m h and merged with scales µ = m h + p hh ⊥ and µ = 2(m h + p hh ⊥ ). The merging parameters were fixed to E Tclus = 50 GeV (or 70 GeV) and clus = 30 GeV. The tt background was generated via aMC@NLO [63,64] along with the decays, and was assumed to have a total cross section of σ tt = 900 pb [65,66]. Showering and hadronization were performed using HERWIG++, and the simulation of the underlying event was included via multiple secondary parton interactions [67]. We follow the basic analysis steps as given in [30]: we assume 80% τ -reconstruction efficiency with negligible fake rate 6 and require two τ -tagged jets with at least p ⊥ > 20 GeV. We require that the taus, taken from the Monte Carlo truth, reproduce the Higgs mass within a 50 GeV window, to account for the reconstruction smearing, as done in [30]. We use the Cambridge-Aachen jet algorithm available in the FastJet package [68,69] with a radius parameter R = 1.4 to search for so-called 'fat jets'. We require the existence of one fat jet in the event satisfying the mass-drop criteria as done in the hV study in ref. [70]. We require the two hardest 'filtered' sub-jets to be b-tagged 7 and to be central (|η| < 2.5) and the filtered fat jet to be in (m h − 25 GeV, m h + 25 GeV). The b-tagging efficiency was taken to be 70%, again with negligible fake rate for the sake of simplicity. We require a loose cut on the transverse momentum of the fat jet (after filtering) that satisfies the above criteria, p fat ⊥ > 100 GeV. This is done to maintain a sufficient number of events to examine the change of efficiencies with respect to other cuts. We also apply a transverse momentum cut on the τ + τ − system of equal magnitude, p τ τ ⊥ > 100 GeV. We wish to examine the stability of the merged samples against that of the un-merged samples with respect to scale variations. It is obvious that sufficiently inclusive quantities should not differ in a way that will impact the analyses. However, there are quantities for which the merged sample and the un-merged sample differ substantially. As an exercise, we examine two such observables here: the distance between the (τ + τ − ) system and the (bb) system (equivalent to the distance between the Higgs bosons), and the transverse momentum of the τ + τ − bb system (equivalent to the transverse momentum of the the hh system). Figure 8 demonstrates that it is conceivable that both of these observables may be used for background rejections. Moreover, it is evident that the largest uncertainties in the hh signal predictions are present in the exact same region that one would wish to place the cuts in: p T (hh) ∼ 100 GeV and ∆R(h, h) ∼ 3. What is also important is the fact that both the 0-jet exclusive and the 1-jet inclusive signal samples contribute in the region of interest, as demonstrated in figure 9 for the case where µ = m h + p hh ⊥ . We can further quantify the effect observed in figure 8. The stability of the samples can be assessed by taking ratios of the efficiencies for different cut values. If the efficiency does JHEP03(2014)126  not vary substantially with the cuts, then the sample (merged or un-merged) is deemed to be stable, and the theoretical systematic uncertainty on the efficiency of the cut can be considered to be low. Figure 10 shows the variation of the ratio of efficiencies for different parameters for the merged and un-merged samples: R = (cut efficiency, sample i)/(cut efficiency, sample j), (4.2) for cuts on the aforementioned observables which we abbreviate as p ⊥,max and ∆R min (h, h). For details of the parameters used for each of the samples {α, β, γ, δ, ι, κ}, see the caption of figure 10. The un-merged samples {κ, λ} exhibit a fairly substantial change in the ratio of efficiencies for the two chosen scales, starting from 10% and going up to ∼ 20% for some values of the cuts. This change can be interpreted as a theoretical systematic uncertainty on the efficiency itself. The merged samples {α, β, γ, δ} perform better, with lower overall variation of the efficiency ratio, with the deviations always < 10% as demonstrated in the figure 10, while, more importantly, on average possessing an efficiency variation of ∼ 5%. The differences are due to the fact that the chosen observables are sensitive to the behaviour of the extra radiation, which, at high transverse momentum or large separations between the Higgs bosons, is not predicted reliably by the parton shower. 8 For completeness, in table 2 we show a set of cuts and resulting cross sections resulting from the analysis. We provide also explicit cuts on the variables we have examined: ∆R(h, h) > 2.8 and p hh ⊥ < 80 GeV. The final result for this basic analysis is optimistic, with S/B ∼ 0.4 − 0.5 for all samples, leading to a reasonable significance at 600 fb −1 of integrated luminosity, expected to be collected at the LHC in the full phase of Run II. It is worthwhile to note that the final cross section prediction for the merged samples (α and β) after all cuts only exhibits a ∼ 2% variation compared to the un-merged samples (κ and λ) which exhibit a ∼ 13% variation. 8 The ratio differs from unity since we expect differences in response to cuts in other observables between the two samples. This can also be seen in the cut flows of table 2. Nevertheless, the merged sample ratios are still closer to unity overall than those of the un-merged samples, demonstrating further the increased accuracy of the calculation.   Table 2. Cross sections for the hh signal and tt aMC@NLO background after series of cuts. The un-merged samples κ and λ have µ = m h and µ = 2m h respectively and the merged signal samples 'α' and 'β' have µ = m h + p hh ⊥ and µ = 2(m h + p hh ⊥ ) respectively, as well as E Tclus = 50 GeV and clus = 30 GeV. The final two cuts were chosen to be ∆R(h, h) > 2.8 and p hh ⊥ < 80 GeV.

Conclusions
We have described the implementation of Higgs boson pair production merged to the onejet matrix elements, generated using OpenLoops, in the HERWIG++ event generator. We have examined the systematic uncertainties associated with the merging. Moreover, we have provided examples of the magnitude of the effects of using the merged samples in a realistic analysis. As was demonstrated in this analysis, using the leading order matrix elements in conjunction with the parton shower can potentially introduce O(20%) systematic uncertainties in the predictions of the efficiencies of experimental cuts. The uncertainty will inexorably propagate to measurements of the Higgs boson self-coupling. The merged samples demonstrate theoretical uncertainties on the efficiencies that are 10% or better for the examined observables. We expect such conclusions to remain valid for a future NLO simulation matched to the parton shower. We thus recommend the use of samples that include the merged exact one-jet matrix elements in all future phenomenological or experimental analyses of the process. The Monte Carlo event generator developed for this project is available as an add-on to the HERWIG++ event generator at http://www.physik.uzh.ch/data/openloops/download/projects/hhmerge/.