Showcasing HH production: Benchmarks for the (HL-)LHC

Current projections suggest that the LHC will have only limited sensitivity to di-Higgs production in the Standard Model (SM), possibly even after the completion of its high luminosity phase. Multi-Higgs final states play a fundamental role in many extensions of the SM as they are intrinsically sensitive to modifications of the Higgs sector. Therefore, any new observation in multi-Higgs final states could be linked to a range of beyond the SM (BSM) phenomena that are not sufficiently addressed by the SM. Extensions of the Higgs sector typically lead to new phenomenological signatures in multi-Higgs final states that are vastly different from the SM expectation. In this work, we provide a range of signature-driven benchmark points for resonant and non-resonant BSM di-Higgs production that motivate non-SM kinematic correlations and multi-fermion discovery channels. Relying on theoretically well-motivated assumptions, special attention is devoted to the particular case where the presence of new physics will dominantly manifest itself in multi-Higgs final states.


I. INTRODUCTION
The Higgs precision spectroscopy program that ensued after the discovery of the Higgs boson in 2012 [1,2] has assumed a central role in particle physics over the past years. One reason why the Higgs couplings and properties measurements move increasingly into the focus of searches for physics beyond the Standard Model (BSM) is the lack of conclusive hints for new interactions in the plethora of BSM searches performed by the ATLAS and CMS experiments. The Higgs boson as the direct implication of electroweak symmetry breaking is typically considered as a harbinger of new physics due to its special role in the unitarization of scattering amplitudes at high energy [3][4][5][6] and its relation to naturalness of the electroweak scale [7], to only name a couple of examples.
Although modifications of Higgs physics at the TeV scale of this size are still well within the limits set by recent 13 TeV LHC measurements, now that the SM can be considered as complete no additional ultra-violet (UV) energy scale can be predicted from the SM alone. This becomes even more pressing when coupling measurements stay consistent with the SM expectation in the future.
Additional requirements are more conveniently imposed in model-specific approaches, which try to mend apparent shortcomings of the SM such as the lack of a viable dark matter candidate or an insufficiently firstorder electroweak phase transition to address criteria of baryogenesis [8]. While model-independent approaches based on effective field theory (EFT) [9] can inform UV completions that address these questions through matching calculations, the appearance of novel phenomenolog-ical signatures such as resonances or thresholds within the LHC's kinematic coverage typically falls outside the region of reliability of these techniques.
A process that highlights the shortcomings of EFT methods in the presence of thresholds already in the SM context is di-Higgs production pp → hh (see e.g. [10,11] for recent summaries). Therefore, EFT methods for ggh n that can be summarized as [12][13][14][15][16] reflect the destructive interference between the top triangle and box diagrams of Fig. 1. While it can be assumed that the Higgs-top coupling can be accessed at the LHC in the range of ∼ 10% and the intermediate top mass scale is under theoretical control at the nextto-leading order (NLO) level * [17][18][19][20][21][22][23], the situation for the trilinear coupling is less clear as it probes a direction in the dimension-6 linear EFT space ∼ c 6 (H † H) 3 /Λ 2 (we denote the SM Higgs doublet with H in the following) [24]. As c 6 is essentially a free parameter unless a matching calculation is performed, its size is only limited by technical considerations related to perturbation theory, on which we need to rely to make phenomenological predictions. This has raised the question of how large trilinear coupling modifications can be [25][26][27][28] to inform di-Higgs investigations. Current CMS projections show that a sensitivity of 0.5 σ at 3 ab −1 [29] seems a realistic target at the LHC. Recent theoretical studies are slightly more optimistic [30][31][32], but the Higgs trilinear coupling will only understood in great detail even at large luminosity, in particular in the light of possible top quark Yukawa coupling modifications [33]. As a consequence, di-Higgs production is a main motivation for considering a major energy upgrade of the LHC or a future hadron collider [34][35][36][37][38][39].
In concrete UV scenarios that address fundamental BSM questions, typically a number of exotic states appear in the spectrum that will not only impact the "standard" single Higgs phenomenology at the LHC but can lead to even more dramatic changes in the production of multiple Higgs bosons. For instance, the triangle and box diagrams probe different aspects of top physics in extended top sectors [40][41][42][43][44].
Interleaving modifications of single Higgs physics with theoretically well-motivated UV considerations can therefore turn di-Higgs production into a strong probe of new physics: new kinematic features can appear that motivate new final states and search strategies, which might currently be unconsidered; di-Higgs production can be enhanced or suppressed. Such phenomenological modifications become particularly relevant when extrapolations of standard single Higgs channels do not show a significant departure from their SM expectation in these scenarios.
We will address these questions in this work using a particular set of models that allow us to contrast precise theoretical and phenomenological requirements with concrete predictions of single-and multi-Higgs production. Imposing, e.g., a strong first-order electroweak phase transition, dark matter constraints, electric dipole measurements and consistency with current Higgs coupling measurements as well as an extrapolation thereof, we discuss the results of a comprehensive scan of the models' parameter space with a particular emphasis on the relevance of multi-Higgs final states. We distill this scan into a number of representative benchmark points of BSM theories that highlight the importance of di-Higgs measurements in the future. In passing, we discuss how non-SM signatures are correlated with modifications of single Higgs physics BSM effects (or a lack of the latter). In turn, this also allows us to formulate an upper limit of the SM-like di-Higgs production cross section in these models when there are no conclusive hints for new physics in single Higgs phenomenology.
This work is structured as follows: In Sec. II, we outline the models that we consider for the purpose of this u-type d-type leptons type I (T1) Φ2  Φ2  Φ2  type II (T2) Φ2  Φ1  Φ1  lepton-specific Φ2  Φ2  Φ1  flipped Φ2 Φ1 Φ2 work. These are the CP-violating 2-Higgs-Doublet Model (C2HDM) and the Next-to-Minimal Supersymmetric Extension (NMSSM). Both models are special in the sense that they feature extended Higgs sectors that allow for Higgs-to-Higgs decays, also into final states with different Higgs bosons, or even cascade Higgs-to-Higgs decays. In Sec. III we outline the details of our scan over the model space and discuss the overall scan results. We present our benchmarks for the C2HDM and NMSSM together with their phenomenological properties in Sec. IV. We conclude in Sec. V.

