Parton Shower and NLOMatching Uncertainties in Higgs Boson Pair Production
Abstract
We perform a detailed study of NLO parton shower matching uncertainties in Higgs boson pair production through gluon fusion at the LHC based on a generic and process independent implementation of NLO subtraction and parton shower matching schemes for loopinduced processes in the Sherpa event generator. We take into account the full topquark mass dependence in the twoloop virtual corrections and compare the results to an effective theory approximation. In the full calculation, our findings suggest large parton shower matching uncertainties that are absent in the effective theory approximation. We observe large uncertainties even in regions of phase space where fixedorder calculations are theoretically well motivated and parton shower effects expected to be small. We compare our results to NLO matched parton shower simulations and analytic resummation results that are available in the literature.
MPP2017249
SLACPUB17190
Introduction
In the Standard Model (SM) the production of a pair of Higgs bosons at a hadron collider proceeds dominantly through the annihilation of gluons. As there is no direct coupling between the Higgs boson and gluons this process is mediated by intermediate massive quark loops. The topquark loop contributions dominate by far, due to the large Yukawa coupling. Bottomquark loops only contribute at the level to the total cross section at leading order (LO) and can thus safely be neglected in most situations. The scattering amplitudes relevant for the calculation of the topquark contributions up to nexttonextto leading order (NNLO) are known in the approximation of an infinitely heavy topquark (commonly referred to as the HEFT approximation) [11, 14, 13, 26, 25, 12]. However, the validity of this approximation is questionable for Higgs boson pair production due to the large momentum transfer required to produce the Higgs bosons. Techniques to systematically improve upon it were extensively studied in [19, 33, 24, 26, 25, 32]. The full result, which is exact in the mass of the topquark, was known only to leading order (LO) [16, 23, 38] until recently. This is due to the complexity of computing the nexttoleading order (NLO) virtual corrections which feature twoloop, fourpoint integrals with both massive internal propagators and massive external lines. They have recently been calculated through the numerical evaluation of all relevant twoloop integrals as part of the full NLO calculation in [6, 5].
At small transverse momenta of the Higgs boson pair, the accuracy of any fixedorder calculation is spoiled by the presence of large logarithms of the form . They can be resummed to all orders using analytical resummation techniques which have been applied to Higgs boson pair production in [18]. Alternatively, parton shower simulations can be employed. In addition to providing a reliable transverse momentum spectrum at small , they also provide results that are fully differential in the kinematics of any soft and collinear QCD radiation. Standard techniques exist for the consistent matching of NLO fixed order calculations to parton shower simulations [20, 34]. They were recently applied to Higgs boson pair production in reference [27], where the MC@NLO and POWHEG matching techniques were used to combine the fixedorder NLO calculation with the Pythia parton shower [41, 40]. The results of [27] suggest that the parton shower matching can have sizeable effects not only in the region of small , but also in the region of large , where one would expect the fixedorder calculation to be reliable and the approximations inherent to parton shower simulations to break down. These effects even exceed the scale uncertainties of the fixedorder calculation.
In this publication we aim to critically assess the origin and size of the aforementioned effects and associated uncertainties. For this purpose we implemented a fully generic and process independent NLO subtraction along with the corresponding parton shower matching techniques for loopinduced processes in the Monte Carlo event generator Sherpa [21]. This allows us to perform our studies using the two different showers that are implemented in Sherpa [39, 42] within the same parton shower matching framework.
This publication is structured as follows. In Section 1 we describe in detail the setup of our calculation along with a brief review of the MC@NLO matching technique and the parton showers we used for our studies. We present the results of our simulations in Section 2, focusing on the origin and size of uncertainties that are inherent to the matching technique applied. We also point out crucial differences that arise when going from the HEFT approximation to the full calculation. Our conclusions are presented in Section 3.
1 Calculational Setup
1.1 FixedOrder NLO Calculation
For the virtual twoloop amplitude we utilize the result of reference [6, 5], retaining the full finite topquark mass effects. This amplitude was obtained by numerically evaluating all relevant 2loop 4point Feynman diagrams with up to 4 scales. We adopt the input parameters of reference [5], with , the mass of the topquark set to , the Higgs boson mass set to , and their widths neglected. We also adopt the choice of reference [5] for the factorization and renormalization scales . Perturbative uncertainties in the fixedorder part of the calculation are estimated by independently varying these scales through factors of and . All studies are performed with hadronic centerofmass energy . The NLO virtual amplitude is provided in the literature in the form of an interpolation grid in two Mandelstam variables, based on a fixed number of precomputed phasespace points [27]. We extract the finite part of the UV renormalized virtual amplitude in the Conventional Dimensional Regularization (CDR) scheme with residual IR singularities subtracted according to the Catani and Seymour scheme [10], as required by the Sherpa event generator, using relations (2.5) and (2.6) of reference [27].
The leading order oneloop squared amplitudes for the Born process and real emission contributions are provided by OpenLoops [9]. For the evaluation of tensor and scalar oneloop integrals, we employ the Collier library [15], CutTools [37], and OneLOop [44, 43].
For the regularization and numerical cancellation of infrared divergences in the realemission part of the calculation we employ the dipole subtraction scheme of Catani and Seymour [10]. We have reimplemented this scheme within Sherpa in a fully generic and processindependent way for loopinduced processes. This implementation is qualitatively equivalent to the implementation in one of Sherpa’s internal matrix element generators Amegic++ [31, 22], apart from the fact that color and spincorrelated amplitudes are to be provided externally through generic interfaces. Through a dedicated interface to OpenLoops and the aforementioned tools, NLO calculations for loopinduced SM processes are thus fully automated (given the availability of the virtual twoloop corrections) in Sherpa and will become available in a public version of the code.
We have validated our implementation in Sherpa by comparing our results for the total cross section, for the differential Higgs boson pair invariant mass distributions, and for the differential single Higgs boson transverse momentum distributions to those published in reference [5].
1.2 Parton Showers
We consider two parton showers for matching to the fixedorder NLO calculation. Both algorithms are dipoletype showers in which QCD radiation is generated coherently off color dipoles spanned by pairs of preexisting partons. Both showers are publicly available as part of the Sherpa event generator package.
Due to the dipole character of the parton showers, their splitting kernels can be used for the purpose of fixedorder NLO subtraction, thus simplifying the implementation of parton shower matching. The CS shower [39] directly uses the splitting kernels of the original CataniSeymour subtraction scheme for parton evolution. The Dire shower [42] uses splitting kernels that are modified in such a way as to reproduce the collinear anomalous dimensions of the DGLAP equations. For NLO matching to the Dire shower, we use a modified version of the original CataniSeymour subtraction scheme that reflects these changes in the splitting kernels [29]. The kernels of both showers approximate real emission amplitudes arbitrarily well in the limit of soft and collinear momenta. Away from the soft and collinear regions, however, they differ.
A further crucial difference between the two algorithms is the choice of evolution variable, which we generically denote by in the following. The choice of evolution variable together with the shower starting scale dictates how much of the phase space away from the soft and collinear regions is available to the parton shower since the starting scale implements the following phase space constraint:
(1) 
For the discussion of the evolution variables we focus on the first (hardest) emission in the production of a colorneutral final state of invariant mass . We illustrate the kinematics of the first emission, producing a final state parton with momentum from the collision of two incoming massless partons with momenta and , in Figure 1. It is useful to consider the variables and , which are closely related to the standard Mandelstam variables , , and :
(2) 
Due to momentum conservation which implies,
(3) 
In terms of and , the evolution variable in Dire is given by
(4) 
The functional form (4) implies that due to equation (3). It follows that for a parton shower starting scale of , Dire behaves like a “power shower” in the sense that it populates the full phase space since and thus (1) is trivially fulfilled.
In the CS shower, the evolution variable is given by
(5) 
This implies that, for a given kinematic configuration, is typically larger than , such that for a given value of the emission phase space of the CS shower is more restricted than that of Dire. It is worth noting that in particular does not correspond to a “power shower” when employing the CS shower. This choice in fact severely constrains the emission phase space, since can get close to and thereby give rise to large values of .
1.3 NLO Parton Shower Matching
In the following we will focus on the MC@NLO matching prescription [20] using the notation of [30] with no distinction between fixedorder NLO subtraction terms and parton shower matching terms since we use the parton shower splitting kernels both for parton evolution and for infrared subtraction and keep all phase space constraints explicit. We thus denote the sum of subtraction terms as a function of the real emission phase space by , where the real emission phase space can be decomposed into the Born phase space and an extra oneparticle emission phase space . We then define the fixedorder differential seed cross sections and in terms of the leading order (Born) term , the UVsubtracted virtual corrections , and the realemission contributions by
(6)  
(7) 
where is the map from a kinematic real emission configuration to the parton shower evolution variable . The Heaviside function in (7) implements the constraint (1). For notational convenience, we will omit the explicit dependence and write in the following. In terms of the quantities introduced above, the fixedorder total NLO cross section is given by
(8) 
In MC@NLO, we generate events according to
(9) 
where is the infrared cutoff scale of the parton shower and the modified Sudakov form factor , which gives the probability for no emission to occur between scales and for the first parton shower step, is given by
(10)  
(11) 
The first line of (9) corresponds to events that have Born kinematics at the level of the fixedorder seed event with weight (Sevents). They either don’t undergo any emission above the infrared parton shower cutoff scale (first term in the square bracket) or they undergo their hardest emission at some scale between and (second term in the square bracket). The second line of (9) corresponds to events with realemission kinematics at the level of the fixedorder seed event and weight (Hevents). All events are treated further by the parton shower precisely as in the leading order case, apart from the Sevents that haven’t undergone any emission, which are kept as they are.
Since the square bracket in (9) integrates to , the total cross section and any observable that is insensitive to QCD radiation is unaltered in MC@NLO compared to the fixedorder NLO result. In fact, it can be shown that a MC@NLO event sample will reproduce the fixedorder NLO result event to order relative to the Born for any infrared safe observable [20]. The parametric NLO accuracy is therefore not spoiled by the parton shower matching.
1.4 Parton Shower Matching Uncertainties
As stated in the previous section, NLO parton shower matching according to the MC@NLO method preserves the parametric accuracy of the fixedorder NLO calculation. Deviations from fixedorder results can numerically be nonetheless significant [30]. Such differences reflect genuine parton shower matching uncertainties, they can be particularly prominent for observables that are sensitive to real emission configurations and thereby to the interplay between parton shower emissions and fixedorder real emission configurations. We will therefore focus on the distribution in the following section, comparing MC@NLO matched parton shower simulations to fixedorder results with both the Dire and the CS shower.
In order to formally compare the MC@NLO result to a fixedorder prediction for this spectrum, we first consider a generic observable that is insensitive to kinematic Born configurations. For such an observable we need to take into account Hevents and parton shower emissions off Sevents. At order relative to the Born we have
(12) 
where the first integral corresponds to Sevents in which the parton shower has generated a nonvanishing value of and the second integral corresponds to Hevents, where a nonvanishing value of is implied by the realemission kinematics of the fixedorder seed event. In the tail of the distribution where we can neglect the Sudakov suppression and set , we obtain after plugging in the definition of :
(13) 
To order we have and the first integral cancels as required by the matching conditions, thus restoring the fixedorder result. This explicitly demonstrates how variations in the parton shower contributions induced by Sevents are subtracted by the MC@NLO subtraction terms in the definition of . Numerically, however, this cancellation can be severely spoiled, potentially leading to large deviations from the fixedorder result. For the deviations to be significant the term on the first line of equation (13) must be similar in size to the fixedorder term on the second line of (13). One can therefore expect large deviations from the fixedorder calculation only if both of the following conditions are met:

