The NMSSM is within Reach of the LHC: Mass Correlations&Decay Signatures

The Next-to-Minimal Supersymmetric Standard Model (NMSSM), the singlet extension of the MSSM which fixes many of the MSSM's shortcomings, is shown to be within reach of the upcoming runs of the Large Hadron Collider (LHC). A systematic treatment of the various Higgs decay channels and their interplay has been lacking due to the seemingly large number of free parameters in the NMSSM's Higgs sector. We demonstrate that due to the SM-like nature of the observed Higgs boson, the NMSSM's Higgs and neutralino sectors have highly correlated masses and couplings and can effectively be described by four physically intuitive parameters: the physical masses of the two CP-odd states and their mixing angle, and $\tan\beta$, which plays a minor role. The heavy Higgs bosons in the NMSSM have large branching ratios into pairs of lighter Higgs bosons or a light Higgs and a $Z$ boson. Search channels arising via these Higgs cascades are unique to models like the NMSSM with a Higgs sector larger than that of the MSSM. In order to cover as much of the NMSSM parameter space as possible, one must combine conventional search strategies employing decays of the additional Higgs bosons into pairs of SM particles with Higgs cascade channels. We demonstrate that such a combination would allow a significant fraction of the viable NMSSM parameter space containing additional Higgs bosons with masses below 1 TeV to be probed at future runs of the LHC.


Introduction
The discovery of the 125 GeV Standard Model (SM)-like Higgs boson [1,2] has prompted the search for additional Higgs bosons at the Large Hadron Collider (LHC). The most straightforward context for such searches is provided by two Higgs Doublet Models (2HDMs) [3], which extend the SM's particle content by a second Higgs SU (2) L doublet. The simplest supersymmetric realization of a 2HDM is the Minimal Supersymmetric Standard Model (MSSM). The collider signatures of such heavy Higgs boson at the LHC have been extensively studied in the literature, see e.g. Refs. [4][5][6][7][8].
The discovery of the SM-like 125 GeV Higgs boson also sparked renewed attention in the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [9,10] since it not only solves the µ-problem [11] of the MSSM but also alleviates the fine-tuning associated with the 125 GeV Higgs boson and the tension implied by the current lack of evidence for superpartners below the weak scale (see e.g. [12][13][14][15]). The NMSSM augments the field content of the MSSM by a SM-singlet chiral superfield S; this extends the particle content by singlet scalar and pseudo-scalar bosons H S and A S , which mix with their corresponding Higgs-doublet counterparts, and a singlet fermion, the singlino S, which mixes with the neutralinos. One of the three CP-even states in the NMSSM must be identified with the 125 GeV SM-like state observed at the LHC. In the following, we reserve the notation h 125 for this SM-like Higgs boson.
The presence of these singlet states introduces new interactions and decay channels, enriching the collider phenomenology of the NMSSM Higgs sector compared to the MSSM. In particular, so-called Higgs cascade decays appear prominently, where a heavy Higgs decays into two lighter Higgs bosons or a light Higgs and a Z boson, see Fig. 3. The presence of Higgs cascade decays warrants the extension of search strategies for additional Higgs bosons beyond the conventional search channels for heavy Higgs bosons developed mostly for models with a Higgs sector consisting of only two Higgs doublets, such as the MSSM [16][17][18][19][20][21][22][23] 1 . Many authors have studied the Higgs and neutralino LHC phenomenology in the NMSSM (see e.g. Refs. [16][17][18][19][20][21][22][27][28][29][30][31] and references therein). However, these studies only cover specific regions of the NMSSM parameter space and typically consider one search channel at a time. A systematic study of the possible signals of the NMSSM Higgs sector and their correlations in parameter space has been perceived to be a challenging task due to the seemingly large number of free parameters controlling the theory.
In this paper, we provide the first systematic approach towards categorizing the NMSSM Higgs sector. We simplify the parameter space of the theory by making use of the SM-like nature of the observed 125 GeV Higgs boson. In the region of NMSSM parameter space where the non SM-like Higgs bosons are light enough to be potentially accessible at the LHC, approximate alignment without decoupling (see e.g. Refs. [18,[32][33][34][35]) must be realized. Such alignment implies correlations between the masses, mixing angles, and couplings in the Higgs and neutralino sector of the NMSSM. We stress that while our analytical understanding of the physically viable parameter space in the NMSSM is guided by assuming perfect alignment in the Higgs sector, we have verified our claims by extensive numerical scans over the NMSSM parameter space using NMSSMTools [36][37][38][39][40] where alignment was not assumed a priori, see Figs. 1 and 2.
The correlations between masses, mixing angles, and couplings become rather clouded when parameterizing the NMSSM in terms of the 7 parameters appearing in the Higgs scalar potential. We show that the region of NMSSM parameter space containing a SMlike 125 GeV Higgs boson and additional Higgs bosons with masses below ∼ 1 TeV, i.e. the region most relevant for Higgs searches at the LHC, can be effectively described by only four physically intuitive parameters: the two physical masses of the CP-odd states, one mixing angle in the CP-odd sector, and tan β. Note that the low tan β regime is of particular relevance for the NMSSM; there, modifying tan β has only minor effects on the NMSSM's phenomenology. Hence, we find that the phenomenology of the entire Higgs and neutralino sectors is governed largely by only three physical parameters in the CP-odd sector, see Eq. (2.34).
Our NMSSM re-parameterization in terms of the masses and mixing angles allows for transparent identification of the most relevant search strategies for different regions of parameter space. These insights allow us to analytically and numerically study, for the first time, the potential of a combination of different search channels arising via Higgs cascade decays. This categorization and coordination of possible Higgs decay channels is important both for extending the coverage of the NMSSM parameter space, as well as identifying the underlying model giving rise to a potential discovery of additional Higgs bosons in the future, e.g. distinguishing the MSSM from the NMSSM. We show that combining Higgs cascade searches with more conventional search modes via decays of non SM-like Higgs bosons into pairs of SM particles, the LHC collaborations will be able to probe ≈ 50 % of the currently remaining viable NMSSM parameter space containing additional Higgs bosons below 1 TeV in future runs of the LHC. We also entertain the scenario that the LHC collaborations can improve the sensitivity of the Higgs cascade decay based searches by an order of magnitude with respect to our projections and the sensitivity in conventional search channels by two orders of magnitude with respect to current limits based on O(30) fb of data. Then, ≈ 90 % of the remaining parameter space containing Higgs bosons below 1 TeV could be probed in the upcoming runs of the LHC, see Figs. 8, 9 and 10. While this latter scenario is optimistic, such sensitivities should be understood as a target for the experimental collaborations which would allow them to probe much of the remaining phenomenologically interesting NMSSM parameter space.
Note that a similar approach can be used to tackle a generic 2HDM+(complex) singlet model [41]. However, the lack of relations in the Higgs sector's parameters prevents making concrete predictions for LHC phenomenology and the interplay of the search modes.
The remainder of this paper is organized as follows. In section 2, we describe the NMSSM parameter space, the correlations in the Higgs and neutralino sector, and our re-parameterization. We validate our analytic claims with extensive numerical parameter scans. In section 3 we discuss various decay channels, their correlations, and their sensitivity at the high luminosity LHC. The coordination of search strategies to cover the parameter space of the NMSSM is presented in section 4. We reserve section 5 for our conclusions. Details regarding the implemented LHC constraints, benchmark points, collider simulations, and analytic expression for the Higgs trilinear couplings are presented in Appendices A-D.

NMSSM Parameter Space
The Next-to-Minimal Supersymmetric Standard Model augments the MSSM particle content with a chiral superfield S uncharged under any of the SM gauge groups. In this paper, we study the scale-invariant NMSSM, where all dimensionful parameters in the superpotential are set to zero. This model enjoys an accidental Z 3 symmetry under which all chiral superfields transform by e 2πi/3 . The additional terms in the superpotential with respect to the MSSM are where H u , H d are the up-and down-type Higgs doublets and λ and κ are dimensionless coefficients. The µ H u · H d term of the MSSM is forbidden in the superpotential of the scale-invariant NMSSM; however, an effective µ-term is generated in the scalar potential when the scalar component of the superfield S gets a vacuum expectation value (vev), µ = λ S . If the vev of the singlet is induced by the breaking of supersymmetry, S is of the order of the supersymmetry breaking scale, thereby alleviating the µ-problem for low-scale supersymmetry. 2 The terms in the scalar potential involving only the Higgs doublets and the singlet are given by [18] V Hu,H d ,S = m 2 S S † S + m 2 where the m 2 i and A i are soft SUSY-breaking parameters of dimension mass squared and mass, respectively, and g 1 and g 2 are the U (1) Y and SU (2) L gauge couplings.
Trading the parameters {m 2 H d , m 2 Hu , m 2 S } for the corresponding vevs via the minimization equations, fixing 3 GeV, and defining tan β ≡ v u /v d , the scalar potential is controlled by the following parameters {λ, κ, tan β, µ, A λ , A κ }. (2.3) Note that all parameters are real in the CP-conserving NMSSM. Of the dimensionless parameters, λ and tan β can be chosen positive without loss of generality, while κ and the dimensionful parameters can have both signs. It it useful to rotate the doublet-like states to the (extended) Higgs basis [18,32,33,[42][43][44][45][46]

5)
where the H 0 i are the neutral components of the corresponding doublet fields H i . The couplings to pairs of SM particles take the particularly simple form

8)
A NSM (down, up, VV) = (g SM tan β, g SM / tan β, g SM ) , (2.9) where "down" ("up") stands for pairs of down-type (up-type) SM fermions, "VV" for pairs of vector bosons, and g SM is the coupling of an SM Higgs boson of the same mass to such particles. The CP-even and CP-odd interaction states from the singlet S do not couple to SM particles and are defined via The charged Higgs is defined by The remaining degrees of freedom make up the longitudinal polarization of the W ± and Z bosons after electroweak symmetry breaking. The elements of the symmetric squared mass matrix for the CP-even Higgs bosons in the extended Higgs basis {H SM , H NSM , H S } at tree-level are where we traded A λ for M 2 A , defined as and used the short-hand notation The tree-level elements of the symmetric squared mass matrix for the CP-odd Higgs boson in the basis {A NSM , A S } are given by 22) and the mass of the charged Higgs boson is where s W ≡ sin θ W , with θ W the weak mixing angle. In this paper, we decouple the gauginos from the collider phenomenology by setting the bino and wino mass parameters {|M 1 |, |M 2 |} |µ|.

