# Top-quark mass effects in double and triple Higgs production in gluon-gluon fusion at NLO

## Abstract:

The observation of double and triple scalar boson production at hadron colliders could provide key information on the Higgs self couplings and the potential. As for single Higgs production the largest rates for multiple Higgs production come from gluon-gluon fusion processes mediated by a top-quark loop. However, at variance with single Higgs production, top-quark mass and width effects from the loops cannot be neglected. Computations including the exact top-quark mass dependence are only available at the leading order, and currently predictions at higher orders are obtained by means of approximations based on the Higgs-gluon effective field theory (HEFT). In this work we present a reweighting technique that, starting from events obtained via the MCNLO method in the HEFT, allows to exactly include the top-quark mass and width effects coming from one- and two-loop amplitudes. We describe our approach and apply it to double Higgs production at NLO in QCD, computing the needed one-loop amplitudes and using approximations for the unknown two-loop ones. The results are compared to other approaches used in the literature, arguing that they provide more accurate predictions for distributions and for total rates as well. As a novel application of our procedure we present predictions at NLO in QCD for triple Higgs production at hadron colliders.

^{1}

## 1 Introduction

Present LHC data already provide convincing evidence that the scalar particle observed at the LHC is the one predicted by the Brout-Englert-Higgs breaking mechanism [1, 2] of the symmetry as implemented in the Standard Model (SM) [3]. Here, the strength of the Higgs boson couplings is uniquely determined by the masses of the elementary particles, including the Higgs boson itself. The measured couplings to fermions and vector bosons are found to agree within 10-20% with the SM predictions [4, 5]. No direct information, however, has been collected so far on the Higgs self-couplings that appear in the potential:

(1) |

The values of the Higgs self-couplings and are fixed in the SM by gauge invariance and renormalisability to , i.e. fully determined by the mass of the Higgs boson and the Higgs field vacuum expectation value . Direct information on the Higgs three-point and four-point interactions would therefore provide key information on the upper scale of validity of the SM when thought of as an effective theory itself, or on the possible existence of a richer scalar sector, featuring additional scalar fields, possibly in other representations.

In this context, multiple Higgs production plays a special role. At the lowest order, Higgs pair production is the simplest production process that is sensitive to the trilinear self-coupling , while to probe the quartic Higgs coupling one would need to consider at least triple Higgs production. Unfortunately, in the SM multiple Higgs production rates at the LHC are quite small [6, 7] and the prospects of making precise enough measurements at the LHC (assuming SM values) are at best challenging [8, 9, 10] for double Higgs production and rather bleak for triple Higgs production [7, 11].

However, multiple Higgs production rates could be enhanced by new physics effects and therefore the process provides one with a wealth of possibilities for probing new physics. These possibilities range from explicit models featuring new particles, such as the 2HDM [12, 13, 14, 15], SUSY [6, 16, 17, 18], extended colored sectors [19, 20, 21], Little Higgs Models [22, 23], Higgs portal [24, 25], flavor symmetry models [26] and Composite Higgs models [27, 28] to model independent higher-dimensional interactions [29, 30, 31]. In any case, precise predictions for rates and distributions will be needed to be able to extract valuable information on or on new physics effects in general. Higgs pair production will be relevant for the high-luminosity LHC and several phenomenological studies [32, 33, 34, 35, 36, 37, 38, 39, 40, 41] have recently investigated the potential of detecting the signal for the process over the backgrounds in various Higgs decay channels. Triple Higgs production rates are smaller by more than two orders of magnitude [7, 11], making the process possibly relevant for a future 100 TeV collider.

Analogously to single Higgs production, several channels can lead to a final state involving two or three Higgs bosons. They involve either the Higgs coupling to the top quark (e.g., gluon-gluon fusion via top loops and associated production), or to vector bosons (e.g., vector boson fusion and associated production), or to both (e.g., single-top associated production). The dominant production mechanism for multiple Higgs production is gluon-gluon fusion via top-quark loops, exactly as in the case of single Higgs production. Cross sections corresponding to the other channels are at least one order of magnitude smaller, even though possibly interesting because of different sensitivity to the ’s or to other couplings, such as [42, 43], or to new physics. In addition, typically, different channels open different possibilities of exploiting a wider range of Higgs decay signatures. Predictions for the six main production channels at NLO in QCD for hadron collider energies ranging from 8 to 100 TeV can be found in [44]. NNLO results for vector boson fusion [45] and vector boson associated production [33] are also available. For triple Higgs production gluon-gluon fusion is again the dominant production channel by at least one order of magnitude. However, the importance of the subdominant production channels is slightly different from the pair production at 14 TeV, with associated production giving the second largest contribution, and VBF the third.

Gluon-gluon fusion being the dominant production mechanism for multiple Higgs-boson production, one would like to have the most accurate and precise predictions possible for this channel. However, as Higgs pair production is a loop induced process at the Born level, higher-order computations become very involved. Computations including the exact top-quark mass dependence, i.e. in the Full Theory (FT), are only available at the leading order, and currently predictions at NLO [46] and NNLO accuracy [47, 48] are obtained by means of approximations that build upon the Higgs-gluon effective field theory. Besides, at variance with the inclusive cross section computation for single Higgs production, top-quark mass effects in the loops cannot be neglected. In this work we present a method that, starting from an NLO calculation matched to parton shower (PS) in the Higgs effective-field theory allows us to systematically include the heavy-quark mass effects coming from the one- and two-loop amplitudes, as long as results for the latter are available. We describe our approach, validate it against single Higgs production, for which all the necessary loop amplitudes are known, and apply it to double and triple Higgs production using all available information, i.e. the exact one-loop (real) amplitudes and the known parts of the two-loop amplitudes. We dub our results FT at NLO(+PS) and compare them to other approximations in the literature. We show that with just the currently available information on the loop amplitudes, our reweighting method already provides more accurate results for the differential distributions at NLO plus parton shower level and argue that that is the case for total rates as well.

