Higgs pair production at the LHC with NLO and parton-shower effects

We present predictions for the SM-Higgs-pair production channels of relevance at the LHC: gluon-gluon fusion, VBF, and top-pair, W, Z and single-top associated production. All these results are at the NLO accuracy in QCD, and matched to parton showers by means of the MC@NLO method; hence, they are fully differential. With the exception of the gluon-gluon fusion process, for which a special treatment is needed in order to improve upon the infinite-top-mass limit, our predictions are obtained in a fully automatic way within the publicly available MadGraph5_aMC@NLO framework. We show that for all channels in general, and for gluon-gluon fusion and top-pair associated production in particular, NLO corrections reduce the theoretical uncertainties, and are needed in order to arrive at reliable predictions for total rates as well as for distributions.


Introduction
Present LHC data provide evidence that the scalar particle observed at the LHC is the one predicted by the Brout-Englert-Higgs symmetry breaking mechanism [1,2] of SU (2) L × U (1) Y as implemented in the Standard Model (SM) [2]. In this case, the strengths of the Higgs boson couplings are uniquely determined by the masses of the elementary particles. The measured couplings to fermions and vector bosons agree within 10-20% with the SM predictions [3,4]. No information, however, has been collected so far on the Higgs self-coupling λ. In the SM the Higgs boson mass itself fixes the value of this self coupling in the scalar potential whose form, in turn, is determined by the global symmetries and the requirement of renormalisability. These conditions, however, have no raison d'être once experimental indications (such as the existence of dark matter) as well as theoretical arguments (such as naturalness) are put forward. In this respect, it is appropriate (and, in fact, advantageous) to consider the SM as the subset of operators of dimension less than or equal to four of an effective field theory (EFT) lagrangian with an SU (2) L ×U (1) Y symmetry. Direct information on the Higgs three-and four-point interactions could therefore provide a key indication of the structure of the scalar potential, and of where the scale Λ characterising such an EFT might lie.
In this context, Higgs pair production could play a key role. Not only it is the simplest production process that is sensitive to the self-coupling λ, but it also provides one with a wealth of possibilities for probing higher-dimensional interactions as well as the existence of heavier states coupled to the Higgs [5,6,7,8,9,10,11]. Unfortunately, in the SM the rates for Higgs pair production at the LHC are quite small [12,13,14,15]. So unless new physics produces sizable enhancements (something quite possible in several scenarios), a measurement of the HH production cross sections will necessitate considerable integrated luminosity even at 14 TeV centre-of-mass energy. In any case, precise predictions for rates and distributions will be needed in order to be able to extract valuable information on λ or on new physics effects in general.
Analogously to single-Higgs production, several channels can lead to a final state involving two Higgs bosons. They entail the Higgs coupling to either the top quark (as in the case of gluon-gluon fusion and of tt associated production), or vector bosons (in VBF, and in W and Z associated production), or both (for single-top associated production). The dominant production mechanism is gluon-gluon fusion via a top loop, 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 λ or to new physics, and because of the possibility of exploiting a wider range of Higgs decay signatures.
In this Letter we present results accurate to NLO in QCD for the six production channels mentioned before, which are the largest in the SM. For all of them our predictions improve upon existing ones in at least one aspect. We shall discuss this point in more details in what fol-H H H H Figure 1: Classes of diagrams for Higgs pair production in hadron hadron collisions: double Higgs production without HHH vertices on the left-hand side, and, on the right-hand side, the contribution due to the Higgs self interaction. Final state particles other than the Higgs bosons are understood. lows. Here, we limit ourselves to pointing out that HH production via gluon-gluon fusion is computed at the NLO in a "loop-improved" EFT approach, using the exact oneloop real-emission and improved one-loop virtual matrix elements; that in the case of ttHH and tjHH production exact NLO QCD results are presented in this paper for the first time; and that by matching NLO computations to parton showers we generate samples of events, also for the first time, for each of the production channels, which can be used for fully realistic simulations, including those at detector level. With the exception of the gluon-gluon fusion process which, being loop-induced, needs an ad-hoc treatment, our results are obtained automatically with the publicly-available version of the Mad-Graph5_aMC@NLO framework [16,17].
In the next section we introduce and review the main features of the Higgs-pair production channels. In section 3 we present the calculation and simulation framework, and in section 4 we collect results for some selected observables together with their uncertainties. We summarise our findings and prospects in the conclusions.