II. THE MODELS
A. The C2HDM The 2-Higgs-Doublet Model (2HDM) [45][46][47] is obtained from the SM by adding a second SU (2) L Higgs doublet. The Higgs potential of a general 2HDM with a softly broken Z 2 symmetry, under which Φ 1 → Φ 1 and Φ 2 → −Φ 2 , can be written as The absence of Flavour Changing Neutral Currents (FCNC) is guaranteed by extending the Z 2 symmetry to the fermions. Depending on the Z 2 charge assignments, there are four phenomenologically different types of 2HDMs that are summarized in Tab. I. Hermiticity of the potential requires all parameters to be real, except for λ 5 and m 2 12 . If they have different unrelated complex phases we are in the framework of the complex or CPviolating 2HDM [48], which depends on ten real parameters. In the description of the C2HDM we follow the conventions of [49]. The vacuum expectation values (VEVs) developed by the Higgs doublets after electroweak symmetry breaking (EWSB) could in principle be complex in the C2HDM. Since the phase can be removed by a basis change [48], we set it to zero without loss of generality. In terms of the complex charged fields φ + i and the real B. The NMSSM Supersymmetric (SUSY) extensions require the introduction of at least a second Higgs doublet. In the NMSSM, the minimal field content with the doublet su-perfieldsĤ u andĤ d is extended by a complex superfield S (for reviews on the NMSSM, see for example [54,55]). The NMSSM Higgs potential is derived from the superpotential, the soft SUSY breaking Lagrangian and the D-term contributions. The scale-invariant NMSSM superpotential reads in terms of the hatted superfields For simplicity, we only included here the third generation fermion superfields, given by the left-handed doublet quark (Q 3 ), and lepton (L 3 ) superfields as well as righthanded singlet quark (t c R ,b c R ) and lepton (τ c R ) superfields. The soft SUSY breaking Lagrangian contains the mass terms m x for the Higgs (x = H u , H d , S) and sfermion (x =Q 3 ,t R ,b R ,L 3 ,τ R ) fields, obtained from the complex scalar components of the superfields. The Lagrangian with the trilinear soft SUSY breaking interactions A λ,κ,t,b,τ between the sfermion and Higgs fields reads The contribution to soft SUSY breaking from the gaugino mass parameters M 1,2,3 of the bino (B), winos (W ) and gluinos (G), respectively, is given by Expanding the tree-level scalar potential around the nonvanishing VEVs of the Higgs doublet and singlet fields, leads to the Higgs mass matrices for the three scalar (h d , h u , h s ), the three pseudoscalar (a d , a u , a s ) and the charged Higgs states (h ± u , h ∓ d ). We choose the VEVs v u , v d and v s to be real and positive. The three CPeven mass eigenstates H i (i = 1, 2, 3) are obtained from the interaction states through rotation with the orthogonal matrix R S that diagonalizes the 3 × 3 mass matrix squared of the CP-even fields, The mass eigenstates are ordered by ascending mass, M H1 ≤ M H2 ≤ M H3 . The CP-odd mass eigenstates A 1 and A 2 are obtained from a rotation R G that separates the Goldstone boson, followed by a rotation R P into the mass eigenstates, They are also ordered by ascending mass, M A1 ≤ M A2 . Altogether, the NMSSM Higgs spectrum consists of seven physical Higgs states, three neutral CP-even, two neutral CP-odd and two charged Higgs bosons. We use the three minimisation conditions of the scalar potential to replace the soft SUSY breaking masses squared for H u , H d and S in L mass by the remaining parameters of the tree-level potential so that the set of six parameters parametrising the tree-level NMSSM Higgs sector is given by The sign conventions are such that λ and tan β are positive, while κ, A λ , A κ and µ eff can take both signs. Note that the Higgs boson masses are not input parameters, but dependent parameters calculated from the input values. The inclusion of higher-order corrections in the Higgs boson masses is crucial here to shift the mass of the SM-like Higgs boson to the observed value of 125 GeV.

III. DETAILS OF THE SCAN
A. The C2HDM Scan The benchmark points ‡ provided in this paper have to satisfy theoretical and experimental constraints. In order to find valid points, we perform a scan in the C2HDM parameter space and additionally require the mass of one of the Higgs bosons, to be identified with the SM-like one and denoted by h, to be m h = 125.09 GeV [58]. The scan ranges are summarized in Tab. II. For simplicity, we only consider the C2HDM Type 1 (T1) and Type 2 (T2), which cover to a large extent the phenomenological effects to be expected in the C2HDM. Since physical parameter points with Re(m 2 12 ) < 0 are extremely rare, though possible, we neglect them in our scan. We test for ‡ 2HDM benchmarks for double Higgs production can be found in Refs. [56,57].  compatibility with the flavour constraints on R b [59,60] and B → X s γ [60][61][62][63][64] as 2σ exclusion bounds in the m H ± − tan β plane. In accordance with [64] we therefore require m H ± to be above 580 GeV in the C2HDM T2, whereas in the C2HDM T1 this bound is much weaker and depends more strongly on tan β. We verify agreement with the electroweak precision data by using the oblique parameters S, T and U -the 2HDM formulae are given in [47,65] -for which we demand 2σ compatibility with the SM fit [66], including the full correlation among the three parameters. We require perturbative unitarity to hold at tree level. The third neutral Higgs boson mass m Hj =Hi,h , which is not an independent input parameter in the C2HDM, but calculated from the other input values, is demanded to lie in the interval 10 GeV ≤ m Hj < 1.5 TeV.
To avoid degenerate Higgs signals, we additionally impose m Hi =h to be 5 GeV away from 125 GeV. For the SM input parameters we use [67,68] In order to perform the scans and find valid parameter points we use the program ScannerS [69,70]. Besides the above mentioned constraints it also tests for the potential to be bounded from below and uses the tree-level discriminant of [71] to enforce the electroweak vacuum to be the global minimum of the treelevel Higgs potential. Agreement with the Higgs exclusion limits from LEP, Tevatron and LHC is checked by using HiggsBounds5.2.0 [72][73][74] and with the Higgs rates by using HiggsSignals2.2.1 [75].
The required decay widths and branching ratios are obtained from the C2HDM implementation C2HDM HDECAY [53] in HDECAY [76,77]. In the production cross sections we included the QCD corrections taken over from the SM and  the Minimal Supersymmetric Extension (MSSM), where available. Electroweak corrections are consistently neglected both in production and decays. For more details on the production cross sections, see [53,78]. Working in the C2HDM, we also make sure to be in agreement with the measurements of the electric dipole moment (EDM), where the strongest constraint originates from the electron EDM [79]. We take the experimental limit given by the ACME collaboration [80]. Finally, we also investigate for the C2HDM if the parameters of the final data set induce a strong first order phase transition, a necessary condition for successful baryogenesis [8,81,82], by using the C++ code BSMPT [83].