The paper is structured as follows. In section 2 we discuss the current status of the computation for multiple Higgs production. In section 3 we describe our method and in particular how the top-quark mass effects can be consistently included in the calculation. In section 4 we present our results for Higgs pair production, comparing our method with previous results in the literature. In section 5 for the first time we present results for triple Higgs production beyond the leading order, obtained using the same method. We summarise our findings and draw some conclusions in the final section.

## 2 Multiple Higgs production in gluon-gluon fusion

In this section we briefly summarise the state-of-the-art in the computation of multiple Higgs production in gluon-gluon fusion. For the sake of conciseness we will focus only on Higgs pair production, most of the features being easily extendable to the general case. In the SM, the diagrams contributing to the gluon-gluon fusion channel can be organised in two classes: those where both Higgs bosons couple to the heavy quarks in the loop and those that feature the Higgs self coupling. The corresponding classes of diagrams appearing at the leading order are shown in fig. 1. At NLO, real one-loop diagrams, fig. 2 a), and virtual two-loop diagrams, both boxes and triangles, fig. 2 b), contribute, each of them being infrared divergent yet overall giving a finite result when combined in the computation of collinear and infrared safe observables. The one-loop matrix elements entering the Born amplitude as well as the real corrections can be obtained using modern loop techniques [49] automatically [50]. The two-loop triangles (such as the first one in fig. 2 b) ) featuring the Higgs self coupling are the same as those entering the single Higgs production and therefore also known for a long time [51, 52, 53] and used in publicly available codes for single Higgs production [54, 55]. The two-loop box amplitudes, however, are not known. This computation is extremely challenging due to the presence of several scales () in the loops and is currently out of reach. As it will become clear in the coming section, our reweighting technique can efficiently provide the NLO result with the full top-quark mass dependence once the corresponding amplitudes become known.

In the meantime, precisely due to the difficulty of including higher-order corrections exactly, the strategy widely employed for single Higgs-boson production has been adopted for Higgs pair production. An effective field theory, where the top quark has been integrated out from the theory and the Higgs boson couples directly to the gluon field, has been introduced, where the corresponding lagrangian reads

(2) |

being the QCD field tensor. The main motivation for using this approximation is that it makes the computation of higher-order corrections feasible. The approximation has been proven to work extremely well for single Higgs production [56]. The HEFT provides accurate predictions for the total rates as well as for the differential distributions when the invariants involved are not much larger than the top quark mass. Unfortunately, in the case of double Higgs production, the relevant scale is at least the invariant mass of the pair which is typically and therefore the HEFT provides only a rough approximation for the total rates and a very poor one for the relevant distributions [19, 34].

Given the fact that the full NLO results are not presently available and that the HEFT gives a poor description of the process, efforts have been made to improve results taking into account heavy-quark loop effects at least in an approximated way. A first step in this direction has been taken in the seminal NLO calculation for Higgs pair production, as implemented in the code HPAIR [6, 46], which provides total cross sections in the SM and in SUSY. In this case, the NLO calculation is performed within the HEFT, yet all contributions (virtual and real) to the short-distance parton-parton cross section are expressed in terms of the LO cross section times an correction. The LO cross section in the HEFT is then substituted by the LO one with the full heavy-quark mass dependence. As the aim of HPAIR is to provide results for the total cross section, this approximation is certainly a first useful step. However, such an approach is obviously not suitable in general and in particular for observables that receive considerable contributions from the real configurations, the most glaring example being the tail of the transverse momentum of the Higgs pair, .

We also mention that the top-quark mass effects at NLO in QCD have been estimated through an expansion [57] (on which we will further comment later). The NNLO computation for total rates in the HEFT [47] has also now been completed with the calculation of the three-loop matching coefficient [48], while recently results have been obtained by merging samples of different parton multiplicities in [37, 58] and including threshold resummation [59].

## 3 The inclusion of heavy-quark loop effects

In this section we discuss how to improve further on the HEFT approximation,
and in particular how one can improve on pure NLO in HEFT by including
heavy-quark loop effects in the Born amplitude as well as in the real emission ones. Results using this method have already
been presented in [44]. They have been achieved by implementing all the ingredients
in an automatically generated code for the HEFT within the MadGraph5_aMC@NLO
framework [60]. In practice, the effective Lagrangian
of eq. (2) is implemented in FeynRules [61, 62]^{2}

To include the heavy-quark amplitudes we make use of the structure of an NLO calculation as performed within MadGraph5_aMC@NLO. For completeness, let us consider the computation of a cross section at NLO, namely

(3) |

where are respectively the Born, virtual and real emission contributions, are the local counterterms (needed in order to render the integral over finite) and is the integrated form of (over the extra parton phase space). The detailed form of the counterterms depends on the subtraction scheme in use for the computation (e.g., FKS[72] or CS[73]). In general they involve a Born matrix element times some process-independent splitting kernel together with a dedicated phase-space mapping.

A very similar formulation as the one above holds in the case of matching NLO computations with parton-shower in the MC@NLO formalism, where on top of the local counterterms one has to also include the so-called Monte-Carlo counterterms [74]. The important difference with respect to the NLO formulation of eq. (3) is that the MC counterterms are such that Born-like (-events) and real-emission (-events) unweighted events can obtained as the corresponding subtracted cross sections are separately finite. The corresponding contributions to the total cross section can be written as

(4) | |||||

(5) |

In the MadGraph5_aMC@NLO framework, one can automatically generate the code corresponding to the Born, virtual, real amplitudes, the counter terms and the phase space [75, 50] in one go in order to compute cross sections and generate events for at NLO in QCD in the HEFT. All the finite heavy-quark one-loop matrix-elements (i.e. those entering the Born and real contributions) needed can also be obtained within MadGraph5_aMC@NLO. Note, however, that two limitations presently make the automatic computation of the exact NLO result not possible. First, the computation of cross sections that have a loop Born matrix-element is not automated yet (even at the LO only). Second, even with the automation for loop-induced processes, the need for the two-loop amplitudes would require an external routine, as this cannot be performed automatically by MadLoop. Therefore, the inclusion of heavy-quark effects needs manipulation that can in principle be performed in two ways.

