Higgs pair production via gluon fusion in the Two-Higgs-Doublet Model

We study the production of Higgs boson pairs via gluon fusion at the LHC in the Two-Higgs-Doublet Model. We present predictions at NLO accuracy in QCD, matched to parton showers through the MC@NLO method. A dedicated reweighting technique is used to improve the NLO calculation upon the infinite top-mass limit. We perform our calculation within the MadGraph5_aMC@NLO framework, along with the 2HDM implementation based on the NLOCT package. The inclusion of the NLO corrections leads to large K-factors and significantly reduced theoretical uncertainties. We examine the seven 2HDM Higgs pair combinations using a number of representative 2HDM scenarios. We show how the model-specific features modify the Higgs pair total rates and distribution shapes, leading to trademark signatures of an extended Higgs sector.


Introduction
Recent LHC data provide evidence that the scalar particle observed at the LHC is the one predicted by the electroweak symmetry breaking mechanism [1,2], as implemented in the Standard Model (SM) [2]. The mechanism predicts the strengths of all Higgs boson couplings to be uniquely determined by the masses of the elementary particles. This also includes the triple and quartic Higgs boson self-couplings, which in the SM are linked to the Higgs mass and the vacuum expectation value (VEV), and are therefore fully fixed after the Higgs mass measurement. Current measurements of the Higgs couplings to fermions and vector bosons [3, 4] agree within 10-20% with the SM predictions, while no direct information is available on the Higgs boson self-interactions. In view of this experimental picture, possible non-minimal Higgs sectors are constrained, but certainly not ruled out. While some new physics models are no longer compatible with the present collider data, there are still many extended Higgs sectors which can accommodate a ∼ 126 GeV Higgs boson with coupling strengths similar to the SM. One prime candidate for such a theory is the Two-Higgs-Doublet Model (2HDM) [5]. This model provides a simple UV-complete perturbative extension of the SM, which we can view as the low-energy Higgs sector of more fundamental theories such as the Minimal Supersymmetric Standard Model (MSSM) [6,7], GUTs [8,9], composite Higgs models [10][11][12], and little Higgs models [13,14]. Aside from its very rich and distinctive phenomenology, the 2HDM also sets the ground for novel approaches to diverse unsettled conundrums, from e.g. the origin of neutrino masses [15] to Naturalness [16], Electroweak Baryogenesis [17][18][19][20] and Dark Matter [21][22][23]. If effectively realized in Nature, the 2HDM could lead to genuine indications of new physics at colliders. For instance, the direct production of additional heavy scalars offers excellent prospects to identify these novel particles through a variety of decay modes and final-state signatures [24][25][26][27][28][29][30][31][32][33]. Not less important are the indirect signatures from such extended Higgs sectors, which could arise via modified Higgs couplings to the fermions and gauge bosons, as well as through modified Higgs self-interactions.
At present, the experimental reconstruction of the Higgs potential is a major step towards the fundamental understanding of the EW symmetry breaking mechanism. This task can only be fully accomplished through the direct measurement of the Higgs boson three-point and four-point interactions [34][35][36]. Multiple Higgs boson production processes are therefore instrumental in this endeavour, not only because they directly depend on the Higgs self-couplings, but also because they are sensitive to possible new heavier states and/or to higher-dimensional operators [37][38][39][40][41][42][43][44]. Extracting the Higgs self-couplings from collider data is known to be an arduous task [45,46]. While the triple Higgs rates lie beyond the reach of the LHC capabilities [47,48], the prospects of measuring Higgs pair production are better, but still challenging. There are multiple production channels which lead to Higgs boson pairs at hadron colliders, including vector boson fusion; associated production with gauge bosons and heavy quarks; and gluon fusion. The latter is dominant in the SM and gives an approximate next-to-leading order (NLO) cross section of approximately 35 fb at the 14 TeV LHC [49]. Recent studies [42,43,[50][51][52][53][54][55] have investigated the potential of digging out the di-Higgs signal over its backgrounds in various Higgs decay channels, among them γγbb, W + W − bb, bbτ + τ − and bbbb, assisted e.g. by jet substructure techniques. The current picture is that, aside from potential new physics effects, the extraction of the coupling g hhh will require copious integrated luminosity. Optimistic estimates point towards values of 3000 fb −1 at 14 TeV in order to reach an accuracy of ∼ 40% [56]. The reason is not only the rather modest total rate, but also the limited sensitivity to the trilinear coupling g hhh and the large theoretical uncertainties. With this research program ahead, accurate predictions for the total and differential di-Higgs rates are in order -for the SM and beyond, while more sophisticated experimental analyses are required to explore the precise reach of the LHC.
In this work we investigate the phenomenological possibilities of Higgs pair production via gluon fusion beyond the SM. We resort to the 2HDM as an illustrative new physics framework, and present results for all seven combinations of 2HDM Higgs pair final states. We compute the total and differential Higgs pair rates at NLO accuracy in QCD, matched to the Pythia8 [57] parton shower (PS) with the MC@NLO method. A dedicated reweighting strategy is employed to improve the NLO calculation beyond the infinite top-mass limit, by also employing the exact real-emission matrix elements [58]. On the practical side, we carry out our computation using the MadGraph5 aMC@NLO framework [59]. The results presented here improve upon, and further extend, the earlier studies in the literature. To the best of our knowledge, our work provides for the first time NLO+PS event samples for these processes, which can be readily used for realistic simulations, in-Type I Type II 1 + ∆ h 0 t cos α sin β = 1 + ξ/ tan β − ξ 2 /2 + O(ξ 3 ) cos α sin β = 1 + ξ/ tan β − ξ 2 /2 + O(ξ 3 ) sin α sin β = −1/ tan β + ξ + ξ 2 /(2 tan β) + O(ξ 3 ) cos α cos β = tan β + ξ − ξ 2 /2 tan β + O(ξ 3 ) cluding those at the detector level once the Higgs bosons are allowed to decay. Accurate predictions for the differential rates are also important, as they can be used to identify the distinctive properties of the signal kinematics and compare them to the backgrounds. This is e.g. instrumental in the context of jet substructure techniques [41,42], which are known to improve upon the more traditional search strategies based on inclusive observables.
The remainder of this paper is organized as follows. In Section 2.1 we begin by outlining the 2HDM setup, field content, coupling structure, and parameter space constraints. Section 2.2 defines a series of 2HDM benchmark scenarios, which are constructed to cover all representative phenomenological features of the model. In Section 2.3 we move on to describe the theoretical structure of Higgs pair production at leading and next-to-leading order. Before closing this first part, the technical setup of the calculation is introduced in Section 3. In the second part of the paper we present a comprehensive numerical analysis of all Higgs pair production channels within the 2HDM. Total rates are documented and discussed in Section 4.1, while in Section 4.2 we focus on the light di-Higgs differential distributions. Finally, we summarize and conclude in Section 5.

The Two-Higgs-Doublet Model
The Two-Higgs-Doublet Model (2HDM) [5] extends the minimal scalar sector of the SM by introducing a second SU (2) L doublet Φ 2 with weak hypercharge Y = +1. This gives rise to an enlarged particle content with five physical Higgs bosons: one light (resp. heavy) neutral, CP-even state h 0 (resp. H 0 ); one neutral, CP-odd state A 0 ; and two charged Higgs bosons H ± . Throughout the paper we identify the state h 0 with the Higgs particle observed at the LHC and fix its mass to m h 0 = 126 GeV. The mixing angle α is introduced to diagonalize the CP-even squared mass matrix. Assuming natural flavor conservation [60], the absence of tree-level flavor-changing neutral-current (FCNC) interactions is protected by a global, flavor-blind, Z 2 discrete symmetry Φ i → (−1) i Φ i . The latter is approximate up to the soft-breaking mass term L soft ⊃ m 2 12 Φ † 1 Φ 2 +h.c. We neglect extra sources of CPviolation by considering real mass terms and self-couplings in the Higgs potential. After electroweak symmetry breaking, the neutral components of the Higgs doublets acquire real The ratio of the two VEVs is given by tan β ≡ v 2 /v 1 . Overall, we are left with 7 free input parameters, which we can sort out as: The convention 0 ≤ β − α < π (with 0 < β < π/2) guarantees that the Higgs coupling to the weak gauge bosons g 2HDM hV V has the same sign in the 2HDM and in the SM. This criterion fixes the possible sign ambiguities in the generic parametrization of the model [61,62].
The different possible choices of fermion field transformations under Z 2 lead to different Yukawa coupling patterns. We will hereafter concentrate on the two canonical setups: i) type-I, in which all fermions couple to just one of the Higgs doublets; and ii) type-II, where up-type (down-type) fermions couple exclusively to Φ 2 (Φ 1 ). The resulting Yukawa couplings deviate from the SM in a way that we parametrize in the notation of [63], Analytic expressions for these coupling shifts are provided in Table 1 for both type-I and type-II 2HDMs. For further reference, also the trilinear Higgs self-couplings involving the CP-even neutral states are quoted below in Table 2.

