Singlet Extensions of the Standard Model at LHC Run 2: Benchmarks and Comparison with the NMSSM

The Complex singlet extension of the Standard Model (CxSM) is the simplest extension that provides scenarios for Higgs pair production with different masses. The model has two interesting phases: the dark matter phase, with a Standard Model-like Higgs boson, a new scalar and a dark matter candidate; and the broken phase, with all three neutral scalars mixing. In the latter phase Higgs decays into a pair of two different Higgs bosons are possible. In this study we analyse Higgs-to-Higgs decays in the framework of singlet extensions of the Standard Model (SM), with focus on the CxSM. After demonstrating that scenarios with large rates for such chain decays are possible we perform a comparison between the NMSSM and the CxSM. We find that, based on Higgs-to-Higgs decays, the only possibility to distinguish the two models at the LHC run 2 is through final states with two different scalars. This conclusion builds a strong case for searches for final states with two different scalars at the LHC run 2. Finally, we propose a set of benchmark points for the real and complex singlet extensions to be tested at the LHC run 2. They have been chosen such that the discovery prospects of the involved scalars are maximised and they fulfil the dark matter constraints. Furthermore, for some of the points the theory is stable up to high energy scales. For the computation of the decay widths and branching ratios we developed the Fortran code, sHDECAY, which is based on the implementation of the real and complex singlet extensions of the SM in HDECAY.


