Exotic Decays Of A Heavy Neutral Higgs Through HZ/AZ Channel

Models of electroweak symmetry breaking with extended Higgs sectors are theoretically well motivated. In this study, we focus on the Two Higgs Doublet Model with a low energy spectrum containing scalars $H$ and a pseudoscalar $A$. We study the decays $A\rightarrow HZ$ or $H\rightarrow AZ$, which could reach sizable branching fractions in certain parameter regions. With detailed collider analysis, we obtain model independent exclusion bounds as well as discovery reach at the 14 TeV LHC for the process: $gg\rightarrow A/H\rightarrow HZ/AZ$, looking at final states bbll, \tau\tau ll and ZZZ(4l2j) for l =e,\mu. We further interpret these bounds in the context of the Type II Two Higgs Doublet Model, considering three different classes of processes: $A\rightarrow h^0Z$, $A\rightarrow H^0Z$, and $H^0 \rightarrow AZ$, in which $h^0$ and $H^0$ are the light and heavy CP-even Higgses respectively. For 100 fb$^{-1}$ integrated luminosity at the 14 TeV LHC, we find that for parent particle mass around 300-400 GeV, $A\rightarrow h^0Z$ has the greatest reach when $H^0$ is interpreted as the 126 GeV Higgs: most regions in the tan\beta-sin(\beta-\alpha) parameter space can be covered by exclusion and discovery. For 126 GeV $h^0$, only relatively small tan\beta<10 (5) can be reached by exclusion (discovery). For $A\rightarrow H^0Z$, the reach is typically restricted to sin$(\beta-\alpha)\sim\pm 1$ with tan\beta<10 in bbll and \tau\tau ll channels. The ZZZ (4l2j) channel, on the other hand, covers a wide range of 0.3<|sin(\beta-\alpha)|<1 for tan\beta<4. $H^0\rightarrow AZ$ typically favors negative values of sin(\beta-\alpha), with exclusion/discovery reach possibly extending to all values of tan\beta. A study of exotic decays of extra Higgses would extend the reach at the LHC and provides nice complementarity to conventional Higgs search channels.


Introduction
The greatest experimental triumph of the Large Hadron Collider (LHC) till date is the discovery of a scalar resonance at 126 GeV with properties consistent with that of the Standard Model (SM) Higgs [1][2][3][4]. The mass of this particle along with its spin [2,4,5] has now been established, and a complete characterization of all its possible decay modes is underway. At the same time, from the theoretical front, we have now known for a while that the SM, though in excellent agreement with experiments, has to be supplanted with other dynamics if it is to explain many puzzles facing particle physics today, viz., the hierarchy problem, neutrino masses, and the nature of dark matter, to name a few. Many beyond the SM scenarios are constructed to explain one or many of these puzzles, and are becoming more constrained by the Higgs observation at the LHC. This is particularly true for theories constructed with an extended Higgs sector. Well known examples are the Minimal Supersymmetric Standard Model (MSSM) [6][7][8], Next to Minimal Supersymmetric Standard Model (NMSSM) [9,10] and Two Higgs Doublet Models (2HDM) [11][12][13][14]. In addition to the SM-like Higgs boson in these models, the low energy spectrum includes other CP-even Higgses, CP-odd Higgses, as well as charged ones.
Models with an extended Higgs sector hold a lot of phenomenological interest. The discovery of extra Higgses would be an unambiguous evidence for new physics beyond the SM. Other than the decay of these extra Higgses into the SM final states γγ, ZZ, W W , bb and τ τ , which have been the focus of the current Higgs searches, the decay of heavy Higgses into light Higgses, or Higgs plus gauge boson final states could also be sizable. Such decays are particularly relevant as the 126 GeV resonance could show up as a decay of a heavier state, opening up the interesting possibility of using the SM-like Higgs to discover its heavier counterparts. It is thus timely to study these exotic Higgs decay channels and fully explore the experimental discovery potential for the enlarged Higgs sector.
In this paper, we focus on the decays H → AZ or A → HZ, with H and A referring to generic CP-even and CP-odd Higgs, respectively 1 . We consider leptonic decays of the Z, with the A/H in the final states decaying to either a pair of fermions (bb or τ τ ) or ZZ and explore the exclusion bounds as well as discovery reach at the LHC for various combinations of (m A , m H ).
In the 2HDM or NMSSM, both decays H i → A j Z and A i → H j Z could appear with large branching fractions as shown in [15][16][17][18]. Ref. [19] also argued that A → h 0 Z could have a sizable branching fraction in the low tan β region of the MSSM with the light CPeven h 0 being SM-like. A brief Snowmass study of A/H → HZ/AZ with bb final state can be found in Ref. [20]. Another Snowmass study of heavy Higgses [21] explored sensitivities in the H 0 → ZZ → 4 and A → Zh 0 → bb , τ τ channels at the 14 TeV and 33 TeV LHC, focusing on the case with h 0 being the 126 GeV Higgs. In our study, we consider a variety of daughter Higgs masses in bb and τ τ channels, and analyze A → H 0 Z → ZZZ in addition. We also interpret the search results in the context of the Type II 2HDM.
The paper is organized as follows. In Sec. 2, we present a brief overview of models and parameter regions where the channels under consideration can be significant. In Sec. 3, we summarize the current experimental search limits on heavy Higgses. In Sec. 4.1, we present the details of the analysis of the HZ/AZ with the bb final states. We also show modelindependent results of 95% C.L. exclusion as well as 5σ discovery limits for σ × BR(gg → A/H → HZ/AZ → bb ) at the 14 TeV LHC with 100, 300 and 1000 fb −1 integrated luminosity. In Secs. 4.2 and 4.3, we present the analysis for the τ τ and ZZZ final states, respectively. In Sec. 5, we study the implications of the collider search limits on the parameter regions of the Type II 2HDM. We conclude in Sec. 6.

