Two-real-scalar-singlet extension of the SM: LHC phenomenology and benchmark scenarios

We investigate the LHC phenomenology of a model where the Standard Model (SM) scalar sector is extended by two real scalar singlets. A Z 2 ⊗ Z 2 (cid:48) discrete symmetry is imposed to reduce the number of scalar potential parameters, which is spontaneously broken by the vacuum expectation values of the singlet ﬁelds. As a result, all three neutral scalar ﬁelds mix, leading to three neutral CP-even scalar bosons, out of which one is identiﬁed with the observed Higgs boson at 125 GeV. We explore all relevant collider signatures of the three scalars in this model. Besides the single production of a scalar boson decaying directly to SM particle ﬁnal states, we extensively discuss the possibility of resonant multi-scalar production. The latter includes decays of the produced scalar boson to two identical scalars (“symmetric decays”), as well as to two diﬀerent scalars (“asymmetric decays”). Furthermore, we discuss the possibility of successive decays to the lightest scalar states (“cascade decays”), which lead to experimentally spectacular three- and four-Higgs ﬁnal states. We provide six benchmark scenarios for detailed experimental studies of these Higgs-to-Higgs decay signatures.


I. INTRODUCTION
The Large Hadron Collider (LHC) at CERN is the first experimental facility that directly probes the mechanism of electroweak symmetry breaking (EWSB), described in the Standard Model of particle physics (SM) by the Brout-Englert-Higgs mechanism [1][2][3][4][5][6]. The milestone discovery of a Higgs boson with a mass of ∼ 125 GeV in 2012 [7,8] and the ongoing measurements of its properties at the LHC (see e.g. Ref. [9]) open the door to a deeper understanding of the structure of EWSB. Indeed, this experimental endeavor may reveal first signs of new physics beyond the SM (BSM), as many well-motivated BSM extensions affect the phenomenology of the observed scalar particle. However, by the end of Run-II of the LHC, with the full collected data of ∼ 150 fb −1 per experiment (ATLAS and CMS) still being analyzed, Higgs signal rate measurements in various production and decay channels [10][11][12][13][14][15][16][17][18][19][20] are so far in very good agreement with the SM predictions.
Extensions of the SM by scalar singlets are among the simplest possible model beyond the SM (BSM). The most general extension of the SM by n real scalar singlet fields φ i (i ∈ [1, . . . , n]) has a scalar potential of the form where with real coefficient tensors. Here, Φ describes the scalar SU (2) L doublet field of the SM and V SM denotes the scalar potential of the SM. An extension by complex singlets can always be brought into this form by expanding fields and coefficients into real and imaginary parts. Since the φ i are pure gauge singlets they have trivial kinetic terms that do not induce any gauge interactions, leading to the following contributions to the electroweak (EW) Lagrangian: In addition, it is not possible to write down gauge invariant and renormalizable interactions between a scalar singlet and any of the SM fermions. The singlets will therefore only interact with the SM Higgs boson through the couplings of the scalar potential and, if they acquire a non-zero vacuum expectation value (vev), mix with the SM Higgs boson and thereby inherit some of its gauge and Yukawa couplings. This is also the reason why -as long as no new interactions of the scalar singlet fields with additional particles occur -there is no physical difference between a parametrization in terms of N complex scalar singlet fields or 2N real scalar singlet fields. Naively, one would expect that imaginary parts of complex scalar fields are CP-odd, and mixing them with the real parts or the SM Higgs boson would lead to CP-violation. However, due to the singlets not having any gauge or fermion couplings it is always possible to find a CP-transformation under which all of them are CP-even [21,22]. Thus any pure singlet extension of the SM is a theory of only CP-even scalars.
Experimentally, singlet extensions can be explored in two complementary ways at the LHC. First, precise measurements of the 125 GeV Higgs signal rates probe the structure of the doublet-singlet mixing, as well as possible new decay modes of the observed Higgs boson to new light scalar states. Second, direct searches for new scalars may reveal the existence of the mostly singlet-like Higgs bosons. For the latter the discovery prospects depend on the singlet-doublet mixing and the new scalar's mass (both governing the production rates), and on the decay pattern of the produced scalar state. In general, decays directly to SM particle final states as well as to two lighter scalar states ("Higgs-to-Higgs decays") are possible. While some of the former decays are already searched for by the LHC experiments, current searches for Higgs-to-Higgs decays focus almost exclusively on the signatures h S → h 125 h 125 (where h S denotes the new Higgs state with mass above 250 GeV) [58][59][60][61][62][63][64][65][66][67][68][69], or h 125 → h S h S (with the h S mass below 62.5 GeV) [70][71][72][73][74][75][76]. The model considered in the following, however, features also Higgs decays to unidentical scalar bosons ("asymmetric decays"), Higgs decays involving only non-SM-like Higgs bosons, as well as the possibility of successive Higgs-to-Higgs cascade decays. All of these decay signatures have not been experimentally explored in detail yet. 1 We will extensively discuss them in this paper and show that they lead to novel collider signatures with sizable signal rates that are experimentally interesting for the analysis of Run-II data as well as the upcoming LHC runs. We provide six dedicated two-dimensional benchmark scenarios, each highlighting a different Higgs-to-Higgs decay signature that has not been probed experimentally. We strongly encourage the experimental collaborations to investigate these novel signatures using current and future collider data. This paper is structured as follows. We introduce the model in Section II and summarize all relevant theoretical and experimental constraints on the parameter space in Section III. In Section IV we discuss the collider signatures of the model and present the impact of current LHC searches on the parameter space. In Section V we propose six benchmark scenarios for LHC searches for Higgs-to-Higgs decay signatures. We conclude in Section VI.