Introduction
The ATLAS [1] and CMS [2] collaborations at the Large Hadron Collider (LHC) have established the existence of a new scalar particle identified with the Higgs boson. The investigations are now focused on pinning down the pattern of electroweak symmetry breaking (EWSB), which already started with the measurements of the Higgs couplings and will continue with the search for physical processes that could give us some direct insight on the structure of the potential. These are the processes involving triple scalar vertices [3][4][5] from which the resonant ones are the most promising to lead to a detection during the next LHC runs. So far no such resonant decay has been observed. Both the ATLAS and the CMS collaboration have searched for double Higgs final states produced by a narrow resonance decaying into bbbb [6,7] and γγbb [8,9]. These and other final states in double Higgs production, will be one of the priorities in Higgs physics during the LHC run 2. Furthermore, the Higgs couplings extracted so far from the 7 and 8 TeV data show no large deviations from the Standard Model (SM) expectations and point towards a very SM-like Higgs sector.
In the SM, the double Higgs production cross section is rather small and even at the high luminosity stage of the LHC it will be extremely challenging to measure the triple Higgs coupling. However, in many beyond the SM (BSM) scenarios, resonant decays are possible due to the existence of other scalar particles in the theory. This increases the prospects for observing a final state with two scalars. The simplest extension where a di-Higgs final state would be detectable at the LHC is the singlet extension of the SM, where a hypercharge zero singlet is added to the scalar field content of the model. When the singlet is real one either obtains a new scalar mixing with the Higgs boson or the minimal model for dark matter [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. The model can also accommodate electroweak baryogenesis by allowing a strong first-order phase transition during the era of EWSB [26][27][28][29][30], if the singlet is complex. In this case the spectrum features three Higgs bosons, either all visible or one being associated with dark matter. In the present work we will focus mainly on the complex version of the singlet extension although we will also discuss the case of a real singlet. The collider phenomenology of the singlet model, which can lead to some distinctive signatures at the LHC, has been previously discussed in [31][32][33][34][35][36][37][38][39][40].
Both the real and the complex singlet extensions, in the minimal versions that we will discuss, can have at least two phases. Namely, a Z 2 -symmetric phase with a dark matter candidate and a broken phase where the singlet component(s) mix with the neutral scalar field of the SM doublet. Therefore these models can be used as simple benchmarks for resonant double Higgs production pp → H → hh, where H generically denotes a new heavy scalar and h the SM-like Higgs boson. Furthermore, the broken phase of the complex model also allows for a decay of a heavy scalar h 3 into two other scalars h 1 , h 2 of different mass, h 3 → h 2 + h 1 . This particular decay is common in other extensions of the SM, such as the Complex Two-Higgs Doublet Model (C2HDM) and the Nextto-Minimal Supersymmetric extension of the SM (NMSSM) [41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56]. In the NMSSM a complex superfieldŜ is added to the Minimal Supersymmetric Model field content allowing for a dynamical solution of the µ problem when the singlet field acquires a non-vanishing vacuum expectation value (VEV). NMSSM Higgs sectors that feature light Higgs states in the spectrum can lead to sizeable decay widths for Higgs decays into a pair of lighter Higgs bosons. Higgs-to-Higgs decays in the NMSSM have been studied in [57][58][59][60][61][62] and several recent studies have investigated the production of Higgs pairs in the NMSSM [57,59,[63][64][65]. These processes involve the trilinear Higgs self-coupling, which is known in the NMSSM including the full one-loop corrections [57] and the two-loop O(α t α s ) corrections [62] both in the real and in the complex NMSSM. A recent study within the C2HDM on all possible CP-violating scalar decays and the chances of detecting direct CP violation at the LHC can be found in [66]. It is interesting to point out that decays of the type h i → h j Z are forbidden in the singlet extension because the model is CP-conserving and has no CP-odd states. In the real NMSSM (and also in the CP-conserving 2HDM) the decay is possible provided CP(h i ) = − CP(h j ) while in CP-violating models, like the C2HDM and the complex NMSSM, h i → h j Z is always possible if kinematically allowed. Note that combinations of scalar decays into pairs of scalars together with h i → h j Z allow us not only to probe the CP nature of the model but can also be used as a tool to identify the CP quantum numbers of the scalars of the theory [66]. Finally there is an interesting scenario where all three scalars h i are detected in h i → ZZ. If the potential is CP-conserving then all three scalars are CP-even. This is exactly what happens in the CxSM and in the NMSSM but not in other simple extension of the the SM like the 2HDM 5 .
In order to calculate the branching ratios of the singlet models we have modified HDECAY [68,69] to include the new vertices. This new tool, sHDECAY, can be used to calculate all branching ratios and widths in both the real and the complex singlet extensions of the SM that we discuss. In each extension, both the symmetric phase and the broken phase were considered leading to a total of four models. The detailed description of sHDECAY can be found in appendix A 6 .
In this study we define benchmark points for the LHC run 2 both for the complex singlet model (CxSM) and the real singlet model (RxSM). We are guided by phenomenological and theoretical goals. The first is to maximize at least one of the Higgs-to-Higgs decays. The second is related to the stability of the theory at higher orders. As discussed in [70], the radiative stability of the model, at two loops, up to a high-energy scale, combined with the constraint that the 125 GeV Higgs boson found at the LHC is in the spectrum, forces the new scalar to be heavy. When we include all experimental constraints and measurements from collider data, dark matter direct detection experiments and from the Planck satellite and in addition force stability at least up to the Grand Unification Theory (GUT) scale, we find that the lower bound on the new scalar is about 170 GeV [70]. Finally, whenever a dark matter candidate exists we force it to fully saturate the dark matter relic density measured by Planck as well as to make it consistent with the latest direct detection bounds.
The renormalisation of the real singlet model has recently been addressed in [71,72]. In [71] the focus was on the modifications of the couplings of the SM-like Higgs to fermions and gauge bosons while in [72] the goal was to determine the electroweak corrections to the heavy Higgs decay, H, into the two light ones, h. In the former it was found that the electroweak corrections are at most of the order of 1 % and this maximal value is attained in the limit where the theory decouples and becomes indistinguishable from the SM. In the latter it was found that the corrections to the triple scalar vertex (Hhh) are small, typically of a few percent, once all theoretical and experimental constraints are taken into account. Furthermore it was concluded that electroweak corrections are stable for changes of the renormalisation scheme. There are so far no calculations of electroweak corrections in the complex singlet extension of the SM.
We end the introduction with a word of caution regarding interference effects in the real singlet model. The problem of interference in BSM scenarios has been raised in [73] for the particular 5 Note however that if CP-violation is introduced via the Yukawa Lagrangian, the pseudoscalar can decay to ZZ although in principle with a very small rate as shown in [67] for a particular realisation of the 2HDM. 6 The program sHDECAY can be downloaded from http://www.itp.kit.edu/∼maggie/sHDECAY/.
case of a real singlet extension of the SM. As shown in [73], the SM-like Higgs contribution to the process gg → h * , H ( * ) → ZZ → 4l is non-negligible outside the non-SM scalar (H) peak region and interference effects can be quite large. These effects were shown to range from O(10%) to O(1) for the integrated cross sections at the LHC at a center-of-mass energy of 8 TeV [74]. Moreover, these effects also modify the line shape of the heavier scalar. However, as shown in [74] judicious kinematical cuts can be used in the analysis to reduce the interference effects to O(10%).
Recently, interference effects were studied in the process gg → h * , H ( * ) → hh at next-to-leading order (NLO) in QCD [75]. It was found that the interference effects distort the double Higgs invariant mass distributions. Depending on the heavier Higgs mass value, they can either decrease by 30% or increase by 20%. Further, it was shown that the NLO QCD corrections are large and can significantly distort kinematic distributions near the resonance peak. This means that any experimental analysis to be performed in the future should take these effects into account. The paper is organized as follows. In Sect. 2 we define the singlet models, CxSM in Sect. 2.1 and RxSM in Sect. 2.2, and we describe, in Sect. 2.3, the constraints that are applied. In the experimental constraints particular emphasis is put on the LHC run 1 data. The NMSSM is briefly introduced in Sect. 2.4 together with the constraints that are used in the scan. In Sect. 3 we present our numerical analysis and results. First, the possibility of SM Higgs boson production from chain decays in the CxSM and RxSM after run 1 is discussed in 3.1. Afterwards Higgs-to-Higgs decays at run 2 are investigated in 3.2 and the related features of the RxSM are compared to those of the CxSM in Sect. 3.2.1 before moving on to the comparison of the CxSM with the NMSSM, which is performed in Sect. 3.2.2. Finally in Sect. 3.3 benchmark points are presented for the CxSM and the RxSM for the LHC run 2. Our conclusions are given in Section 4. In appendix A the implementation of the singlet models in sHDECAY is described and sample input and output files are given. In appendix B we list the expressions for the scalar vertices in the singlet extensions used in the scalar decays into pairs of scalars.

Models and Applied Constraints
In this section we define the reference complex and real singlet models. These will be analysed to define benchmarks for the various allowed kinematical situations at the LHC run 2, with a centerof-mass (c.m.) energy of 13 TeV. We review how EWSB proceeds consistently with the observed Higgs boson and define the couplings of the scalars of each theory with the SM particles as well as within the scalar sector. We then describe the various constraints that we apply to these singlet models. Subsequently the NMSSM is briefly introduced, and the ranges of the parameters used in the NMSSM scan are described, together with the applied constraints.

The CxSM
The main model discussed is this work is a simple extension of the SM where a complex singlet field with hypercharge zero, is added do the SM field content. All interactions are determined by the scalar potential, which can be seen as a model with a U (1) global symmetry that is broken softly.
The most general renormalisable scalar potential, with soft breaking terms with mass dimension up to two, is given by where the soft breaking terms are shown in parenthesis and the doublet and complex singlet are, respectively, Here v ≈ 246 GeV is the SM Higgs VEV, and v S and v A are, respectively, the real and imaginary parts of the complex singlet field VEV.
In [76], the various phases of the model were discussed. In order to classify them it is convenient to treat the real and the imaginary components of the complex singlet as independent, which is equivalent to building a model with two real singlet fields. We focus on a version of the model that is obtained by requiring a Z 2 symmetry for the imaginary component A (this is equivalent to imposing a symmetry under S → S * ). As a consequence of this symmetry, the soft breaking couplings must be both real, i.e. a 1 ∈ R and b 1 ∈ R. Observe that m, λ, δ 2 , b 2 and d 2 must be real parameters for the potential to be real.
By analysing the minimum conditions one finds two possible phases that are consistent with EWSB triggering the Higgs mechanism. They are: • v A = 0 and v S = 0, in which case mixing between the doublet field h and the real component s of the singlet field occurs, while the imaginary component A ≡ a becomes a dark matter candidate. We call this the symmetric or dark matter phase.
• v S = 0 and v A = 0, which we call the broken phase, with no dark matter candidate and mixing among all scalars.
The model phases are summarized in table 1.

Phase Scalar content
VEVs at global minimum Symmetric (dark) 2 mixed + 1 dark S = 0 and A = 0 Broken ( Z 2 ) 3 mixed S = 0 and A = 0 As discussed in [76], simpler models can be obtained with the same field content by imposing extra symmetries on the potential. The exact U (1)-symmetric potential has a 1 = b 1 = 0, leading to either one or to two dark matter candidates depending on the pattern of symmetry breaking. One can also impose a separate Z 2 symmetry for S and A that forces a 1 = 0 and b 1 ∈ R and, again, gives rise to the possibility of having one or two dark matter candidates. However, from the phenomenological point of view, the model presented here covers all possible scenarios in terms of the accessible physical processes to be probed at the LHC in a model with three scalars. In fact, even simpler models, like the real singlet extension of the SM discussed below, have some similarities with this version of the CxSM in what concerns the planned searches for the next LHC runs. However, typically, the allowed parameter space of different models can also be quite different, implying that the number of expected events may allow the exclusion of one model but not of another qualitatively similar one.

Physical states and couplings
To obtain the couplings of the scalars to the SM particles, we define the mass eigenstates as h i (i = 1, 2, 3). They are obtained from the gauge eigenstates h, s and a through the mixing matrix R with and m 1 ≤ m 2 ≤ m 3 denoting the masses of the neutral Higgs particles. The mixing matrix R is parametrized as with s i ≡ sin α i and c i ≡ cos α i (i = 1, 2, 3) and In the dark matter phase, α 2 = α 3 = 0, the a field coincides with A and it is the dark matter candidate, which does not mix with h nor with s. The couplings of h i to the SM particles are all modified by the same matrix element R i1 . This means that, for any SM coupling λ (p) h SM where p runs over all SM fermions and gauge bosons, the corresponding coupling in the singlet model for the scalar h i is given by The coupling modification factor is hence independent of the specific SM particle to which the coupling corresponds.
In the dark matter phase the same arguments apply with R i1 = (R 11 , R 21 , 0). The state i = 3 then corresponds to the dark matter candidate, which does not couple to any of the SM fermions and gauge bosons.
The couplings that remain to be defined are the ones of the scalar sector. They are read directly from the scalar potential, Eq. (2), after replacing the expansion about the vacua, Eq. (3). For the purpose of our analysis we only need the scalar triple couplings. We define them according to the following normalisation (i, j, k = 1, 2, 3): These couplings determine, up to the phase space factor, the leading order expressions for the Higgs-to-Higgs decay widths. In appendix B.1 we provide the expressions for these couplings as well as for the quartic couplings, for completeness. −π/2 π/2 a 1 (GeV 3 ) −10 8 0

Parameters of the model
In Sect. 3 we will present scans over the parameter space of this model and we will identify benchmark points that represent various qualitatively different regions in this space. Regardless of the phase of the model under consideration, the model always has seven parameters. The particular choice of independent parameters is just a matter of convenience. For the broken phase we choose the set {α 1 , α 2 , α 3 , v, v S , m 1 , m 3 } and for the dark phase the set {α 1 , v, v S , a 1 , m 1 , m 2 , m 3 ≡ m A }.
In our scans, all other dependent parameters are determined internally by the ScannerS program [76,77] according to the minimum conditions for the vacuum corresponding to the given phase. This program provides tools to automatise the parameter space scans of generic scalar extensions of the SM. It also contains generic modules for testing local vacuum stability, to detect symmetries, and library interfaces to the codes we have used to implement the constraints (described later in Sect. 2.3). Apart from v, which is obtained internally from the Fermi constant G F , sHDECAY uses the same sets of input parameters and all other dependent parameters are calculated internally from these. In Table 2 we present the values and ranges of parameters that were allowed in the scan, before applying constraints. Note that in the broken phase the mass that is not an input is obtained from the following relation

The RxSM
In this study our main focus is the CxSM model in which all possibilities exist for the decays of a scalar into two different or identical scalars, in a theory with three CP-even scalars. The special case when the two scalars in the final state of the decay are identical is, however, also allowed in simpler models with just two scalar degrees of freedom. The simplest such model is the real singlet model (RxSM). In fact, the broken phase of the CxSM, the dark phase of the CxSM and the broken  RxSM all share this possibility. In order to compare these models in our analysis we now briefly summarise the RxSM. This model is obtained by adding a real singlet S with a Z 2 symmetry (S → −S) to the SM. Then the most general renormalisable potential reads where m, λ, λ HS , m S and λ S are real. Electroweak symmetry breaking, with a vacuum consistent with the Higgs mechanism occurs when the following vacuum expectation values are chosen: Again v ≈ 246 GeV is the SM Higgs VEV, and v S is the singlet VEV. In this model, we are primarily interested in the broken phase, v S = 0, where we use the notation m 1 and m 2 for the scalar states, ordered in mass, that mix two-by-two. The mixing matrix is precisely the same as the sub-block responsible for the mixing in the dark phase of the complex singlet when α 1 ≡ α and α 2 = α 3 = 0.
In the symmetric phase, only the observed Higgs boson with mass m 1 ≡ m h 125 is visible. For the dark matter candidate we have the mass m 2 ≡ m D . In the broken phase, on the contrary, both scalars are visible with one of them corresponding to the SM-like Higgs boson. In the following we will focus on the broken phase since it allows for the Higgs-to-Higgs decays that we want to use in the comparison with the CxSM.
Regarding couplings everything is similar to the CxSM, i.e. the couplings of the various scalars are controlled by the same rule, Eq. (8), with the reduced two-by-two mixing matrix. The couplings among the scalars are again directly read from the potential of the model, Eq. (11), after expanding about the vacua, Eq. (12), and they are given in Appendix B.2.
Finally, for the scans over the parameter space of this model we note that the model has five independent parameters, which are chosen as {α 1 ≡ α, m 1 , m 2 , v, v S }. In the implementation of the decays for the broken phase in sHDECAY we choose the same set, with the SM VEV, v, determined from the Fermi constant G F . All other dependent parameters are determined internally by sHDECAY, and by ScannerS according to the minimum conditions for the vacuum corresponding to the given phase. In Table 3 we present the values and ranges of parameters that were allowed in the scan before applying constraints.
For completeness we note that the input parameters for the implementation of the dark matter phase in sHDECAY, which we do not analyse here, are chosen as m 1 ≡ m h 125 , m 2 ≡ m D , λ S , m 2 S , v (v is again determined from G F ).

Constraints Applied to the Singlet Models
To restrict the parameter space of the singlet models we apply various theoretical and phenomenological constraints. These have been described in [76] and recently updated in [70]. We will only discuss them briefly expanding mostly on the constraints from collider physics, which is the focus of our discussion on benchmarks for the LHC run 2.

Theoretical constraints
We apply both to the CxSM and the RxSM the constraints that: i) the potential must be bounded from below; ii) the vacuum we have chosen must be a global minimum and iii) that perturbative unitarity holds. The first two are implemented in the ScannerS code [76,77] using the relevant inequalities for the CxSM and the RxSM. Perturbative unitarity is tested using a general internal numerical procedure in ScannerS, which can perform this test for any model, see [76,77].

