Anatomy of Exotic Higgs Decays in 2HDM

Large mass splittings between new scalars in two-Higgs-doublet models (2HDM) open a key avenue to search for these new states via exotic heavy Higgs decays. We discuss in detail the different search channels for these new scalars at the LHC in the presence of a sizable mass splitting, i.e. a hierarchical 2HDM scenario, taking into account the theoretical and experimental constraints. We provide benchmark planes to exploit the complementarity among these searches, analyzing their potential to probe the hierarchical 2HDM parameter space during LHC Run 2.


Introduction
Analyses of the results from the LHC 7-8 TeV run by both ATLAS and CMS show that the properties of the Higgs particle at m h ∼ 125 GeV are close to those expected for the Standard Model (SM) Higgs boson h SM [1,2]. The complete nature of the scalar sector responsible for electroweak (EW) symmetry-breaking, however, remains to be determined, and it is particularly interesting to ascertain whether the Higgs sector consists of only one SU(2) L scalar doublet or has a richer structure containing additional states. Addressing this question is a key task for present and future studies at the Large Hadron Collider (LHC).
Two Higgs doublet models (2HDM) constitute the prime example of a well-motivated extended Higgs sector, appearing in many extensions of the SM such as the MSSM [3], composite Higgs models [4] and viable EW baryogenesis scenarios [5]. In addition to the SM-like CP-even Higgs boson, the 2HDM spectrum contains one more CP-even Higgs, a CP-odd Higgs and a pair of charged ones * . In recent years, its allowed parameter space has been scrutinized in light of ATLAS/CMS Higgs coupling measurements and searches for extra Higgses at the LHC [6][7][8][9][10][11][12][13][14][15].
A key avenue to probe the 2HDM heavy Higgs bosons at the LHC which has started to attract attention recently is the search for exotic decays of the heavy Higgses in the presence of a sizable mass splitting among them [16][17][18][19][20][21] (see also [15]). These sizable splittings are difficult to realize in the MSSM, while in more general 2HDM scenarios they may lead to important physical consequences † . While the conventional decay channels of a heavy Higgs into two SM quarks, leptons or gauge bosons have been the focus of most of the existing searches, the exotic (non-SM) modes of a heavy Higgs decaying into two light Higgses, or one light Higgs with one SM gauge boson quickly dominate once they are kinematically open. The current exclusion bounds on extra Higgses based on their conventional decays only will be therefore significantly relaxed. On the other hand, the exotic decay modes offer new discovery channels, which have already shown exclusion power during the 8 TeV LHC run [22,23], and yield very promising prospects for the 13 TeV LHC run. In this work, we aim to provide a comprehensive categorization and analysis of the exotic search channels for the new 2HDM scalars, highlighting the complementarity among them, and provide guiding benchmark planes for Run 2 of the LHC at 13 TeV.
After a review of the 2HDM in Section 2, we present the constraints on the 2HDM parameter space coming from theoretical considerations (stability of the EW minimum, perturbativity and tree-level unitarity) and experimental measurements in Section 3, where we also introduce the salient features of our benchmark scenarios for exotic 2HDM Higgs decays (Section 3.6) motivated by the theoretical and experimental constraints. In Section 4 we discuss the production and decay of non-SM Higgses at the LHC, and then analyze in depth our different benchmark scenarios in Section 5, before concluding in Section 6.

2HDM Lagrangian and Higgs Potential
In the 2HDM, we introduce two SU(2) L doublets Φ i (i = 1, 2): where v i are the vacuum expectation values (vev) of the neutral components, satisfying v 2 1 + v 2 2 = v 2 , with v = 246 GeV. The ratio of vevs is defined as tan β ≡ v 2 /v 1 . The 2HDM Lagrangian for Φ i can be written as where the first term denotes the kinetic term for the two Higgs doublets, V (Φ 1 , Φ 2 ) is the Higgs potential and the last term denotes the Yukawa interactions between Φ i and the SM fermions. Assuming CP conservation and a soft Z 2 symmetry breaking, the 2HDM Higgs potential can be written down as ‡ :

(2.3)
After EW symmetry breaking, the physical 2HDM scalar spectrum consists of five states: two CP-even Higgses h, H with m h < m H , a CP-odd scalar A and a charged scalar pair H ± [24], which may be written as with the angle α parametrizing the mixing between the neutral CP-even components (we use the shorthand notation s x ≡ sin x, c x ≡ cos x, t x ≡ tan x). The Goldstone bosons G and G ± are absorbed as longitudinal components of the Z and W ± bosons. In the limit c β−α = 0 (the alignment limit for h), the state h can be identified with the SM Higgs, its couplings to fermions and gauge bosons being precisely those predicted by the SM § . It is thus convenient to describe the model in terms of t β , c β−α , the physical scalar masses m h , m H , m A , m H ± , the soft Z 2 symmetry breaking parameter m 2 12 and the vev v. The quartic couplings in Eq. (2.3) can be expressed in terms of the physical masses and mixing angles as (see e.g. [26]) s β c β . (2.5)