A. Scalar potential and model parameters
The two real singlet model (TRSM) adds two real singlet degrees of freedom to the SM. These are written as two real singlet fields S and X. In order to reduce the number of free parameters two discrete Z 2 symmetries are introduced. The most general renormalizable scalar potential invariant under the Z 2 S ⊗ Z 2 X symmetry is All coefficients in Eq. (5) are real, thus the scalar potential contains nine model parameters in total. We provide a translation of these coefficients to the scalar potential parameters in the complex scalar singlet parametrization in Appendix A.
Depending on the vevs acquired by the scalars different phases of the model can be realized. We decompose the fields (in unitary gauge) as leading to the tadpole equations These have solutions for any values of v, v S , v X . However, to achieve electroweak symmetry breaking v = v SM ≈ 246 GeV is required. If v S , v X = 0 the Z 2 symmetries are spontaneously broken, and the fields φ h,S,X mix into three physical scalar states. This is called the broken phase.
If v X = 0 the field φ X does not mix with φ h,S , does not acquire any couplings to SM particles, and is stabilized by the Z 2 X symmetry. 2 This makes it a candidate particle for dark matter (DM). The phenomenology of the two visible scalar states is very similar to the 2 The case of v S = 0 is equivalent by renaming S ←→ X.
real singlet extension [28][29][30][31][32]. If both singlet vevs vanish, φ h is the SM Higgs boson, and the two singlets both form separate dark sectors stabilized by their respective Z 2 symmetries. In this case collider phenomenology is (at tree-level) only impacted by possible invisible decays of h 125 to the DM particles.
In this work, we focus on the broken phase as it leads to the most interesting collider phenomenology. The mass eigenstates h 1,2,3 are related to the fields φ h,S,X through the 3 × 3 orthogonal mixing matrix R  We assume the mass eigenstates to be ordered by their masses and parametrize the mixing matrix R by three mixing angles θ hS , θ hX , θ SX . Using the short-hand notation it is given by Where the angles θ can be chosen to lie in without loss of generality. In the TRSM it is possible to express the nine parameters of the scalar potential through the three physical Higgs masses, the three mixing angles, and the three vevs. These relations are given by where a sum over i is implied. Fixing one of the Higgs masses to the mass of the observed Higgs boson, M a 125 GeV, and fixing the Higgs doublet vev to its SM value, v 246 GeV, leaves seven free input parameters of the TRSM: This is an important practical difference between the TRSM and another well-studied extension of the SM with two real singlet degrees of freedom, the CxSM [35,38]. The CxSM is expressed in terms of a complex singlet with a softly broken U (1) symmetry of the singlet phase imposed on the scalar potential. This more stringent symmetry leaves only seven model parameters such that one of the physical scalar masses and one of the singlet vevs are dependent parameters. In contrast, the TRSM is consistent for any combination of masses, mixing angles, and vevs, and therefore allows to cover the full possible kinematic phase space of Higgs-to-Higgs decay signatures, as we will do when defining the benchmark scenarios in Section V.
As in all pure singlet extensions the couplings of the scalar boson h a (a ∈ {1, 2, 3}) to all SM particles are given by the SM prediction rescaled by a common factor where R a1 denotes the doublet admixture of the mass eigenstate h a . Due to the orthogonality of the mixing matrix these obey the important sum rule

B. Collider Phenomenology
The triple Higgs couplings are of special phenomenological interest in the TRSM. Using Eq. (15) they can be expressed directly through the input parameters of Eq. (16). The Similarly, the coupling of three different scalars is given by and the triple Higgs self-couplingλ aaa reads With these definitions the tree-level partial decay width of a scalar h a into two scalars h b and h c (where b = c is allowed) is then given by with With this information, the phenomenology of a TRSM Higgs boson h a can be fully obtained from the predictions for a SM-like Higgs boson h SM of the same mass. Throughout this work we employ the narrow width approximation to factorize production cross sections and branching ratios (BRs).
For a certain production process (e.g. gluon gluon fusion) the cross section, σ, for h a with mass M a can be obtained from the corresponding SM Higgs production cross section, σ SM , by simply rescaling Since κ a rescales all Higgs couplings to SM particles, Eq. (24) is exact up to genuine EW corrections involving Higgs self-interactions. In particular, this holds to all orders in QCD. The scaling factor κ a also rescales universally the partial widths of h a decays into SM particles, which in turn leads to a rescaling of the SM total width as where Γ(h a → SM; M a ) denotes the sum of all partial widths of h a into SM particle final states. Note that this alone can never change the BR predictions of h a into SM particles.
Using the results of Eq. (22) we can obtain the BRs of h a decays to other scalar bosons, Denoting the sum of these "new physics" (NP) decay rates to scalar boson final states as the BRs of h a decays into any final state F SM composed entirely of SM fermions and gauge bosons are given by One important special case is that in the absence of BSM decay modes -which is always the case for the lightest Higgs bosons h 1 -h a has BRs identical to a SM-like Higgs boson of the same mass. Figure 1 shows the decay branching ratios of a SM-like Higgs boson h SM as a function of its mass. As long as BR(h a → NP) = 0, i.e. if no Higgs-to-Higgs decays are possible for h a , the scalar boson h a has exactly the BRs shown in Fig. 1. The numerical values are taken from Ref. [77], based on state-of-the-art evaluations using HDECAY [78][79][80] and Profecy4F [81][82][83].

