Two Higgs Doublets and a Complex Singlet: Disentangling the Decay Topologies and Associated Phenomenology

We present a systematic study of an extension of the Standard Model (SM) with two Higgs doublets and one complex singlet (2HDM+S). In order to gain analytical understanding of the parameter space, we re-parameterize the 27 parameters in the Lagrangian by quantities more closely related to physical observables: physical masses, mixing angles, trilinear and quadratic couplings, and vacuum expectation values. Embedding the 125\,GeV SM-like Higgs boson observed at the LHC places stringent constraints on the parameter space. In particular, the mixing of the SM-like interaction state with the remaining states is severely constrained, requiring approximate alignment without decoupling in the region of parameter space where the additional Higgs bosons are light enough to be accessible at the LHC. In contrast to 2HDM models, large branching ratios of the heavy Higgs bosons into two lighter Higgs bosons or a light Higgs and a $Z$ boson, so-called Higgs cascade decays, are ubiquitous in the 2HDM+S. Using currently available limits, future projections, and our own collider simulations, we show that combining different final states arising from Higgs cascades would allow to probe most of the interesting region of parameter space with Higgs boson masses up to 1\,TeV at the LHC with $L=3000\,{\rm fb}^{-1}$ of data.


Introduction
With the discovery of a Higgs boson at the Large Hadron Collider (LHC) in 2012 [1,2] with mass m h 125 ≈ 125 GeV [3,4] and couplings compatible with those of a Standard Model (SM) Higgs [4][5][6][7][8], the final ingredient of the SM has been found. However, despite the SM holding up to all laboratory test it has been subjected to so far, it fails to explain the behavior of the Universe at large scales. In particular, the SM does not contain a suitable candidate for the observed Dark Matter and fails to explain the matter-antimatter asymmetry. Beyond such phenomenological problems, the SM also suffers from some issues more theoretical in nature, e.g. the hierarchy problem or the lack of explanation of the SM's flavor structure. During the past 50 years, a multitude of Beyond the Standard Model (BSM) particle physics models have been developed to address the aforementioned problems. The vast majority of BSM models, in particular those valid to energy scales much larger than the electroweak scale, feature a scalar sector extended beyond the SM's one SU (2)-doublet Higgs. For example, consistent supersymmetric extensions of the SM require a Higgs sector containing at least two Higgs doublets. Furthermore, while the parameters of the Higgs sector are those least well measured in the SM, many of the SM's shortcomings are intimately related to the scalar sector e.g. the hierarchy problem, the matter-antimatter asymmetry, and the flavor structure.
There are two avenues for studying BSM physics: One can either take a top down approach, starting from well-motivated SM extensions valid at energy scales much larger than the electroweak scale, or, one can take a bottom up approach by parameterizing our ignorance of high-scale physics by writing down the most general form of BSM models at the weak scale. In the case of extensions of the SM's Higgs sector, this may be exemplified by the well studied case of BSM models containing two Higgs doublets. One can either choose to study more complete models containing such a Higgs sector, e.g. the Minimal Supersymmetric Standard Model (MSSM), or, one can take a more model independent approach and study the most general form of a two Higgs doublet Model (2HDM).
From a bottom up perspective, extensions of 2HDMs with an additional singlet may be motivated from the well known fact that they facilitate baryogenesis [9][10][11] to explain the matter-antimatter asymmetry, from flavor considerations [12,13] and from being useful for constructing Dark Matter models [14]. In addition, well motivated top down BSM models exist containing such a scalar sector, for example the Next-to-Minimal Supersymmetric Standard Model (NMSSM) [15,16]. In this work, we take the bottom up approach to study extensions of the SM's scalar sector with one additional doublet and one complex SM gauge singlet S, which we refer to as 2HDM+S in the following.
While a large body of literature exists on general 2HDM models, see e.g. Refs. [17][18][19][20][21] for recent works, no systematic study of 2HDM+S models exists in the literature to the best of our knowledge. Some aspects of 2HDM+S models have been discussed in Ref. [22] and in the appendix of Ref. [23]. The extension of 2HDMs with a real singlet has been discussed in Refs. [24][25][26][27][28]. In this work, we provide a first systematic study of the 2HDM+S parameter space. The physical Higgs sector contains 5 real neutral Higgs bosons and one charged Higgs boson after electroweak symmetry breaking. Assuming CP-conservation, the 5 neutral Higgs bosons can be subdivided into 3 CP-even states and 2 CP-odd states. In order to be compatible with the observed phenomenology, one of the CP-even states must have couplings to pairs of SM particles similar to those of a SM Higgs with a mass of 125 GeV, such that it can be identified with the 125 GeV Higgs boson observed at the LHC. As we will see, the presence of this SM-like state has important implications for the behavior of the remaining Higgs bosons.
Extended Higgs sectors such as the 2HDM+S can be tested at the LHC in two complimentary ways: one can study the behavior of the observed approximately SM-like h 125 , constraining its possible mixing with additional Higgs bosons, or, one may directly search for the remaining Higgs bosons in the model. The first way has thus far constrained a number of couplings of h 125 to SM particles to be within O(10 %) of the SM values. In most cases, future data expected to be collected at the LHC during the next 20 years will allow only for marginal improvements on these bounds [29]. Note however, that a Higgs factory such as the International Linear Collider would allow for precision tests of the SMlike Higgs boson, probing some of its couplings below the O(1 %) level [30]. This would yield powerful constraints on any BSM Higgs sector.
Such indirect searches are complimented by direct searches for the additional Higgs bosons as carried out at the LHC in a plenitude of final states. Most commonly, one searches for direct production of the additional Higgs bosons which then decay into pairs of SM particles. Such search strategies are similar to those employed in the hunt for the 125 GeV Higgs boson prior to its discovery. For a 2HDM Higgs sector, searches for additional Higgs bosons are complicated by three different issues: 1) While CP-even states are allowed to decay into pairs of SM vector bosons, the SM-like nature of h 125 suppresses the couplings of the additional states to pairs of vector bosons. 1 2) In addition to suppressing the branching ratio into pairs of SM vector bosons, the SM-like nature of h 125 also suppresses the coupling of heavy Higgs bosons to pairs of SM-like Higgs boson or a Z boson and a SM-like Higgs, reducing the power of the corresponding search channels. 3) For large regions of parameter space the decays of any Higgs boson with mass m Φ 350 GeV are dominated by decays into pairs of top quarks. Due to the interference of this signal with the QCD background [31][32][33][34][35][36][37][38], such decays are challenging to use for Higgs searches at the LHC with current search strategies.
All three of these complications are still present in models with enlarged Higgs sectors such as the 2HDM+S. In addition, the presence of the singlet, which can mix with the doublets, complicates conventional searches further since the production cross section of any Higgs state at the LHC is suppressed with a growing singlet fraction. However, the presence of additional Higgs bosons allows for new decay channels of the heavy Higgs bosons, so-called Higgs cascades, where the heavy Higgs bosons decay into two lighter Higgs states or a light Higgs and a Z boson [23,[39][40][41][42][43][44]. Such decays can also be present in the generic 2HDM without the singlet if the CP-even and CP-odd Higgs bosons have non-degenerate masses, although large mass splittings are difficult to achieve in consistent 2HDMs [45]. However, in this case the direct production cross section of the lighter state has no suppression from a singlet component, generally making direct searches for such a state more powerful. We briefly discuss the limiting case of decoupling the singlet from the doublets in Sec. 2.3.
As mentioned above, the branching ratios of any of the heavy Higgs bosons to pairs of SM-like h 125 's, pairs of vector bosons, or a Z boson and a h 125 , are suppressed by the SM-like nature of h 125 . Of the remaining decay modes with possibly large, O(10 %), branching ratios, those consisting of a SM-like Higgs boson or a Z boson and one additional non SM-like Higgs are most useful for searches at the LHC. This is because one can use the decay products of the SM-like Higgs or the Z boson with known mass and branching ratios to tag such events. Furthermore, final state search signatures will also contain the decay products of the non SM-like Higgs bosons: if they decay into visible states, e.g. a pair of SM particles, one searches for a resonance in the invariant mass spectrum of the decay products. However, the non SM-like Higgs bosons can also decay into additional new states, for example pairs of new stable neutral particle playing the role of a Dark Matter candidate. In this case, the decay products of the non SM-like Higgs bosons produced in the Higgs cascades would leave the detector without depositing energy, manifesting as missing transverse energy (E miss T ) in the detector, giving rise to so-called mono-Z and mono-Higgs signatures.
As shown in the context of the NMSSM in Refs. [23,43,44], the branching ratios of Higgs cascades are sizable in large parts of the parameter space, leading to observable signatures at the LHC. As we will see below, this is also the case in the general 2HDM+S.
The remainder of this work is organized as follows: In Sec. 2, we discuss the 2HDM+S. We focus on the mass spectrum, and the implications of the SM-like nature of the 125 GeV Higgs boson in Sec. 2.1, and the couplings between the Higgs bosons and the relevant model parameters for the Higgs cascades in Sec. 2.2. In Sec. 2.3 we discuss the production cross sections and decay modes of the 2HDM+S Higgs bosons, and in Sec. 2.5 we discuss the Higgs cascades in more detail. The more experimentally inclined reader may skip directly to Sec. 3 where we discuss the future reach of the LHC for the full 2HDM+S using Higgs cascades. In Sec. 3.2 we demonstrate that the 2HDM+S gives rise to cross sections within reach of the Higgs cascade searches at the LHC for large regions of parameter space. We reserve Sec. 4 for our conclusions. Appendices A, B, C, D, and E contain details about the mass matrices, trilinear couplings, quartic couplings, the collider simulation performed for estimating the reach of the mono-Z channel, and the mapping of the 2HDM+S parameters onto the NMSSM, respectively.

2HDM+S
The most general 2HDM scalar potential is given by [17] where Φ 1 , Φ 2 are SU (2) doublets with hypercharge Y = 1/2. The m 2 ij parameters have dimension mass squared, while the λ i are dimensionless. In this work, we study the CPconserving case, where all parameters can be chosen manifestly real. Extending the field content by a complex scalar SM gauge singlet S, the most general scalar potential for the additional terms is The parameter ξ has dimensions mass cubed, the parameters {µ Si , µ ij } have dimension mass, and the {λ i , λ i } are dimensionless. For the CP-conserving case, all parameters can again be chosen to be manifestly real. The most general scalar potential for a 2HDM+S model is then given by Note that the NMSSM has the same Higgs sector, although supersymmetry severely restricts the parameters. We provide a mapping of the 2HDM+S onto the Higgs sector of the NMSSM in Appendix E.
Out of the 29 parameters in the scalar potential only 27 are physical. One can perform a global SO(2) rotation in (Φ 1 , Φ 2 ) space [17], and a real shift of S. These transformations render one of the parameters in the doublet sector and one parameter in the singlet sector unphysical, respectively. In the following, we will use the basis where the shift of S has been used to eliminate the tadpole term ξ from the scalar potential. We do not further define the particular basis for the doublets yet, as such a choice is intimately linked to the choice of couplings to the SM fermions, discussed below.
One can furthermore use the minimization conditions For the purposes of this work, we will assume that the vevs lie along the neutral direction of the doublets and that the vacuum is CP-conserving. Note that these assumptions may be broken even when using the vevs {v 1 , v 2 , v S } as input parameters and choosing them manifestly real, as we will do in the following, see for example the discussion in Ref. [17]. Further, deeper minima, including CP-violating ones, may be present. An analysis of the vacuum structure of the 2HDM+S is beyond the scope of this work and we leave it for the future.
In addition to the field transformations discussed above, one can redefine each of the doublets as well as the singlet by a phase. Demanding the potential to remain manifestly CP-conserving under such transformations, i.e. all parameters remain real, constrains these transformation to U (1) ⊗ Z 2 in the doublet sector and Z 2 for the singlet. Although these transformations cannot be employed to absorb additional parameters, one can use them to enforce all vevs positive, {v 1 , v 2 , v S } ≥ 0.
As customary, we define The where {H SM , H NSM , H S } and {A NSM , A S } are the neutral CP-even and CP-odd real Higgs basis interaction states and G 0 (G ± ) is the neutral (charged) Goldstone mode. In this basis, of the states coming from the doublets, only H SM = √ 2v acquires a vev, and it is straightforward to work out the coupling of SM fermions to the Higgs basis states to e.g. study possible flavor changing neutral currents (FCNCs). Potentially dangerous FCNCs can be omitted if the groups of right-handed SM fermions with the same quantum numbers couple to only one of the doublets, respectively, as in the so-called Type I, Type II, flipped, and lepton specific 2HDM versions 3 [17]. Note that couplings of the fermions to the singlet S are forbidden by gauge invariance. Restricting ourselves to the case where the values of the Yukawa matrices are chosen such that the observed SM fermion masses are obtained, the couplings of pairs of SM particles to the Higgs basis states can be written as (2.10) H NSM (f 1 , f 2 , VV) = (g SM / tan β, −g SM tan β, 0) , (2.11)

12)
A NSM (f 1 , f 2 , VV) = (g SM / tan β, g SM tan β, 0) , (2.13) where "f 1 " ("f 2 ") stands for SM-fermions coupling to Φ 1 (Φ 2 ), "VV" for pairs of W or Z gauge bosons, and g SM is the coupling of a SM Higgs boson to such particles. Note that CP-odd states couple to pseudoscalar fermion bilinears (f γ 5 f ), instead the CP-even states couple to the scalar bilinear (f f ). For concreteness, in the following we assume a Type II Yukawa structure, where Φ ≡ iσ 2 Φ * , the Y i are the 3 × 3 Yukawa matrices, and the left-handed quarks Q and leptons L as well as the right-handed up-type (down-type) quarks u R (d R ) and the right-handed leptons e R should be understood as vectors in generation space. We define all parameters of the 2HDM+S scalar potential, in particular the vevs, in the interaction basis as determined by the Yukawa structure [cf. the discussion below Eq. (2. 3)], and where the singlet is shifted such that ξ = 0. In particular, this means that we keep all parameters of the doublet sector in our expressions, but one should keep in mind that one of these parameters is unphysical due to the SO(2) symmetry in (Φ 1 , Φ 2 ) space. We stress that the scalar potential of our model allows more general Yukawa structures than the one chosen in Eq. (2.15), including ones leading to tree-level FCNCs. Consistency with a desired Yukawa structure, e.g. Type I or Type II, can be ensured by imposing an ad-hoc Z 2 symmetry, as is done in generic 2HDM's. For example, the Type II structure in Eq. (2.15) and its radiative stability can be ensured by a Z 2 symmetry under which Φ 1 → −Φ 1 , u R → −u R , and all other fields transform trivially. 4 In this work, we remain agnostic about the mechanism which ensures the desired Yukawa structure and its radiative stability. Hence, we will not impose additional symmetries on the model. The inclined reader may choose their favorite mechanism; imposing the corresponding restrictions on the 2HDM+S parameter space is straightforward. Note also that our results in the remainder of this paper will in general hold for a different Yukawa structure. Some quantitative details may change because of the change of the Yukawa enhancement/suppression of the fermion couplings. However, such modifications will be small since we mostly consider the low tan β = O(1) regime in this work.