Alignment
The (neutral) interaction states of the Higgs basis mix into three CP-even and two CP-odd mass eigenstates. We denote the CP-even mass eigenstates h i = {h 125 , H, h}, where again m A > m a , and P j a i denotes the j = {NSM, S} component of the a i mass eigenstate. The S j h i and P j a i are obtained by diagonalizing the squared mass matrices for the CP-even states, Eqs. (2.12)-(2.17), and CP-odd states, Eqs. (2.20)-(2.22), respectively.
The measured branching ratios of the 125 GeV mass eigenstate observed at the LHC are compatible with those of a SM Higgs boson, although current experimental precision allows for O(10 %) deviations [47][48][49][50][51]. Thus, in order to be compatible with the observed phenomenology, the h 125 eigenstate we identify with the observed Higgs boson must have a mass of ∼ 125 GeV and be dominantly composed of the interaction eigenstates H SM , whose couplings are identical to the SM Higgs boson's.
As is well known, the (squared) mass of the h 125 mass eigenstate receives an additional contribution λ 2 v 2 s 2 2β relative to the MSSM case, and at tree-level is given by For small to moderate values of tan β the λ 2 contribution to the mass is sizable and allows for a tree-level mass of h 125 close to 125 GeV. This makes the low tan β region in the NMSSM particularly interesting because there the observed mass of the SM-like Higgs boson can be obtained without the need for large radiative corrections, as for example are required in the MSSM. There are two possibilities to achieve approximate alignment of H SM with h 125 [18,34,[52][53][54][55][56][57]: Either, the remaining mass eigenstates H and h are much heavier than 125 GeV, the so-called decoupling limit, or, the parameters of the model conspire to (approximately) useful to re-parameterize the physically relevant region of parameter space by approximate alignment and the physical masses of the CP-odd Higgs bosons. The remaining freedom of the parameter space can be described by the mixing angle in the CP-odd sector P S A [Eq. (2.26)] and the value of tan β. Hence the basis is physically much more intuitive than the usual parameterization in terms of the parameters appearing in the scalar potential, cf. Eq. (2.3). While current experimental constraints allow for λ and κ to be slightly shifted from the values expected from perfect alignment, in practice, we can use the alignment conditions, Eqs. (2.29) and (2.30), to fix λ and κ to a very good approximation [18,21]. Further, it is well known (and easily seen from the mass matrices) that the precise value of tan β is a small effect in the low tan β regime, which is of most interest in the NMSSM. Hence, the phenomenology of the Higgs and neutralino sectors is, to a large degree, governed by the three parameters The value of µ can be obtained from M 2 P,12 as given in Eq. (2.32), and using the relationship for κ/λ as dictated by alignment, Eq. (2.30), with the corresponding . (2.37) Finally, from the matrix element M 2 , and the alignment relations have been assumed for κ and λ. Note that for each set of input parameters {m a , m A , P S A }, there are two sets of correlated solutions for µ, κ and A κ . In our analytical formulae and figures, we will denote these by µ ± . We also note that P S A ↔ (−P S A ) corresponds to µ ± ↔ (−µ ∓ ).
We stress that these relations define the masses as well as all couplings between the NMSSM Higgs bosons and between Higgs bosons, neutralinos 5 , and SM particles from the input parameters {tan β, m A , m a , P S A }, assuming (approximate) alignment as dictated by h 125 phenomenology.
A comment about radiative corrections is in order here. In general, sizable corrections are present in the NMSSM, in particular via stop loops due to the large top Yukawa couplings as well as via Higgs loops via potentially large quartic couplings between the Higgs bosons, see e.g. Refs. [9,18,29,58]. Since our re-parameterization of the parameter space is obtained at tree level except for the first alignment condition, Eq. (2.29), not all such corrections are explicitly included. This may somewhat cloud the relation of our parameter basis, which uses physical masses, a mixing angle, tan β, and the couplings λ and κ, with the usual parameterization of the Higgs sector in terms of the parameters appearing in the scalar potential. Our relations in Eqs. (2.35)-(2.38) should strictly be understood as relations to obtain parameters shifted with respect to the bare parameters after absorbing relevant radiative corrections. Note also that we did not include obtaining m h 125 = 125 GeV as a condition on our parameter basis, rather, the required mass of the SM-like eigenstate should be understood as setting the size of the stop corrections. In the NMSSM, a 125 GeV mass for the SM-like Higgs mass eigenstate can be obtained without large radiative corrections. Thus, the phenomenology of the Higgs and neutralino sector can be studied in a region of parameter space where the radiative corrections from the stops are small and the relation of the parameters obtained from our Eqs. (2.35)-(2.38) with input parameters for numerical tools is rather direct. Thus, even though the parameters obtained from Eqs. (2.35)-(2.38) cannot generally be directly used as input in spectrum generators like NMSSMTools, SOFTSUSY, NMSSMCALC, etc, in practice, this is a minor problem as discussed further in the following section. Finally, the main advantage of our re-parameterization is that it allows for the transparent understanding of the Higgs sector in the physically viable region of parameter space. While a precision study would require one to carefully incorporate radiative corrections, here we are interested in mapping the qualitative behavior and in identifying search strategies to cover as much of the NMSSM's parameter space as possible. Radiative corrections shift the parameters, but do not impact the qualitative behavior of the NMSSM.

Mass Correlations
In the previous section we re-parameterized the NMSSM parameters governing the Higgs and neutralino sectors of the NMSSM. We showed that in the alignment limit, only four free parameters remain in the Higgs sector. Most of the phenomenology is controlled by {m a , m A , P S A }, while tan β plays a minor role. This leads to strong correlations between the masses in the Higgs sector as well as between the Higgs masses and the Higgsino and singlino parameters.
From the form of the mass matrices, it is straightforward to see that the scale of all masses, except for the SM-like Higgs state, are controlled by the parameter |µ|, which is in turn highly correlated with M A due to the requirement of approximate alignment [see Eq. (2.36)]. Because of the large numerical factor in the square root in Eq. (2.36), regardless of the mixing angles and the mass splitting in the CP-odd sector, Combining this with the alignment condition given in Eq. (2.30) dictates that |κ|/λ should be small for most of the region under consideration, however it can be driven to larger values even for small deviations of µ from Eq. (2.39), where the χ i with 1 ≤ i ≤ 5 are the neutralino mass eigenstates in ascending order of their masses m χ i , and the N ij denote the interaction eigenstate components of the Since tan β is small, the neutralino sector is mostly controlled by the mass splitting between m χ i and µ, and hence by the value of κ/λ given by Eq. (2.37). 6 Thus, we expect the masses of the Higgsino-like neutralinos to be correlated with the mass of the doublet-like Higgs states, while the mass of the singlino-like state is much more strongly effected by m a and P S A . Considering the Higgs sector, we first note from Eq. (2.30) that 2s 2β |κ|/λ ∝ /|µ| 1, and hence the mixing angles and the heavy masses in the CP-odd and even sectors are generally expected to be correlated, cf. Eqs. (2.31) and (2.32), with masses approximately given by M A . The singlet-like states are less tightly correlated; taking into account firstorder mixing effects, their masses can be approximated as [18] (2.43) 6 This holds under the assumption that the absolute value of the bino and wino mass parameters |M1| and |M2| are much larger than |µ|. However, note that allowing the bino and the wino to be light does not add new parameters beyond M1 and M2 since the bino and wino mix only with the Higgsinos, and the mixing is controlled only by tan β. 7 For compactness, we denote the masses of the singlet-like CP-even and CP-odd mass eigenstate by m h and ma here, respectively. A priori, the singlet-like states can be heavier than the doublet-like states, in such a case the singlet-like states should be identified with A or H in our notation. It should be noted however that for typical choices of parameters the singlet-like masses are lighter than MA, such that they comprise the lighter CP-odd and non SM-like CP-even mass eigenstates.  Table 1. NMSSM parameter ranges used in NMSSMTools scans. We decouple the remaining supersymmetric partners by setting all sfermion mass parameters (except the stop parameters M Q3 = M U3 ) to 3 TeV, the bino and wino mass parameters to M 1 = M 2 = 1 TeV, and the gluino mass to 2 TeV. The stop and sbottom mixing parameters are set to X t ≡ (A t − µ cot β) = 0 and See the text for a discussion of these choices as well as the ranges of the parameters we scan over.
Note that the sum (3m 2 h + m 2 a ) is independent of A κ [18]. The opposite sign contribution from A κ to m 2 h and m 2 a induces an anti-correlation in their masses for fixed values of the remaining parameters. Further, compared to the CP-even state, the singlet-like CP-odd state receives a factor of 3 larger contribution from A κ , and no large contribution from either M A or µ. Thus, m a has the smallest correlation with M A (or µ) of the non SM-like Higgs states, justifying our choice of parameterization for the Higgs sector in Eq. (2.33). On the other hand, apart from its anti-correlation with m a , the mass of the singlet-like CP-even state receives large contributions proportional to κ 2 µ 2 /λ 2 , and is thus expected to be quite correlated with the masses of the doublet-like states, as well as the CP-odd mixing angle. It can be further shown that in the parameter region of interest, the maximal value for the CP-even state is obtained for the smallest values of m a , and generally obeys m h M A /2 [18].
In order to demonstrate these correlations, we will show the statistical properties of the masses of NMSSM spectra obtained from a parameter scan with NMSSMTools 4.9.3 [36][37][38][39][40] in Figs. 1 and 2. NMSSM parameters are drawn from linear flat distributions over the ranges listed in Tab. 1. Note that for our numerical scans NMSSMTools requires us to use the parameters appearing in the scalar potential, {tan β, λ, κ, A λ , A κ , µ}, and the stop mass parameter, M Q 3 , as input, not the parameters of our more physical re-parameterization given in Eq. (2.33). In addition to the standard scan we also perform a scan over a narrower range of parameters focused on producing lighter Higgs spectra accessible at the LHC which we label the light subset. The range 1 ≤ tan β ≤ 5 is motivated by obtaining m h 125 = 125 GeV without the need for large radiative corrections; recall that the contribution λ 2 v 2 s 2 2β to M 2 S,11 is suppressed for larger values of tan β. We decouple the remaining supersymmetric partners from our study by setting all sfermion mass parameters (except the stop parameters) to 3 TeV, the bino and wino mass parameters to M 1 = M 2 = 1 TeV, and the gluino mass to 2 TeV. Since large third generation squark mixing is not necessary to obtain the correct SM-like Higgs mass in the NMSSM, we set the stop and sbottom  mixing parameters X t ≡ (A t − µ cot β) = 0 and X b ≡ (A b − µ tan β) = 0. Parameter points from the scan are kept if they satisfy a subset of the standard constraints implemented in NMSSMTools, in particular, the Higgs spectra must contain a SM-like Higgs boson with mass and couplings compatible with the SM-like 125 GeV state observed at the LHC, as well as evade constraints from searches for additional Higgs bosons and sparticles at the Large Electron-Positron Collider (LEP), the Tevatron, and the LHC. Furthermore, we require the lightest neutralino χ 1 to be the lightest supersymmetric particle (LSP). Beyond the constraints implemented in NMSSMTools, we require compatibility with direct searches for Higgs bosons at the LHC listed in Tab. 2 located in the Appendix. As discussed in Ref. [21] we find that parameter points with a Higgs boson with mass and couplings compatible with the 125 GeV state observed at the LHC approximately satisfy the alignment conditions, although our chosen parameter ranges do not a priori impose these conditions. Before we validate our analytical claims with numerics, let us highlight a few points  regarding the coverage of the NMSSM parameters in terms of the physical basis we have chosen. First we note that the requirement of non-tachyonic m 2 h means that not all values of {m a , m A , P S A } are physically allowed. Second, as discussed above, depending on the choice of the CP-odd mass parameters, the alignment conditions may lead to large values of |κ|/λ. However, perturbative consistency generally demands |κ|/λ 1/2 [9,18]. Generically, this implies that for a fixed value of m A , large mixing angles would demand too large values of |κ|/λ for small values of m a , whereas values of m a close to degeneracy with m A tend to drive m h tachyonic. We also note that large mixing angles for light m a can be in tension with direct searches for doublet-like scalars at the LHC, cf. Tab. 2, further reducing the allowed range of mixing angles in such situations. Therefore, even though we started off with a linear flat distribution for our scans of the NMSSM parameter space, the resultant physically viable regions predominantly correspond to CP-odd masses with small mixing angles. The sum of these requirements leads to a rather small dependence of the mass spectra on the mixing angles beyond the ordering of the neutralino masses and the mass of m h , despite our random scan a priori allowing for large mixing angles. Hence, we present our numerical results for the mass spectra in the m A − m a plane.
As commented at the end of the previous section, radiative corrections affect the relations between our parameter basis and the inputs used in the NMSSMTools scan. Since we set the stop and sbottom mixing parameters to zero in our numerical scan, the stop corrections to the Higgs mass matrices are relatively simple, cf. Refs. [9,18] for the relevant expressions. However, depending on the value of the stop mass parameter M Q 3 as well as the size of the quartic couplings between the Higgs bosons, sizable radiative corrections may still be present. Hence, care must be taken when performing precision studies of the NMSSM parameter space to ensure that radiative corrections are properly incorporated when comparing NMSSMTools numerical output to our analytical alignment conditions. The phenomenologically most interesting region of parameter space is where the additional Higgs bosons have masses below 1 TeV and are hence accessible at the LHC. In this region, excellent agreement is obtained between analytic alignment conditions and the full numerical output from NMSSMTools as shown in Fig. 1 of Ref. [21]. We have further checked that our Eqs. (2.35)-(2.38) yield broad agreement when comparing the NMSSMTools input parameters with the corresponding quantities obtained from these equations for the points in our numerical scan.
In Fig. 1 we show the masses of the non SM-like CP-even states H and h, and in Fig. 2 the masses of the singlino-and Higgsino-like neutralinos χ i , i = {1, 2, 3}, together with the value of |µ| in the m A − m a plane. In these figures, the color scale shows the mean of the respective mass (or |µ|) binned in the m A − m a plane. In addition, we show the standard deviation of the entries in each bin in units of the mean value with the blue error-bars. The error-bars are normalized such that the they would span the height of the bin if the standard deviation is equal to the mean. Note that the scale of the error-bar is linear with respect to the bin height, while the bin widths as well as the color scale are logarithmic.
From the left panel of Fig. 1, we see that as expected the heavy (usually doublet-like) CP-even state H and the CP-odd state A are approximately mass degenerate and tightly correlated. The mass of H is virtually independent of m a . From the right panel of Fig. 1 we observe that the mass of the light (usually singlet-like) CP-even state h is also correlated with the mass of the heavy doublet-like states, with larger m A leading to heavier m h . Further, as expected the mass satisfies m h m A /2. From the same panel, we also observe that the lightest m h is obtained for the heaviest m a and vice versa, showing the expected anti-correlation in their masses. The large standard deviation ∼50 % for most values of {m a , m A } shows the weaker dependence on the mixing angle. Fig. 2 shows the correlation of the µ parameter and the masses of the three lightest neutralinos with the Higgs masses. We first note that the value of |µ| is very tightly correlated with m A , with correspondingly small error-bars, particularly in the low m A region. This is in agreement with Eq. (2.39), and stems from viable h 125 phenomenology requiring approximate alignment. We stress again that these parameter points are selected from a random parameter scan by requiring compatibility with the observed h 125 phenomenology without a priori imposing alignment conditions. Thus, these results justify our use of the alignment conditions to reduce the number of free parameters, which facilitates the analytic understanding of the parameter space. The |µ| parameter not only controls the Higgsino masses, but the singlino mass is also proportional to µ. In our parameter scans we decoupled the bino and wino mass parameters {|M 1 , |M 2 |} |µ|, hence the three lightest neutralinos are dominantly composed of the Higgsinos and the singlino. Since the Higgsinos are mass degenerate before taking into account mixing effects, we expect either the lightest or the third-lightest neutralino to be singlino-like, while the remaining two of the three lightest neutralinos are Higgsino-like. The effect on the masses can be seen in Fig. 2: The second lightest neutralino is usually Higgsino-like and its mass thus quite tightly correlated with µ. The mass-scales of the lightest and third-lightest neutralino on the other hand are also correlated with µ, but we see both larger standard deviations as well as masses smaller than the mean of |µ| for χ 1 and larger than |µ| for χ 3 . Both effects are due to either χ 1 or χ 3 being singlino-like.
In summary, we find that the mass spectra of the Higgs sector as well as the associated neutralinos can be described by a simple parameterization in terms of four quantities: the two physical masses and the mixing angle in the CP-odd sector, and tan β (the latter plays a minor role). The presence of a SM-like state with a mass of 125 GeV requires approximate alignment for Higgs spectra accessible at the LHC, determining the preferred values of tan β, λ and κ/λ. The non SM-like CP-even doublet-like state and the CP-odd doubletlike state are approximately mass degenerate and heavier than the singlet-like states. We use the masses of the CP-odd states m a and m A as an input parameters. Together with the value of tan β, one can then obtain the value of µ which controls the Higgsino masses. Since the value of κ/λ is given by the alignment condition one also directly obtains the singlino mass. In the limit where the bino and the wino are heavy, {|M 1 |, |M 2 |} |µ|, we find that the second-lightest neutralino is Higgsino-like with mass given by |µ|. Either the lightest (if 2κ/λ 1) or third lightest neutralino (if 2κ/λ 1) is singlino-like with mass ∼ |2κµ/λ|, and the remaining state is again Higgsino-like with its mass pushed away from |µ| due to the mixing effects with the singlino. The remaining Higgs states are the singlet-like CP-odd and CP-even states. The mass of the CP-even state is mostly governed by m A and satisfies m h m A /2. The mass of the CP-odd singlet-like state is only weakly correlated with the mass of the remaining states. Most prominently, the mass of the singlet-like scalar is (weakly) anti-correlated with the mass of the singlet-like pseudo-scalar. Figure 3. Illustration of NMSSM-specific Higgs decay channels, where the Φ i,j,k stand for one of the five NMSSM Higgs bosons. For channel (a), either one or all three of the Φ i,j,k must be CPeven. For channel (b), if Φ i is CP-even, Φ j must be a CP-odd state, and vice-versa. For channel (c), the final state can be χ 1 χ 1 h i , χ 1 χ 1 a i , or χ 1 χ 1 Z, and Φ i can be CP-even or -odd. As discussed further in the text, the most important channels considered in this work are (gg → H → hh 125 )

Higgs Decays
The Higgs bosons in the NMSSM can decay into a variety of final states, however the bulk of the experimental searches at the LHC have been focused on Higgs bosons decaying into pairs of SM particles, see Tab. 2. The presence of the singlet-like states in the NMSSM both poses a challenge and offers new opportunities for Higgs searches at the LHC when compared to the MSSM. On the one hand, since the singlet does not directly couple to any SM particle, production cross sections of the NMSSM Higgs bosons at colliders are suppressed by the respective singlet component of the Higgs boson in question. On the other hand, the additional singlet-like states offer new decay modes for the doublet-like Higgs bosons, illustrated in Fig. 3. As discussed e.g. in Refs. [18,21], branching ratios into pairs of lighter Higgs bosons or a light Higgs and a Z boson can be sizable and even compete with decays into pairs of top quarks. Note that decays into pairs of SM-like Higgs bosons or a SM-like Higgs and a Z boson are suppressed, since the corresponding couplings vanish in the alignment limit. Therefore, of all the decays into bosons, the experimentally most promising channels are cascade decays into a SM-like Higgs and an additional non SM-like Higgs boson, or into a Z boson and an additional non SM-like Higgs. The corresponding couplings are not suppressed by the presence of the SM-like h 125 , and the Z or h 125 in the final state allows for tagging of such events due to their known masses and branching ratios. We will discuss such decays in some detail below.

Cascade Decays
In order to study which of the different final states is most relevant for the different regions of NMSSM parameter space, it is useful to start by studying the ratios of σ(gg → Φ 1 → ZΦ 2 ) and σ(gg → Φ 1 → h 125 Φ 2 ) at the LHC, where Φ i stands for any of the non SM-like NMSSM Higgs mass eigenstates. The branching ratio BR(Φ 1 → h 125 Φ 2 ) in particular is intimately related to the couplings λ and κ, while the (Φ 1 → ZΦ 2 ), (Φ → SM SM) and (Φ → χ i χ j ) branching ratios depend mostly on the mass spectrum and the respective mixing angles. The dependence of these Higgs cascade decays on the relevant masses and mixing angles have been studied in great detail in the context of a generic 2HDM+singlet model in Ref. [41]. However, unlike the generic 2HDM+S model, in the NMSSM many parameters are correlated as discussed in the previous section. In the following we will discuss how these parameter correlations dictate the behavior of the Higgs cascade decays.
In terms of the mixing angles and masses and assuming alignment, the most relevant ratios can be written as [41] where σ ggh (m) is the gluon fusion production cross section of a SM Higgs boson of mass m, and the form factor is defined as In addition to the masses and the mixing angles, these ratios depend on the decay widths of the parent state Γ Φ , and one combination of trilinear couplings between the involved states, which are given bỹ For values in proximity of alignment λ ∼ 0.65, moderate values of |κ| < 1, and low tan β, these combinations of couplings are at most ∼ O(v). Hence, unless the channel in the denominator of Eqs.
, these couplings play no important role for the ratios.
Note that the ratio [σ(gg when exchanging all quantities referring to CP-even states to the corresponding quantity for CPodd states, and vice versa. Likewise, the ratio [σ(gg 3) keeping in mind that exchanging the ratio of the gluon fusion production cross sections entails replacing In Figs. 4 and 5 we show these ratios in the plane of the light CP-odd mass m a vs. its NSM fraction (P NSM a ) 2 = (P S A ) 2 for fixed values of m A and tan β. As discussed in section 2.2, in the alignment limit all parameters controlling the NMSSM Higgs sector are fixed in terms of these inputs, except for the choice of the two different solutions for µ ± , cf. Eq. (2.36), corresponding to the two panels in each of the figures. As discussed in section 2.2, due to the correlation of parameters, not all of the parameter region shown in Figs. 4 and 5 is allowed. In particular, κ takes large absolute values for large (P S A ) 2 0.3 and small values of m a 150 GeV. In order to prevent a Landau pole close to the SUSY breaking scale, we constrain |κ| < 1 in the figures. Furthermore, for one of the solutions for µ ± , the light non SM-like CP-even Higgs state becomes tachyonic for large values of m a as indicated in the corresponding panels.
In some regions of parameter space only one of the decay modes appearing in the respective ratio is kinematically allowed. Recalling the correlation of the masses discussed    In the region of parameter space where both channels appearing in the respective ratios are allowed, we find the two cross sections to generally be of the same order of magnitude. However, in particular regions of parameter space one of the channels can dominate, even far from regions where kinematic suppression is effective. This occurs for example in the top panels of Fig. 4 and the bottom panels of Fig. 5 for (P S A ) 2 ∼ 0.5; for this value the cross section in the denominator σ(gg → A → h 125 a) is strongly suppressed, cf. Eqs. (3.1) and (3.4). Therefore, since no single decay channel is dominant throughout parameter space, it is important to consider all of them in order to fully cover the parameter space of the NMSSM at the LHC. We have verified the analytical results presented in this section by computing and comparing these ratios from the output of our NMSSMTools scan.
The final decay products of the daughter Higgs and Z bosons produced from the decay of the heavy parent Higgs in the Higgs cascade decays discussed above will dictate the sensitivity of the LHC to such channels. Higgs and Z bosons decay into pairs of SM particles such as τ + τ − , bb, ZZ, or W + W − , or, if kinematically accessible, they might also decay into pairs of neutralinos. In general, in the low tan β regime, the branching ratios of the light non SM-like Higgs bosons are similar to those of a SM Higgs boson of the same mass, with the exception that they can have sizable branching ratios into pairs of neutralinos if kinematically accessible, and that CP-odd Higgs bosons do not decay into pairs of gauge bosons at tree-level. Note that the decays into pairs of SM fermions are controlled by the (tan β suppressed/enhanced) Yukawa couplings. The decays into pairs of Higgsino/singlino-like neutralinos are controlled by λ and κ instead. In the alignment limit we find λ ∼ 0.65, much larger than all Yukawa couplings but the top Yukawa. Thus, if kinematically accessible, NMSSM Higgs bosons typically have large branching ratios into pairs of neutralinos below the top threshold, i.e. m Φ i 2m t ∼ 350 GeV. This qualitative behavior is unchanged when allowing for light binos or winos. The couplings of the doublet-like Higgs bosons to binos and winos are controlled by the U (1) Y and SU (2) L gauge couplings, respectively, which are again larger than all Yukawa couplings but the top Yukawa. Expressions for the couplings of the NMSSM Higgs bosons to pairs of SM particles as well as pairs of neutralinos can e.g. be found in Refs. [9,29,59], and a more detailed discussion of the branching ratios can be found in Refs. [21,41].
Incorporating the final state decays of the non SM-like daughter Higgs bosons, we classify the Higgs cascade decay channels leading to different final states as follows: If the heavy Higgs boson decays into a SM-like Higgs and a light Higgs [ Fig. 3 (a) Mono-Higgs or mono-Z signatures can also arise if the heavy Higgs decays directly into neutralinos where one of the neutralinos is not the lightest one, and subsequently decays into the lightest neutralino and a SM-like Higgs or a Z boson [ Fig. 3 (c)]. However, as discussed in Refs. [21,41], such decays are kinematically unfavorable for collider searches since the neutralinos might conspire to be produced approximately back-to-back in the transversal plane yielding small E miss T . Note that the categorization above misses some final states, such as when both the h 125 (or the Z) and the light additional Higgs boson decay into invisible states, or if the heavy Higgs decays to two heavier neutralinos which subsequently decay into the lightest neutralino and additional particles. The former type of decay channels may e.g. be probed via monojet-type searches. The latter decay channel may in principle be probed with strategies similar to what is discussed here, although it will in general be more challenging since they yield softer final states.

LHC Prospects for Cascade Decays
Not all final states are equal -the sensitivity of the LHC is very channel dependent. To determine the coverage of the NMSSM parameter space at the LHC we need to compare the cross sections for each channel to the sensitivity of the LHC. To this end we compare the cross sections for our NMSSM parameter scan to the projected sensitivity in the different channels at the 13 TeV LHC assuming L = 3000 fb −1 of data. For the first time, we exploit all of the mono-Z, mono-Higgs, Z+visible, and Higgs+visible classes of final states for probing the NMSSM at the LHC, whereas the previous literature considered one class of final states at a time.
The sensitivity of the mono-Higgs and mono-Z channels has been extensively discussed in Refs. [21,41]. The sensitivity of Higgs+visible channels in the bbbb and bbγγ final states has been discussed in Ref. [22]. The importance of the Z+visible channel has been discussed in [18,21], but to date no estimate of the sensitivity at the 13 TeV LHC is available. For the purposes of this work, we extrapolate the sensitivity at the 13 TeV LHC for L = 3000 fb −1 of data, σ 13 TeV; 3000 fb −1  In Fig. 6 we show the signal strength for our NMSSM parameter scan where σ is the cross section of the parameter point and σ 3000 fb −1 Proj.
is the cross section expected to be probed with L = 3000 fb −1 of data in the respective channel. The different panels in Fig. 6 correspond to the different Higgs cascade channels (mono-Higgs, For the right panels, the x-axis corresponds to m h (m A ) and the y-axis to m H (m a ). The hard cutoff at masses of the parent state below ∼ 350 GeV in the top right panel (Higgs+visible) is due to the mass ranges for which the sensitivity in these finals states is available in Ref. [22]. Note also that in the bottom left panel the mass range extends up to 1.5 TeV, while in the other panels we show only the mass range up to 1 TeV.
Higgs+visible, mono-Z and Z+visible) arising from the decay of a heavy Higgs into a pair of lighter Higgses or a light Higgs and a Z boson, corresponding to the diagrams (a) and (b) in Fig. 3. The signal strength for mono-Z and mono-Higgs final states arising via decays of a heavy Higgs into neutralinos, where one of the neutralinos in turn radiates off a Z or a Higgs boson, cf. diagram (c) in Fig. 3, is shown in Fig. 11 located in Appendix D. 8 For completeness, we present four benchmark points in Appendix B, chosen to have large signal cross sections in the mono-Higgs, Higgs+visible, mono-Z, and Z+visible classes of Higgs cascades, respectively.
All four search channels shown in Fig. 6 are able to probe sizable regions of the NMSSM parameter space. Note that the comparison between different channels should not be taken at face value but as a qualitative comparison, since the extrapolations of the future LHC sensitivity assume different systematic uncertainties for each channel; in particular, the assumptions in Ref. [22] for the Higgs+visible channel are more optimistic than those in the extrapolation of the remaining channels. Most notably, all searches maintain sensitivity in the {m H , m A } 350 GeV region which is difficult to probe with conventional Higgs searches. For {m H , m A } 350 GeV, decays into top quarks {H, A} → tt dominate the decays of heavy Higgs bosons into pairs of SM particles; this decay channel is very challenging to probe at the LHC due to interference with the QCD background [6,[62][63][64][65][66][67][68].
Comparing various final states arising in individual triangles portrayed in each panel of Fig. 6, we clearly find the effects of the correlation of masses discussed in section 2.3. For the mono-Z and Z+visible (mono-Higgs and Higgs+visible) final states, the decay chains induced by the parent CP-odd state A (CP-even state H) contain only the light CP-even state h and neutralinos with tightly correlated masses. This leads to the behavior of the cross sections of the respective channels with respect to the extrapolated LHC sensitivity being relatively uniform, cf. the bottom right triangles in the left panels and the top left triangles in the right panels of Fig. 6 respectively. On the other hand, mono-Z and Z+visible (mono-Higgs and Higgs+visible) final states induced by the parent CP-even state H (CP-odd state A) involve the light pseudo-scalar a, whose mass is much less tightly correlated with the remaining Higgs states. This leads to somewhat less regular behavior, as can be seen from top left triangles in left panels and bottom right triangles in right panels of Fig. 6 respectively. In particular, while one may find mass spectra with large mass gaps leading to readily observable final states, one can also easily end up in situations where the a is too heavy such that this decay of the heavy parent state is kinematically suppressed or may not yield sufficiently hard decay products required for large E miss T in the mono-Z and mono-Higgs final states. On the other hand, one may also end up in the situation where a is lighter than the pair of lightest neutralinos, such that a → χ 1 χ 1 decays required for mono-Z and mono-Higgs final states are kinematically forbidden.
Comparing the Z+visible and Higgs+visible to the mono-Z and mono-Higgs final states, we find that the Z+visible and Higgs+visible channels are usually more effective as long as the light Higgs involved in the decay chain is below the top threshold, m h/a Figure 7. Largest Z+visible or Higgs+visible signal strength (x-axis) vs. the largest mono-Z or mono-Higgs signal strength (y-axis) for points from our NMSSM parameter scan. The color coding denotes the Higgs cascade channel with the largest signal strength as indicated in the legend. The Φ in the legend can be either the light non SM-like CP-even state h or the light CP-odd state a, depending on the CP properties of the parent state in the cascade, and whether it is accompanied by an h 125 or a Z boson, cf. Fig. 3. The displayed parameter points satisfy all current constraints from conventional searches at the LHC as listed in Tab. 2 and feature at least one of the heavy Higgs bosons H or A lighter than 1 TeV. Points in the L-shaped region either above or to the right of the solid lines, indicating a signal strength µ 3000 fb −1 Proj. = 1, are within our projected sensitivity of the LHC with 3000 fb −1 of data. The different Higgs cascade channels are clearly complimentary such that one must employ all of them in order to probe as large a portion of the parameter space as possible. Note that as discussed further in the text, the LHC collaborations may improve the sensitivities by at least one order of magnitude compared to our estimates. Points below (to the left of) the dashed lines have signal strengths µ 3000 fb −1 Proj. < 10 −10 in the mono-Z and mono-Higgs (Z+visible and Higgs+visible) channels, rendering such points difficult to detect in the respective channels even at the high energy LHC [69], and may require the 100 TeV collider [70]. 350 GeV. Once the light state is allowed to decay to a pair of top quarks, such decays will usually dominate, rendering searches in (h/a → bb/τ + τ − /γγ) final states less effective. This effect is particularly visible in the (gg → H → Za) and (gg → A → h 125 a) decay channels, while it is somewhat less pronounced in the (gg → A → Zh) and (gg → H → h 125 h) channels because the kinematic cutoff for a → tt decays is much harder than for h → tt decays and because decays of CP-odd Higgs bosons into pairs of vector bosons are forbidden at tree-level. Mono-Z and mono-Higgs final states remain sensitive above the top threshold for the light Higgs; as discussed above, the branching ratios of the light Higgs bosons into pairs of neutralinos can be comparable with the branching ratios into pairs of top quarks. Thus, mono-Z and mono-Higgs final states are particularly powerful in the parameter region hard to probe with conventional searches where all of the non SM-like Higgs bosons are above the top threshold. This region may also be accessible with Z+visible or Higgs+visible final states when using final states arising from decays of the light Higgs bosons into top quarks or W or Z bosons. Sizable cross sections into such final states have been demonstrated in Ref. [21], however, no estimate of the LHC sensitivity for such final states exists to date.
Finally, in Fig. 7 we compare the signal strengths of the mono-Z and mono-Higgs channels to the Higgs+visible and Z+visible channels for all points from our NMSSM parameter scans passing all constraints, in particular evading all current bounds from conventional searches at the LHC as listed in Tab. 2. The different Higgs cascade channels are clearly complimentary such that one must employ all of them in order to cover as large a portion of the parameter space as possible. Recall that the comparison of the different channels should be understood qualitatively and not be taken at face value, since the extrapolation of the sensitivity for the different channels assume e.g. different systematic errors of the background.
We stress that the sensitivity extrapolations we have used are somewhat conservative, in particular in the mono-Higgs, mono-Z, and Z+visible final states, where systematic errors may be reduced significantly by the experimental collaborations compared to what was assumed when estimating the sensitivity. Hence we expect that the true sensitivity of LHC searches may be up to approximately one order of magnitude better than what is shown in Figs. 6 and 7. This renders in particular the mono-Z final state very promising, allowing the LHC to probe heavy Higgs boson with masses larger than 1 TeV.

Combining Searches to Cover the NMSSM Parameter Space
In the previous section we discussed the sensitivity of the LHC for the different final states arising from Higgs cascade decays. In this section we will demonstrate how by combining these searches with conventional searches (utilizing direct decays of the heavy Higgs bosons into pairs of SM particles), significant progress towards coverage of the NMSSM parameter space can be made. We find that the NMSSM parameter space which realizes non SM-like Higgs bosons lighter than ∼ 1 TeV could be almost completely probed by the 13 TeV LHC with 3000 fb −1 of data.
In Fig. 8 we compare the projected signal strength of the Higgs cascade channels, µ 3000 fb −1

Proj.
, with the current signal strength for the conventional Higgs searches, µ <37 fb −1 Curr. Lim. . Once more we note the complementarity of the different channels. In particular, when  Fig. 7 but that the x-axis shows the largest signal strength of all conventional Higgs searches listed in Tab. 2 arising through the production of the light states h and a, and the y-axis shows the largest signal strength of all the Higgs cascade searches. Note that for the Higgs cascades modes we use the projected sensitivity for L = 3000 fb −1 of data while for the conventional searches we use the best current limit. We cut off the x-axis at µ <37 fb −1 Curr. Lim. = 1 since points to the right of that are already excluded. Right: Same as the left panel, but the x-axis shows the best signal strength of the conventional channels including those arising via the direct production of the heavy non SM-like CP-even and CP-odd states H and A. In Sec. 4, we entertain the scenario that the LHC collaborations will be able to improve their sensitivities by one order of magnitude for the Higgs cascade decays compared to our projections (i.e. to µ 3000 fb −1 Proj = 0.1 on the y-axis) and two orders of magnitude compared to current conventional limits (i.e. to µ <37 fb −1 Curr. Lim. = 10 −2 on the x-axis). Then, all points except those in the bottom-left quadrangle bounded by the dash-dotted lines may be probed at the LHC with 3000 fb −1 of data. This quadrangle encloses only ≈ 10 % of the points shown -thus, such an improvement would allow future runs of the LHC to cover almost all (≈ 90 %) of the phenomenologically viable NMSSM parameter space containing additional Higgs bosons with masses below 1 TeV. Note that the scales of the x-axes differ between the panels. considering the detectability of the lighter states h and a via conventional searches, cf. the left panel, we find that for parameter points where one class of searches becomes ineffective, the other one usually fares well. If the lighter states h and a evade constraints from conventional Higgs searches, they are usually quite singlet-like such that their production cross section at the LHC is suppressed. However, mostly singlet-like light states are readily produced via Higgs cascade decays: For example, (Φ 1 → h 125 Φ 2 ) decays, where Φ i stands for a non SM-like Higgs mass eigenstate, are mostly controlled by the coupling λ if Φ 2 is singlet-like, and Φ 1 doublet-like. Since λ takes large values λ ∼ 0.65 in the alignment limit, the corresponding branching ratios are large such that searches utilizing Higgs cascades remain sensitive. From the right panel of Fig. 8, we see that direct searches for the heavy states H and A provide an additional handle for the Higgs cascade decays. Combining Higgs cascade decays with (conventional) direct searches for all the NMSSM Higgs bosons, the entire parameter space of the NMSSM with heavy Higgs bosons H and A lighter than ∼ 1 TeV is at most 2 − 3 orders of magnitude below current limits. Note the different scale for the x-axes between the left and right panels in Fig. 8.
In Fig. 8 we used our projected sensitivity for the 13 TeV LHC with 3000 fb −1 of data for the Higgs cascade channels, while we used current limits based, depending on the channel, on at most 37 fb −1 of data for the conventional Higgs searches. The increased statistical power of the future 3000 fb −1 data set should allow the bounds in the conventional searches to improve by approximately one order of magnitude. This would allow ≈ 50 % of our parameter points with masses of the additional Higgs bosons below 1 TeV to be probed by the LHC. Note that all of these parameter points satisfy current constraints. Hence, by combining all search channels, the LHC can make significant progress towards complete coverage of the NMSSM parameter space.
In the remainder of this section, we entertain the scenario that the LHC collaborations will be able to improve the sensitivity of their searches in the Higgs cascade channels by one order of magnitude compared to our projections (i.e. µ 3000 fb −1 Proj = 0.1) and two orders of magnitude compared to current limits in the conventional channels (i.e. µ ,37 fb −1 Curr.Lim. = 0.01). For the Higgs cascade decay based searches, such improvements could be realized by a combination of better rejection of reducible backgrounds and reduced systematic uncertainty of the remaining backgrounds. Note that improvements of comparable size have been demonstrated by both the ATLAS and the CMS collaboration in the past when updating analyses with increased statistics. For the conventional searches, we estimate that increasing the luminosity from the current L = O(30) fb −1 of data to L = 3000 fb −1 could yield one order of magnitude better sensitivity as discussed above, while another order of magnitude of improvement may be possible by improved background rejection/systematics and search strategies. While this is an optimistic scenario, it presents a clear target for the experimental collaborations which would allow the LHC to probe almost all of the remaining phenomenologically viable NMSSM parameter space featuring additional Higgs bosons with masses below 1 TeV.
In Figs. 9 and 10 we show the coverage of the parameter space when combining conventional searches with searches utilizing Higgs cascades under the assumptions that the sensitivity of the searches will be improved to µ 3000 fb −1 Proj.
(Higgs cascades) = 0.1 and Curr. Lim. (conventional) = 0.01 as discussed in the previous paragraph. These figures demonstrate that such a combination will allow the LHC to probe most of the NMSSM parameter space where at least one of the heavy Higgs states H or A has a mass below 1 TeV. The left panel of Fig. 9 presents the fraction of points scanned in each bin which will be probed at the LHC with 3000 fb −1 of data in the m A vs. (P S A ) 2 plane. The left panel of Fig. 9 and Fig. 10   region in the left panel of Fig. 9); in this region of parameter space a retains relatively large production cross sections at the LHC such that it can be searched for with conventional Higgs searches. The region with the least sensitivity in the left panel of Fig. 9 (heavy m A with sizable mixing (P S A ) 2 , colored pale yellow) is correlated with the heavy m h region in the right panel, particularly clustered around small (P S h ) 2 1, implying a dominantly doublet-like m h and a mostly singlet-like m H . Further, it can be seen from the left panel of Fig. 10 that this region corresponds to m a 350 GeV, where because of sizable mixing (as seen from the left panel of Fig. 9), the a is expected to have large branching ratios into pairs of top quarks, degrading search sensitivities. We again note the hard cut-off in the sensitivity for a at the top threshold visible in the lower triangle in the left panel of Fig. 10. As noted previously, while suppressed by alignment, m h can decay to pairs of gauge bosons, whereas such decays are forbidden for the CP-odd states at tree-level.
In the left panel of Fig. 10 we can most clearly see the complementarity between conventional searches and the Higgs cascade channels: If the non SM-like CP-even states (h and H) and the CP-odd states (a and A) are approximately mass degenerate, i.e. close to the diagonal, conventional searches are most powerful since all of the Higgs bosons could be heavily mixed and thus be copiously produced and decay to pairs of SM particles. On   the other hand, if the Higgs bosons are not mass degenerate, the Higgs cascade channels are most powerful, as we have already seen in Fig. 6. In the right panel of Fig. 10 we shows these sensitivities in the conventional m A -tan β plane for comparison with results generically presented in the MSSM. We see that the possibility of Higgs cascade decays in the NMSSM will allow the LHC to probe the low tan β region up to m A ∼ 1 TeV.
We note again that the distribution of points portrayed in Figs. 9 and 10 does not reflect the actual density of points scanned. As mentioned in section 2.3, the viable NMSSM parameter space we have analyzed corresponds to predominantly doublet-like m A and singlet-like m a , with tan β heavily clustered around ∼ 2.5. Combining conventional and Higgs cascade search channels, we expect ≈ 50 % of the points from our NMSSM parameter scan consistent with h 125 phenomenology, current direct search constraints, and featuring spectra with the additional Higgs bosons min(m A , m H ) ≤ 1 TeV to be probed by the LHC with 3000 fb −1 of data. Under the optimistic assumption that the LHC collaborations are able to improve their reach in the Higgs cascade channels by one order of magnitude with respect to our projections, and two order of magnitude in the conventional search channels with respect to current limits, almost all (≈ 90 %) of the points could be probed at the future LHC.

Conclusions
In this work, we have studied the collider phenomenology of the Z 3 -invariant NMSSM. We have focused on the Higgs and neutralino sector of the model, which is usually described in terms of the parameters appearing in the scalar potential. However, in order to be compatible with the observed Higgs phenomenology, the model must contain a 125 GeV Higgs mass eigenstate with SM-like couplings. This leads to strong correlations between the physical parameters, in particular the masses of the additional Higgs bosons and their supersymmetric partners, which are part of the neutralino sector. We have demonstrated that the Higgs and neutralino sector of the NMSSM can be effectively described by four physically intuitive parameters: the physical masses of the two CP-odd Higgs bosons, the mixing angle in the CP-odd Higgs sector, and tan β, all of which are quite transparently connected to the couplings of the physical Higgs and neutralino states. This reduction in parameters due to h 125 phenomenology, and the induced correlation in the physical masses and couplings, makes the NMSSM much more tractable than previously thought. We stress that we verified our conclusions with intensive numerics using NMSSMTools. Without implementing alignment a priori as an input to our scans, the correlated parameters and masses were obtained as an output.
Most search efforts for an extended Higgs sector at the LHC have been focused on direct searches, looking for resonant decays of heavy scalars/pseudo-scalars into SM particles as is generically predicted in supersymmetric models like the MSSM or generic 2HDMs. While such strategies have led to the exclusion of large mass ranges for the heavy doublet-like non SM-like states for large values of tan β, the low tan β region is challenging to probe for masses of the non SM-like Higgs bosons above ∼ 350 GeV where decays into pairs of top quarks are kinematically allowed. The presence of the singlet sector in the NMSSM allows for Higgs cascade decays, where a heavy Higgs boson decays into two lighter Higgs bosons or one light Higgs and a Z boson. As has been previously discussed e.g. in Refs. [18,21,22], these Higgs cascades can play an important role in the phenomenology of the NMSSM and provide a promising means to probe the low tan β regime. Of such decays, the branching ratios into pairs of SM-like Higgs boson or one SM-like Higgs and a Z boson are suppressed by the presence of the SM-like 125 GeV state observed at the LHC. Thus, the most relevant Higgs cascade modes are those into one SM-like Higgs and one non SM-like (light) Higgs boson, or into a Z boson and one non SM-like (light) Higgs boson; the corresponding branching ratios are not suppressed and the presence of the SM-like Higgs or the Z boson allows for the tagging of such processes at the LHC. If the additional Higgs bosons decay dominantly into pairs of SM particles, such Higgs cascades lead to Higgs+visible and Z+visible final states. If the dominant decay mode is into neutralinos, then the Higgs cascades lead to mono-Higgs and mono-Z final states.
Previously, the potential of such Higgs cascade decays modes for probing the NMSSM parameter space has been studied on a channel-by-channel basis in the literature. Here, we have provided a systematic comparison of the different Higgs cascade modes, gaining analytical understanding of the phenomenology due to our re-parameterization of the NMSSM parameters in terms of the physical parameters of the CP-odd sector. Most importantly, we have demonstrated that it is crucial to use as many different final states arising through Higgs cascades as possible, since no single class of final states dominates throughout the NMSSM parameter space. Further, we note that Higgs cascade modes may play a crucial role for differentiating models, e.g. the MSSM from the NMSSM, if additional Higgs bosons are discovered at the LHC. Higgs cascade decays usually involve the singlet-like states characteristic of the NMSSM, which are challenging to probe via conventional Higgs searches at the LHC since their production cross sections are suppressed with respect to doublet-like Higgs bosons.
In closing, we have demonstrated that the combination of Higgs cascade searches with conventional strategies to search for additional Higgs bosons via their decay into pairs of SM particles will allow ≈ 50% of the phenomenologically viable NMSSM parameter space with masses 1 TeV to be probed by the upcoming runs of the LHC. Under the optimistic assumption that the LHC collaborations are able to improve their reach in the Higgs cascade channels by one order of magnitude over our projections, and in the conventional search channels by two orders of magnitude with respect to current limits, ≈ 90 % of this NMSSM parameter space may be accessible to the LHC. While this is an optimistic scenario, it sets a target for the sensitivity required to probe most of the remaining interesting parameter space of the NMSSM.  Table 2. Direct Higgs searches at the LHC used for this work. The second column indicates the NMSSM Higgs bosons which can take the place of the generic scalar Φ i in the first column, recall that H/h (A/a) are the heavier/lighter non SM-like CP-even (CP-odd) states and h 125 is the observed SM-like 125 GeV Higgs boson. In the last row, the second column indicates possible pairs of (Φ i , Φ j ) in the corresponding process in the first column.

B Benchmark Points
In Tab. 3 we present the NMSSM parameters and mass spectra, in Tab Table 3. NMSSM parameters and mass spectra for our benchmark points categorized according to max µ 3000 fb −1 Proj.
(Higgs cascades) , BP1: Mono-Higgs, BP2: Higgs+visible, BP3: Mono-Z, and BP: Z+visibles. The first block from the top shows the parameters used as input parameters in NMSSMTools {λ, κ, tan β, µ, A λ , A κ , M Q3 } where the first 6 parameters are those appearing in the scalar potential, cf. Eq. (2.3), and M Q3 = M U3 is the stop mass parameter which controls the radiative corrections to the scalar mass matrices. For the convenience of the reader we also record the value of M A , defined in Eq. (2.18). As noted in Section 2.3, the remaining parameters are fixed to M 1 = M 2 = 1 TeV, M 3 = 2 TeV, A t = µ cot β, A b = µ tan β, and all sfermion mass parameters (except M Q3 = M U3 ) are fixed to 3 TeV. The second block shows the mass spectrum of the Higgs sector and the third block values of the singlet components of the non SM-like Higgs bosons. In particular, these blocks contain the masses of the CP-odd states a and A and the mixing angle in the CP-odd sector P S A . Recall that these three quantities were used in our physical reparameterization of the NMSSM (cf. Section 2.2) together with the value of tan β, which plays a minor role, and the values of λ and κ, which are approximately fixed by alignment. In the fourth block we record the masses of the three lightest neutralinos. Since we set the bino and wino mass parameters to M 1 = M 2 = 1 TeV, the two heaviest neutralinos χ 4 and χ 5 are bino-and wino-like with masses m χ4 ≈ m χ5 ≈ 1 TeV, while the three lightest neutralinos, χ 1 , χ 2 , and χ 3 , are Higgsinoand singlino-like.  (gg → H → Za → + − χ 1 χ 1 ) 1.73 -2.14 -   Regarding the branching ratios important for Higgs cascade decays, we first note that the branching ratio of heavy Higgs bosons into pairs of SM-like Higgs bosons or a SMlike Higgs and a Z boson is suppressed due to the proximity to alignment as discussed in Section 3, see also Refs. [18,21,41]. For all benchmark points, we find Additionally we note that in agreement with our discussions in Sec. 3.1, branching ratios of the heavy non-SM like doublets into either h 125 or Z and an additional singlet like state are generally comparable. This leads to multiple channels that may be probed at the LHC for each benchmark point, as discussed in detail below.

BP1 -Mono-Higgs
This benchmark point features a Higgs spectrum with comparable masses of the singletlike states a and h, m a = 205 GeV and m h = 165 GeV. The heavier states A and H are mostly composed of A NSM and H NSM , respectively, and are approximately mass degenerate with m A ≈ m H ≈ 650 GeV. The Higgsino mass parameter has a value of µ = −193 GeV, and κ = −0.281, leading to 2|κ|/λ = 0.93. Thus, the lightest neutralino χ 1 is mostly singlino-like but has sizable Higgsino components. Its mass is m χ 1 = 102 GeV, allowing for (a → χ 1 χ 1 ) decays but not for (h → χ 1 χ 1 ) decays. The second-lightest neutralino χ 2 is dominantly Higgsino-like with a mass of m χ 2 = 212 GeV ≈ |µ|, while χ 3 is mostly Higgsino-like but has a sizable singlino component and a mass of m χ 3 = 292 GeV.
Due to their singlet-like nature, the direct production cross sections of a and h are much smaller than those of a SM Higgs boson of the same mass, rendering them beyond the reach of conventional search channels at the LHC which rely on direct production of a or h. The dominant decay modes of h are into pairs of b-quarks, and, facilitated by its (small) doublet component, into pairs of W -bosons. The singlet-like pseudo-scalar on the other hand is kinematically allowed to decay into pairs of neutralinos (a → χ 1 χ 1 ). Because χ 1 has sizable singlino as well as Higgsino components, such decays proceed via both of the NMSSM's large couplings λ and κ, rendering the corresponding branching ratio large, BR(a → χ 1 χ 1 ) = 0.985.
The heavier (doublet-like) CP-even state H mostly decays into pairs of top quarks, neutralinos, and, most relevant for Higgs cascade channels, via (H → hh 125 ) and (H → Za). The cross section (gg → H → hh 125 ) is not large enough for it to be within reach of the Higgs+visible search modes. However, facilitated by the sizable branching ratios of (H → Za) and (a → χ 1 χ 1 ), this benchmark point is within the projected reach of mono-Z searches, µ 3000 fb −1 Proj.
(gg → H → Za → + − χ 1 χ 1 ) = 1.73. The heavier (doublet-like) CP-odd state A mostly decays into pairs of top quarks, neutralinos, and through the (A → h 125 a) channel. The sizable branching ratio of the latter decay mode, BR(A → Zh 125 ) = 0.304, together with the large branching ratio corresponding to (a → χ 1 χ 1 ) decays leads to a large projected signal strength in mono-Higgs searches via the corresponding decay chain, µ 3000 fb −1 Proj.
Neither H nor A have large branching ratios into pairs of SM states except into pairs of top quarks, rendering them very difficult to detect by conventional searches. Thus, the best chances to detect BP1 would be in mono-Z searches via (gg → H → Za) (the projected signal strength in this channel is µ 3000 fb −1 Proj. = 1.73) and particularly in mono-Higgs searches via (gg → A → h 125 a) with a projected signal strength µ 3000 fb −1 Proj. = 2.36.

BP2 -Higgs+visible
Benchmark point BP2 features a Higgs spectrum with a large split between the masses of the singlet-like states a and h, m a = 168 GeV and m h = 561 GeV. The heavier doubletlike states A and H are almost mass degenerate, m A ≈ m H ≈ 750 GeV. The Higgsino mass parameter takes much larger absolute value than for BP1, µ = −466 GeV. Further, κ also has a larger absolute value than for BP1, κ = 0.347, leading to 2|κ|/λ = 1.15. Thus, the two lightest neutralinos, χ 1 and χ 2 , are mostly Higgsino like and approximately mass degenerate, m χ 1 = 486 GeV and m χ 2 = 494 GeV. The third-lightest neutralino, χ 3 , is mostly composed of the singlino and has a mass of m χ 3 = 572 GeV. Note that because |2|κ|/λ − 1| is larger than for BP1, the Higgsino and singlino mass parameters are split further for BP2 than for BP1, leading to much smaller singlino-Higgsino mixing. Further, because of the relatively large masses of the neutralinos, none of the Higgs bosons are kinematically allowed to decay into pairs of neutralinos.
Similar to BP1, the large singlet components of a and h lead to direct production cross sections at the LHC much smaller than those of SM Higgs bosons of the same mass. Thus, they are out of reach of conventional search strategies. The dominant decay mode of the CP-even state h is into pairs of W -bosons and pairs of the much lighter singlet-like CP-odd state, BR(h → aa) = 0.740. The CP-odd state a decays mostly into pairs of b-quarks with a branching ratio of BR(a → bb) = 0.797. It also has a sizable branching ratio into τ -leptons, BR(a → τ + τ − ) = 0.0905.
The heavier (doublet-like) CP-even state H predominantly decays into pairs of top quarks. Because of the small value of tan β compared to BP1, (H → h 125 h) decays, which are mostly controlled by the (H SM H NSM H S ) coupling given in Tab. 6, are suppressed. The largest branching ratio of H relevant for Higgs cascade searches is BR(H → Za) = 0.0269. However, this branching ratio is not sufficiently large to put BP2 within reach of Z+visible searches where the projected signal strength is only µ 3000 fb −1 Proj.
Similar to H, the CP-odd doublet-like state A mostly decays into pairs of top quarks. The largest branching ratio relevant for Higgs cascade searches is BR(A → h 125 a) = 0.0188. This decay mode is mostly controlled by the (H SM A NSM A S ) coupling, which becomes largest for values of tan β = 1, but is suppressed for sgn(κ) = +1 see Tab. 6. Nonetheless, together with the large branching ratio of a into pairs of b-quarks, this branching ratio is sufficiently large to render BP2 within reach of Higgs+visible searches with a projected signal strength µ 3000 fb −1 Proj.
(gg → A → h 125 a → bbbb) = 1.22. Since both H and A decay dominantly into pairs of top quarks, BP2 is very challenging to discover with conventional search strategies. The most promising channel to discover this benchmark point at the LHC are Higgs+visible Higgs cascade searches, particularly in the bbbb final state with a projected signal strength of µ 3000 fb −1 Proj.
(gg → A → h 125 a → bbbb) = 1.22. If the sensitivity of the LHC can be improved by an order of magnitude over our projections, BP2 could also be probed via Z+visible searches through the (gg → H → Za) channel, µ 3000 fb −1 Proj.

BP3 -Mono-Z
Similar to BP2, the benchmark points BP3 also features a Higgs mass spectrum with a sizable split between the masses of the singlet-like states a and h, m a = 290 GeV and m h = 66.8 GeV. However, note that we have the inverted hierarchy for BP3 compared to BP2: for BP3 the CP-even state h is much lighter than the CP-odd state a, while for BP2 h is much lighter than a. The doublet like states are A and H are again approximately mass degenerate, m A = 696 GeV and m H = 678 GeV. The Higgsino mass parameter µ takes a moderate absolute value, µ = 251 GeV, and κ = −0.203 takes somewhat smaller absolute values than for BP1 and particularly BP2. The value of κ implies 2|κ|/λ = 0.640, implying a singlino mass parameter much smaller than the Higgsino masses. Correspondingly, we find that the lightest neutralino, χ 1 , is mostly singlino-like with a mass of m χ 1 = 97.6 GeV. The second-and third-lightest neutralino, χ 2 and χ 3 , are mostly Higgsino-like and have masses of m χ 2 = 248 GeV and m χ 3 = 323 GeV, respectively. The singlino-Higgsino mixing is smaller than for BP1, but χ 3 still has a singlino fraction of ∼ 0.3. Note that due to the mass spectra, the singlet-like CP-odd state a is allowed to decay into pairs of χ 1 's, while such decays are kinematically forbidden for the CP-even state h.
Due to their mostly-singlet nature, the direct production cross sections of a and h are much smaller than those of SM Higgs bosons of the same mass, rendering them difficult to detect with conventional search strategies. The CP-even state h dominantly decays into pairs of b-quarks with a branching ratio of BR(h → bb) = 0.909. The CP-odd state a is kinematically allowed to decay into pairs of neutralinos χ 1 . Similar to the case of BP1, the (somewhat smaller but still sizable) Higgsino components of χ 1 and its large singlino component renders the corresponding branching ratio large, BR(a → χ 1 χ 1 ) = 0.994, since such decays proceed through both of the NMSSM's large couplings λ and κ.
The heavier (doublet-like) CP-even state H has large branching ratios into pairs of top quarks, neutralinos, and into a Z boson and an a. The latter decay mode is particularly relevant for Higgs cascade searches, the corresponding branching ratio is BR(H → Za) = 0.249. Together with the large branching ratio of a into pairs of neutralinos χ 1 , we find a large projected signal strength in the corresponding mono-Z final state, µ 3000 fb −1 Proj.
(gg → H → Za → + − χ 1 χ 1 ) = 2.14. Further, although the branching ratio for (H → hh 125 ) decays is rather small, the large branching ratio of h into pairs of b-quarks renders BP2 in reach of Higgs+visible searches in the bbbb final state with projected signal strength of µ 3000 fb −1 Proj.
(gg → H → h 125 h → bbbb) = 1.64. The doublet-like CP-odd state A mostly decays into pairs of top quarks or a SM-like Higgs and the light pseudo-scalar state (A → h 125 a). The branching ratio of the latter decay is large, BR(A → h 125 a) = 0.212, and together with the large branching ratio of (a → χ 1 χ 1 ) decays puts this point within the projected reach of mono-Higgs searches with a signal strength of µ 3000 fb −1 Proj.
(gg → A → h 125 a → γγχ 1 χ 1 ) = 1.77. The branching ratio for (A → Zh) decays is small, BR(A → Zh) = 5.56 × 10 −2 , such that despite a rather large branching ratio of h into pairs of τ -leptons, BR(h → τ + τ − ) = 8.78 × 10 −2 this benchmark point would only be probed in Z+visible Higgs cascade searches via (gg → A → Zh) if the sensitivity of the search is improved by approximately one order of magnitude beyond our projections.
In summary, as was the case for BP1 and BP2, benchmark point BP3 is challenging to probe via conventional search strategies. This is because the heavier doublet-like states H and A have no large branching ratios into pairs of SM particles except into pairs of top quarks, and because the production cross sections of the lighter states a and h are suppressed due to their singlet-like nature. The most promising probe of BP3 are Z+visible Higgs cascade searches with a projected signal strength of µ 3000 fb −1 Proj.

BP4 -Z+visible
This benchmark points features somewhat lighter doublet-like Higgs states A and H than BP1−BP3. Furthermore, the mixing angle of the CP-odd Higgs bosons is sizable, P S A = 0.550. This leads to a larger split between the masses of the doublet-like states, m A = 533 GeV and m H = 460 GeV. The singlet-like states a and h have masses of m a = 157 GeV and m h = 298 GeV, respectively. The Higgsino mass parameter takes a comparatively small absolute value, µ = 141 GeV, while κ has a large absolute value, κ = −0.734. Thus, 2|κ|/λ = 2.20, corresponding to a singlino mass parameter much larger than the Higgsino mass parameter. Accordingly, we find that the two lightest neutralinos, χ 1 and χ 2 , are mostly Higgsino-like with masses of m χ 1 = 96.6 GeV and m χ 2 = 142 GeV, respectively. The third-lightest neutralino, χ 3 , is mostly singlet-like with a mass of m χ 3 = 376 GeV. These masses imply that the singlet-like CP-even state h can decay into pairs of χ 1 's, while such decays are kinematically forbidden for the singlet-like CP-odd state a.
The direct production cross sections of the singlet-like states a and h are again suppressed by their singlet-like nature. Note that due to its sizable A NSM component, the direct production cross section of a is only suppressed by a relatively small factor of σ(gg → a)/σ ggh (m a ) = 0.08 with respect to the gluon fusion production cross section of a SM Higgs boson of the same mass, σ ggh (m a ). However, this state is still very challenging to detect at the LHC because its couplings to pairs of W and Z bosons vanish at tree-level for a CP-odd state. Thus, its branching ratios into pairs of W and Z bosons as well as into pairs of photons are much reduced compared to a SM Higgs boson of the same mass. The most promising conventional search channels for a are via its dominant decay modes (a → bb) and (a → τ + τ − ), but the current limits in the corresponding search channels are relatively weak. The largest signal strength in conventional searches is in the (gg → a → τ + τ − ) search mode, µ <37 fb −1 Curr. Lim. (gg → a → τ + τ − ) = 0.117. Recall that this signal strength is calculated with respect to the current limit. We expect the limit in this channel to improve in future runs of the LHC such that BP4 may be in reach of conven-tional search channels via the (gg → a → τ + τ − ) channel. The CP-even singlet-like state h on the other hand still has only a small doublet fraction and thus a small production cross section at the LHC. Moreover, its dominant decay mode is into pairs of neutralinos with a corresponding branching ratio of BR(h → χ 1 χ 1 ) = 0.680, rendering it virtually impossible to discover with conventional search strategies.
The heavier (doublet-like) CP-even state H mostly decays into pairs of top quarks and via the (H → Za) channel. The latter decay has a large branching ratio of BR(H → Za) = 0.569. Together with the sizable branching ratio of a into pairs of τ -leptons, BP4 is rendered within reach of Z+visible Higgs cascade searches via the (gg → A → Zh) channel. The projected signal strength is µ 3000 fb (gg → A → h 125 a → γγbb) = 2.17. This benchmark point is difficult to probe in mono-Higgs and mono-Z channels despite the sizable branching fraction of h into pairs of neutralinos. The most promising mode is via the (gg → A → Zh) channel; however, the corresponding branching ratio of A, BR(A → Zh) = 0.06 is not sufficiently large to put the point in the projected reach of the LHC. The projected signal strength in this channel is µ 3000 fb −1 Proj.
(gg → A → Zh → + − χ 1 χ 1 ) = 0.189. In summary, although passing all current limits arising through conventional search strategies, BP4 may be probed by conventional searches if their sensitivity can be improved significantly in future runs of the LHC. The (gg → a → τ + τ − ) mode is the most promising decay channel with a current signal strength of µ <37 fb −1 Curr. Lim. (gg → a → τ + τ − ) = 0.117. Furthermore, this benchmark point can be probed by Higgs cascade searches. We find the largest projected signal strength in Z+visible searches, µ 3000 fb (gg → A → h 125 a → γγbb) = 2.17. Figure 11. Same as Fig. 6, but for the mono-Z (left) and mono-Higgs (right) final states arising through decays of the parent heavy Higgs into a pair of neutralinos where one of the neutralinos subsequently radiates off a Z or a Higgs, cf. diagram (c) in Fig. 3.