Leptonic cascade decays of a heavy Higgs boson through vectorlike leptons at the LHC

We demonstrate the potential of fully leptonic cascade decays of a heavy neutral Higgs boson through vectorlike leptons as a simultaneous probe for extended Higgs sectors and extra matter particles at the LHC. The processes we explore are unique in that their event topologies lead to di-boson-like leptonic final states with a lepton pair which does not reconstruct the mass of a gauge boson. By recasting existing $2\ell + E_T^{\rm miss}$ and $3/4\ell$ searches channels using run2 data from the LHC we obtain $\textit{model independent}$ bounds on the masses of heavy scalars and vectorlike leptons and use these results to explore future prospects at the HL-LHC. Our results can be directly applied to any kind of new physics scenarios sharing the final states and the event topology. For concreteness, we apply our results to a benchmark scenario: a two Higgs doublet model type-II augmented with vectorlike leptons. Remarkably, even with current data the sensitivity of our analysis shows a reach for masses of a heavy neutral Higgs and vectorlike leptons up to 2 TeV and 1.5 TeV, respectively. Even for low $\tan\beta \gtrsim 1,$ the analysis retains sensitivity to heavy Higgs masses slightly above 1 TeV. The future sensitivities at the HL-LHC extend the reach for heavy Higgses and new leptons to 2.7 TeV and 2 TeV, respectively.