III. SETUP OF THE PARAMETER SCAN
In order to assess the phenomenologically viable regions of the parameter space we apply all relevant theoretical and experimental constraints, which are discussed in the following. We furthermore provide details of our numerical setup.

A. Theoretical Constraints
Unitarity constraints provide important upper bounds on the multi-scalar couplings and the scalar masses. In the TRSM we have derived perturbative unitarity bounds in the high energy limit by requiring the eigenvalues M i of the 2-to-2 scalar scattering matrix M to fulfill The resulting constraints on the parameters of the scalar potential are where a 1,2,3 are the three real roots of the cubic polynomial Closed form conditions for boundedness of the scalar potential, Eq. (5), have been derived in [84,85]. In our notation they read It has been proven in Ref. [39] that at tree-level a vacuum of the form of Eq. (6) is always the global minimum of the scalar TRSM potential. Therefore no additional constraints from vacuum decay need to be considered.

B. Experimental Constraints
We use the oblique parameters S, T and U [86][87][88][89] to parametrize constraints from electroweak precision measurements, which are compared to the latest fit results [90]. The results of Refs. [91,92] are applicable to the TRSM to obtain model predictions for S, T and U . 3 Flavor constraints are not relevant as the singlets do not change the Yukawa sector. We use HiggsBounds-5.4.0 [98][99][100][101][102][103] to check for agreement with the bounds from searches for additional Higgs bosons.
Important constraints on the model parameter space arise from the signal rate measurements of the observed 125 GeV Higgs boson, which we denote by h 125 in the following. These constraints are especially relevant in singlet extensions as there are effectively only two BSM parameters that enter the phenomenology of the h 125 : its coupling scale factor κ 125 and its decay rate BR(h 125 → NP) into new particles (see Section II B).
We use HiggsSignals-2.3.0 [104][105][106][107] to test for agreement with the observations at the 2σ level using a profiled likelihood ratio test with the SM as alternative hypothesis. In practice, the likelihood ratio test statistic is calculated via the difference between the loglikelihoods, which in turn is approximated as ∆χ 2 = χ 2 − χ 2 SM within HiggsSignals. As we have two relevant statistical degrees of freedom that can influence the Higgs signal rate predictions (see above), we obtain the 2σ confidence region by demanding ∆χ 2 ≤ 6.18. [10][11][12][13][14][15][16][17]19, 20] with up to ∼ 137 fb −1 of data collected during Run-II at a center-of-mass energy of 13 TeV, as well as the ATLAS and CMS combined measurements from Run-I [9].

HiggsSignals-2.3.0 contains the latest measurements from ATLAS [108] and CMS
A further complication may arise in this model in case that two or even all three scalar bosons have a mass around 125 GeV. HiggsSignals then automatically takes into account a possible superposition of their signals in the test against the Higgs rate measurements by incoherently summing the contributions of all scalars. This approach neglects any possible interference effects, see Ref. [105] for details. A similar approach is employed in HiggsBounds to combine multiple scalars that lie within the experimental mass resolution.
Assuming only one scalar boson is responsible for the observed signal at 125 GeV, we show the constraints from Higgs signal rate measurements in the (simplified) two-dimensional parameter plane (κ 125 , BR(h 125 → NP)) in Fig. 2. 4 If no BSM decay modes of h 125 exist, a lower bound on κ 125 > 0.963 at 95 % C.L. is obtained. For the other limiting case of exactly SM-like couplings, κ 125 = 1, we find a limit of BR(h 125 → NP) < 7.3 %. The 2σ limit between these two limiting cases follows approximately a linear slope. The region κ > 1 is only included for completeness in Fig. 2 but cannot be realized in the TRSM, see Eq. (18). Note that this analysis is applicable to any model where a singlet scalar mixes with the Higgs boson. This bound can e.g. be applied to Higgs portal models, where it gives a stronger constraint than direct measurements of BR(h 125 → invisible) [110,111].