The factor dressed with the parton shower splitting kernels ( in (13)) is comparable in magnitude relative to the realemission matrix elements in . This depends on the size of the NLO corrections that enter and on the splitting kernels in the phase space region of hard emissions.

The phase space of interest is accessible to the parton shower so that the first integral in (13) has support in that region. This depends on the choice of and on the shower (through the functional form of ).
The formally subleading contributions originating from the parton shower matching in the first integral in (13) are, to a large extent, controlled by the choice of . To access the matching uncertainties we will therefore vary this parameter by factors and . With two different parton showers at our disposal we have an additional handle on these uncertainties through the functional form of . The nominal choice for in the CS shower will be , in line with and . As outlined above, such a choice would open up the entire emission phase space in case of the Dire shower. Our nominal choice for the Dire shower will therefore be , which allows us to perform both up and downwards variations.
Based on the argument presented above, one might expect to see large parton shower contributions in the high tails of other processes with large Kfactors. However, it is important to note that for such effects to be visible the factor must remain large, relative to the realemission matrix element, also when multiplied by the parton shower splitting kernels. In single Higgs boson production through gluon fusion, for example, one might anticipate a large shower contribution in the tail of the Higgs transverse momentum spectrum due to the large NLO Kfactor. However, in this case the parton shower splitting kernels underestimate the realemission matrix elements significantly, such that the parton shower contributions in the tail are very small in an MC@NLOmatched calculation [1, 30]. For the parton showers considered in this work, this holds even when taking into account the full top mass dependence in the realemission matrix elements, which reduce the size of by more than an order of magnitude in the tail of the transverse momentum distribution.
2 Results
2.1 Leading Order Results
We start our discussion with predictions obtained in the most simple setup, using leading order matrix elements for inclusive Higgs boson pair production supplemented by a parton shower. This type of simulation will be referred to as “LO+PS” in what follows. Since the transverse momentum of the Higgs boson pair is zero at leading order, any nonzero value of this observable is entirely generated by the parton shower. As a reference, we use a fixedorder prediction obtained by simulating the process with leading order matrix elements.
Figure 2 and 3 show the result of our comparison both in the HEFT approximation and in the full theory, respectively. Comparing the full SM and the HEFT approximation, we observe qualitatively different parton shower effects. In the HEFT approximation, both parton showers significantly underestimate the fixed order spectrum in the tail of the distribution. Even if the full phase space is made available to the showers they do not reproduce the slowly falling transverse momentum spectrum predicted by the fixedorder HEFT matrix elements. To show this for the CS shower we display also LO+PS results obtained by setting the parton shower starting scale to the hadronic centerofmass energy . In the case of the Dire shower, the full phase space is already available for , which corresponds to the upper edge of the uncertainty band.
In the full SM, by contrast, for large enough values of the parton shower starting scale both parton showers overestimate the fixedorder prediction. For the CS shower, this effect is restricted to smaller transverse momenta, due to the choice of evolution variable. If we lift any phase space restriction in the CS shower, by setting , we observe that in the tail of the distribution the shower overestimates the fixedorder predictions by more than an order of magnitude. The upper edge of the Dire shower uncertainty band also overestimates the fixedorder prediction, although this feature is not as pronounced as for the CS shower.
It is therefore evident that the Born matrix elements dressed with splitting kernels can strongly overestimate the real emission matrix elements in the full SM and strongly underestimate the real emission matrix elements in the HEFT. The former effect is, however, to a certain extent limited by the phase space constraint implemented through the parton shower starting scales.
Naively, the large differences between the HEFT and the full SM simulations may seem surprising. However, the highenergy behaviour of the HEFT real emission amplitudes is unphysical because the momentum transfers vastly exceed the topquark masses which have been integrated out in the HEFT approximation. As a result, the spectrum calculated at fixedorder falls off extremely slowly in the HEFT and the parton shower kernels thus tend to underestimate the spectrum in the tail. A similar slow fall off can been observed in the HEFT approximation for the Higgs boson transverse momentum in single Higgs boson production [4, 17, 7, 36, 8].
2.2 NLO Results
We start the discussion of NLOmatched parton shower simulations with the results of a HEFT treatment, shown in Figure 4. As discussed in Section 2.1 (and shown in Figure 2) the combination of Born matrix elements and parton shower splitting kernels strongly undershoot the full realemission matrix elements when employing the HEFT approximation. We therefore expect the tail of the distributions to converge to the fixedorder result both for the CS shower and the Dire shower. As shown in Figure 4, this is indeed the case.
In addition to a more precise description of the tail (compared to the LO+PS type simulations) we observe a reduction of the parton shower starting scale uncertainties. The individual variations of H and Sevent contributions are of order one for some values of but cancel to a large extent in the sum.
Moving on to the discussion of results in the full SM, we remind the reader of our findings in the corresponding LO+PS type simulations. In the full SM, the parton shower splitting kernels in combination with Born matrix elements overestimate the realemission matrix elements (see Figure 3). The parton shower effects in the tail of the distribution are therefore large. As shown in Figure 5, this also holds at NLO. The parton shower matched results converge to the fixed order result in the tail for nominal choices for . Upward variation of , however, (indicated by the upper edge of the blue uncertainty bands) lead to parton shower effects of up to + even in the tail of the distribution. As shown in the lower panels of Figure 5, the excesses in the tail are indeed generated by parton shower emissions off Sevents. The extent of these effects is limited by the phase space available to the parton shower, as determined by the choice of and the functional form of the evolution variable. We observe that results generated using the Dire shower, particularly for larger values of , have a different shape than those generated with the CS shower.
For the MC@NLO algorithm, by construction, the large parton shower effects in the tail should be cancelled to first order in . As outlined in Section 1.4, any miscancellation is due to a numerically large discrepancy between and . We demonstrate this explicitly in Figure 6, where we show modified Dire MC@NLO with substituted for , leading to a complete cancellation of the first integral in (13). This procedure eliminates large parts of the excess in the tail independently of , as anticipated. Variations in S and Hevent contributions remain large, as shown in the lower panels of Figure 6, but they cancel in the sum. The procedure of replacing with would, of course, spoil the NLO accuracy of any inclusive observable but allows us here to demonstrate the origin of the discrepancy between the showered and fixedorder results in the tail of the distribution.
In the HEFT approximation one may naively expect effects of similar size in the tail of the distributions. As demonstrated using a LO+PS simulation and as shown in Figure 2, however, the fixedorder real emission contributions completely dominate in this region. The bulk of the contributions in the tail are hence generated by the second integral of (13). As a result, the relative impact of parton shower effects in the tail remains small as shown in Figure 4.
2.3 Comparison to the Literature
In Figure 7 we compare our results for the spectrum to the NLO parton shower matched results presented in reference [27]. These results were obtained with the Pythia 8 shower [41, 40] interfaced to MadGraph5_aMC@NLO [3, 28] and POWHEG BOX [2] for matching according to the MC@NLO method and the POWHEG method [34], respectively. In MadGraph5_aMC@NLO the nominal value of is set randomly in the interval , where is the sum of the transverse masses of the Higgs bosons. For the simulations based on MadGraph5_aMC@NLO and Pythia we show uncertainty bands that were obtained by varying the nominal parton shower starting scale by factors of and . The POWHEG matching prescription can be recovered within the MC@NLO framework by setting the parton shower starting scale to the collider energy and by setting [30]. Therefore, lacking a natural equivalent to in the POWHEG framework, we compare only to nominal POWHEG predictions produced with the hdamp parameter set to , as described in [27].
Focusing on the region , we note that all MC@NLO predictions considered here are generally compatible within the uncertainty bands. However, the agreement between the nominal results of our simulations and the fixedorder result is much better in the tail. The POWHEG results exhibit a very large excess in the tail that is not covered by the uncertainty bands of our Sherpa predictions. Similar discrepancies between MC@NLO and POWHEG have been observed in the context of other processes [1, 30] and can be attributed to the large phase space available to the parton shower as a result of setting and the numerically large discrepancy between and [35]. As described in Section 1.2, the former can be achieved in Dire by setting , which is represented by the upper edge of the uncertainty band around the Dire prediction. We note that the shape of this curve is in fact most similar to the POWHEG prediction in the tail of the distribution.
Comparing the different uncertainty bands themselves, we observe large differences. The shape of uncertainty bands obtained with MadGraph5_aMC@NLO and with the CS shower are somewhat similar with a peak around but the size of the uncertainty band around the MadGraph5_aMC@NLO result is much larger throughout. The uncertainties on the Dire prediction describe a more evenly shaped band.
Differences in the region of small transverse momenta are not fully covered by the variation bands. Although we expect these variations to be indicative of NLOmatching uncertainties, we do not expect them to cover all parton shower uncertainties. These include, but are not limited to, ambiguities in the choice of a kinematic recoil scheme and ambiguities in the choice of the renormalization scales for the strong coupling in the splitting kernels.
In Figure 8 we show a comparison to the calculation of reference [18] which employed analytic nexttoleadinglog (NLL) resummation techniques instead of parton showers. We observe good agreement within the uncertainties except near the peak region and the region around where we find discrepancies of about that are not fully covered by our uncertainty bands. Taking into account the resummation scale uncertainties on the analytic results (not shown in Figure 8), which are of the order of in the peak region and near [18], we consider the observed agreement satisfactory.
2.4 Other Observables
Having discussed the Higgs boson pair transverse momentum at length, we briefly discuss parton shower effects on a number of other observables.
Lorentz invariant observables that depend only on the momenta of the Higgs bosons are not affected by the parton shower. The kinematics of the Higgs boson pair is only altered by emissions off dipoles spanned by two initial state partons. The recoil generated by such an emission affects the Higgs boson pair only through a Lorentz boost. Lorentzinvariant quantities like the Higgs boson pair invariant mass are therefore not affected by the parton shower. Our simulations were checked by inspecting the Higgs boson pair invariant mass distributions, we observe agreement within the statistical uncertainties of well below one percent.
Figure 9 shows the leading jet transverse momentum . As opposed to the Higgs boson pair, the parton emitted in the hardest emission off an Sevent and the final state parton in an H event is affected strongly by secondary emissions because the kinematics are altered by finalfinal and finalinitial dipoles. Such emissions decorrelate the leading jet and the Higgs boson pair momenta. Parton shower effects are therefore qualitatively different. The effect of the parton shower are generally moderate and the spectrum remains compatible with the fixedorder prediction within the scale uncertainties. It is worth noting that even the full variation bands remain within the uncertainty bands of the fixedorder calculation.
This picture drastically changes when considering observables that are more sensitive to high multiplicity final states. As an example we consider the differential distribution, defined as the scalar sum of jet transverse momenta:
where the index labels all jets in the respective event. In a parton shower event, the total energy is typically distributed among a larger number of jets than in a fixedorder calculation. Their scalar contributions to are added, giving rise to larger values of than for the fixedorder NLO calculation. We show this effect in Figure 9. Comparing the Dire and CS shower predictions, we note that the uncertainty bands overlap but that the shower starting scale variations are much larger for the CS shower.
In Figure 10 we show the azimuthal separation between the Higgs bosons . At leading order, the momenta of the Higgs bosons are perfectly correlated due to momentum conservation. Only in events with additional radiation can one observe a nontrivial distribution of the azimuthal separation between the Higgs bosons. As shown in Figure 10, parton shower corrections to the fixedorder result are mostly covered by the fixedorder uncertainties except in the region of which corresponds to backtoback configurations and which is sensitive to soft QCD emissions.
Also shown in Figure 10 is the transverse momentum of a randomly chosen Higgs boson. The effect of a parton shower emission on the transverse momentum of a given Higgs boson is random, either decreasing or increasing its value. If the distribution was completely flat, any parton shower effects would therefore average out. Since the distribution is falling, the parton shower effects of increasing the transverse momenta of low Higgs bosons is not counterbalanced by the effect of decreasing the transverse momenta of high Higgs bosons, thus inducing a slope relative to the fixedorder result. This effect is small but clearly visible in Figure 10.
3 Conclusions
We have presented a study of NLO parton shower matching uncertainties in Higgs boson pair production through gluon fusion at the LHC. We assessed these uncertainties by matching the fixedorder NLO calculation to two dipole shower algorithms in the Sherpa event generator according to the MC@NLO framework. The interplay between fixedorder real emission contributions and parton shower emissions was studied in detail through variations of the parton shower starting scale. We find large matching uncertainties that exceed the fixedorder uncertainties even in regions of phase space where the fixedorder calculation is well motivated and where parton shower matching effects are expected to be small. Our nominal predictions are in good agreement with the fixedorder result in these regions, however. A comparison to MC@NLO matched results from the literature revealed qualitative differences which are, nevertheless, compatible within the large uncertainties. We observe larger differences in a comparison to POWHEG predictions in the tail of the transverse momentum spectrum, where POWHEG overestimates the fixedorder spectrum by a factor of . We find reasonable agreement throughout between our results and those obtained through analytic resummation techniques.
Acknowledgements
We would like to thank Stefan Höche for providing analytic results for the subtraction and matching terms required in the interface to the Dire shower, for technical assistance in the implementation, and for helpful discussions. We thank Eleni Vryonidou for evaluating the MadGraph5_aMC@NLO + Pythia parton shower uncertainties and Matthias Kerner for his advice regarding the use of the fixedorder NLO results. We would also like to thank Jonas Lindert and Philipp Maierhöfer for their support with OpenLoops. We are indebted to Gudrun Heinrich, Joey Huston, Matthias Kerner, and Frank Krauss for their valuable comments on the draft. This work is supported in part by the U.S. Department of Energy under contract number DEAC0276SF00515 and by the Research Executive Agency (REA) of the European Union under the Grant Agreement PITNGA2012316704 (HiggsTools).
References
 [1] Simone Alioli, Paolo Nason, Carlo Oleari, and Emanuele Re. NLO Higgs boson production via gluon fusion matched with shower in POWHEG. JHEP, 04:002, 2009.
 [2] Simone Alioli, Paolo Nason, Carlo Oleari, and Emanuele Re. A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX. JHEP, 06:043, 2010.
 [3] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro. The Automated Computation of TreeLevel and NexttoLeading Order Differential Cross Sections, and their Matching to Parton Shower Simulations. JHEP, 07:079, 2014.
 [4] U. Baur and E.W.N. Glover. Higgs boson production at large transverse momentum in hadronic collisions. Nuclear Physics B, 339(1):38 – 66, 1990.
 [5] S. Borowka, N. Greiner, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk, and T. Zirke. Full top quark mass dependence in Higgs boson pair production at NLO. JHEP, 10:107, 2016.
 [6] S. Borowka, N. Greiner, G. Heinrich, S.âP. Jones, M. Kerner, J. Schlenk, U. Schubert, and T. Zirke. Higgs Boson Pair Production in Gluon Fusion at NexttoLeading Order with Full TopQuark Mass Dependence. Phys. Rev. Lett., 117(1):012001, 2016. [Erratum: Phys. Rev. Lett.117,no.7,079901(2016)].
 [7] Malte Buschmann, Dorival Goncalves, Silvan Kuttimalai, Marek Schonherr, Frank Krauss, and Tilman Plehn. Mass Effects in the HiggsGluon Coupling: Boosted vs OffShell Production. JHEP, 02:038, 2015.
 [8] Fabrizio Caola, Stefano Forte, Simone Marzani, Claudio Muselli, and Gherardo Vita. The Higgs transverse momentum spectrum with finite quark masses beyond leading order. JHEP, 08:150, 2016.
 [9] Fabio Cascioli, Philipp Maierhofer, and Stefano Pozzorini. Scattering Amplitudes with Open Loops. Phys. Rev. Lett., 108:111601, 2012.
 [10] S. Catani and M.H. Seymour. A general algorithm for calculating jet crosssections in nlo qcd. Nuclear Physics B, 485:291–419, 1997.
 [11] S. Dawson, S. Dittmaier, and M. Spira. Neutral Higgs boson pair production at hadron colliders: QCD corrections. Phys. Rev., D58:115012, 1998.
 [12] Daniel de Florian, Massimiliano Grazzini, Catalin Hanga, Stefan Kallweit, Jonas M. Lindert, Philipp Maierhöfer, Javier Mazzitelli, and Dirk Rathlev. Differential Higgs Boson Pair Production at NexttoNexttoLeading Order in QCD. JHEP, 09:151, 2016.
 [13] Daniel de Florian and Javier Mazzitelli. Higgs Boson Pair Production at NexttoNexttoLeading Order in QCD. Phys. Rev. Lett., 111:201801, 2013.
 [14] Daniel de Florian and Javier Mazzitelli. Twoloop virtual corrections to Higgs pair production. Phys. Lett., B724:306–309, 2013.
 [15] Ansgar Denner, Stefan Dittmaier, and Lars Hofer. Collier: a fortranbased Complex OneLoop LIbrary in Extended Regularizations. 2016.
 [16] Oscar J. P. Eboli, G. C. Marques, S. F. Novaes, and A. A. Natale. Twin higgs boson production. Phys. Lett., B197:269–272, 1987.
 [17] R.K. Ellis, I. Hinchliffe, M. Soldate, and J.J. Van Der Bij. Higgs decay to  a possible signature of intermediate mass higgs bosons at high energy hadron colliders. Nuclear Physics B, 297(2):221 – 243, 1988.
 [18] Giancarlo Ferrera and Joao Pires. Transversemomentum resummation for Higgs boson pair production at the LHC with topquark mass effects. JHEP, 02:139, 2017.
 [19] R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, P. Torrielli, E. Vryonidou, and M. Zaro. Higgs pair production at the LHC with NLO and partonshower effects. Phys. Lett., B732:142–149, 2014.
 [20] S. Frixione and B. R. Webber. Matching nlo qcd computations and parton shower simulations. Journal of High Energy Physics, 0206:029, 2002.
 [21] T. Gleisberg, Stefan. Hoeche, F. Krauss, M. Schonherr, S. Schumann, et al. Event generation with SHERPA 1.1. Journal of High Energy Physics, 0902:007, 2009.
 [22] T. Gleisberg and F. Krauss. Automating dipole subtraction for qcd nlo calculations. European Physical Journal C, 53:501–523, 2008.
 [23] E. W. Nigel Glover and J. J. van der Bij. Higgs boson pair production via gluon fusion. Nucl. Phys., B309:282–294, 1988.
 [24] Jonathan Grigo, Jens Hoff, Kirill Melnikov, and Matthias Steinhauser. On the Higgs boson pair production at the LHC. Nucl. Phys., B875:1–17, 2013.
 [25] Jonathan Grigo, Jens Hoff, and Matthias Steinhauser. Higgs boson pair production: top quark mass effects at NLO and NNLO. Nucl. Phys., B900:412–430, 2015.
 [26] Jonathan Grigo, Kirill Melnikov, and Matthias Steinhauser. Virtual corrections to Higgs boson pair production in the large top quark mass limit. Nucl. Phys., B888:17–29, 2014.
 [27] G. Heinrich, S. P. Jones, M. Kerner, G. Luisoni, and E. Vryonidou. NLO predictions for Higgs boson pair production with full top quark mass dependence matched to parton showers. JHEP, 08:088, 2017.
 [28] Valentin Hirschi and Olivier Mattelaer. Automated event generation for loopinduced processes. JHEP, 10:146, 2015.
 [29] Stefan Höche. Unpublished, personal communications.
 [30] S. Hoeche, F. Krauss, M. Schonherr, and F. Siegert. A critical appraisal of nlo+ps matching methods. Journal of High Energy Physics, 1209:049, 2012.
 [31] F. Krauss, R. Kuhn, and G. Soff. AMEGIC++ 1.0: A Matrix element generator in C++. JHEP, 02:044, 2002.
 [32] Philipp Maierhöfer and Andreas Papaefstathiou. Higgs Boson pair production merged to one jet. JHEP, 03:126, 2014.
 [33] F. Maltoni, E. Vryonidou, and M. Zaro. Topquark mass effects in double and triple Higgs production in gluongluon fusion at NLO. JHEP, 11:079, 2014.
 [34] Paolo Nason. A new method for combining nlo qcd with shower monte carlo algorithms. JHEP, 11:040, 2004.
 [35] Paolo Nason and Bryan Webber. NexttoLeadingOrder Event Generators. Ann. Rev. Nucl. Part. Sci., 62:187–213, 2012.
 [36] Tobias Neumann and Ciaran Williams. The Higgs boson at high . Phys. Rev., D95(1):014004, 2017.
 [37] Giovanni Ossola, Costas G. Papadopoulos, and Roberto Pittau. CutTools: A Program implementing the OPP reduction method to compute oneloop amplitudes. JHEP, 03:042, 2008.
 [38] T. Plehn, M. Spira, and P. M. Zerwas. Pair production of neutral Higgs particles in gluongluon collisions. Nucl. Phys., B479:46–64, 1996. [Erratum: Nucl. Phys.B531,655(1998)].
 [39] S. Schumann and F. Krauss. A parton shower algorithm based on cataniseymour dipole factorisation. Journal of High Energy Physics, 0803:038, 2008.
 [40] Torbjorn Sjostrand, Stefan Ask, Jesper R. Christiansen, Richard Corke, Nishita Desai, Philip Ilten, Stephen Mrenna, Stefan Prestel, Christine O. Rasmussen, and Peter Z. Skands. An Introduction to PYTHIA 8.2. Comput. Phys. Commun., 191:159–177, 2015.
 [41] Torbjorn Sjostrand, Stephen Mrenna, and Peter Z. Skands. A Brief Introduction to PYTHIA 8.1. Comput. Phys. Commun., 178:852–867, 2008.
 [42] Höche Stefan and Stefan Prestel. The midpoint between dipole and parton showers. Eur. Phys. J., C75(9):461, 2015.
 [43] A. van Hameren. OneLOop: For the evaluation of oneloop scalar functions. Comput. Phys. Commun., 182:2427–2438, 2011.
 [44] A. van Hameren, C. G. Papadopoulos, and R. Pittau. Automated oneloop calculations: A Proof of concept. JHEP, 09:106, 2009.