The first option is to generate the code for an NLO computation in the HEFT and then replace the matrix-elements (for and ) with the corresponding ones in the FT. Even though this is the simplest option, it features several drawbacks. First, this method is very inefficient as the (computationally expensive) one-loop and two-loop matrix elements routines would then be called many times to probe and map all regions of phase space. In addition, it requires the evaluation of the real one-loop matrix elements in the FT in regions of phase space very close to the soft/collinear limits, i.e. where they might feature unstable configurations. For such points, multiple precision needs to be employed at the cost of a growth of the running time by a factor of a hundred.

The second option is to include the top-quark mass effects by reweighting after having generated the short-distance events and before these are passed to a parton shower program. In order for this procedure to be applied, all the weights corresponding to the separate contributions (events and counter events) and the corresponding kinematics, which is in general different between events and each of the counter events, need to be saved in an intermediate event file. With this information it is then possible to recompute the total event weight by reweighting each contribution by the matrix-elements in the FT. The weights corresponding to are rescaled by the ratio , while those corresponding to by the ratio . When unweighted events are generated, this amounts into rescaling the whole weight of -events with Born matrix-elements, and the different terms corresponding to -events as written above. This solution has the advantage of requiring the FT matrix-elements to be evaluated in significantly fewer phase space points than those used while integrating it directly. In addition, it is completely general and only assumes that there are no regions in phase space where the HEFT gives a vanishing contribution while the full theory does not. In our case this condition is satisfied. Moreover, being non-renormalisable, the HEFT leads to distributions that are in general harder in the high-energy tails than the FT. Therefore using the HEFT to generate events efficiently populates phase space regions that are later suppressed by the FT matrix elements, implying no large degradation of statistical accuracy.

In this work we follow the latter strategy, which we find efficient and general. Such technique has been proposed and employed to include heavy-quark corrections at LO in multi-parton merged samples for single Higgs [76] and [37] production. The implementation of the reweighting is simplified considerably by the fact that MadGraph5_aMC@NLO already features the needed structure and automatically stores all the relevant weights (though in a slightly different form than the one presented in eq. (3)), in order to evaluate the uncertainties associated with scale variation and PDF on an event-by-event basis as described in [77].

If the two-loop matrix elements were known, our procedure would have allowed to exactly compute cross sections at NLO accuracy and perform event generation at NLO+PS in the FT. As the finite part of the two-loop box virtual correction is unknown, this is not yet possible. In our work, we therefore employ and study the effect of using different approximations for the only unknown terms, while including all one-loop (Born and real) known terms. Once the two-loop results for the virtual corrections become available, it will be possible to include them straightforwardly, for example following the conventions of the Binoth Les Houches Accord [78].

The method described here can be applied to other loop-induced processes, for which suitably defined effective interactions can be used to automatically generate events while the FT one–loop matrix elements can then be employed to correctly reweight the events. As already mentioned, our method has been already used in the NLO calculation of Higgs pair production within the SM in [44] and more recently in the 2HDM in [13]. In this work, we will also apply the same procedure to triple Higgs production. We stress that it could also be applied to other loop-induced processes, such jet and associated production in gluon-gluon fusion.

## 4 Higgs pair production

In this section we present the results obtained by applying the method described in the previous section to production. The first effects we wish to discuss are those coming from the inclusion of the top-quark width. The relevance of the top-quark width in the context of single Higgs production in gluon-gluon fusion has been discussed in [54] where it has been verified that for a light Higgs below the threshold the effect is negligible but rises to the percent level for heavier Higgs bosons, i.e., closer to the threshold. Such virtualities are exactly those relevant for Higgs pair production and one can therefore expect the top-quark width effects not to be negligible.

As already mentioned, in our calculation the and one-loop amplitudes are obtained via MadLoop [50]. The computation is performed using the complex mass scheme for the top quark [79, 80] as implemented in MadLoop and ALOHA in MadGraph5_aMC@NLO. In practice, for the top quark this simply amounts to performing the replacement

(6) |

everywhere in the computation, i.e. to modify the kinematical mass as well as the Yukawa coupling.
The effect of including a non-zero top-quark width is shown in fig. 3, where the LO matrix
element squared for is plotted as a function of the invariant mass of the Higgs pair for and 1.5 GeV.^{3}

We now consider the inclusion of the finite top mass in the NLO computation. In what we dub NLO FT, the Born and real configurations are reweighted with the corresponding Born and real emission finite top-quark mass matrix elements and for the virtual configurations, the HEFT result, yet rescaled by the Born in the FT, is used. We stress again that the only approximation made in this procedure is that coming from the absence of the exact results for the two-loop virtual terms. As a check we have applied this method to single Higgs production in gluon-gluon fusion where all elements of the calculation are available and found excellent agreement with the corresponding FT implementation in MC@NLO v4.0 [81].

We also compare our NLO FT results to two other approximations. The first of these two approximations corresponds to what in the following we refer to as “Born-improved” HEFT results and amounts to rescaling all events, Born and real-emission like, by the Born in the FT, i.e., to always use the weight . In the case of the real-emission events the kinematics used is that determined by the counter terms mapping. This approach mimics the one used in HPAIR, the only (practically irrelevant) difference being that in HPAIR the reweighting is performed after the integration over the polar angle, while our HEFT Born reweighting is fully differential. Indeed a detailed check with the results of HPAIR shows that once the same parameters are chosen, and in particular the top-width is set to zero, an excellent agreement is found.

In addition to the Born-improved HEFT results, we compare our reference results to a second approximation,
dubbed NLO FT, by which we assess
the possible shifts of the central value of the NLO cross section due to the
inclusion of the exact two-loop corrections in the FT. As already mentioned several times, the two-loop box contributions, e.g., those coming from the second and third diagrams of fig. 2b) are not known. However, the two-loop triangle contributions, e.g., first diagram of fig. 2b), are known (and form a gauge independent subset) as they are the same as those entering single Higgs production. To this aim, we use the finite contribution of the two–loop corrections for single Higgs production, which are publicly available as a function of the Higgs mass, e.g. in SusHi [55],
and assume that the same form factor would hold for the two-loop box amplitude. In other words, we proceed under the assumption that these corrections factorise out of the matrix element globally, i.e. for both the box and the triangle diagrams. Needless to say, such an assumption is *ad hoc* and while we do not claim that the triangle form-factor should resemble the box one, this test is still a useful one. The main reason comes from the fact that it allows us to study the extent of possible cancellations taking place between the top-quark mass corrections from the real (which are included) and virtual (which are only approximate) matrix elements. It can be observed^{4}