B. The NMSSM Scan
In order to find benchmark points § that are compatible with the recent experimental constraints we proceed as described in [84][85][86], where also further details can be found. We perform a scan in the NMSSM parameter range with scan ranges summarized in Tab The soft SUSY breaking trilinear couplings of the first two generations are set equal to the corresponding values of the third generation. In order to ensure perturbativity we applied the rough constraint In accordance with the SUSY Les Houches Accord (SLHA) format [87,88] the soft SUSY breaking masses and trilinear couplings are understood as DR parameters at the scale µ R = M s = mQ 3 mt R . § For NMSSM benchmarks for double Higgs production from Higgsto-Higgs decays, see [84].
The SM input parameters have been chosen as in the C2HDM scan, with the exception of the top quark mass which has been set to m t = 173.5 GeV. The small difference of 1 GeV has no effect on the scan results.
The spectrum of the Higgs and SUSY particles including higher-order corrections is calculated with NMSSMTools5.2.0 [89][90][91][92][93][94] which also checks for the constraints from flavour and low-energy observables. It provides the input for HiggsBounds5.2.0 [72][73][74] to check for compatibility with the exclusion bounds from the Higgs searches. The mass of one of the neutral CP-even Higgs bosons, identified with the SM-like Higgs boson denoted by h, has to lie in the range and the masses of all other Higgs bosons are demanded to be separated by at least 1 GeV in order to avoid two overlapping signals. The signal strengths of this Higgs boson have to be in agreement with the signal strength fit of [95] at the 2 × 1σ level. For the computation of the signal strengths we need the production cross section and branching ratios for the NMSSM Higgs bosons. To compute production through gluon fusion and bb annihilation, we take the SM cross sections and multiply them with the effective couplings obtained from NMSSMTools. The SM values are calculated with SusHi [96,97] and include in gluon fusion the NLO corrections with the full top quark mass dependence [98] as well as the next-to-next-to-leading order (NNLO) corrections in the heavy quark effective theory [99][100][101][102][103]. The next-to-next-to-next-to-leading order (N 3 LO) corrections are taken into account in a threshold expansion [104][105][106][107] for Higgs masses below 300 GeV. For masses above 50 GeV, bb annihilation cross sections that match between the five-and four-flavor scheme are used, obtained in the soft-collinear effective theory [108,109]. They are in accordance with the results from [110,111]. For masses below 50 GeV, cross sections obtained in the Santander matching [112] are used, with the fiveflavor scheme cross sections from [113] and the fourflavor scheme ones from [114][115][116]. The branching ratios are taken from NMSSMTools and cross-checked against against NMSSMCALC [117].
The parameter points also have to satisfy the bounds from SUSY searches at the LHC. The gluino mass and the lightest squark mass of the second generation are demanded to lie above 1.85 TeV, respectively, see [118]. The stop and sbottom masses are required to be above 800 GeV, respectively, [118,119], the slepton masses above 400 GeV [118] and the absolute value of the lightest chargino mass above 300 GeV [120].
Through an interface with micrOMEGAS [94] we obtain the relic density that we require not to exceed the value measured by the PLANCK collaboration [121]. The spinindependent nucleon-dark matter direct detection cross section, also provided by micrOMEGAS, is required not to violate the upper bound from the LUX experiment [122]. We furthermore test for compatibility with the direct detection limits from XENON1T [123] and check the Dark Matter annihilation cross section against the results provided by FermiLat [124].

C. Extrapolations
We include a range of extrapolations of single Higgs measurements to our discussion to identify an approximate "exclusion luminosity" (see below) at which single Higgs measurements will start to become sensitive to a particular scenario and spectrum. This notion will allow us to put multi-Higgs final states in direct comparison with single Higgs measurement expectations and identify interesting regions of parameter space.
In particular, we include projections for the m h 125 GeV standard single Higgs production modes, gg → h (gluon fusion), qq → hjj (weak boson fusion), qq → V h, V = W ± , Z (Higgs radiation), pp → tth (associated production with top quarks), and consider decays h → ZZ, h → W W , h → γγ, h → bb and h → τ + τ − . For the decays h → γγ, ZZ we use the CMS projections provided in Ref. [29]; these include the production modes gluon fusion, weak boson fusion and tth. We interpolate between different luminosities using a √ L luminosity dependence at all times.
Projections for h → W W are obtained using Ref. [125] and we rescale these results taking into account the cross section differences between 13 TeV and 8 TeV using the results provided by the Higgs cross section working group [126]. For h → bb we consider extrapolations based on V h production [127], tth production [128] as well as weak boson fusion [129]. h → τ τ is based on [130], which agrees with the ECFA results of [29] upon projecting to 3 ab −1 .
The improved determination of the SM-like Higgs boson with m h 125 GeV needs to be contrasted with additional coverage of Higgs-like searches for masses m h = 125 GeV. We include projections of existing resonance searches in γγ [131], τ τ [132], W W [133] and ZZ [134] final states. By far the most constraining exotic searches result from tt resonance searches, and we extrapolate the results of Ref. [135]. This analysis is performed in the context of a Z model and can therefore be interpreted as only a rough estimate of the search potential of tt Higgs resonances. To our knowledge no comprehensive analysis of exotic heavy Higgs masses is publicly available. This is mostly due to the dedicated interference between the background and the Higgs signal that also depends on the CP character of the produced state [136]. This has a significant impact on the sensitivity of tt final states. Including such effects is beyond the scope of this work.