Dark matter constraints
In the dark phase of the CxSM, we compute the relic density Ω A h 2 with micrOMEGAS [78] and use it to exclude parameters space points for which Ω A h 2 is larger than Ω c h 2 +3σ, where Ω c h 2 = 0.1199± 0.0027 is the combination of the measurements from the WMAP and Planck satellites [79,80] and σ denotes the standard deviation. As for bounds on the direct detection of dark matter, we compute the spin-independent scattering cross section of weakly interacting massive particles (WIMPS) on nucleons also with micrOMEGAS with the procedure described in [76] and the points are rejected if the cross section is larger than the upper bound set by the LUX2013 collaboration [81].

Electroweak precision observables
We also apply a 95% exclusion limit from the electroweak precision observables S, T, U [82,83]. We follow precisely the procedure reviewed in [70] and implemented in the ScannerS code and therefore do not repeat the discussion here.

Collider constraints
The strongest phenomenological constraints on models with an extended Higgs sector typically arise from collider data, most importantly from the LHC data. On one hand, one must ensure that one of the scalar states matches the observed signals for a Higgs boson with a mass of 125 GeV. The masses of the other two scalars can either be both larger, both smaller and one larger and one smaller than that of the discovered Higgs.
On the other hand, all other new scalars must be consistent with the exclusion bounds from searches at various colliders, namely the Tevatron, LEP and LHC searches. The strongest constraints on the parameter space come from the measurements of the rates of the 125 GeV Higgs performed by ATLAS and CMS during the LHC run 1. The experimental bounds are imposed with the help of ScannerS making use of its external interfaces with other codes. We have applied 95% C.L. exclusion limits using HiggsBounds [84]. As for consistency with the Higgs signal measurements we test the global signal strength of the 125 GeV Higgs boson with the latest combination of the ATLAS and the CMS LHC run 1 datasets, i.e. µ 125 = 1.09 ± 0.11 [85]. In the remainder of the analysis, for the singlet models, we fix the SM-like Higgs mass to 125.1 GeV, which is the central value 7 of the ATLAS/CMS combination reported in [86]. Note, also, that in our scans we required the non-SM-like Higgs masses to deviate by at least 3.5 GeV from 125.1 GeV, in order to avoid degenerate Higgs signals.
HiggsBounds computes internally various experimental quantities such as the signal rates Here σ New (h i ) denotes the production cross section of the Higgs boson, h i , in the new model under consideration and σ SM (h i ) the SM production cross section of a Higgs boson with the same mass. Similarly, BR New is the branching ratio for h i in the new model to decay into a final state with SM particles X SM , and BR SM the corresponding SM quantity for a Higgs boson with the same mass. For the computation of the rates HiggsBounds requires as input for all scalars: the ratios of the cross sections for the various production modes, the branching ratios and the total decay widths. In singlet models, at leading order (and also at higher order in QCD), the cross section ratios are all simply given by the suppression factor squared, R 2 i1 , see Eq. (8). Regarding the branching ratios and total decay widths we have used the new implementation of the CxSM and the RxSM in sHDECAY as described in appendix A.
In the analyses of Sects. 3.2 and 3.3 the h i production cross sections through gluon fusion, gg → h i , are needed at 13 TeV c.m. energy. We have computed these both with the programs SusHi [87] and HIGLU [88] at next-to-next-to-leading order (NNLO) QCD and with higher order electroweak corrections consistently turned off. The results of the two computations agreed within the numerical integration errors and within parton density function uncertainties. In the analysis we have used the results from HIGLU.

The NMSSM
Supersymmetric models require the introduction of at least two Higgs doublets. In the NMSSM this minimal Higgs sector is extended by an additional complex superfieldŜ. After EWSB the NMSSM Higgs sector then features seven Higgs bosons. These are, in the CP-conserving case, three neutral CP-even, two neutral CP-odd and two charged Higgs bosons. The NMSSM Higgs potential is derived from the superpotential, the soft SUSY breaking terms and the D-term contributions. The scale-invariant NMSSM superpotential in terms of the hatted superfields reads where we have included only the third generation fermions. While the first term replaces the µterm µĤ uĤd of the MSSM superpotential, the second term, cubic in the singlet superfield, breaks the Peccei-Quinn symmetry so that no massless axion can appear. The last three terms represent the Yukawa interactions. The soft SUSY breaking Lagrangian gets contributions from the scalar mass parameters for the Higgs and the sfermion fields. In terms of the fields corresponding to the complex scalar components of the superfields the Lagrangian reads The contributions from the trilinear soft SUSY breaking interactions between the sfermions and the Higgs fields are comprised in The soft SUSY breaking Lagrangian that contains the gaugino mass parameters is given by We will allow for non-universal soft terms at the GUT scale. After EWSB the Higgs doublet and singlet fields acquire non-vanishing VEVs. By exploiting the three minimisation conditions of the scalar potential, the soft SUSY breaking masses squared for H u , H d and S in L mass are traded for their tadpole parameters. The Higgs mass matrices for the three scalar, two pseudoscalar and the charged Higgs bosons are obtained from the tree-level scalar potential after expanding the Higgs fields about their VEVs v u , v d and v s , which we choose to be real and positive, The 3 × 3 mass matrix squared, M 2 S , for the CP-even Higgs fields is diagonalised through a rotation matrix R S yielding the CP-even mass eigenstates H i (i = 1, 2, 3), The H i are ordered by ascending mass, In order to obtain the CP-odd mass eigenstates, A 1 and A 2 , first a rotation R G is performed to separate the massless Goldstone boson G, which is followed by a rotation R P to obtain the mass eigenstates which are ordered such that The tree-level NMSSM Higgs sector can be parametrised by the six parameters We choose the sign conventions for λ and tan β such that they are positive, whereas κ, A λ , A κ and µ eff can have both signs. In order to make realistic predictions for the Higgs masses we have taken into account the higher order corrections to the masses. Through these corrections also the soft SUSY breaking mass terms for the scalars and the gauginos as well as the trilinear soft SUSY breaking couplings enter the Higgs sector.
In the NMSSM several Higgs-to-Higgs decays are possible if kinematically allowed. In the CP-conserving case the heavier pseudoscalar A 2 can decay into the lighter one A 1 and one of the two lighter Higgs bosons H 1 or H 2 8 . The heavier CP-even Higgs bosons can decay into various combinations of a pair of lighter CP-even Higgs bosons or into a pair of light pseudoscalars. Hence, in principle we can have, Not all of these decays lead to measurable rates, however. Furthermore, if the decay widths of the SM-like Higgs boson into lighter Higgs pairs are too large, the compatibility with the measured rates is lost. In order to investigate what rates in the various SM final states can be expected from Higgs-to-Higgs decays an extensive scan in the NMSSM parameter spaces needs to be performed. We have generated a new sample of NMSSM scenarios, following a procedure similar to an earlier publication [59], to obtain the decay rates into Higgs pairs. The scan ranges and the applied criteria shall be briefly discussed here and we refer the reader to [59] for further details.
For the scan, tan β, the effective µ parameter and the NMSSM specific parameters λ, κ, A λ and A κ have been varied in the ranges Perturbativity constraints have been taken into account by applying the rough constraint We have varied the trilinear soft SUSY breaking couplings of the up-and down-type quarks and the charged leptons, A U , A D and A L with U ≡ u, c, t, D ≡ d, s, b and L ≡ e, µ, τ , independently in the range For the soft SUSY breaking right-and left-handed masses of the third generation we choose For the first two generations we take The gaugino soft SUSY breaking masses are chosen to be positive and varied in the ranges The chosen parameter ranges entail particle masses that are compatible with the exclusion limits on SUSY particle masses [89][90][91][92][93][94][95][96][97][98][99][100][101][102] and with the lower bound on the charged Higgs mass [103,104]. With the help of the program package NMSSMTools [105][106][107] we checked for the constraints from lowenergy observables and computed the input necessary for HiggsBounds to check for consistency with the LEP, Tevatron and the LHC run 1 latest exclusion limits from searches for new Higgs bosons. Details are given on the webpage of the program 9 . Via the interface with micrOMEGAS [108][109][110][111] also the compatibility with the upper bound of the relic abundance of the lightest neutralino as the NMSSM dark matter candidate with the latest PLANCK results was verified [79]. Among the points that are compatible with the NMSSMTools constraints, only those are kept that feature a SM-like Higgs boson with mass between 124 and 126 GeV. This can be either H 1 or H 2 and will be called h in the following. In addition, according to our criteria, a Higgs boson is SM-like if its signal rates into the W W , ZZ, γγ, bb and τ τ final states are within 2 times the 1σ interval around the respective best fit value. We use the latest combined signal rates and errors reported by ATLAS and CMS in [85], in a 6-parameter fit. In this fit, the five rates mentioned above are for fermion mediated production, so in their computation we approximate the inclusive production cross section by the dominant production mechanisms, gluon fusion and bb annihilation. The sixth parameter of the fit is the ratio between the rate for vector mediated production and the fermion mediated production, which we also check to be within 2σ. The cross section for gluon fusion is obtained by multiplying the SM gluon fusion cross section with the ratio between the NMSSM Higgs decay width into gluons and the corresponding SM decay width at the same mass value. The two latter rates are obtained from NMSSMTools at NLO QCD. The SM cross section was calculated at NNLO QCD with HIGLU [88]. For the SM-like Higgs boson this procedure approximates the NMSSM Higgs cross section, computed at NNLO QCD with HIGLU, to better than 1%. As for the cross section for bb annihilation we multiply the SM bb annihilation cross section with the effective squared bb coupling obtained from NMSSMTools. For the cross section values we use the data from [112], which was produced with the code SusHi [87]. Note that, similarly to the CxSM analysis, we have excluded degenerate cases where other non-SM-like Higgs bosons would contribute to the SM-like Higgs signal, by requiring that the masses of the non-SM-like Higgs bosons deviate by at least 3.5 GeV from 125.1 GeV. Finally, note that in [59] the branching ratios obtained with NMSSMTools were cross-checked against the ones calculated with NMSSMCALC [113,114]. There are differences due to the treatment of the radiative corrections to the Higgs boson masses 10 , as well as due to the more sophisticated and up-to-date inclusion of the dominant higher order corrections to the decay widths and the consideration of off-shell effects in NMSSMCALC. The overall picture, however, remains unchanged. 9 http://www.th.u-psud.fr/NMHDECAY/nmssmtools.html 10 In NMSSMTools the full one-loop and the two-loop O(αs(α b +αt)) corrections at vanishing external momentum [115] are included. NMSSMCALC provides both for the real and the complex NMSSM the full one-loop corrections including the momentum dependence in a mixed DR−on-shell renormalisation scheme [116][117][118] and the two-loop O(αsαt) corrections at vanishing external momentum. For the latter, in the (s)top sector the user can choose between on-shell and DR renormalisation conditions. For a comparison of the codes, see also [119].

