Benchmarking Di-Higgs Production in Various Extended Higgs Sector Models

We present a comprehensive study on Higgs pair production in various archetypical extended Higgs sectors such as the real and the complex 2-Higgs-Doublet Model, the 2-Higgs-Doublet Model augmented by a real singlet field and the Next-to-Minimal Supersymmetric extension of the Standard Model. We take into account all relevant theoretical and experimental constraints, in particular the experimental limits on non-resonant and resonant Higgs pair production. We present the allowed cross sections for Standard Model (SM)-like Higgs pair production and the ranges of the SM-like Yukawa and trilinear Higgs self-coupling that are still compatible with the applied constraints. Furthermore, we give results for the pair production of a SM-like with a non-SM-like Higgs boson and for the production of a pair of non-SM-like Higgs bosons. We find that di-Higgs production in the models under investigation can exceed the SM rate substantially, not only in the non-resonance region but also due to resonant enhancement. We give several benchmarks with interesting features such as large cross sections, the possibility to test CP violation, Higgs-to-Higgs cascade decays or di-Higgs production beating single Higgs production. In all of our benchmark points, the next-to-leading order QCD corrections are included in the large top-mass limit. For these points, we found that, depending on the model and the Higgs pair final state, the corrections increase the leading order cross section by a factor of 1.79 to 2.24. We also discuss the relation between the description of Higgs pair production in an effective field theory approach and in the specific models investigated here.


Introduction
The discovery of the Higgs boson in 2012 by the ATLAS [1] and CMS [2] collaborations at the Large Hadron Collider (LHC) structurally completes the Standard Model (SM) of particle physics. The subsequent investigation of its properties revealed the Higgs boson to be very SMlike [3][4][5][6]. In order to verify that the Brout-Englert-Higgs mechanism [7][8][9][10] is indeed responsible for the generation of elementary particle masses we not only need to measure the Higgs couplings to massive SM particles at the highest precision but we also need to reconstruct the Higgs potential itself to establish experimentally whether indeed its Mexican hat shape is responsible for the Brout-Englert-Higgs mechanism of spontaneous symmetry breaking [11,12]. Moreover, despite the success of the SM, which has been tested very accurately at the quantum level, some open problems remain that cannot be solved by the SM, like e.g. the nature of Dark Matter (DM) or why there is more matter than antimatter in the universe. Extensions of the Higgs sector beyond the SM can provide a DM candidate [13][14][15][16][17][18][19][20][21][22]. Electroweak baryogenesis provides a mechanism to dynamically generate the matter-anti-matter asymmetry if the three Sakharov conditions [23] are fulfilled. This may be possible within extended Higgs sectors with additional scalar degrees of freedom and new sources of CP violation. In summary, the investigation of the Higgs potential itself will provide deep insights into the mechanism underlying electroweak symmetry breaking and into the possible landscape of new physics extensions of the Higgs sector. Understanding the Higgs potential will help us to answer some of the open questions of contemporary particle physics.
In the SM the trilinear and quartic Higgs self-couplings are given in terms of the Higgs boson mass. This is not the case in extended Higgs sectors. Experimentally, the Higgs selfinteractions are accessible in multi-Higgs production, the trilinear Higgs self-coupling in double and the quartic Higgs self-coupling in triple Higgs production. At the LHC the dominant di-Higgs production process is given by gluon fusion into Higgs pairs [24][25][26]. At leading order the process is mediated by heavy quark triangle and box diagrams [27][28][29], implying a small SM cross section of 31.05 fb at FT approx [30] and for a c.m. energy of √ s = 13 TeV. At FT approx , the cross section is computed at next-to-next-to-leading order (NNLO) QCD in the heavy-top limit with full leading order (LO) and next-to-leading order (NLO) mass effects and full mass dependence in the one-loop double real corrections at NNLO QCD. Note that the NLO cross section in the heavy top limit was first presented in [31] supplemented by a large-top-mass expansion [32,33] and the inclusion of the full real corrections [34,35]. The NLO result including the full top-quark mass dependence has been calculated by [36][37][38][39][40] and confirmed in [41,42] by applying suitable expansion methods. The NNLO corrections in the large m t limit have been calculated in [43][44][45], the results at next-to-next-to-leading logarithmic accuracy (NNLL) became available in [46,47]. And the corrections up to next-to-next-to-next-to leading order (N 3 LO) were presented in [48][49][50] for the heavy top-mass limit. Note that these corrections apply to the SM case and cannot necessarily be taken over to extended Higgs sectors where e.g. the bottom loops may play an important role or additional diagrams are involved.
In the SM the process suffers from a destructive interference between the triangle and box diagrams which makes its observation very challenging at the LHC. In extended Higgs sectors SM-like di-Higgs production can be enhanced because Yukawa and trilinear Higgs self-couplings are modified relative to the SM. This can result in altering the interference structure between the various diagrams. Additional Higgs bosons can enhance the cross section resonantly and new colored particles that run in the loops can also lead to an increased cross section. Therefore, extended Higgs sectors not only answer some of the open questions but they can possibly also facilitate access to di-Higgs production. Furthermore, the additional Higgs states lead to a large variety of di-Higgs final states implying a plethora of multi-particle final state signatures. To get a comprehensive picture of which final state signatures in Higgs pair production processes are possible and to be able to give a meaningful guideline to experimentalists in their searches for Higgs pair production, we have to take into account all available experimental and theoretical constraints on beyond-the-SM (BSM) extensions. Since the discovered Higgs boson behaves very SM-like, new physics extensions become increasingly constrained. The trilinear Higgs selfcoupling, however, is not as constrained as single Higgs couplings yet [51,52] and we may still expect some distinctive signatures from Higgs pair production.
The goal of this paper is to investigate Higgs pair production in some archetypical BSM extensions. Imposing all available constraints -also from recent di-Higgs searches -we derive the parameter space of these models that is still allowed. Based on this data set we will derive limits on the involved couplings in Higgs pair production, namely the trilinear Higgs self-coupling and the Higgs-Yukawa coupling, investigate the possible enhancement (or also suppression) of SM-like Higgs pair production and investigate what kind of non-SM-like signatures might appear. We provide benchmark points so that experimentalists can match their derived limits on specific models for further interpretation. Since the experiments derive limits from non-resonant and resonant Higgs pair production and our extended Higgs sector models include both effects, we derive a strategy on how we can apply the available di-Higgs limits on our models. We included for all benchmark scenarios and parameter points fulfilling the applied constraints the information on the resonant part of the cross section, where applicable. Our aim is to give a global and comprehensive overview on Higgs pair production in BSM Higgs sectors.
The models that we consider are both non-supersymmetric and supersymmetric ones. Supersymmetry (SUSY) [53][54][55][56][57][58][59][60][61][62][63][64][65][66] is able to solve many of the open problems of the SM. The nonminimal supersymmetric extension (NMSSM) [67][68][69][70][71][72][73][74][75][76][77][78] solves the little hierarchy problem and more easily complies with the discovered SM-like Higgs mass after inclusion of the higher-order corrections [79]. The Higgs sector consists of two Higgs doublets to which a complex singlet superfield is added so that after electroweak symmetry breaking (EWSB) we have three neutral CP-even, two neutral CP-odd and two charged Higgs bosons in the spectrum. Supersymmetric relations constrain the Higgs potential parameters in a different way than non-SUSY models. Therefore, we also investigate non-SUSY Higgs sector extensions where the trilinear couplings are less constrained from a theoretical point of view. This way we make sure not to miss some possibly interesting di-Higgs signatures. We start with one of the most popular extensions complying with ρ = 1 at tree level, the CP-conserving 2-Higgs doublet model (R2HDM) [80,81] where a second Higgs doublet is added to the SM sector. Incorporating a minimal set of BSM Higgs bosons (five in total, three neutral and two charged ones) allows for resonant di-Higgs enhancement. We additionally take into account the possibility of CP violation (which is required for electroweak baryogenesis) by investigating the CP-violating 2HDM (C2HDM) [81][82][83][84] which consists of three CP-mixed and two charged Higgs bosons. In this case the SM-like Higgs couplings can be diluted by CP admixture, the same happens through singlet admixture. Thus, light Higgs bosons may not be excluded yet because they may have escaped discovery through small couplings to the SM particles. Such a singlet admixture is realized in the next-to-2HDM (N2HDM) [18,85,86]. By adding a real singlet field to the 2HDM Higgs sector the Higgs spec-trum then consists of three neutral CP-even Higgs bosons, one neutral CP-odd and two charged Higgs bosons, allowing for the possibility of Higgs-to-Higgs cascade decays. This is also possible in the C2HDM and the NMSSM. For simplicity, we will focus on the type I and II versions of the R2HDM, C2HDM and N2HDM. With the models investigated in this paper 1 we cover a broad range of interesting new physics and in particular a large variety of possible new physics signatures in di-Higgs production. In turn, we provide guidelines to the experiment. We also take the occasion to confront our models with a simple effective field theory (EFT) approach and investigate to which extent this model-independent parametrisation of new physics, becoming effective at high scales, can describe the effects in our investigated specific UV-complete models.
The paper is organised as follows. In Sec. 2, we briefly introduce our models. In Sec. 3, we present the regions which we scanned for each model and the theoretical and experimental constraints that we take into account. After a brief re-capitulation of the di-Higgs production process through gluon fusion, we explain in detail in Sec. 4 how we apply the experimental limits from resonant and non-resonant di-Higgs searches on our models. We then investigate the impact of the di-Higgs constraints on our parameter sample. We present the distributions of the Higgs mass spectra in the different models and give the ranges for the SM-like Higgs top-Yukawa and Higgs self-couplings that are still allowed after considering all constraints. Subsequently, we present scatter plots for all models showing the cross section values for SM-like Higgs pair production that are compatible with the constraints, and we list the maximum values for resonant and non-resonant Higgs pair production possible in each model. In Sec. 5 we present the maximum values from resonant SM-like di-Higgs production and present the corresponding benchmark points along with their specific features. In Sec. 6 we investigate to which extent di-Higgs production can constrain the parameter values of the models. Section 7 is devoted to the comparison of the EFT description of BSM Higgs pair production with the results in specific UV-finite models. The last two sections, 8 and 9, are devoted to the pair production of a SM-like Higgs together with a non-SM-like Higgs, of a pair of non-SM-like Higgs bosons and to the cascade decays leading to multi-Higgs final states. We also present benchmark points where di-Higgs production beats single Higgs production. We conclude in Sec. 10. In the Appendix we present cross sections for resonant and non-resonant production and discuss the conditions for alignment in the C2HDM and the N2HDM.

The Models
In this section we provide a very brief description of the different models we will be studying, highlighting the diverse scalar spectra of each model as well as the different input parameter sets for each of them.