D. Results
Both in the C2HDM and the NMSSM the enlarged Higgs sector leads to a plethora of di-Higgs production processes. In particular, they feature processes with two different Higgs bosons in the final state. Compared to the SM, the cross sections can be largely enhanced in case of resonant production of a heavy Higgs boson that subsequently decays into a pair of lighter Higgs bosons, provided the Higgs self-coupling is not too small. Moreover, the different Higgs self-couplings themselves can enhance the cross section in view of the well-known fact that in the SM the triangle and box diagrams interfere destructively. Additionally, loops with bottom quarks may play a role in scenarios with enhanced down-type Yukawa couplings for large tan β in the NMSSM or the C2HDM T2. In the NMSSM loops with stop and sbottom quarks also contribute to Higgs pair production, and we furthermore have the possibility to produce a di-Higgs final state with pseudoscalars. These processes can yield even larger rates as has been discussed in detail in [85]. However, due to supersymmetry, the Higgs self-couplings are given in terms of the gauge couplings limiting to some extent deviations from the SM. This is not the case for the C2HDM, so that effects different from the NMSSM may be expected here. In the NMSSM, on the other hand, we have the possibility of the final state Higgs bosons to decay into non-SM final states like e.g. neutralinos, inducing signatures with phenomenologically interesting features. Altogether, both models provide a large playground for possible BSM effects in Higgs pair production that can be rather different.
Furthermore, we restrict ourselves to SM final states. Most of our results show the leading order (LO) Higgs pair production cross sections. For the benchmark points we also computed the NLO QCD corrections in the limit of heavy loop particles. They typically increase the cross section by about a factor two. We have implemented the NLO QCD corrections both for the NMSSM ¶ [85] and C2HDM [138] in the Fortran code HPAIR that was originally designed to compute the SM and MSSM Higgs pair production at NLO QCD. All Higgs pair production processes have been computed at a c.m. energy of √ s = 14 TeV, and we have adopted the CT14 parton densities [139] for the LO and NLO cross sections with ¶ They only include the corrections to the top quark loops. For NLO QCD corrections including also the squarks in the limit of vanishing external momenta, see [137]. See M. Spira's website, http://tiger.web.psi.ch/proglist.html. α s (M Z ) = 0.118 at LO and NLO. The renormalisation scale has been set equal to M HH /2, where M HH generically denotes the invariant mass of the final state Higgs pair. Consistent with the application of the heavy top quark limit in the NLO QCD corrections, we neglect the bottom quark loops in the LO cross section.
In view of the possibility of (largely) enhanced production of a pair of SM-like Higgs bosons, in the selection of valid scenarios, we also took into account limits set by LHC 4b [140][141][142], (2b)(2τ ) [143,144,160] and (2b)(2γ) [145] final states from the production of a heavy scalar resonance that decays into two 125 GeV Higgs bosons.

C2HDM
We start by discussing the possible sizes of Higgs pair production that are compatible with all present experimental constraints. Table IV summarizes the maximum cross section values found in the sample of valid parameter points where we additionally applied the extrapolations of Subsection III C. Taking these into account, we only kept the points that are not excluded at 64 fb −1 , which corresponds approximately to the present luminosity acquired by the LHC experiments. We will come back to the role of the extrapolations below. In the following, we denote the SM-like Higgs boson with a mass of 125 GeV h, the lighter of the non-SM like neutral Higgs bosons is called H ↓ , the heavier one H ↑ . All cross sections are calculated at LO QCD and hence still increase by approximately a factor two when QCD corrections are included.
The SM Higgs pair production cross section in gluon fusion amounts to 19.72 fb at LO and 38.19 fb at NLO * * with a K-factor of K = 1.95. Inspection of the Tab. IV shows, that both in the T1 and T2 scenarios, the maximum attained cross section for hh production can exceed the SM value, in T1 by a factor of about 40 HiHj), in the C2HDM T1 and T2, with an exclusion luminosity ≥ 64 fb −1 that satisfy all theoretical and experimental constraints described above. * * This value differs from the one given in [18]. This is because we do not include top quark mass effects here and use a different pdf set.
by about a factor 1.7. This is also the case at NLO, with NLO cross sections and K-factors for hh production in the T1 and T2 models given by The reason for the large enhancement in T1 is the resonant production of the heavier Higgs bosons H ↓ and H ↑ with a mass of 285 GeV and 287 GeV, respectively, that subsequently decay into a pair of SM-like Higgs bosons. This is also the reason for the enhancement in T2, where the masses of the non-SM-like Higgs bosons amount to 794 and 798 GeV. The reason for the much smaller enhancement in hh production in T2 compared to T1 is the overall heavier Higgs spectrum in T2. In particular, the intermediate heavy resonances in the T2 scenarios that can produce hh in their decay usually fall into the heavy mass range where the ATLAS and CMS limits on the upper cross section for (4b) production drop rapidly, cf. [140,141]. In T1, furthermore, the maximum values of di-Higgs production processes involving H ↓ can compete with SM Higgs pair production or even largely exceed it. Thus the production of a SM-like Higgs boson and H ↓ can be larger by a factor 2.5. This final state is interesting as it is clearly a non-SM-like signature where the experiments can use the SM-like Higgs boson to calibrate or "tag" this exotic configuration. This does not apply for T2, however, where due to the experimental constraints, the non-SM-like Higgs bosons are in general heavier than in T1, inducing small di-Higgs production processes due to a much smaller phase space.