Numerical Analysis
In this section, we analyse the consequences of the Higgs-to-Higgs decays on the allowed parameter space after applying the constraints described in Sects. 2.3 and 2.4. First we analyse the importance of chain decays for the production of the SM-like Higgs boson (Sect. 3.1), then we discuss the prospects of using such Higgs-to-Higgs decays to distinguish between the various models (Sect. 3.2) and, finally, we provide a set of benchmark points for the singlet models in Sect. 3.3.

Higgs Boson Production from Chain Decays in the CxSM and RxSM
In models with extended Higgs sectors, where additional scalar states exist that can interact with the SM Higgs boson, the signal rates for a given SM final state X SM , Eq. (13), do not account for the total process leading to X SM through Higgs boson decays. This is the case for the singlet models, where the same final state may be reached through an intermediate step with a heavy scalar, h i , decaying into two lighter scalars h j , h k (which may be different or not), that finally decay into a final state X SM .
In the following we will investigate the question: Can such resonant decays of a heavy Higgs boson into a final state containing a reconstructed SM-like Higgs boson (which is the one being observed at the LHC) compete with the direct production of the observed Higgs state [120]? We will call the former production mechanism chain production. For a given channel X SM in which the SM-like Higgs is observed, the total rate is then defined as where Here N h i ,h 125 is the expected number of h 125 Higgs bosons produced in the decay of h i , which, for a model with up to 3 scalars, is given by where P h i ,n×h 125 are the probabilities of producing n of the observed Higgs boson h 125 from the decay of h i , with i = j and m h j < m h i . Replacing in Eq. (31) the branching ratios by the partial and total widths and exploiting the fact that the production cross section and the direct decay widths of h i into SM final states are both modified by R 2 i1 , compared to the SM, we arrive at It is important to note, from Eq. (34), that µ C h 125 is in fact independent of the SM final state X SM . This simplifies the analysis because we can study this global quantity instead of each final state individually.
In the remainder of this section we compare the relative importance of direct production to chain production in the singlet models through the ratio µ C h 125 /µ T h 125 . This measures the fraction of Higgs events from chain production with respect to the total number of Higgs events. The allowed parameter space for the CxSM model after LHC run 1 was recently analysed in detail in [70]. In this study it was found that, in the broken phase, all kinematical possibilities are still allowed, i.e. the observed Higgs boson can be any of the three scalars and the corresponding couplings may deviate considerably from the SM. Similarly, in the dark phase of the CxSM or the broken phase of the RxSM, all kinematically different situations are still allowed for the two visible scalars.
Here we point out that even for the interpretation of the LHC run 1 data, it is important to take into account the contributions from chain decays in the measurements of the SM Higgs signal rates. In some cases, these contribution can be up to 15% of the total signal rate. In singlet models, where the direct decay signal rate is simply suppressed relative to the SM limit, this can lead to viable points that would naively be excluded if this contribution is not included. On the other hand it is interesting to investigate which points are still allowed by the LHC run 1 data and, simultaneously, have the largest possible chain decay contributions. In this way, a new heavy scalar may be found indirectly at run 2 in decays into a SM Higgs pair (in any of the singlet models we discuss) or into a SM Higgs with a new light scalar (in the broken phase of the CxSM).
In Fig. 1 we present a sample of points generated for the broken phase of the CxSM in the scenario where the observed Higgs boson is the lightest scalar (top panels) and the next-to-lightest scalar (bottom panels). We present projections against the masses of the two new scalars that can be involved in the chain decay contributions. We overlay three layers of points for which the total signal rate µ T h 125 is, respectively, within 3σ (red), 2σ (blue) and 1σ (green) of the LHC run 1 best fit point for the global signal strength. The vertical axis shows the fraction of the signal rate that is due to the chain decay contribution. In the upper plots, where h 1 ≡ h 125 , chain decays become possible above m h 3 = 250 GeV, cf. Fig. 1 (upper left). In the upper right plot points exist when m h 2 129 GeV due to the minimum required distance of 3.5 GeV from m h 125 . In the lower two plots, h 2 ≡ h 125 , so that with the lowest h 1 mass value of 30 GeV in our scan, cf. Fig. 1 (lower right), chain decays come into the game for m h 3 155 GeV. It can be inferred from the top panels, that when the Higgs boson is the lightest scalar, the chain decay contribution can be up to ∼ 6% (∼ 12%) at 2σ (3σ). We have also produced plots where we apply the 2σ (3σ) cut on the rate for direct h 125 production and decay, i.e. on µ h 125 instead of µ T h 125 . In that case we observe a reduction of the maximum allowed µ C h 125 /µ T h 125 and the numbers will change to ∼ 4% (∼ 9%). This means that imposing a constraint on the parameter space without including the chain decay contribution would be too strong and points with larger chain decays at the LHC run 2 would in fact still be allowed. The kink observed in Fig. 1 (upper right) at m h 2 = 250 GeV is due to the opening of the decay channel h 2 → h 125 +h 125 which now also adds to the chain decay contributions from h 3 → h 125 + h 125 and h 3 → h 125 + h 2 . The observed tail for large masses stems from the constraints from electroweak precision observables (EWPOs). In this region both m h 2 and m h 3 are large so that EWPOs force the factors R 2 i1 for the SM couplings of the two heavy Higgs bosons to be small, thus suppressing the contribution from chain decays. In the upper left plot on the other hand, the points for large m h 3 also include the cases where m h 2 can be small, which can then have a larger modification factor R 2 21 without being in conflict with EW precision data. Overall, the shape of the three different σ regions in the various panels is the result of an interplay between the kinematics and the applied constraints. To a certain extent the structure can be directly related to the exclusion curves from the collider searches imposed by HiggsBounds. This is e.g. the case for the peaks at m h 2 ∼ 170 GeV and m h 2 ∼ 280 GeV, as we have checked explicitly by generating a sample of points without imposing collider constraints. 11 In the bottom plots of Fig. 1, where h 125 is the next-to-lightest Higgs boson, the relative fraction of the chain contribution on the total rate is much smaller with at most ∼ 2% (∼ 4%) for the 2σ 11 A similar behaviour was found in Fig. 8 (top and bottom right) of Ref. [70]. Each layer corresponds to points that fall inside the region where the total signal rate µ T h125 is within 3σ (red), 2σ (blue) and 1σ (green) from the best fit point of the global signal rate using all LHC run 1 data from the ATLAS and CMS combination [85]. (3σ) region. The reason is that here only the decays h 3 → h 1 + h 2 , h 2 + h 2 contribute to the chain decays, while in the scenario with h 125 ≡ h 1 also the h 2 → h 1 + h 1 decay is possible. The strong increase at m h 3 = 250 GeV in the lower left panel is due to the opening of the decay h 3 → h 2 + h 2 , and the decrease in the number of points for large values of m h 3 is due to EWPOs. The remainder of the shapes of the three σ region again is explained by the interplay between kinematics and applied constraints, namely the exclusion curves from the Higgs data and their particular shape. In Fig. 2 we display the fraction of chain decays for the scenarios where the Higgs is the lightest visible scalar in the dark phase of the CxSM (left) and in the broken phase of the RxSM (right). Both models show a similar behaviour. Above the kinematic threshold for the chain decay, at m h 2 = 250 GeV, the fraction of chain decays increases to about 14% in the CxSM (left plot) and 17% in the RxSM (right plot) for m h 2 ≈ 270 GeV and a total signal rate within 3σ of the measured value of the global signal rate µ. For the total rate within 2σ these numbers decrease to ∼ 7% (left plot), respectively, ∼ 9% (right plot). As before, the suppression of the fraction of chain decays at large masses m h 2 can be attributed to EW precision constraints, and the second peak at larger h 2 masses together with the overall shape is explained by the combination of kinematics and constraints, in particular the Higgs exclusion curves.
Note, that in the CxSM the density of points is lower simply because the parameter space is higher dimensional and the independent parameters that were used to sample it do not necessarily generate a uniform distribution of points in terms of the fraction of chain decays. The samples have been generated in both cases with a few million points.