Scenarios with large H → AZ or A → HZ
In the 2HDM, we introduce two SU(2) doublets Φ i , i = 1, 2: where v 1 and v 2 are the vacuum expectation values of the neutral components which satisfy the relation: v 2 1 + v 2 2 = 246 GeV after electroweak symmetry breaking. Assuming a discrete Z 2 symmetry imposed on the Lagrangian, we are left with six free parameters, which can be chosen as four Higgs masses (m h , m H , m A , m H ± ), the mixing angle α between the two CP-even Higgses, and the ratio of the two vacuum expectation values, tan β = v 2 /v 1 . In the case in which a soft breaking of the Z 2 symmetry is allowed, there is an additional parameter m 2 12 . The mass eigenstates contain a pair of CP-even Higgses: h 0 , H 0 , one CP-odd Higgs, A and a pair of charged Higgses H ±2 : Two types of couplings that are of particular interest are ZAH 0 /h 0 couplings and H 0 /h 0 V V couplings, with V being the SM gauge bosons W ± and Z. Both are determined by the gauge coupling structure and the mixing angles. The couplings for ZAH 0 and ZAh 0 are [22]: with g being the SU(2) coupling, θ w being the Weinberg angle and p µ being the incoming momentum of the corresponding particle.
The H 0 V V and h 0 V V couplings are: . (2.4) Note that A always couples to the non-SM-like Higgs more strongly. If we demand h 0 (H 0 ) to be SM-like, then | sin(β − α)| ∼ 1 (| cos(β − α)| ∼ 1) is preferred, and the ZAH 0 (ZAh 0 ) coupling is unsuppressed. Therefore, in the h 0 -126 case, A is more likely to decay to H 0 Z than h 0 Z, unless the former decay is kinematically suppressed. H 0 → AZ could also be dominant once it is kinematically open. Particularly for a heavy H 0 , as we will demonstrate later in Sec. 5, H 0 → AZ can have a large branching fraction in the sin(β − α) = ±1 regions. On the contrary, for H 0 being SM-like with | cos(β − α)| ∼ 1, A → h 0 Z dominates over H 0 Z channel. For very light m A , h 0 → AZ could also open. The detectability of this channel, however, is challenging given the soft or collinear final decay products from a light A. Therefore, for our discussion below, we will focus on the cases A → h 0 Z, H 0 Z and H 0 → AZ only.
In the generic 2HDM, there are no mass relations between the pseudoscalar and the scalar states. Thus, the decays A → h 0 Z, H 0 Z and H 0 → AZ can happen in different regions of parameter spaces. It was shown in Ref. [23] that in the Type II 2HDM with Z 2 symmetry, imposing all experimental and theoretical constraints still leaves sizable regions in the parameter space. In those parameter spaces, such exotic decays can have unsuppressed decay branching fractions. It was also pointed out in Ref. [11] that in the Type I 2HDM, for cos 2 (α − β) > 1/2, the decay h 0 → AZ will actually dominate the W W decay for a light A. Results obtained in this study can also be applied to the CP-violating 2HDM in which H i → H j Z could be sizable with H i,j being mixtures of CP-even and CP-odd states. Appropriate rescaling of the production cross sections and decay branching fractions is needed to recast the results.
The Higgs sector in the MSSM is more restricted, given that the quartic Higgs couplings are fixed by the gauge couplings and the tree-level Higgs mass matrix only depends on m A and tan β. In the usual decoupling region with large m A , the light CP-even Higgs h 0 is SMlike while the other Higgses are almost degenerate: m H 0 ∼ m A ∼ m H ± . Thus, A → ZH 0 or H 0 → ZA is not allowed kinematically. A → Zh 0 is typically suppressed by the small coupling: cos(β − α) ∼ 0, and is only relevant for small tan β. In the NMSSM, the Higgs sector of MSSM is enlarged to include an additional singlet. It was shown in Ref. [17] that there are regions of parameter space where the decay A i → H j Z can be significant.

Current experimental limits
Searches for the non-SM like Higgses, mainly in the bb, µµ, τ τ or W W/ZZ channels have been performed both by ATLAS and CMS. No evidence for a neutral non-SM like Higgs was found.
Searches for the neutral Higgs bosons Φ of the MSSM in the process pp → Φ → µ + µ − /τ + τ − have been performed by the ATLAS [24], and in the τ + τ − channel at CMS [25]. Limits in the µµ channel are much weaker given the extremely small branching fraction in the MSSM. The production mechanisms considered were both gluon fusion and bb associated production, and the exclusion results were reported for the MSSM m max h scenario. The ATLAS study was performed at √ s = 7 TeV with 4.7 -4.8 fb −1 integrated luminosity looking at three different possible τ τ final states, τ e τ µ , τ lep τ had , and τ had τ had . The ATLAS search rules out a fairly sizable portion of the MSSM parameter space, extending from about tan β of 10 for m A ∼ 130 GeV, to tan β ≈ 60 for m A = 500 GeV. The corresponding exclusion in σ Φ × BR(Φ → τ τ ) extends from roughly 40 pb to 0.3 pb in that mass range. The CMS study was performed with 19.7 fb −1 integrated luminosity at 8 TeV and 4.9 fb −1 at 7 TeV in the τ e τ µ , τ µ τ µ , τ lep τ had , and τ had τ had final states. The search excludes roughly between tan β of 4 for m A = 140 GeV and tan β ≈ 60 for m A = 1000 GeV. The corresponding exclusion in σ Φ × BR(Φ → τ τ ) extends from roughly 2 pb to 13 fb in that mass range. In Fig. 1, we recast the current 95% C.L. limit of pp → Φ → τ + τ − in the (m A , tan β) parameter space of the Type II 2HDM [25] (left panel) and the projected 5σ reach at the 14 TeV LHC with 30 fb −1 luminosity [26] (right panel). In both plots, the solid black curves correspond to the limits in the MSSM, when m A ≈ m H 0 with both A and H 0 contributing to the signal. The solid red curves correspond to the limits in the type II 2HDM, when only contribution from A is included and H 0 is decoupled. The reach is considerably weaker: the current exclusion is about tan β ∼ 12 at m A = 160 GeV, and tan β ∼ 46 for m A = 600 GeV. At the 14 TeV LHC with 30 fb −1 luminosity, the 5σ reach extends beyond the current exclusion for large m A . Dashed lines indicate the reduced reach in the τ τ channel once A → h 0 Z mode opens, for a benchmark point of sin(β − α) = 0, m h 0 = 50 GeV and m H 0 = 126 GeV.
Searches with bb final states have also been performed for the MSSM Higgs in the associated production pp → bΦ + X. The CMS search, done with 2.7 − 4.8 fb −1 of data at √ s = 7 TeV excludes tan β values between 18 and 42 in the mass range 90 GeV < m A < 350 GeV [27].
The ATLAS collaboration has also looked for the heavier CP-even Higgs in the Type I and Type II 2HDM, assuming the lighter CP-even Higgs is the discovered 126 GeV boson [28]. The study was performed with 13 fb −1 integrated luminosity at 8 TeV and considered both gluon fusion and vector boson fusion production. Searches in the process H 0 → W W → eµν e ν µ exclude a significant region of the m H 0 − cos α parameter space in the mass range 135 GeV< m H 0 < 200 GeV for the Type II 2HDM. The excluded region shrinks for higher tan β due to the reduced branching ratio to W W . This would serve as a useful constraint if we were to look at decays of the relatively light H 0 to light A's. In this paper, we consider values of m H outside this mass range so this constraint does not apply.
The CMS collaboration has also searched for the heavier CP-even Higgs H 0 and a heavy CP-odd Higgs A in 2HDM via the processes gg → A → h 0 Z and gg → H 0 → h 0 h 0 , assuming the lighter Higgs h 0 is the discovered 126 GeV boson [29]. The study was performed with 19.5 fb −1 integrated luminosity at 8 TeV. Various possible decays of the SM-Higgs were taken into account. Assuming SM branching ratios for h 0 , this study gives an upper bound on σ × BR(A → h 0 Z) of roughly 1.5 pb for m A between 260 and 360 GeV and σ × BR(H 0 → h 0 h 0 ) between 8 pb and 6 pb for masses m H 0 between 260 GeV and 360 GeV. The corresponding excluded parameter space for the Type II 2HDM in the tan β − cos(β − α) plane was also analyzed. In the analysis presented in this paper, we do not necessarily require that the daughter Higgs in A → HZ to be the SM-like Higgs or have SM-like branching ratios. Furthermore we also analyze the process H → AZ for light A and its implication in the Type II 2HDM.