We quantify the differences between the approximations discussed above by calculating the total cross section for Higgs pair production in gluon-gluon fusion at 14 TeV in tab. 1. In this computation, we have set the Higgs mass to GeV and the top-quark mass to GeV. We note that the top-quark mass dependence of the cross section around the reference top mass of 172.5 GeV can be parametrised approximately (at LO FT as well as at NLO FT) by GeV. Parton distribution functions (PDFs) are evaluated by using the MSTW2008 (LO and NLO) parametrisation in the five-flavour scheme [84]. The renormalisation and factorisation scales are set to . The dependence of the predictions on scale and PDF variations can be estimated at no extra computational cost via a reweighting technique [77]. Scales are varied independently in the range and PDF uncertainties at the 68% C.L. are obtained following the prescription given by the MSTW collaboration [84]. Even though -quark loops can be computed in our setup, -quark masses as well as their tiny (0.3%) contribution to the cross section are neglected in the following.

Table 1 collects our results. We first verify that the effect of the non–zero top-quark width on the total cross section at LO, a % decrease, directly follows from the results shown in fig. 3 and the fact that the invariant mass distribution peaks at 400 GeV. We also note the well-known fact that the process receives large QCD corrections as well as the expected reduction of the theoretical uncertainties for the NLO computations. We then show three NLO results: i) the Born-improved HEFT result through a local event-by-event reweighting, ii) the NLO FT result, obtained by combining the exact real emission matrix elements, with the Born-rescaled HEFT results for the virtual corrections and iii) the NLO FT result obtained by combining the exact real emission matrix elements, with the exact results of single Higgs production for the virtual corrections, as described previously. For all NLO results we keep the finite top-quark width of 1.5 GeV.

production in gluon-gluon fusion at 14 TeV | Cross section [fb] | |
---|---|---|

HEFT | 19.2 | |

LO | FT, GeV | 23.2 |

FT, GeV | 22.7 | |

NLO | HEFT | 32.9 |

HEFT Born-improved | 38.5 | |

FT (virtuals: Born-rescaled HEFT ) | 34.3 | |

FT (virtuals: estimated from single Higgs in FT) | 35.0 |

We can now compare the different approximations of the FT NLO result. The first important observation is that the Born-improved result is 11% larger than our baseline one. We also note here that the Born-improved result obtained by a local event-by-event rescaling is within 1% of what one would obtain from a global Born rescaling obtained using the total cross section numbers, i.e., . The difference from the Born-improved result only slightly reduces (9%) when an estimate for the finite top-quark mass terms from the two-loop contributions is included, see last line of tab. 1. Our NLO FT result is rather stable in that respect. This is related to the fact that the cancellation we discussed earlier for single Higgs production is only relevant very close to the threshold, with the Born-rescaled result rapidly rising over the exact one above the threshold. In the case of single Higgs production, we have indeed checked that for Higgs masses above 400 GeV the NLO FT result (only including the exact real emission matrix element but not the known two–loop virtual results) is closer to the exact result than the corresponding Born-improved one. In the case of Higgs pair production, one could also argue that even if a similar cancellation of the top-quark mass effects between the real and virtual corrections occurred at the threshold, it would not have a very pronounced effect on the total cross section, as for Higgs pair production the peak of the invariant mass distribution is located at higher mass values.

At this point it is worth to recall the results of ref. [57], where the top-quark mass effects at NLO in QCD were estimated by computing the first few terms in the expansion for the factor. The expansion is known not to converge well at LO [19] and is not supposed to work beyond or even close to the threshold, around and beyond which the bulk of the cross section resides. However, in ref. [57] an attempt was made by combining the exact Born cross section with the expanded factors, as a “taming” technique for the expansion. A +10% increase with respect to the Born-rescaled HEFT result was found, i.e., an effect similar in size but opposite in sign to our estimate. Combined with our findings, the estimate of ref. [57] implies that the difference between the finite part of the Born-rescaled HEFT virtuals and the exact ones should account for a +20% increase of the total cross section, a quite large effect indeed, especially considering that by including top-mass effects in the virtual corrections estimated via the known two-loop triangles, leads only to a couple of percent increase. Besides, we note that the results of the same expansion approach applied to the production of a single heavy Higgs of mass between 400 and 500 GeV, are known to overestimate the exact results in the FT when no high-energy matching is performed [85, 55, 86].

While only an exact calculation of the missing two-loop amplitudes will finally settle this issue, the NLO FT approach provides central values for the cross sections that appear rather robust, predicting a correction of about -10% with respect to those obtained by means of the Born-improved HEFT. In addition, together with the results of ref. [57], our study provides an estimate of about 10% for the uncertainty to be associated with the HEFT calculation due to the missing top-quark mass effects. Such an uncertainty should be quoted along with the other theoretical uncertainties in the HEFT calculations, at NLO but also at NNLO.

Finally, we note that including the exact one–loop matrix elements provides a more accurate description of the tails of the distributions where hard parton emissions take place, and the factorisation of the real–emission amplitudes into the LO amplitudes, as implicit in the Born-improved approach, cannot accurately describe the hard parton kinematics. To emphasize this point, we compare in fig. 4 the differential distributions of the invariant mass and the transverse momentum of the Higgs pair obtained by NLO FT and Born-improved HEFT. Pythia8 [87] has been used for parton shower and hadronisation.

By studying the two distributions, we see that including the exact real emission matrix elements affects the two observables in different ways. On the one hand, for the invariant mass distribution the effect is a uniform modification of the cross section by about 10%. On the other hand, for the transverse momentum of the Higgs pair, at low values where the distribution is dominated by the parton shower, there is no visible difference between the two results, while at high values where the distribution is dominated by hard parton emission, coming from the real matrix elements, we see that there is a significant deviation from the Born-improved result. In that region we trust our FT prediction more as it describes better the kinematics.