Experimental accessibility and exclusion luminosity
In order to assess the experimental accessibility of these cross sections, however, we need to look at their decay products. We therefore applied the narrow width approximation and multiplied the produced Higgs bosons with their branching ratios in various SM final states. The most promising final states for the investigation of Higgs pair production at the LHC are the (bb)(γγ) [146], the (bb)(τ τ ) [147][148][149] and (bb)(bb) [147,150,151] final states (for other final states see also [152][153][154]). In Fig. 2 we show for all parameter points that pass our applied constraints for the C2HDM T1, the cross section values of SM Higgs pair production in the (bb)(γγ) final state (left) and for H ↓ H ↓ production in the (bb)(bb) final state (right) normalized to the corresponding SM values as a function of the exclusion luminosity. By the latter we define the luminosity at which this process would be excluded experimentally, based on the extrapolations described in Subsection III C. First of all we notice that the figures contain parameter points at lower luminosity that should have been excluded by HiggsBounds. The reason why they are there is that HiggsBounds relies on the published experimental results and cannot check for certain signatures that become relevant in BSM Higgs sectors. Thus there exist Higgs spectra with heavy Higgs bosons that dominantly decay into top quark pairs. These would induce exotic four-top final states in heavy Higgs pair production. Such signatures compete, however, with single heavy Higgs production and subsequent decay into a top-quark pair. Applying our rough estimate on the exclusion power of the experiments for this process, based on the Z data, such scenarios are excluded already, although they have been let through by HiggsBounds due to the lack of a dedicated experimental analysis for this. This shows the importance of experimental analyses investigating top pair final states from heavy Higgs production in order to properly assess the exclusion limits for BSM Higgs sectors -with dramatic effects on possible Higgs pair production signatures. While our rough extrapolation excludes about 0.6% of the T1 points for a luminosity of about 36 fb −1 , the effect is much larger for the T2 sample allowed by HiggsBounds † † . Here about 22% of the points would be excluded. This is because of the overall heavy non-SM-like Higgs bosons in T2 and their prominent decays into top-quark pairs.
As can be inferred from the figures in the C2HDM T1, the production of a SM-like Higgs pair with subsequent decay into (bb)(γγ) can exceed the SM rates by up to a factor 60. This maximum enhancement factor is the same for all final states, as the branching ratios of the SM-like Higgs boson h are almost the same as in the SM. In the following, we will use the quantity to classify whether a Higgs boson X has a sizable non-SM branching ratio and decay phenomenology. If Σ X 1 then the exotic states can be dominantly discovered in "standard" SM-Higgs-like decay channels, e.g. X → bb or tt if the mass of X permits such a decay.
In the H ↓ H ↓ final state with both H ↓ 's decaying into bottom quarks the enhancement can even be up to a factor of about 200. The point with the maximum enhancement corresponds to the one quoted in Tab. IV and the enhancement is due to the large di-Higgs production process of 3.2 pb and a slightly enhanced branching ratio into b-quarks as compared to the SM. The same factor is found for the (bb)(ττ ) final state. Due to a smaller branching ratio into photons, however, the maximum allowed enhancement in the (bb)(γγ) final state only amounts up to a factor of 40. The H ↓ in this scenario has a mass of m H ↓ = 131 GeV, and the mass of H ↑ is m H ↑ = 313 GeV. Its main branching ratios are BR(H ↑ → ZH ↓ )= 0.53 and BR(H ↑ → H ↓ H ↓ ) = 0.46. The maximum branching ratios of the charged Higgs boson with a mass of m H + = 312 GeV are BR(H + → W + H ↓ )=0.65 and BR(H + → tb)=0.34. With its large di-Higgs production cross section and the large non-SM-like branching ratios, this parameter point is an interesting scenario for studying new physics effects (also beyond the Higgs pair events that we consider here).
All remaining di-Higgs production processes are less promising. Thus the enhancement factor for hH ↓ production remains below 3 in the 4b and 2b2τ final state and below 2 in the 2b2γ final state. All other final states range below the SM values.
As can already be inferred from the maximum di-Higgs production values in T2, given in Tab. IV the situation looks much less promising in the C2HDM T2. There are very few points in hh production with subsequent decay into the (2b)(2τ ) and 4b final state that exceed the SM rate, and only by a factor of about 2.4. The maximum enhancement found in the (2b)(2γ) final state is about 2.4. All other final states lead to smaller rates than in the SM.
From these considerations we can conclude that there are promising di-Higgs signatures with large rates in the C2HDM T1 both for SM-like Higgs pair production but also for final states with non-SM-like Higgs bosons. The exotic Higgs bosons appear in SM-like final states, however, with different kinematic correlations due to different masses. This highlights the need to conduct Higgs pair analyses in a broad range of kinematic possibilities. Furthermore, the strict constraints on T2 scenarios, would exclude the model if di-Higgs signatures much larger than in the SM are found.

NMSSM
In Table V we summarize for the NMSSM the maximum di-Higgs production cross section values found in the sample of valid parameter points that are not excluded at a luminosity of 64 fb −1 . All cross sections are calculated at LO QCD and hence still increase by approximately a factor two when QCD corrections are included. By A ↓ we denote the lighter of the two pseudoscalar Higgs bosons.  The reason for the large enhancement of σ(gg → hh) is the intermediate resonant production of heavy Higgs bosons H ↓ and H ↑ with subsequent decay into a SMlike Higgs pair. The H ↓ H ↓ production cross section is so large because of the smallness of the H ↓ mass, m H ↓ = 39.52 GeV. The enhancement in hA ↓ production is due to the resonant A 2 ≡ A ↑ production with subsequent decay into hA ↓ . The huge enhancement in A ↓ A ↓ production is on the one hand due to the smallness of the A ↓ mass of m A ↓ = 37 GeV, on the other hand due to the resonant H ↑ production with subsequent decay into A ↓ A ↓ (the resonant H ↓ production plays a minor role). Searches for relatively low mass states are performed in the γγ [155] and τ τ [156] channels, however, with rather limited sensitivity.