C. Numerical Evaluation
Based on these constraints we performed a large scan of the TRSM parameter space using an updated private version of the code ScannerS [35,38,112,113]. For the determination of viable regions in the parameter space, we apply all of the constraints described above. Note that bounds from signal strength measurements are evaluated with HiggsSignals for each point individually. This guarantees that the possibility that two or even all three Higgs bosons may have masses close to 125 GeV and therefore contribute to the observed signal is correctly accounted for.
We parametrize the model via the input parameters given in Eq. (16). For the numerical results presented in Section IV we independently sample from uniform distributions for each parameter. We allow for the non-h 125 Higgs masses and the singlet vevs to lie within and vary the mixing angles throughout their allowed range, Eq. (14). We only keep parameter points that pass all constraints. For the benchmark scenarios in Section V, we fix all parameters apart from the non-SM scalar masses, and scan the two-dimensional parameter space in a grid within the defined parameter ranges. The singlet vevs are only mildly constrained by current experimental results while theoretical constraints -in particular perturbative unitarity -only require them to not be substantially smaller than the scalar masses. On the other hand, as they enter the triple scalar couplings, Eqs. (19) and (20), they can influence the relative importance of different Higgs-to-Higgs decay modes without changing the remaining phenomenology. We therefore expect that future results from LHC searches for Higgs-to-Higgs decays will be able to constrain the vacuum structure of the singlet fields.
For the SM Higgs production cross sections and decay rates we use the predictions from Refs. [77,114]. The h SM production cross sections and total width are rescaled according to Eqs. (24) and (25) and combined with leading-order decay widths for the Higgs-to-Higgs decays from Eq. (22). For the h 125 production rates, we use the results of the N 3 LO calculation in the gluon gluon fusion (ggF) channel [115]. This calculation uses an effective description of the top-induced contributions. For scalar bosons with masses M a = 125 GeV we instead employ results from the NNLO+NNLL calculation [114] that accounts for top-quark mass induced effects up to NLO. Indeed, we find that the predictions of these calculations differ sizably for scalar boson masses M a 2m t , for instance, at M a = 400 GeV, In the following discussion of collider signatures we assume the production of a single scalar state via the dominant ggF process. In some cases, though, it might be worthwhile to investigate the discovery potential of the subdominant Higgs production processes of vector boson fusion or Higgs-Strahlung, pp → V φ (V = W ± , Z), as these give additional trigger options and may help to reduce the background. We leave the detailed exploration of the prospects for various production modes to future studies. Scalar pair production can proceed through a top-quark box diagram in addition to single Higgs production followed by a Higgs-to-Higgs decay. For pair production of h 125 these diagrams and their interference effects with the resonant production are known to be important (see e.g. Refs. [116][117][118][119]). For cases other than h 125 -pair production the box diagrams are significantly less important as they are always suppressed by the small κ factors of the non-h 125 scalars. Signal-signal interference effects between different resonant scalars of similar mass have also been shown to significantly impact di-Higgs production cross sections [119]. However, such mass configurations play no important role for most of the scenarios in the following discussion.

A. Signatures of New Scalars decaying to SM particles
The additional scalar bosons h a (h a = h 125 ) can decay directly to SM particles. The branching ratios of the various SM particle final states (F SM ) are obtained according to Eq. (28), and their relative rates (i.e. the ratios of branching ratios for different F SM decay modes) are identical to the corresponding SM predictions for a Higgs boson with mass M a . The rate for h a signal processes leading to F SM , normalized to the corresponding SM prediction, can therefore be expressed as This quantity is shown in Fig. 3 as a function of M a in the low mass region (left panel ) and high mass region (right panel ) for the sampled parameter points that pass all relevant constraints (see Section III). For M a roughly between 12 and 85 GeV LEP searches for e + e − → h a Z → bbZ [120] lead to an upper limit on the possible signal rate, as shown by the red lines in Fig. 3 (left). At larger mass values 190 GeV, the upper limit originates from LHC searches for pp → h a → W + W − and ZZ. The latest ATLAS [121] and CMS [122] limits are overlaid as green and orange lines, respectively, in Fig. 3 (right). For very large mass values 700 GeV direct LHC searches are not yet sensitive to probe the parameter space. In addition, we include in Fig. 3

B. Signatures of Resonant Scalar Pair Production
The model allows for resonant scalar pair-production at the LHC, or, in other words, the direct production of a single scalar h a followed by the "symmetric" or "asymmetric" decay into identical or different scalar states, respectively. Specifically, where, in the symmetric case, Eq. (38), a = 2, b = 1 or a = 3, b ∈ {1, 2}, and X denotes not further defined objects that may be produced in association with the scalar state (e.g., jets, vector bosons, etc.). The h 125 can be either of the three scalar states h a (a ∈ {1, 2, 3}). Processes of the symmetric type, Eq. (38), leading to pair production of h 125 are already being investigated, see e.g. Refs. [58][59][60][61][62][63][64][65][66][67][68][69] for recent LHC Run-II searches. Figure 4 (left) shows the 13 TeV LHC signal rate for the resonant scalar pair production process pp → h a → h 125 h 125 (a ∈ {2, 3}) as a function of the h a mass, M a . Figure 4 contains the complete sample of allowed parameter points generated according to Section III C. Overlaid are the most recent experimental limits on this process from the ATLAS [69] and CMS [66] collaborations. Figure 4 (left) illustrates that experimental searches in this channel are beginning to directly constrain the TRSM for resonance masses between around 380 GeV and 550 GeV. In contrast, LHC searches [70][71][72][73][74][75][76] for the inverted signature of single-production of h 125 which then decays into a pair of light h a (a ∈ {1, 2}) are not yet sensitive, as the indirect constraints from Higgs signal rates on the possible new decay modes, BR(h 125 → NP) ≤ 7.3% (see Fig. 2), are much stronger than the direct limits from these searches. 5 Both of these processes are under active experimental investigation and we expect the bounds to improve in the future.
We will now turn to the more exotic signatures resulting from Eqs. (38) and (39) that are not yet under active investigation. Following the processes in Eqs. (38) and (39), the two produced scalar states may further decay directly to SM particles. Alternatively, an h 2 final state may decay into the two lightest scalars: h 2 → h 1 h 1 . This can lead to interesting decay cascades leading to three or four scalar states that eventually decay to SM particles. The possible decay patterns within our model are depicted in a generic form in Fig. 5. Here and in the following we denote final states from Higgs decays composed of SM particles (i.e. gauge bosons or fermions) generically F SM , unless otherwise specified. For the more complicated final states we will use F n SM to denote an n-particle SM final state, where we count the SM particles before their decay (i.e., W ± , Z, and t are counted as one particle). As discussed above, the relative fractions of their decay rates solely depend on the mass of the decaying Higgs state.
We find that all possible Higgs-to-Higgs decay signatures, Eqs. (38) and (39), can appear at sizable rates in the allowed TRSM parameter space. In the next section we therefore design six two-dimensional benchmark scenarios that highlight these signatures in detail, and are tailored to initiate dedicated experimental studies and facilitate the design of corresponding searches. As a final remark, we briefly want to comment on the possible size of the total width of the resonantly-produced scalar state h a . Figure 4 (right) shows the ratio of the total width over the mass, Γ a /M a , in dependence of M a and the sum of h a decays to scalar states, BR(h a → NP). Parameter points with larger values of Γ a /M a overlay parameter points with smaller values. We can clearly see that parameter points with larger Γ a /M a tend to have sizable decay rates to scalar states. However, overall, Γ a /M a never exceeds values greater than around 18% in the considered mass range up to 1 TeV. In the discussion of the benchmark scenarios below we will comment on cases where Γ a /M a 1 %.