Collider analysis
In this section, we will present model independent limits on the σ × BR for both 95% C.L. exclusion and 5σ discovery for A/H → HZ/AZ in the various final states of bb , τ τ and ZZZ(4 2j). In this study we focus on the leptonic decay of the Z, which allows precise mass reconstruction and suppresses the background sufficiently. Other decay modes of the Z, for example Z → τ τ , might be useful in studying this channel as well. In the discussion of the analyses and results below, we use the decay A → HZ for m A > m H + m Z as an illustration. Since we do not make use of angular correlations, the bounds obtained for A → HZ apply to H → AZ as well with the values of m A and m H switched.

A/H → HZ/AZ → bb
We start our analysis by looking at the channel A/H → HZ/AZ → bb for = e, µ, focusing only on the gluon fusion production channels. We use H to refer to either the light or the heavy CP-even Higgs. Since the only allowed couplings are of the type H − A − Z, if the parent particle is a scalar H, the daughter particle is necessarily a pseudoscalar A and vice versa.
The dominant SM backgrounds for bb final states are Z/γ * bb with leptonic Z/γ * decay, tt with leptonically decaying top quarks, ZZ → bb , and H SM Z [30][31][32][33]. We have ignored the subdominant backgrounds from W Z, W W , H SM → ZZ, W bb, Multijet QCD Background, Zjj, Z as well as tW b. These backgrounds either have small production cross sections, or can be sufficiently suppressed by the cuts imposed. We have included H SM Z here even if the cross section is very small because it has the same final state as the process under consideration, especially for the A → H SM Z case. The total cross sections for these backgrounds can be found in Table 1.
We use Madgraph 5/MadEvent v1.5.11 [34] to generate our signal and background events. These events are passed to Pythia v2.1.21 [35] to simulate initial and final state radiation, showering and hadronization. The events are further passed through Delphes 3.09 [36] with the Snowmass combined LHC detector card [37] to simulate detector effects.
For the signal process, we generated event samples at the 14 TeV LHC for gg → A → HZ with the daughter particle mass fixed at 50, 126, and 200 GeV while varying the parent particle mass in the range of 150 − 600 GeV. We applied the following cuts to identify the signal from the backgrounds 3 : 1. Two isolated leptons, two tagged b's: For jet reconstruction, the anti-k T jet algorithm with R = 0.5 is used.
2. Lepton trigger [38]: 3. Dilepton mass m : We require the dilepton mass to be in the Z-mass window: 80 GeV < m < 100 GeV. 4. m bb versus m bb : We require the dijet mass m bb to be close to the daughter-Higgs mass m H and the mass m bb to be close to the parent-Higgs mass m A . These two invariant masses are correlated, i.e., if we underestimate m bb we also underestimate m bb . To take this into account we apply a two-dimensional cut: (4.4) where w bb × m H is the width of the dijet mass window. Note that the slightly shifted reconstructed Higgs mass m bb (0.95 m H instead of m H ) is due to the reconstruction of the b-jet with a small size of R = 0.5. The second condition describes two lines going through the points (m A ± w bb , m H ) with slope (m Z + m H )/m A . We choose a width for the m bb peak of w bb = Max(Γ H SM | m A , 0.075m A ) where Γ H SM | m A is the width of a SM Higgs with mass m A [39]. This accounts for both small Higgs masses for which the width of the peak is caused by detector effects and large Higgs masses for which the physical width dominates.
The cuts given in Eq. (4.5) follow from simple relativistic kinematics applied to the process as applicable to the entire momenta, i.e., b jets assuming that the parent Higgs A is at rest. We have chosen to specialize this formula to the transverse part alone, including an optimization factor of 0.6. In Fig. 3, we show how this p T cut helps in extracting the signal over the backgrounds for the case where the parent mass is 500 GeV and the daughter mass is 126 GeV. The regions of the plot to the left of the two lines are excluded. It can be seen that while the signal is largely intact, a good portion of the backgrounds gets cut out.  GeV and m H = 126 GeV at the 14 TeV LHC. We have chosen a nominal value for σ × BR(gg → A → HZ → bb ) of 100 fb to illustrate the cut efficiencies for the signal process. In the last column, S/ √ B is shown for an integrated luminosity of L = 300 fb −1 .
In Table 1, we show the signal and background cross sections with cuts for signal benchmark point of m A = 300 GeV and m H = 126 GeV at the 14 TeV LHC. We have chosen a nominal value for σ × BR(gg → A/H → HZ/AZ → bb ) of 100 fb to illustrate the cut efficiencies for the signal process. In the last column, S/ √ B is shown for an integrated luminosity of L = 300 fb −1 . Note that for both the signal and the backgrounds, the biggest reduction of the cross sections arises upon demanding exactly two isolated leptons and b jets. In fact, the signal cross section drops from 100 fb to 6.35 fb at this stage. The two b tag efficiencies bring down the cross section by 0.7 2 ≈ 50%. Other contributing factors are leptons and b jets that are either soft or in the forward direction, or non-isolated leptons and b jets. We also remark that the m cut does not have a significant effect on either the signal or the bb and H SM Z backgrounds since these are dominated by the leptons coming from Z, but does have a pronounced effect on the tt background. The second to last row clearly demonstrates the efficacy of the two dimensional cut in the m bb − m bb plane.  In Fig. 4, we display the results at the 14 TeV LHC for 95% C.L. exclusion (left panel) and 5σ discovery (right panel) limits for σ × BR(gg → A → HZ → bb ), which applies for H → AZ as well with m A and m H switched. The blue, red, and green curves correspond to the daughter particle being 50 GeV, 126 GeV, and 200 GeV, respectively. The masses of the daughter particle are chosen such that they represent cases with a light Higgs, a SMlike Higgs, as well as a heavy Higgs that can decay to W W/ZZ. For each mass, we have displayed the results for three luminosities: 100 fb −1 (dashed), 300 fb −1 (solid), and 1000 fb −1 (dot-dashed), with 10% systematic error included [40]. Better sensitivity is achieved for larger m A since the mass cuts on m bb and m bb have a more pronounced effect on SM backgrounds for larger masses. The limit, however, gets worse for the m H = 50 GeV case when m A 400 GeV (blue curves). This is due to the decrease of the signal cut efficiency for a highly boosted daughter particle with two collimated b jets. For the interesting case where the daughter particle is 126 GeV, it is seen that the discovery limits for a 300 fb −1 collider fall from about 0.7 pb for m A of 225 GeV, to less than 20 fb for a 600 GeV parent particle. These numbers do not change appreciably between the three chosen luminosity values, except for the case of m H = 50 GeV and m A 400 GeV. This is because we have chosen a uniform 10% systematic error on the backgrounds, which dominates the statistical errors for most of the parameter region. For a given parent particle mass m A , limits are better for smaller m H = 50 GeV. This is because the m bb distribution for the dominating Zbb and tt backgrounds peaks around higher masses m bb ≈ 70 -200 GeV and therefore the background rejection efficiency for m bb ≈ 50 GeV is high. We reiterate here these exclusion and discovery limits are completely model independent. Whether or not discovery/exclusion is actually feasible in this channel should be answered within the context of a particular model, in which the theoretically predicted cross sections and branching fractions can be compared with the exclusion or discovery limits. We will do this in Sec. 5 using Type II 2HDM as a specific example.