Experimental accessibility and exclusion luminosity
In Fig. 3 we show for all parameter points that pass our applied constraints, the NMSSM cross section values of SM Higgs pair production in the (bb)(bb) final state (left) and for A ↓ A ↓ production in the (bb)(bb) final state (right) normalized to the corresponding SM values as a function of the exclusion luminosity.
As can be inferred from Fig. 3 (left), the 4b final state rates from SM-like Higgs pair production exceed the SM reference value by less than a factor 10 and only for lower exclusion luminosities. Large enhancement factors are basically limited by the LHC upper limits on heavy resonant scalar production with subsequent decay into a SM-like Higgs pair. The situation looks even less promising in the production of an SM-like Higgs boson together with the lighter of the CP-even non-SM-like Higgs bosons, where only an enhancement factor slightly above 2.3 at most is found. This is the case for high exclusion luminosities beyond 1 ab −1 so that nevertheless this process might be accessible at high luminosities. The situation is different in the production of hA ↓ . Because the lighter pseudoscalar can be relatively light and decays dominantly into (bb) ‡ ‡ we can have enhancement factors above 10 up to about 45 in the 4b state. This makes it particularly interesting, moreover in view of the exotic final state with two different Higgs masses in di-Higgs production. The enhancement factors can become huge in A ↓ A ↓ production, which is mainly due to the lightness of A ↓ . In 4b production it can be up to 1000. For larger exclusion luminosities the enhancement factor can still be a factor up to 10, cf. Fig. 3 (right). In the (bb)(γγ) final state the enhancement can be larger than 100 up to about 360. Figure 4 nicely summarizes the situation of the enhanced di-Higgs cross sections involving very light Higgs bosons. It shows the production of A ↓ A ↓ with subsequent decay in the 4b final state normalized to the value of the corresponding process for the SM-like di-Higgs production, as a function of the mass of the light pseudoscalar. The color code denotes the exclusion luminosity. For very light masses below 125 GeV (note the gap around 125 GeV is due to our scan procedure) the rates are largely enhanced because of the large di-Higgs production cross sections. With increasing mass the rates decrease. The exclusion luminosities are high for exotic Higgs masses above 125 GeV and below the top-pair threshold. Above the top-pair threshold the exclusion luminosities are much lower due to the exclusion limits in the top-pair final state. For masses below the SM-like ‡ ‡ Note that typical trigger criteria are too selective to directly observe pp → A ↓ → bb. Higgs mass, however, there are parameter points where the exclusion luminosities can exceed the 100 fb −1 and even 1 ab −1 while still featuring large rates. The reason is that these points are not excluded from single Higgs searches as light Higgs states with dominant decays into bb final states are difficult to probe. On the other hand this enhancement combined with the large di-Higgs production cross section implies huge 4b final state rates, that may be tested at the high luminosities, but with associated experimental difficulties. This is a nice example of the interplay between difficult single-Higgs searches and large exotic di-Higgs rates, where new physics may be found.

A. Type 1 Benchmarks
We describe a representative set of benchmarks of the C2HDM T1 model and their associated (exotic) multi-Higgs phenomenology. The input parameters, the derived third neutral Higgs boson mass, the CP-odd admixtures in terms of the squared mixing matrix elements R 2 i3 and the exclusion luminosity L excl are given in Tab. VI. We also give the NLO QCD gluon fusion hh production cross section at √ s = 14 TeV together with its K-factor, given by the ratio of the NLO cross section to the LO one. In Table VII we present the 4b, (2b)(2τ ) and (2b)(2γ) rates from Higgs pair production normalized to the rate expected in the SM from Higgs pairs relevant for the discussion of the various benchmark points.
• T1BP1 -Highest exclusion lumi: The point with the highest exclusion luminosity in the complete sample.
The exclusion luminosity for this point is found to be 11.5 ab −1 , i.e. well above the LHC design luminosity even after the high luminosity phase. All di-Higgs cross sections involving non-SM-like Higgs bosons have values below the SM reference point. Altogether this benchmark point behaves very SMlike as expected for such a high exclusion luminosity. The neutral Higgs mass spectrum is relatively degenerate with all masses in the vicinity of the SM-like Higgs boson at approximately 125 GeV. The SM-like Higgs pair production cross section for this point is SM-like; bbγγ is about 10% below the SM expectation while bbτ τ and bbbb are slightly enhanced by 5%, cf. Tab. VII. The cases when a SM Higgs is accompanied by an additional     exotic Higgs are around 30% of the SM expectation in bbτ τ . With the latest improvements in hadronic tau tagging [157][158][159] such a signature might be phenomenologically accessible at the LHC with large luminosity.
• T1BP2 -HighLMaxHsmHl: Among the points with exclusion luminosities ≥ 1 ab −1 the point with the maximum cross section gg → hH ↓ .
We have σ(gg → hH ↓ ) = 28.47 fb at LO. As summarized in Tab. VII, for the final states involving b's, τ 's and γ's we find for this di-Higgs final state relative to the production of a SM Higgs pair with subsequent decay into the same final state: The neutral Higgs mass spectrum in this case is slightly split while the pair production of the SMlike Higgs bosons follows largely the SM paradigm. The mass of the lighter neutral state is rather close to the SM boson, which allows us to compare the rates with the SM itself. H ↓ has a decay phenomenology that is SM-like. At 3 ab −1 we can therefore expect around 26k exotic bbτ τ events while the bbbb rate is enhanced by 70% over the SM expectation.
This point will have a highly non-SM decay phenomenology and cascade decays are an interesting avenue to look for such a scenario. All di-Higgs production cross sections involving non-SM-like Higgs bosons lie below 5 fb. While the point has interesting signatures for non-SM-like single Higgs production, exotic di-Higgs production is not a very promising avenue. The cross sections are far below the SM value, as no resonance-enhancement is possible as the Higgs mass values are too close. Additionally the branching ratios in SM-like Higgs states are very small, as decays into non-SM-like final states dominate. We have Σ H ↓ = 0.180, Σ H ↑ = 0.149. The SM-like Higgs pair production is consistent with the SM expectation within 10% for this point. Exotic production of H ↓ together with the SM Higgs leads to around 220 exotic (bb)(τ τ ) events and around 1900 (bb)(bb) at 3 ab −1 .
• T1BP4 -MaxLEnhancedHsmHl : Among the points with gg → hH ↓ ≥ 10 fb the point with the maximum exclusion luminosity. The SM-like Higgs pair production modes fall again within ∼ 10% of the SM expectation. The mass of the lightest neutral Higgs boson of about 120 GeV is again reasonably close to the SM-like Higgs to allow a direct comparison of expected rates, which are slightly smaller than the SM. Specifically, the exotic bbτ τ and bbbb modes are 35% smaller than what we would expect for the SM mode with m h 125 GeV.
• T1BP5 -EnhancedHsmHsm: Among the points with exclusion luminosities ≥ 1 ab −1 the point with the maximum cross section gg → hh.
• T1BP6 -EnhancedHlHl: Among the points with exclusion luminosities ≥ 1 ab −1 the point with the maximum cross section gg → H ↓ H ↓ .
The exclusion luminosity is 2.58 ab −1 and we have a light H ↓ with a mass just above half the SMlike mass so that the branching ratios of the latter remain in accordance with the LHC data. The di-Higgs production of H ↓ H ↓ amounts to 1.249 pb at LO QCD. Comparing to the SM expectation, we have σ(gg → H ↓ H ↓ → 4b)/SM = 145 , σ(gg → H ↓ H ↓ → (2b)(2τ ))/SM = 124 , σ(gg → H ↓ H ↓ → (2b)(2γ))/SM = 0.29 , due to the suppressed decay H ↓ → γγ. The light Higgs almost exclusively decays into b pairs at a branching fraction of 85%. This means that such a state is difficult to observe in single Higgs production as trigger criteria typically remove such events from the busy hadronic LHC environment. There is a possibility to observe this state in its τ modes ( 8%). This point is a nice example how the 4b mode can be an important BSM discovery tool when SM mass correlations are relaxed.
We note that none of these points features a strong first order electroweak phase transition. In general we observe that points with a strong first order phase transition do not lead to enhanced rates in the four-particle final states; the cross section for SM-like di-Higgs production is close to the SM value while other di-Higgs production cross sections are smaller than the SM expectation.
In addition to these points which are all characterized by relatively light exotic states we include benchmarks with heavy neutral exotics. This is achieved by adding the additional requirement m H ↑ ≥ 400 GeV to the defining criteria of the benchmark points quoted above. We restrict ourselves to three benchmark points highlighting special features. They are called T1BP1 H, T1BP2 H and T1BP5 H in analogy to the their lighter mass spectrum counterparts T1BP1,2,5. Detailed information on the benchmark points is summarized in Tab. VIII.
• T1BP1 H -Highest exclusion lumi: The exclusion luminosity for this point is found to be 2.46 ab −1 . All di-Higgs cross sections involving non-SM-like Higgs bosons have values below the SM reference as the Higgs pair production cross section steeply falls for heavy Higgs production. The additional resonant structures that are sourced in the SM-like Higgs pair production amount to an increase above the SM expectation by a factor of ∼ 5.8 across the standard search channels 4b, (2b)(2τ ), (2b)(2γ).
• T1BP2 H -HighLMaxHsmHl : We have σ(gg → hH ↓ ) = 2.34 fb at LO, which is rather large given the mass of the exotic Higgs. The decay phenomenology of the additional Higgs is completely dominated by decays into top-final states. In this sense the single Higgs production and exotic hH ↓ production are fully correlated. The exclusion luminosity is 1.8 ab −1 and results from the extrapolation of the tt resonance search. This point, although not relevant for di-Higgs analyses shows how single Higgs measurements in the tt channel influence multi-Higgs final states. Such a benchmark could be adopted to further clarify the role of single Higgs measurements for exotic multi-Higgs final states.
This point hence gives access to SM-like Higgs pair production even for a heavy Higgs spectrum, and is an immediate sign of BSM physics as the di-Higgs cross section is enhanced.