Higgs pair production channels
In the SM, the diagrams contributing to Higgs pair production can be organised in two classes (see fig. 1): those where both Higgs bosons couple only to vector bosons or to heavy quarks, and those that feature the Higgs self coupling.
The dominant channel for Higgs pair production is gluon-gluon fusion via virtual top quarks, i.e., box and triangle diagrams. This process therefore starts at the leading order with a loop, exactly as single-Higgs production. In contrast with the latter, however, the effective field theory approach (where Higgs-gluons vertices are included in the lagrangian L HEFT = α S /(3πv 2 )(φ † φ)GG, G being the QCD field tensor) provides only a rough approximation for total rates, and a very poor one for distributions [18,10]. Better predictions, which take loop effects into account exactly, need therefore to be employed in actual phenomenological and experimental studies. Results have been available for some time, and implemented in the code HPAIR [12,13], which deals with both the SM and the MSSM, but is only capable of computing total cross sections. In HPAIR the NLO calculation is essentially performed with EFT techniques; the exact oneloop Born amplitudes are however employed as leadingorder contribution to the NLO cross section, and used to reweight (after the integration over the polar scattering angle) the HEFT virtual-and real-emission matrix elements. In this work we improve on the HPAIR approach on several counts. Firstly, we include the exact one-loop results not only for the 2 → 2 Born amplitudes, but also for the 2 → 3 real-emission processes, which we compute with MadLoop [19]. In other words, the only approximation made at the level of matrix elements is that for the finite part of the two-loop virtual corrections which, being presently unknown, is approximated by the corresponding one-loop HEFT result reweighted (without any intermediate integration) with the exact one-loop Born amplitude. Secondly, in the loops we make use of the complex mass (and Yukawa) scheme for the top quark [20,21]. Thirdly, our results are fully differential, and can be used to obtain any distribution after matching with parton shower. In summary, our predictions improve both on the total cross sections that can be obtained with HPAIR, and on the differential, hadron-level (i.e., showered) observables recently presented in ref. [22,23] (which do not include virtual effects, and are therefore akin to tree-level merged results). We also stress we do not make use of the recentlyderived 1/m t effects at the NLO accuracy [24], of the NNLO HEFT results for total rates [25] and of threshold resummation [26]. More details on the procedure employed in this work will be presented elsewhere [27].
The second-largest production channel is vector boson fusion (VBF). In this case the NLO QCD corrections are trivial, as they involve the same contributions as for single-Higgs production. In VBF we compute only vertex loop corrections, i.e., the finite part of the pentagon and hexagon loop diagrams are discarded for simplicity. These contributions only affect interferences between diagrams that feature identical quarks, which are negligibly small already at the LO. NLO results have been presented in the literature (see e.g. [15]) only for total rates. In this paper we study, for the first time, differential observables for VBF in the SM at fixed NLO and matched to parton showers, showing distributions for the latter. Distributions at fixed NLO in the two Higgs doublet model have appeared in ref. [28]. We point out that, although NNLO corrections to the total VBF cross sections are not known, they could be easily computed following the approach of ref. [29].
At variance with single-Higgs production, the production of a Higgs pair in association with a tt pair is the third most important process and, in fact, it is even larger than VBF at high Higgs-pair transverse momenta, or for collider centre-of-mass energies higher than that of the LHC. The inclusion of NLO QCD corrections in this process has never been achieved prior to this work, even at a fully inclusive level, as it involves thousands of Feynman dia-grams of high complexity, such as pentagon and hexagon loops. Our framework, however, has no problems in handling it in a fully automatic way. For instance, the total (sequential) CPU time required to generate one million unweighted events and to obtain a cross section accurate at the per-mil level is about one hundred and sixty CPU hours on a 2.3GHz machine. This renders the computation feasible on a medium-size (30 core) cluster in a few hours.
The channels of vector boson associated production are technically the easiest ones, as all QCD corrections factorise and are relevant only to the initial state. As in the case of VBF, we improve upon existing NLO results by giving one the possibility of studying fully-differential observables; in this work, we do not include the finite one-loop, gg-initiated contributions to ZHH production, which however can also be handled by MadLoop. Our results correspond to on-shell final state vector bosons. NNLO QCD corrections to total cross sections are known to be small [15].
Finally, in order to provide the complete set of possibly interesting final states, we also compute for the first time at the NLO the single-top associated production, by including both s-and t-channel contributions and by considering both top and anti-top in the final state. The corresponding cross sections are tiny at the LHC, and of very limited phenomenological relevance in the SM. However, this process is at least of academic interest because it is sensitive to couplings to both vector bosons and top quarks, and to their relative phases. In addition to that, given that it has the largest sensitivity to the self-coupling λ, it might become relevant at a future proton-proton 100 TeV collider.