Benchmarks
Phenomenologically viable 2HDM scenarios satisfying all model constraints and compatible with the LHC data have been extensively scrutinized in the literature [70,. These studies highlight a preference for a low-mass, SM-like Higgs field h 0 along with heavier Higgs companions. In the decoupling limit [69], the 2HDM can be mapped onto an effective theory, whose expansion parameter ξ ≡ cos(β − α) ∼ v 2 /M 2 heavy determines the hierarchy between the light m h 0 = O(v) and the heavy scalar masses [138,139]. The decoupling condition ξ = cos(β − α) ≪ 1 correlates the two mixing angles through cos β ∼ sin α, which in the limit of large tan β can be expressed by 3) The behavior of the relevant Higgs interactions in the decoupling limit is explicitly shown in Tables 1-2. Notice that even for ξ ≪ 1 some of the couplings may be substantially shifted. This behavior appears for instance in the tan β ≫ 1 (tan β ≪ 1) regimes as a reflect of the so-called delayed decoupling [140]. As for the Yukawa couplings, these shifts may be more (in type-I) or less (in type-II) correlated within each fermion generation and can lead to enhanced, suppressed, or even sign-flipped couplings [141]. In turn, the triple     Table 3. All couplings are normalized to their SM counterparts.
Higgs self-coupling g h 0 h 0 h 0 may be enhanced by up to 100% above the SM in type-I models -while for type-II, the LHC data favour g h 0 h 0 h 0 g SM HHH with allowed suppressions up to O(50)% [142]. Let us also note the particular decoupling limits g h 0 h 0 h 0 → g SM HHH and g H 0 h 0 h 0 → 0 for ξ ≪ 1. The potentially large Higgs self-coupling deviations constitute a genuine trait of the 2HDM, with no counterpart in e.g. the Higgs sector of the MSSM. In the latter case, Supersymmetry relates all Higgs self-couplings to the gauge couplings, implying that their size becomes restricted. The fact that the conventional decoupling limit ξ ≪ 1 is not the unique 2HDM setup consistent with a SM-like ∼ 126 GeV resonance is another remarkable feature of the model. In the so-called alignment limit [69,130,143,144], a SM-like Higgs state can still be made compatible with additional Higgs bosons as light as ∼ 200 GeV [70,. Interestingly, this low-mass region m H 0 250 GeV is also elusive to direct searches -mainly because of the problematic background isolation [145].
These rich phenomenological possibilities are captured by the set of 2HDM benchmark scenarios which we introduce in Table 3. We employ them further down in Section 4 to examine the distinctive 2HDM signatures on the Higgs pair production observables. They have been constructed in agreement with all up-to-date parameter space constraints, which we have included through an in-house interface of the public tools 2HDMC [146], HiggsBounds [147,148], SuperIso [149,150] and HiggsSignals [151,152] along with additional routines of our own. For the most recent direct heavy Higgs searches [100,102,103] not yet available from HiggsBounds, it has been checked explicitly that the benchmarks evade the exclusion bounds. In order to better illustrate certain model features, in some scenarios we tolerate deviations slightly beyond 1σ in the averaged Higgs signal strength.
In Table 4 we quote the numerical values for sample Yukawa and Higgs self-couplings (a selection which is relevant to the light Higgs pair production) for all seven 2HDM benchmarks defined in Table 3. All couplings are normalized to their SM counterparts, as denoted byĝ hxx ≡ g 2HDM hxx /g SM Hxx , where H stands for the SM Higgs boson. The heavy Higgs trilinear coupling is normalized asĝ The key physics properties of the different 2HDM scenarios can be summarized as follows: • B1: Moderate mass hierarchy -taken from benchmark H1 in Ref. [153]. It corresponds to a type II 2HDM with moderate (viz. 300 − 500 GeV) heavy Higgs masses. The small values of tan β and ξ ≡ cos(β − α) guarantee that all light Higgs boson couplings remain very close to the SM -with an O(5)% suppression in the triple Higgs self coupling and an O(10)% enhancement in the bottom Yukawa. This is also the reason why this benchmark evades the recent CMS [100] and ATLAS limits [102]. As we will show in Section 4, the major new physics effects in this case originate from the resonant heavy Higgs-mediated contribution gg → H 0 → h 0 h 0 . These effects are particularly enhanced here due to the dominant heavy Higgs cascade decay H 0 → h 0 h 0 , whose branching fraction is BR h 0 h 0 ≃ 0.6.
• B2: Large mass hierarchy -taken from benchmark a.1 again in Ref. [153]. At variance with the B1 scenario above, all of the additional Higgs bosons have in this case masses of O(700) GeV, so that they effectively decouple. The very large Z 2 soft-breaking term m 2 12 is responsible for the O(40)% suppression of the trilinear Higgs self-coupling g h 0 h 0 h 0 , while at the same time it enlarges g H 0 h 0 h 0 .
• B3: Non-resonant SM limit -in which the relatively light CP-even Higgs com- The parameter choice is inspired by the benchmarks of class [c] from Ref. [153]. They are characterized by a rather light Higgs spectrum, along with moderate values for tan β and m 2 12 . The Yukawa interactions are once again fixed according to a type-II setup. The resulting coupling patterns in this case are SM-like for the lightest Higgs field, with mild deviations not larger than 1%. For the heavy neutral field H 0 we find a suppressed (resp. enhanced) Yukawa coupling to the top (resp. bottom) with opposite signs, and a strongly reduced trilinear coupling g H 0 h 0 h 0 .
• B4: Enhanced triple Higgs self-coupling -we assume i) type-I Yukawa couplings, for which the LHC Higgs signal strength imposes weaker constraints; ii) a small tan β value; iii) a value of ξ for which the Higgs coupling to the gauge bosons g hV V ∼ sin(β − α) ∼ 1 − ξ 2 /2 is reduced by ∼ O(1)%; iv) a soft-breaking mass term close to the unitarity limit. Such parameter setup leads to an O(30)% enhancement for the light Higgs triple self-coupling, and of O(10)% for the heavy quark Yukawas. The negative sign m 2 12 < 0 protects the stability of the vacuum. The relatively low mass m H 0 < 2m h 0 prevents the resonant heavy Higgs production and allows to better access the genuine model-specific imprints on the di-Higgs production observables.
• B5: Fermiophobic heavy Higgs -a situation which can only be accomplished in type-I models for sin α = 0. On the one hand, compatibility with the LHC Higgs signal strength enforces ξ ≪ 1, meaning that these scenarios are only viable at large tan β, according to Eq. (2.3). The value of m 2 12 is tailored to fulfill the vacuum stability and unitarity bounds, which become very tight at large tan β. Notice that for ξ ≪ 1 the (relatively low-mass) heavy Higgs H 0 hardly couples to the gauge bosons and the light Higgs h 0 , while by construction it cannot couple to any fermion either. Consequently, in this case there is no heavy Higgs contribution to the light di-Higgs production gg → H 0 → h 0 h 0 .
• B6 and B7: Enhanced and sign-flipped bottom Yukawa -both instances are possible in type-II models at large tan β, where the bottom Yukawa coupling g h 0 bb = (1 − ξ tan β) g SM Hbb may yield either g h 0 bb ≃ −g SM Hbb or g h 0 bb > g SM Hbb , in correspondence to the two branches sin α > 0 and sin α < 0 along the decoupling condition of Eq. (2.3). The possibility of a strongly modified bottom Yukawa, with all of the remaining couplings being SM-like, is a trademark property of the 2HDM [141], and relies on the delayed decoupling behavior mentioned earlier [140]. A similar mechanism for tan β < 1 should in principle lead to a sign-reverted top Yukawa, although this situation is in practice disfavored by the flavor constraints.