V. BENCHMARK SCENARIOS
In this section we define six benchmark scenarios in order to motivate and enable dedicated experimental studies of Higgs-to-Higgs decay signatures. Each scenario focusses on one (or more) novel signatures and features a (close-to) maximal signal yield that can be expected within the model while obeying the constraints described in Section III. The bench-benchmark scenario h 125 candidate target signature possible successive decays  Fig. 5).
The model parameters for all scenarios as well as the coupling scale factors κ a are given in Table II. All cross section values given in the following refer to production of the initial scalar through ggF at the 13 TeV LHC.
We employ a factorized approach relying on the narrow width approximation. For each benchmark scenario we will show both the BR(h a → h b h c ) (a, b, c ∈ {1, 2, 3}, a = b, c) and the cross section In all scenarios where either b = c or h b,c ≡ h 125 there is only one unknown BSM mass in the final state h b h c . In this case we will employ a further factorization where we present the BR(h b h c → F SM ) as a function of the remaining mass parameter. In this case the full cross section into a given SM final state can be obtained by where potential cascades, Fig. 5 parameters that may lead to different cross sections, branching ratios, or regions excluded by some constraints. As such, the regions of parameter space that are excluded by some constraint in a benchmark scenario should under no circumstances discourage experimental searches in this parameter region.
In the first benchmark scenario, BP1, we identify the heaviest scalar state h 3 with h 125 , and focus on the asymmetric decay h 125 → h 1 h 2 . The parameter values (see Table II) are chosen such that the couplings of h 3 to SM particles are nearly identical to the SM predictions, κ 3 1. At the same time, the parameter choice maximizes -within the experimentally allowed range -the branching ratio BR(h 125 → h 1 h 2 ), which is shown in Fig. 6 (top left) as a function of M 1 and M 2 . In Fig. 6 (top right) we show the corresponding signal rate for inclusive production via gluon gluon fusion. We find that the BR for h 3 → h 1 h 2 reaches up to 7 − 8 % which translates into a signal rate of around 3 pb. These maximal branching ratios are reached in the intermediate mass range for h 2 , M 2 ∼ 60 − 80 GeV. This feature is caused by the fact that the triple Higgs couplings are proportional to the masses (see Eq. (20)). Therefore, although phase space opens up significantly for light decay products, the branching ratios become smaller for M 2 < 40 GeV. In the hatched region in Fig. 6 the decay rate slightly exceeds the 2σ upper limit inferred from the LHC Higgs rate measurements (using HiggsSignals). We stress again that this excluded area is dependent on our parameter choices and strongly encourage experimental searches to cover the whole mass range.
Due to the sum rule, Eq. (18), the coupling scale factors κ 1,2 have to be very close to zero in order to achieve κ 3 ∼ 1. This means that the couplings of h 1 and h 2 to SM particles are strongly suppressed. As a result, if the decay channel h 2 → h 1 h 1 is kinematically open, M 2 > 2M 1 , it is the dominant decay mode leading to a significant rate for the h 1 h 1 h 1 final state. In BP1 we find that BR(h 2 → h 1 h 1 ) 100 % in this kinematic regime (i.e. above the red line in Fig. 6) with a very sharp transition at the threshold. If in addition M 1 10 GeV the h 1 decays dominantly into bb leading to a sizeable rate for the bbbbbb final state as shown in Fig. 6 (bottom right). If the h 2 → h 1 h 1 decay is kinematically closed, M 2 < 2M 1 , both scalars h 1 and h 2 decay directly to SM particles, with BRs identical to a SM-like Higgs boson with the corresponding mass (see Fig. 1). Therefore, for masses M 1 , M 2 10 GeV, the bbbb final state dominates, as shown in Fig. 6 (bottom left), while at smaller masses, combinations with τ -leptons and eventually final states containing charm quarks and muons become relevant.
In the second benchmark scenario, BP2, we identify h 125 ≡ h 2 and consider the production of h 3 followed by the asymmetric decay h 3 → h 1 h 125 . The scenario is defined in the (M 1 , M 3 ) parameter plane, and the remaining parameters are fixed to the values given in Table II. The mixing angles are chosen such that the production rate of h 3 is maximized, while the h 2 properties remain consistent with the measured Higgs signal rates. This results in a h 3 production rate of roughly 4% of the production cross section for a h SM at the same mass.
The phenomenology of BP2 is illustrated by Fig. 7. The BR(h 3 → h 1 h 2 ) shown in Fig. 7 (left) mostly stays above 20 % for M 3 350 GeV, reaching maximal values of around 50 − 55 % in the low mass region, M 3 ∼ 150 − 170 GeV. In this region, the corresponding signal rate in Fig. 7 (right) is about 0.6 pb. It remains above 50 fb as long as M 3 450 GeV. The shaded region in Fig. 7 is excluded by boundedness of the scalar potential. Again, this constraint depends strongly on the values of the model parameters and should not discourage experimental efforts to perform model-independent searches in this mass range. The total width of h 3 can reach maximal values of Γ 3 /M 3 ∼ 1.1 % in this benchmark scenario for M 3 480 GeV. The branching ratios for decays to SM final states originating from the h 1 h 125 two-particle state are shown in Fig. 8 for BP2 as a function of M 1 . In most of the mass range, the bbbb final state dominates, followed by bbW + W − and bbτ + τ − final states.
The cascade decay h 125 ≡ h 2 → h 1 h 1 is in principle possible if kinematically allowed and in compliance with the observed h 125 properties. However, we chose κ 2 2 small in order to maximize κ 3 within the experimental constraints. From Fig. 2 we see that, at the corresponding value of κ 2 , BR(h 125 → h 1 h 1 ) must not exceed ∼ 2.5 %. In BP2 this decay rate is always below 0.1 %.
Besides the asymmetric decay h 3 → h 1 h 2 the symmetric decays h 3 → h 1 h 1 and h 3 → h 2 h 2 are also present in this scenario. The decay h 3 → h 1 h 1 has a rate 25 % in the mass range M 3 250 GeV. The decay mode h 3 → h 2 h 2 only becomes kinematically open for M 3 2M 2 = 250 GeV, and reaches rates up to ∼ 28 %. Although these rates are not negligible in BP2, we shall define dedicated benchmark scenarios BP5 and BP6 below where these decay modes clearly dominate.
In benchmark scenario BP3 we identify h 125 ≡ h 1 and consider the production of h 3 followed by the asymmetric decay h 3 → h 125 h 2 . Similar to the BP2 scenario the mixing angles are chosen to maximize κ 2 3 5.7 % and BR(h 3 → h 1 h 2 ). The benchmark plane corresponding to the parameters given in Table II is shown in Fig. 9.
The BR(h 3 → h 125 h 2 ) shown in Fig. 9  If M 2 < 250 GeV BSM decay modes of h 2 are prohibited and its decay rates are identical to an h SM of the same mass (see Fig. 1). In this region the h 125 h 2 state dominantly decays into final states involving b-quarks and heavy gauge bosons as shown in Fig. 10. As soon as M 2 > 250 GeV the decay h 2 → h 125 h 125 becomes dominant, quickly reaching a rate of ∼ 70 %. Above threshold this rate remains largely independent of M 2 . The decay BRs of the resulting h 125 h 125 h 125 state to the most important six particle SM final states, F 6 SM , are given in Table III. The first row lists the direct branching ratios of h 125 h 125 h 125 while the second row includes the factor BR(h 2 → h 125 h 125 ) ≈ 68 %, which is an approximation obtained in the mass region 260 GeV < M 2 < 500 GeV. The resulting values can thus be compared directly to the BRs of the four particle F 4 SM in Fig. 10. For instance, we find that rates for bbbbbb, bbbbW + W − and bbW + W − final states are of comparable size for M 2 270 GeV.
The maximal production rates for the SM signatures amount to around 0.3 pb and 0.14 pb, respectively, where the latter is found when both decays are just above threshold, M 3 380 GeV and M 2 255 GeV.
In BP3, the competing symmetric decay h 3 → h 2 h 2 reaches rates of 20 % if kinematically allowed. Otherwise the decay h 3 → h 125 h 125 reaches similar values (and becomes dominant in the threshold region, We now turn to the symmetric Higgs-to-Higgs decay signatures. In benchmark scenario BP4 we identify h 125 ≡ h 3 and focus on the production of h 2 followed by its decay h 2 → h 1 h 1 . In order to avoid constraints from the Higgs rate measurements on the possible decays h 125 → h a h b (a, b ∈ {1, 2}), the relevant couplings must be tuned to rather small values while keeping |κ 2 | relatively large to ensure sizeable direct production of h 2 . The parameter choices for BP4 are listed in Table II. Fig. 11 shows the collider phenomenology of BP4. The branching ratio BR(h 2 → h 1 h 1 ) is larger than 50 % throughout the allowed parameter plane, as shown in Fig. 11 (left). For M 2 40 GeV it is by far the dominant decay mode of h 2 with a BR of more than 90 %. As the produced scalar boson is light, the signal rates shown in Fig. 11 (right) are enhanced by the large ggF cross section for light scalars. Even though h 2 is only produced with a rate of about κ 2 2 ∼ 5% of the SM Higgs cross section at the same mass, we still obtain signal rates of O(100 pb) in the low mass region M 2 20 GeV. However, this parameter region is partly constrained by LEP searches for e + e − → Zh 2 → Z(bb) [120]. For M 2 ≥ 20 GeV, where this limit is no longer sensitive, the signal rate can still reach 60 pb. Still, the signature remains experimentally challenging as the decay products for these low M 1 will be very soft.
The BRs for the decay modes of the h 1 h 1 state into SM particles are shown in Fig. 12. For M 1 10 GeV the decay into bbbb is dominant, followed by bbτ + τ − . For even lighter M 1 the predominant decay is into charm quarks. The h 125 h 1 h 1 coupling is very small in this scenario. Still -due to the large κ 125the process pp → h 125 → h 1 h 1 can enhance the total h 1 h 1 production cross section by up to ∼ 15% for large M 2 ∼ 125 GeV. On the other hand, interference effects between the resonant h 2 and h 125 contributions -similar to those discussed in Ref. [119] -remain negligible.
In the benchmark plane BP5 we identify h 125 ≡ h 2 and consider the production of the heavier scalar h 3 followed by its symmetric decay to the lightest scalar, h 3 → h 1 h 1 . In our parameter scan of the TRSM (see Section III C) we found that parameter points exhibiting a sizeable pp → h 3 → h 1 h 1 rate also tend to be strongly constrained by the Higgs signal strength measurements if 2M 1 < 125 GeV. In addition, if kinematically accessible, the decay modes h 3 → h 2 h 2 and/or h 3 → h 1 h 2 tend to dominate over the decay h 3 → h 1 h 1 . In order to define a suitable benchmark scenario for the pp → h 3 → h 1 h 1 process it is therefore necessary that all triple Higgs couplings except forλ 113 are small while not overly suppressing κ 3 . The chosen parameter values of BP5 are given in Table II.
The phenomenology of BP5 is shown in Fig. 13. Throughout the parameter plane BR(h 3 → h 1 h 1 ) -shown in Fig. 13 (left) -exceeds 85 % and approaches 100 % for low values of M 3 . The heavy scalar h 3 is produced at a rate of around κ 2 eter region at M 1 120 GeV and M 3 350 GeV is constrained by LHC Higgs searches for resonant double Higgs production [63,66]. These are applied under the assumption that h 1 cannot be experimentally distinguished from h 125 ≡ h 2 if they are close in mass and thus contributes to the predicted signal rate for this process.
The BRs of the h 1 h 1 two-particle state can again be found in Fig. 12. They are identical to those discussed for BP4 since the BRs of h 1 are always identical to those of a SM Higgs boson of the same mass (see Section II B). However, now the scenario extends to M 1 values up to 125 GeV and with increasing M 1 the final state bbW + W − becomes sizable. In contrast to BP4, the two h 1 may be boosted if M 3 2M 1 , leading to collimated h 1 decay products. This may provide an additional experimental handle, enabling the reduction of combinatoric background, and leading to a potential increase of the trigger sensitivity as well as the applicability of jet substructure techniques. Indeed, a recent ATLAS search for highly collimated photon-jets [123] probes the signature pp → h 3 → h 1 h 1 → γγγγ in the mass range M 3 ≥ 200 GeV, 0.1 GeV ≤ M 1 ≤ 10 GeV. However, the currently obtained limit is still several orders of magnitude larger than the predicted rate in BP5.
In benchmark plane BP6 we identify h 125 ≡ h 1 and consider the production of the heaviest scalar h 3 followed by its symmetric decay h 3 → h 2 h 2 . This constrains the mass range for h 3 to values M 3 > 250 GeV. This, in combination with the suppression of κ 3 due to the sum rule, Eq. (18), leads to relatively low production cross sections. The input parameters for BP6 are listed in Table II. 400 GeV, where h 2 decays directly to SM particles. The signal rates in the mass range M 3 600 GeV, which is interesting for cascade decays, can reach up to 100 fb for pp → h 3 → h 2 h 2 at the 13 TeV LHC. In BP6 the total width of of h 3 can reach up to Γ 3 /M 3 ∼ 14 % without violating the unitarity constraint. Therefore, it may be important to include finite width effects in experimental analyses of this scenario.
The shaded region at large masses, M 3 800 GeV, indicates that the parameter region is in conflict with perturbative unitarity. Additionally, experimental searches [66]  M 2 160 GeV, M 3 330 GeV, as shown in Fig. 14. We expect this search analysis to sensitively probe this benchmark scenario in the future, in particular, if the currently considered mass range is extended. The ATLAS search only considers h 3 masses up to the tt threshold, M 3 ≤ 340 GeV. However, as we discuss here, the W + W − W + W − final state remains the dominant four-particle SM final state even beyond the tt threshold. Figure 15 shows the BRs of the decays of the h 2 h 2 state resulting from the h 3 decay in BP6. At low M 2 < 250 GeV only h 2 → F SM decays are kinematically allowed. As shown in Fig. 15  h 2 to four-particle F SM are given in Fig. 15 (left). The dominant final states of this class are W + W − W + W − and W + W − ZZ, with W + W − tt becoming comparable at high M 2 . Figure 15 (left) also shows the "inclusive" branching ratio for the single cascade h 125 h 125 F SM (summed over all possible h 2 → F SM ) and the double cascade decay rate to the h 125 h 125 h 125 h 125 final state.
The branching ratios of h 2 h 2 into various six-particle SM final states via the single cascade decay are shown in Fig. 15  comparable to the four-particle final states. The decay h 2 h 2 → h 125 h 125 h 2 → bbbbW + W − is the third most likely decay mode of h 2 h 2 for 250 GeV < M 2 < 350 GeV.
The branching ratios of the h 125 h 125 h 125 h 125 → F 8 SM decays via a double cascade are independent of the model parameters. They are given in Table IV. Since the BR for the double cascade shown in Fig. 15 (left) is almost independent of M 2 we include in the second row of Table IV an approximate branching ratio for the decay of h 2 h 2 into an eight-particle SM final state through the double cascade. For this we use the averaged BR(h 2 h 2 → h 125 h 125 h 125 h 125 ) = 14.5 %, evaluated in the mass range 260 GeV < M 2 < 500 GeV. The most important eight-particle final states are all combinations of decays into b quarks and W bosons -the most likely decay products of h 125 . Due to combinatorial factors their overall branching fractions are, again, in some cases comparable to the four-and six-particle final states. For example, the bbbbbbW + W − is similar to the ZZZZ final state rate for masses M 3 ∼ 300 −350 GeV. Near the kinematic threshold, M 3 500 GeV and M 2 250 GeV, the signal cross section for pp → h 3 → h 2 h 2 → h 125 h 125 h 125 h 125 amounts to around 14 fb.

VI. CONCLUSION
We presented the collider phenomenology of a simple extension of the SM Higgs sector, where two real scalar singlet fields are added to the particle content. In this two-real-singlet model we imposed a discrete Z 2 symmetry for each scalar singlet field that is spontaneously broken by the singlet field's vacuum expectation value. Consequently, all scalar fields mix, leading to three neutral CP-even Higgs states h a (a = 1, 2, 3). Any of these states can be identified with the Higgs boson with mass 125 GeV observed at the LHC.
The model leads to an interesting collider phenomenology for searches for the additional Higgs states. Following the single production of one of the Higgs states, h a , this state can either decay directly to SM particles, or it can decay into two lighter Higgs states, h a → h b h c , where the lighter states can either be identical ("symmetric" Higgs-to-Higgs decays with b = c = 1, 2), or different ("asymmetric" Higgs-to-Higgs decays with b = 1, c = 2). In the latter case, successive decays of the second lightest Higgs state to the lightest Higgs state, h 2 → h 1 h 1 may be possible if kinematically allowed. This leads to interesting Higgs-to-Higgs cascade decay signatures, in particular, h 3 → h 1,2 h 2 → F SM h 1 h 1 ("single cascade") and h 3 → h 2 h 2 → h 1 h 1 h 1 h 1 ("double cascade"), as shown in Fig. 5. We find that rates for all these possible Higgs-to-Higgs decays can in general be sizable, easily dominating the direct decay modes to SM particles.
Many of these Higgs-to-Higgs decay signatures have not been investigated experimentally to date. We therefore presented six two-dimensional benchmark scenarios to facilitate the design of dedicated experimental searches. Each scenario is defined such that one of the novel signatures has a nearly-maximal signal rate, while still obeying all theoretical and experimental constraints on the model. Moreover, as the model can be parametrized conveniently in terms of the relevant physical parameters, i.e. the three Higgs masses, three mixing angles (governing the Higgs coupling strengths to SM particles) and the three vevs, the benchmark scenarios can cover the entire kinematical phase space for the decay signatures, thus rendering them as ideal references for experimental searches.
For each benchmark scenario, we discussed in detail the rates of the relevant decays, as well as the expected signal rates in the TRSM at the 13 TeV LHC. We furthermore provided an overview of the most relevant SM particle final states, as a function of the relevant mass parameters. This should provide a first step for experimental analyses to estimate the discovery potential of corresponding searches. We expect that some of the presented signatures can already be probed sensitively at the LHC with the data of ∼ 150 fb −1 per experiment collected during Run-II.
It should be kept in mind that the Higgs-to-Higgs decay signatures (and potentially also the cascade decays) discussed here can generically appear also in other BSM models that feature three (or more) Higgs states. In that case, however, the Higgs coupling properties do not necessarily agree with those of the TRSM. This may result in different production rates of the resonantly-produced Higgs state, as well as different decay rates, in partic-ular concerning the Higgs decays to SM particles. It is therefore important that future experimental searches present their results as limits -or ideally measurements -of the model-independent signal rate, as a function of the relevant kinematical quantities (Higgs masses and, possibly, total widths). Furthermore, Higgs-to-Higgs decays to possible SM particle final states that are not dominant in the TRSM may still be worthwhile to probe experimentally, as the anticipated rates may be different in other models. In the case of a future discovery of an additional scalar state, signal rate measurements in various complementary production and decay modes will be crucial to probe its coupling structure and, in turn, to discriminate between the possible BSM interpretations.
The exploration of the scalar sector -leading to a better understanding of the mechanism of electroweak symmetry breaking -is one of the most important scientific goals of the LHC program. This endeavor requires an open-minded and unbiased view on the potential collider signatures of new scalars. Our discussion of the TRSM and the presented benchmark scenarios demonstrate that there is a plethora of currently unexplored collider signatures involving Higgs-to-Higgs decays, and we hope that this work will initiate and facilitate the design of corresponding LHC searches in the near future.

ACKNOWLEDGMENTS
We thank the LHC Higgs Cross Section Working Group WG3 conveners for encouraging this study, as well as Massimiliano Grazzini and Georg Weiglein for useful comments. We also thank Claudia Seitz for helpful discussions regarding experimental questions and Liang Li for helpful explanations regarding the experimental search presented in Ref. [61]. The work of TS and JW is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2121 "Quantum Universe" -390833306. TR was supported by the European Union through the Programme Horizon 2020 via the COST action CA15108 -Connecting insights in fundamental physics (FUN-DAMENTALCONNECTIONS), and additionally wants to thank the DESY Theory group for repeated hospitality while this work was completed.
In this model soft U (1) breaking terms -such as a 1 , b 1 = 0 -are required to avoid a massless Goldstone boson. With these terms included, the resulting model is no longer a special case of the TRSM.