A/H → HZ/AZ → τ τ
We now turn to the process gg → A/H → HZ/AZ → τ τ . Since we want to reconstruct the final state particles unambiguously, we will employ τ tags and thus will only consider fully hadronic τ decays. While the signal is typically suppressed compared to the bb case due to the smaller H → τ τ branching fraction, the SM backgrounds [32,33] are much smaller due to the absence of b jets in the final states. The dominant background is ZZ. We have also included H SM Z background even though it is negligible for most cases.
Here, we list the cuts employed: 1. Two isolated leptons and two tagged τ 's: We do not impose jet veto.
2. Lepton trigger: 3. Dilepton mass m : 80 GeV < m < 100 GeV. 4. m τ τ versus m τ τ : The expected Higgs mass is shifted more towards smaller values compared to the bb case. This is because of the hadronic decay of τ with missing energy carried away by neutrinos. Our 2-D cuts are modified as follows: We show the normalized 2-D distribution as well as cuts imposed as indicated by red lines in Fig. 5 for the signal (left panel) and the backgrounds (right panel). The cut filters out most of the backgrounds while retaining the signal, yielding a good S/ √ B value.

5.
Transverse momentum: The looser cut on τ p T compared to the bb case is again due to the extra missing E T in the τ decay.
In Table 2, we present the cross sections after the individual cut is imposed sequentially. We take a nominal signal cross section of 10 fb to illustrate the efficiency of the chosen cuts. Again, the 2-D m τ τ − m τ τ cut improves the S/ √ B value significantly. In Fig. 6, we show the 95% C.L. exclusion and 5σ discovery reach in σ × BR(gg → A → HZ → τ τ ) for the 14 TeV LHC. The general feature of these plots follows that of Fig. 4, particularly with highly boosted daughter particles making τ identification more challenging, as shown by the blue curves for 50 GeV daughter particle mass, which exhibit worse limits for m A > 400 GeV. The exclusion limits are lowest for small m H = 50 GeV and also for high m H = 200 GeV since the dominating ZZ background peaks at m τ τ ≈ 90 GeV and therefore our m τ τ mass cut leads to a high background rejection for lower or higher m H . Since the statistical error dominates the 10% systematic error, the σ × BR limits scale roughly with 1/ √ L, as indicated by the dashed, solid and dot-dashed lines for different luminosities.
Compared to the bb case, the σ × BR reach in τ τ case is better due to significantly lower SM backgrounds. For the 126 GeV daughter particle case with 300 fb −1 , the 5σ discovery reach varies from about 20 fb for parent mass of 225 GeV to about 3 fb for 600 GeV. Thus, given the typical ratio of Br(H/A → bb) : Br(H/A → τ τ ) ∼ 3m 2 b /m 2 τ , the reach in τ τ can be comparable or even better than bb channel, in particular, for smaller parent Higgs masses.