Higgs pair production in the 2HDM
The pair production of Higgs bosons at hadron colliders can proceed via weak gauge boson fusion [154][155][156][157][158][159], double Higgs-strahlung off the W and Z bosons [160], and gluon fusion [158,[161][162][163]. Because of the large gluon luminosity in the high-energy proton beams, the gluon gluon fusion channel dominates. Predictions for the SM are known at NLO [164] and NNLO [165] in QCD, both in the infinite top-mass effective theory. Further studies have reported on the subleading O(1/m 2 t ) terms [166], threshold resummation [167], as well as on gluon fusion results merged to one jet [168]. More recently, predictions for all Higgs pair production channels at NLO and matched to parton showers have been presented in Ref. [49]. These studies conclude that higher-order effects are large, especially for the dominant gluon fusion channel, and that including them significantly reduces the theoretical uncertainties.
In the context of new physics, double Higgs production has been addressed in a variety of models, for instance the MSSM [34,162], the NMSSM [169][170][171], the 2HDM [37,172,173], Figure 1: Generic Feynman diagrams describing the production of neutral Higgs boson pairs (h = h 0 , H 0 , A 0 ) in the 2HDM through gluon fusion at leading order. The Feynman diagrams have been generated using FeynArts.sty [185]. extended colored sectors [38,40,174], Little Higgs [175,176], Higgs portal [41,177], a flavor symmetry model [178] and Composite Higgs models [179,180]. Dedicated studies on the charged Higgs pair case have been presented e.g. in Refs. [181,182]. Model-independent approaches based on anomalous couplings and effective operators have been considered as well [39,183,184].
In this study we focus on the Higgs pair production in the 2HDM, and consider the seven possible final-state double Higgs combinations. These may be sorted out into the following three categories:

Charged Higgs boson pairs: H
There are two leading-order (LO) mechanisms contributing to the gluon fusion di-Higgs channels pp(gg) → hh (with h = h 0 , H 0 , A 0 ), whose generic Feynman diagrams we display in Figure 1. These correspond to: 1. Triangle topologies, which give rise to O(G F α sĝq ) contributions through the schannel exchange of (at least) one neutral Higgs boson. The Higgs boson couples to the gluons via the usual heavy-quark loops.
2. Box topologies, which contribute at O(G F α sĝ 2 q ) through the virtual heavy quark exchange.
The normalized heavy quark Yukawasĝ q ≡ g q /g SM q are related to the 2HDM coupling shifts through Eq. (2.2) and Table 1. While the triangle topologies have a linear dependence inĝ q and in the Higgs self-coupling, the boxes are quadratic in the former. In the SM, the two topologies interfere destructively. This effect is particularly strong near the threshold [186] and explains in part why the total rates are quite modest. One additional destructive interference arises between the top and bottom-mediated loops, although the bottom quark effects are very small in the SM.
For category 2), we have additional tree-level contributions for which a quark-antiquark pair annihilates into a virtual Z-boson (cf. the right-most diagram in Figure 2). For the charged Higgs boson pairs of category 3), also the photon exchange from the qq-annihilation contributes. The relative size of these qq-initiated subchannels is quantified in Section 4.1.
To gain further insight into the structure of the (loop-induced) gluon fusion mechanism, let us focus on the neutral CP-even Higgs pairs of category 1). The partonic cross-section at leading order may be written as The functions F △ , F and G denote the one-loop gauge-invariant form factors from the triangle and the box contributions, and depend on the Mandelstam kinematical invariants and the quark and Higgs masses. The variablet corresponds to the momentum transfer squared from one of the incoming initial-state gluons to one of the Higgs bosons in the final state. The sum over the neutral CP-even Higgs fields h = h 0 , H 0 accounts for the respective Higgs-mediated triangle topologies. Notice that the s-channel exchange of the CP-odd Higgs A 0 is forbidden by CP-conservation. The model-dependent Yukawa couplings are related to their SM counterparts through the coupling shifts quoted in Table 1. As for the Higgs self-interactions (cf. Table 2), we normalize them as λ hhh ≡ i v g hhh , with the SM Higgs self-coupling being g SM HHH = −3im 2 H /v. Finally, the symmetry factor c ij = 1/2 (1) for i = j (i = j) properly accounts for the cases with identical (different) final-state particles. The overall normalization of Eq. (2.4) is consistent with the notation of Ref. [162]. The corresponding hadronic cross section is obtained by convoluting Eq. (2.4) with the gluon luminosities.
The two box form factors F , G correspond to the two possible S-wave and D-wave contributions. These are linked to equal (resp. opposite) gluon helicities, which respectively add to a total angular momentum J z = 0 (resp. J z = 2) along the collision axis. Instead, the triangle diagrams only contribute for J z = 0 and hence give rise to a single form factor F △ . The large (resp. low) mass limits of the loop form factors read [162]: The opposite signs for these form factors reflect the negative interference patterns between i) boxes and triangles; ii) top and bottom-mediated loops.
The Higgs pair total rates and distribution shapes are determined by the size and relative signs of these different loop-induced contributions, and the way they interplay according to Eq. (2.4). The results are thereby sensitive to the underlying 2HDM dynamics. For instance, sign-flipped couplings may revert the partial cancellation between the triangle and the box loops. The possible variations in the Yukawa and the Higgs self-couplings may either enhance or suppress these interference patterns. If kinematically allowed, the onshell production of a heavy neutral CP-even Higgs through the heavy top triangle, followed by the cascade decay H 0 → h 0 h 0 , will overwhelm the SM expectations. All these cases will be examined in Section 4 with the help of the different 2HDM benchmark scenarios defined in Table 3.