Higgs-to-Higgs Decays at the LHC Run2
Many extensions of the SM allow for the decay of a Higgs boson into two lighter Higgs states of different masses. Such a decay is not necessarily a sign of CP violation. In fact, in the broken phase of the CxSM, although all three scalars mix, they all retain the quantum numbers of the scalar in the doublet, i.e. they are all even under a CP transformation. On the contrary, in the C2HDM it could be a signal of CP violation. Furthermore, even in CP-conserving models, the CP number of the new heavy scalars may not be accessible at an early stage if they are to be found at the LHC run 2. Thus, models that have a different theoretical structure and Higgs spectrum may in fact look very similar at an early stage of discovery.
In section 3.2.2 we will perform a comparison between the broken phase of the CxSM and the NMSSM, which contains similar possibilities in terms of scalar decays into scalars, if the CP numbers are not measured. We focus on resonant decays, which allow to probe the scalar couplings of the theory, and, in many scenarios, may provide an alternative discovery channel for the new scalars. First, however, we compare the Higgs-to-Higgs decay rates at the LHC run 2 for the real singlet extension in its broken phase with the CxSM. The results of this comparison will then allow us to draw meaningful conclusions in the subsequent comparison between the NMSSM and CxSM-broken.

Comparison between the RxSM-broken and the Two Phases of the CxSM
We discuss the broken phase of the RxSM because it is the only one allowing for Higgs-to-Higgs decays in this model. Depending on the mass of the non-SM-like additional Higgs boson, we have two possible scenarios for the comparison of RxSM-broken and CxSM in its symmetric (dark) and its broken phase: 2) SM-like Higgs boson decaying into two identical Higgs bosons: Here the non-SM-like Higgs is lighter than the SM-like Higgs boson, and we denote it by ϕ. If it is light enough the decay h 125 → ϕ + ϕ is possible. This is the only case that can be realized in RxSM-broken and CxSM-dark. In contrast, the decay possibilities of CxSM-broken are given by For simplicity again 4b final states are investigated. Figure 3 (left) shows the decay rates for case 1) in the broken phase (blue points) and the dark matter phase (green) of the CxSM, and in the RxSM-broken (red). It demonstrates, that the maximum rates in the RxSM-broken exceed those of the CxSM, although the differences are not large.
The results for case 2) are displayed in Fig. 3 (right). Again the maximum rates in RxSMbroken are larger than those in the CxSM. In the dark matter phase the maximum rates are not much smaller, while the largest rates achievable in CxSM-broken lie one order of magnitude below those of the RxSM. We have verified that the larger rates allowed for the models with two-by-two mixing (CxSM-dark and RxSM) result from their different vacuum structure. The larger rates can be traced back to the branching ratio for the Higgs-to-Higgs decay, BR(h 125 → ϕ + ϕ), which differs, from model to model, in its allowed parameter space for the new scalar couplings of the theory.

Comparison of CxSM-broken with the NMSSM
We now turn to the comparison of the Higgs-to-Higgs decay rates that can be achieved in the broken phase of the complex singlet extension with those of the NMSSM. We will focus, in our discussion, on the comparison of final state signatures that are common to both models. We will not consider additional decay channels, that are possible for the NMSSM Higgs bosons, such as decays into lighter SUSY particles, namely neutralinos or charginos, or into a gauge and Higgs boson pair. These decays would add to the distinction of the CxSM-broken from the NMSSM, if they can be identified as such. Depending on the NMSSM parameter points, various such decay possibilities may be possible and will require a dedicated analysis that is beyond the scope of this paper. Our approach here is a different one. We assume we are in a situation where we have found so far only a subset of the Higgs bosons common to both models, where we have not observed any non-SM final state signatures yet and where we do not have information yet on the CP properties of the decaying Higgs boson. Additionally, we assume that we do not observe any final state signatures that are unique in either of the models. 12 We then ask the question: If one focusses on Higgs-to-Higgs decays only in final states that are common to both models, will it be possible to tell the CxSM-broken from the NMSSM based on the total rates? The discovery of additional non-CxSM Higgs bosons and the observation of non-CxSM final state signatures would add new information but not limit our findings with respect to the distinction of the models based on the signatures that we investigate here. Our analysis, provided the answer is positive, can therefore be seen as a trigger for further studies in the future taking into account other decay processes.
To organise the discussion we distinguish four cases, which are in principle different from the point of view of the experimental searches. The cases that are possible both in the CxSM-broken and in the NMSSM are: For the new scalar produced in association, in the low mass region the most important decay channel is the one into b-quarks followed by the decay into τ leptons. 14 Decays into photons might become interesting due to their clean signature. In the high mass regions, if kinematically allowed, the decays into pairs of massive vector bosons V ≡ W, Z and of top quarks become relevant. In the NMSSM the importance of the various decays depends on the value of tan β and the amount of the singlet component in the Higgs mass eigenstate, whereas in the singlet models it only depends on the singlet admixture to the doublet state. In order to simplify the discussion, we resort to the 4b final state. We explicitly checked other possible final states to be sure that they do not change the conclusion of our analysis in the following.