Higgs Mass Eigenstates and Alignment
The mass eigenstates are obtained from the diagonalization of the squared-mass matrix for the Higgs basis states. The neutral and charged Goldstone modes G 0 and G ± are massless by construction and do not mix with the other interaction states. In the following, we remove the Goldstone modes from the theory by choosing the unitary gauge. Furthermore, there is no mixing between CP-even and CP-odd states in the CP-conserving 2HDM+S.
We denote the three CP-even mass eigenstates where S j h i with j = {SM, NSM, S} denotes the components of the mass eigenstates in terms of the interaction basis. Likewise, we denote the two CP-odd mass eigenstates where again m A > m a , and where the components are similarly denoted by P j a i . The S j h i (P j a i ) are obtained from diagonalizing the (symmetric) squared mass matrix for the CP-even (CP-odd) Higgs bosons, M 2 S (M 2 P ). The values of the entries of the mass matrices and the mass of the charged Higgs boson are recorded in Appendix A.
The observation of a m h 125 ≈ 125 GeV mass eigenstate with couplings to SM particles compatible with that of a SM Higgs boson at the LHC implies that our model must contain a mass eigenstate with Hence we see that there are two possibilities to achieve alignment: either the left hand sides of Eqs. (2.21) and (2.22) go to zero, or, the right hand sides become large while the left hand sides remain non-zero and sizable. The latter possibility is the so-called decoupling limit, corresponding to implying {m H , m h } m h 125 . The first option is the so-called alignment without decoupling limit, and is of particular interest for LHC phenomenology. This is because the additional CP-even mass eigenstates H and h are not necessarily much heavier than h 125 and hence may be directly accessible at the LHC [18,19,23,[52][53][54][55]. In this case, must be satisfied in order to ensure approximate alignment. Neglecting radiative corrections, perfect alignment is achieved for 5 where the first condition ensures M 2 S,12 = 0 and the second condition M 2 S,13 = 0. Here and in the following we employ a shorthand notation,

Couplings and Parameters
The trilinear couplings of the Higgs bosons can be obtained from the scalar potential as (2.28) We first note that CP-conservation forbids couplings such as g h i h j a k or g a i a j a k . In the following, we discuss the remaining trilinear couplings in the extended Higgs basis. From these, the couplings for the mass eigenstates can be obtained via (2.29) Some of the couplings allowed by CP-conservation, vanish due to gauge invariance. For example, the coupling is identically zero since it would arise only through terms proportional to S(Φ i Φ j ), i, j = {1, 2}, in the scalar potential which are forbidden by U (1) Y invariance. A number of the trilinear couplings are proportional to entries in the mass matrices: Categorizing the Higgs basis states as SM-like (H SM ), NSM-like (H NSM , A NSM , H ± ), and singlet-like (H S , A S ), we find that couplings involving only SM-like states are proportional to the corresponding diagonal entry of the CP-even mass matrix Couplings involving two SM-like states are proportional to the corresponding entry in the CP-even squared mass matrix mixing such states, Couplings involving one SM-like state, one NSM-like state, and one singlet-like state are proportional to the singlet-doublet mixing entries in one of the mass matrices This implies that the couplings of two SM-like states to one NSM-like or one singlet-like state are suppressed in the alignment limit. Indeed, for perfect alignment one finds and hence the couplings of the mass eigenstates, are suppressed in proximity to alignment, and vanish for perfect alignment.
Besides the trilinear couplings proportional to the entries of the mass matrices, a few of the couplings are identical or can be written as linear combinations of entries of the mass matrices. Accounting for such degeneracies, we find that there are 10 independent trilinear couplings, For completeness we provide a list of the trilinear couplings and their values not a priori forbidden by CP or charge conservation in Table 1 located in Appendix B, which match the results obtained in Ref. [23].
It is straightforward to show from the scalar potential that there are 4 independent quartic couplings between the interaction states which cannot be written as linear combinations of the entries of the mass matrices or the trilinear couplings, (2.39) We list the values of these couplings in Appendix C. All together, the 9 entries of the mass matrices, which we parameterize by the 5 physical masses the 10 independent linear couplings listed in Eq. 2.38, and the 4 independent quartic couplings given in Eq. 2.39 yield a set of 27 independent parameters. In the remainder of this work we will describe the 2HDM+S parameter space in terms of these parameters since they are more closely related to physical observables instead of the 27 parameters in the scalar potential. That said, we note that of these 27 parameters, only a subset are relevant for the analysis of the heavy Higgs decay topologies, especially given the constraints due to the SM-like nature of h 125 . In particular, the quartics play no role (however, these may be interesting when considering di-Higgs production), and as we shall see, only a few of the trilinear couplings listed in Eq. (2.38) will be needed.

Production Cross Section and Partial Decay Widths
For moderate/low values of tan β = O(1), the production cross section for non-standard Higgs bosons at the LHC is dominated by gluon fusion. For larger values of tan β, bbH/bbA associated production may become relevant. Note that in proximity to the alignment limit the vector boson fusion production cross sections of additional Higgs boson is suppressed since only H SM couples to pairs of vector bosons.
To first approximation, the gluon fusion production cross section can be parameterized by the contribution from top quarks in the loop. At one-loop order, the ratio of the gluon fusion production cross sections of a CP-even Higgs boson h i and a CP-odd Higgs boson a i is then given by [56] where σ ggh (m) is the gluon fusion production cross section of a SM Higgs boson with mass m, τ a i ≡ (m a i /2m t ) 2 , m h i is the mass of h i , m a i is the mass of a i , m t is the top quark mass, and As long as the narrow width approximation is valid, the production cross section of any final state arising from gluon fusion production of a Higgs boson is given by the product of the gluon fusion production cross section of the Higgs boson of interest, and its branching ratio into the relevant final state. The relevant branching ratio is given in terms of the partial decay widths by where the width Γ Φ i of the state Φ i is obtained by summing over all partial widths. In the remainder of this section, we discuss partial widths of the most relevant decay modes of the additional Higgs bosons.
The partial decay width for a CP-even mass eigenstate h i into pairs of vector bosons is given by whereas such decays are forbidden for CP-odd mass eigenstates at tree level. Decays into pairs of SM fermions are allowed for both CP-even and CP-odd Higgs bosons. The corresponding partial widths can be written , for down − type quarks and leptons f , Here and in the following, we use C J Φ to refer to the mixing angles for both the CP-even and the CP-odd mass eigenstates Φ with mass m Φ ; for example, the corresponding partial width is given by In addition to the above, the extended Higgs sector of the 2HDM+S model allows potentially large decay widths of the heaviest Higgs bosons into pairs of lighter Higgs bosons or a Higgs and a Z boson. CP-conservation allows decays into pairs of Higgs bosons only if they are of the type (h i → h j h k ), (h i → a j a k ), or (a i → h j a k ), while decays into a Z and a Higgs boson must be of the type (h i → Za j ) or (a i → Zh j ). The corresponding partial widths are given by where the trilinear couplings between the Higgs mass eigenstates are given in Eq. (2.29).

Branching Ratios: SM Fermions, Misalignment and the 2HDM Limit
Considering the decays of the heavy Higgs into pairs of SM particles, it is worth repeating that in the alignment limit only the SM-like mass eigenstate h 125 decays into pairs of vector bosons, while the corresponding branching ratio for the non SM-like mass eigenstates such as (H → ZZ) vanish. The dominant decay mode into SM particles will thus be into the pair of kinematically accessible SM fermions m Φ > 2m f with the largest tan β suppressed (enhanced) Yukawa coupling where γ = −1 (γ = +1) for up-type (down-type) fermions accounts for the tan β suppression (enhancement) of the Yukawa couplings. For moderate values of tan β, the dominant decay mode will thus be into pairs of top quarks if m Φ > 2m t . However, decays into bottom quarks will dominate over those into top quarks even for m Φ > 2m t if tan β m t /m b ≈ 6.4. In the following, we identify the regions of parameter space where other decay modes may dominate over decays into pairs of SM fermions.
First, we note that if kinematically accessible, any of the Higgs bosons can decay into possible additional singlet fermions χ i χ j . Ignoring kinematic factors, such decays would compete with decays into SM fermions if  2v, which would significantly reduce the production cross section of such a Higgs boson at the LHC.
If kinematically accessible, we find from Eqs. (2.53) and (2.49), that decays of a non SM-like Higgs boson Φ i into a Z and a lighter Higgs state Φ j will compete with decays into pairs of SM fermions if From the form of the trilinear couplings tabulated in Tab. 1, we see that respectively, where we have defined and assumed perfect alignment, implying S NSM H = S S h etc. Since the couplings g H SM H NSM H S and g H SM A NSM A S are proportional to the entries of the mass matrices, these couplings can be rewritten as Hence, such decays will compete with ff if Since there is no additional suppression of this due to either the NSM or the S component of h 125 , it can be achieved much more easily, even for small values of C S Φ i and m Φ i ∼ 500 GeV.
The trilinear couplings relevant for decays of the heavy Higgs into daughter Higgs bosons without h 125 , i.e. (H → hh), (H → aa), and (A → ha), are not necessarily suppressed by small mixing angles, and the corresponding branching ratios can be large. However, these branching ratios are mostly governed by parameters in the scalar potential which are not related to either the mixing angles or physical masses. Hence, these branching ratios can essentially be taken as free parameters in a generic 2HDM+S.
Observe that while (Φ i → Φ j Φ k ) cascade decays are controlled by the trilinear couplings g Φ i Φ j Φ k which are a function of the mixing angles, masses, and the Higgs basis trilinear couplings between the Higgs basis states listed in Eq. are not suppressed, and if kinematically accessible may be the most relevant cascade decays for probing such models at the LHC. This is because if compared to decays into pairs of new light Higgs boson, such as (H → hh) or (A → ha), they offer a final state with known mass and branching ratios into pairs of SM particles, which can be employed to tag such signatures. Finally, we would like to comment on the limit where the singlet states decouple from the doublet-like Higgs states. Such a situation is achieved if the mixing angles satisfy 69) where the first (second) option corresponds to the heavier (lighter) non SM-like states H and A (h and a) to be comprised of the singlets H S and A S . The first option corresponds to the 2HDM limit of the 2HDM+S.
The latter case can be achieved by a conspiracy of parameters yielding {M 2 S,23 , M 2 P,12 } → 0 in addition to the second alignment condition implying M 2 S,13 → 0. Due to the vanishing mixing of the singlet-like and doublet-like states, the singlet-like states do not couple to SM particles and hence can neither be directly produced at the LHC nor decay, since all couplings to other particles vanish. Recall that in the alignment limit the g H SM H SM H S coupling vanishes as well such that even (h → h 125 h 125 ) decays are forbidden. Thus, the singlet-like states h and a can play the role of DM candidates [57][58][59][60][61]. If kinematically allowed, the singlets may be produced at the LHC in the decays of h 125 . The couplings corresponding to (h 125 → hh) and (h 125 → aa) decays are g H SM H S H S and g H SM A S A S to leading order, respectively. These couplings are unsuppressed, and hence these decays are only constrained by the invisible decay width of h 125 [4]. If at least one of the doublet-like states H and A is sufficiently light such that it can readily be produced at the LHC, are the only other production channels of the singlet-like states at the LHC.
In addition, if the doublet-like states A and H have sufficiently different masses and are light enough such that the heavier state is readily produced at the LHC, depending on the mass ordering, one of the Higgs cascade channels remain effective, cf. the discussion in the previous paragraph. All other Higgs cascades modes are kinematically forbidden, suppressed by alignment, or suppressed by the vanishing production cross section of the singlet-like states at the LHC.