Next-to-leading order corrections
The NLO QCD corrections arise at O(α 3 s ) and are obviously linked to the color charges of the initial partons. They originate from i) virtual gluon exchange; and ii) light parton radiation. The NLO virtual corrections to gg → hh are intrinsically a two-loop effect, not accessible by current loop calculation techniques. As an alternative, Higgs pair production studies at NLO and NNLO [164,165] resort to an Effective Field Theory framework. In the so-called Higgs Effective Field Theory (HEFT) approach, the effective three-point and four-point Higgs boson interactions with the gluons are obtained by integrating out the topquark. These effective interactions, which are connected to the low-energy theorems [187][188][189], can be written in terms of higher-dimensional operators and cast into the Lagrangian: where G A µν denotes the gluon field-strength tensor while h = h 0 , H 0 . Similar expressions may be derived for the gluon couplings to the pseudoscalar Higgs, cf. e.g. [164]. The (1 + ∆ h t ) factor (cf. Table 1) accounts for the model-dependent shifts to the effective ggh and gghh interactions. The relative sign between the two couplings constitutes the low-energy equivalent of the destructive interference between the triangle and the box topologies.
A selection of the Feynman diagrams describing the relevant NLO QCD effects are illustrated in Figure 3. The NLO virtual corrections in the HEFT picture are obtained from one-loop topologies, as shown by the first two Feynman diagrams of Figure 3. The shaded blobs denote the effective three-point and four-point Higgs couplings to the gluons. In our calculation, these HEFT virtual corrections are combined with the exact 2 → 3 real-emission amplitudes (cf. the right-most diagrams in Figure 3) as an alternative to the more traditional approach, which only uses the exact LO result to improve upon the infinite top-mass limit. Further details are provided in Section 3. Aside from the leading-order gg-initiated partonic subchannel, notice that for the NLO real-emission contributions also the mixed qg fusion channels open up.

Calculation setup
In this section we describe in detail the technical setup of our calculation. As already mentioned, one important aspect is the treatment of the NLO corrections to the gluon fusion channels. Given that this production mechanism is loop-induced at LO, two ingredients would be needed for an exact NLO calculation: i) the one-loop 2 → 3 real-emission amplitudes; ii) the two-loop 2 → 2 virtual correction amplitudes. While the former can be calculated by means of standard techniques, the latter are still beyond the reach of the state-of-the-art calculations. On the other hand, the HEFT description from Eq. (2.6) relies on the infinite top-mass approximation and thereby has a limited validity. It is in fact well-known that it provides only a rough estimate of total rates [38], while it poorly reproduces the kinematical distributions [42,190]. The usual approach in higher-order studies is to extract the QCD corrections from the HEFT, and then employ the exact one-loop LO amplitudes to reweigh the HEFT virtual-and real-emission matrix elements. This conventional strategy, as implemented in HPAIR [162,164], will be hereafter referred to as the "born-improved" approach. One important point in this procedure is that the reweighting of the 2 → 3 real-emission part is based on its factorization into the 2 → 2 LO matrix element times a global factor.
As an alternative, in this work we include in addition the "loop-improved" approach which was first presented in [49]. For completeness, here we describe the method briefly, while a detailed discussion is provided in [58]. In the "loop-improved" calculation we include the exact one-loop results not only for the 2 → 2 LO amplitudes, but also for the NLO 2 → 3 real-emission matrix elements. Therefore, the only approximation we make at the amplitude level concerns the finite part of the NLO virtual corrections which, in the absence of the exact two-loop calculation, is first taken from the one-loop HEFT results, and then reweighted with the exact one-loop LO matrix elements. Including the exact one-loop 2 → 3 matrix elements provides a more accurate description of the tails of the distributions, e.g. at large Higgs boson pair transverse momentum p T (hh). In these phase space regions, in which hard parton emissions take place, the factorization of the 2 → 3 real-emission amplitudes into the 2 → 2 LO amplitudes, as implicit in the "bornimproved" approach, cannot accurately describe the hard parton kinematics. In Section 4, we compare the Higgs pair total rate predictions obtained within the "born-improved" and the "loop-improved" methods, while for the differential rates we concentrate on the "loop-improved" results. Despite the upgraded treatment of hard real emission, our "loopimproved" method still relies on the one-loop HEFT virtual corrections in place of the exact two-loop results. The impact of these exact two-loop results cannot be quantified until they become available. The validity of the approximation is nonetheless limited, not only because the top quark mass is not parametrically large as compared to the other relevant scales of the process, but also because of the enhanced bottom-quark contributions which are possible in the 2HDM. Given that the low-energy theorems do not hold for bottommediated loops, we expect the HEFT virtual results for certain 2HDM configurations to be less reliable than in the SM case. Let us also note that, for the gg-induced production of the mixed CP-even/CP-odd pairs, the NLO virtual corrections must again be approximated by the one-loop HEFT results. As for the s-channel Z-mediated contributions in the latter cases, we use the same factorized form as for the A 0 -mediated ones, which is exact in the m t → ∞, m b → 0 limit as discussed in Ref. [164]. Concerning the NLO realemission corrections, in all these channels we include the full set of non-zero Feynman diagrams contributing to gg → gh 0 A 0 /H 0 A 0 , among them the box diagrams of the form gg → gZ * → gh 0 A 0 /H 0 A 0 . The latter were not included in the analysis of Ref. [164], as they cannot be factorized into the LO amplitude and a universal correction factor. These contributions are infrared and collinear-finite and are expected to be small, but we nevertheless calculate them for completeness. Similar approximations are made for the charged Higgs pair case.
From the technical point of view, our results are obtained within the MadGraph5 aMC@NLO framework, which allows the calculation of LO and NLO cross sections at a fully differential level. In the MadGraph5 aMC@NLO setup, NLO computations are carried out using two independent modules: i) MadFKS [191], which takes care of the Born and the real-emission amplitudes, subtracts the infrared singularities according to the FKS prescription [192,193], and generates the parton shower subtraction terms required by the MC@NLO method [194]; ii) MadLoop [195], which handles the calculation of the one-loop matrix elements using the OPP integrand-reduction method [196] (as implemented in CutTools [197]) and the OpenLoops method [198]. For the 2HDM Higgs pair processes which are tree-level at LO, such as the qq → h 0 A 0 subchannels, the NLO results may be obtained automatically within MadGraph5 aMC@NLO. At variance, for the loop-induced gluon fusion channel events are generated at NLO using the HEFT results and then reweighted using the exact 2 → 2 and 2 → 3 one-loop amplitudes, by means of a separate reweighting routine. All one-loop matrix elements are computed by MadLoop [195]. The final results are in all cases fully differential, so that they can be used to obtain any distribution at will, after matching with the parton shower.
As for the 2HDM implementation, we use the 2HDM@NLO model obtained with the package NLOCT [199]. The model is based on the FeynRules [200] and UFO [201] frameworks. It includes all relevant UV counterterms and R2 vertices needed for the MadLoop calculation, and allows to compute tree-level and one-loop processes within a completely general 2HDM setup. The input parameters for the different benchmarks are submitted to MadGraph5 aMC@NLO [59] through parameter cards in the standard MadGraph setup. The latter are constructed with the help of an in-house modification of the public calculator 2HDMC [146]. Parton distribution functions (PDFs) are evaluated using the MSTW2008 sets with four active flavors at LO and NLO consistently [202]. Therefore, the bb-initiated channels are not included in our calculation. The MadGraph5 aMC@NLO framework allows fully flexible renormalization and factorization scale choices. For all results that follow in this study, we set common factorisation and renormalisation scales, fixing their central value to half the invariant mass of the Higgs pair µ 0 R = µ 0 F = m h i h j /2, where ij account for all possible final-state combinations. This scale choice has been proved to yield perturbatively stable results for the Higgs boson pair production at NLO [49], and to behave very similarly with respect to alternative scale settings used in the literature (cf. e.g. Refs. [164,203]). The strong coupling constant is evaluated at the renormalization scale α s (µ 0 R ) along with the PDFs. The scale and PDF uncertainties are generated automatically and at no extra computational cost within MadGraph5 aMC@NLO, following the reweighting prescription of [204]. In our analysis we vary independently the scales in the range µ 0 /2 < µ 0 R , µ 0 F < 2µ 0 . PDF uncertainties at the 68% C.L. are extracted by following the prescription given by the MSTW collaboration [202]. Pythia8 [57] is used for parton shower and hadronisation. The matching to the Pythia8 parton shower (virtually ordered, and p T -ordered for processes with no final-state radiation, such as the production of Higgs pairs) is also automated within MadGraph5 aMC@NLO . For conciseness, only observables linked to the final-state Higgs pair are shown in the following. This implies that the Higgs bosons are kept stable at the parton shower stage.
The dedicated codes for the calculation of the gluon fusion channels can be downloaded from [205]. This website contains in addition a selection of standalone codes for different Higgs pair production processes in the SM and beyond.