Interactions in the 2HDM
The couplings of the CP-even scalars to a pair of gauge bosons, arising from the Higgs kinetic term in Eq. (2.2), are [24] g hZZ = 2im 2 (2.6) The CP-odd scalar A does not couple to pairs of vector bosons, while the charged scalar H ± only couples to pair of vector bosons at loop level. In addition, the couplings of two scalars and one vector boson read in which p µ are the outgoing momentum for the corresponding particle. The hHZ-coupling is absent due to CP conservation. We note that, considering h (H) to be the SM-like 125 GeV Higgs with c β−α = 0 (s β−α = 0), gauge boson couplings to two non-SM like Higgses are unsuppressed, while the gauge boson couplings to h (H) and one non-SM like Higgs are suppressed by c β−α (s β−α ).
Regarding the cubic couplings among scalars arising from the 2HDM scalar potential Eq. (2.3), the relevant ones for our analysis are  [24] assume the MSSM relation m 2 12 = m 2 A s β c β ). We also stress that for a light CP-odd scalar A with m A < m h /2, the decay channel h → AA could be open, being however very constrained experimentally ¶ (see [28] for a discussion of this region of the 2HDM parameter space).

State
Up-type fermions Finally, as is well-known the couplings of the 2HDM scalars to SM fermions, contained in L Yuk in Eq. (2.2) are not univocally determined by the gauge structure of the model. In the presence of a Z 2 symmetry guaranteeing the absence of tree-level flavour changing neutral currents [29], four possible 2HDM types exist (see [30] for a discussion). The couplings of the neutral scalar states to SM fermions, normalized to their SM values, can be expressed in terms of functions of α and β, shown in Table 1 for the particular case of a Type II 2HDM (one Higgs doublet Φ 2 couples to the up-type quarks, while the other Higgs doublet Φ 1 couples to the down-type quarks and leptons).