Setup
As was mentioned above, apart from the gluon-gluon fusion channel, all results presented in this work have been obtained in a fully automatic way with MadGraph5_aMC-@NLO [16,17]. This program is designed to perform the computation of tree-level and NLO cross sections, including their matching to parton showers and the merging of samples with different parton multiplicities. A user can generate a given process through a simple shell interface (in a manner fully analogous to that of MadGraph5 [30]), with the corresponding self-contained code being generated on the fly. While it is possible to obtain predictions at the ME+PS level (i.e., with the MLM-k T tree-level merging technique of refs. [31,32,33] and its analogues) in this work we limit ourselves to NLO+PS results. This is because the smallness of the Higgs-pair cross sections rather emphasises observables which are inclusive with respect to extra radiation, and for which NLO-level results have to be preferred to tree-level merged ones, since they provide one with better predictions for absolute normalisations and for theoretical uncertainties.
Within MadGraph5_aMC@NLO, any NLO computation is performed by means of two independent mod-ules: MadFKS [34] takes care of the Born and of the realemission amplitudes, and it also carries out the subtraction of the infrared singularities according to the FKS prescription [35,36] as well as the generation of the Monte Carlo subtraction terms required by the MC@NLO method [37].
In our simulations we set the Higgs mass equal to m H = 125 GeV. Parton distributions functions (PDFs) are evaluated by using the MSTW2008 (LO and NLO) set in the five-flavour scheme [42]. b-quark masses as well as their coupling to the Higgs are neglected. For the sake of brevity, we only show observables related to the Higgs bosons and therefore we have left the latter stable in the simulations. We stress, however, that the top quarks and the vector bosons that appear in the final states can be decayed with the built-in MadSpin package [43], which allows one to include all spin-correlation effects. On the other hand, Higgs decays can be handled correctly also by the Monte Carlos, thanks to the Higgs being a spin-0 particle.
The code allows full flexibility as far as the choice of the renormalisation and factorisation scales µ R,F is concerned.
The central values of these scales have been chosen as follows. For gluon-gluon fusion, VBF, and V HH production we set µ 0 = m HH /2, m W and m V HH , respectively. For 1/4 , m T being the transverse energy of the corresponding particle, as we find that in this way the cross section displays a rather stable behaviour. For single-top associated production tjHH we simply use the fixed value µ 0 = m H + m t /2. Scale and PDF uncertainties can be evaluated at no extra computational cost thanks to the reweighting technique introduced in ref. [44], the user deciding the range of variation. In addition, such information is available on an event-by-event basis and therefore uncertainty bands can be plotted for any differential observable of interest. In our analysis we vary independently the scales in the range 1/2µ 0 < µ R , µ F < 2µ 0 . PDF uncertainties at the 68% C.L. are obtained by following the prescription given by the MSTW collaboration [42].
For the studies shown in this paper we employ HER-WIG6 [45] and Pythia8 [46] for parton shower and hadronisation. The matching to HERWIG++ [47] and Py-thia6 [48] (virtuality ordered, plus p T ordering for processes with no final-state radiation) is also available in MadGraph5_ aMC@NLO.