Total rates
We report on the 2HDM predictions for the Higgs pair total cross sections at the LHC. Our results are documented in Tables 5-11. Characteristic model-dependent features are highlighted by the different benchmarks defined in Section 2.2. For the sake of brevity, we select a representative subset of them (B1 -B3) to compute the total rates for all seven combinations of 2HDM Higgs pair final states. The corresponding results are shown in Tables 5-10. For the remaining parameter choices (B4 -B7), we concentrate on the light di-Higgs production h 0 h 0 (cf. Table 11). In addition to the total rate predictions for the (leading) gluon fusion mechanism, Tables 6, 8 and 10 include the contribution from the additional qq-initiated channels. The latter only feature for the mixed neutral Higgs pairs (h 0 A 0 , H 0 A 0 ) and the charged Higgs pairs H + H − , and occur already at tree-level through the s-channel Z 0 boson exchange -as well as through photon exchange in the charged Higgs case. We show in all cases the total rate predictions at LO and NLO accuracy in QCD. For the gluon fusion channels, we explicitly distinguish between the two reweighting schemes in use, as described earlier in Section 3: i) the NLO-"loop-improved", with both the exact 2 → 2 LO matrix elements and the exact one-loop 2 → 3 real-emission matrix elements used to reweight the HEFT results; ii) the NLO-"born-improved", in which only the exact one-loop 2 → 2 LO matrix elements are used for reweighting. The total rate central values are shown along with the scale and PDF uncertainties, whose upper (lower) bands are readily evaluated in MadGraph5 aMC@NLO. The size of the NLO corrections is quantified in terms of the K-factor, provided in the right-most column of the tables. The latter is defined consistently through K ≡ σ NLO /σ LO , where the NLO cross section for the gg-initiated channels stands for the "loop-improved" result. For the qq-initiated channels, the NLO prediction corresponds to the exact result.
Before we discuss the characteristic model-dependent features, some general comments are in order. First of all, we find large NLO corrections, with typical K-factors in the range of 1.5 − 1.7 for the gg-initiated channels. These sizable QCD effects are primarily due to the NLO real emission. This is a well-known fact for processes with color-singlet final states, where there is plenty of phase space to accommodate the radiation of initial-state light partons [42,[206][207][208]. Notice that these K-factors do not change significantly as we compare the different gg-initiated di-Higgs final-states. Likewise, they do not depend on the chosen benchmark and they are quantitatively similar to the SM prediction. We note here for completeness that the corresponding predictions for the SM give a LO cross section of 23.0 fb, a NLO "loop-improved" result of 34.9 fb and a NLO "Born-improved" one of 38.9 fb [49] and therefore a K-factor of 1.52. These traits can be again explained by the structure of the QCD corrections (cf. Section 2.4), which is common to all Higgs pair channels and unlinked to the underlying 2HDM dynamics. Actually, all genuine 2HDM imprints (viz. the modified couplings and the heavy resonances) modify the LO and NLO predictions in exactly the same way, leaving the K-factor values unaltered.   Table 6: Total Higgs pair cross sections via quark-antiquark annihilation pp(qq) → h i h j (in fb) for the mixed CP-even/CP-odd and charged Higgs pair combinations within the 2HDM, in the same setup as Table 5. The 2HDM parameters are fixed to benchmark B1 in Table 3. when we compare the total rates at LO and NLO accuracy. For gluon fusion, we find typical scale uncertainties spanning ∆σ/σ ≡ [σ(µ 0 /2) − σ(2µ 0 )]/σ ∼ 30 − 40% at LO, while they shrink down to ∆σ/σ ∼ 15 − 20% at NLO. The PDF uncertainties lie at the per-cent level and increase with heavier Higgs masses. This behavior can be attributed to the uncertainty in the gluon parton density, which increases with the Bjorken variable x in the kinematically relevant regions for these processes. On the other hand, the qq-initiated subprocesses exhibit an overall milder scale uncertainty, in line with the fact that they do not depend on α s at LO. The corresponding K-factors are also smaller in these cases (cf. Tables 6, 8 and 10) and remain in the range K ∼ 1.2 − 1.3, as expected for Drell-Yan-like processes. Notice also that, unlike the scale uncertainties, the PDF uncertainties slightly grow from the LO to the NLO predictions. This is explained by the additional   Table 8: Total Higgs pair cross sections via quark-antiquark annihilation pp(qq) → h i h j (in fb) for the mixed CP-even/CP-odd and charged Higgs pair combinations within the 2HDM, in the same setup as Table 5. The 2HDM parameters are fixed to benchmark B2 in Table 3. initial-state partons which become active at NLO. The differences between the "loopimproved" and "born-improved" results are typically smaller than 10% and therefore lie within the theoretical uncertainties. Another aspect to mention is the relative size of the (tree-level) qq-initiated channels versus the (loop-induced) gg-fusion mechanism. This varies significantly with the different Higgs pair combinations. For instance, while gluon fusion prevails for h 0 A 0 , in the case of H 0 A 0 we find that the bulk contribution is qq induced. This difference can be traced back to the coupling g h 0 A 0 Z 0 (resp. g H 0 A 0 Z 0 ), which is proportional to cos(β − α) (resp. sin(β − α)) and therefore vanishes (resp. maximizes) in the decoupling limit ξ = cos(β − α) → 0.

1.22
While the K-factors barely depend on the specific 2HDM scenario, the total rates critically rely on the chosen benchmark. In the light Higgs pair case, these can vary from σ NLO ∼ 30 fb up to σ NLO ∼ 2 pb. The latter pb-level rates are roughly two orders of magnitude above the SM expectations, and are linked to the resonant heavy Higgs contribution gg → H 0 → h 0 h 0 . This trait is common to all scenarios in which the cascade decay H 0 → h 0 h 0 is kinematically accessible (m H 0 > 2m h 0 ) and reflects the fact that the resonant production effectively involves one power of G F less compared to the continuum pair production. Such enhancements are tamed if the position of the resonance is shifted away from the di-Higgs threshold 2m h 0 and vanish consistently in the SM limit of ξ → 0 and m H 0 ≫ m h 0 . The fact that the resonant contribution to σ(h 0 h 0 ) decreases with increasing m H 0 values is explained by i) the larger phase space required to produce the intermediate heavy state H 0 ; and ii) the correspondingly lower gluon luminosity, which is probed at larger x-values, the more massive the heavy field becomes. It is also remarkable that these resonant situations are almost insensitive to the actual value of the triple Higgs self-couplings. This is due to the fact that i) the dependence of the resonant production gg → H 0 on the coupling g H 0 h 0 h 0 is cancelled by the dominant decay mode H 0 → h 0 h 0 and ii) the light Higgs-mediated diagrams gg → h 0 * → h 0 h 0 are subdominant with respect to the resonant part, in such a way that the dependence on the coupling g h 0 h 0 h 0 is overshadowed. A number of model-specific fingerprints can be unraveled by comparing the total rate predictions for the different 2HDM benchmarks. These can be ultimately traced back to the characteristic coupling patterns and the extra Higgs boson masses in each scenario, given in Tables 3 and 4. For the resonant scenarios, B1 and B2 , the on-shell production of the heavy Higgs field H 0 is responsible for the total rate enhancement in σ(h 0 h 0 ), by a factor 70 (resp. 4) above the SM. The difference of one order of magnitude between the two scenarios is linked to the fact that, while in B1 the heavy Higgs mass is relatively low and lies right above the di-Higgs threshold m H 0 ≃ 2m h 0 , in scenario B2 the H 0 field is substantially heavier. The total h 0 h 0 rates in the first case thus benefit from the larger on-shell single H 0 rates; as well as from the overly dominant decay mode H 0 → h 0 h 0 . For scenarios B6 and B7 in Table 11 we find milder (∼ 30%) resonant enhancements of the h 0 h 0 cross section as compared to B1, even though the H 0 mass is also quite low in these cases. This is because the H 0 contribution through gg → H 0 → h 0 h 0 is suppressed by the reduced heavy Higgs top Yukawa g H 0 t . Finally, for the non-resonant scenarios (viz. B3, B4 and B5) our predictions for σ(h 0 h 0 ) fall down to values close, or even slightly below the SM. This is explained naturally by the SM-like coupling pattern of B3 and B5 (cf. Table 4). For B4, the reduced rate is a consequence of the stronger destructive interference between i) the triangle-mediated contributions gg → h 0 * → h 0 h 0 and gg → H 0 * → h 0 h 0 , which are enhanced by the larger trilinear couplings; and ii) the box-mediated diagrams. Barring possible on-shell resonances, we conclude that significant deviations from the SM Higgs pair predictions are not possible rate-wise. This is in agreement with the findings reported in Ref. [153]. Similarly, we see that the distinctive model-specific features have very little impact on the total rates and on the size of the QCD corrections.
Before closing this discussion, let us devote one last word to the heavier di-Higgs combinations, shown in Tables 5, 7 and 9 for the B1, B2 and B3 scenarios respectively. The reported total rates span a wide range from σ ∼ O(10 −2 ) fb to σ ∼ O(10) fb and lie in all cases below the light di-Higgs pair predictions. The reason is twofold: i) the relative phase space suppression and ii) the lower initial-parton luminosities involved in the production of these heavier states. The cases in which different rates are obtained for different Higgs  1.67 Table 9: Total Higgs pair cross sections via gluon fusion pp(gg) → h i h j (in fb) for all seven finalstate combinations within the 2HDM, in the same setup as Table 5. The 2HDM parameters are fixed to benchmark B3 in Table 3.