Comparison of Higgs Cascades
In the previous section we compared the decay widths corresponding to the Higgs cascades to the decays into pairs of SM particles. In the following, we will compare the final state cross sections for the different Higgs cascades to each other, i.e. we identify which Higgs cascade mode is the most relevant for different regions of the 2HDM+S parameter space.
In particular, we want to compare cascade decays into lighter Higgs bosons to cascade decays into a Z and a light Higgs boson. The ratio of the production cross sections is given by in the narrow width approximation, assuming that the production cross section of the heavy Higgs bosons is dominated by gluon fusion.
In the rest of our analyses and results, for simplicity, we will assume perfect alignment for h 125 . While some misalignment is allowed by current experimental observations, h 125 can generically only be contaminated by a NSM fraction of few % and a singlet fraction of 20% (see e.g. Fig. 1 in Ref. [23]). In turn, this contamination is precisely the SM component the other Higgs bosons in the model can have. Since we have so far presented all results in terms of mixing angles and masses, the effect of misalignment on the various cross sections can be deduced by assuming a small SM component for the additional Higgs bosons. From Eqs. (2.44)-(2.54) it is easy to see that although the allowed misalignment quantitatively affects the final state cross sections, the qualitative behavior of the model will not be affected. We stress that the phenomenology of the extended Higgs sector of the 2HDM+S will thus primarily be dictated by the results we present in the alignment limit.
We first compare the relevance of the two Higgs cascades for the same parent state Φ i = Φ l . Using Eqs. (2.52) and (2.53) we find where we have implicitly assumed perfect alignment, of particular relevance for estimating the trilinear couplings for the mass eigenstates relevant for the denominators in the ratios above (see discussion in Sec. The ratios of the final state production cross sections for processes involving different parent Higgs boson are slightly more involved, since one has to take into account the ratio of the gluon fusion production cross sections as well as the ratio of the decay widths. From Eqs. (2.52), (2.53), and (2.74), we find in the alignment limit   78) also depend on the ratio of the gluon fusion production cross sections and the total decay widths of the parent heavy Higgs bosons. The total decay widths are well approximated by the summation of the partial widths listed in Eqs. (2.47)-(2.53) and in general are functions of the Higgs mass spectrum, the mixing angles, tan β, and additional free trilinear couplings between the Higgs bosons. Observe that since these free trilinear couplings can be sufficiently large such that the corresponding decay modes dominate the decay width, in principle they can dramatically affect the ratios in Eqs. (2.77) and (2.78). However, in certain limits the ratio of the total width can be well approximated rather simply, e.g. if the decay width of both A and H is dominated by the partial width into pairs of top quarks, we can approximate the ratio in the alignment limit for (2.80) In Figs. 1-3 we show the ratios of the final state cross sections for the cascade decays in the most relevant parameter planes, fixing the remaining parameters to fiducial values. Note that, as discussed previously, the complimentary ratios not shown can be obtained by swapping all parameters referring to the CP-even sector with the corresponding quantities in the CP-odd sector. We have chosen the mass scales of the involved particles in the range where the LHC would be most sensitive: heavy Higgs bosons with masses above ∼ 1 TeV are challenging because the gluon fusion production cross section drops sharply with mass, while much smaller masses of mostly doublet-like additional Higgs bosons are constrained by direct Higgs searches. shown in the right panel can be understood from the fact that, ignoring the dependence on the ratio of the total decay widths, but incorporating the effect of the mixing angles on the respective production cross sections, this ratio is proportional to (P NSM cf. Eq. (2.77).
In Fig. 2 we show the same ratios in the planes of the parameters which, apart from the mass spectrum, control the mass eigenstates couplings g h 125 Hh and g h 125 Aa relevant for   of the corresponding cascade topologies to decrease with growing mass of the parent Higgs boson, even though the relevant widths for the decay of the parent state increase with growing mass. In addition, the ratio shown in the right panel of Fig. 3 depends on the ratio of the total decay width of the involved parent Higgs bosons Γ H /Γ A . For the range of masses displayed, no strong effects are expected from the total decay widths. However, one finds more pronounced effects for different choices of the masses; in particular threshold effects, e.g. when m A or m H falls below 2m t ≈ 350 GeV. Then, the corresponding decay mode, which may have had a large branching ratio, becomes kinematically forbidden.
In summary, we observe from Figs. 1-3 that while there may well be a conspiracy of parameters such that one decay mode is severely suppressed compared to the others, for large parts of the parameter space the cross sections for Higgs cascades with a SM-like Higgs boson and a light additional state (gg → Φ i → h 125 Φ j ) are of the same order of magnitude as the cross sections for the cascades with a Z boson and an additional light Higgs in the final state (Φ i → ZΦ j ). Hence, one needs to study all possible decay topologies at the LHC in order to comprehensively cover the 2HDM+S parameter space.

Current Limits and Future Projections
We want to compare the cross section for the cascade decays in the 2HDM+S with the projected reach at the future LHC, assuming a collected luminosity of L = 3000 fb −1 . In the preceding section we discussed the production of pairs of scalars (Φ j Φ k ) or a Z boson and a scalar (ZΦ j ) from an intermediate scalar Φ i produced by gluon fusion (gg → Φ i ). However, at the LHC one would not observe such states directly, but only the decay products of the Higgs and Z bosons. Since the branching ratios into pairs of SM particles as well as the masses of the Z boson and the observed SM-like Higgs state h 125 are known, they can be used to tag such processes. On the other hand, neither the branching ratios nor the mass of the additional scalars a and h are known. Depending on the mass spectrum and matter content of the model, additional Higgs bosons may decay either in pairs of SM particles (visible), or, if present and kinematically allowed in pairs of new fermions χ which, if they are neutral and stable, would leave the detector without depositing energy (invisible) yielding missing transverse energy (E miss T ). In the following, we classify the final states arising from such Higgs cascade decays as • h 125 + invisible, or mono-Higgs, • Z + visible, • Z + invisible, or mono-Z, and discuss each signature in detail below.

h 125 + visible:
Experimental searches in the [(h 125 → bb) + (Φ j → qq )] and [(h 125 → bb) + (Φ j → bb)] final states have been performed by the ATLAS collaboration [62,63]. The future reach at the LHC in the (h 125 + visible) channel has been presented in Ref. [44] for the bbbb and bbγγ final states.

Z + visible:
For the (Z + visible) channel, the CMS collaboration has performed a search in the (Z + bb/τ + τ − ) final states [80]. No estimate for the future reach at the LHC exists, although the relevance for BSM searches has been pointed out in Refs. [23,43].
We extrapolate the reach at the 13 TeV LHC for 3000 fb −1 of data, σ 3000 fb −1

Z+vis
, from the limits set by the experimental search carried out by the CMS collaboration at √ s = 8 TeV with L = 19.8 fb −1 of data [80]. We rescale the reported limit σ 8 TeV; 19.8 fb −1
where σ √ s ggh (m) is the gluon fusion production cross section of a SM Higgs boson with mass m at the LHC with center-of-mass energy √ s. Note, that this is a conservative extrapolation of the reach relying purely on the increased statistics, while the ATLAS and CMS collaborations have demonstrated in the past significant improvements on background rejection as well as increased control of the systematic errors when updating searches.
We performed a dedicated collider simulation to estimate the future reach for the mono-Z channel, the details of which are presented in Appendix D. We find that the most relevant topology giving rise to the mono-Z signature proceeds from the decay [gg → Φ i → Z + (Φ j → χ 1 χ 1 )], cf. left panel of Fig. 4. The reach is mainly controlled by the mass splitting ∆m Φ = m Φ 1 − m Φ j + m Z , which controls the E miss T in the process since the (Φ j → χ 1 χ 1 ) system and the Z boson are produced back-to-back from the schannel resonance Φ i . As long as (2m χ 1 < m Φ j ), such that the (Φ j → 2χ 1 ) decay is kinematically allowed, the reach is virtually independent of the mass of χ 1 . Note that it is then straightforward to map the obtained reach to any BSM model containing the required particles by simply computing the corresponding cross section in that model. The reach for this topology from our simulation is presented in Fig. 4. We see that the reach of the mono-Z search at L = 3000 fb −1 can be quite powerful. Particularly for ∆m Φ 300 GeV, this channel will allow to probe cross sections as small as σ[gg → Φ i → Z +(Φ j → χ 1 χ 1 )] = 10 fb.

Future Prospects of the 2HDM+S
In Sec. 2 we compared the production cross sections of the different Higgs cascade decays in the 2HDM+S, in particular, in Sec. 2.5 we discussed the relations between the 2HDM+S parameters and the prevalence of different Higgs cascade channels.
In the following, we compute the reach at the future LHC for the cascade topologies in the 2HDM+S discussed above. To this end, we need to include the final state branching ratios of the daughter Higgs bosons produced in the Higgs cascades. For each point in 2HDM+S parameter space, these branching ratios may be obtained from the partial decay widths listed in Sec. 2.3. The relevant final states are Φ → {bb, τ + τ − , χ 1 χ 1 }, where Φ = {h, a} is the mass eigenstate produced as the daughter state in the Higgs cascades. The first two final states are relevant for the (Z/h 125 + visible) search modes, while the last final state gives rise to mono-Higgs/Z signatures. Note, that the additional fermion χ 1 is not necessarily part of the 2HDM+S model, however if such a fermion were to be stable on collider time scales, its presence would produce the mono-Higgs and mono-Z signatures. If such a fermion is stable on cosmological time scales and is produced in the early Universe via some mechanism, it may yield a viable Dark Matter candidate. This possibility has been discussed in detail in Ref. [14]. In general, the couplings of such a fermion to the different 2HDM+S Higgs bosons are interrelated, cf. Ref. [14], however, for the purposes of this work we treat them as free parameters.
In Fig. 5 we show representative regions of 2HDM+S parameter space within reach of the different Higgs Cascade search modes at the future LHC with 3000 fb −1 of data. We have fixed the parent Higgs boson masses to 750 GeV, high enough that there is no kinematic suppression for any of the SM channels, and standard searches (in particular tt final states) are challenging. 7 Other relevant parameters are chosen as shown in the legend. In particular, we choose tan β = 2, the mass of the additional fermion is set to m χ 1 = 50 GeV, and its couplings to the light additional Higgs bosons to g hχ 1 χ 1 = 0.01 and g aχ 1 χ 1 = 0.05. These benchmark values of the couplings are chosen such that as long as decays of the h/a into pairs of top quarks are kinematically forbidden, g hχ 1 χ 1 {g hff , g aff } g aχ 1 χ 1 , where g Φff stands for the relevant couplings of Φ to pairs of SM fermions,  [91,92], but no longer significant [93,94]) 750 GeV di-photon excess purely by accident. Since except for the charged Higgs no additional charged particles beyond the SM are present in the model considered here, it would be very challenging, if not impossible, to obtain the required di-photon cross sections to match the previously observed excess.
while cascade decays of the type (A → h 125 a) and (H → Za) will in general lead to mono-Higgs and mono-Z signals. This behavior is most pronounced when comparing the (Z + visible) with the mono-Z region in the left panel of Fig. 5: Large values of (P S A ) 2 suppress the gluon fusion production cross section of A and hence the cross section of (A → Zh) cascades giving rise to (Z + visible) final states, while large values of (S S H ) suppress the gluon fusion production cross section of H and hence (H → Za) cascades giving rise to mono-Z final states. Note, that since the decay width for In order to understand the behavior of the mono-Higgs and (h 125 + visible) regions, recall that for the chosen values of the Higgs basis trilinear couplings ( / tan β is much larger than the couplings to other SM states or the benchmark values chosen for the g Φχ 1 χ 1 couplings unless tan β is large or Φ is very singlet like, corresponding to C NSM Φ → 0. In this region, the signal would be tth 125 or ttZ. Such decays may be probed by SM tth 125 searches in the near future [95,96]. Observe that this region would be within reach of the mono-Z and mono-Higgs searches for large values of the g Φχ 1 χ 1 1 couplings. Alternatively, if there is some departure from the alignment limit such that h has some non-negligible H SM component, one may search for final states arising from (h → ZZ/W W ) decays.
Note that the cut-off for low masses of h or a for the mono-Higgs and mono-Z searches is due to the (h/a → χ 1 χ 1 ) searches being kinematically forbidden for m h /m a < 2m χ 1 = 100 GeV. For the (h 125 /Z + visible) final states the cutoff mainly comes from the fact that Refs. [44,80] provide the reach only down to a minimal value of m h /m a . Observe that the region for the (Z + visible) final state is not as smooth as for the other channels, in particular in the region (150 m h /GeV 250, 100 m a /GeV 300) because the reach for this final states is extrapolated from the experimental result in Ref. [80]. Hence, the reach as a function of the involved masses suffers from statistical fluctuations, while the reach in the remaining final states is evaluated from theoretical predictions based on Monte Carlo simulations, which suffer much less from statistical fluctuation and give rise to smooth contours in the mass plane.
The remaining behavior portrayed in the right panel of Fig. 5 Higgs Casacades and are hence controlled by the mass of h. As discussed above, this separation is much clearer when comparing the (Z + visible) and mono-Z signatures than when comparing the (h 125 + visible) and mono-Higgs channels due to the values chosen for the g Φχ 1 χ 1 couplings.
Note that while we have fixed the mass of the parent Higgs bosons at 750 GeV when presenting our results, and assumed perfect alignment, the effect of varying these quantities is easy to deduce. First, different masses for the parent Higgs bosons would primarily affect the gluon fusion production cross section, whose scaling with mass is well known. The affect of misalignment would be quantitatively negligible on the reach of the Higgs Cascades discussed here. However, various decay chains not considered in the above analysis, such as (H → h 125 h 125 ) or (A → Zh 125 ) would be present. As discussed in Sec. 2.4, such decays are suppressed by either the NSM or S component of h 125 compared to decays into h. The reach for such decays can be extrapolated from those presented by convoluting the relevant decay widths with the misalignment of h 125 and identifying m h = 125 GeV in the right panel of Fig. 5. Observe that our results can also be mapped to the case of the decoupled singlet and non-degenerate CP-odd and CP-even NSM-like Higgs bosons by appropriately choosing the NSM and S components for the parent and daughter Higgs bosons in the decay chain of interest (c.f. discussion in Sec 2.4).
Finally, note that Fig. 5 is meant to portray only the prospects of exploring the 2HDM+S parameter space using Higgs cascades. The regions displayed do not take into account existing bounds from searches for additional Higgs boson beyond h 125 . In particular, the charged Higgs does not play a role in any of the searches shown. Existing constraints on charged Higgs bosons, e.g. from flavor physics observables, can be satisfied by choosing a sufficiently large mass of the charged Higgs. Recall that in this work we treat the masses of the physical Higgs bosons as free parameters. Even without considering effects of mixing, mass splittings of order of a few 100 GeV between the mostly doublet-like pseudo-scalar and the charged Higgs are easily achievable; from the equations given in Appendix A we find for example [M 2 11 is the element of the squared mass matrix corresponding to the A NSM interaction eigenstate. Furthermore, in more complete models with larger particle content than the 2HDM+S considered here such as the NMSSM, indirect observables such as those from flavor physics receive additional contributions beyond those from the charged Higgs which may loosen the bounds on the mass of the charged Higgs, cf. Refs. [97,98].
In summary, as evident from Fig. 5, there are large regions of parameter space in reach of the different Higgs Cascade search modes. In particular, the Higgs Cascades enable the LHC to probe regions of parameter space challenging to access with traditional searches for the direct decays of additional Higgs states: Singlet-like light states are difficult to directly produce due to the small couplings to pairs of SM particles. On the other hand, doublet-like states are readily produced, but if their mass is above the kinematic threshold allowing for decays into pairs of top quarks, Φ → tt decays will dominate over the decays into other SM states. Pairs of top quarks produced from an s-channel resonance are very difficult to detect at the LHC due to interference effects with the QCD background, which makes the m Φ 350 GeV, low tan β region extremely challenging to probe at the LHC through direct Higgs decays with current search strategies [31][32][33][34][35][36][37][38].

Conclusions
The 2HDM+S is well motivated both by phenomenological considerations such as facilitating baryogenesis or DM model building, and since interesting high energy models with such a Higgs sector exist, e.g. the NMSSM. Here, we presented the first systematic study of this model. The seemingly large number of parameters complicates analyzing the model space. By mapping the parameters of the generic Higgs potential to physically relevant quantities like masses and mixing angles, and defining the remaining set of independent trilinear and quartic couplings, we obtain a much more intuitive understanding of the 2HDM+S parameter space. This mapping also makes transparent how embedding the 125 GeV Higgs boson observed at the LHC and its required SM-like couplings constrain the 2HDM+S parameter space.
The extended Higgs sector of the 2HDM+S allows for a rich phenomenology beyond what is encountered in 2HDMs. In particular, the presence of the singlet allows for interesting cascade decays of the heavy Higgs bosons into two lighter Higgs bosons or a light Higgs and a Z boson. Using such Higgs cascades in addition to the conventional searches for additional Higgs bosons allows to cover a large part of the model space despite the seemingly large number of free parameters. Although short of direct evidence (such as observing all the 2HDM+S Higgs bosons), Higgs cascades are also an important handle for differentiating such a model from e.g. the 2HDM, where Higgs cascades are hard to realize.
We compute the various production cross sections and branching ratios of the Higgs bosons (both CP-even and odd), allowing us to analytically understand which model parameters control the phenomenology. We show that cascade decays of a parent Higgs bosons into final states involving a single 125 GeV Higgs boson or a Z can have relevant cross sections. In particular, we show that while in certain regions of parameter space these two can be complimentary, for the most part, both signatures are comparable.
For collider searches, one must also consider the decays of the daughter Higgs bosons produced in the Higgs cascade decays. Besides decays into pairs of SM particles, these light Higgs bosons may also decay into new neutral states which, if stable on collider times scales, give rise to missing energy signatures. In order to include models with such additional states, which may serve as a Dark Matter candidate, we add an arbitrary fermion to the 2HDM+S. The final states arising from Higgs cascades can then be categorized as h 125 /Z+visible or h 125 /Z+invisible (mono-Higgs/Z). Convoluting the production cross sections and branching ratios and comparing to the estimated sensitivity of the LHC with L = 3000 fb −1 of data in the various final states, we show that most of the interesting region of 2HDM+S parameter space can be probed by combining the different Higgs cascade modes.
The NMSSM is arguably the best motivated high scale model with a 2HDM+S Higgs sector. The Z 3 -invariant NMSSM usually discussed in the literature represents a very constrained version of the 2HDM+S with only 7 free parameters controlling the Higgs sector. The general NMSSM is much more complicated, but can be mapped trivially to a generic 2HDM+S. Without any direct evidence at the LHC for non-SM like Higgs bosons, it is important and timely to consider all the different types of signatures that may arise in the generic situation, in particular since the low tan β region of parameter space is challenging to probe experimentally above the top threshold.

A Mass Matrices
The entries of the (symmetric) squared mass matrices are obtained from the scalar potential via (A.1) By construction, the neutral and charged Goldstone modes G 0 and G ± , respectively, are massless and do not mix with the remaining states. The mixing of the charged Higgs H ± with the neutral states is forbidden by gauge invariance. Finally, CP conservation forbids mixing of CP-even and CP-odd states. Thus, the non-trivial entries of the mass matrix can be split into 6 entries concerning the CP-even Higgs bosons, 3 entries concerning the CP-odd states, and the (squared) mass of the charged Higgs.
The entries of the CP-even Higgs squared mass matrix in the extended Higgs basis {H SM , H NSM , H S } at tree level are 8 where and we used a shorthand notation s β ≡ sin β, c β ≡ cos β, etc.
The entries of the (symmetric) squared mass matrix for the CP-odd Higgs bosons in the basis {A NSM , A S } are at tree level given by The mass of the charged Higgs boson is given by (A.14)

B Trilinear couplings
The trilinear couplings are obtained from the scalar potential via We list them for the interaction states of the extended Higgs basis in Table 1 in terms of the parameters in the scalar potential 9 . Couplings which are strictly vanishing are marked green, while couplings which are proportional to other couplings are marked blue. Besides those couplings marked in the tables, the couplings can be written in terms of the entries of the mass matrices.
C Independent quartic couplings The quartic couplings are obtained from the scalar potential via (C.1) We list the independent quartic couplings which can not be obtained as linear combinations of the entries of the mass matrix or the trilinear couplings in Table 2.  Figure 6: Illustration of resonant decay channels yielding mono-Z signals at the LHC. For the left diagram, Φ i = H and Φ j = a, or Φ i = A and Φ j = h. Decays of the type A → Zh 125 are suppressed by alignment, see the discussion in Sec. 2.3. This topology is henceforth referred to as the scalar topology. For the right diagram, Φ = H or Φ = A, and we refer to this topology as the fermion topology.
In this Appendix, we obtain the future reach in the mono-Z channel in as model independent a fashion as possible. We consider two different topologies giving rise to signals at the LHC, depicted in Fig. 6: The scalar topology [gg → Φ i → Z + (Φ j → χ 1 χ 1 )] where the intermediate states are part of the Higgs sector and the fermion topology [gg → Φ → χ 1 + (χ j → χ 1 + Z)] where the scalar Φ in the s-channel decays into a pair of new fermions χ 1 χ j and the χ j in turn decays into a Z-boson and a χ 1 fermion. One may understand the result as derived in a simplified model framework, where in the case of the scalar topology the SM is extended by two scalars Φ i , one of them CP-even and one CPodd, and a new fermion χ 1 , and in the case of the fermion topology by one scalar Φ and two new fermions χ j and χ 1 . In both cases, χ 1 must be stable on collider time scales and carry no color or electric charge such that it would not deposit energy in the detector, as would naturally be the case if χ 1 is a Dark Matter candidate. It is then straightforward to map the obtained reach to any BSM model containing the required particles by simply computing the corresponding cross section in that model. In particular, in the 2HDM+S the scalar topology arises if one new (stable) fermion χ 1 with mass m χ 1 m h /2 (m χ 1 < m a /2) is added and in addition the scalar mass spectrum is arranged such that the primary decay is allowed, i.e. m A > m Z + m h (m H > m Z + m a ). Likewise, the fermion topology arises if two additional fermions χ j and χ 1 are added to the model and the mass spectrum is arranged such that m χ j + m χ 1 < m Φ and m χ 1 + m Z < m χ j , where Φ can be any of the non SM-like Higgs bosons in the 2HDM+S.
Before going any further we can already note that, in general, for similar mass spectra, the reach in the scalar topology will be better than in the fermion topology, analogously to what discussed in Ref. [43] for the mono-Higgs signal: In the scalar topology, the E miss T is given by the transverse momentum of the intermediate scalar Φ j produced in the (Φ i → Z + Φ j ) decay back-to-back with the Z boson in the transverse plane. In the fermion topology on the other hand there is no such correlation, it is possible that the χ 1 's in the final state are orientated such that the (vectorial) sum of their respective transverse momenta approximately cancels and hence no E miss T is observed in the event. In order to project the reach at the future LHC, we simulate signal samples in both the scalar and the fermion topology using the Monte Carlo generator Madgraph5 aMC@NLO v2.5.4 [99] for the generation of the hard event, pythia8 [100,101] for showering, and the fast detector simulation Delphes3 [102]. We compare the signal sample for each set of masses with the background taken from Ref. [103] in the four different signal categories employed in [103], {(Z → µ + µ − ) + E miss T , (Z → e + e − ) + E miss T } ⊗ {HM, LM}, where HM stands for the high mass and LM for the low mass search performed in [103] by employing the same set of cuts 10 .
In the LM search, the final discriminating variable is E miss T . In the HM search, the discriminating variable is the transverse mass where p T is the transverse momentum of the lepton system. It is worth noting that the LM search employs a lower E miss T cut, E miss T > 90 GeV, than the HM search, E miss T > 120 GeV. 10 Note, that our Monte Carlo configuration differs from that used in Ref. [103]; in particular, we use the fast detector simulation Delphes3 (with the default ATLAS card in Madgraph5 v2.5.4) instead of a full detector simulation. This may have numerical impact, although beyond the scope of this work. We then obtain the future reach at the LHC, σ min , in both the HM and LM searches by adding the signals in the Z → e + e − and Z → µ + µ − categories, comparing to the background taken from Ref. [103] rescaled to L = 3000 fb −1 , and performing a cut-and-count analysis in the signal region where we add an additional lower cut in the discriminating variable of the respective search. We consider a point to be within reach if where S (B) is the number of signal (background) events after all cuts. In the first condition, the two terms in the denominator parameterize the statistical and systematical uncertainty in the background, respectively. In Ref. [103], the systematic uncertainty of the background was approximately ∆ ≈ 8 %. Assuming moderate improvement of the systematic uncertainty of the background, we assume ∆ = 5 % for our estimate of the reach.
We show the estimated reach at the LHC for L = 3000 fb −1 of data in Fig. 4, located in the main text, for the scalar topology and in Fig. 7 for the fermion topology, where for each mass hypothesis we show the better reach out of the HM and LM search. Note, that unless the mass spectrum is close to the kinematic edges indicated in Figs. 4 and 7, the HM search usually has better reach since it is optimized for larger E miss T . As discussed in the main text, in the scalar topology shown in Fig. 4, the reach is mainly controlled by the mass splitting ∆m Φ = m Φ 1 − m Φ j + m Z , which controls the E miss T in the process since the (Φ j → χ 1 χ 1 ) system and the Z boson are produced back-to-back from the s-channel resonance Φ i . As long as (2m χ 1 < m Φ j ), such that the (Φ j → 2χ 1 ) decay is kinematically allowed, the reach is virtually independent of the mass of χ 1 . For ∆m Φ 300 GeV, the reach of the mono-Z search at L = 3000 fb −1 is better than σ[gg → Φ i → Z + (Φ j → χ 1 χ 1 )] = 10 fb.
In the fermion topology, the reach depends on m Φ , which controls the overall energy scale, as well as the mass splittings at the two vertices, [m Φ − m χ j + m χ 1 ] and [m χ j − (m χ 1 + m Z )]. Since the reach is controlled by three masses, we show it in the m χ j − m χ 1 plane for different values of the mass of the parent Higgs boson Φ in Fig. 7. The reach is significantly weaker than in the scalar topology because the Z boson is produced at the secondary vertex which leads to softer E miss T spectra than in the scalar topology where the Z boson is produced back-to-back with the Φ j at the primary vertex.

E Mapping to the NMSSM
The Next-to-Minimal Supersymmetric Standard Model (NMSSM) has the same Higgs sector as the 2HDM+S, albeit its parameters are much more constrained due to supersymmetry. In this Appendix, we provide a mapping of the parameters of the 2HDM+S's scalar potential to those of the NMSSM.
The scalar potential of the general CP-conserving NMSSM can be written as (cf.
where M 4 is a new parameter of dimension mass 4 , which has no effect on the model's phenomenology, we can map the parameters of the 2HDM+S onto the general NMSSM as given in Tab. 3. Note that the mapping requires use of the identity  Table 3: Mapping of the 2HDM+S parameters to the general NMSSM. The Z 3 NMSSM is obtained by setting µ = µ = m 2 3 = m 2 S = ξ F = ξ S = 0.