A → HZ → ZZZ → 4 + 2j
We now consider the case where the daughter particle decays to a pair of Z bosons, which only applies to A → HZ → ZZZ. This process involves a trade-off between having a clean final state with suppressed backgrounds and suppressed signal cross section for detection. We find that the best final states combination that yields signal cross sections that are not too suppressed in realistic models with controllable backgrounds is the 4 + 2j final state: A → HZ → ZZZ → 4 + 2j. The SM backgrounds for this process come from the single, double and triple vector boson processes including additional jets as well as tt background [37,41,42].
Note that the Z's from the H decay could be either on-shell or off-shell depending on m H . We will display our results for two cases: one where one of the final state Z's is necessarily off-shell, and another where both are on-shell. We will find that the latter case leads to much better discovery prospects.
• Three Z-candidates: We reconstruct the hadronically decaying Z using the 2 hardest jets. To reconstruct the leptonically decaying Z's: -4e or 4µ: If we have 4e or 4µ, we first find the combination of electrons or muons with opposite charge that is closest to the Z-mass. The other 2 electrons or muons are combined to find the last Z.
-2e2µ: Here, we combine the same flavored leptons in a straightforward manner.
• Z masses: We require the hardonically decaying Z 1 , the well reconstructed leptonically decaying Z 2 and the final reconstructed leptonically decaying Z 3 to be in the following windows: 60 GeV <m Z 1 < 115 GeV.
m min <m Z 3 < 115 GeV. (4.12) Here, we assume Z 1 to be on-shell. However, we allow for the possibility that Z 3 could be far off-shell. The m min employed here mimics the LHC search strategy for the SM Higgs, and its value depends on the Higgs mass and can be found in Table 2 of Ref. [43].
• m H and m A : The Z produced in the A decay typically has a higher p T than the Z's produced in H decay. Therefore we assume that the lower p T Z's are coming from the H. For the reconstructed H with mass m ZZ and A with mass m ZZZ we require: (4.14) In Table 3 We note that the nominal value for the cross section that is used in Table 3 can, in typical BSM scenarios, be enhanced at small tan β, due to the top loop contributions to the gluon fusion production, as well as the suppression of the H → bb branching fraction. Fig. 7 shows the 95% C.L. exclusion and 5σ discovery at the 14 TeV LHC for different integrated luminosities: L = 100 fb −1 , 300 fb −1 , and 1000 fb −1 . Even for L = 300 fb −1 ,   the discovery limits vary only between about 3 fb and 1.5 fb with 200 GeV m H for m A between 300 GeV and 600 GeV. Thus, the only challenge in this channel is to have high enough signal cross sections, as the SM backgrounds prove to be less of a threat compared to the bb final state.

Implications for the Type II 2HDM
The decays A/H → HZ/AZ appear in many models that have an extension of the SM Higgs sector. In this section, we illustrate the implications of the exclusion or discovery limits of bb , τ τ and ZZZ(4 2j) searches on these models using Type II 2HDM as an explicit example.
In the Type II 2HDM, one Higgs doublet Φ 1 provides masses for the down-type quarks and charged leptons, while the other Higgs doublet Φ 2 provides masses for the up-type quarks. The couplings of the CP-even Higgses h 0 , H 0 and the CP-odd Higgs A to the SM gauge bosons and fermions are scaled by a factor ξ relative to the SM value, which are presented in Table 4. Table 4. The multiplicative factors ξ by which the couplings of the CP-even Higgses and the CPodd Higgs to the gauge bosons and fermions scale with respect to the SM value. The superscripts u, d, l and V V refer to the up-type quarks, down-type quarks, leptons, and W W/ZZ respectively.
The implication of the current Higgs search results on the Type II 2HDM has been studied in the literature [15,16,18,23,44]. In particular, a detailed analysis of the surviving regions of the Type II 2HDM was performed in [23], considering various theoretical constraints and including the latest experimental results from both the ATLAS and the CMS. Either the light or the heavy CP-even Higgs can be interpreted as the observed 126 GeV SM-like Higgs, with very different preferred parameter regions. In the h 0 -126 case, we are restricted to narrow regions with sin(β − α) ∼ ± 1 with tan β up to 4 or an extended region in 0.55 < sin(β −α) < 0.9 with 1.5 < tan β < 4. The masses m H 0 , m H ± , and m A are, however, relatively unconstrained. In the H 0 -126 case, we are restricted to a narrow region of sin(β − α) ∼ 0 with tan β up to about 8, or an extended region of sin(β − α) between −0.8 to −0.05, with tan β extending to 30 or higher. m A and m H ± are nearly degenerate due to ∆ρ constraints. Imposing the flavor constraints in addition further narrows down the preferred parameter space.
Given the different parameter dependence of the gluon fusion cross section for A and H 0 , the branching fractions of h 0 , H 0 and A, as well as the coupling difference between h 0 AZ and H 0 AZ, we can identify three different classes of processes: gg → A → h 0 Z, gg → A → H 0 Z, and gg → H 0 → AZ when interpreting the exclusion and discovery limits from the previous sections. We do not consider the decay of h 0 → AZ since this channel is experimentally challenging given that both h 0 and A are relatively light. Table 5. Benchmark points shown for illustrating the discovery and exclusion limits in the processes considered in the context of the Type II 2HDM. The checkmarks indicate kinematically allowed channels. Also shown are the typical favored region of sin(β − α) for each case (see Ref. [23]).
In Table. 5, we list the benchmark points that we use for the interpretation of the exclusion and discovery bounds in the Type II 2HDM.  Table 5 are the preferred regions in sin(β − α) once all the theoretical and experimental constraints are imposed, following Ref. [23].
Note that in our study, we have decoupled the charged Higgs so that it does not appear in the decay products of A or H. For a light charged Higgs that is accessible in the decays of A/H → H ± W ∓ , H + H − , decay branching fractions of A/H → HZ/AZ will decrease correspondingly, which reduces the reach of this channel. However, the new decay channels involving the charged Higgs might provide new discovery modes for A or H, which have been explored elsewhere [45][46][47]. In particular, for A/H → H ± W ∓ , H + H − with H ± → τ ± ν, the spin correlation in the τ decay can be used to identify the signal from the SM backgrounds. The sensitivity of this channel involving H ± in the intermediate to large tan β region provides a nice complementarity to the A/H → HZ/AZ channels [45].
To be more general, in the discussion below when we interpret the search results of bb , τ τ and ZZZ(4 2j) channels in the model parameter space, we do not restrict ourselves to the narrow preferred parameter regions for h 0 -126 or H 0 -126 case as shown in Ref. [23]. In particular, we consider the broad range of −1 ≤ sin(β − α) ≤ +1 and 1 ≤ tan β ≤ 50. This is because the allowed regions would change if a soft Z 2 symmetry breaking is incorporated which Ref. [23] did not deal with. Furthermore, the Higgs sector of 2HDM and the subsequent symmetry breaking structure is rather general and the results presented in this section can be interpreted in the context of any such model if the Higgs couplings to the fermions follow a similar pattern. We do, however, point out the interplay between the exotic Higgs decay channels and the SM-like Higgs search results at the end of each discussion.