[fb] | TeV | TeV | TeV |
---|---|---|---|

LO FT | 22.7 | 145 | 1080 |

NLO FT | 34.3 | 199 | 1250 |

For completeness, we also show the results for the total cross section at LO and NLO for Higgs pair production in gluon-gluon fusion as a function of the hadronic centre-of-mass energy in fig. 5, including the uncertainties due to scale variation and PDF added linearly. The uncertainty bands demonstrate nicely the reduction of the scale uncertainties for the NLO calculation. Results for a selected number of hadronic energies are also shown in tab. 2.

## 5 Higgs triple production

We now apply the same procedure followed for Higgs pair production to obtain results beyond the LO for triple Higgs production through gluon-gluon fusion. Representative topologies of diagrams that contribute to the process at LO are shown in fig. 6. These include pentagons where the Higgs bosons couple only to the heavy quarks of the loop, boxes in combination with one power of the trilinear Higgs coupling and triangles with either two powers of the trilinear coupling or the quartic Higgs coupling. It turns out that at LO, the hierarchy of the various contributions in decreasing size order is first pentagons, then boxes and finally triangles [7]. With the triangles contribution already being the smallest of the three, and the fact that the quartic coupling itself has only a very small contribution to that (it only appears in one of the four diagrams), any attempt to access the quartic coupling through the measurement of this process becomes not feasible at the LHC and extremely challenging even at a future 100 TeV collider. In any case, measuring the quartic Higgs coupling in this process would require a very precise knowledge of the triple Higgs coupling [7]. The situation could possibly change in the presence of physics beyond the SM. In this respect having predictions as accurate as possible for production cross sections at hadron colliders in the SM is still be very useful.

We therefore propose to follow the same procedure as for Higgs pair production and improve the NLO HEFT results by systematically including the information from the FT one-loop amplitudes, the Born and the real contributions. In this case, the information which will be included in an approximated way is that related to the two-loop boxes and two-loop pentagons, whose calculation is extremely challenging (note that even the two-loop boxes are more complicated than those in production as they feature one more scale) and probably will not be available for some time. The method and the setup follows exactly that of , even though the calculation is more complicated and the resulting reweighting procedure is substantially slower. For this process, we find it necessary to use a small computing cluster and fully parallelise the reweighting on an event-by-event basis.

We show our results for the production cross sections as a function of the centre-of-mass energy in fig. 7 and a few representative results in tab. 3. The NLO calculation leads to factors and uncertainties which are similar to Higgs pair production. The most important information conveyed by fig. 7 is that the cross section remains very small even at 100 TeV collisions.

We conclude by showing the NLO+PS effects in two key distributions, i.e. that of the invariant mass of the three Higgs-bosons and the transverse momentum of the triplet. The latter distribution features two important damping effects: at small the spectrum is softened by the soft resummation performed by the shower and at high where the top-quark loop effects matter and the HEFT is not reliable.

[fb] | TeV | TeV | TeV |
---|---|---|---|

LO FT | 0.0557 | 0.438 | 3.78 |

NLO FT | 0.0894 | 0.677 | 5.09 |

## 6 Conclusions

The observation of multiple Higgs production at hadron colliders is a very challenging task, yet a crucial one to obtain key information on the form of the Higgs potential. Rates for these processes are rather low and the accessible signatures swamped by large backgrounds. In any case, any effort to gather information from measurements or bounds on these processes requires accurate predictions for the SM total cross sections. In addition, differential distributions are needed not only to calculate acceptances but also to improve the potential of disentangling these processes from the various backgrounds by selecting the most sensitive regions in phase space. As in single Higgs production, the largest rates for multiple Higgs production come from gluon-gluon fusion mediated by a top-quark loop. Loop-induced processes pose one with the difficulty of obtaining higher order predictions, as these require the computation of multi–loop Feynman diagrams.

In this work we have first focused on Higgs pair production by considering different approximations to improve the HEFT NLO calculation upon the infinite top-quark mass limit. We have introduced a reweighting procedure that allows the inclusion of the top-quark mass and width effects exactly. We have then applied it to production using the available information, i.e. the exact (Born and real) one-loop amplitudes and the approximated two-loop matrix elements to appropriately reweight events generated automatically by means of the MCNLO method in the effective field theory. As a first result and at variance with single Higgs production, we have found that including a non-zero top-quark width reduces the cross section by a couple of percents, the largest effect reaching -4% just above the threshold.

We have then performed a study to assess the relevance of various corrections and the accuracy of other approaches used in the literature to approximate NLO results in the FT. In particular we have compared to a Born reweighting, where only the exact Born results are used to improve upon the HEFT results. At the total cross section level our best estimate gives a result about 10% smaller than the Born-improved HEFT. We have then considered the difference between the two approaches for differential distributions, and found that including the exact real emission matrix elements provides a better description of the phase space regions where hard emissions take place. We have then argued that total rates are improved too: by using an estimated form of the unknown virtual corrections in the FT using available results from single Higgs production, we have shown that our results are rather stable under variations of the unknown finite terms. Even though the effect of the missing two-loop virtual corrections on the total cross section cannot be quantified until these become available, comparing the different approximations allows one to conservatively associate an uncertainty of order 10% with the calculation due to the missing top-quark mass effects. Note that this uncertainty should be quoted along with others until the exact NLO result becomes available. Finally, we have applied our method to triple Higgs production, providing for the first time predictions for hadron colliders at NLO (+PS) accuracy in the SM. We have found a reduction of the theoretical uncertainties and enhancements of the cross section similar to those of production over a large range of collision energies.

## Acknowledgements

We are grateful to Claude Duhr, Giampiero Passarino, and Michael Spira for stimulating discussions, and to Marius Wiesemann, Tobias Neumann and Michael Rauch for their collaboration. We thank Stefano Frixione for many discussions and for comments on the manuscript. This work has been supported in part by the Research Executive Agency (REA) of the European Union under the Grant Agreement number PITN-GA-2010-264564 (LHCPhenoNet) and PITN-GA-2012-315877 (MCNet). The work of MZ is partially supported by the ILP LABEX (ANR- 10-LABX-63), in turn supported by French state funds managed by the ANR within the “Investissements d’Avenir” programme under reference ANR-11-IDEX-0004-02.