The Real and Complex 2HDM
The 2HDM is one of the simplest extensions of the SM, where instead of a single Higgs doublet we now have two, carrying identical hypercharges. The model was first proposed by Lee in 1973 [80] to provide an extra source of CP violation via spontaneous symmetry breaking, and has a rich phenomenology (for a review, see [141]). We considered the 2HDM version with a softly broken discrete Z 2 symmetry of the form Φ 1 → Φ 1 and Φ 2 → −Φ 2 . In terms of the two Table 1: Four left rows: The four Yukawa types of the Z2-symmetric 2HDM, stating which Higgs doublet couples to the different fermion types. Five right columns: Corresponding Z2 assignment for the quark doublet Q, the up-type quark singlet uR, the down-type quark singlet dR, the lepton doublet L, and the lepton singlet lR.
SU (2) L Higgs doublets Φ 1,2 with hypercharge Y = +1, the most general scalar potential which is SU (2) L × U (1) Y invariant and possesses a softly broken Z 2 symmetry is given by The Z 2 symmetry is introduced in the model in order to avoid dangerous flavour-changing neutral currents (FCNCs) mediated by the neutral scalar. Since the Z 2 symmetry is extended to the fermion sector, it will force all families of same-charge fermions to couple to a single doublet which eliminates tree-level FCNCs [141,142]. This implies four different types of doublet couplings to the fermions listed in Tab. 1.
The dimension-2 coefficient m 2 12 which breaks the Z 2 symmetry softly, is introduced to allow for the existence of a decoupling limit, in which all scalars other than the SM-like one have very large masses, with suppressed couplings to fermions and gauge bosons. Furthermore, in the case where all parameters in Eq. (2.1) are real we are in the CP-conserving 2HDM, which we will call R2HDM from now on. In the case of complex m 2 12 and λ 5 parameters 2 the model explicitly breaks the CP symmetry. This is the CP-violating version of the 2HDM, called C2HDM. The two complex doublet fields can be parametrised as with v 1,2 being the vacuum expectation values (VEVs) of the two doublets Φ 1,2 . After EWSB three of the eight degrees of freedom initially present in Φ 1,2 are taken by the Goldstone bosons to give masses to the gauge bosons W ± and Z, and we are left with five physical Higgs bosons.
In the CP-conserving case, these are two neutral CP-even Higgs bosons, h and H, where by convention m h < m H , one neutral CP-odd A and a pair of charged Higgs bosons H ± . In the following we will denote h and H by H 1 and H 2 , respectively, in order to standardise the notation for all considered models. The CP-even neutral Higgs mass matrix is diagonalised by a mixing angle α, whereas both the neutral CP-odd and charged mass matrices are diagonalised by a mixing angle β, such that The phases of these two parameters must be such that arg(m 2 12 ) = arg(λ5)/2, otherwise a trivial field rephasing would render the potential real.
In the C2HDM, the three neutral Higgs bosons mix, resulting in three neutral Higgs mass eigenstates H i (i = 1, 2, 3) with no definite CP quantum number and which by convention are ordered as m H 1 ≤ m H 2 ≤ m H 3 . The rotation matrix R diagonalising the neutral Higgs sector can be parametrised in terms of three mixing angles α i (i = 1, 2, 3) as where s i ≡ sin α i , c i ≡ cos α i , and, without loss of generality, the angles vary in the range (2.5) The CP-conserving limit of the C2HDM is obtained by setting α 2 = α 3 = 0 and α 1 = α+π/2 [83]. The shift by π/2 in this limit is necessary to match the usual 2HDM convention. We discuss the alignment limit, where the C2HDM approaches the SM in Appendix B.

By identifying
where v is the SM VEV, v ≈ 246 GeV, and using the two minimisation conditions, the scalar sector of the 2HDM can be described by eight independent input parameters. In the C2HDM, as stated above, the parameters λ 5 and m 2 12 can be complex so that the C2HDM Higgs sector at tree level is described by ten parameters. Notice that it is always possible to perform a basis change to make one of these phases vanish so that we end up with nine independent parameters. In the C2HDM the three neutral Higgs boson masses are not independent. The third neutral Higgs mass is a dependent quantity and is obtained from the input parameters, cf. [143]. We hence choose two of the three neutral Higgs boson masses as input values and calculate the third one. The chosen input masses are called m H i and m H j with H i per default denoting the lighter one, i.e. m H i < m H j . They denote any two of the three neutral Higgs bosons among which we take one to be the 125 GeV SM-like scalar. We furthermore replace the three mixing angles α 1,2,3 by two coupling values of H i and by a matrix element of our rotation matrix. These are the squared H i couplings to the massive gauge bosons V and to the top quarks t, c 2 H i V V and c 2 H i tt , respectively, and the neutral mixing matrix entry R 23 . We furthermore fix the sign of R 13 , sg(R 13 ), to either +1 or -1 in order to lift the degeneracy that we introduce by specifying only the squared values of the H i couplings. This choice of input parameters complies with the input parameters of the program code ScannerS that we will use for our parameter scans as explained below. We hence have the input parameter set v , tan β , c 2

The N2HDM
In the following, we give a brief introduction to the N2HDM and refer to [85] for more details. The scalar potential of the N2HDM is obtained from the 2HDM potential by adding a real singlet field Φ S . In terms of the two SU (2) L Higgs doublets Φ 1 and Φ 2 , defined in Eq. (2.2), and the singlet field, defined as the N2HDM potential is given by The above scalar potential is obtained by imposing two Z 2 symmetries, The first (softly-broken) Z 2 symmetry is the extension of the usual 2HDM Z 2 symmetry to the N2HDM which, once extended to the Yukawa sector, will forbid FCNCs at tree level, implying four different N2HDM versions just like in the 2HDM, cf. Tab. 1. The second Z 2 symmetry is an exact symmetry which will be spontaneously broken by the singlet VEV and as such does not allow the model to have a DM candidate. Other versions of the model choose parameters such that v S = 0 yielding very interesting DM phenomenology, but in the current work we will not consider these possibilities.
After EWSB, we have three neutral CP-even Higgs bosons H 1,2,3 with masses ranked as m H 1 < m H 2 < m H 3 , one neutral CP-odd boson A and a pair of charged Higgs bosons H ± . The physical states H 1,2,3 are obtained from the weak basis (ρ 1 , ρ 2 , ρ S ) by an orthogonal transformation R which is defined by 3 mixing angles α 1,2,3 that are in the same range as in the C2HDM. After exploiting the minimisation conditions, we are left with twelve independent input parameters for the N2HDM. For the scan, we will again replace the three mixing angles α 1,2,3 by the squared H 1 couplings to massive gauge bosons V and the top quarks t, c 2 H 1 V V and c 2 H 1 tt , respectively, and the neutral mixing matrix element R 23 , so that our input parameters read tan β , c 2 Like in the 2HDM, we fix sg(R 13 ) to either +1 or -1 in order to lift the introduced degeneracy through the squared values of the H 1 couplings. The limit of the real 2HDM with an added decoupled singlet field is obtained from the N2HDM spectrum by letting α 2,3 → 0 and α 1 → α + π/2. Again, the shift by π/2 in this limit is necessary to match the usual 2HDM convention. The alignment limit for the N2HDM is discussed in Appendix C.

The NMSSM
As a supersymmetric benchmark model, we consider the Next-to Minimal Supersymmetric SM (NMSSM) [69][70][71][72][73][74][75][76][77][78]. It extends the two doublet fieldsĤ u andĤ d of the MSSM by a complex superfieldŜ. When the singlet field acquires a non-vanishing VEV, this not only solves the µ problem [144] but, compared to the MSSM, it also relaxes the tension on the stop mass values that need to be large for the SM-like Higgs boson mass value to be compatible with the measured 125.09 GeV. Indeed in supersymmetry the neutral Higgs masses are given in terms of the gauge parameters at tree level so that there is an upper mass bound on the lightest neutral scalar which, in the MSSM, is given by the Z boson mass. Substantial higher-order corrections to the Higgs boson mass are therefore required to obtain phenomenologically valid mass values for the SM-like Higgs boson. The additional singlet contribution to the tree-level mass of the lightest neutral Higgs boson shifts its mass to larger values compared to the MSSM prediction, thus no longer requiring large radiative corrections. The scale-invariant NMSSM superpotential that is added to the MSSM superpotential W MSSM reads where for simplicity we only included the third generation fermion superfields, given by the left-handed doublet quark ( Q 3 ) and lepton ( L 3 ) superfields, and the right-handed singlet quark ( t c R , b c R ) and lepton ( τ c R ) superfields. The NMSSM-type couplings λ and κ are dimensionless and taken real since we consider the CP-conserving NMSSM. The Yukawa couplings y t , y b , y τ can always be taken real. The scalar part ofŜ will develop a VEV v S / √ 2, which dynamically generates the effective µ parameter µ eff = λv S / √ 2 through the first term in the superpotential. The second term, cubic inŜ, breaks the Peccei-Quinn symmetry and thus avoids a massless axion, and W MSSM contains the Yukawa interactions. The symplectic product x · y = ij x i y j (i, j = 1, 2) is built by the antisymmetric tensor 12 = 12 = 1. The soft SUSY breaking Lagrangian reads where again only the third generation of fermions and sfermions have been taken into account. The tilde over the fields denotes the complex scalar component of the corresponding superfields. The soft SUSY breaking gaugino parameters M k (k = 1, 2, 3) of the bino, wino and gluino fields B, W and G, as well as the soft SUSY breaking trilinear couplings A x (x = λ, κ, t, b, τ ) are in general complex, whereas the soft SUSY breaking mass parameters of the scalar fields, m 2 Since we consider the CP-conserving NMSSM, they are all taken real. In what follows, we will use conventions such that λ and tan β are positive, whereas κ, A λ , A κ and µ eff are allowed to have both signs.
After EWSB, we expand the Higgs fields around their VEVs v u , v d , and v S , respectively, which are chosen to be real and positive After applying the minimisation conditions, we choose as independent input parameters for the tree-level NMSSM Higgs sector the following, The sign conventions are chosen such that λ and tan β are positive, while κ, A λ , A κ and µ eff are allowed to have both signs. Further parameters will become relevant upon inclusion of the higher-order corrections to the Higgs boson mass that are crucial to shift the SM-like Higgs boson mass to the measured value.