gg → A → h 0 Z
We compute the production cross section for the CP-odd Higgs A by a simple rescaling of the SM Higgs cross section as follows: where τ f = 4m 2 f /m 2 A and the scalar and pseudoscalar loop factors F h 1/2 and F A 1/2 are given by: [22] F and with η ± ≡ 1 ± √ 1 − τ . We have ignored the contribution from other Higgses in the loop, which is typically small. The left panel of Fig. 8 shows the contour plot of the σ(gg → A) normalized to that of the SM Higgs with the same mass. The tan β dependence is due to the Att and Abb couplings, while the mass dependence comes from the different dependence of F 1/2 (τ f ) on τ f for pseudoscalar compared to a scalar. Enhancements over the SM value is possible for large tan β at small m A due to the bottom loop, or small tan β for all values of m A due to the top loop. The bump in the plot for m A around 350 GeV corresponds to top threshold effects. Note that for A, the production cross section only depends on tan β and is independent of α. Also shown in the right panel of Fig. 8 are contours of σ(gg → A) in the m A − tan β plane for the 14 TeV LHC, with the cross sections for the SM Higgs production obtained from Ref. [39,48]. Significant cross sections of 10 pb or more are possible for large m A up to 500 GeV for small tan β. Cross sections of similar magnitude are also possible at large tan β due to the bottom loop enhancement effects, albeit only for relatively small m A .
In Fig. 9, we show contour plots of BR(A → h 0 Z) for BP1 (left panel) and BP3 (right panel). BR(A → h 0 Z) always maximizes at sin(β − α) = 0, and decreases for larger | sin(β − α)|, since g ZAh 0 ∼ cos(β − α). For BP1 with (m A , m H 0 , m h 0 ) = (400, 126, 50) GeV, both A → h 0 Z and A → H 0 Z open, with the coupling of the latter process proportional to sin(β − α). Therefore, BR(A → h 0 Z) decreases more rapidly when | sin(β − α)| gets bigger. BR(A → h 0 Z) decreases at large tan β as A → bb becomes more and more important. For BP3 with (m A , m H 0 , m h 0 ) = (300, 400, 126) GeV, only A → h 0 Z opens with no competitive process from A → H 0 Z and A → tt. Therefore, comparing to BP1, BR(A → h 0 Z) decreases much slower as sin(β − α) approaches ±1. BR(A → h 0 Z) is also maximized at smaller tan β due to both the absence of A → tt and the suppression of A → bb.
To compare with the exclusion and discovery limits in the bb , τ τ channels, it is also important to know the branching fractions of h 0 → bb, τ τ , which depend mostly on m h 0 . For BP1 with m h 0 = 50 GeV, we used BR(h 0 → bb)= 82% and BR(h 0 → τ τ )= 8%. For the other benchmark points with h 0 being the SM-like 126 GeV Higgs, the branching fraction is obtained by rescaling the SM value of the BR with relevant coupling coefficients as given in Table. 4. We show a contour plot of BR(h 0 → bb) in Fig. 10 for h 0 being the 126 GeV Higgs. While h 0 → bb reaches 80% and saturates in most of the parameter space, there is a wedge shaped region around 0.5 < sin(β − α) < 1 at small tan β in which h 0 → bb could be suppressed.
In Fig. 11, we show the LHC 100 fb −1 discovery/exclusion reach for gg → A → h 0 Z in the bb (red curves) and τ τ (blue curves) channels for BP1 (left panel), BP2 (middle panel) and BP3 (right panel). 95% Exclusion regions are shown as yellow regions enclosed For BP2 with m A = 400 GeV and m h 0 = 126 GeV, regions of tan β < 10 or tan β > 32 will be excluded if no signal is detected, and regions of tan β < 4 can be discovered if there are positive signals. For BP3 with m A = 300 GeV and m h 0 = 126 GeV, the exclusion and discovery regions shrink further at small tan β. The wedge-shaped region toward sin(β − α) = 1 corresponds to the wedge region in Fig. 10. Our results agree with that of Ref. [21] for A → h 0 Z with h 0 being the SM-like Higgs.
We note the interesting feature that the bb limits are better than the τ τ ones for BP1 and BP2, while the behavior flipped for BP3. This is because τ τ typically has better reach than bb process at small m A , while bb does better at large m A , when the BR(h 0 → bb)/BR(h 0 → τ τ ) ∼ 3m 2 b /m 2 τ is taken into account. Given the smallness of the branching fraction of h 0 → ZZ for the m h 0 values chosen, the ZZZ channel will not be useful in probing the parameter space with gg → A → h 0 Z. We also note that for the H 0 -126 case (BP1) with the favored region to interpret H 0 as the SM-like Higgs being around sin(β − α) ∼ 0, gg → A → h 0 Z will be extremely useful in probing this region. For the h 0 -126 case (BP2 and BP3), the favored region to interpret h 0 as the SM-like Higgs is around sin(β − α) = ±1. Even though the A → h 0 Z branching ratio is typically suppressed when sin(β − α) approaches ±1, we could still have reach in sin(β − α) extending fairly close to ±1. In the left panel of Fig. 12, we show the reach in tan β versus m A plane for m h 0 = 50 GeV and sin(β − α) = 0, with 95% C.L. exclusion (yellow regions enclosed by the solid curves) and 5σ discovery (cyan regions enclosed by dashed curves) given for bb channel (red lines) and τ τ channel (blue lines). While τ τ is more sensitive at low m A , bb extends the reach at large m A . In general, small tan β (lower region) or large tan β (top region) are within reach due to the enhancement of the top and bottom Yukawa couplings in those regions. For small tan β ∼ 1, almost all values of m A up to 600 GeV can be covered, with regions of m A shrink for increasing tan β. At large tan β 10, small m A can not be approached due to the weakening of the experimental limit, while large m A can not be approached due to the decreasing of the signal cross sections.
In the right panel of Fig. 12, we show the reach in m A − tan β plane for m h 0 = 126 GeV and sin(β − α) = 0.6. Note that we have chose a value for sin(β − α) that is consistent with the current Higgs search results [23] of a 126 GeV h 0 while still allowing a sizable branching fraction for A → h 0 Z. We have decoupled the heavy CP-even Higgs H 0 so that A → H 0 Z does not occur. Given the reduced branching fraction for A → h 0 Z, as well as the worse exclusion/discovery limits, the exclusion and discovery regions are smaller, compared to the left panel with m h 0 = 50 GeV, sin(β − α) = 0. In particular, only regions with tan β 8 or a small region in tan β 50 around m A ∼ 450 GeV are viable.