3) SM-like Higgs boson decaying into two light
In the NMSSM also h 3 → A 1 + A 1 is possible as well as In search channels where a non-SM-like Higgs boson is produced in the decay of a heavier non-SM-like Higgs boson we have a large variety of final state signatures. Leaving apart the decay channels not present in the CxSM-broken, the final state Higgs bosons can decay into massive gauge boson or top quark pairs, if they are heavy enough. Otherwise final states with b-quarks, τ leptons and photons become interesting. As we found that the various final states do not change our following conclusions, for simplicity we again focus on the 4b channel in the comparison with the NMSSM.
Note that, in the CxSM, we do not have the possibility of a heavier Higgs decaying into a pair of non-SM-like Higgs bosons that are not identical, so that this case does not appear in the above list. Figure 4 shows the rates corresponding for case 1). A heavier Higgs boson Φ is produced and decays into a pair of SM-like Higgs bosons h 125 that subsequently decay into a b-quark pair each. The red (blue) points represent all decays that are possible in the NMSSM (CxSM-broken), i.e. Φ = h 3,2 for h 125 ≡ h 1 and Φ = h 3 for h 125 ≡ h 2 . As can be inferred from the plot, such decay chains do not allow for a distinction of the two models. The maximum possible rates in the NMSSM can be as high as in the CxSM. The differences in the lowest possible rates for the two layers are due to the difference in the density of the samples. Such small rates, however, are not accessible experimentally. Note finally that here and in the following plots the inclusion of all possible Higgsto-Higgs decay channels in each of the models is indeed essential. Otherwise a fraction of points might be missed and a possible distinction of the models might be falsely mimicked.
The rates for case 2) are displayed in non-SM-like lighter Higgs boson ϕ, where the plot strictly refers to scenarios with m ϕ < m h 125 . Thus the blue points (CxSM) and red points (NMSSM) represent the cases Φ = h 3 and ϕ = h 1 with h 2 ≡ h 125 . Additionally, in the NMSSM, the green points have to be included to cover the case Φ = A 2 , ϕ = A 1 with h 125 ≡ h 1 or h 2 . The possible rates in the two models are shown as a function of m ϕ (left plot) and as a function of m Φ (right plot). The overall lower density of points in the NMSSM is due to its higher dimensional parameter space which, combined with its more involved structure, limits the computational speed in the generation of the samples. The performed scan starts at 30 GeV for the lightest Higgs boson mass, so that the Higgs-to-Higgs decays set in at m Φ = 155 GeV (blue points, right plot). There are no points above m ϕ = 121.5 GeV due to the imposed minimum mass difference of 3.5 GeV from the SM-like Higgs boson mass to avoid degenerate Higgs signals. The left plot shows that the masses of the A 1 , in the NMSSM, are mostly larger than about 60 GeV. While a more extensive scan of the NMSSM could yield additional possible scenarios, light Higgs masses allow for h 125 decays into Higgs pairs, that move the h 125 signal rates out of the allowed experimental range. This explains why there are less NMSSM points for very small masses m ϕ . The figure clearly demonstrates that the maximum achieved rates of the NMSSM in this decay chain can be enhanced by up to two orders of magnitude compared to the CxSM over the whole mass ranges of m ϕ and m Φ where they are possible. The observation of a much larger rate than expected in the CxSM-broken in the decay of a heavy Higgs boson into a SM-like Higgs and a lighter Higgs state would therefore be a hint to a different model, in this case the NMSSM.
The enhancement of these points can be traced back to larger values for the production cross sections of the heavy Higgs boson, for the Higgs-to-Higgs decay branching ratio and for the branching ratios of the lighter Higgs bosons into b-quark pairs in the NMSSM compared to the CxSM. In the NMSSM the larger production cross sections are on the one hand due to pseudoscalar Higgs broken CxSM vs NMSSM: m ϕ < m h125 broken CxSM vs NMSSM: m ϕ < m h125 production. In gluon fusion these yield larger cross sections than for scalars, provided that the top Yukawa couplings are similar, in scenarios where the top loops dominate. Additionally, the investigation of the top Yukawa coupling, which is the most important one for gluon fusion for not too large values of tan β, shows that in the NMSSM for the enhanced points it is close to the SM coupling or even somewhat larger. In comparison with these NMSSM top-Yukawa couplings, in the CxSM the top-Yukawa couplings are suppressed for all parameter points. This is because in the CxSM all Higgs couplings to SM particles can at most reach SM values and are typically smaller, for the new Higgs bosons, due to the sum rule i R 2 i1 = 1. This rule assigns most of the coupling to the observed Higgs boson and therefore its coupling factor can not deviate too much from 1. We also have NMSSM scenarios where the bottom Yukawa coupling is larger by one to two orders of magnitude in the NMSSM compared to the CxSM, This enhancement is due to tan β, for which we allow values between 1 and 30. In this case, where the top Yukawa coupling is suppressed, the Higgs bosons are dominantly produced in bb annihilation. The behaviour of the branching ratios can be best understood by first recalling that in the CxSM the branching ratios for each Higgs boson are equal to those of a SM Higgs boson with the same mass, in case no Higgs-to-Higgs decays are present. This is because all Higgs couplings to SM particles are modified by the same factor. If Higgs-to-Higgs decays are allowed the branching ratios even drop below the corresponding SM value. In the NMSSM, however, the branching ratio BR(ϕ → bb) is close to 1. Despite ϕ being the singlet-like A 1 or H 1 it dominantly decays into bb as the SUSY particles of our scan are too heavy to allow for ϕ to decay into them. As for the branching ratio BR(h 125 → bb), in the CxSM it can reach at most the SM value of around 15 0.59, while in NMSSM this branching ratio can be of up to about 0.7 due to enhanced couplings to the b-quarks. Note, in particular, that the uncertainty in the experimental value of the SM Higgs boson rates into bb in the NMSSM still allows for significant deviations of the SM-like Higgs coupling to bottom quarks from the SM value. Finally, the branching ratio BR(Φ → h 125 ϕ) for the enhanced points is larger in the NMSSM compared to the singlet case. This can be either due to the involved trilinear Higgs self-coupling or due to a larger phase space. The self-coupling λ Φh 125 ϕ depends on a combination of the NMSSM specific parameters λ, κ, A κ , A λ , v s , the Higgs mixing matrix elements and tan β. Through the Higgs mixing matrix elements it also depends on soft SUSY breaking masses and trilinear couplings, as we include higher order corrections in the NMSSM Higgs masses and mixing matrix elements. Therefore a distinct NMSSM parameter region or specific combination of parameters that is responsible for the enhanced self-couplings compared to the CxSM cannot be identified. The same holds for a possible different mass configuration of the involved Higgs bosons, with the higher-order corrected masses depending on NMSSM specific and SUSY breaking parameters. Figure 6 also refers to case 2), but this time the mass of the non-SM-like Higgs boson ϕ is larger than m h 125 . With this mass configuration the blue and red points of the CxSM and NMSSM, respectively, represent the cases Φ = h 3 and ϕ = h 2 with h 1 ≡ h 125 . In the NMSSM we also have the possibilities Φ = A 2 , ϕ = A 1 and h 125 ≡ h 1 or h 2 , which are covered by the green points. In the left plot the points set in at m ϕ = 128.5 GeV due to the imposed minimum mass distance from m h 125 . The upper limits on m ϕ in both models are due to the bounds on the input parameters chosen for the scans. In the right plot the points start at the kinematic lower limit of 253.5 GeV. Both plots clearly demonstrate, that the NMSSM Higgs-to-Higgs decay rates exceed those of the CxSM-broken in the whole mass range of m ϕ , respectively m Φ , by up to two orders of magnitude and allow for a distinction of the models if the largest possible signal rates in the NMSSM are discovered. The enhancement can be understood by looking at the production cross section of the heavy Higgs boson and at the various involved branching ratios. The enhanced production compared to the CxSM case is again most important in pseudoscalar NMSSM Higgs production, but also the H 3 production can be somewhat enhanced, because in the singlet case the heavy Higgs couplings to tops and bottoms are suppressed compared to the SM, while in the NMSSM we can have points where the top Yukawa coupling can be close to the SM value or even a bit larger and we can also have other points where the bottom Yukawa coupling can be much larger than the corresponding SM coupling while the top Yukawa coupling is suppressed. We have verified that there are scenarios where the top Yukawa coupling provides the dominant contribution to the cross section enhancement relative to the CxSM and other scenarios where it is the bottom Yukawa coupling. In the region where the NMSSM points yield larger rates, the branching ratio BR(Φ → h 125 ϕ) can be larger, but there are also cases, where the singlet case and the NMSSM lead to similar branching ratios. However, it turns out that the value of BR(ϕ → bb) in the NMSSM clearly exceeds that of the singlet case. For masses below the top pair threshold, it can reach values close to 1 in the NMSSM. In the singlet case it reaches at most the value of a SM Higgs boson with the same mass, or lower, when the decay into h 125 h 125 is kinematically possible. In the NMSSM below the top pair threshold the only important decay is the one into b-quark pairs, as the SUSY particles from our scan are too heavy and decays into massive gauge bosons are forbidden for ϕ = A 2 because it is CP-odd, and for ϕ = H 3 because of the coupling sum rules for the scalar Higgs couplings to gauge bosons. Above the top pair threshold ϕ can then also decay into tt.
For large values of tan β, however, the branching ratio into bb can reach values close to 1. In the cases where also decays into lighter SUSY particles, Higgs pairs or Higgs and gauge boson pairs broken CxSM vs NMSSM: m ϕ > m h125 are kinematically possible, the branching ratio into bb drops below 1, but is still larger than the corresponding branching ratio in the CxSM. As for the larger h 125 branching ratio into bb we have discussed this already in the context of Fig. 5.
The clear difference in the maximal values of the rates between CxSM-broken and NMSSM is less pronounced in the 2b2W final state, in particular for small masses m ϕ and m Φ . This is obvious, as here it is the decay A 2 → h 125 + A 1 that leads to larger rates (green points) than the CxSM counterpart. As the pseudoscalar A 1 , however, does not decay into massive gauge bosons, these rates go down in the 2b2W final state.
The investigation of case 3) (not plotted here), where a SM-like Higgs decays into two light Higgs states with subsequent decays into SM particles, shows that the maximum possible NMSSM rates are slightly larger than the maximum rates achieved in the CxSM-broken. This decay chain also takes place in RxSM-broken with maximum rates somewhat exceeding those of the NMSSM. The maximum possible rates are clearly limited by the observed rates for h 125 . An enhanced rate of the light ϕ into SM particles leaves some room to further increase the rate without violating the experimental measurements of the h 125 rates. In the NMSSM, however, the light scalar states h 1 and A 1 are mostly singlet-like with suppressed couplings to SM particles. The maximum possible rates in the NMSSM therefore cannot be expected to exceed those of RxSM-broken by far, so that the distinction of these two models (and hence also of the NMSSM and CxSM-broken) based on SM-like Higgs decays into lighter Higgs states is not possible.
In case 4) (also not plotted here) with a heavy Higgs boson decaying into a pair of non-SM-like Higgs bosons ϕ with m ϕ > m h 125 the maximum possible NMSSM rates exceed those of the CxSM in the 4b final state for certain kinematic configurations where m ϕ > ∼ 170 GeV. With maximum values of about 1 fb, they are, however, too small to be exploited. On the other hand the rates of the 4W/4Z final state yields much larger rates in both models. However, the maximum rates in the CxSM and NMSSM are for these final states comparable making impossible to distinguish them. In the opposite case, where m ϕ < m h 125 , in the dominant 4b final state the maximum rates are again comparable and no distinction of the models is possible.
To finalise the discussion we remark that our conclusions do not depend on the approximation of Higgs pair production from the decay of a resonantly produced heavy Higgs boson. In order to verify this we computed for all possible di-Higgs final states the complete Higgs pair production processes from gluon fusion which take into account both the resonant and the continuum diagrams. For the computation of these processes we implemented our models (i.e. the NMSSM and the complex singlet extension CxSM) into the existing code HPAIR 16 . This Fortran code, initially designed for the SM and the MSSM, calculates Higgs pair production at NLO QCD accuracy. We found that the difference between the resonant approximation and the full calculation, including the continuum, is typically of the order of ∼ 10% in the region of the plots with large rates where we can distinguish the CxSM from the NMSSM. However, it should be noted that in a full experimental analysis, the resonant production gives an extra handle in reducing the background relative to the continuum. This extra handle will surely make the signal significance in the resonant case much larger than the continuum one. In fact, although we are presenting here total cross sections, the signal events are mainly concentrated in the bins around the peak of the resonant scalar. This kinematical cut around the peak, even if broad, is only possible in the resonant case. Thus our conclusions are unaffected by the use of the approximation.