B. Type 2 Benchmarks
As already visible from Tab. IV, the C2HDM T2 model gives rise to less spectacular signatures than the C2HDM T1. We give two representative scenarios below. In general, the spectrum is much heavier than for T1. There is no scenario with the SM-like Higgs boson being H 2 ; the SM Higgs is always the lightest state H 1 . The input parameters for these points as well as further relevant information are summarized in Tab. IX.
• T2BP1 -Highest exclusion lumi: The exclusion luminosity for this point is found to be 2.66 ab −1 . All di-Higgs cross sections involving non-SM-like Higgs bosons have values below the SM reference value. Altogether this benchmark point behaves very SM-like as expected for such a high exclusion luminosity where tt resonance searches become sensitive to this scenario.
Despite a lower exclusion luminosity this point behaves also very SM like and barely exceeds the rates into SM-like final states of T2BP1 which has a higher exclusion luminosity.
Overall, it will be difficult to probe the C2HDM T2, which features a heavy Higgs spectrum, in di-Higgs production. This is partly also due to the fact that enhanced SM-like Higgs pair production cross sections are already excluded by the LHC limits on resonant heavy scalar production with decay in a SM Higgs boson pair.

C. NMSSM Benchmarks
Let us finally turn to the NMSSM. The criteria for selecting the benchmark points are as follows.
• NMBP1: The point with the largest 4b rate from SMlike Higgs boson pair production with an exclusion luminosity above 1 ab −1 .
• NMBP2: The same as BP1 but with an exclusion luminosity beyond 100 fb −1 .
• NMBP3: The point with the largest 4b rate from the production of an SM-like Higgs boson and the lighter of the CP-even non-SM-like Higgs bosons, H ↓ , with an exclusion luminosity above 1 ab −1 .
• NMBP4: The point with the largest 4b rate from H ↓ H ↓ production with an exclusion luminosity above 100 fb −1 .
We also provide benchmark points for di-Higgs final states involving a light pseudoscalar A ↓ in the final state: • NMBP5: The point with the largest 4b rate from hA ↓ production with an exclusion luminosity above 1 ab −1 .   • NMBP6: The same as BP5 but with an exclusion luminosity beyond 100 fb −1 .
• NMBP7: The point with the largest 4b rate from A ↓ A ↓ production and an exclusion luminosity  above 1 ab −1 . It turns out that NMBP7 is identical to NMBP5.
• NMBP8: The same definition as for NMBP7 but with an exclusion luminosity beyond 100 fb −1 .
The input values, the derived Higgs boson masses, the exclusion luminosity and the NLO QCD cross section for hh production with its K-factor for the various benchmark points are listed in Tabs. X and XII. Note, that we include benchmark points with an exclusion luminosity around 100 fb −1 in case the rates are much enhanced compared to the SM as in this case a luminosity of 100 fb −1 might be enough to test this parameter point. As the rates for di-Higgs production involving a heavy scalar or pseudoscalar are low we do not present benchmarks for these cases.
From Tab. XI we can read off that the 4b rates from hh production for an exclusion luminosity above 100 fb −1 (NMBP2) can be up almost a factor of about 4.8, so that this process will be difficult to be accessed at the lower luminosity. Assuming an exclusion luminosity above 1 ab −1 (NMBP1) the enhancement compared to the SM rate is around 2.4. Again, larger enhancements in the final state with an SM-like Higgs boson pair are excluded by the limits provided from ATLAS and CMS [141,142] . In the hH ↓ final state the enhancement factor is only 2.3 as the di-Higgs production cross section NMBP3 : σ(hH ↓ ) = 49.13 fb (K = 1.92) , is not much larger than in the SM. For H ↓ H ↓ production we have, however, leading to an enhancement factor of up to 540. The H ↓ mass is much lower here than in the former case inducing dominant branching ratios into bb. The large cross section is mainly due to the small H ↓ mass. The resonant enhancement plays a minor role here. Both benchmark points are special in the sense that the SM-like Higgs boson is not the lightest but the second lightest CP-even Higgs boson in the spectrum. We hence have a light CP-even Higgs boson in these scenarios. The H ↓ is very singlet-like in both cases and decays with a branching ratio of about 0.9 into bb. From Tab. XIII we can read off that the 4b, (2b)(2τ ) and (2b)(2γ) final states from both hA ↓ and A ↓ A ↓ production can be enhanced above the SM rate even for an exclusion luminosity of 1 ab −1 . For the lower exclusion luminosity the enhancement can be huge, in particular in the 4 fermion final state from A ↓ A ↓ production (NMBP8). These enhancements are due to large di-Higgs production cross sections which at NLO QCD amount to NMBP5 : σ(hA ↓ ) = 56.04 fb (K = 1.93) , NMBP6 : σ(hA ↓ ) = 988 fb (K = 1.99) , NMBP8 : σ(A ↓ A ↓ ) = 34.59 pb (K = 2.30) .
The enhanced di-Higgs cross section values are on the one hand due to the light pseudoscalar masses and on the other hand due to resonant scalar production for A ↓ A ↓ or pseudoscalar production for hA ↓ production. As already noted in the discussion of Fig. 4 this is an example where new physics may lead to huge measurable effects in Higgs pair production while the single Higgs process, here A ↓ production in gluon fusion, is difficult to access. This is a prime example that demonstrates that despite the very SM-like nature of the 125 GeV Higgs boson Higgs pair production can be far from being SM-like.
We finally remark that in all NMSSM scenarios the stop masses are quite large, of the order of 1 TeV and larger.