gg → A → H 0 Z
A → H 0 Z opens once it is kinematically accessible. Since m h 0 < m H 0 , A → h 0 Z is always accessible and more favorable in phase space. Whether A → H 0 Z dominates or not depends largely on sin(β − α), which controls the coupling of ZAH 0 as well as ZAh 0 . Fig. 13 shows the contours of BR(A → H 0 Z) in the parameter space of tan β versus sin(β − α), for BP1 in the left panel. Contrary to the A → h 0 Z case as shown in Fig. 9, the branching ratios become larger for larger | sin(β − α)|, which is maximized at sin(β − α) = ±1, consistent The branching fraction is less than 10% over almost the entire parameter space. It is also evident that unlike BP1 and BP2, there is no suppression of the branching fractions at small tan β due to the absence of the tt decay mode.  0.2 for small tan β. There is also a small additional bump around sin(β − α) = −0.6, mainly due to the increasing of BR(H 0 → bb), as shown in the left panel of Fig. 14. The reach is greatly reduced for BP2 due to the suppression of H 0 → bb, except for sin(β − α) ∼ ±1. Only thin slices of parameter region near sin(β − α) ∼ ±1 can be covered, which extends to tan β 8 for the exclusion, and tan β 4.5 for discovery.
Note that for BP1 with (m A , m H 0 , m h 0 ) = (400, 126, 50), both A → h 0 Z and A → H 0 Z open. The former is more sensitive to the sin(β − α) ∼ 0 region, as shown in the left panel of Fig. 11, while the latter is more sensitive to sin(β − α) ∼ ±1, as shown in the left panel of Fig. 15. Searches in these two channels are complementary to each other. When combined, they could cover the entire region of sin(β − α), in particular, for tan β 10. Note that when combined with the current experimental search results for the 126 GeV Higgs being the H 0 , the region with sin(β − α) ∼ 0 is favored, with a thin slice of extended region at negative −0.8 < sin(β − α) < −0.05 as well [23].
Similar complementarity between A → h 0 Z and A → H 0 Z can be found for BP2 with (m A , m H 0 , m h 0 ) = (400, 200, 126) GeV, for the entire region of sin(β − α). Interpreting h 0 being the 126 GeV observed Higgs boson, furthermore, favors sin(β − α) ∼ ±1 or a thin slice of extended region at 0.55 sin(β − α) 0.9 [23]. In the left panel of Fig. 16, we present the exclusion and discovery reach in the tan β versus m A plane for A → H 0 Z with m H 0 = 126 GeV, m h 0 = 50 GeV and sin(β −α) = −0.8. We have chosen the value of sin(β − α) such that the branching faction of A → H 0 Z is sizable while still consistent with the experimental Higgs search results [23] with a 126 GeV H 0 . We see that tan β up to about 6.5 can be reached for exclusion, and tan β up to about 3.5 can be reached for discovery.
In the right panel of Fig. 16, we present the exclusion and discovery reach in the tan β versus m A plane for m H = 200 GeV, m h 0 = 126 GeV and sin(β − α) = 1. For 350 GeV m A 600 GeV, tan β up to about 6 can be excluded, and up to about 3 can be discovered in the bb channel. τ τ channel does better in the low m A region.
For BP2 with m A = 400 GeV, m H 0 = 200 GeV, we can also study the parameter reach of A → H 0 Z with H 0 → ZZ. In Fig. 17, we show BR(H 0 → ZZ) in the left panel, which reaches a maximum of 25% for | sin(β − α)| 0.2. It gets larger for small tan β when H 0 → bb is further suppressed. In the right panel of Fig. 17, we show the discovery and exclusion contours in the tan β versus sin(β − α) plane for 100 fb −1 luminosity at the LHC. While H 0 → ZZ maximizes at sin(β − α) ∼ 0, A → H 0 Z is minimized in this region. As a result, regions of 0.3 | sin(β − α)| 1 with tan β up to 4.7 can be excluded while the discovery regions are 0.5 | sin(β − α)| 1 with tan β 2.8. Note also that this channel is complementary to A → H 0 Z → bb/τ τ as shown in Fig. 15, which is sensitive to sin(β − α) ∼ ±1 region.  In Fig. 18, we present the exclusion and discovery in tan β versus m A plane for gg → A → H 0 Z → ZZZ(4 2j) with m H = 200 GeV, sin(β − α) = 0.9. We have chosen the value of sin(β − α) such that the branching fractions of both A → H 0 Z and H 0 → ZZ is sizable while still consistent with the experimental Higgs search results [23] with a 126 GeV h 0 . We see that the whole region of 300 GeV < m A < 600 GeV can be covered at small tan β, with the maximum reach in tan β obtained for m A ∼ 350 GeV: tan β 3 for discovery and tan β 5 for exclusion.

gg → H 0 → AZ
For this process, we restrict to the m h 0 = 126 GeV case with a heavier H 0 . We use BP4 with (m A , m H 0 , m h 0 ) = (50, 400, 126) GeV and BP5 with (m A , m H 0 , m h 0 ) = (200, 400, 126) GeV as an illustration. The gluon fusion production cross section for H 0 can be rescaled from the SM cross section: where the loop factors F 's are defined in Eq. (5.2). We note that in contrast to the production of A in Eq. (5.1), the production of H 0 involves both α and β. In the left panel of Fig. 19, we show contours of the production cross section of H 0 normalized to the SM value in the sin(β − α) − tan β plane for m H 0 = 400 GeV. We see that for positive sin(β − α), the cross section is always relatively more suppressed than that for negative sin(β − α), introduced by the interference between the top and bottom loops in Eq. (5.4).
For sin(β − α) = ±1, which is preferred by the interpretation of h 0 being the SM-like Higgs, the cross section receives the strongest suppression: only 10% of the corresponding SM value. In the right panel of Fig. 19, we show contours of the production cross section at 14 TeV LHC in the m H 0 − tan β plane. We see that cross sections of 10 pb or more is possible for m H 0 up to 425 GeV for small tan β -slightly lower than the corresponding numbers for σ(gg → A) as shown in Fig. 8. However, the bottom loop enhancement plays a slightly more significant role in this case at large tan β, compared to the A case.  Since g ZAH 0 ∝ sin(β − α), the branching fraction gets bigger for larger | sin(β − α)|, and is maximized at sin(β − α) = ±1. Branching fractions in BP4 is larger than that of BP5 due to the bigger phase space for H 0 → AZ. For A → bb and τ τ , the branching fraction is about 94% and 6% respectively, which does not vary much for BP4 with m A = 50 GeV and BP5 with m A = 200 GeV. tan β 26 can be excluded and a smaller region in −1 sin(β − α) −0.6 with tan β 5 can be discovered. While bb channel has better reach for BP4, τ τ channel has a slightly better sensitivity for BP5. The reach is also much better for negative sin(β − α) because of the less suppressed cross sections of gg → H 0 . In the left panel of Fig. 22, we show the exclusion and discovery each with 100 fb −1 luminosity at 14 TeV LHC in tan β versus m H plane, for gg → H 0 → AZ with bb (red) and τ τ (blue) final states. We have chosen m A = 50 GeV and sin(β − α) = −1. Discovery is possible for small values of tan β 5 or larger values of tan β 20. The exclusion reach, however, is much more extended. All values of tan β can be covered for m H 0 up to 450 GeV, with reach extended further at larger and smaller values of tan β. The reach with daughter particle mass m A = 200 GeV is shown in the right panel of Fig. 22. Both the exclusion and discovery regions shrink greatly. Only very small tan β 4 or very large tan β 44 can be excluded. Note that while sin(β − α) = ±1 is preferred by the interpretation of the h 0 being the SM-like Higgs, the suppression of gg → H 0 in that region results in a reduced exclusion/discovery reach. Even a small deviation of sin(β − α) away from ±1 would introduce a much larger reach in gg → H 0 → AZ.