Benchmark Points for the LHC Run 2
In this section we provide a set of benchmark points with focus on the CxSM. At the end, we also discuss the RxSM which is simpler with respect to the possibilities for Higgs-to-Higgs decays.
As discussed in section 2, in these singlet models, at leading order in the electroweak corrections, the production cross sections and decay widths into SM particles are given by multiplying the corresponding SM result with the mixing matrix element squared, R 2 i1 , cf. Eq. (8). One can show that the signal strength factor is approximately given by There is hence only one common signal strength factor regardless of the decay channel into SM particles. This provides a useful measure for the deviation of the model point from the SM, and it is given by the squared overlap of the scalar state with the SM Higgs fluctuation (i.e. the mixing matrix element) multiplied by the sum of the branching ratios into SM particles. If there are no Higgs-to-Higgs decays the latter factor is one.
In the next subsections we present several benchmark points. These were chosen as to cover various physical situations. From a phenomenological perspective, we are interested in maximising the visibility of the new scalars in the LHC run 2 and in covering, simultaneously, many kinematically different possibilities. In particular we are interested in scenarios where µ C h 125 /µ T h 125 is as large as possible, so that we can observe the new scalars produced in chains containing the SM-like Higgs boson at the LHC run 2, while preserving consistency with the LHC run 1 measurements. Thus, in many of the points presented below we have tried to maximise the cross section for Higgs-to-Higgs decays. At the same time we required the rates of the SM-like Higgs to be within 2σ or at most 3σ from the global signal strength provided by the combination of the ATLAS and CMS data from the LHC run 1. We require at least 3σ consistency but in many cases we find points that satisfy the required properties within 2σ.
From the theoretical viewpoint, whenever possible, we choose points for which the model remains stable up to a large cutoff scale µ. Here µ is the cutoff scale at which the couplings of the theory reach a Landau pole or the scalar potential develops a runaway direction. This scale was obtained through a renormalisation group analysis as detailed in [70] and it will be indicated by the quantity log 10 (µ/GeV). In the dark matter phase of the CxSM we also require that the dark matter relic density predicted by the model, Ω A h 2 , is consistent with the combination of the measurements of the cosmic microwave background (CMB) from the WMAP and Planck satellites [79,80]. We choose points where Ω A h 2 is within 3σ of the central value of the combination, which is Ω c h 2 = 0.1199 ± 0.0027.