Results
We start by presenting in fig. 2 the predictions for the total rates at proton-proton colliders with up to 100 TeV c.m. energy. The thickness of the curves corresponds to the  Firstly, contrary to what happens in single-Higgs production, the top-pair associated channel is the thirdlargest starting at about √ s =10 TeV, and becomes the second-largest when c.m. energies approach √ s =100 TeV. Secondly, the theoretical uncertainties due to scale variations in the three most important processes (gluon-gluon fusion, VBF, and tt associated production) are sizably reduced by the inclusion of the NLO corrections. Thirdly, the K-factor is always slightly larger than one, except for gluon-gluon fusion where it is of order two, and for the toppair associated channel where it is smaller than one. Finally, PDF uncertainties are comparable to NLO scale uncertainties, except in the case of gluon-gluon fusion, where the latter are dominant. In the case of V HH and tjHH production it is manifest that the standard procedure of determining uncertainties due to missing higher orders by varying the scales does not give a reliable estimate, as NLO corrections for these processes are much larger than the LO scale dependence band. This is due to two facts: these processes are purely electro-weak processes at the LO, and therefore the scale uncertainties are artificially small; furthermore in the kinematic region probed by these processes, the quark-gluon initiated channel which opens up at the NLO can be important.
In fig. 3 we display total LO and NLO cross sections for the six dominant HH production channels at the LHC with √ s =14 TeV, as a function of the self-interaction coupling λ. The dashed (solid) lines and light-(dark-)colour bands correspond to the LO (NLO) results and to the scale and PDF uncertainties added linearly. The SM value of the cross section corresponds to λ/λ SM = 1. The sensitivity of the total cross sections to the actual value of λ depends in a non-trivial way on the relative couplings of the Higgs to vector bosons and top quarks, and on the kinematics in a way that is a difficult to predict a priori, i.e., without an explicit calculation. The reduction of the scale uncertainties that affect the gg → HH, VBF, and ttHH rates, due to the inclusion of NLO corrections, and pointed out in table 1 for the SM, is seen here also for values of λ = λ SM .
We then plot typical distributions for all channels and at the 14 TeV LHC, which we obtain by generating samples of events at parton level, which are then showered with Pythia8 (solid) and HERWIG6 (dashes). Being tiny at the 14 TeV LHC, we do not show the results for single-top associated production. We present observables at the NLO+PS accuracy in the main frames of the plots: the transverse momentum of the hardest (softest) Higgs in fig. 4 (fig. 5), and the transverse momentum ( fig. 6)    NLO effects appear as overall rescaling factors only in some distributions and on a channel-dependent basis. Moreover, differences between results obtained with the two different shower programs are very mild for all observables and anyway decreasing when going from LO to NLO. In addition, we have checked that differences in the distributions between NLO+Pythia8/NLO+HERWIG6 and NLO fixed-order results are quite small (typically less than a few percent), and that in general the NLO+PS results are slightly closer to each other than to the corresponding NLO fixed-order results. The only exception to this general behaviour is seen in some distributions relevant to the ttHH channel, where the NLO curves lie between the NLO+PS ones, which can differ up to 10% (still within the scale uncertainties). However, in all these cases, the corresponding differences at the LO+PS level are systematically larger than for NLO+PS, which thus confirms the stabilisation trend usually seen when higher-order corrections are included.
NLO corrections in the gluon-gluon fusion channel are important rate-wise, yet the shapes are not strongly affected, as is apparent from the rather flat K-factors. NLO corrections in VBF production are of order 10-20%, and modify the shape of the distributions towards low massscale values. NLO effects in ttHH production lead to a drastic reduction of the scale uncertainties, and to minor changes in the shapes, except for m HH . The associated vector boson channels display very similar features: rather flat K-factors for all the distributions studied, except for p T (HH) where the NLO corrections become more and more important at high p T .

Conclusions
Assessing the nature of the newly discovered boson will need a campaign of measurements to be performed at the LHC at an unprecedented accuracy. One of the key processes in this endeavour is Higgs-pair production. Not only it gives one the possibility of measuring the value of the Higgs self-coupling λ, but also of putting constraints on several, still viable, new-physics scenarios. All such measurements will need accurate SM predictions for total cross sections (in order to extract information on the couplings) and differential distributions (in order to establish acceptances and identify optimal selection cuts), including reliable estimates of the theoretical uncertainties.
In this Letter we have presented the first predictions at the NLO accuracy matched with parton shower for all the relevant Higgs-pair production channels in the Standard Model. We find that, as expected, including NLO corrections leads to a reduction of the theoretical uncertainties, especially significant in the gluon-initiated channels, and provides one with reliable predictions for the kinematic distributions of final state particles. With the exception of the gluon-gluon fusion process, which needs an ad-hoc treatment and for which a dedicated procedure and code have been developed [27], all the results presented here can be easily and automatically reproduced with the publicly available version of the MadGraph5_aMC@NLO code [17].
The extension of our study to models that feature physics beyond the SM is in progress.