Conclusion
Given the discovery of a 126 GeV SM-like Higgs boson at the LHC, it is now time to use the experimental data to constrain new physics models while also exploring the detectability of extra Higgs bosons in the extensions of the SM. In this spirit, we explored the production and decay of heavy scalar and pseudoscalar states via the processes gg → H 0 → AZ and gg → A → h 0 Z/H 0 Z with both fermionic (bb, τ τ ) and possible bosonic (ZZ) decays of the daughter Higgs. This channel provides nice complementarity to the conventional search channel pp → A/H → τ τ , which is mostly sensitive to the large tan β region. We presented model independent limits on the 95% C.L. exclusion and 5σ discovery in those channels at the 14 TeV LHC. The possibilities include the interesting case of having the 126 GeV SM-like Higgs as a decay product of a heavy pseudoscalar.
For the 14 TeV LHC with 300 fb −1 integrated luminosity, the 95% C.L. limits on σ×BR for the bb final state (where the b's come from the Higgs in the final state) for a 126 GeV daughter Higgs particle vary between 200 fb to a few fb for the parent heavy Higgs mass in the range of 200 GeV to 600 GeV, while the limit for 5σ discovery is about 3−5 times larger. For the τ τ channel with the same range of A mass, the exclusion bounds are around 5 − 1 fb and the discovery reach is about 20 fb − 3 fb. While the σ × BR reach in the τ τ channel is in general much better than the bb channel, owing mostly to more suppressed backgrounds, it is comparable to bb mode once the branching fraction difference between bb and τ τ modes are taken into account in a given model. gg → A → H 0 Z → ZZZ → 4 2j is useful for heavy Higgses with m H 0 > 2m Z . For m H 0 = 200 GeV and m A = 400 GeV, exclusion in this channel with 300 fb −1 integrated luminosity requires as little as 1 fb in σ × BR while 5σ discovery needs about 3 fb.
We then discussed the implication of the exclusion and discovery bounds of bb , τ τ and ZZZ channels in the Type II 2HDM, studying three classes of processes: gg → A → h 0 Z, gg → A → H 0 Z, and gg → H 0 → AZ. We find, in general, that there is a significant portion of the tan β versus sin(β − α) plane that allows discovery/exclusion possibilities in the bb and τ τ final states. bb and τ τ have comparable reach, with τ τ being slightly better for low parent Higgs masses and bb being better for higher parent Higgs masses.
Specifically, in the channel gg → A → h 0 Z when H 0 is identified as the SM-like Higgs, 95% exclusion covers most of the tan β versus sin(β − α) plane for m A around 400 GeV. tan β < 5 can also be covered by 5σ discovery. On the other hand, the exclusion/discovery range is more restricted when h 0 is identified as the SM Higgs. Typically, we find that for m A = 400 GeV, discovery region lies between −1 < sin(β − α) 0.8 and tan β ≤ 5, while the exclusion region extends to tan β 10 or 30. Note also that even though the reach is always maximized at sin(β − α) ∼ 0, it extends to larger values of | sin(β − α)| close to ±1 as well. A wide range of m A can be covered at low tan β 10, while high tan β can only be approached for m A 500 GeV.
The case where A decays to H 0 Z is complementary to A → h 0 Z in that the discovery and exclusion regions split into two distinct regions around sin(β − α) ∼ ±1. We find that in both the bb and τ τ channels, the discovery reach covers tan β up to about 4, while the exclusion region extends to about 7 for m A up to about 600 GeV. Moreover, for m H 0 ≥ 200 GeV, this channel also allows for an exclusion reach with ZZZ final states with 0.3 < | sin(β − α)| < 1, and tan β up to 4.5 for m A around 400 GeV. For small values of tan β, a wide range of m A can be covered either by exclusion or discovery.
In the last class gg → H 0 → AZ, we find that discovery/exclusion regions favor the negative sin(β − α) regions, largely due to the parameter dependence of gluon fusion production σ(gg → H 0 ). For m H 0 = 400 GeV and m A = 50 GeV, a wide range of tan β versus sin(β−α) space can be covered, except for a small stripe around −0.15 < sin(β−α) < 0.2. For m A = 200 GeV, the regions −1 sin(β − α) −0.5 can be excluded for all values of tan β, while only a smaller region at low tan β can be discovered. For m A = 50 GeV and sin(β − α) = −1, the exclusion reach in m H can be as large as 450 GeV for tan β around 10, which extends even further for smaller and larger tan β.
While extra Higgs bosons other than the observed 126 GeV SM-like Higgs exist in many extension of the SM, the searches for those Higgses in unconventional decay channels have just started. Compared to conventional search channels of bb, τ τ , W W , ZZ and γγ, these exotic decay modes of heavier Higgses decaying into two light Higgses or one Higgs with one gauge boson can be dominant in certain regions of parameter space. In this paper, we explored A/H → HZ/AZ in bb , τ τ and ZZZ modes. Other channels, in particular, those involving charged Higgses can be very promising as well [45][46][47].