Introduction
Since the discovery of the Higgs boson, unveiling the Higgs sector involved in electroweak symmetry breaking (EWSB) has become a strong focus of the LHC program and is essential to understanding possible features of new physics beyond the Standard Model (SM). In particular, there are no symmetry arguments prohibiting the existence of an additional isodoublet field in the Higgs sector. Two Higgs doublet models (2HDMs) are generally considered to be among the simplest and theoretically well-motivated frameworks for possible extensions of the Higgs sector, for a review see, e.g. Refs. [1,2]. Popular examples of new theories beyond the SM (BSM) sharing this 2HDM-like Higgs sector include the Minimal Supersymmetric SM (MSSM) [3,4], the DFSZ axion model [5,6], and Twin Higgs model [7].
Another possibility which is not forbidden by any symmetry arguments is the existence of extra fermions beyond the three generations of the SM. Since the existence of new chiral fermions is strongly disfavored both by theoretical and experimental reasons [8], extra matter fermions, if they exist, are likely to be vectorlike with masses which can be generated independently of the Higgs mechanism. Vectorlike fermions have also been considered in many compelling BSM theories including both non-supersymmetric and supersymmetric models with complete vectorlike families [9][10][11][12][13], composite/little Higgs model [14,15], the KSVZ axion model [16,17], and gauge mediated supersymmetry breaking model [6,[18][19][20][21][22][23]. Furthermore, models that include both an extension of the Higgs sector and vectorlike leptons were recently studied as possible solutions for the discrepancy in the measured value of the muon anomalous magnetic moment [24][25][26][27].
Existing searches at the LHC for extra Higgs bosons or vectorlike fermions have been developed independently. On the other hand, it has been proposed that signatures involving both types of new particles simultaneously have many advantages over the typical searches [28][29][30][31][32][33][34][35][36]. In particular, heavy Higgs cascade decays through a vectorlike lepton mixing with the SM fermion, such as those depicted in Fig. 1, provide relatively clean and distinctive signals due to their unique event topology. The final states of the processes in Fig. 1 are fully leptonic, possibly with missing energy, and hence we dub these signatures leptonic cascade decays throughout this paper. An important feature of our process is that the signature is "di-boson-like" but with one lepton pair whose momentum does not reconstruct the mass of a gauge boson in any way.
In this paper, we obtain new, model-independent constraints on heavy scalars and vectorlike leptons by recasting recent results of the dilepton with missing energy (2 +E miss T ) search [39] and 3 or 4 lepton (3/4 ) search [40] at the LHC run2. We show that, remarkably, even with current data assuming a 50% cascade branching ratio the sensitivity to masses of new scalars and heavy leptons can reach well above 1 TeV, proving the effectiveness of our leptonic cascade processes. Our analysis results channel. The diagrams are drawn by TikZ-FeynHand [37,38].
are quite promising compared to the typical consensus that the sensitivities for new leptons are weaker than those for vectorlike quarks considering the conventional production processes through gauge bosons. 1 For concreteness, we interpret our results in the context of a reference model: a two Higgs doublet model (2HDM) type-II, where Φ = H (A) is the heavy CP even (odd) neutral Higgs boson, augmented by vectorlike leptons which mix with the muon with e 4 (ν 4 ) being the lightest charged (neutral) vectorlike lepton. This is exactly the scenario which can explain the muon anomalous magnetic moment, even with multi-TeV Higgs bosons and new leptons [26,27]. This is the main motivation to pursue signals with muons in the final state. However, results generated from vectorlike leptons mixing with the first generation would be almost identical. We find that the projected sensitivity of the analysis for the HL-LHC extends the reach for heavy Higgses and new leptons to 2.7 TeV and 2 TeV, respectively. We emphasize that the search strategies and results in this paper can be readily applied to other new physics scenarios with the same kinematic topology and final states such as models of vectorlike leptons with a Z boson [43][44][45][46], singlet-extended 2HDM or Next-to-MSSM type models [47]. This paper is organized as follows. In Sec. 2, we outline search strategies at the LHC for the leptonic heavy Higgs cascade decays. In Sec. 3, we obtain current and projected upper limits on the proposed cross sections and compare them to the corresponding limits from conventional heavy Higgs search channels. As part of the comparison we also obtain sensitivities to masses of heavy Higgses and new leptons in our reference model. Section 4 is devoted to our conclusions and outlooks. Details of the 2HDM we consider are explained in Appendix A.

Search strategies
In this section, we provide a detailed description of the analysis strategies we propose to obtain new constraints for heavy Higgses and new leptons based on the full LHC run2 data and projected expectations for the HL-LHC. As part of our analysis we recast the conventional searches for heavy Higgses based on decays to τ lepton pairs, as well as existing searches for leptonic signals in 2 + E miss T and 3/4 final states initially motivated in supersymmetric and seesaw models. Importantly, our recast of the 2 + E miss T and 3/4 searches can be generically applied to other BSM scenarios with the same event topology and final states, whereas our results for the heavy Higgs to di-τ channel are relevant for our reference model of the two Higgs doublet model type-II augmented by vectorlike leptons in the alignment limit.
For simplicity, events for SM background processes are taken from data available from existing results, giving the relevant references along the way. Nevertheless, as we will see our approach provides powerful enough experimental sensitivities on our leptonic cascade processes, as will be shown in Sec. 3. Although not implemented here, further investigation for multiple lepton resonances in the 3/4 channel would increase the sensitivities even more, as already demonstrated in another cascade process in Ref. [31].

Heavy neutral Higgs to τ τ
Searches for heavy BSM Higgs bosons have been considered as important tasks in the LHC collaborations. Among the conventional searches for a heavy neutral Higgs boson decaying to the SM fermions, we impose constraints from the Φ → τ τ searches [48,49] which provides the most stringent limits in a wide range of parameter space (for example, see Ref. [50]).
The neutral Higgs bosons dominantly decay to the SM fermions in the third generation. In our numerical analysis, we calculated H → cc, bb, tt, τ τ , γγ, gg, hh and A → cc, bb, tt, τ τ , γγ, gg, assuming the MSSM Higgs potential for H → hh.
Although not pursued in this paper, our results could be weakened upon further consideration of the Higgs potential where BR(H → hh) is sizable. The branching fractions involving weak gauge bosons are vanishing in the alignment limit and we do not consider those decay modes. Note, however, that our final states nevertheless mimic certain di-boson-like signals and hence can be potentially constrained using the corresponding searches proposed in [28][29][30].
We consider the limit on σ × Br(Φ → τ τ ) using the latest ATLAS data [49] for gluon-gluon fusion (ggΦ) and b-annihilation (bbΦ). In our analysis, the production of the neutral Higgs bosons, H and A, via these processes are calculated using SuShi [51,52] and we further assume m H = m A to avoid bounds from custodial symmetry breaking. Values of the production cross sections at √ s = 13 TeV are shown in  For the constraints discussed later in our reference model, we use these values of production cross sections and consider a parameter set to be excluded if either is larger than the experimental bounds on σ × BR (Φ → τ τ ), where BR (Φ → τ τ ) is calculated in context of the model. In estimating the future sensitivity at the HL-LHC, we assume the statistical uncertainty dominates and simply rescale the upper bound on the cross section by √ R L , where R L := L run2 /L HL , to estimate the future sensitivity at the HL-LHC, where the integrated luminosities are L run2 = 139 fb −1 and L HL = 3 ab −1 .

2 + E miss
T channel All three leptonic cascade processes in Fig. 1 contribute to the 2 + E miss T channel. We recast bounds for slepton and chargino masses in SUSY models from Ref. [39]   whose signal region includes a lepton pair which does not reconstruct the mass of a weak gauge boson.
Let us first explain how we obtain the sensitivities. We calculate the 95% CL s limits by requiring [53] CL s (µ = 1) : where F is the cummulative distribution function of the normal distribution.
Here, the asymptotic formula for the distribution of test statistic q µ , defined as is used [54], where the likelihood function is defined as: (2.4) In the above definition, ν i = µs i + b i with s i , b 0 i , σ i being the number of signal events, the number of background events and the error of background in the i-th bin, respectively. Here, we assume that only background events have uncertainties and obey the Gaussian distribution. (μ,b) are the values of µ andb = {b i } which maximize the likelihood, andb = {b i } is the value of b i for a given µ. Note that n i corresponds to the observed data (central values of backgrounds) for q obs µ (q A µ ). For future limits, the significance for exclusion is given by [54] Z excl q n=b µ=1 , n i = b i , (2.5) while that for discovery is given by For the future sensitivity with 3 ab −1 data, the number of events (error) are rescaled by R L ( √ R L ). In Ref. The number of signal events in each m T 2 bin is calculated as where L is the integrated luminosity. Here, Φ, P and J run over scalar fields, production process and decay modes, respectively. The cut acceptance times the detector efficiency, J,P i , is a function of the scalar and new fermion masses and is the fraction of events that pass the cuts in a given signal region. From now on, we simply call the J,P i values the efficiencies. These are calculated by means of Monte Carlo simulations. In our analysis, we generated events using MadGraph5 2 8 2 [56] based on a UFO [57] model file generated with FeynRules 2 3 43 [58,59]. In order to boost up the speed of the simulation, decays of the neutral heavy Higgs boson are handled using MadSpin [60]. Showering and hadronization are controlled by PYTHIA8 [61], and detector simulation by Delphes3.4.2 [62]. We used the default ATLAS card for the fast detector simulation. Jets are reconstructed using the anti-k T algorithm [63,64] with ∆R = 0.4. For the decay modes J = EZ, EW and NW, which are depicted from the left-to-right panels respectively in Fig. 1, we simulate the processes where V (V = Z, W ) indicates leptonic decays including τ leptons. 3 Production of Φ is simulated via the CP-even Higgs boson production from gluon-gluon fusion in the 4 flavor scheme with the 5-dimensional effective interaction and the b-annihilation process in the 5 flavor scheme separately. In the simulation up to two additional partons are included and these are matched to the showered events by the MLM matching [65] with xqcut = m H /10. The m T 2 distributions from our leptonic cascade decays are shown in Fig. 3. We see that m T 2 has a relatively sharp edge for lighter e 4 masses in the EZ and EW decays where all the E miss T comes wholly from the decay products of e 4 , while the behavior is the other way around in the NW decay. Thus, the m T 2 cuts can also be efficient in a wide range of parameter space in discriminating the leptonic cascade decays from the SM backgrounds, in particular for W W and tt events whose m T 2 distributions are mostly shifted below ∼ 160 GeV in all the SRs (See Fig. 5 of Ref. [39].).
Among the m T 2 bins, the highest bins with m T 2 > 260 GeV are the most important to discriminate the signals from SM backgrounds. From the ggΦ (bbΦ) production followed by the EZ decay, the efficiencies to the m T 2 > 260 GeV bins are roughly 1-5% (5-10%) in SF-0J, and 5-10% (5-10%) in SF-1J, where nJ (n = 0, 1) is the number of jets in the SRs.

3/4 channel
Leptonic cascade decays involving the Z boson, i.e., the left-most panel of Fig. 1 (the EZ mode), can produce fully leptonic final states when including Z → , which yields clear multiple leptonic resonance signals. The process is hence constrained by the searches for 3/4 events depending on the cuts for the softest lepton. For this analysis, we recast the search results of Ref. [40].
Among the SRs defined in Ref. [40], the 4 SRs with one Z-pair and small E miss T have the strongest sensitivities, although there could be sub-dominant contributions from the other SRs due to the mis-reconstruction of leptons. There are two bins with m inv < 400 GeV and m inv > 400 GeV, where m inv is an invariant mass of the four leptons, which corresponds to the mass of the neutral Higgs Φ in our reference model. This is shown in the lower-right panel of Fig. 3, where the events are clearly clustered around m Φ = 1 TeV. Since all SRs in Ref. [40] are mutually exclusive, we combine them all in the same way as in the 2 + E miss T search. The efficiencies of the EZ decay in the SR requiring four leptons, with one Z-like lepton pair with E miss T < 50 GeV and m inv > 400 GeV are roughly 5-10% for both ggΦ and bbΦ productions.

Current constraints and future prospects
In the following subsections, we present the current limits based on the LHC run2 data with 139 fb −1 integrated luminosity and future prospects at the HL-LHC with luminosity 3 ab −1 based on the analysis strategies explained in the previous section. We emphasize that our analysis results can be generically applied to other BSM scenarios which share the same kinematic topology and final states; thereby modelindependent upper limits and prospects on the total cross sections will be shown first. For scenarios with a resonant particle production like a Higgs boson, i.e., via gluon-gluon fusion and/or b-annihilation, our results on the branching ratios can be readily applied. Finally, we provide the model-dependent constraints and expected sensitivities for the vectorlike lepton mass and heavy neutral Higgs mass in a 2HDM type-II.

Model independent limits on cross sections
and respectively 4 . Here, ( 4 , , V, ) = (e 4 , µ, Z, µ), (e 4 , µ, W, ν µ ), (ν 4 , ν µ , W, µ) refers to the EZ, EW and NW decays, respectively. Tables of the values used in these figures are attached in supplemental material. Note that the production and decays of Φ are handled separately by MadGraph and MadSpin respectively to boost up the speed of our event generations; thereby our method does not technically include the cases of an off-shell production of Φ or the breakdown of the narrow width approximation, which might provide different values of the cut acceptances. The number of signal events in the i-th m T 2 bin is calculated by 2 + E miss T , Run-2 2 + E miss T , Run-2  without summations for each choice of Φ, J, P , and we obtain upper limits on σ P Φ × Br Φ,J . The current limits are O (1−10) fb for m Φ 1 TeV for the EZ and EW decay modes. The search is more sensitive for larger m e 4 due to the increasing number of events passing the m T 2 cut for larger lepton masses, as displayed in Fig. 3. The limits on the NW mode is of the same order, but the search is more sensitive to the smaller m ν 4 since the m T 2 distribution shows the opposite behavior with the vectorlike lepton mass. In both cases the sensitivities increase with increasing m H as a larger Higgs mass tends to produce events with larger values of m T 2 . The bottom-right panels in Figs. 4 -7 show limits on the EZ mode, where Z → , which is constrained by the 3/4 search channel. The limit is almost independent with respect to the vectorlike lepton mass because m inv is determined by m Φ , so the only requirement is that four leptons should be reconstructed with a sufficient p T . At the HL-LHC, we expect the experimental sensitivities to the cross sections would be increased by about √ R L 4.7 and hence they can be improved to be O (0.2 − 1) fb. 2 + E miss T , Run-2

Constraints on the branching fractions
We now discuss the constraints on the branching fractions of the neutral heavy Higgs bosons. As stated earlier, our bounds on the branching fractions can be applied to a wide range of other scenarios with a resonant particle production like a Higgs boson. Figures 8 and 9 (10 and 11) show the upper bounds on the branching fractions from the current (future) data assuming a production cross section for Φ as in a 2HDM type-II for tan β = 10 and 50, respectively. The number of signal events is calculated in a combined way by where Φ = H, A and P = gg, bb, but there is no summation over the decay modes J = EZ, EW, NW. The production cross sections from gluon fusion and b-annihilation are calculated using SuShi [51,52]. In these figures, we neglect the difference in the branching fractions between the CP-even and CP-odd Higgs bosons, which is expected to be small for m 2 t /m 2 H 1, i.e., we set Br H = Br A . We assume that only one of the three decay modes contributes to the signal regions to extract the corresponding limit. Here, we also show the limits on the NW decay mode for completeness, although this branching fraction for a mostly singlet-like ν 4 is constrained 2 + E miss T , HL-LHC to be below ∼ 5 % (0.1 %) for m H 340 GeV (above 340 GeV) in our reference model. This decay is constrained mainly by the Φ → τ τ search results but also from electroweak precision measurements (see the discussion in Appendix A.). The limits obtained by combining all decay modes are presented in the next section.
The 2 +E miss T search excludes the Higgs boson mass up to about 1.3 (2.1) TeV for m e 4 500 GeV and tan β = 10 (50) if the branching fraction is 50%. The sensitivity of this search becomes weaker for smaller m e 4 , because the m T 2 distribution drops off around m e 4 as shown in the top panels of Fig. 3, and hence cannot pass the m T 2 cut. The limits for the EZ decay mode are tighter than those for the EW decay mode because of the broader m T 2 distribution. For the NW decay, although the branching fraction cannot be sizable in our reference model, a similar range of Higgs masses can be probed unless the mass difference is small, implying a smaller m T 2 . At the HL-LHC with 3 ab −1 data, the limits will be strengthened to about 1.8 (2.8) TeV for tan β = 10 (50) and 50% branching fraction.
The limits on the branching fraction of the EZ decay mode from the 3/4 search are shown in the lower-right panels of Figs. 8, 10, 9 and 11. The current (future) sensitivities on m H with the 50% total branching fraction are about 1.2 and 1. of the masses of the vectorlike leptons because the key kinematic cut is given by m inv = m Φ .

Model dependent constraints
Finally, we apply our results to the reference model: a 2HDM type-II augmented by vectorlike leptons. The details of our reference model are described in Refs. [28,30] but we simply show the Lagrangian and explain the field contents briefly to aid with our discussion. The most general Lagrangian of Yukawa interactions and mass terms relevant for our processes include:     leptons including a SM singlet lepton N L,R (denoted by κs). The final line denotes the vectorlike mass terms of the new leptons. In order to avoid strong bounds in the Higgs sector, we take the alignment limit [66][67][68][69], i.e., α = β − π/2 where α is the neutral Higgs mixing angle and tan β = H u / H d is the ratio of the vacuum expectation values of the two Higgs doublets. More details on the model, the field contents, mass mixing, and the interactions among the mass eigenstates are explained in Ref. [28]. In our numerical analysis, we scan the parameters to vary in the ranges: where we consider c max = 1 and 3.5 ∼ √ 4π, with the latter being motivated by the upper limit of couplings near the weak scale from perturbativity. Since there are many parameters in the Lagrangian above, we consider an optimized parameter scan strategy; two representative cases maximizing BR(H → e 4 µ) are picked to focus on emphasizing the current and future sensitivities of our processes: The "light-L" denotes the case where the lightest new leptons e 4 , ν 4 are almost isodoublet, i.e., (e − 4 , ν 4 ) ∼ (L − , L 0 ). The "light-E" denotes the case where e 4 is almost isosinglet, i.e., e 4 ∼ E. In order to focus on our leptonic cascade processes, we require the other vectorlike leptons to be heavier than the neutral Higgses, H/A. Note that these two representative cases correspond to the simple scenarios where the branching ratios of e 4 typically follow the pattern expected by the Goldstone boson equivalence theorem, i.e., BR(e 4 → W µ):BR(e 4 → Zµ):BR(e 4 → hµ) = 2:1:1 (for isosinglet) and 0:1:1 (for isodoublet), and the approximation (λ L , λ E , λ,λ)v/(m L , m E ) 1 is valid in most of the parameter space; the couplings among the mass eigenstates are easily expressed with the Lagrangian parameters as in Ref. [30,70]. As pointed out in Ref. [34] for the vectorlike quarks, general vectorlike lepton scenarios can include the possibilities of small λ L /λ E as well as sizable mixing between the isodoublet and isosinglet, which allow all the values between 0 and 1 for the e 4 branching ratios.
In the parameter scan, we do not study the case in which ν 4 ∼ N when Φ → ν 4 ν µ → W µν µ can dominate for κ N > 0.5. This is because, in the limit where (λ,λ, λ E , λ L , κ,κ, κ N )v/(m N , m L ) 1 is valid, large values of κ N maximizing the decay width Φ → ν 4 ν µ are strongly constrained by the electroweak precision mea-surements, especially from the Fermi constant G F at large tan β. 5 Moreover, the values of BR(Φ → ν 4 ν µ ) are limited by the competing decay modes, Φ → τ + τ − and Φ → tt, and hence we find the total values of BR(Φ → ν 4 ν µ → W µν µ ) are preferred to be smaller than 4%. See Appendix A for more details.
In the light-L (light-E) case, we scan the absolute value of the Yukawa coupling constant λ L (λ E ) in the [0.5, c max ] range, so that the branching fraction BR(Φ → e 4 µ) covers the full range of possible values. For each point, we calculate the contributions to the electroweak precision observables (EWPOs), cross sections to the di-tau channels, CL s , Z excl , and Z disc .
Note that the heavy Higgs Φ can decay into the vectorlike lepton pair, e.g., Φ → e 4 e 4 , but the decay width is suppressed by (λ L , λ E )v u /(m L , m E ) compared to that of Φ → e 4 µ in our simple scenarios "light-L" and "light-E". Nevertheless, in general scenarios, such a double e 4 production can be sizable and provide another interesting signature, which is beyond the scope of this paper. Recall also that the vectorlike leptons can be pair produced through gauge boson interactions which can be subjected to robust bounds, as in Ref. [41]. The most recent analysis by ATLAS in Ref. [73] shows that a nominal bound of m e 4 ,ν 4 800 GeV can be obtained assuming our total cross section of pp → e 4 ν 4 → W W ν µ µ is similar to what is expected in the type-III seesaw model. We do not include this bound (after rescaling the production cross section) in our parameter scan but we emphasize that some light e 4 region can be constrained further by this complementary search.
The upper (lower) panels of Fig. 12 show the sensitivities of our reference model in the 2 + E miss T analysis for the light-E (light-L) case when c max = 1. In the left panels, where the sensitivities are represented in the m H − m e 4 plane, the pink (2 run2), green (2 disc.) and blue (2 excl.) points correspond to the regions which can be covered by our 2 + E miss T analysis with current data at 95% C.L., the corresponding discovery region at the HL-LHC, and the future exclusion region, respectively. Note that here exclusion regions mean that some choices of parameters would be excluded. We include the contributions from the three decay modes, EZ, EW and NW, as defined in Eq. (2.7).
The points are plotted in the order given in the legend, i.e. pink points are plotted over green points, and green points over blue ones. In the right panels, we classify the blue points (2 excl.) by whether they are expected to have sensitivity to the Φ → τ τ search or not at the HL-LHC. The blue (2 only) and cyan (2 & τ τ ) points are both within the future sensitivity of the 2 + E miss T channel but we expect the searches for Φ → τ τ to also have sensitivity in the cyan region. The brown (τ τ run2) points are excluded by the current τ τ search, and the purple (τ τ excl.) Figure 12. Sensitivities of the 2 + E miss T search to the light-E (light-L) scenario with c max = 1 in the upper (lower) panels. In the labels, "2 " is a short for 2 + E miss T search and "τ τ " is the di-tau search. In addition, "run2", "disc." and "excl." indicate that the points are constrained by the current run-2 data, expected to be discovered and constrained by the future HL-LHC data, respectively. See the texts for more details.   points correspond to the region we expect to be possibly excluded by the future τ τ search at the HL-LHC but not by the future 2 + E miss T search. The white points (no excl.) are not excluded by any of the search channels discussed in this paper. The brown solid (purple dashed) line corresponds to the nominal exclusion limit from the current (future) sensitivity of Φ → τ τ searches without the presence of vectorlike leptons. These nominal bounds can be pushed back to the brown or purple scattered dots in the presence of vectorlike leptons depending on the parameter choices. The (small) purple region in the top-right corner denotes the region within the reach of the τ τ search but not of our leptonic cascade process.
We find the current search for the 2 + E miss T channel is sensitive to m Φ ≤ 1.7 TeV and m e 4 ≤ 1 TeV when c max = 1 and the sensitivity to Higgs masses close to 1 TeV remains for any tan β > 1. This is quite promising even compared to the conventional heavy Higgs search Φ → τ τ . Moreover the future sensitivities with 3 ab −1 extend to m Φ 2.2 TeV and m e 4 1.5 TeV covering a much wider range of tan β than the conventional searches. Figure 13 shows the corresponding sensitivities for the 3/4 search. The labels 4 indicate that these are limits or sensitivities from the 3/4 search instead of the 2 + E miss T search. In the light-E case, the limit is much weaker than that from the 2 + E miss T search because the EW decay mode does not contribute to the SRs of the 3/4 search. While in the light-L case, the limit is similar because BR (e 4 → Zµ) is larger than in the light-E case. We may be able to test whether the lightest charged vectorlike lepton e 4 is almost singlet-like or doublet-like by using both of 2 + E miss T and 3/4 channels in a complementary way. For example, if we were able to discover the signal at m Φ = 1.7 TeV and m e 4 = 1 TeV for c max = 1 in both channels, the discovered e 4 should be almost doublet-like.
In order to demonstrate the maximal capability of our analysis, we also allow larger Yukawa coupling constants up to c max = 3.5. The corresponding results are shown in Figs. 14 and 15. The branching fractions can be larger than those in the c max = 1.0 case, and thus the searches can probe heavier mass regions: m Φ 2 TeV and m e 4 1.5 TeV with current data, and m Φ 2.5 TeV and m e 4 1.8 TeV at the HL-LHC. Note that the presence of vectorlike leptons with large couplings can significantly dilute the typical search for H → τ τ as seen by the absence of cyan points and purple triangles in Fig. 14 (where they are completely covered by the pink and green points) as compared to Fig. 12. We found that the τ τ search at the HL-LHC can lose its sensitivity up to about the run2 limit, i.e. the brown solid line, due to the large Yukawa coupling constant λ L or λ E ∼ 3.5. In that case, the leptonic cascade search strategy presented here is necessary to probe the details of the Higgs sector.

Conclusion
In this paper, we have studied the potential of leptonic cascade decays of a heavy neutral Higgs boson through vectorlike leptons as a simultaneous probe of an extended Higgs sector and extra matter particles at the LHC. The fully leptonic final states contribute to multi-lepton signals such as 2 + E miss T and 3/4 which are already under investigation motivated by BSM scenarios such as SUSY and seesaw models. We have found that by simply recasting the existing searches we can obtain new strong constraints to any BSM scenarios sharing the same event topology and final states as our reference scenario, depicted in Fig. 1. Therefore, some of our analysis results are shown in model independent ways. The 95% C.L. model independent upper limits on the total cross sections are found to be O(1 − 10 fb) for a heavy neutral Higgs (or a resonant particle producing the same topology) in the mass range between 1 -3 TeV, using the current run2 data with 139 fb −1 . Naively rescaling the size of the data to an integrated luminosity of 3 ab −1 , assuming the statistical uncertainty dominates, the future sensitivities at the HL-LHC extend to O(0.2 − 1 fb).
The model-independent bounds could be transformed into a useful form, e.g. the total branching fraction of the resonant particle in Fig. 1, in a wide class of BSM scenarios where the resonant particle production cross sections are the same as (or simply rescaled from) the neutral heavy Higgs production cross sections in 2HDM type-II. If the branching fraction of one's interest is 50%, the heavy Higgs mass up to about 1. Although not implemented here, we expect that further investigation of an additional lepton resonance from the decay of e 4 would increase the sensitivity of the 3/4 channel. Therefore, we conclude that searches for our leptonic cascade processes can shed light on probing both an extended Higgs sector and extra matter, and are generally more advantageous than conventional heavy Higgs searches in these scenarios.

A Approximated expressions in the reference model
In this appendix, we summarize analytical formulae of the important EWPOs and decay widths at leading order in: Here the Lagrangian parameters are given in Eq. (3.5). Further assuming enough splittings between m L and m E (m N ), the lightest charged (neutral) vectorlike lepton e 4 (ν 4 ) is mostly isodoublet or isosinglet (SM singlet). Detailed approximate expressions of gauge and Yukawa couplings in the mass eigenstate basis are in Refs. [30,70]. Here, we focus on the formulae useful in applying constraints from electroweak precision measurements and several important branching fractions.

A.1 Electroweak precision measurements
The vectorlike leptons in our reference model are constrained by various electroweak precision measurements listed in Table 1, which are also discussed in Refs. [30,70]. We calculate these observables at tree-level for the first 6 observables and at the one-loop level for the last 4 observables.
The contributions by the vectorlike leptons to the oblique corrections are almost the same as the formulae in Ref. [74] given for the vectorlike quarks except the definitions of the following functions: due to the difference in hypercharges. Here, y a := m 2 νa /m 2 Z , y β := m 2 e β /m 2 Z (a, β = 1, 2, · · · , 5) and the function f (y 1 , y 2 ) is defined in Ref. [74]. 6 The ratio R γγ := Γ (h → γγ) /Γ (h → γγ) SM is calculated with the formula shown in Ref. [76].
In our reference model, the contributions of vectorlike leptons to the Fermi constant G F , determined from the muon decay, and R µ := Γ (Z → had)/Γ (Z → µµ) are the most strongly constrained observables in Table 1. Hence, we provide the approximate expressions (when Eq. (A.1) is valid) of the leading contributions to G F and R µ : where t β = tan β. For t β 1, the leading contribution to G F is κ 2 N /m 2 N and the other contributions are suppressed by t 2 β . From Eq. (A.4), the upper bound on κ N is estimated as Therefore, κ N should be less than O (10 −3 ) to be consistent with EWPOs. The limits on λ L and λ E are weaker by a factor tan 2 β, and hence λ L,E ∼ O (1) is allowed for sufficiently large tan β and vectorlike lepton masses.

A.2 Branching fractions
The partial widths for the heavy netural Higgs decays to the charged leptons, e + a e − b (e a = µ, e 4 , e 5 ), are given by where β(x, y) := 1 − 2(y + x) + (y − x) 2 . Here, Y H e is the heavy CP-even Higgs Yukawa matrix of the charged leptons in the mass basis. Those for the CP-odd Higgs boson can be obtained by replacing H → A, and those for the neutrinos are given by replacing e → ν and µ → ν µ .
When e 4 and ν 4 are mostly doublet-like, i.e. (e − 4 , ν 4 ) ∼ (L − , L 0 ), the decay widths to a vectorlike lepton and a SM lepton are approximately given by and for mostly singlet-like leptons, i.e. e − 4 ∼ E, ν 4 ∼ N , Here the sub-dominant contributions suppressed by the new Yukawa couplings times v 2 /m 2 VLL , VLL = E, L, N , and the SM lepton masses are neglected. Note that H → ν µ L 0 appears only at sub-leading order and vanishes for κ N = 0. Thus, the Higgs boson decays mostly to a charged vectorlike lepton and a SM lepton when λ L (λ E ) is large and e − 4 L − (E). In the numerical analysis, the heavy Higges decay to SM particles, H → cc, bb, tt, τ τ , γγ, gg, hh and A → cc, bb, tt, τ τ , γγ, gg are calculated using the formulas presented in Ref. [1]. We assume the MSSM for the triple Higgs coupling. The decays to gauge bosons, H → ZZ, W W , and A → hZ vanish in the alignment limit.
The partial decay widths of the charged vectorlike lepton e 4 is given by where Y h e , g Z e L,R and g W L,R are the coupling matrices to the SM bosons in the mass basis, see Refs. [30,70] for the details. The mass squared ratios are defined as In our analysis, we assume only one of the vectorlike leptons is lighter than the heavy Higgs for simplicity. For the "light-E" scenario, i.e., e 4 ∼ E, the partial widths are approximately given by For the "light-L" scenario, i.e., e − 4 ∼ L − , the partial decay widths are approximately given by If ν 4 is mostly isodoublet-like, then we have while if it is mostly singlet-like, then (A. 26) In the limit of κ N λ L , λ E , the lepton ν 4 can dominantly decay to a W boson.  Figure 16 shows the values of the branching fractions of the leptonic cascade decay H → e 4 µ → V 2 µ, with V 2 = Zµ, W ν µ . In the upper (lower) panels, the hierarchy m E < m H < m L , m N (m L < m H < m E , m N ) is assumed such that e 4 E (L). The Higgs mass is 1 TeV and the heavier vectorlike lepton masses are set to be 3 TeV. The Yukawa coupling constant λ E (λ L ) is set to 1.0 or 3.5 ∼ √ 4π on the left and right panels, respectively, for the case of e 4 E (L). The other Yukawa coupling constants are set to 10 −3 for simplicity. In addition, we fix κ N = 0 to suppress the contribution to G F . The regions below the cyan (gray) line is 2σ away from the measured value of G F (R µ ). G F provides the stronger constraint for an isosinglet-like e 4 , while only that of R µ provides a constraint for an isodoublet-like e 4 due to the assumption κ N = 0. When λ E (λ L ) = 1, the branching fraction can be as large as 50% (40%) for a singlet-like (doublet-like) e 4 . If we allow larger coupling constants, BR (H → e 4 µ) ∼ 1 is possible, thus, the total branching fraction can be as large as BR (e 4 → V ), i.e., 50% (75%) for the doublet-like (singlet-like) e 4 . The decay to the neutral component L 0 vanishes because of κ N = 0. Figure 17 shows the branching fraction H → ν 4 ν µ → W µν µ where ν 4 is the SM singlet-like. The constraints from the electroweak precision measurements and the H → τ τ searches are also displayed with the brown and yellow lines, respectively. The value of G F is more than 2σ away from the measured value to the right of the red lines for m N = 250 (2500) GeV in the left (right) panel. In the left panel, the region above the yellow dash-dot line is excluded by the H → τ τ search. The red dashed line in the right panel shows the G F constraint for m N = 500 GeV. In the left panel, we set m H = 300 GeV < 2m t , where the branching fraction can be as large as 10 %. However, a wide range of parameter space tan β 0.75 is already excluded by the recent H → τ τ search. Below the yellow line, the CP-even (CP-odd) Higgs decay is dominated by H → hh (A → gg) because we assume a MSSM-like Higgs potential. Hence the decay to tau leptons are suppressed. This region can be further excluded by the searches for H → hh depending on the quartic couplings of Higgs potential [77,78]. The value of BR (H → ν 4 ν µ ) is even smaller for larger tan β because the decay widths to bottom quarks and tau leptons are enhanced by tan 4 β compared with the decay to the decay mode into ν 4 . On the right panel, we assume m H = 3 TeV which is well above the limit from the di-tau decay channel. Even if m N = 2.5 TeV, κ N 0.5 is required to be consistent with EW precision measurements. The Higgs width is dominated by decays to SM fermions, top quarks for small tan β while bottom quarks for large tan β, and BR (H → ν 4 ν µ ) is at most O (10 −3 ). The values of the branching fraction would be similar (or slightly smaller) in the regions m H < 3 TeV where our analysis has sensitivities and hence we conclude the leptonic process H → ν 4 ν µ → W µν µ would not be constrained in most of our parameter space.