The Alignment Limit and the Role of m 2 12
It is useful to cast the relations between the quartic couplings and the physical masses Eq. (2.5) in terms of c β−α , which characterizes the departure from the alignment limit for h Current data from LHC Run 1 favour the alignment limit c β−α = 0 [31] (see also [6-8, 10-12, 14] The combination m 2 12 −m 2 H s β c β in Eq. (2.10) will play a key role in the following discussion: the value of m 2 12 is not fixed by the mass spectrum or the scalar couplings to gauge bosons and fermions, only entering the trilinear scalar couplings Eq. (2.8). Its possible allowed values are dictated by theoretical constraints on the 2HDM parameter space, namely the boundedness from below of the scalar potential Eq. (2.3) and the stability of the EW minimum, and the requirements of perturbativity and tree-level unitarity on the quartic couplings λ i , as shown in the next section. These have a large impact on the allowed values of masses m H , m A , m H ± , m 2 12 and t β (and c β−α away from alignment), as the absence of a value of m 2 12 satisfying the theoretical constraints for a given set of values for m H , m A , m H ± and t β , indicates that such set of values is not physically viable (see e.g. [15]).

Vacuum Stability
In order to have a stable vacuum, the following conditions need to be fulfilled [26] For c β−α = 0, satisfying the first two conditions requires m 2 12 − m 2 H s β c β 0 for either t β > 1 or t β < 1, as seen from Eq. (2.10). Moreover, Eq. (2.9) shows that a departure from alignment generically has a negative impact on the first two stability conditions. Focusing on the alignment limit, the first two requirements are automatically satisfied for m 2 12 − m 2 H s β c β = 0, with the last two given by This implies that for m H > m A , m H ± , the mass splittings between the heavy CP-even Higgs H and the other heavy scalars A and H ± have to be small, such that the decays of H into AZ, AA, H + H − or H ± W ∓ are not kinematically allowed. For m 2 12 = 0 all four stability conditions of Eq. (3.1) are automatically fulfilled. The allowed region in the m 12 vs. t β plane is shown in the left panel of Figure 1 for m A = m H ± = 400 GeV and m H = 200, 300, 400 GeV as an illustration. As seen from Figure 1, the regions m 2 12 < m 2 H s β c β are generically allowed by the vacuum stability requirement.

Perturbativity and Unitarity
Upon imposing the perturbativity condition |λ i | ≤ 4π, the strongest constraints in the alignment limit come respectively from Thus, perturbativity requires m 2 12 − m 2 H s β c β v 2 unless t β ∼ 1. Moreover, even for m 2 12 = m 2 H s β c β , perturbativity of λ 3−5 imposes constraints on the size of the mass splittings among the new scalars.
Even stronger constraints are found when requiring tree-level unitarity of the scattering matrix in the 2HDM scalar sector [33]. The eigenvalues of the scattering matrix read Performing a partial-wave expansion of the scattering amplitudes yields limits on the partial wave amplitudes, which for the J = 0 case translate into the constraint |Re(Λ i )| < 8π (see e.g. [34]), which we consider here. A quick inspection of Eq. (3.3) shows that for t β 1 the scattering matrix eigenvalues scale as Λ 7,9,11 ∼ λ 1 (particularly Λ 11 3 λ 1 ), which An analysis of unitarity constraints at one-loop level has been performed in [34].
again imposes m 2 12 − m 2 H s β c β v 2 (and yields an even stronger constraint than the perturbativity one). A similar argument follows for t β 1, this time with Λ 7,9,11 ∼ λ 2 . As a result, m 2 12 ≈ m 2 H s β c β is strongly preferred unless t β ∼ 1, as shown explicitly in the right panel of Figure 1 (for m A = m H ± ). In the limit m 2 12 = m 2 H s β c β , the scattering matrix eigenvalues from Eq. (3.3) become independent of t β (in alignment c β−α = 0) and read such that |Λ i | < 8π (note that Λ i are real) impose upper limits on the mass splittings (although not on the masses themselves). We also note that for m 2 12 = 0, Λ 1−6 are independent of t β (depending only on the scalar masses) while Λ 7−12 do depend on t β , which once again results in t β ≈ 1 being the only accessible region for large mass splittings in this case.

Electroweak Precision Measurements
Measurements of EW precision observables (EWPO) impose strong constraints on the 2HDM mass spectrum. Adopting the current 95% C.L. constraints on the S and T oblique parameters (with U = 0) [35], the allowed region of parameter space in the (m A , m H ± ) plane is shown, for c β−α = 0 (neither t β nor m 2 12 affect S and T ), in the left panel of Figure 2 respectively for m H = 400 GeV (red), m H = 300 GeV (blue) and m H = 200 GeV (green). Satisfying EWPO constraints requires the charged scalar mass to be close to one of the heavy neutral scalar masses: m H ± ≈ m H or m H ± ≈ m A .  Away from the alignment limit, additional contributions to S and T proportional to c β−α appear [24] (see also [36]), such that the scenario m H = m A = m H ± is only allowed for small |c β−α | once m H v is realized, as shown in the right panel of Figure 2. The departure from alignment also allows for mild mass splittings among all the new scalars (e.g. m A > m H + m Z and m H m H ± + m W ), which however does not significantly alter the phenomenology of exotic Higgs decays at the LHC, discussed in detail in Section 5.

Flavour Constraints
Various flavour measurements [37] provide indirect constraints on the charged scalar mass m H ± as a function of t β . The different limits are computed for the case of a Type II 2HDM with SuperIso [38,39], and shown in the left panel of Figure 3. The most stringent of these * * comes from the measurement of the branching fraction (BR) of b → sγ (B 0 d → X s γ), which sets a limit m H ± > 480 GeV at 95 % C.L. [43] (we note that the limit is even stronger for t β < 2). For large t β 20, the lower limit on m H ± set by the measurement of the branching fraction B + d → τ + ν is significantly stronger, with m H ± 700 GeV for t β = 30. Similarly, the region t β 1 is very strongly constrained by B 0 s → µµ and ∆m B d . For m H = 125 GeV and s β−α = 0, when the heavy CP-even scalar H is the SM-like Higgs, the mass of the light CP-even Higgs h is constrained by flavour measurements as well, as shown in the right panel of Figure 3 for Type II 2HDM. The strongest constraint in this case comes from B 0 s → µ + µ − , which can exclude up to m h < 100 GeV (the precise bound depending on m 2 12 and t β ) for masses m H ± satisfying the b → sγ constraint. Note that flavour constraints are typically very model dependent. Contributions from additional sectors in the model could relax the constraints, as has e.g. been studied in the MSSM framework for b → sγ [44]. Being mostly focused on the collider aspects of 2HDM Higgses, we will not consider flavour as a hard constraint in the following, however indicating its effect on the parameter space under consideration.

LHC and LEP Constraints
We now review the constraints from direct searches of the new scalars. Besides the LEP bound m H ± > 80 GeV (72 GeV) for Type II (I) 2HDM [27], LEP searches for e + e − → * * We note here that the recent measurement from the BaBar Collaboration of the ratios of B → D * τ ν to B → D * ν decays and B → Dτ ν to B → D ν decays cannot be accommodated within the Type II 2HDM [40]. However, a new measurement of the former ratio by the Belle Collaboration [41,42] is in tension with this conclusion. Since this matter is not settled yet, we choose not to include these flavour measurements in our discussion.
A H (H → bb/τ τ, A → bb/τ τ ) constrain the sum of the masses m A + m H 209 GeV [45]. At the LHC, the searches for A/H in bb-associated production and decaying to τ τ by ATLAS/CMS [46,47] constrain the high t β region in the Type II 2HDM. Away from alignment, searches by ATLAS/CMS for [51,52] and H → hh → bbγγ, bbbb [53][54][55] yield strong constraints on the (c β−α , t β ) plane as a function of the respective mass m H /m A (see e.g. [13,15]). We however stress that the limits summarized above can be significantly weakened once exotic Higgs decay channels are open [15-17, 19, 21]). Searches for these new channels, e.g. via A/H → HZ/AZ are then crucial for probing 2HDM scenarios with large mass splittings among the new states (i.e. hierarchical 2HDM scenarios), and there is already ongoing effort by CMS in this direction [22,23].
Finally, ATLAS/CMS searches impose constraints on the charged scalar [56,57] beyond those of LEP. A light charged scalar m H ± m t is mostly excluded by the non-observation of the decay t → H + b → τ νb where the top is produced in top pair production. For m H ± > m t , the current limit is very weak and only constrains the high t β region for m H ± not much above the top mass (see [19] for a detailed discussion).