V. CONCLUSIONS
Multi-Higgs final states are statistically limited at the LHC, but are key processes to gain a precise understanding of the mechanism of electroweak symmetry breaking. Phenomenologically, they are highly correlated with measurements in single Higgs final states. The question of how much additional information can be gained from the investigation of multi-Higgs final states is therefore best addressed using concrete BSM extensions.
Another particularly relevant question when considering di-Higgs final states is whether they could play a key discovery tool for BSM interactions. This could happen at the LHC in situations when single Higgs analyses are simply not competitive due to, e.g., trigger thresholds that might be mitigated in more complex multi-Higgs final states. Large branching ratios of additional scalars into SM-like Higgs bosons serve as an additional avenue to observe resonantly enhanced SM-like Higgs production. In these scenarios, the kinematic correlations are often significantly modified compared to the SM.   In this work we have performed a comprehensive scan over the complex 2HDM and the NMSSM, with a particular emphasis on the expected di-Higgs phenomenology in these models, taking into account a variety of current constraints and future projections. We find that in particular in the C2HDM Type 1 models, the di-Higgs phenomenology can significantly differ from the SM expectation. The differences range from new signatures in SM-like search channels for light Higgs bosons all the way to new resonant structures in di-Higgs final states for larger exotic Higgs masses. In particular, in final states involving light Higgs bosons the multi-fermion final states can be significantly enhanced compared to the SM case. This is also the case in the NMSSM where we can have light scalars or pseudoscalars in the spectrum. In the C2HDM Type 2 models with its heavy Higgs spectrum, the multi-Higgs final states play a more subdominant role as new physics discoveries typically occur in tt resonance searches before (exotic) di-Higgs production becomes relevant.
We have distilled our scans into a representative number of benchmark points that not only reflect the phenomenological possibilities that present themselves in non-standard Higgs sectors, but also point to a particular range of phenomenological situations: Firstly, the expected sensitivity of the tt resonance search is a crucial factor in deciding on the relevance of di-Higgs searches. The multi-Higgs signal is typically driven by top-mediated gluon fusion. Therefore, the decay to top final states is directly correlated with a large di-Higgs cross section for resonant production as well as enhancements in the decay to photons for non-resonant production.
Secondly, when exotic Higgs masses fall below the tt threshold, di-Higgs final states typically follow the SM decay rates with compressed neutral Higgs masses. This highlights the necessity to achieve a high mass resolution in the standard search channels (2b)(2b), (2b)(2τ ), (2b)(2γ), even when SM-mass correlations are abandoned. While di-Higgs production would be enhanced in this instance, providing clear evidence of the presence of BSM interactions, their precise nature would remain elusive to some extent.
Alternatively, additional Higgs exotics can create multiple resonant features leading to a large enhancement of the total SM Higgs pair production rate. The extrapolated signal-strength constraints locate viable candidates for enhanced di-Higgs production in bb final states, which are difficult to access experimentally in single Higgs production. Here di-Higgs production can play a significant role as a discovery tool for BSM interactions due to smaller backgrounds and better kinematical handles.
Thirdly, relatively light Higgs bosons with significant branching ratios can lead to a strong enhancement in multi-fermion final states. Such signatures are already studied by the experimental collaborations. Our results indicate the importance of these analyses in the future.
While we have specifically focused on di-Higgs production, it is clear that these scenarios can have interesting non-SM signatures that can be exploited to observe or constrain a certain parameter point in a more targeted, yet parameter point-dependent way. We will leave this for future work.