### Footnotes

- preprint: CP3-14-61, MCnet-14-18, LPN14-105
- More complete implementations including interactions coming from the full set of dimension-six operators are also available [63, 64, 65, 66].
- We note that here we assumed a 90 scattering for all points included in fig. 3, but as the matrix element has an extremely weak angular dependence [19], this provides a perfectly good demonstration of the effect also at the level of the Higgs pair invariant mass distribution.
- We are grateful to Michael Spira for bringing up this point.

### References

- F. Englert and R. Brout, Broken Symmetry and the Mass of Gauge Vector Mesons, Phys.Rev.Lett. 13 (1964) 321–323.
- P. W. Higgs, Broken Symmetries and the Masses of Gauge Bosons, Phys.Rev.Lett. 13 (1964) 508–509.
- S. Weinberg, Implications of Dynamical Symmetry Breaking, Phys.Rev. D13 (1976) 974–996.
- CMS-HIG-13-003. CMS-HIG-13-004. CMS-HIG-13-006. CMS-HIG-13-009. (2013).
- ATLAS-CONF-2013-009. ATLAS-CONF-2013-010. ATLAS-CONF-2013-012. ATLAS- CONF-2013-013. (2013).
- T. Plehn, M. Spira, and P. Zerwas, Pair production of neutral Higgs particles in gluon-gluon collisions, Nucl.Phys. B479 (1996) 46–64, [hep-ph/9603205].
- T. Plehn and M. Rauch, The quartic higgs coupling at hadron colliders, Phys.Rev. D72 (2005) 053008, [hep-ph/0507321].
- U. Baur, T. Plehn, and D. L. Rainwater, Measuring the Higgs boson self coupling at the LHC and finite top mass matrix elements, Phys.Rev.Lett. 89 (2002) 151801, [hep-ph/0206024].
- U. Baur, T. Plehn, and D. L. Rainwater, Examining the Higgs boson potential at lepton and hadron colliders: A Comparative analysis, Phys.Rev. D68 (2003) 033001, [hep-ph/0304015].
- U. Baur, T. Plehn, and D. L. Rainwater, Probing the Higgs selfcoupling at hadron colliders using rare decays, Phys.Rev. D69 (2004) 053004, [hep-ph/0310056].
- T. Binoth, S. Karg, N. Kauer, and R. Ruckl, Multi-Higgs boson production in the Standard Model and beyond, Phys.Rev. D74 (2006) 113008, [hep-ph/0608057].
- J. Baglio, O. Eberhardt, U. Nierste, and M. Wiebusch, Benchmarks for Higgs Pair Production and Heavy Higgs Searches in the Two-Higgs-Doublet Model of Type II, Phys.Rev. D90 (2014) 015008, [arXiv:1403.1264].
- B. Hespel, D. Lopez-Val, and E. Vryonidou, Higgs pair production via gluon fusion in the Two-Higgs-Doublet Model, arXiv:1407.0281.
- E. Asakawa, D. Harada, S. Kanemura, Y. Okada, and K. Tsumura, Higgs boson pair production in new physics models at hadron, lepton, and photon colliders, Phys.Rev. D82 (2010) 115002, [arXiv:1009.4670].
- A. Arhrib, R. Benbrik, C.-H. Chen, R. Guedes, and R. Santos, Double Neutral Higgs production in the Two-Higgs doublet model at the LHC, JHEP 0908 (2009) 035, [arXiv:0906.0387].
- J. Cao, Z. Heng, L. Shang, P. Wan, and J. M. Yang, Pair Production of a 125 GeV Higgs Boson in MSSM and NMSSM at the LHC, JHEP 1304 (2013) 134, [arXiv:1301.6437].
- U. Ellwanger, Higgs pair production in the NMSSM at the LHC, JHEP 1308 (2013) 077, [arXiv:1306.5541].
- D. T. Nhung, M. Muhlleitner, J. Streicher, and K. Walz, Higher Order Corrections to the Trilinear Higgs Self-Couplings in the Real NMSSM, JHEP 1311 (2013) 181, [arXiv:1306.3926].
- S. Dawson, E. Furlan, and I. Lewis, Unravelling an extended quark sector through multiple Higgs production?, Phys.Rev. D87 (2013) 014007, [arXiv:1210.6663].
- G. D. Kribs and A. Martin, Enhanced di-Higgs Production through Light Colored Scalars, Phys.Rev. D86 (2012) 095023, [arXiv:1207.4496].
- C.-Y. Chen, S. Dawson, and I. Lewis, Top Partners and Higgs Production, arXiv:1406.3349.
- C. O. Dib, R. Rosenfeld, and A. Zerwekh, Double Higgs production and quadratic divergence cancellation in little Higgs models with T parity, JHEP 0605 (2006) 074, [hep-ph/0509179].
- L. Wang, W. Wang, J. M. Yang, and H. Zhang, Higgs-pair production in littlest Higgs model with T-parity, Phys.Rev. D76 (2007) 017702, [arXiv:0705.3392].
- M. J. Dolan, C. Englert, and M. Spannowsky, New Physics in LHC Higgs boson pair production, Phys.Rev. D87 (2013), no. 5 055002, [arXiv:1210.8166].
- J. M. No and M. Ramsey-Musolf, Probing the Higgs Portal at the LHC Through Resonant di-Higgs Production, Phys.Rev. D89 (2014) 095031, [arXiv:1310.6035].
- E. L. Berger, S. B. Giddings, H. Wang, and H. Zhang, Higgs-flavon mixing and LHC phenomenology in a simplified model of broken flavor symmetry, arXiv:1406.6054.
- R. Grober and M. Muhlleitner, Composite Higgs Boson Pair Production at the LHC, JHEP 1106 (2011) 020, [arXiv:1012.1562].
- M. Gillioz, R. Grober, C. Grojean, M. Muhlleitner, and E. Salvioni, Higgs Low-Energy Theorem (and its corrections) in Composite Models, JHEP 1210 (2012) 004, [arXiv:1206.7120].
- R. Contino, M. Ghezzi, M. Moretti, G. Panico, F. Piccinini, et al., Anomalous Couplings in Double Higgs Production, JHEP 1208 (2012) 154, [arXiv:1205.5444].
- A. Pierce, J. Thaler, and L.-T. Wang, Disentangling Dimension Six Operators through Di-Higgs Boson Production, JHEP 0705 (2007) 070, [hep-ph/0609049].
- J. Liu, X.-P. Wang, and S.-h. Zhu, Discovering extra Higgs boson via pair production of the SM-like Higgs bosons, arXiv:1310.3634.
- A. Papaefstathiou, L. L. Yang, and J. Zurita, Higgs boson pair production at the LHC in the channel, Phys.Rev. D87 (2013) 011301, [arXiv:1209.1489].
- J. Baglio, A. Djouadi, R. Gröber, M. Mühlleitner, J. Quevillon, et al., The measurement of the Higgs self-coupling at the LHC: theoretical status, JHEP 1304 (2013) 151, [arXiv:1212.5581].
- M. J. Dolan, C. Englert, and M. Spannowsky, Higgs self-coupling measurements at the LHC, JHEP 1210 (2012) 112, [arXiv:1206.5001].
- A. J. Barr, M. J. Dolan, C. Englert, and M. Spannowsky, Di-Higgs final states augMT2ed – selecting events at the high luminosity LHC, Phys.Lett. B728 (2014) 308–313, [arXiv:1309.6318].
- M. Gouzevitch, A. Oliveira, J. Rojo, R. Rosenfeld, G. P. Salam, et al., Scale-invariant resonance tagging in multijet events and new physics in Higgs pair production, JHEP 1307 (2013) 148, [arXiv:1303.6636].
- Q. Li, Q.-S. Yan, and X. Zhao, Higgs Pair Production: Improved Description by Matrix Element Matching, Phys.Rev. D89 (2014) 033015, [arXiv:1312.3830].
- F. Goertz, A. Papaefstathiou, L. L. Yang, and J. Zurita, Higgs Boson self-coupling measurements using ratios of cross sections, JHEP 1306 (2013) 016, [arXiv:1301.3492].
- D. E. Ferreira de Lima, A. Papaefstathiou, and M. Spannowsky, Standard model Higgs boson pair production in the ( )( ) final state, JHEP 1408 (2014) 030, [arXiv:1404.7139].
- V. Barger, L. L. Everett, C. Jackson, and G. Shaughnessy, Higgs-Pair Production and Measurement of the Triscalar Coupling at LHC(8,14), Phys.Lett. B728 (2014) 433–436, [arXiv:1311.2931].
- M. Slawinska, W. v. d. Wollenberg, B. van Eijk, and S. Bentvelsen, Phenomenology of the trilinear Higgs coupling at proton-proton colliders, arXiv:1408.5010.
- R. Contino, C. Grojean, D. Pappadopulo, R. Rattazzi, and A. Thamm, Strong Higgs Interactions at a Linear Collider, JHEP 1402 (2014) 006, [arXiv:1309.7038].
- D. Pappadopulo, A. Thamm, R. Torre, and A. Wulzer, Heavy Vector Triplets: Bridging Theory and Data, arXiv:1402.4431.
- R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, et al., Higgs pair production at the LHC with NLO and parton-shower effects, Phys.Lett. B732 (2014) 142–149, [arXiv:1401.7340].
- L. Liu-Sheng, Z. Ren-You, M. Wen-Gan, G. Lei, L. Wei-Hua, et al., NNLO QCD corrections to Higgs pair production via vector boson fusion at hadron colliders, Phys.Rev. D89 (2014) 073001, [arXiv:1401.7754].
- S. Dawson, S. Dittmaier, and M. Spira, Neutral Higgs boson pair production at hadron colliders: QCD corrections, Phys.Rev. D58 (1998) 115012, [hep-ph/9805244].
- D. de Florian and J. Mazzitelli, Higgs Boson Pair Production at Next-to-Next-to-Leading Order in QCD, Phys. Rev. Lett. 111, 201801 (2013) [arXiv:1309.6594].
- J. Grigo, K. Melnikov, and M. Steinhauser, Virtual corrections to Higgs boson pair production in the large top quark mass limit, arXiv:1408.2422.
- G. Ossola, C. G. Papadopoulos, and R. Pittau, CutTools: A Program implementing the OPP reduction method to compute one-loop amplitudes, JHEP 0803 (2008) 042, [arXiv:0711.3596].
- V. Hirschi et al., Automation of one-loop QCD corrections, JHEP 05 (2011) 044, [arXiv:1103.0621].
- D. Graudenz, M. Spira, and P. Zerwas, QCD corrections to Higgs boson production at proton proton colliders, Phys.Rev.Lett. 70 (1993) 1372–1375.
- M. Spira, A. Djouadi, D. Graudenz, and P. Zerwas, Higgs boson production at the LHC, Nucl.Phys. B453 (1995) 17–82, [hep-ph/9504378].
- R. Harlander and P. Kant, Higgs production and decay: Analytic results at next-to-leading order QCD, JHEP 0512 (2005) 015, [hep-ph/0509189].
- C. Anastasiou, S. Buehler, F. Herzog, and A. Lazopoulos, Total cross-section for Higgs boson hadroproduction with anomalous Standard Model interactions, JHEP 1112 (2011) 058, [arXiv:1107.0683].
- R. V. Harlander, S. Liebler, and H. Mantler, SusHi: A program for the calculation of Higgs production in gluon fusion and bottom-quark annihilation in the Standard Model and the MSSM, Computer Physics Communications 184 (2013) 1605–1617, [arXiv:1212.3249].
- R. V. Harlander and K. J. Ozeren, Finite top mass effects for hadronic Higgs production at next-to-next-to-leading order, JHEP 0911 (2009) 088, [arXiv:0909.3420].
- J. Grigo, J. Hoff, K. Melnikov, and M. Steinhauser, On the Higgs boson pair production at the LHC, Nucl.Phys. B875 (2013) 1–17, [arXiv:1305.7340].
- P. MaierhÃ¶fer and A. Papaefstathiou, Higgs Boson pair production merged to one jet, JHEP 1403 (2014) 126, [arXiv:1401.0007].
- D. Y. Shao, C. S. Li, H. T. Li, and J. Wang, Threshold resummation effects in Higgs boson pair production at the LHC, JHEP 1307 (2013) 169, [arXiv:1301.1245].
- J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, et al., The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 1407 (2014) 079, [arXiv:1405.0301].
- N. D. Christensen and C. Duhr, FeynRules - Feynman rules made easy, Comput.Phys.Commun. 180 (2009) 1614–1641, [arXiv:0806.4194].
- A. Alloul, N. D. Christensen, C. Degrande, C. Duhr, and B. Fuks, FeynRules 2.0 - A complete toolbox for tree-level phenomenology, Comput.Phys.Commun. 185 (2014) 2250–2300, [arXiv:1310.1921].
- P. Artoisenet, P. de Aquino, F. Demartin, R. Frederix, S. Frixione, et al., A framework for Higgs characterisation, JHEP 1311 (2013) 043, [arXiv:1306.6464].
- A. Alloul, B. Fuks, and V. Sanz, Phenomenology of the Higgs Effective Lagrangian via FEYNRULES, JHEP 1404 (2014) 110, [arXiv:1310.5150].
- F. Maltoni, K. Mawatari, and M. Zaro, Higgs characterisation via vector-boson fusion and associated production: NLO and parton-shower effects, Eur.Phys.J. 74 (2014) 2710, [arXiv:1311.1829].
- F. Demartin, F. Maltoni, K. Mawatari, B. Page, and M. Zaro, Higgs characterisation at NLO in QCD: CP properties of the top-quark Yukawa interaction, arXiv:1407.5089.
- G. Ossola, C. G. Papadopoulos, and R. Pittau, Reducing full one-loop amplitudes to scalar integrals at the integrand level, Nucl.Phys. B763 (2007) 147–169, [hep-ph/0609007].
- C. Degrande, Automatic evaluation of UV and R2 terms for beyond the Standard Model Lagrangians: a proof-of-principle, arXiv:1406.3030.
- B. Page and R. Pittau, vertices for the effective ggH theory, JHEP 1309 (2013) 078, [arXiv:1307.6142].
- C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, Mattelaer, et al., UFO - The Universal FeynRules Output, Comput.Phys.Commun. 183 (2012) 1201–1214, [arXiv:1108.2040].
- P. de Aquino, W. Link, F. Maltoni, O. Mattelaer, and T. Stelzer, ALOHA: Automatic Libraries Of Helicity Amplitudes for Feynman Diagram Computations, Comput.Phys.Commun. 183 (2012) 2254–2263, [arXiv:1108.2041].
- S. Frixione, Z. Kunszt, and A. Signer, Three jet cross-sections to next-to-leading order, Nucl.Phys. B467 (1996) 399–442, [hep-ph/9512328].
- S. Catani and M. Seymour, A General algorithm for calculating jet cross-sections in NLO QCD, Nucl.Phys. B485 (1997) 291–419, [hep-ph/9605323].
- S. Frixione and B. R. Webber, Matching NLO QCD computations and parton shower simulations, JHEP 06 (2002) 029, [hep-ph/0204244].
- R. Frederix, S. Frixione, F. Maltoni, and T. Stelzer, Automation of next-to-leading order computations in QCD: the FKS subtraction, JHEP 10 (2009) 003, [arXiv:0908.4272].
- J. Alwall, Q. Li, and F. Maltoni, Matched predictions for Higgs production via heavy-quark loops in the SM and beyond, Phys.Rev. D85 (2012) 014031, [arXiv:1110.1728].
- R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, R. Pittau, et al., Four-lepton production at hadron colliders: aMC@NLO predictions with theoretical uncertainties, JHEP 1202 (2012) 099, [arXiv:1110.4738].
- S. Alioli, S. Badger, J. Bellm, B. Biedermann, F. Boudjema, et al., Update of the Binoth Les Houches Accord for a standard interface between Monte Carlo tools and one-loop programs, Comput.Phys.Commun. 185 (2014) 560–571, [arXiv:1308.3462].
- A. Denner, S. Dittmaier, M. Roth, and D. Wackeroth, Predictions for all processes e+ e- 4 fermions + gamma, Nucl.Phys. B560 (1999) 33–65, [hep-ph/9904472].
- A. Denner, S. Dittmaier, M. Roth, and L. Wieders, Electroweak corrections to charged-current e+ e- 4 fermion processes: Technical details and further results, Nucl.Phys. B724 (2005) 247–294, [hep-ph/0505042].
- S. Frixione, F. Stoeckli, P. Torrielli, B. R. Webber, and C. D. White, The MCaNLO 4.0 Event Generator, arXiv:1010.0819.
- R. Harlander, Supersymmetric Higgs production at the Large Hadron Collider, Eur.Phys.J. C33 (2004) S454–S456, [hep-ph/0311005].
- C. Anastasiou, S. Bucherer, and Z. Kunszt, HPro: A NLO Monte-Carlo for Higgs production via gluon fusion with finite heavy quark masses, JHEP 0910 (2009) 068, [arXiv:0907.2362].
- A. Martin, W. Stirling, R. Thorne, and G. Watt, Parton distributions for the LHC, Eur.Phys.J. C63 (2009) 189–285, [arXiv:0901.0002].
- R. V. Harlander and W. B. Kilgore, Next-to-next-to-leading order Higgs production at hadron colliders, Phys.Rev.Lett. 88 (2002) 201801, [hep-ph/0201206].
- R. V. Harlander, H. Mantler, S. Marzani, and K. J. Ozeren, Higgs production in gluon fusion at next-to-next-to-leading order QCD for finite top mass, Eur.Phys.J. C66 (2010) 359–372, [arXiv:0912.2104].
- T. Sjostrand, S. Mrenna, and P. Z. Skands, A Brief Introduction to PYTHIA 8.1, Comput.Phys.Commun. 178 (2008) 852–867, [arXiv:0710.3820].