Benchmark B3
qq-channels LO NLO K-factor  Table 10: Total Higgs pair cross sections via quark-antiquark annihilation pp(qq) → h i h j (in fb) for the mixed CP-even/CP-odd and charged Higgs pair combinations within the 2HDM, in the same setup as Table 5. The 2HDM parameters are fixed to benchmark B3 in Table 3. 1.54 Table 11: Total light Higgs pair cross sections via gluon fusion pp(gg) → h 0 h 0 (in fb) within the 2HDM, in the same setup as Table 5. The 2HDM parameters are fixed to benchmarks B4 -B7 in Table 3. B3) can be understood by considering the size of the Higgs couplings in the given scenario.

Differential distributions
As we have seen so far, the most apparent 2HDM imprints rate-wise are of resonant nature. These appear when the on-shell subprocess gg → H 0 → h 0 h 0 is kinematically available and adds to the continuum production. This possibility crucially depends on the heavy Higgs spectrum and is almost insensitive to the distinctive coupling patterns of the 2HDM. The mere observation of an enhancement in the total rate would thus not be undisputed evidence of an underlying 2HDM structure.
Instead, a richer landscape opens up as we move on to the kinematical distributions. These are particularly helpful to track down the model-specific effects which are otherwise smeared once we integrate over the whole phase space. For the remainder of this section we shall focus on the light neutral CP-even Higgs pair channel h 0 h 0 and study representative di-Higgs distributions for the LHC at 14 TeV. Our results are displayed in Figures 4-10, in which we show the light Higgs pair rates as a function of the di-Higgs invariant mass m h 0 h 0 (left panels) and the hardest Higgs transverse momentum p h 0 T (right panels). We concentrate on these two for the sake of conciseness, even though any alternative distribution can be obtained at will, as our setup is fully differential. All histograms are shown at both LO+PS and NLO+PS accuracy, where the NLO results follow from the "loop-improved" approach. The SM prediction (also at NLO+PS) is overlayed for comparison. In the lower subpanels we display the bin-by-bin ratio of the 2HDM NLO+PS results over the SM values.
We begin by pointing out a number of features common to all 2HDM scenarios. Let us first recall that, unlike the SM, there are two types of triangle topologies that contribute in the 2HDM; one of them is linked to the s-channel exchange of the light Higgs boson h 0 , while the other one proceeds through the heavy Higgs exchange H 0 . The first subprocess gg → h 0 * → h 0 h 0 follows the shape of the virtual Higgs boson propagator and would peak aroundŝ ≃ m 2 h 0 . At larger invariant masses, the light Higgs propagator is probed off-shell, which means that for these kinematical configurations, the light Higgs triangles become subleading. The virtual heavy Higgs counterpart gg → H 0 * → h 0 h 0 becomes relevant at larger invariant mass values m h 0 h 0 ∼ m H 0 and may either add to, or partially cancel, the light Higgs triangle amplitudes in this region -depending on the relative signs of the different Yukawa and trilinear Higgs self-couplings. If m H 0 > 2m h 0 , the heavy Higgs is produced on-shell and its resonant peak takes over. As alluded to above, these resonant situations overshadow all other possible new physics effects. In particular, neither the total nor the differential rates are sensitive to a possible enhancement or suppression of the trilinear Higgs self-couplings.
The interplay between the triangles and box topologies, which have different phase space dependence, generates a variety of model-specific signatures which are reflected in the histograms. For instance, a modified trilinear coupling g h 0 h 0 h 0 will mostly reveal itself at low invariant masses. Instead, the boxes are unresponsive to variations in the triple Higgs self-interactions. Shifted Yukawa couplings will typically become more apparent in the larger m h 0 h 0 and p h 0 T regions, through their influence on the box contributions. The latter topologies have a slower decrease with m h 0 h 0 and p h 0 T compared to the triangles, and hence dominate in the hard Higgs tails. Independent changes in the top and bottom-quark Yukawas, which are possible for type-II models, can in addition modify the interference between the top and the bottom loops.
We can see from all histograms that, barring the cases with low-mass H 0 resonances, the majority of the di-Higgs events concentrates on the invariant masses well above the Higgs pair thresholdŝ 4m 2 h 0 . This is after all the reason for the breakdown of the infinite top-mass effective theory (HEFT), which is meant to hold forŝ ≪ m 2 t . The HEFT fails to correctly reproduce the exact distributions not only at large invariant masses but also for moderate values. The transverse momentum distributions are problematic for the same reason [180,190].
We also notice that the QCD NLO effects, while quantitatively important rate-wise, have very little effect on the distribution shapes. Notice that both the LO and the NLO histograms in all figures vary in parallel, which implies a fairly constant K-factor. This is certainly not unexpected, in view of the structure of the NLO QCD corrections, as we have described in Section 2.4. Neither the exchange of virtual gluons between the incoming partons nor the light parton radiation off the initial-state colored legs can significantly influence the leading kinematical features appearing in the m h 0 h 0 and p h 0 T distributions, which rely fundamentally on ii) the relative sizes and signs of the heavy-quark Yukawa and the triple Higgs self-couplings; and iii) the potential enhancement due to a resonant (or close-to-resonant) intermediate heavy Higgs exchange -all these mechanisms being already present at the LO.
Let us now move on to the different model-specific features which can be appreciated in the plots. We begin in Figure  T regions. The dip in the signal right after the resonant peak is due to the interference between the boxes and the heavy Higgs-mediated triangles. The small deviations from the SM away from the resonant peak are caused by the modified trilinear and Yukawa couplings. As expected, at large invariant masses m h 0 h 0 600 GeV the cross section is dominated by the box diagrams. Thereby we explain the 2HDM/SM ratio through the rescaled top Yukawa (g h 0 tt /g SM Htt ) 4 ∼ 0.85. A qualitatively similar situation is encountered for benchmark B2, as shown in Figure 5. Again, the heavy Higgs resonant peak is manifest for m h 0 h 0 ≃ 700 GeV and p h 0 T ≃ 330 GeV. Given that the heavy Higgs mass is now larger, its on-shell single production via gg → H 0 is suppressed by phase space and by the lower gluon luminosity. This accounts for the smaller rates with respect to the B1 scenario discussed above. Close to the light Higgs pair threshold m h 0 h 0 ≃ 2m h 0 , we find an enhanced differential rate with respect to the SM. The interplay between several effects is responsible for this behavior. On the one hand, the heavy Higgs-mediated triangles are small in these bins, while the light Higgs-mediated and the boxes -which again leads to a partial cancellation in this case, because the g H 0 tt and g H 0 h 0 h 0 couplings have the same (negative) sign. In the large m h 0 h 0 tail, instead, the differential rates mostly depend on the box contributions, and hence lie roughly ∼ 20% below the SM yields.
In contrast to the previous cases, the differential distributions in scenario B3 (cf. Figure 6) barely depart from the SM. The reason is twofold: i) the absence of the on-shell heavy Higgs contribution; ii) the SM-like pattern of all light Higgs couplings in the limit ξ = cos(β −α) → 0. In this case, also the trilinear coupling g H 0 h 0 h 0 is extremely suppressed, so that the heavy Higgs-mediated process gg → H 0 → h 0 h 0 barely contributes. The flat 2HDM/SM cross-section ratio is a further indication that no distinctive 2HDM imprints arise in any region of the phase space. The results for benchmark B5, which we do not show explicitly, are also featureless. The SM-like profile in the latter case also results from the fermiophobic nature of the heavy Higgs boson, which cannot couple to the quarks and hence has no influence on the light di-Higgs production.
A variety of non-standard imprints can be appreciated in Figure 7 for benchmark B4. These genuine 2HDM effects have in this case a non-resonant origin and can be better interpreted by analysing the individual components from the different topologies. The latter are shown separately in Figure 8. One first observation is the increased rates right next to the di-Higgs threshold m hh ≃ 2m h 0 . These are primarily caused by the additional heavy Higgs contribution, and are particularly strong in this case not only because of the low heavy Higgs mass m H 0 = 200 GeV, but also due to the unsuppressed trilinear coupling g H 0 h 0 h 0 -in contrast to scenario B3 discussed above. In fact, it turns out that in this region (viz. the first bin of Figure 8 masses and in the boosted Higgs tails (viz. m h 0 h 0 800 GeV and p T 400 GeV) the rates are dominated by the box contributions. Correspondingly, the O(10%)-enhanced Yukawa coupling accounts for the increased 2HDM rates in these tails, which exceed the SM results by up to O(50)%, as expected from the overall rescaling (g h 0 tt /g SM Htt ) 4 ≃ 1.5.
Last but not least, in Figures 9 and 10 we study the impact of a strongly modified bottom-quark Yukawa, as realized by benchmarks B6 and B7. While the two scenarios share a number of similar features, their differences highlight some remarkable properties of the 2HDM. In both cases we observe the expected on-shell peak from the heavy Higgs cascade decay, although here it is softer as compared to the previous resonant scenarios. This milder effect is visible not only in the distributions but also rate-wise and can be explained by the strongly suppressed top-quark Yukawa coupling to the heavy Higgs g H 0 tt (cf. Table 4). One further common trait to both benchmarks is the asymptotic behavior at large di-Higgs invariant masses and in the boosted Higgs tails, where we obtain results very close to the SM predictions. All that said, we can also appreciate some relevant differences. On the one hand, the O(20)% enhanced bottom Yukawa in B6 reinforces the destructive interference between the top and bottom-mediated triangles. In addition, the amplitude of the individual triangle diagrams is pulled down further by the (slightly suppressed) trilinear self-coupling g h 0 h 0 h 0 . Both effects cooperate to reduce the interference between the triangle and the box topologies in the lowest m h 0 h 0 and p h 0 T bins, in such a way that we obtain rates slightly above the SM expectations.
A similar mechanism operates in the B7 scenario, albeit in the opposite direction. In this case, the sign-flipped bottom-quark Yukawa causes the top and bottom-mediated triangles to interfere constructively. As a result, we are left with slightly enhanced triangle amplitudes, which thus reinforce the interference with the boxes. In agreement with this fact, the predicted number of events for m h 0 h 0 m H 0 falls slightly below the SM expectations. Unlike the previous resonant scenarios, in this case we find no dip right after the resonance peak. Instead, we obtain slightly enhanced rates for m h 0 h 0 m H 0 because the fact that (g H 0 h 0 h 0 )(g H 0 tt ) < 0 (cf. Table 4) leads to a constructive interference with the boxes. A dip does appear instead right below m h 0 h 0 m H 0 due to the additional negative sign from the heavy Higgs propagator 1/(s − m 2 H 0 ). At this point, let us mention that we are aware of possible caveats in matching the NLO prediction to the parton shower in the case of enhanced bottom Yukawas. These have been discussed in the context of single Higgs production [31,209] and are related to the heavy-quark mass dependence of the Higgs transverse momentum distributions, which is not known exactly beyond the LO. These issues mostly concern the low Higgs p T region and become more relevant if the bottom-quark Yukawa increases. This is indeed the reason why we do not expect them to matter in our case, as in all of the scenarios that we have explored only the bottom-quark coupling to the heavy Higgs boson g H 0 bb is significantly enhanced. In view of this fact, and since the bottom-mediated effects to the gg-induced heavy Higgs production are very small, we do not find it necessary to adjust our setup to include a more dedicated treatment.

Conclusions
The production of Higgs boson pairs at hadron colliders provides a direct handle on the trilinear Higgs self-coupling. Its direct experimental measurement plays a paramount role in the reconstruction of the Higgs potential and represents a key step towards unraveling the fundamental structure of electroweak symmetry breaking. While the profile of the ∼ 126 GeV resonance discovered at the LHC largely agrees with a SM Higgs profile, there is still experimental room, and also strong theoretical motivation, for a non-minimal Higgs sector. Accurate predictions for Higgs production cross sections and distributions, including reliable estimates on the theoretical uncertainties, are thus necessary to identify and characterize eventual deviations from the SM. The 2HDM constitutes an exemplary framework in which to describe the possible signatures of an extended Higgs sector. The model-specific imprints on the Higgs pair production phenomenology at the LHC can have a resonant or non-resonant origin and may be categorized as follows: i) direct production of additional Higgs boson pairs; ii) virtual or resonant heavy Higgs-mediated contribution to the production of light Higgs pairs; iii) modified Higgs couplings.
In this paper we have examined the production of Higgs boson pairs via gluon fusion in the 2HDM. We have provided predictions for the 14 TeV LHC at NLO accuracy in QCD, matched to the Pythia8 parton shower. We have considered all seven possible finalstate Higgs pair combinations in the model and explored representative up-to-date 2HDM benchmarks. The Higgs bosons in the final state have been kept stable in our simulated event samples. Our calculations have been based on the publicly available tool MadGraph5 aMC@NLO together with the 2HDM@NLO model implementation obtained with the NLOCT package. The NLO QCD corrections to the (loop-induced) gluon fusion channels have been handled via a reweighting procedure which includes the exact one-loop matrix elements. Dedicated codes for this computation are available in [205].
We have reported large QCD corrections, reflected in the sizable K-factors in the ballpark of K ∼ 1.5 − 1.7 (for the gg-initiated channels) and K ∼ 1.2 − 1.3 (for the qq-initiated ones). These QCD effects are dominated by the initial-state light-parton radiation. They remain almost constant over the phase space and barely depend on the model parameters. Once they are taken into account, the theoretical uncertainties on the Higgs pair rate predictions are significantly reduced.
We have examined a variety of characteristic 2HDM features and evaluated their effect on the total Higgs pair rates and kinematical distributions. The underlying model structure influences the light di-Higgs production in different ways i) the virtual (real) heavy Higgsmediated contribution gg → H 0 * (H 0 ) → h 0 h 0 , which may enhance the total rates by up to roughly 2 orders of magnitude above the SM expectations; ii) the diverse possible combinations of enhanced, suppressed and/or sign-flipped Higgs couplings, which lead to increased or reduced rates, particularly apparent in the differential distributions.
One further step in this direction is to combine our current results with the corresponding Higgs branching ratios, and to evaluate the signal-over-background significances for the respective decay modes. Identifying the channels with the best prospects of being measured requires dedicated studies which should involve detector effects and selection cuts to optimize the signal extraction over the backgrounds. Complementary information can be obtained by including the additional double Higgs mechanisms besides the dominant gluon fusion mode, namely vector boson fusion; associated production with gauge bosons; and Higgs radiation off heavy quarks. The ultimate goal is to identify the most promising experimental opportunities for these extended Higgs sector searches at the LHC.