Scans and Theoretical and Experimental Constraints
Our goal is to investigate the landscape of Higgs pair production in extended Higgs sectors by taking into account the relevant theoretical and experimental constraints. In order to do so, we performed a scan in the various parameter spaces of the models and checked each parameter point for compatibility with our applied constraints. In this section, we briefly describe the constraints used in our study. For details, we refer to our previous papers on the different models [85,94,[145][146][147][148][149].
We performed the scans with the help of the program ScannerS [150][151][152] for all models except for the NMSSM. In Tables 2, 3, and 4, we list the scan ranges of the R2HDM, the C2HDM, and the N2HDM, respectively. We give them for the various set-ups with respect to which neutral Higgs boson takes the role of the SM-like Higgs which we will denote H SM from now on. We distinguish the cases "light" where the lightest of the neutral Higgs bosons is SM-like (H 1 ≡ H SM ), "medium" with H 2 ≡ H SM , and "heavy" with the heaviest being SM-like (H 3 ≡ H SM ). In the R2HDM, with only two neutral Higgs bosons, we have the cases "light" (H 1 ≡ H SM ) and "heavy" (H 2 ≡ H SM ) only. Note that in the C2HDM, only two of the three neutral Higgs boson masses are independent input quantities whereas the third one is dependent and computed from the input parameters. Therefore, in the generation of the data points two cases are considered, namely H 1 /H 2 ≡ H SM . The third one is calculated and all three masses are subsequently ordered by ascending mass so that all three set-ups, with either of the H i (i = 1, 2, 3) respresenting the H SM , are covered by the two scans described in Tab. 3. Note also that we restrict ourselves to the type I and II models. For all these models, the R2HDM, the C2HDM and the N2HDM, we apply the same theoretical constraints, which have different expressions for each model, requiring that all potentials are bounded from below, that perturbative unitarity holds and that the electroweak vacuum is the global minimum. In the R2HDM we use for the latter the discriminant from [153] and for the C2HDM the one from [154].
As for experimental constraints, we impose compatibility with the electroweak precision data by demanding the computed S, T and U values to be within 2σ of the SM fit [155], taking into account the full correlation among the three parameters. We require one of the Higgs bosons to have a mass of [156] m H SM = 125.09 GeV ,     and to behave SM-like. Compatibility with the Higgs signal data is checked through HiggsSignals version 2.6.1 [157] which is linked to ScannerS. We furthermore suppress interfering Higgs signals by forcing any other neutral scalar mass to deviate by more than ±2.5 GeV from m H SM . Scenar-ios with neutral Higgs bosons that are close in mass are particularly interesting for non-resonant di-Higgs production as they may have discriminating power with respect to the SM case. The appearance of non-trivial interference effects requires, however, a dedicated thorough study that is beyond the focus of this study and is left for future work. We require 95% C.L. exclusion limits on non-observed scalar states by using HiggsBounds version 5.9.0 [158][159][160]. Additionally, we checked our sample with respect to the recent ATLAS analyses in the ZZ [161] and γγ [162] final states that were not yet included in HiggsBounds. Consistency with recent flavour constraints is ensured by testing for the compatibility with R b [163,164] and B → X s γ [164][165][166][167][168][169] in the m H ± − tan β plane. For the non-supersymmetric type II models, we imposed the latest bound on the charged Higgs mass given in [169], m H ± ≥ 800 GeV for essentially all values of tan β, whereas in the type I models this bound is much weaker and is strongly correlated with tan β. Lower values for m H ± allow, via electroweak precision constraints, different ranges for the masses of the neutral Higgs bosons, which will therefore affect our predictions for di-Higgs production.
In the C2HDM, we additionally have to take into account constraints on CP violation in the Higgs sector arising from electric dipole moment (EDM) measurements. Among these, the data from the EDM of the electron imposes the strongest constraints [170], with the current best experimental limit given by the ACME collaboration [171]. We demand compatibility with the values given in [171] at 90% C.L.
In the NMSSM, we use the program NMSSMCALC [172,173] and compute the Higgs mass corrections up to O((α t +α λ +α κ ) 2 +α t α s ) [174][175][176] with on-shell renormalisation in the top/stop sector. We demand the computed SM-like Higgs boson mass to lie in the range 122 GeV...128GeV which accounts for the present typically applied theoretical error of 3 GeV [79]. We use HiggsBounds and HiggSignals to check for compatibility with the Higgs constraints. Furthermore, we omit parameter points with the following mass configurations for the lightest charginõ χ ± 1 and the lightest stopt 1 , to take into account lower limits on the lightest chargino and the lightest stop mass. The experimental limits given by the LHC experiments ATLAS and CMS rely on assumptions on the mass spectra and are often based on simplified models. The quotation of a lower limit therefore necessarily requires a scenario that matches the assumptions made by the experiments. For our parameter scan we therefore chose a conservative approach to apply limits that roughly comply with the recent limits given by ATLAS and CMS [177,178]. For further details of the Higgs mass computation and of the input parameters as well as their scan ranges, we refer to [176]. Figure 1: Generic diagrams contributing to leading-order C2HDM Higgs HiHj (i, j, k = 1, 2, 3) pair production in gluon fusion.

Gluon Fusion into Higgs Pairs
All our di-Higgs production cross sections through gluon fusion have been computed by adapting the public code HPAIR [179] to the R2HDM, the C2HDM [106,109], the N2HDM and the NMSSM [109,131]. The process is mediated through heavy quark loops already at leading order (LO). Generic diagrams are shown for the example of H i H j (i, j = 1, 2, 3) production in the C2HDM in Fig. 1. The diagrams that contribute are triangle diagrams and box diagrams. The first triangle diagram contains a Higgs boson H k (k = 1, 2, 3) in the s-channel that couples to the final state Higgs bosons H i and H j through the trilinear coupling λ H i H j H k . The box diagrams (third diagram) are proportional solely to the Yukawa couplings of H i and H j to the top and bottom quarks. In the C2HDM, which is CP-violating, we additionally have a Z boson exchange in the s-channel (second diagram) which couples to the CP-mixed Higgs boson final states. This diagram also has to be taken into account in the CP-conserving models when the production of a mixed pair of one CP-even and one CP-odd Higgs boson in the final state is considered. Note that the contribution of the Z boson exchange diagram to the overall cross section is small. Furthermore, the QCD corrections from the SM cannot be taken over here. Our implementation of the BSM models in HPAIR allows us to take the QCD corrections (in the heavy top limit) correctly into account also for this diagram.
In extended Higgs sectors we have several modifications compared to the SM. The additional Higgs bosons H k can lead to resonant enhancement of the di-Higgs cross section compared to the SM in case m H k > m H i + m H j . In Higgs pair production we will call parameter configurations where the resonant rates makes up for a significant part of the cross section "resonant production". For mediator masses of m H k < m H i + m H j resonant enhancement is kinematically not possible. This is a clear case of "non-resonant" production. However, note that, for parameter configurations with m H k > m H i + m H j , the resonance contribution may be very suppressed if the involved couplings are small, the mediator mass is very heavy, its total width is large, or if there are destructive interferences between different diagrams. From an experimental point of view, the cross section would not be distinguishable from "non-resonant" production then. The transition between "resonant" and "non-resonant" is of course fluid. We will address this in detail in the discussion of our application of the experimental limits from resonant and nonresonant di-Higgs production. Further differences from the SM case arise from Higgs-Yukawa and trilinear Higgs couplings deviating from those of the SM Higgs boson and from additional particles running in the loop. The latter is the case for the NMSSM where supersymmetric partners of the top and bottom quark contribute to the loop. An interesting feature is that in the SM we have a destructive interference between the triangle and the box diagrams, implying possible enhancements in extended Higgs sectors where the couplings differ from the SM case. This can be inferred from Fig. 2, where we show the LO Higgs pair production cross section when we vary the SM Higgs top-Yukawa coupling (upper left), the trilinear Higgs self-coupling (upper right) and both couplings (lower) while keeping all other couplings fixed to the SM values. Note, that for the sake of illustration we varied the top-Yukawa coupling in ranges beyond the experimental exclusion limits. 3 We see the destructive interference which becomes largest for λ HHH /λ SM HHH = 2.48. The cross section drops to zero (modulo the small bottom quark contribution) for the top-Yukawa coupling y t = 0 as the Higgs does not couple to the top quarks any more. Note finally that the di-Higgs cross section values through the s-channel exchange triangle diagrams are sensitive to the total widths of the exchanged Higgs bosons as well, that have to be provided for the computation of the cross section. where m HH denotes the invariant Higgs pair mass and where we take the corresponding CT14 pdf set [180] for the LO and the next-to-leading order (NLO) computation. The NLO QCD corrections which are of two-loop order are computed in the limit of heavy loop particle masses. When we explicitly present NLO results we consistently set the bottom-quark mass to zero both at LO and at NLO. For presented K-factors, i.e. the ratio of the NLO to the LO cross section, we consistently set the pdf sets to LO and NLO for the LO and the NLO calculation, respectively.
Note that, in order to keep our scans economic in time, we compute the Higgs pair production cross sections for the results covering the whole scanned parameter space of the models at leading order. The NLO QCD corrections are roughly taken into account by applying a factor of two. This rough approximation works reasonably well for SM-like Higgs pair production [31, 36-40, 106, 131]. When we present specific benchmark scenarios, however, we compute the cross sections explicitly at NLO QCD in the heavy loop particle limit for the considered model by using our adapted HPAIR codes. Be aware, however, that the QCD corrections should be taken with caution for scenarios with large values of tan β. In this case the bottom loop contributions become more important so that the limit of heavy loop particle masses cannot be applied any more.

The Impact of Di-Higgs Constraints
We first want to discuss how we take into account the already available results from the experimental SM-like di-Higgs searches, resonant and non-resonant. Since the experimental limits in the two cases rely on different topologies, the question arises how these limits can be applied in our cases where both resonant and non-resonant production are included in the evaluated cross sections. Moreover, in some of the models we can have more than one resonance decaying into a SM-like Higgs pair. In the following, we will outline our applied strategy.
For the generation of the data we turned off in HiggsBounds the experimental limits from resonant di-Higgs production while applying all other constraints described above. We then computed, for each data point that fulfils the kinematic constraint m H k > 2m H SM , the production cross section σ(H k ), for all possible intermediate resonances H k . The latter has been obtained with the code SusHi v1.6.1 [181][182][183] at NNLO QCD and it includes both production in gluon fusion and in association with a b-quark pair. Note, however, that associated production with a b-quark pair does not play an important role for our scenarios. This production cross section is subsequently multiplied with the branching ratio of the H k decay into H SM H SM . The branching ratio has been obtained for the R2HDM, C2HDM, N2HDM, and NMSSM by using the public codes HDECAY [184][185][186], C2HDM HDECAY [84], N2HDECAY [85,86], and NMSSMCALC [172], respectively. We then compare the cross section σ(H k ) × BR(H k → H SM H SM ) for each possible intermediate resonance 4 with the experimental limits on resonant di-Higgs production. Since these limits are obtained from the LHC run at √ s = 13 TeV, we also compute the SusHi cross sections at this c.m. energy (in contrast to the Higgs pair production cross sections that are all evaluated at √ s = 14 TeV). We took into account the experimental limits that were the stringest ones at the time of the production of our plots. They are given in Refs. [187] for the 4b, [188,189] for the (2b)(2τ ), [190] for the (2b)(2γ), [191] for the (2b)(2W ), [192] for the (2b)(ZZ), [193] for the (2W )(2γ) and [194] for the 4W final state. 5 We furthermore included the recently published results on 4b by ATLAS [196] and CMS [197], (2b)(2τ ) [198] and (2b)(2γ) [199]. Parameter points where at least for one possible resonance H k the experimental limit for any of the final state signatures is exceeded, are rejected. Nevertheless, there is one exception. Since the experimental limits are given assuming narrow resonances, we do not reject points where the ratio of the total width Γ tot (H k ) of H k and its mass m H k exceeds the value (Γ tot (H k )/m H k ) limit = 5%. 6 For illustration, we show the impact of the limits from resonant searches for the N2HDM type I (N2HDM-I) where the lightest CP-even scalar H 1 is the SM-like Higgs boson H SM . The yellow points in Fig. 3 (left) show for the points of our scan in the N2HDM-I parameter space that pass the constraints described in Sec. 3, the single production cross sections of the heavy Higgs boson H 2 , computed with SusHi, and multiplied with the branching ratio into a pair of two SM-like Higgs bosons H 1 . In other words, they represent one of the resonant production modes of a SM-like Higgs pair. The dashed line and the dot-dashed line are the experimental limits obtained from resonant di-Higgs production searches in the 4b final state [196] and the (2b)(2τ ) final state [198], respectively. These limits are now applied on all yellow points. Note, however, that we not only apply them on resonant H 2 but also on resonant H 3 production. The right plot in Fig. 3 shows the situation after applying the aforementioned experimental constraints plus the bounds from (2b)(2γ) [199] and the CMS bounds from 4b [197], which only affect the very low and heavy mass region, respectively, and, due to better visibility, were not added in the plot. All previous experimental results are weaker in the whole heavy resonance mass range and thus automatically satisfied. We see that some of the yellow points above the experimental limits are left over. Here we do not fulfil our criteria of (Γ tot (H k )/m H k ) limit < 5% so that the experimental limits cannot be applied and, thus, no statement about the validity of these points can be made.
The red points in Fig. 3 (left) show for all parameter scenarios passing the constraints of Sec. 3, the cross sections for SM-like Higgs pair production as a function of the mass of the non-SM-like Higgs boson H 2 . As described above, the cross section is calculated at LO and multiplied by a factor two to approximately take into account NLO QCD corrections. The constraints from resonant di-Higgs searches are taken into account by referring to the yellow points. Only those scenarios where the yellow points passed the resonant search limits are retained for the di-Higgs cross sections and result in the allowed red points presented in Fig. 3 (right). The comparison of the left and right plot in Fig. 3 clearly shows that the present di-Higgs searches are already sensitive to the N2HDM parameter space and exclude parts of it beyond single Higgs data constraints.
From the right plot, we infer that there are many points left after application of the resonant search limits. In many of them, the contribution from resonant diagrams is suppressed or kinematically forbidden.
Looking only at the total cross section values, as we do it here, and not at distributions, the sizes of the resonant Higgs pair production cross sections in the suppressed cases are similar to the non-resonant ones or even smaller, so that they cannot be distinguished based on the total rates. Also, invariant mass distributions barely change as we explicitely verified, since the contributions from the resonances are largely suppressed w.r.t. to those from the non-resonant parts. We therefore applied on these points the corresponding non-resonant limits. We come back to this point below. One can also see from this plot that for m H 2 < ∼ 2m H 1 the cross section can be suppressed relative to the SM value. This is to be attributed to destructive interferences between the various diagrams contributing to the di-Higgs cross section. As for the experimental limits from non-resonant searches, they mostly do not constrain our models. The latest results in (bb)(γγ) [199], however, start cutting on the N2HDM-I (with Fig. 4 for type I (left) and type II (right). For all points passing our constraints, we plot the NLO QCD (approximated by a factor 2) gluon fusion SM-like Higgs pair production cross sections for the N2HDM-I (left) and II (right), versus the NNLO QCD gluon fusion production cross section of a heavy non-SM-like Higgs H i that subsequently decays into the SM-like Higgs pair. For H 2 ≡ H SM , we have H i = H 3 , for H 1 ≡ H SM we sum over the two possibilities H i = H 2 , H 3 . From the plot, we can infer that for parameter points where H SM H SM production from resonant heavy Higgs production dominates the di-Higgs process, both cross sections, di-Higgs and single Higgs times Higgsto-Higgs decay, approach each other (see diagonal line in the plot). 7 For the smaller cross sections, resonant production stops playing a significant role and the experimental limits from non-resonant di-Higgs searches can be applied. The most stringent one among the various final states is presently given by the (2b)(2γ) final state [199] 8 , which is visualized in the plots by the horizontal lines. We see that in the N2HDM-I this already cuts into the parameter space so that non-resonant Higgs search constraints start to play a role for certain models. 9 The transition between non-resonant and resonant production is fluid of course. In order to be able to apply the experimental limits we are forced to define a separation between the two cases which is arbitrary. We define a cross section to be resonantly dominated when the single non-SM Higgs production with subsequent decay into SM-like Higgs bosons makes up for more than 10% of the di-Higgs result and accordingly apply the resonant limits. This region separation is shown by the diagonal dashed line in each plot. The shaded region is hence the region where we apply the non-resonant search limits. Apart from the N2HDM-I case, we found that non-resonant searches do not constrain the investigated models at present. The previous definition is arbitrary and a sophisticated experimental analysis taking into account distributions would be required. This is beyond the scope of this paper. Since at present the non-resonant searches are not very sensitive, this approach is good enough for our purpose of drawing an overall phenomenological picture.

Parameter Dependences
The size of the cross section for SM-like Higgs pair production depends on the SM-like Higgs values of the trilinear Higgs self-couplings and of the top and bottom Yukawa couplings. The influence of the latter is only relevant for BSM models in the case of large values of tan β. The cross section value furthermore depends on the masses of all additionally involved non-SM-like Higgs bosons H k , their total widths, their Yukawa couplings and the trilinear Higgs self-couplings λ H k H SM H SM . In this subsection we present the distributions of the mass spectra of the non-SMlike neutral Higgs bosons and discuss the allowed sizes of the trilinear Higgs self-coupling and top-Yukawa coupling of the SM-like Higgs boson in our investigated BSM models. 7 Note that the di-Higgs and single Higgs cross sections are not exactly the same for several reasons. The SusHi single Higgs gluon fusion results are computed at NNLO QCD and √ s = 13 TeV whereas HPAIR for di-Higgs production is run at √ s = 14 TeV for LO QCD and multiplies the result afterwards by a factor of two. (The SM Higgs pair production cross sections at FTapprox [30] differ by 18% at √ s = 13 TeV and 14 TeV, respectively.) Furthermore, different pdfs were used in the computation. Also, HPAIR includes all s-channel Higgs exchange and the box diagrams in the computation of the cross section. The impact of the difference between the cross sections w.r.t. to the application of the experimental limits is negligible, however, as we explicitly verified. 8 Apart from the combined limit which we do not apply. 9 We remind the reader that our di-Higgs cross sections are computed at √ s = 14 TeV while the limits are obtained at 13 TeV. With future experimental results at 14 TeV, even more points will be excluded. For our rough analysis, however, this is a good enough approach. Figure 5 displays the mass distributions of the two neutral non-SM-like Higgs bosons contributing to the cross section for SM-like Higgs pair production in the C2HDM (upper), the N2HDM (middle) and NMSSM (lower) for type I (left) and type II (right) in the former two models. 10 We denote the masses of the heavier one by m ↑ and the one of the lighter Higgs by m ↓ , respectively. In the N2HDM and the NMSSM these Higgs bosons are CP-even, in the C2HDM they have no definite CP quantum number. We found in our scans that the largest freedom in the distributions, after applying all considered constraints, is found in the N2HDM (middle row) where in type I all three neutral Higgs bosons H 1,2,3 can be SM-like (red, green, blue points). The overall lightest Higgs mass spectrum is realized for H 3 ≡ H SM and becomes increasingly heavier if instead H 2 or H 1 are SM-like, respectively. For the latter two cases large mass gaps can occur between m ↑ and m ↓ . In the NMSSM large mass gaps are also possible. In the C2HDM-I, where all three SM-like cases can be realized, this is not the case any more. In particular for H 1 = H SM the masses m ↑ and m ↓ quickly become nearly degenerate with increasing values. This is even more pronounced for the C2HDM-II. Large differences in the mass hierarchies allow for Higgs-to-Higgs cascade decays which is particularly interesting for non-SM-like single Higgs searches or Higgs pair production and is a clear sign of non-minimally extended Higgs sectors in contrast to the R2HDM or MSSM e.g. which feature only two CP-even neutral Higgs bosons. We can hence expect such decays in the N2HDM-I, II, and the NMSSM, and to a less extent in the C2HDM-I because of the more compressed mass spectrum. We further found that in the C2HDM-II only the case H 1 ≡ H SM is realized in our scan after including the constraints. In the N2HDM-II and the NMSSM, H 1 and H 2 being SM-like is still possible but not H 3 . In the R2HDM, not shown in the plots, in type I, for H 1 ≡ H SM , the heavier H 2 mass ranges between 130 GeV and, the upper scan limit of 3 TeV. For H 2 ≡ H SM , the lighter H 1 mass varies between 30 and 122.5 GeV. 11 In type II, m H 2 ranges between 800 GeV and the upper scan limit. Also for all other models, the lightest Higgs masses found to be allowed are betweeen around 30 and 122.5 GeV, where the lower bound is due to our scan limit. In the type-II models the overall mass spectrum is pushed to higher values because of the lower limit on the charged Higgs mass.

Higgs Couplings
In the SM, the triangle and box diagrams contributing to Higgs pair production interfere destructively. Deviations of the trilinear Higgs self-coupling λ 3H SM and Yukawa coupling to top quarks y t,H SM of the SM-like Higgs H SM with respect to the SM values λ 3H and y t,H , respectively, will mitigate this destructive interference. In turn, it will allow for larger cross sections in (the non-resonant) SM-like di-Higgs production. Experimental limits on the trilinear Higgs self-coupling derived from Higgs pair production results often assume the top-Yukawa coupling to be SM-like. 12 We therefore want to answer the following questions: • After applying constraints from single Higgs data, what are the allowed ranges for λ 3H SM and y t,H SM in our models? What is the impact of di-Higgs constraints?
• How do di-Higgs cross sections behave as a function of λ 3H SM and y t,H SM ?
10 For a previous phenomenological investigation of the models, cf. [147]. Due to less strict constraints at that time, the mass spectra were less constrained. 11 Note that we applied a gap of ±2.5 GeV around 125.09 GeV. 12 Results where this assumption has been dropped are given in [51,[199][200][201][202]. • What is the relation of the λ 3H SM and y t,H SM limits to the effective coupling parameters in SM effective field theory (SMEFT)?  red points finally also comply with the di-Higgs constraints. The comparison of the red with the blue points hence shows the impact of the di-Higgs measurements on the coupling values.

Allowed Coupling Ranges and Impact of Di-Higgs Constraints
As we can see from the left plot, in the R2HDM-I, after applying single Higgs constraints, the Yukawa coupling is restricted to values close to the SM, 0.87 ≤ (y R2HDM t,H SM /y t,H ) ≤ 1.07 (blue points). The trilinear Higgs self-coupling is less restricted with values in the range −0.2 ≤ (λ R2HDM 3H SM /λ 3H ) ≤ 1.15 (blue points). Note that the restriction of the trilinear Higgs self-coupling at this stage basically stems from the constraints on the Yukawa coupling from the single Higgs data. Di-Higgs constraints (red points) slightly reduce these values, to about −0.1 ≤ (λ R2HDM 3H SM /λ 3H ) ≤ 1.1. Departures from the SM value in the trilinear couplings come along with non-SM-like Yukawa couplings (cf. wedge region). Turning to the N2HDM-I (right plot) we see that the usual plus single Higgs constraints reduce the Yukawa coupling to about the same values as in the R2HDM-I (blue points). The trilinear coupling on the other hand can vary in a larger negative range with values between -16 and 1 times the SM value. After including the di-Higgs constraints the trilinear coupling range is reduced to about -1 to 1 times the SM value (red points). This is primarily caused by the unitarity constraint. In addition to the perturbative unitarity check performed by ScannerS we applied the following approximate limit on all trilinear Higgs self-couplings combinations in the respective model. 13 We found that this additional constraint only affects the N2HDM. In all other models the inclusion of the constraints through ScannerS left over only scenarios that already respect this unitarity constraint. Besides the latter, we also found that the di-Higgs searches cut on the allowed trilinear Higgs self-coupling values, though to a lesser extent.
The unitarity constraints are responsible for the wedge regions in the plots. Comparing 13 The value is derived by assuming a rough perturbative limit on the Higgs mass of MH = 700 GeV, implying a limit on the trilinear coupling of λ perturb  the shape of the wedge regions in the R2HDM-I and the N2HDM-I we see that an increased precision in the Yukawa coupling will affect the allowed deviation in the trilinear coupling in the N2HDM-I more than in the R2HDM-I. Overall, we find that the trilinear coupling gets more and more restricted but significant deviations are still possible and that they come along with a non-SM-like Yukawa coupling. The present (observed) limits on the trilinear Higgs self-coupling assuming a SM top-Yukawa coupling are -1.0 to 6.6 times the SM trilinear Higgs self-coupling at 95% CL as derived by ATLAS [51] and -3.3 to 8.5 as given by CMS [52]. These experimental sensitivities to the trilinear Higgs self-coupling of the SM start to constrain the parameter space of our models, namely the N2HDM. 14 This can be inferred from Tab. 5 where we list the allowed ranges for the top-Yukawa and the trilinear Higgs self-couplings in our investigated models after applying all described constraints. For all models, due to the single Higgs constraints, the top-Yukawa coupling is bounded to a range of at most ±0.1 around the SM case, with the exception of the NMSSM where it can deviate by up to 17%. 15 The trilinear couplings are less constrained. For the N2HDM-I with H 1 or H 2 being SM-like they are outside the lower ATLAS limit; however, only assuming SM-like Yukawa couplings which is not the case as can be inferred from Fig. 6. Note also that a vanishing trilinear SM-like Higgs self-coupling is also still allowed in some of the models.
There is one caveat to be made on the values given in Tab. 5. These limits have been obtained from the scans in the chosen parameter space with application of all constraints. Hence, they depend on the constraints that we apply, and they also depend on our scanning procedure and sampling. More extended scan ranges and scans adapted to specific parameter regions could possibly find more points and extend these allowed coupling values somewhat. With the given coupling values, however, we are on the conservative side. Furthermore, also note that the C2HDM contains per definition the limit of the R2HDM. This is not reflected, however, in the coupling ranges (and will not be in the plots shown below either). The reason is, that the scan in the C2HDM is performed in different input parameters than in the R2HDM and for finite scan ranges necessarily leads to differences. We explicitly checked that larger R2HDM ranges than in the C2HDM indeed coincide with the CP-conserving limit in the C2HDM and that larger C2HDM ranges compared to the R2HDM are due to truly CP-violating points. We chose not to merge the C2HDM sample with the R2HDM as it allows us to investigate CP-violating effects. As a side remark we add that for the values of our scan the SM-like Higgs boson in the C2HDM-I can still have a CP-violating admixture 16 of up to 16%, 20% and 10% for H 1 , H 2 and H 3 being SM-like, respectively, and of up to 2% in the C2HDM-II with H 1 ≡ H SM .

SM-Like Higgs Pair Production
In this section we discuss the production of two SM-like Higgs bosons H SM in our four considered models R2HDM, C2HDM, N2HDM, and NMSSM. We display, for all valid parameter points passing the applied constraints (including those from di-Higgs searches), the cross sections for SM-like Higgs pair production in the R2HDM, C2HDM, and N2HDM in Fig. 7 and for the NMSSM in Fig. 8. 17 The plots are shown as a function of the heavier of the neutral non-SM-like Higgs bosons, denoted by m ↑ . Red points depict the results where the SM-like Higgs boson is given by the lightest Higgs boson H 1 , green points are those for H SM ≡ H 2 , and blue ones correspond to H SM ≡ H 3 . In the R2HDM, we of course only have red and green points. The overall lower density of the points in the C2HDM compared to the R2HDM is an artifact of the scan, as in the C2HDM one of the neutral masses is not an input parameter but computed from the other two neutral input masses in contrast to the R2HDM so that the same coverage of the parameter space would require a dedicated scan in specific regions, cf. also the remark made above on the comparison of the coupling ranges of the R2HDM and C2HDM. Qualitatively, a larger scan would not change the overall picture, however.
From the plots we can infer from the steep rise of the cross section values, once m ↑ ≥ 2m H SM , the resonant enhancement of the cross sections. In the case of the R2HDM this is due to resonant H 2 production and for the C2HDM, N2HDM or NMSSM this can be resonant H 2 and/or H 3 production. For the latter models, in the case of H SM ≡ H 3 , and for the R2HDM in the case of H SM ≡ H 2 we can only have non-resonant di-Higgs production. We can also see that the cross section values can be suppressed compared to the SM case which is given by the horizontal line in the plots. The R2HDM and C2HDM cross sections approach the SM-like cross section for large mass values (red points) whereas this is not the case for the N2HDM with H 1 ≡ H SM . In 16 It is defined by the rotation matrix element squared |Ri3| 2 , the index i denotes the SM-like Higgs boson in the mass basis, the index 3 the CP-violating degree of freedom in the interaction basis. 17 The values of the signal rates for the various final states after the HSM decays can be obtained by multiplying the cross section by the SM branching ratio into the final state of interest. To a very good approximation the values of the branching ratios of the SM Higgs boson can be applied here as the present LHC results push any SM extension strongly towards the SM limit. the R2HDM, we approach the decoupling limit when m H 2 becomes large. For large H 2 mass the production cross section decreases, and furthermore with increasing m H 2 the trilinear coupling λ H 2 H 1 H 1 that is relevant for resonant H 1 H 1 production goes to zero as it is proportional to cos(β − α) which approaches zero in the decoupling limit. In the N2HDM, for H 1 ≡ H SM , we   can have resonant enhancement both from H 2 and H 3 production so that although m H 3 = m ↑ grows in Fig. 7  In Tab. 6 (left), we compare the maximum cross section values in the different models for all parameter scenarios where m H i =H SM < 2m H SM , in case the SM-like Higgs boson is not given by the heaviest neutral one. For these scenarios resonance enhancement is kinematically forbidden. We have included the NLO QCD correction in the large top-mass limit. In the R2HDM-II and C2HDM-II for H SM ≡ H 1 , in the N2HDM-II and NMSSM for H SM ≡ H 1,2 we did not find any points where m H i =H SM < 2m H SM , because of the constraints. Furthermore, in the C2HDM-II the application of the experimental constraints excludes scenarios with H SM ≡ H 2 or H 3 , and in the N2HDM-II and the NMSSM for H SM ≡ H 3 . In the table, we put dashes for all these cases. We observe that the non-resonant cross sections for all models (where they are available) are above the SM value which is given by σ SM = 38 fb. 18 The largest value is obtained for the R2HDM-I with H 1 ≡ H SM . The enhancements are due to a combination of Yukawa and self-coupling values deviating from the SM. Altogether the cross sections deviate by not more than a factor 1.5 from the SM value in the defined non-resonant regions. It can also be that, although kinematically possible, the di-Higgs production cross section is not much enhanced compared to the SM value through resonance production. This is the case for suppressed Yukawa and/or trilinear coupling values of the s-channel exchanged heavy Higgs boson or because it is very heavy or its total width is large. Also destructive interferences can lead to a suppression. Therefore there is no clear correlation between the computed resonantly enhanced cross sections and the Higgs mass values of the resonantly produced Higgs boson as can also be inferred from Figs. 7 and 8. Only the obtained maximum values show the correlation with the mass value, namely that these maximum possible values decrease with increasing mass of the non-SM-like heavy Higgs as it is expected in models with a decoupling limit. In the suppressed case, from an experimental point of view, these scenarios would be interpreted as non-resonant di-Higgs production. As discussed above, the transition between resonant and non-resonant cases is not trivial. In accordance with our choice of separation given above, we also give the maximum cross section values for those scenarios where the resonance contribution makes up for less than 10% of the total di-Higgs cross section. 19 The NLO QCD values are summarised in Tab. 6 (right) for all models with H SM = H 1,2 (only H 1 for the R2HDM). 20 As expected, they exceed the values given in Tab. 6 (left) but by at most a factor 3.
If not suppressed, the cross sections can be much more enhanced compared to the SM case in the resonance case. This can nicely be inferred from Fig. 9 which depicts the cross section values for SM-like Higgs pair production in the N2HDM-I as a function of the SM-like trilinear Higgs self-coupling ratio for all possible H SM scenarios. The full line in the plot shows the SM di-Higgs production cross section as a function of the trilinear Higgs self-coupling variation, the dot-dashed lines additionally depict the change if we allow for a ±10% variation of the top Yukawa coupling as is still compatible with the data. The comparison with these lines allows us to estimate the resonant enhancement effects beyond the cross section change due to modified trilinear Higgs self-and top-Yukawa-couplings. Close to the SM case λ 3H SM /λ 3H = 1 and y t,H SM /y t,H = 1, we observe enhancements of up to 10 (9) times the SM expectation for    The numbers have to be taken with some care, however, as a more complete scan could find additional values. They should rather be taken as a rough guideline. Note furthermore, that in the NMSSM with H 2 ≡ H SM the resonant production makes up only a small contribution of the full cross section meaning that it is the same point as the one quoted in Tab. 6 (right). The enhancement of the cross section is rather due to coupling deviations of the SM-like trilinear Higgs coupling, while the H 2 top-Yukawa coupling is close to the SM value, and from an experimental point of view this would look like non-resonant di-Higgs production.
For completion, we include in Appendix A separately the information on resonant and nonresonant cross section values corresponding to Figs. 7 and 8. In the transition region, where resonant and non-resonant cross sections are similar in size, they should be taken, however, with a grain of salt.

Benchmarks for SM-Like Higgs Pair Production
In the following, we give the input parameters and basic features of all scenarios listed in Tab. 7. The various models and mass hierarchies have different interesting features: large di-Higgs cross sections, large rates for di-Higgs associated production with gauge bosons but also large triple Higgs production rates from Higgs-to-Higgs cascade decays that are only possible in non-minimal Higgs sectors. Additionally, in the C2HDM, we have the possibility to test CP violation through the measurement of Higgs-to-Higgs or Higgs-to-gauge plus Higgs boson decays. These features will be specified for the various benchmarks. We will also give the dominant decay channels of the non-SM-like Higgs bosons. Due to space limitations, we cannot give all decay channels. They can be obtained, however, from the publicly available programs HDECAY [184][185][186], C2HDM HDECAY [84], N2HDECAY [85,86], and NMSSMCALC [172] by using the input values that we give for each benchmark point. The benchmark points comply with the present and past Higgs searches as explicitly checked by applying the code HiggsBounds where we made sure to take into account the latest search limits. Note in this context that light Higgs bosons still escape detection due to suppressed couplings. Moreover, for some light mass regions there are no experimental analyses available yet. This emphasizes once again the necessity to investigate the low BSM Higgs mass region to constrain extended Higgs sectors. Also the heavy Higgs bosons given in the following benchmarks are escaping all present exclusion limits with rates lying below the current experimental sensitivity. All Higgs pair cross sections have been obtained with our codes based on HPAIR [179], that has been extended to include the various models, at NLO QCD in the heavy loop particle limit. The single Higgs cross sections are calculated at NNLO with the program SusHi v1.6.1 [181][182][183]. The single and double Higgs production cross sections are given for √ s = 14 TeV. Further information on all benchmarks can be provided on request. Please write an email to the authors of the paper in this case.

Resonant SM-like H 1 H 1 Production in the R2HDM-I
This benchmark scenario provides maximum SM-like H 1 H 1 production through resonant enhancement in the R2HDM-I, and is defined by the input parameters given in Tab  With the values given we see that indeed H 2 resonant production is responsible for the enhancement of the cross section as The fact that the HPAIR result for H 1 H 1 production is somewhat lower is due to the finite width of H 2 that is taken into account in the HPAIR calculation when we integrate across the H 2 resonance in the s-channel. We find that for associated H 1 H 1 production with a Z boson we have This leads to a Z + (4b) final state with a signal rate of 84 fb so that this discovery channel for A exceeds A production in the tt final state which amounts to 43 fb.

Resonant SM-like H 1 H 1 Production in the R2HDM-II
This is the benchmark scenario for maximum SM-like H 1 H 1 production through resonant enhancement in the R2HDM-II, with the input parameters given in Tab. 9 (upper) and Higgs pair and single Higgs production information in Tab. 9 (lower). As we are in the type 2 model, the overall Higgs spectrum is heavier than BP1. Here we also have, apart from a large SM-like Higgs    With the single Higgs production cross sections given in Tab. 10 (lower) the large cross section is due to resonant enhancement from both H 2 and H 3 , and we arrive at the following Higgs-to-Higgs, Higgs-to-gauge bosons and Higgs-to-gauge+Higgs rates (5.28) 22 Note that deviations between the full HPAIR result and the SusHi result are to be expected as HPAIR takes into account all s-channel (including the total widths of the s-channel particles) and box contributions and their interferences.
These large rates allow for the test of CP violation through Higgs decays. The decays of H 2/3 each into W W/ZZ and the SM-like Higgs boson pair H 1 H 1 assuming H 1 is CP-even 23 attribute the CP quantum number +1 to H 2/3 . On the other hand, the decays into ZH 1 of H 2/3 identify them to be CP-odd. The simultaneous measurement of all these decays with substantial rates would clearly identify H 2/3 to be states with mixed CP quantum numbers. Note that BR(H 1 → bb) = 0.592.

Resonant SM-like H 2 H 2 Production in the C2HDM-I
Information on this benchmark point is gathered in Tab. 11. It features an overall light Higgs mass spectrum. In contrast to the previous C2HDM scenario, the corresponding rates are too small to allow for the test of CP violation through Higgs decays. This is also the case for the other C2HDM scenario presented in the following. The dominant branching ratios of the

Resonant SM-like H 1 H 1 Production in the C2HDM-II
For the C2HDM-II the maximum resonant production of a SM-like Higgs pair with H SM ≡ H 1 is given by BP5 with the input parameters defined in Tab. 12 (upper) and additional information related to double and single Higgs production in Tab. 12 (lower). Overall, the Higgs spectrum is rather heavy as expected in type-2 models. The dominant branching ratios of the non-SM-like 23 Note that although it is by now clear that the SM-like Higgs cannot be a pure CP-odd state, we are far from excluding large CP-odd components in its Yukawa couplings. In fact, there are so far only direct measurements of thettHSM and τ + τ − HSM couplings. Both ATLAS and CMS [204,205] were able to exclude the purely CP-odd hypothesis in the process pp →tt(HSM → γγ) with 3.9 standard deviations and to establish a 95% CL observed (expected) exclusion upper limit for the mixing angle of 43 o (63 o ). Recently CMS [206] has performed the first measurement of the CP mixing angle of the tau lepton Yukawa coupling, using 13 TeV data, and an integrated luminosity of 137 fb −1 The CP mixing angle was found to be 4 o ± 17 o , allowing to set an observed (expected) exclusion upper limit for the mixing angle of 36 o (55 o ). This angle is defined as arctan(b/a), if the generic Yukawa coupling is written as a + ibγ5.   The dominant contribution to resonant production stems from H 2 . More specifically, we have The sum of the two resonant contributions makes up for about half of the total cross section. The remaining part is given by non-resonant production which is enhanced compared to the SM case because the trilinear coupling between three SM-like Higgs bosons deviates from the SM case so that the destructive interference between box and triangle diagrams is not effective.

Resonant SM-like
leading to a Z + (4b) rate of 288 fb.

Resonant SM-like H 2 H 2 Production in the N2HDM-I
Here we have the scenario BP7 where H 2 ≡ H SM . The point is interesting not only because of its large H 2 H 2 cross section but also because it allows for significant production of a Higgs pair with a Z boson in the final state and it leads to a large pair production rate for the SM-like Higgs together with a non-SM-like lighter one. All ZH 1 H 2 , ZH 1 H 1 and H 1 H 2 production channels have significant rates in the 4b (plus Z) final states. The input parameters are listed in Tab. 14 (upper) together with double and single Higgs production related information in Tab. 14 (lower). The CP-even Higgs bosons are rather light.
This leads to substantial production rates for Higgs pair plus gauge boson final states, namely

Resonant SM-like H 1 H 1 Production in the NMSSM
This NMSSM benchmark point features, besides a large H 1 H 1 production cross section, large rates for ZH 1 H 1 and triple H 1 production. As stated above, in the NMSSM, the Higgs boson masses are computed from the input parameters of the model and higher-order corrections are important to shift the SM-like Higgs mass to the observed 125 GeV. We have computed these masses using the new version of NMSSMCALC which includes the two-loop Higgs mass corrections at O((α λ +α κ +α t ) 2 +α t α s ) [176]. In Tab. 17, we list all input parameters for this benchmark point 24 24 In accordance with the SUSY Les Houches Accord (SLHA) [207,208] the soft SUSY breaking masses and trilinear couplings are understood as DR parameters at the scale MSUSY = √ mQ 3 mt R which is also the renormalisation scale used in the computation of the higher-order corrections to the Higgs masses. The soft SUSY breaking parameters of the first two generations are not listed as their influence is negligible. The remaining SM input parameters are given in Ref. [176]. The Higgs mass corrections are computed with on-shell renormalisation in the top/stop sector and on-shell renormalised charged Higgs mass, cf. [176] for details. For completeness, we also list the corresponding value of A λ . that are required by NMSSMCALC to compute the Higgs spectrum. Higgs masses, mixing angles and the total widths of the Higgs bosons, are given in Tab. 18. The table contains additional information related to double and single Higgs production. The given CP-even mixing elements h ij (i, j = 1, 2, 3) comply with the SLHA [207,208]       The dominant branching ratios of the non-SM-like Higgs bosons are given by Since the branching ratio BR(H 3 → H 2 H 2 )= 6.5 × 10 −5 is tiny, it is not resonance H 3 production that enhances the cross section but rather the SM-like trilinear coupling deviation from the SM as stated above.

Constraining Model Parameters
With the SM-like Higgs pair production cross section in the non-resonant case being three orders of magnitude smaller than single Higgs production, deriving constraints on the parameter spaces of the models from di-Higgs production may not be very efficient. This picture changes of course in case of resonance enhancements. To get a rough picture of what can be learnt from di-Higgs production, we present a few selected heat plots for the SM-like Higgs pair production cross sections as a function of relevant model parameters.
In the non-resonant case, it is the top Yukawa and trilinear Higgs self-couplings of the SMlike Higgs boson that determine the size of the SM-like Higgs pair production cross section. In Fig. 10, we show as colour code the size of the cross section for SM-like Higgs pair production normalized to the SM value in the N2HDM-I for the case where H 3 ≡ H SM , as a function of its top-Yukawa and trilinear Higgs self-coupling normalized to their SM values, respectively. By choosing H 3 ≡ H SM , we make sure that the contributions from lighter s-channel Higgs boson exchanges to the cross section are subdominant and that it is the SM-like couplings that determine its size. As can be inferred from the plot, the size of the cross section is mostly insensitive to the value of the trilinear Higgs self-coupling in the limited range that is still allowed for it, while it shows a strong dependence on the top-Yukawa coupling. Taking off the resonance contribution for the case where H 1 ≡ H SM and where the allowed range for the trilinear Higgs self-coupling is larger we find a significant dependence of the cross section on the trilinear coupling. Let us investigate if di-Higgs production can possibly contribute to constraining the parameter space of the model. For this, we resort to the simpler R2HDM whose tree-level Higgs couplings are described by only two mixing angles. In Fig. 11, we depict in the tan β-cos(β − α) plane, through the colour code, the NLO QCD-corrected (by including a factor 2) di-Higgs cross section normalized to the SM value for the R2HDM-I, where the lighter Higgs H 1 ≡ H SM . On the left, the complete cross section is plotted. On the right plot, we consider only points where the resonant contribution from H 2 → H 1 H 1 makes up for less than 10% of the total cross section. More specifically, we only include points where σ NNLO From an experimental point of view this would correspond to non-resonant H 1 H 1 production (according to our definition). From the right plot, we could infer that cross section values deviating from the SM value allow us to constrain cos(β −α). The true picture is more complex though, as shown by the left plot: the possible resonance enhancements in BSM di-Higgs production also allow for larger cross sections very close to the SM Higgs alignment limit cos(β − α) = 0. Indeed we found by comparing the constraints on cos(β − α) before and after applying the di-Higgs constraints that the impact of di-Higgs constraints is vanishingly small.
With increasing complexity of the Higgs sectors, the superposition of the various Higgs contributions to the cross sections make it more and more complicated to derive conclusive statements on the various trilinear and Yukawa couplings and require the combination of different cross sections to extract all involved couplings (cf. [11] for a discussion in the MSSM).

Effective Field Theory versus Specific Models
In this paper so far, we discussed mainly Higgs pair production in specific models. In this section, we switch gears towards a different approach to describe new physics effects. This is given by the effective field theory (EFT) framework where BSM physics is expected to appear at some high new physics scale Λ. In the linear approach called SMEFT [209][210][211][212] new physics is formulated as a power series in the dimensionful parameter 1/Λ. The non-linearly realized EFT, on the other hand, can be viewed as organised by chiral dimension [213][214][215][216][217][218][219][220][221][222][223][224][225][226]. If we choose to describe Higgs pair production in the EFT approach this means that effects from additional non-SM-like light Higgs bosons cannot be described. 25 A discussion of the higher-dimensional operators relevant for Higgs pair production can be found in [230][231][232][233][234]. The QCD corrections in the infinite top mass limit, m t → ∞ have been provided at NLO QCD in [97] and also extended to the CP-violating case in [106]. At NNLO QCD they have been calculated in [235]. The authors of [236] presented the NLO QCD corrections including the full top quark mass effects in a non-linearly realized EFT. An interface with POWHEG [237][238][239] has been provided in [125].
For the models that we considered, only the new physics operators that modify the trilinear Higgs self-coupling and the Yukawa coupling are relevant. The induced effective couplings of one or two Higgs bosons to two gluons only appear in the NMSSM, where integrating out heavy stops and sbottoms in the Higgs-to-gluon loop couplings would induce such couplings. We neglect that effect in the present discussion for simplicity, by setting the associated couplings (c g and c gg in the notation of the non-linear Lagrangian of Ref. [97]) to zero. In the R2HDM, C2HDM and N2HDM, these couplings do not appear as long as we do not include additional heavy coloured particles beyond the SM. Furthermore, we do not consider effects from the chromomagnetic operator as they are of different order in the chiral expansion. 26 Finally, integrating out a possible heavy Higgs boson exchange in the s-channel leads to an effective two-Higgs-two-fermion coupling. Denoting by c 3 the trilinear coupling modification and by c t the top-Yukawa coupling modification with respect to the SM and by c tt the effective two-Higgs-two-fermion coupling coefficient, i.e. adapting the notation of [97], our considered correction ∆L non-lin to the SM Lagrangian reads, where h denotes the physical Higgs boson. Since the single and double-Higgs coefficients c t and c tt to the top-quark pair are taken to be independent, we adapt the non-linear effective Lagrangian approach here. In SMEFT, they are correlated (as well as c g and c gg ). In Fig. 12, we show the generic diagrams that contribute to our EFT approach to Higgs pair production and indicate the EFT coupling modifiers. In the notation of [29] and [97], we have the following Figure 12: Diagrams contributing to Higgs pair production in the EFT approach (with cg = cgg = 0 and neglecting the chromomagnetic operator). The blue, red and green blobs denote the modified Higgs trilinear, Higgs top-Yukawa and the new two-Higgs-two-top-quark couplings, respectively. 25 For an extension of the EFT approach to include an extended particle content, an EFT for the 2HDM, cf. [227]. Also for composite Higgs models a concrete model with two Higgs doublets has been proposed, cf. [228,229]. 26 For a discussion on the chromomagnetic operator, see [236].
two-Higgs-two-top quark coupling : Here g φ t (α i , β) denotes the dimensionless function of the mixing angles α i and β that specifies for each model under consideration the modification of the Yukawa coupling of a Higgs boson φ of the model with respect to the SM Yukawa coupling. The function g H SM H SM H SM 3 (p i ) denotes the dimensionful trilinear coupling of three SM-like Higgs bosons H SM in our BSM model that in the SM case would approach 3M 2 H SM /v. We denote by p i the various parameters on which the trilinear coupling depends in the respective model. The third matching relation to c tt is obtained by assuming a possible heavy Higgs H k s-channel exchange (cf. the first diagram in Fig. 1 for H k = H SM ) where the mass of the exchanged Higgs boson, denoted by m H k , is very large. By g H k H SM H SM 3 we denote the corresponding H k trilinear coupling to two SM-like Higgs bosons. Note that, in non-minimal models, we would have two such contributions. Hence k max = 1 in the R2HDM and 2 in the C2HDM, N2HDM and NMSSM. 27 Table 5 gives us an overview of the c t and c 3 values that are allowed by the bulk of the parameter points.
We have chosen a few benchmark points from our samples in order to investigate the validity of the EFT approach. In Tab. 21, we present the benchmark point SMEFTBP1 for the R2HDM-II with the heavy scalar Higgs M H 2 mass above 1 TeV so that the EFT approach should be justified. The corresponding SMEFT coupling coefficients are   values together with the R2HDM-II and the SMEFT results for LO H2H2 production including the resonance contribution.
the results in the SMEFT approach, as well as the ratios of these two cross sections. Note that  Note that in this subsection all gluon fusion cross sections are given at LO. From the second line, we read off that in our scenario the SMEFT approach approximates the cross section in the full model by only 86% for a Higgs mass m H 2 of the order of 1 TeV. When we turn off the H 2 resonance and compare the results with the one in the SMEFT approach where we accordingly set c tt = 0 we get Both cross sections agree as expected in contrast to the case with the resonance included. Since in the di-Higgs cross section we integrate √ s in the s-channel exchange across the resonance, the SMEFT approach is not a good approximation. We want to investigate the minimum mass values from which the SMEFT rate is close to the full R2HDM result. For this, we gradually increase m H 2 and calculate for the corresponding c tt and trilinear Higgs coupling values the SMEFT cross section and also the full R2HDM cross section. 28 The values are given in the third and fourth line in Tab. 22. We clearly see that with increasing m H 2 , and hence decreasing contribution of the resonance to the cross section, the SMEFT and the full R2HDM results approach each other. Starting from about m H 2 = 1200 GeV the deviation is less than 10%, continuously decreasing with increasing m H 2 .
We perform the same investigation but now for the N2HDM-I with H 1 = H SM where we have two resonance contributions. As benchmark point SMEFTBP2 we take the N2HDM benchmark point BP6 given in Tab. 13. For convenience, we repeat the input parameters in Tab. 23  in GeV together with the N2HDM-I and the SMEFT result for H1H1 production at LO including the resonance contribution. The cross sections in the two approaches agree as expected. Starting from our original N2HDM-I scenario we then gradually increase the m H 2 mass which hence changes c H 2 tt and thereby c tt in order to investigate when the SMEFT result starts to reproduce the full result. The corresponding values are given in Tab. 24 from the third line onwards. The SMEFT and the N2HDM results start to deviate by less than 10% for H 2 masses above about 440 GeV. This agreement also depends on the total width Γ H 2 of the s-channel resonance. Keeping e.g. the total width at the value Γ H 2 = 0.075 GeV (corresponding to the mass m H 2 = 269 GeV) the agreement between full and EFT approach within 10% would be reached around M H 2 = 465 GeV. 29 We hence find in this scenario with two possible heavy resonances that the Higgs mass limit, from which the SMEFT approach starts to approximate the full result, ranges at lower values.
The investigation of these two benchmarks with additional Higgs bosons has shown that first, the results calculated in the full theory and in the EFT approach can differ severly. Second, the agreement between the full theory and the EFT approach depends on the masses of the additionally present Higgs bosons and their total widths. The total width plays an important role when we integrate across the resonance in the s-channel within HPAIR. A priori, one cannot predict to which extent the full theory and the EFT approach agree, as this depends on the parameters of the full model.

Mixed Higgs Pair Final States -H SM + Φ
In all presented models it is possible to produce Higgs pair final states that consist of a SMlike Higgs boson H SM plus a non-SM-like one Φ. In the R2HDM, N2HDM and the NMSSM the non-SM-like Higgs boson can be a scalar or a pseudoscalar. In the C2HDM it would be a CP-mixed state. There is a plethora of final states with substantial production rates possible depending on the major decay modes of the non-SM-like Higgs boson. They range from pure multi-fermion, multi-photon, mixed fermion-photon to multi-Higgs boson, multi-gauge boson, mixed Higgs-plus-gauge boson or mixed fermion-plus-Higgs or gauge boson final states. In the NMSSM, we can additionally have supersymmetric particles in the final state which we will not discuss in this paper. In the following, we present some selected benchmark points. All of them have the common feature that their rates exceed 10 fb at NLO. We have many more final state signatures beyond the presented ones that we can provide on request.

The (bb)(bb) Final State
The largest 4b final state from mixed Higgs pair production is found in the N2HDM-I. We present a benchmark point for  The same channel in the N2HDM-II reaches 18 fb. We present in Tab. 26 for the 4b final state the maximum NLO QCD values in the heavy loop particle limit together with the respective K-factors, for these channels and other mixed Higgs pair combinations in the models under investigation. In the case of resonantly enhanced di-Higgs cross sections, we give the mass of the "resonant" Higgs boson and the corresponding resonant production cross section. It is obtained by calculating the production cross section of the "resonant" particle with SusHi at NNLO QCD and subsequently multiplying it with the branching ratio into the investigated Higgs pair final state.

The (bb)(W W ) Final State
If the SM-like Higgs boson decays into W W then the rates are easily obtained from those of the previous subsection in the 4b final state by multiplying them with BR(H SM → W W )/BR(H SM → bb) ≈ 1/3. However, we can also have the case that the non-SM-like Higgs boson decays into W W , which are the benchmark points that we list here. The maximum rate (at NLO) is obtained for    Table 28: Maximum rates at NLO QCD in the (bb)(W W ) final state for different mixed Higgs pair final states in the investigated models with the non-SM-like Higgs decaying into W W ; the corresponding K-factor is given in the last column. In case of resonantly enhanced production, we give in the third and fourth column, respectively, the mass of the resonant Higgs boson and the resonant cross section as defined in the text. The fifth column contains the mass of the non-SM-like final state Higgs boson. More details on these points can be provided on request.

The (bb)(tt) Final State
As the SM-Higgs decay into tt is kinematically forbidden, it is always the non-SM-like Higgs that decays into tt. We find the maximum rate for   Table 30: Maximum rates at NLO QCD in the (bb)(tt) final state at NLO for different mixed Higgs pair final states in the investigated models with the non-SM-like Higgs decaying into tt; the corresponding K-factor is given in the last column. In case of resonantly enhanced production, we give in the third and fourth column, respectively, the mass of the resonant Higgs boson and the resonant cross section as defined in the text. The fifth column contains the mass of the non-SM-like final state Higgs boson. More details on these points can be provided on request.

Multi-Higgs Final States
In non-minimal Higgs models like the C2HDM, N2HDM, and NMSSM we can have multi-Higgs final states from cascade Higgs-to-Higgs decays. In the production of a SM-like plus non-SM-like Higgs final state, H SM Φ, we found that both the Higgs-to-Higgs decay of the SM-like Higgs or the non-SM-like one can lead to substantial final state rates. The largest NLO rates that we found above 10 fb, in the multi-Higgs final state, are summarised in Tab. 31. In the C2HDM, we did not find NLO rates above 10 fb. We maintain the ordering of particles with regards to their decay chains, so that it becomes clear which Higgs boson decays into which Higgs pair. We give the rates in the (6b) final state as they lead to the largest cross sections for all shown scenarios. In the following, we highlight a few benchmark scenarios from the  229 119 518 26 Table 31: Upper: Maximum rates for multi-Higgs final states given at NLO QCD. The K-factor is given in the last column. In the third and fourth column we also give the mass values mΦ 1 and mΦ 2 of the non-SM-like Higgs bosons involved in the process, in the order of their appearance. Lower: In case of resonantly enhanced production the mass of the resonantly produced Higgs boson is given together with the NNLO QCD production rate. More details on these points can be provided on request.

Non-SM-like Higgs Search: Di-Higgs beats Single Higgs
In the following we present N2HDM-I and NMSSM scenarios with three SM-like Higgs bosons in the final states with H 1 being SM-like and with NLO rates above 10 fb. These benchmark points are special in the sense that the production of the non-SM-like Higgs boson H 2 from di-Higgs states beats, or is at least comparable, to its direct production. 30 This appears in cases where the non-SM-like Higgs is singlet-like and/or is more down-than up-type like. The latter suppresses direct production from gluon fusion. The former suppresses all couplings to SM-like particles. In these cases the heavy non-SM-like Higgs boson might rather be discovered in the di-Higgs channel than in direct single Higgs production.  Note that the H 2 branching ratio into (bb) is tiny. The second lightest Higgs boson H 2 has a significant down-type and large singlet admixture but only a small up-type admixture so that its production in gluon fusion is not very large 31 and also its decay branching ratios into a lighter Higgs pair are comparable to the largest decay rates into SM particles. In this case, the non-SMlike Higgs boson H 2 has better chances of being discovered in di-Higgs when compared to single Higgs channels. Note, that the W bosons still need to decay into fermionic final states where additionally the neutrinos are not detectable so that the H 2 mass cannot be reconstructed.
The input parameters for the first NMSSM scenario that we discuss here are given in Tab. 33. We also specify in Tab. 34 the parameters required for the computation of the Higgs pair production cross sections through HPAIR.  Since H 2 is rather singlet-like, its production cross section through gluon fusion is small and also its decay branching ratios into SM-final states. The gluon fusion production cross section amounts to σ NNLO (H 2 ) = 13.54 fb .
(8.62) 31 The production in association with b quarks is very small for the small tan β value of this scenario.    The SM-like Higgs is H 1 and the K-factor for the NLO QCD production of H 2 H 2 is 1.82. Also in the NMSSM and C2HDM we can have multiple Higgs production but the rates are below 10 fb after the decays of the Higgs bosons. In the N2HDM, we can even produce up to eight Higgs bosons in the final states but the rates are too small to be measurable.

Conclusions
In this paper, we have performed a comprehensive analysis of Higgs pair production in some archetypical BSM models, namely the R2HDM, the C2HDM, and the N2HDM as non-SUSY representatives, and the NMSSM as a SUSY model. After applying the relevant theoretical and experimental constraints, in particular limits from non-resonant and resonant di-Higgs searches, we explore the ranges of the parameter spaces of these models that are still allowed. We find that while the SM-like Higgs top-Yukawa couplings are constrained to within about 10% of the SM model value, there is still some freedom on the trilinear Higgs self-coupling. In particular, zero values for the SM-like trilinear Higgs self-coupling are still allowed in all models. Interestingly, the experimental searches start to constrain the trilinear couplings of the N2HDM. In general, in order to derive limits on the couplings both resonant and non-resonant searches will be required. Overall, the delineation of the parameter space from di-Higgs production is difficult as in BSM models we have the Yukawa and trilinear couplings of various Higgs bosons involved and also their total widths play a role in the size of the cross section.
As for the maximum possible sizes of the resonantly enhanced cross sections for SM-like Higgs pair production, we find that they can be quite different across the models studied, so that the cross section value itself might exclude certain models provided it exceeds a specific limit. We presented benchmark scenarios for the maximum cross sections. They not only feature cross sections that can exceed the SM value by up to a factor of 12 but they are also interesting because they can lead to measurable rates of triple Higgs production or the production of a Higgs boson pair in association with a Z boson. In the C2HDM, we presented a scenario where the simultaneous measurement of Higgs-to-bosons decays would allow for the test of CP violation.
We also investigated to which extent an EFT approach can reproduce the Higgs pair results if the resonances are rather light. Moreover, we found that the value of the resonance mass from which on the EFT approach starts to approximate the result in the specific UV-complete model by better than 10% depends on the benchmark scenario itself and that the total width of the intermediate resonance plays an important role. This again emphasizes the importance of investigating specific models besides a more general EFT approach in order to get a complete picture of the new physics landscape.
We also presented benchmark points for the production of a SM-like Higgs together with a non-SM-like one, leading to a plethora of different final state signatures. The presented benchmark points represent those with significant rates in certain final state signatures, namely (bb)(bb), (bb)(W W ), and (bb)(tt). But also multi-photon final states can have important rates as shown by our results, in particular in non-SM-like Higgs pair final states. We highlighted additionally scenarios that involve further Higgs-to-Higgs decays, of the mixed SM-like plus non-SM-like final state, leading to three-Higgs final states with significant rates. On top of that, three benchmark scenarios were presented where new heavy Higgs bosons might rather be discovered through resonant di-Higgs production than from direct single Higgs production. This is the case for the singlet extended models N2HDM and NMSSM, where the couplings of the searched heavy Higgs boson to SM particles are suppressed due to large singlet admixtures but not its trilinear coupling to other Higgs bosons. Finally, we gave a short overview of the maximum possible rates in the production of two non-SM-like Higgs bosons and also for the production of multiple Higgs bosons.
For all benchmark points, we provided the di-Higgs production cross sections at NLO QCD in the heavy loop particle limit. We found that the K-factors range betwee 1.79 and 2.24 depending on the model and the final Higgs pair state.
We collected a large amount of viable parameter points in the R2HDM, C2HDM, N2HDM, and NMSSM with interesting features in the context of Higgs pair production and necessarily had to restrict ourselves on specific scenarios to highlight some prominent features. We emphasize, however, that we can provide benchmarks with specific features on request and invite the readers to contact us.
With this work, we hope to have given a close to complete and comprehensive overview of possible signatures in di-Higgs or even multi-Higgs production in representative BSM models and what can be learnt from them. It can be a starting point for further investigations in many different final state signatures based on the data sample that we have generated. Having a guideline at hand of what can be expected may help us find our way through the vast new physics landscape and get deeper insights in the mechanism behind electroweak symmetry breaking. Ultimately helping us to answer some of our most pressing open questions in the world of elementary particle physics. acknowledges support the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762 -TRR 257. J. El F. would like to thank the Abdus Salam International Centre for Theoretical Physics (ICTP) for hospitality and financial support where part of this work has been done. R.S. and P.F. are supported by CFTC-UL under FCT contracts UIDB/00618/2020, UIDP/00618/2020, and by the projects CERN/FISPAR/0002/2017 and CERN/FIS-PAR/0014/2019. P.F. is supported by the project CERN/FIS-PAR/0004/2019.

A Resonant and Non-Resonant Production Cross Sections
In this section, we re-present the cross sections of Figs. 7 and 8 by, first, showing all the points that have a resonant contribution and, second, displaying only points that are considered nonresonant by the definition given in Sec. 4.2. Notice that the points displayed for each part are not complementary in the full sample, i.e. some point for which we show its resonant contribution might be defined as non-resonant and appear in the non-resonant plots as well. This is because a point being defined as dominantly non-resonant (by our definition) is not exempt from being compared to experimental resonant constraints. On the other hand, it makes no sense to compare points where resonance production is dominant with experimental non-resonant limits, hence only points that we define as non-resonant ones can be compared to non-resonant experimental limits. To be specific, points with resonant contribution means parameter points where at least one of the heavier neutral Higgs bosons 32 has a mass large enough to decay into a pair of two SMlike Higgs bosons. And the resonant production cross sections are obtained by calculating with SusHi the NNLO QCD production cross section of the "resonant" Higgs boson and subsequently multiplying it with its branching ratio into the two SM-like Higgs bosons. The points included in the non-resonant plots on the other hand are all those scatter points of Figs. 7 and 8 where resonant production is kinematically not possible, or where the resonant production cross section accounts for less than 10% of the total di-Higgs cross section. Figure 13 shows the resonant production cross section values for the R2HDM for the points presented in Fig. 7 (upper). The color code denotes the ratio of the total width of the resonant Higgs boson and its mass.
In Fig. 14 we display the resonant production cross sections for the C2HDM for the different cases w.r.t. to which of the H 1,2 is the SM-like Higgs boson. Note that in the C2HDM T2 we do not have scenarios compatible with all constraints with resonantly enhanced production in case H 2 is the SM-like Higgs boson. In case H 1 is the SM-like Higgs boson, both H 2 and H 3 can in principle lead to resonant enhancement. We therefore show in separate plots their resonant cross sections as a function of their mass. The resonant production cross sections for the N2HDM for the different cases w.r.t. to which of the H 1,2 is the SM-like Higgs boson are shown in Fig. 15 and those for the NMSSM in Fig. 16.
In Fig. 17 we display the di-Higgs cross sections for points considered non-resonant by our definition, for the R2HDM, C2HDM, and N2HDM type 1 and type 2. We show them as a function of the trilinear Higgs self-coupling of the SM-like Higgs boson in the respective model normalized to the value of the trilinear Higgs self-coupling of the SM, λ 3H SM /λ 3H . In Fig. 18 we show the corresponding NMSSM plot. We also include in the plots as a full line the change of the SM Higgs pair production cross section as a function of λ 3H SM /λ 3H . The dashed lines show its change if additionally the top-Yukawa coupling is varied by 10% away from the SM value. From these plots we can infer the importance of interference effects. We see e.g. in Fig. 17 (lower left) for the N2HDM-I with H 2 ≡ H SM green points that are well below the full and dashed lines. The suppression of the N2HDM cross section cannot be caused by the variation of the trilinear or top-Yukawa coupling away from the SM values (indicated by the full and dashed lines). It is caused by the negative interference between the triangle diagram contributions of H 1 and H 2 , as we explicitely verified.

B Alignment Limit in the C2HDM
In the following, we derive the limits that are required to achieve the alignment limit in the C2HDM. In terms of the matrix elements of the mixing matrix, defined in Eq. (2.4), the H i (i = 1, 2, 3) Higgs couplings to massive gauge bosons (V = Z, W ) and to fermions f in the C2HDM-I and C2HDM-II, read respectively g C2HDM-I/II where u denotes up-type quarks, d down-type quarks, l leptons and g SM HXX are the corresponding SM couplings of the SM Higgs H to XX.
We define the conditions to get alignment between the C2HDM and the SM for the SM-like  with m H i = m H . The implications of these limits for the various possibilities of H 1 , H 2 , or H 3 being SM-like for the mixing angles are as follows:  as required in the alignment limit. For the latter, we used the formulae given at the webpage [240].  We will call the first, the second and the third solution S 21 , S 22 and S 23 , respectively. With these conditions, the SM-like Higgs gauge and fermion couplings, as well as the trilinear Higgs self-coupling approach the corresponding SM values. In the following, we will call the first, the second and the third solution S 31 , S 32 and S 33 , respectively. For solution S 3 we obtain the SM values of the H 3 couplings to gauge bosons and fermions as well as for the trilinear Higgs self-coupling.