CxSM Broken Phase
In Tables 4 and 5 we show a sample of various kinematically allowed set-ups for the three mixing scalars of the broken phase of the CxSM. The first, Table 4, contains the parameters that define the chosen benchmark points and the production rates of the lightest and next-to-lightest Higgs bosons h 1 and h 2 in the various final states. The corresponding values for h 3 are listed in Table 5. For h 2 and h 3 the tables contain, in particular, the Higgs-to-Higgs decay rates. We also give the signal rates µ h i (i = 1, 2, 3) as defined in Eq. (13). We have chosen two points where the SM-like Higgs is the lightest Higgs boson (CxSM.B1 and CxSM.B2); two points where it is the next-to-lightest Higgs boson (CxSM.B3 and CxSM.B4); and one point where it is the heaviest (CxSM.B5). An interesting feature is that there are points for which the model remains stable up to a large cutoff scale µ, as shown by log 10 (µ/GeV) in the last row of Table 5. In particular the benchmark points CxSM.B3 and CxSM.B4 are theoretically very interesting, since the new heavy scalar here stabilises the theory up to a scale close to or above the GUT scale of 10 16 GeV.
Most points were chosen such that the cross sections for the indirect decay channels of the new scalars can compete with the direct decays. In particular, in most cases, we have tried to maximise h 3 → h 2 + h 1 where all three scalars could be observed at once. We have furthermore chosen points with large cross sections for the new scalars, so that they can also be detected directly in their decays. Since each point represents a different kinematic situation we discuss the main features of each of them separately: • CxSM.B1: For this point the SM-like Higgs is the lightest of the three Higgs bosons and all Higgs-to-Higgs decay channels are open apart from h 3 → h 2 + h 2 . 17 The presented point has been chosen such that it has a maximal ratio µ C h 125 /µ T h 125 for this kinematical situation within the imposed bounds. The production rates for Higgs-to-Higgs decays of the non-SM-like Higgs bosons at the LHC run 2 are also maximised. With the direct production of the heavier Higgs bosons in SM-like final states being on the lower side 18 , the additional Higgs-to-Higgs decays with rates of e.g. h 2,3 → h 1 + h 1 → bb + τ τ of a few fb suggest, that all new scalars can be expected to be observed. • CxSM.B3: For this point, where h 125 ≡ h 2 and m h 1 < m h 125 /2, all kinematic situations for the scalar decays are available while the spectrum remains light. As can be inferred from Fig. 1 (lower  state, so that Higgs-to-Higgs decays present an interesting discovery option for the heavy Higgs states. Furthermore, large production cross sections have been required for the new light scalar h 1 so that it will be visible in its direct decays in addition to chain production from heavier scalars.
• CxSM.B4: This scenario differs from the previous one in the larger h 1 mass so that the channel h 2 → h 1 + h 1 is kinematically closed. At the same time the direct h 1 production

CxSM Dark Phase
In the phase of the CxSM that contains dark matter, the first requirement is that the relic density predicted by the model agrees with the measurements from the CMB. In Table 6 we have selected four benchmark points that all obey this requirement.
For the first two points, CxSM.D1 and CxSM.D2, the lightest of the two visible scalars is the SM-like Higgs. In point CxSM.D1, both visible scalars can decay invisibly into the dark matter candidate A whereas in point CxSM.D2 h 1 cannot. Both benchmark points feature large invisible branching ratios, with the largest one reaching 29% for h 1 → A+A in CxSM.D1. Also the branching ratios for the Higgs-to-Higgs decays h 2 → h 1 + h 1 are large with 18.4% for CxSM.D1 and 28% for CxSM.D2. The signal rates for these decays in the bbτ τ final state are 13 fb and 33 fb, respectively, so that chain decays contribute to the discovery of h 2 , and also h 1 . Furthermore, both benchmark points have large cross sections for direct production of h 2 so that it can also be discovered in its direct decays into SM particles. Another attractive feature of these points is that the new heavy scalar h 2 can stabilise the theory up to a high scale close to the GUT scale, as can be inferred from log 10 (µ/GeV) in the last row of the table for CxSM.D1 and CxSM.D2.
In the scenarios CxSM.D3 and CxSM.D4, where the SM-like Higgs boson is the heaviest of the two visible Higgs bosons, the overall spectrum is lighter and the theory must have a UV completion above ∼ 10 3 TeV. Point CxSM.D3 represents a case with no invisible decays allowed, and in CxSM.D4 decays of the SM-like Higgs h 2 into a lighter Higgs pair are forbidden at the expense of allowing for a large invisible decay into the dark matter state A. In CxSM.D3 the light Higgs state h 1 can either be discovered directly or in the chain decay of the SM-like Higgs h 2 into an h 1 pair. The production rates of 9.0 pb in the 4b final state and 1.6 pb in the bbτ τ final state are rather large and complement the discovery of h 1 in direct production. At the same time the rather important branching ratio BR(h 2 → h 1 h 1 ) of 29%, which is responsible for these rates, drives the overall rate of the SM-like Higgs, µ h 2 , to a value at the edge of compatibility with the LHC data. The discovery of the non-SM-like light h 1 in CxSM.D4, that has a mass close to the SM-like Higgs, is extremely challenging. With no contribution from chain decays, it is only accessible in its decays into SM particles, with rates in the bb and τ τ final states that are comparable to those of the SM-like Higgs, but much smaller rates in the gauge boson final states. Table 7 contains four benchmark points for the two possible kinematic configurations in this model. The points RxSM.B1 and RxSM.B2 correspond to the case where the SM-like Higgs boson is the lightest of the two scalar states. Benchmark RxSM.B1 allows for the decay h 2 → h 1 + h 1 and we have chosen a point with a relatively large cross section σ 2 × BR(h 2 → h 1 + h 1 ) for such a chain decay. It is comparable to the direct h 2 production cross section. The bbτ τ final state in the former even reaches 72 fb. Thus h 2 could in principle be found directly or in its chain decays. Note that the fraction of chain decays in SM-like Higgs production amounts to 5%. For complementarity, in RxSM.B2 we selected a point where the decay into scalars is kinematically closed, but instead various direct decay channels of h 2 are enhanced compared to RxSM.B1, most notably the W W final state but also bb, τ τ and γγ.

RxSM Broken Phase
The benchmarks RxSM.B3 and RxSM.B4 feature a SM-like Higgs boson that is the heaviest of the two scalars. RxSM.B3 was again chosen such that the non-SM-like Higgs h 1 can be found directly or in the decay h 2 → h 1 + h 1 . In particular note the large rates for the direct h 1 production and decay into the bb and also τ τ final states and compare with the indirect processes h 2 → h 1 + h 1 → bb + bb and bb + τ τ , where the magnitude of the latter two is comparable to the former two. For RxSM.B4, on the contrary, we have chosen a point to represent the situation where the  indirect channel is closed. Here we have required a larger cross section for direct h 1 production allowing for larger rates into the bb and τ τ final states. Still the discovery of h 1 with a mass very close to the Z boson peak will be challenging.

Conclusions
We have analysed in detail the phenomenology of two Higgs bosons final states in a real (RxSM) and a complex (CxSM) singlet extensions of the SM. Both models contain phases with one or two new Higgs bosons, which may be found either directly or in Higgs-to-Higgs decays in the next runs of the LHC. We have performed a comparison of the achievable rates for Higgs-to-Higgs decays in the various singlet models and with respect to the NMSSM. Finally, we have presented benchmark points for the singlet models at the LHC run 2. We started by presenting the models and the phenomenological constraints that were imposed.
For this purpose a new code based on HDECAY was developed for the calculation of the decay widths and branching ratios of the scalar particles present in the RxSM and the CxSM, both in the symmetric and in the broken phase. In this publicly available tool, sHDECAY, the SM higher order electroweak corrections are consistently turned off. Including the state-of-the-art higher order QCD corrections, the code will allow a rigorous interpretation of the data in these singlet extensions of the SM.
In the numerical analysis, we first investigated how important the contribution from Higgs-to-Higgs decays to the signal of the 125 GeV Higgs boson can be compared to its direct production and decay. Depending on the singlet model and on the scenario, the indirect production of the SM-like Higgs boson, through the decay of a heavier scalar, can attain a maximum value ranging between 2%(4%) and 9%(17%) while remaining compatible with the Higgs signal measurements within 2σ(3σ).
Subsequently, we performed a systematic comparison of Higgs-to-Higgs decay rates at the LHC run 2 within the singlet models and with respect to the NMSSM. Among the singlet extensions, the rates achieved in the broken phase of the RxSM and the dark matter phase of the CxSM, in Higgs-to-Higgs decays with two identical scalars in the final states, are larger than those in the broken phase of the CxSM. In particular in the case where the non-SM-like scalar state is lighter than m h 125 /2, the maximum rates of the former two exceed those of the latter by up to two orders of magnitude and can reach about 10 pb in the 4b-quark final state.
The comparative analysis between the maximum allowed rates in the singlet extensions and those of the NMSSM addressed the question: to which extent Higgs-to-Higgs decays allow for a distinction of different Higgs sectors with more than two neutral Higgs bosons? The broken phase of the CxSM features such a Higgs spectrum. We found that a clear distinction from the NMSSM is possible based on Higgs-to-Higgs decay rates, but only in final states with two different scalars, i.e. in pp → h i → h j +h k , with i = j = k. The maximum rates obtained in the NMSSM significantly exceed those of CxSM-broken. This means that, if no direct hints of new physics outside the Higgs sector are found, two-Higgs final states with different masses can play the role of smoking gun signatures for the distinction of extended Higgs sectors based on pure (complex) singlet extensions and those of non-minimal SUSY sectors. Their analysis should therefore play a prime role at the LHC run 2.
Finally, we proposed benchmark points for the RxSM and CxSM at the LHC run 2. Our guiding principles in defining these benchmarks were: i) to maximize the rates for different Higgs-to-Higgs processes, so that the new scalars can be found in many different channels, ii) to look for points where the theory is stable up to a high scale and, iii) that in the dark phase of the CxSM the points conform with dark matter observables and null searches. We provide full information on the input parameters for each point together with all relevant cross sections at the LHC run 2, their compatibility with dark matter observables and the high-energy stability cutoff scale. The benchmarks cover different kinematic situations and vary in their prospect of discovering the related Higgs spectrum. They can therefore serve as guidelines for the experiments in tuning their experimental analyses.
Higgs self-couplings needed for the decay widths. For the computation of the decay widths in the singlet extensions, the parameter isinglet has to be set to 1. The model is chosen by setting the input parameter icxSM to 1 (RxSM-broken), 2 (RxSM-dark), 3 (CxSM-broken) or 4 (CxSM-dark).
The input values COUPVAR, HIGGS, SM4, FERMPHOB and ON-SH-WZ have to be set to zero. Even if the user does not set them to zero, by choosing isinglet to be 1 they are set internally to zero in shdecay.f. The input parameters of the models that can be chosen in shdecay.in are the same as the ones specified in section 2.1 for the CxSM and in 2.2 for the RxSM. Only v is not an input value but obtained internally from the Fermi constant.
All files necessary for the program can be downloaded at the url: http://www.itp.kit.edu/∼maggie/sHDECAY The tar file provided at the webpage contains, besides the main routine shdecay.f and the input file shdecay.in, several help files, stemming already from the original code HDECAY, and a makefile for compilation. The program is compiled with the file makefile by typing make, which provides an executable file called run. Typing run executes the program and all computed branching ratios for the various Higgs bosons are written out together with their masses and total widths. For the CxSM-broken the output files are called br.cbij with i=1,2,3 denoting one of the three Higgs bosons and j counting the three output files for each of the Higgs bosons. The first one contains the fermionic decays, the second the decays into gauge bosons and the total width, and the third one the Higgs-to-Higgs decays. The corresponding files for CxSM-dark are called br.cdkj (j=1,2,3) for the two visible Higgs bosons k=1,2. The output files for RxSM-broken are named br.rbkj with k=1,2 for the two Higgs states and j=1,2,3. In RxSM-dark we have the output files br.rd1j (j=1,2,3) for one visible Higgs state. Besides these files all input parameters are written out in a file called br.input. On the webpage, sample output files can be found for a given input. Updates will be made available there as well.

B Feynman Rules for the Triple Higgs Vertices
The widths for the decays h i → h j + h j and h i → h j + h k (k = j) are given by and where g ijk is the coupling between the scalars i, j, k defined in Eq. (9) and m j is the mass of the scalar state h j . 20

B.1 CxSM
The triple scalar vertices are shown below for the broken phase, assuming the normalisation of Eq. (9). Note that the g ijk are completely symmetric under interchange of indices. g 111 = 3 2 d 2 R 2 13 + R 2 12 (R 13 v A + R 12 v S ) + δ 2 R 11 v R 2 13 + R 2 12 + δ 2 R 2 11 (R 13 v A + R 12 v S ) + λR 3 11 Similarly we can define the quartic couplings Then for the broken phase one obtains: These expressions are all valid for the dark matter phase by replacing α 2 = α 3 = 0 in the mixing matrix elements and setting v A = 0.

B.2 RxSM
For this model, the relevant cubic couplings are The dark phase is obtained by setting α = 0 and v S = 0.