From Constraints to 2HDM Benchmarks
The combination of previous constraints provides a key guideline to the design of simplified 2HDM benchmark scenarios for LHC Run 2 searches at 13 TeV. EWPO measurements require the mass of the charged scalar to be close to the mass of one of the neutral scalars, and so we fix m H ± = m H or m H ± = m A in the following. In addition, measurements of Higgs signal strengths at the LHC favour the alignment limit (c β−α = 0 if h is the 125 GeV SM-like Higgs), particularly for Type II 2HDM. We then focus our analysis mostly on the alignment limit, and only consider deviations from alignment when discussing possible decays of the new scalars into the SM-like Higgs h.
Regarding the impact of theoretical constraints on the 2HDM parameter space, the previous discussion shows that satisfying unitarity/perturbativity and vacuum stability bounds (close to the alignment limit) for arbitrary values of t β requires m 2 12 = m 2 H s β c β and m H m A , m H ± . Alternatively, stability is satisfied for any mass ordering if m 2 12 = 0, while unitarity requires in this case a low value of t β . We thus consider these two scenarios as benchmark cases for our analysis: , vacuum stability requires m H m A and m H m H + , and thus the exotic decays H → A Z and H → H ± W ∓ are not kinematically allowed.
-Unitarity requires |Λ i | < 8π, constraining the mass differences among the new scalar states (but not the absolute mass values). In particular, using Eq. (3.4) we obtain the bound 5(m 2 The authors of [58] have shown that this case is preferred when requiring the 2HDM to be stable up to the Planck scale. In the alignment limit c β−α = 0 all these couplings vanish, and therefore the decays H -Unitarity imposes an upper bound on the scalar masses (not only on the mass splittings). This bound scales as t −2 β for t β > 1 and as t 2 β for t β < 1, such that only the region t β ∼ 1 is allowed (we recall that in Type II 2HDM, at least one of the neutral scalars needs to be heavy due to the combination of EWPO and flavour constraints).
-The cubic scalar couplings Eq. (2.8) now read In the alignment limit the coupling g Hhh vanishes and thus the decay H → hh is absent. However, the couplings g HAA and g HH + H − are non-vanishing as long as t β = 1.
For our analysis of benchmark scenarios away from alignment, which focus on the decays of A, H, H ± into the SM-like Higgs h, we consider the same two cases above for consistency (even though these cases are motivated by theoretical constraints for c β−α = 0).

LHC Production and Decay of 2HDM Higgses
We now discuss the salient features of the production and decay of the new 2HDM scalars at the LHC. The production of the CP-even and CP-odd neutral scalars H, A at the 13 TeV LHC occurs via gluon fusion (gg → H/A) and bb-associated production. Gluon ‡ ‡ This case has also been analyzed in [59].
fusion is the dominant production mechanism for small and moderate values of t β , while for Type II 2HDM, bb-associated production dominates at large t β . In both cases, we compute the production cross section for H and A at NNLO in QCD via SusHi [60] (for H, the cross section does depend on c β−α , and in that case we consider the alignment limit c β−α = 0). For the charged scalar H ± , the dominant production mode for m H ± > m t is in association with a tb pair, and we use the NLO cross section values provided by the Higgs Cross Section Working Group (HXSWG) for m H ± > 200 GeV [61]. A light charged scalar (m H ± < m t ) is mainly produced through top quark decays t → H + b, and we use Top++2.0 [62] to compute the top pair production to NNLO in QCD, assuming a topquark mass m t = 172.5 GeV. The LHC production cross sections for H, A and H + at 13 TeV are shown in Appendix A.
Regarding the decays of the new scalars, in the alignment limit c β−α = 0 the conventional (SM-like) decays of A and H are into tt (if kinematically accessible), bb, cc, τ τ , and with a highly suppressed branching fraction into gg, γγ and µµ. When open, the decay into tt is dominant for low and moderate t β , followed by the decay into bb. At high t β , for Type II 2HDM, the decay into τ τ becomes important, where the decay into bb can dominate even above the tt threshold. For the CP-even Higgs H, the decay into massive gauge bosons W + W − and ZZ is present away from the alignment limit, and dominates as soon as the departure from alignment is sizeable. For the charged scalar, the decay H ± → tb dominates once it is kinematically open, followed by H ± → τ ν, cs and cb. In the following, we compute all 2HDM branching fractions using 2HDMC [63]. We here stress that the above conventional decays of the new 2HDM scalars become suppressed once exotic (non SM-like) decay modes open up. These can be decays involving several states among H, A, H ± , in the presence of a large mass splitting among the new scalars (see e.g. [16][17][18][19][20][21] for existing studies on individual channels), and/or decays into the SM-like Higgs boson h, namely H → hh, A → hZ, H ± → hW ± , which are possible for c β−α = 0 and are also considered in the following as exotic (despite involving SM decay products) as they don't occur in the SM. In the former case, we can further distinguish between the decay of a new scalar into another one and a gauge boson, and the potential decays of H into either AA or H + H − . The different types of exotic decay modes for the 2HDM are summarized in Table 2.

Parent Scalar Decay Possible Final States Channels in 2HDM
The impact of the presence of exotic Higgs decay modes on the branching ratios is shown in Figure 4 for c β−α = 0. The top two panels show the relevant branching fractions The decay H ± → HW ± (solid blue) dominates with BR(H ± → HW ± ) 50%, particularly for a not too heavy state H ± . In that case, H ± → tb is suppressed to be about few percent for intermediate t β , and only reaches about 50% in the very small and very large t β region. Finally, the bottom right panel in Figure 4 shows the branching fractions of H for m H > m H ± = m A (with m H = 500 GeV and m A = 200 GeV) and m 2 12 = 0. In this case the decays H → AA and H → H + H − are allowed and dominate over most of the t β region, except for t β ∼ 1, where H → AZ and H → H ± W ∓ become dominant due to the accidental suppression of the HH + H − and HAA couplings at t β ∼ 1. Note however that for m 2 12 = 0 the theoretical constraints do not allow a significant departure from t β ∼ 1, such that a large branching fraction for H → AZ and H → H ± W ∓ is expected. Decays to SM fermions are highly suppressed in this scenario.

2HDM Planes for Exotic Higgs Decays
Our analysis of exotic Higgs decays in the 2HDM focuses on a few key benchmark planes which show the complementarity among different LHC search channels for the new scalars. We first focus on the alignment limit: c β−α = 0 for m H > m h = 125 GeV and s β−α = 0 for m h < m H = 125 GeV. In this context, we consider two possible mass planes: m A vs. m H = m H ± (Plane I) and m H vs. m A = m H ± (Plane II). These two choices are motivated by EWPO constraints (recall the discussion in Section 3.3). This is in contrast to a potential m H ± vs. m H = m A plane, highly constrained by EWPO to a small mass splitting m H ± − m H/A which closes the phase space needed for on-shell exotic Higgs decays § § , so that we don't consider such benchmark plane in our current study. Finally let us remark that, while we do not impose the flavour bounds as hard constraints on our 2HDM benchmark planes (recall the discussion in Section 3.4), we do show them as indicative in the following. Our 2HDM benchmark plane (BP) scenarios in alignment are then: As discussed in Section 3.6, this mass ordering is allowed for Case 1 (m 2 12 = m 2 H s β c β and all t β values) and Case 2 (m 2 12 = 0 and t β ∼ 1). We thus consider four scenarios: Case 1 with t β = 1.5, 7, 30 and Case 2 with t β = 1.5.
This mass ordering is not compatible with Case 1 due to vacuum stability (see Sections 3.1 and 3.6). Thus, we only consider Case 2 with t β = 1.5 ¶ ¶ .
As for BP IB, this mass ordering is not compatible with Case 1, and so we only consider Case 2 with t β = 1.5. § § As discussed in Section 3.3, a sizable departure from alignment could allow for a mass hierarchy mA > mH + mZ (such that A → HZ is kinematically allowed) and mH m H ± + mW (such that H → H ± W ∓ is kinematically allowed, but nevertheless phase space suppressed). The phenomenology of this kind of scenario is however largely contained in Planes I-II, and so we do not consider it separately.
¶ ¶ t β = 1 is chosen for the exotic decays into two lighter new scalars (H → AA in this case) not to vanish.
• BP IIB: m H < m A = m H ± .
In order to study the decays of the new scalars into the SM-like Higgs, we also consider a plane in which the departure from alignment is explored, assuming h is the 125 GeV SMlike Higgs (Plane III). We set m H = m A = m H ± for simplicity, and define the plane as m H = m A = m H ± vs. c β−α : We consider Case 1 with t β = 1.5, 7, 30 and Case 2 with t β = 1.5.
A summary of the different benchmark planes considered and the relevant exotic decay modes is shown in Table 3. In all cases, we present the σ× BR of each characteristic decay channel at the 13 TeV LHC, together with a detailed analysis of the regions disfavoured by theoretical and experimental constraints (including flavour constraints, shown for reference only). The results for Planes I and II (c β−α = 0) are presented in Section 5.1, while the results for decays to SM-like Higgs away from alignment, corresponding to Plane III, are presented in Section 5.2. Further details on the cross sections and decay branching fractions for the non-SM like Higgses can be found in Appendix A.  Before we move on to discuss in detail our different 2HDM planes for LHC searches at 13 TeV, let us comment on the comparison of these benchmark scenarios with others proposed in the literature. In particular, our Planes I and II have a substantial overlap with the 2HDM "short cascade" scenario D from [64], while our specific BP IA and BP IIB have similarities with the A → HZ benchmarks for c β−α = 0 in [15] (see also [18]). As compared to [64], the present analysis explores the full mass plane, not restricted to specific benchmark lines with fixed relations * * * among m H , m A and m H ± . We also explore the dependence on t β , which has a significant impact on the allowed 2HDM parameter space for Planes I and II. Moreover, our analysis includes the 8 TeV experimental constraint from the CMS H → AZ/A → HZ search [22,23], precisely tailored to probe these 2HDM scenarios and thus a key ingredient in a study of 2HDM exotic Higgs decays.  It is worth discussing here the extent to which a departure from our Benchmark Plane assumptions (Case 1 with m 2 12 = m 2 H s β c β or Case 2 with m 2 12 = 0, with certain mass degeneracy relations between m A , m H and m H ± ) may lead to modifications of the 2HDM phenomenology w.r.t. the scenarios we consider in this work:

Mass
• Negative m 2 12 : Extending m 2 12 to be negative does not introduce further constraints regarding vacuum stability. In addition, the variation of unitarity constraints for m 2 12 < 0 w.r.t. the m 2 12 = 0 case is mild, such that m 2 12 = 0 effectively captures this scenario. * * * In particular, we note that the fixed relations in [64] result in the exotic Higgs decays being largely subdominant above the tt threshold, which may not be the case in general (see e.g. Figure 4).
• Deviation from m 2 12 = m 2 H s β c β : Assuming the alignment limit c β−α = 0, at high t β the soft Z 2 breaking term m 2 12 is fixed to be around m 2 H s β c β by unitarity. Only when departing significantly from the alignment limit are deviations from m 2 12 = m 2 H s β c β at high t β possible, yielding a larger allowed region in parameter space.
At low values of t β deviations from m 2 12 = m 2 H s β c β can be easily accommodated in the alignment limit. We however stress that changing the value of m 2 12 only affects cubic and quartic Higgs couplings. In particular, HH + H − and HAA coupling vanish for m 2 12 = m 2 H s β c β in the alignment limit (see Section 3.6). Deviations from m 2 12 = m 2 H s β c β will then increase the branching fractions of H → AA and H → H + H − (if allowed by kinematics) with the branching fractions to other final states decreasing accordingly. For m 2 12 far away from m 2 H s β c β , the branching fraction dependence on m 2 12 is mild, with Case 2 (m 2 12 = 0) being a representative scenario.
• Deviations from the alignment limit c β−α = 0: In our discussion we mainly focus on the alignment limit, preferred by Higgs coupling measurements at the LHC. In particular, we find that in the alignment limit, unitarity and vacuum stability cannot be satisfied simultaneously if m H > m A or m H ± at high t β and therefore decays of the type H → AZ, H → H ± W ∓ are not permitted for large t β . This statement can be relaxed when deviating from the alignment limit. For moderate values of |c β−α | around 0.2−0.6, there are regions of parameter space with m H > m A , m H ± , which are allowed by both unitarity and vacuum stability. Note however that once away from the alignment limit, channels like H → AZ, H ± W ∓ will generically have reduced branching fractions. There have been studies in the literature for exotic Higgs decays in those non-alignment region [18,65,66].
• Deviation from mass degeneracy: The mass degeneracy relation is mainly imposed by the electroweak precision measurements. As shown in the left panel of Fig. 2, small deviations from m H ± ≈ m A/H are allowed. In principle, there are small regions of parameter space corresponding to a mass hierarchy m H = m A < m H ± . However, these region of parameter space are already in almost 2σ tension with observation. Other deviations from mass degeneracy will not lead to a phenomenology that would differ significantly from that of BP I and BP II. and t β described in Table 3. Note that for t β = 1.5, 7 we consider the dominant ggA production, while for t β = 30 the bbA production dominates and is considered instead. For each panel, contour lines of 10, 10 2 , 10 3 and 10 4 fb are drawn as light grey dashed curves to guide the eye. Large cross sections σ × BR 1 pb are possible for t β ∼ 1 and t β 1, respectively due to the enhanced top and bottom Yukawa coupling contribution, even for large CP-odd scalar masses m A ∼ 500 − 600 GeV. The shaded areas enclosed by an irregular curve in Figures 5 and 6 are excluded by the CMS A → HZ search [22,23], which already constrains a sizable portion of parameter space and highlights the potential of such a search at LHC 13 TeV in the bb and τ τ final states, as a probe of both A and H.

Exotic Decays in the Alignment
Hatched regions show the parameter space excluded by other experimental searches, as well as unitarity constraints. The former exclusions are mainly due to t → H + b searches, which yield a limit m H ± > m t , as well as H → τ τ searches for large t β , which rule out m H < 600 GeV for t β = 30. We also show the flavour bound m H ± > 480 GeV as a horizontal grey line for indicative purposes.
Regarding unitarity, for Case 1 (m 2 The latter imposes the strongest constraint, which rules out regions with a very large mass splitting m A − m H (as indicated by the hatched region in the lowerright corner of each panel in Figures 5 and 6). For Case 2 (m 2 12 = 0) with m A > m H = m H + , the strongest unitarity constraints come from |Λ 6 |v 2 = 5m 2 In particular, |Λ 11 | rules out the large m H region (upper hatched region in the lower right panel in Figures 5 and 6).
Taking  of parameter space allowed, with signal cross sections σ× BR 100 fb. For Case 2, the lower right panel of Figure 7 shows that the CMS A → HZ search constrains most of the viable parameter space, which may in turn be probed completely by LHC 13.
While the generic features for H ± → HW ± are similar to those of A → HZ, the signal cross sections are about two order of magnitude smaller, due to the suppressed production cross section of pp → H ± tb. This, in addition to the complicated final state HW + W − bb which results, makes this channel challenging for LHC studies at 13 TeV.

BP IB: m A < m H ± = m H
In this scenario, only Case 2 (m 2 12 = 0) is viable. The σ× BR for the three possible exotic decay channels H → AZ, H ± → AW ± and H → AA is shown in Figure 9 for our benchmark t β = 1.5. The m H > 460 GeV region is excluded by unitarity, the strongest unitarity constraint coming from |Λ 11 |v 2 = 1 2 m 2 H ( 1 Below the unitarity limit on m H , the ATLAS/CMS limits on A → τ τ at low t β (ggA production) [46,47] combined with the bounds from the CMS H → AZ search [22,23] rule out m A > 40 GeV down to m H 350 GeV. As can be seen from Figure 9, only a small region of parameter space survives the unitarity and LHC 8 TeV constraints. We  Figure 10. Comparing to BP IIB, the additional collider search limit m H ± > m t applies, which overlaps with the 8 TeV LHC exclusion from A → τ τ . This results in only a small stripe in parameter space, corresponding to 200 GeV < m A = m H ± < 240 GeV and 300 GeV < m H < 450 GeV, being viable. Moreover, we note that the decays H → AA and H → H + H − are essentially not kinematically allowed in the viable region, as shown in the lower panels of Figure 10. This benchmark scenario should indeed be possible to probe completely at LHC 13 TeV. c β−α . In Figures 11, 12, and 13 we respectively show the σ × BR for A → hZ, H ± → hW ± and H → hh, in each case for Case 1 with t β = 1.5, 7, 30 and Case 2 with t β = 1.5 in the (c β−α vs m A = m H = m H ± ) plane.

Exotic Decays into
For Case 1 with t β = 1.5, only the region |c β−α | 0.2 (close to the alignment limit) is viable as a result of Higgs signal strength measurements (mainly driven by the g hV V couplings) considering all the theoretical and experimental constraints. The allowed range for c β−α shrinks as the masses of the heavy 2HDM scalars grow due to stability constraints, being already restricted to −0.02 < c β−α < 0.06 for m A = m H = m H ± 500 GeV. At the same time, LHC bounds on H → ZZ and A → τ τ rule out m A = m H = m H ± < 350 GeV. For significantly higher values of t β (as our t β = 7, 30 scenarios) the stability constraints rule out almost completely the region c β−α < 0, while unitarity imposes a strong constraint on c β−α > 0 for high scalar masses m A = m H = m H ± > 600 GeV. In addition, for t β = 7 the vacuum stability constraint rules out the region c β−α > 0.3 while Higgs signal strengths rule out the region 0.05 < c β−α < 0.24. For t β = 30, Higgs signal strengths rule out c β−α 0.01, while A → τ τ searches restrict the allowed parameter space to m A = m H = m H ± > 650 GeV, leaving only a very narrow stripe as viable parameter space. For Case 2 with t β = 1.5, satisfying the constraints from H → ZZ and A → τ τ requires m A = m H = m H ± > 350 GeV and |c β−α | 0.2, while unitarity imposes an upper bound on the scalar masses in the range 450 GeV -550 GeV depending on c β−α . As shown in Figure 11, the cross sections for A → hZ in the allowed region of parameter space could reach 1 pb or higher for t β = 1.5 both in Case 1 and 2. For t β = 7 (Case 1) the cross section for A → hZ is still sizable in the allowed region 0.24 < c β−α < 0.3, reaching values ∼ 100 fb. For t β = 30 the signal cross section is however very small due to the suppressed branching ratio BR(A → hZ) close to the alignment limit. The signal cross sections for H ± → hW ± shown in Figure 12 follow a trend similar to those for A → hZ, but being typically a factor 10 -100 smaller due to the suppressed production cross section for H ± above m t (see Appendix A.1). Finally for H → hh the signal cross sections, shown in Figure 13, are about factor of 10 smaller than those of A → hZ, and an additional suppression of the branching ratio BR(H → hh) occurs for certain values of c β−α (e.g. c β−α ∼ 0.22 for t β = 7 and c β−α ∼ 0.052 for t β = 30, as seen from Figure 13).

Conclusions
In the 2HDM, other than decaying to pairs of SM quarks, leptons, and gauge bosons, the exotic decays of heavy Higgses into two lighter Higgses or one light Higgs and a SM gauge boson are likely to dominate once they are kinematically open. While the collider search bounds for heavy Higgses based on conventional search modes W W , ZZ, γγ, bb and τ τ for neutral Higgses, and τ ν and cs modes for charged Higgses would be relaxed once those exotic modes are open, the exotic decay modes offer new discovery channels in large regions of the 2HDM parameter space. Away from the 2HDM alignment limit, exotic decays into the SM-like 125 GeV Higgs boson h, namely H → hh, A → hZ and H ± → hW ± , are potentially important, and there is already an ongoing ATLAS and CMS search programme for A → hZ [51,52] and H → hh [53][54][55]. In contrast, close to the alignment limit, as favoured by measurements of Higgs signal strengths, exotic decays among the new 2HDM scalars become particularly relevant. The experimental searches based on those channels, however, have just started with H/A → AZ/HZ [22,23]. In this work, we carefully examine the exotic Higgs decay channels in the 2HDM, both in the presence of a hierarchy between Higgses and away from alignment when this hierarchy is not present. By taking into account the various theoretical and experimental constraints, we propose 2HDM benchmark plane scenarios for LHC searches at 13 TeV: • BP IA: m A > m H = m H ± , with A → HZ, H ± W ∓ .
• BP IB: m A < m H = m H ± , with H → AZ, AA and H ± → AW ± .
• BP IIB: m H < m A = m H ± , with A → HZ and H ± → HW ± .
In each case, we analyze the allowed regions of parameter space and the LHC 13 TeV σ × BR for the relevant exotic Higgs decay modes in those regions.  To summarize, exotic Higgs decays provide new discovery avenues for heavy Higgses. In turn, the exploration of the proposed benchmarks via these decays could help to understand the structure of the electroweak symmetry breaking sector beyond the SM.

A.1 2HDM Production Cross Sections
In Figure 14, we show the gluon fusion production cross section for H (upper left panel) and A (upper right panel), bb-associated production cross section for H (middle left panel) and A (middle right panel), and tbH ± production cross section (bottom) for the charged scalar (details are given in Section 4).

A.2 2HDM Branching Ratios for Exotic Higgs Decays
For illustration, we show in Figure 15 the branching ratios of H a → H b V (with H a,b = H, A, H ± and V = W ± , Z) for Plane I (top) and Plane II (bottom) for Case 2 (m 2 12 = 0) with t β = 1.5 (being the scenario allowed for the four benchmarks BP IA, BP IB, BP IIA and BP IIB). The decay branching ratios for H/A → AZ/H(h)Z are shown on the left panels of Figure 15, while those for A/H ± → H ± W ∓ /AW ± (Plane I) and H/H ± → H ± W ∓ /HW ± (Plane II) are shown on the right panels.      Figure 14. Production cross section for H, A and H + at LHC 13 TeV. The contour lines indicate the cross section of 1, 10, 10 2 , 10 3 , 10 4 , 10 5 and 10 6 fb. [GeV]