Resonant heavy Higgs searches at the HL-LHC

We show the importance of searches for heavy resonant scalars. Taking cue from the present searches, we make projections for searches in an extended scalar sector at the high luminosity $($HL$)$ run of the Large Hadron Collider $($LHC$)$. We study the three most relevant search channels, $viz.$, $H \rightarrow hh, H/A \rightarrow t\bar{t}$ and $b\bar{b}H/A$, where $H$ $(A)$ is a heavy scalar $($pseudoscalar$)$. Upon studying multifarious final states for the resonant double Higgs production, we find that the $b\bar{b}\gamma\gamma$ and $b\bar{b}b\bar{b}$ channels impose the strongest limits. One can probe $\sigma(pp \rightarrow H \rightarrow hh)$ between $(104,17)$ fb for $m_H$ lying between $(300,600)$ GeV and between $(8,4)$ fb for $m_H$ lying between $(800,1000)$ GeV in the $b\bar{b}\gamma\gamma$ and $b\bar{b}b\bar{b}$ channels, respectively, at $95\%$ CL. For the $b\bar{b}H$ channel, we can exclude $\sigma(pp \rightarrow b\bar{b}H)$ between $(15,3)$ fb for $m_H$ lying between $(300,500)$ GeV. Our analyses will hold for several scenarios with heavy Higgs extensions. However, we consider the phenomenological Minimal Supersymmetric Standard Model $($pMSSM$)$ as an example and impose the present constraints from the SM$-$like Higgs measurements and our future direct search limits and obtain strong constraints on the $m_A-\tan{\beta}$ parameter space, where $m_A$ and $\tan{\beta}$ are respectively the mass of the pseudoscalar and the ratio of the vacuum expectation values of the two Higgs doublets. Assuming that the heavy Higgs boson decays only to Standard Model $($SM$)$ states, we find that the $H \rightarrow hh \rightarrow b\bar{b}\gamma\gamma$ $(H \rightarrow t\bar{t})$ channel excludes $\tan{\beta}$ as low as $4$ ($m_A$ between $(400,800)$ GeV) at $95\%$ CL. This weakens up to $6.5$ when the $b\bar{b}H$ dominates. However, upon allowing for non$-$SM decay modes, the limits change.


Introduction
The Higgs boson discovered in 2012, was the last missing piece in the Standard Model of particle physics (SM). The SM, however, is inadequate to explain the nature and existence of dark matter, the small but non-negligible masses of the neutrinos, the asymmetry between baryons and anti-baryons, to name a few. Besides, SM can not explain the hierarchy problem which is inherent in the theory. Well motivated theories including supersymmetry have the potential to solve some of these limitations. There are additional fundamental theoretical requirements that the SM can not satisfy. The aforementioned experimental observations and theoretical requirements compel us to look for physics beyond the Standard Model (BSM). However, the possibilities being innumerable, it is extremely difficult to ascertain the nature of such new physics. Since the discovery of the Higgs boson, and a growing convergence of its properties with the SM expectations [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16], the new physics possibilities are gradually getting strongly constrained. Searches for BSM are being performed at the Large Hadron Collider (LHC) by the CMS, ATLAS, ALICE and LHCb collaborations. Except for some excitement in the flavour physics sector, there have not been any strong hints for new physics in the form of new particles or significant deviations in couplings with respect to the SM. Even though supersymmetry is perhaps one of the most elegant theories of our time, it comes with additional new particles, which need to be discovered at some stage. Even though searches are being conducted for a considerable region of parameter space for the Minimal Supersymmetric Standard Model (MSSM), there are more non-traditional searches which need to be performed. The MSSM parameter space has been extensively studied in light of the constraints from cosmology, flavour physics, and Run-I plus Run-II data from LHC . Implications on the MSSM parameter space from ∼ 36 fb −1 dataset of Run-II has also been studied [39][40][41]. However, there are simple extensions of the MSSM that can weaken the present bounds considerably. On the positive side, the LHC can potentially pin down the Higgs couplings to weak bosons and most of the fermions fermions at the level of O(5 − 10%) [42,43]. However, as has been shown in numerous experimental [44][45][46][47][48][49][50][51][52][53][54][55][56][57] (including future extrapolations [58][59][60]) and phenomenological studies , the measurement of the elusive triple Higgs coupling (λ hhh ) is an almost impossible feat at the LHC. Studies show that future colliders are expected to be more adept in constraining or even measuring this coupling to a great precision [84,[103][104][105][106][107][108][109].
In the following sections, we focus on the various production and decay processes of a resonant scalar, viz., a resonant decay to a pair of SM-like Higgs bosons, to a pair of top quarks, a heavy pseudoscalar, A, decaying to an SM-like Higgs boson and a Zboson and the production of a heavy scalar in association with a pair of bottom quarks. The final theme of this work is in the context of the phenomenological MSSM (pMSSM). However, our results are presented in such a way that they can be mapped onto most models with an extended scalar sector. Table 1 summarises the various bounds set on the double-Higgs production cross-section by CMS and ATLAS in the non-resonant and resonant categories. The resonant production results are mostly interpreted in terms of spin-0 and spin-2 hypotheses. Besides, there are many supersymmetric interpretations for the resonant scalar searches. As an example, for the bbτ + τ − resonant search performed by CMS [110], the MSSM parameters m A (mass of the CP -odd scalar, A) and tan β (ratio of the vacuum expectation values of the two Higgs doublets in the model, viz., H u and H d ) are excluded in the range 230 GeV < m A < 360 GeV and tan β 2, at 95% CL.
In order to be completely sure whether or not there is any extended Higgs sector, it is of utmost importance to measure the Higgs quartic coupling, λ hhhh and the Higgs trilinear coupling, λ hhh = λ hhhh v, where v is the vacuum expectation value of the SM Higgs boson. Now, independent measurements of the Higgs couplings to the gauge bosons and fermions will constrain v and we already have a precise Higgs mass measurement. To confirm this sector of the SM, one needs to measure λ hhh or λ hhhh . The only way to probe the Higgs quartic self-coupling directly is via triple Higgs production [111][112][113][114] and it is an impossible feat to achieve at the LHC because of a signal rate of the order O(0.1) fb. Upon using the above relation, one can indirectly measure λ hhh from Higgs pair production [115][116][117]. Models with extended scalar sectors can modify the Higgs self-coupling [69,79,80,118,119] and thus it is very important to precisely measure the Higgs self-couplings and reconstruct the Higgs potential. Any significant deviations in this sector will give us hints into the kind of new physics that we are looking for. Besides measuring deviations to the Higgs self-coupling, there are other possible channels to look for in order to establish an extended scalar sector. Some of these new channels include the production of the SM-like Higgs in association with a Z-boson reconstructing a resonance. Another possible channel is the production of a pair of top-quarks. Now, the first of these channels can be via a heavy pseudoscalar resonance [120,121], whereas the tt production can be either through a heavy scalar or pseudoscalar [122,123]. However, both these channels can also come about from a heavy Z [121,[124][125][126][127]. Now, in order to be sure whether the Zh or tt production is via a spin-0 or spin-1 resonance, one needs to delve deeper into the angular observables. Lastly, we also study the effects of the high tan β regime for a heavy scalar produced in association with a pair of b-quarks and decaying to a pair of τ -leptons [128].
Our paper is organised as follows. We study the reach of the HL-LHC for the H → hh channel in various final states, in section 2, by showing the 95% CL bounds on σ(pp → H → hh) as functions of the heavy Higgs mass, m H . Following the Higgs pair production,  [136,137] (300-900) (500-3000) W W * W W * 160 9300-2800 [138] (260-500) Table 1: Bounds obtained on the di-Higgs cross-sections (in fb) from CMS and ATLAS studies dedicated to the search for non-resonant (NR) and resonant (R) double Higgs production in various channels. The numbers in brackets show the range of the heavy scalar mass considered in that particular study.
we address the couplings of the CP -even heavy Higgs to a pair of top quarks and to a pair of bottom quarks in sections 3 and 4 respectively. In section 5, we briefly use the previous results to recast our limits in the purview of the phenomenological MSSM (pMSSM) and show the future reach of these searches in the m A − tan β parameter space. We finally summarise our results and present our future outlook in section 6.

The pp → H → hh Channel
As discussed in the introduction, the objective of this work is to scrutinise the viable scalar extensions of the Standard Model (SM). In this section, we focus on a heavy CP -even scalar produced via the gluon fusion mode and subsequently decaying to a pair of SMlike Higgs bosons. In the next five subsections, we study the reach of the HL-LHC in constraining the resonant Higgs pair production cross-section, σ(pp → H → hh), upon studying multifarious channels, viz., bbγγ, bbbb, bb, τ + τ − , bbW + W − and γγW + W − . Many of these channels with τ s and W give different signatures upon considering leptonic or hadronic modes. We study all possible final states giving importance to the total rate as well the cleanliness.

The bbγγ Channel
The bbγγ final state is the golden channel when it comes to studying the non-resonant double Higgs production. The cleanliness of this channel, owing to smaller backgrounds, triumphs over the reduced rate (Br(h → γγ) ∼ 0.2%). Here however, we turn to a resonant scalar production which decays to a pair of SM-like Higgs bosons. Our goal is to ascertain the reach of the HL-LHC in measuring σ(pp → H → hh) for a range of heavy Higgs masses (m H ). One of the reasons for this final state being a favourite is that the reconstruction and identification precision of photons at the LHC is very high.
Even though the signal seems to have a clean final state, there are several backgrounds at play which need to be dealt with carefully. The major backgrounds typically have the form of hh + X which includes the SM double Higgs production, h + X which includes Zh, hbb and tth, and null-Higgs processes like tt + ttγ where jets or leptons may fake photons, bbγγ + ccγγ + jjγγ (henceforth termed as bbγγ * ) where for the latter two, the light-jets may fake b-jets. Other fake backgrounds include bbjγ + ccjγ (we will refer to it as Fake 1 category), bbjj (referred to as the Fake 2 category), where the c-jets may pose as b-jets and the light jet may mimic a photon, and the single Higgs processes, viz., hjj + hcc (classified henceforth as the hjj * category), where the light-jets and c-jets may mimic b-jets. One of the major difference between most of the backgrounds and the signal lies in the invariant mass distribution of the b-jets, m bb . However, even when the m bb distribution of the signal (as well as the non-resonant SM di-Higgs background) peaks around the SM-like Higgs mass, m h , it is broad and can have considerable overlap with the m bb distribution either ensuing from a Z-boson or from a continuum. It should be noted that the most dominant backgrounds come from the QCD-QED bbγγ * , tth and SM-like di-Higgs processes. The former being a continuum covers a large part of the kinematic variable space with the signal. The SM-like di-Higgs background also has very large overlap with the signal. One of the easiest way to break this degeneracy is to utilise the m bbγγ or reconstructed m hh distribution which has a clear peak around the heavy scalar mass, m H , for the signal.
We generate the QCD-QED bbγγ and Zγγ → bbγγ backgrounds upon merging with one additional jet. We employ the MLM merging scheme [139] where the extra jet contains gluon, light quarks, c-as well as b-quarks. Among the h + X category, the Zh is generated with Higgs boson decaying to a pair of photons and the Z-boson decaying to a pair of bottom quarks. Furthermore, the tth and bbh backgrounds are generated with h → γγ. The major fake backgrounds with jets in the final state are generated with the aforementioned jet definition with one exception. We define both of the jets in the jjγγ channel in a way as to have no overlap with the bbγγ background. In case of the tt+X backgrounds, we generate the tt events with both of the top quarks decaying leptonically which ultimately fakes as photons. However, for the ttγ background, we require one of the tops to decay leptonically and the other hadronically. Finally, we generate separate single Higgs backgrounds via gluon fusion in association with c-quarks and also with light jets. The separation between the hcc and hjj backgrounds are necessary in order to appropriately take into account the different fake rates for c → b and j → b. All of these backgrounds are generated with specific cuts at the generation level which we summarise in Appendix A.
The idea of this section is to understand the reach of the HL-LHC in constraining models with extended scalar sectors. We thus employ optimised search strategies for a varied range of scalar masses. We vary m H in the mass range 275 GeV and 1 TeV. Specifically, we consider the following benchmark points, viz., m H = 275, 300, 350, 400, 450, 500, 550, 600, 800 GeV and 1 TeV. In line with our previous work [100], we first perform a classical cut and count analysis to gauge the sensitivity of various benchmark points. We closely follow various cuts from the ATLAS projection study [59]. Namely, we require exactly two btagged jets and two photons in the final state. The photons are required to have transverse momenta, p T > 10 GeV and a pseudorapidity coverage of |η| < 2.5. Moreover, the two photons are also required to lie within the pseudorapidity range, |η γ | < 1.37 (barrel region) or 1.52 < |η γ | < 2.37 (endcap region). After imposing these basic requirements, we apply some stronger selection cuts in order to enhance the signal to background ratio, S/B. We require the invariant mass of the pair of photons, m γγ , to reconstruct sharply about the SM-like Higgs mass in the range (122,128) GeV. We also impose lower bounds on the transverse momenta of the leading and sub-leading b-jets and photons. Moreover, upon inspecting the ∆R γγ and ∆R bb , we find that with larger values of m H , the SM-like Higgs bosons are more boosted yielding more collimated final states. We thus require ∆R γγ and ∆R bb to lie in the range (0.4,3.0), (0.4,2.0) and (0.4,1.5) for m H = 275 and 300 GeV, m H = 400, 450, 500, 550 and 600 GeV, and m H = 800 GeV and 1 TeV respectively. Besides, we require the invariant mass of the b-jets, m bb to lie in the range (90,130) GeV. This choice is related to account for the jet-energy correction and has been described in Ref. [100]. We summarise these cuts in Table 2. As the next logical step, we delve deeper into the kinematics. We reconstruct the invariant mass of the m bbγγ system and the transverse momentum of the SM-like Higgs decaying to a pair of photons, p T,γγ . As can be seen from the p T,γγ distribution in Fig. 1, the spectrum is harder for heavier values of m H . We choose p T,γγ > 50 GeV for m H = 275 and 300 GeV and for all other masses, we choose the transverse momentum of this reconstructed Higgs to be larger than 100 GeV. Thus, after these fixed cuts, we perform a simplified optimisation with the m bbγγ variable in order to enhance the signal over background. These cuts are finally tabulated in Table 3 where we also present the signal efficiency and the background yield at an integrated luminosity, L = 3 ab −1 with m H being varied. The signal efficiency, , here points to the ratio of the total number of signal events remaining after all the cuts applied in sequence to the generated number of events. The second column refers to the range of m bbγγ that optimises the signal and the third column denotes the minimum p T for the diphoton system. Finally, we provide a detailed cut-flow table for m H = 400 GeV in Table 4 with a complete information Fixed cuts 122 GeV < m γγ < 128 GeV p T,b > 40 (30) GeV, p T,γ > 30 (30) GeV 0.4 < ∆R γγ < (3.0/2.0/1.5), 0.4 < ∆R bb < (3.0/2.0/1.5), ∆R γb > 0. 4 90 GeV < m bb < 130 GeV   Upon utilising these results, we are able to set an upper limit on the production cross section of the heavy Higgs in a model independent manner 6 . We calculate the cross-section reach by using the significance formula viz. S/ √ B > N , where N denotes the number of confidence intervals. Here, the signal yield, S, is defined as σ(pp → H → hh → bbγγ)×L× and B represents the total background yield after the cut-based analysis. With this, we derive σ(pp → H → hh) at the N σ level, with N = 2 and 5 respectively corresponding to a 95% and 99.7% confidence level (CL) upper limit. We show the final results in Fig. 2 with the upper limit on σ(pp → H → hh) as a function of m H and we display the 2σ and 5σ lines. The 2σ upper limit is strong between 400 GeV and 1 TeV, varying between ∼37 fb and ∼7 fb. We do not show the results for a multivariate analysis as the sensitivity does not improve significantly.

The pp → A → Zh Channel
With the accumulation of more data, we are on the brink of accepting the fact that in MSSM or in generic two Higgs doublet models, the SM-like Higgs is in the decoupling regime with its coupling to the SM gauge bosons being proportional to sin(β − α), with α and tan β being the mixing angle in the neutral CP -even sector and the ratio of the two vacuum expectation values of the two doublets, respectively. For the up and down type fermions, the Yukawa couplings for the SM-like Higgs boson are proportional to cos α/ sin β and sin α/ cos β respectively. In the decoupling regime sin(β − α) ∼ 1 and hence the decay width A → Zh which is proportional to the coupling cos(β − α), is small. In a non-2 bbγγ + ccγγ + jjγγ. 3 bbjγ + ccjγ. 4 bbjj. 5 (gg → hjj) + (gg → hcc) 6 We consider the cut-based optimisations as final as we did not obtain any observable improvement with a multivariate analysis. decoupling regime, pp → A → Zh can give us deep insight into two scalars simultaneously. In this small subsection, we will focus on the aforementioned channel in the bbγγ final state. We will remain agnostic to the fact that the prospects of observing A → Zh in the decoupling regime may be extremely small. We follow a similar search strategy as before. The main difference here is the fact that in the previous analysis both the diphoton and the bb pairs are required to peak around the SM-like Higgs boson mass, whereas in the present case, the b-jets are required to peak around the Z-boson mass. We follow a cutbased analysis as before and optimise the m bbγγ and p T,γγ cuts for different values of m A . These variables are shown to have substantial discriminatory power and are shown in Fig. 3 for m A = 220 GeV and 400 GeV. Details of the optimised cuts are presented in Table 5. In Table 6 we show the cut-flow table for m A = 400 GeV. After a full optimisation, we show the 95% and 99.7% CL exclusion for σ(pp → A → Zh) in Fig. 4. The bounds are weaker than their H → hh counterpart mainly because of a larger overlap with the Zh background.

The bbbb Channel
After having studied the cleanest possible di-Higgs final state, we turn our attention to the one with the largest rate, viz., pp → H → hh → 4b. Several searches have already been conducted in this channel [56,130,131,143] and provide some of the strongest bounds 7 tt + ttγ 8 bbγγ + ccγγ + jjγγ. 9 bbjγ + ccjγ.     The dominant backgrounds include the multijet production from QCD processes and the top pair production. For the multijet backgrounds, we dissect our sample generation into three different categories each having at least two b-quarks, viz., bbbb, bbcc and bbjj, in order to have sufficient statistics to take into account the different tagging efficiencies and fake rates. We do not generate the h + jets and Z + jets backgrounds separately but we include their tree-level diagrams while generating the multijet backgrounds as they have negligible contributions [56]. We must mention here that we do not consider other possible sources of multijet production viz., cccc, ccjj etc. as these processes will be highly suppressed (with respect to bbbb) upon multiplying by the fake efficiency factors, in succession. We generate the tt background with the top quark decaying to a b-quark and a W -boson. The W -bosons are then further decayed to cs orcs. We avoid the W → ud mode as the probability of a light jet faking a b-jet is ∼ 30 times smaller than that of a c-jet. Lastly, we also consider the subdominant background coming from the non-resonant di-Higgs production (gg → hh) and also from ttbb (including ttZ/tth).
We select events containing exactly 4 b-tagged jets with the requirement of p T,b > 60 GeV and |η b | < 2.5. The scalar sum, H T , of the transverse momenta of all the visible particles in a event must fulfil, H T > 300 GeV. Finally, we form two di-jet systems from these four b-jets. The two jets within a dijet system must satisfy 0.4 < ∆R bb < 1.5. We choose the leading (sub-leading) di-jet system to have p T > 200 (150) GeV. Furthermore, to reduce the contamination from the tt background, we reconstruct the top by combining extra jets in an event with the di-jet systems. These jets must be within ∆R < 1.5 in the η − φ plane with the di-jet system. If an event contains exactly one extra jet, then we choose the di-jet system which is closest to it and combine to form a top quark system, m t 1 . However, when there are two such jets, we compute the minimum of all possible ∆R combinations between these two jets and the two di-jets before reconstructing two other top masses, m t 2 and m t 3 . Because for our signal, we do not expect any proper top quark reconstruction, we thus veto events if the reconstructed mass of any of these possible choices for the top quark exceeds 120 GeV. After imposing this cut the tt background reduces to half with more than 80% of the signal events still to spare. We detail these cuts one by one alongside the signal efficiency and cross-sections for the background processes in Table 7  Finally, after all the aforementioned cuts are applied in succession, we check for any possible improvement upon performing a multivariate analysis. We utilise the BDT algorithm for our purposes and choose the following nine kinematic variables with maximal potency, Here we use the kinematic variables reconstructed from the two di-jet systems viz., invariant mass (m di-jet ), transverse momentum (p T,di-jet ) and azimuthal angle separation between the b-jets forming the dijet systems (∆φ bb,di-jet ). The subscript k = 1, 2 refers to the p T ordering of the di-jets. We also take the separation in the η and η-φ plane between the two di-jets, viz., ∆η di-jets and ∆R di-jets respectively. m bbbb is the invariant mass of the four b-jet system. The top four variables with the highest discriminatory power are shown in Fig. 5. We can see that the lower masses have significantly longer tails while performing the mass reconstructions.
Finally, in Table 8, we present the background yields after the BDT optimisation has been completed. Like in the previous section, we translate these results into an exclusion diagram showing the upper limits on σ(pp → H → hh) as a function of the heavy Higgs mass. We show these in Fig. 6. The limit is very strong between 600 GeV and 1 TeV with the 95% CL upper limit varying between ∼ 23 fb and ∼ 4 fb.        Table 8: Respective background yields for the bbbb channel after the BDT analyses optimised for a heavy Higgs mass of (a) 400 GeV, (b) 600 GeV, (c) 800 GeV and (d) 1 TeV. The tables also list the perturbative order at which the cross-sections are considered.

Discussion about m H = 400 GeV
The 95% and 99.7% CL upper limits on σ(pp → H → hh) for the heavy Higgs with a mass around 400 GeV is very large (∼ 440 fb and ∼ 1100 fb respectively) as compared to the other mass points, even after the BDT optimisation. The reason for this is the following. The signal efficiency for m H = 400 GeV reduces by ∼ 95% after imposing the ∆R bb selection as can be seen from Table 7. Since the heavy Higgs mass (400 GeV) is near the threshold of the non-resonant di-Higgs production, the SM-like Higgs bosons for the resonant case are produced with low p T . This further leads to the Higgs decay products being widely separated in the η − φ plane and thus obviously does not satisfy our di-jet selection criteria of ∆R bb < 1.5 within each di-jet system. With the sole intention of improving the sensitivity, we adopt a χ 2 minimisation technique as described below. We define a new kinematic variable χ 2 hh for the events which do not satisfy the ∆R bb < 1.5 selection criteria as follows where m h = 125 GeV and σ hj = 0.1 × m di-jet,j with j = 1, 2 marks the p T ordering. Thus, in addition to the events satisfying ∆R bb < 1.5, we also consider those events which contain di-jet pairs with ∆R separation between the b-jets to be more than 1.5. Following this, we construct the aforementioned χ 2 hh variable for each possible pair of reconstructed di-jet. The event is finally selected if the non-zero minimum value of the χ 2 hh variable is less than 50 12 . Upon using this modification, the signal efficiency increases by ∼ 24% at the di-jet selection level while simultaneously increasing the dominant backgrounds like bbbb by ∼ 5% and tt by ∼ 13%. However, the limit on the upper limit of the cross-section improves to ∼ 340 fb and ∼ 855 fb at 95% and 99.7% CL respectively.

The bbτ τ Channel
Next, we turn our attention to one of the best probes for the di-Higgs searches, viz., the bbτ τ channel. The intricacy and potential of this channel lies in our ability to reconstruct the τ -leptons as these come with neutrinos which show up as missing transverse energy in the detector. This channel gives rise to three phenomenologically different final states, viz., bb In this work, we will only consider the last category, i.e., the one with the fully hadronic τ decays. The hadronically decayed τ -leptons are termed as τ -hadrons or τ -jets which may either contain one (one-pronged) or three (three-pronged) charged particle(s) inside the jet cone. Thus, it is essential to tag these tau-jets in order to segregate them from regular QCD jets ensuing from the various backgrounds that we will discuss below. We will not discuss the fully leptonic case here as from our previous analysis [100] we know that the sensitivity is extremely low even at the HL-LHC.
We generate the signal samples with Pythia 6. For the background samples, we employ MG5 aMC@NLO [144]. We generate two different samples for the dominant tt background, where either both the W -bosons decay to τ leptons or where one decays into a τ and the other to a pair of jets. The QCD-QED background bb also contributes significantly. Besides, we also generate the subdominant backgrounds which include tth, ttW , ttZ, bbh, Zh and the non-resonant Higgs pair production i.e., gg → hh. We simulate the Zh background upon considering two processes where in one case the Z-boson decays to a pair of bottom quarks and the Higgs boson decays to a pair of τ -leptons and in the other the decays are reversed. Moreover, we also generate the dominant fake background for the hadronic channel in the form of bbjj, where the light jets can be fake τ -tagged jets. We detail the generation level cuts for these various backgrounds in Appendix A. Following the generation level cuts, we further apply some basic cuts on the signal and background samples in order to ensure a common kinematic phase space. The b/τ -jets and the leptons (e, µ) are required to lie within |η| < 2.5 and have p T,b/τ ( ) > 20 (10) GeV. The light jets must satisfy p T,j > 20 GeV and |η j | < 4.5. The minimum distance in the η − φ plane between the b-jets and the leptons, and also among themselves is required to be ∆R > 0.2. The reconstructed invariant mass of the bottom pair and the visible τ pair must obey m bb/τ τ > 50 GeV.

The bbτ h τ h Channel
In this sub-section we briefly outline the prospects of searching for the heavy Higgs in the bbτ h τ h final state. In doing so, we select events containing exactly two b-tagged jets and two τ -tagged jets alongside the cuts described above. Having seen the strength of the multivariate analyses for this channel, in Ref. [100], rather than opting for the classical cut-based analysis, we perform a BDT analysis with the following 13 variables with the maximal discerning capability.
where M T is the transverse mass of the h → τ τ system 13 , p T,tot and m tot are respectively the transverse momenta and mass of the full visible system and m eff is the scalar sum of the transverse mass of all the visible products plus / E T . The rest of the variables have usual definitions. The top five variables are shown in fig. 7. As can be seen, the m T 2 variable is particularly useful for heavier Higgs masses as it can be used to completely eradicate the tt background. We train the signal and background samples and they are optimised for each benchmark signal point. We list the background events after optimising the BDT and imposing the cut for four values of m H , in Table 9. Finally, we show the upper limit on the heavy Higgs production cross-section (assuming BR(H → hh) = 100%) in Fig. 8. The 95% CL upper limit on the cross-section between m H = 600 GeV and 1 TeV varies between ∼ 100 fb and ∼ 40 fb.

The bbW W * Channel
In this section, we consider the situation where a heavy scalar decays to a pair of SM-like Higgs bosons with one of them decaying to a pair of b-quarks and the other to W W * , leading to three possible final states depending on the decays of the W -bosons. We perform our analyses for the fully leptonic (leptons at this stage include e, µ, τ ) and the semi-leptonic channels. We avoid studying the fully hadronic mode as the signal will be overwhelmed by the huge QCD background.
The dominant contribution to the background for both the channels mentioned above comes from top pair production because of its large production cross-section. We generate this background where either or both the W -bosons decay leptonically. The fully hadronic  Table 9: Signal and background yields for the bbτ h τ h channel after the BDT analysis for heavy Higgs mass of (a) 400 GeV, (b) 600 GeV, (c) 800 GeV and (d) 1000 GeV for bbτ h τ h channel.  tt mode is not considered as a potential background as the fake rate for j → is negligible for all practical purposes. The fully leptonic tt background contributes to the fully leptonic channel final state whereas for the semi-leptonic scenario, the contribution comes from both the fully-leptonic as well as the semi-leptonic tt. The second most dominant background for the semi-leptonic channel is W bb + jets, where the W -boson decays leptonically (e, µ, τ ). We generate this background upon merging with two additional jets by exploiting the MLM merging scheme [139]. While generating the W bb + jets background we ensure that there is no double counting ensuing from the semi-leptonic tt background. Besides the aforementioned backgrounds, a significant contribution also comes from the + − bb production where refers to e, µ and τ . Finally, we also consider the subdominant backgrounds viz., tth, ttZ, ttW and the non-resonant gg → hh.
In this subsection and the following section (section 3), the top-pair production is the dominant background. Thus, the reconstruction of the top quarks is a very powerful tool in order to reduce the contribution from this background. For the semi-leptonic case, the only source of missing transverse energy, / E T , arises from the neutrino of the leptonically decaying W -boson from the top decay. We reconstruct the top from its decay products 14 . The quadratic equation gives two possible solutions for the neutrino p z . Besides, because there are two b-jets in the final state, we get four possible choices for the mass of the leptonically decaying top. We use these variables during our analysis. After reconstructing both the tops, we reconstruct the total system from all the final state particles. We also use this later in section 3 which exhibits the same final state. These variables help us greatly in reducing the semi-leptonic tt background for high values of m H .
Before embarking on the final analysis, we impose a common set of trigger cuts for both the leptonic and the semi-leptonic channels. The p T , |η| and ∆R cuts for the various objects are discussed in subsection 2.3 and also in Appendix A. Furthermore, we require generation-level cuts on some invariant masses, viz., m bb > 50 GeV and m > 30 GeV (only for the fully leptonic channel). The cut m > 30 GeV is also applied at the analysis level and applies to lepton pairs of different flavours as well. The selected events are also require to have / E T > 40 GeV upon scrutinising the distribution. The / E T distribution for the 1 and 2 cases are shown in Fig. 9. Finally, we perform separate multivariate analyses for the two final states upon using the BDTD algorithm. While training samples for both the leptonic and semi-leptonic analyses, we only consider the tt background since it constitutes the bulk of the total background. This training is used for testing all other backgrounds which are subdominant in front of tt.
14 First the W -boson mass is reconstructed in order to attain the pZ component of the neutrino.

The 2 2b + / E T Channel
For the fully leptonic final state, we select events with exactly two b-tagged jets, and two isolated leptons having opposite charge meeting the trigger criteria as mentioned above. We choose the following set of kinematic variables in order to perform the multivariate analysis.
p T,bb , η bb , φ bb , m bb , ∆R bb , ∆φ bb , p T, , η , φ , m , ∆R , where M T is the transverse mass of the SM-like Higgs decaying to W -bosons. The rest of the variables have either been defined before or have usual meaning. The top four variables are shown in Fig. 10. The signal distributions are significantly different from the various backgrounds.
Finally, in Table 10, we summarise the number of background events after imposing the optimised cut on the BDT variable. Like in the other channels, we impose an upper limit on σ(pp → H → hh) as a function of the heavy Higgs mass. This is shown in Fig. 11. The 95% CL upper limit varies between ∼ 105 fb and ∼ 60 fb within 600 GeV < m H < 1 TeV and is somewhat weaker than the channels discussed earlier owing to smaller S/B.

The 1 2b2j + / E T Channel
Finally, we discuss the potential of the semi-leptonic final state as well. We require events with exactly two b-tagged jets, one isolated lepton and at least two light jets satisfying      the trigger criteria discussed earlier. Besides, we consider the same set of cuts as for the dileptonic channel (except the m cut) before performing the multivariate analysis. We find the following kinematic variables to have the best discriminatory power and use them for our multivariate analysis. η bb , m bb , m t , m jj , ∆R jj , ∆R ,jj , M T , m T 2 , m bbj 1 , m t11 , m t12 , p T, ν , p T,b 1 , p T, 1 , p T,j 1 , where m t is the transverse mass of the leptonically decaying W -boson. ∆R ,jj is the distance in the η − φ plane between the system comprising of the two hardest jets and the lepton. m bbj 1 refers to the invariant mass of the two b-tagged jets and the hardest p T jet. The reconstructed transverse momentum of the leptonically decaying W -boson is denoted as p T, ν . m tij is the mass of the leptonically decaying top quark with the reconstruction procedure outlined before. The first index i = 1, 2 indicates the p T ordering of the bjet. The second index j = 1, 2 refers to the choice of the z-component of the neutrino momentum. The other variables have usual definitions. The best discriminatory variables are listed in Fig. 12. However, we can see that the separation power for the 1 category is significantly less compared to its 2 counterpart.
Coming to the results, table 11 summarises the background yields after the BDT cut. The upper limit on σ(pp → H → hh) as a function of m H are shown in Fig. 13. The limits are considerably weak in this channel.

The γγW W * Channel
After the bbγγ channel this is the second most cleanest channel in terms of the final state particles but with the pitfall of having very low event rate. In this channel, one of the SMlike Higgs decays to a pair of photons and the other to lepton(s) through h → W W * . Similar to the bbW W * analysis in subsection 2.4, here also we divide the channel into the leptonic and semi-leptonic category. Because of the relatively clean final states, these channels have low contaminations due to backgrounds. We simulate the Zh and W h backgrounds upon merging with two additional jets (the definition of jet is given in subsection 2.1). Here we decay the Z-and the W -bosons leptonically (e, µ, τ )).  Table 11: Respective background yields for the 1 + 2j + 2b + / E T channel after the BDT analyses optimised for m H = (a) 400 GeV, (b) 600 GeV, (c) 800 GeV and (d) 1 TeV. The various orders of the signal and backgrounds are same as in Table 10.  an additional jet and using the same scheme as before. Next, we also consider the tth background with Higgs-boson decayed to a pair of photons. Finally, we also consider the SM Higgs pair production which is subdominant.
Before performing the multivariate analyses, we impose the generic trigger cuts. The p T , |η| and ∆R 15 cuts are the same as has been defined in subsection 2.3. The above cuts for the photons are the same as those on the leptons. Owing to an excellent resolution for the diphoton invariant mass, we require 122 GeV < m γγ < 128 GeV. Finally, we also require m > 30 GeV because we generate the γγ background with this invariant mass cut at the generation level (the details of these cuts are mentioned in Appendix A). We now describe the results of the multivariate analyses for these two channels in the following two subsections.

The γγ1 2j + / E T Channel
Before performing the BDT analysis, we select events with exactly two isolated photons, one isolated lepton and at least two jets in the final state, which fulfils all the aforementioned trigger requirements. Like all the other channels, we consider the following variables to train our signal and background samples for the multivariate analysis.
where the variables carry their usual meaning. The five best variables are shown in Fig. 14. The background yields after the BDT optimisation are shown in Table 12. In Fig. 15, we show the upper limit on σ(pp → H → hh) as a function of m H . The 95% CL upper limit changes from ∼ 205 fb for m H = 400 GeV to ∼ 135 fb for m H = 1 TeV.

The γγ2 + / E T Channel
This is the final channel that we study for the pp → H → hh case. We choose events with exactly two isolated photons, and two isolated leptons with opposite charge, following the trigger cuts mentioned above. Finally, we choose the following kinematic variables for the multivariate analysis.
p T,γγ , ∆R γγ , ∆φ γγ , m , ∆R , M T , m tot , m eff , ∆R γγ, , p T, 1 , with the usual definitions for the variables. Some of the variables of interest are plotted in Fig. 16. The background yields after the BDT cut are tabulated in Table 13 whereas the upper limit on σ(pp → H → hh) as a function of heavy Higgs mass is shown in Fig. 17. The 95% CL upper limit for the leptonic scenario is stronger than its semi-leptonic counterpart in the heavy Higgs mass range of 600 GeV and 1 TeV. The upper limit varies in between ∼ 155 fb and ∼ 75 fb at 95% CL, in the aforementioned range. 15 ∆R γγ/γ > 0.4 and ∆R >0.2 .

Summarising the H → hh channel
Having studied five different channels with more than one sub-processes in three instances, we summarise the results in this subsection. The 95% CL upper limits on σ(pp → H → hh)     for all these channels is shown in Fig. 18. We find that the strongest limits come from the bbγγ and 4b channels. The bbγγ is strongest up to m H ∼ 660 GeV. From 660 GeV onward, the 4b channel is more constraining owing to its larger cross-section.  TeV) from the 4b channel. We find an order of magnitude improvement in the sensitivity.

The pp → H → tt Channel
After having studied the H → hh in multifarious channels in detail, we now turn our attention to a heavy scalar (or pseudoscalar) resonance being produced predominantly by gluon  Table 13: Respective background yields for the γγ2 + / E T channel after the BDT analyses optimised for m H = (a) 400 GeV, (b) 600 GeV, (c) 800 GeV and (d) 1 TeV. The various perturbative orders for the backgrounds are the same as in Table 12. The branching ratio of t → bW being close to 100% makes the channel essentially become a search for H → bbW + W − . However, unlike the H → hh → bbW W * channel studied in subsection 2.4, where one of the W -bosons is off-shell, here both of them are on-shell. This is the first essential difference between the two channels and the reason why one requires a completely different search strategy for the two cases. In the previous section 2,  Figure 18: 95% CL upper limit on σ(pp → H → hh) (fb) as a function of m H (GeV) for the bbγγ, bbbb, bbτ + τ − , bbW W * and W W * γγ channels.

95% C.L Upper Limit
we required BR(H → hh) = 100%. However, in realistic scenarios, if the heavy scalar is produced predominantly via gluon fusion (top/bottom loops), it should also decay to a pair of top quarks (and also bottom quarks) if it is above the tt threshold. Similar to the H → hh → bbW W * channel, here also we divide the analysis into two parts, viz., the leptonic and the semi-leptonic channels. We apply the same trigger-level cuts to the various objects as sketched in subsection 2.4. The backgrounds are the same as before. As before, we implement the production and decay of the heavy scalar in the Pythia 6 framework.

The leptonic Channel
Like in section 2.4.1, here also we select events with two oppositely charged isolated leptons and two b-tagged jets. Without performing a classical cut-based analysis, we optimise our results to obtain the best-possible sensitivity by employing a boosted decision tree analysis. The set of variables which discriminate the signal from the backgrounds are as follows.
where all the variables have their usual meaning as mentioned earlier. The top four discriminatory variables are shown in Fig. 19. In Table 14, the number of background events at an integrated luminosity of 3000 fb −1 , optimised to maximise the sensitivity for various values of m H and after imposing a cut on the BDT observable, are presented. Like in all the other channels, we present the 95% and 99.7% upper limit on σ(pp → H → tt) as a function of m H , in Fig. 20. We find that the 95% upper limit on the cross-section lies between ∼ 460 fb and ∼ 160 fb for m H varying between 400 GeV and 1 TeV.

The semi-leptonic Channel
We end this section by analysing the semi-leptonic final state ensuing from the semi-leptonic decays of tt. We select events which contain a single isolated lepton, two b-tagged jets and at least two light jets after applying the same set of trigger cuts as discussed in section 2.4. Finally, we perform a multivariate analysis with the following set of kinematic variables.   Table 14: Respective background yields for the 2 + 2b + / E T channel after the BDT analyses optimised for m H = (a) 400 GeV, (b) 600 GeV, (c) 800 GeV and (d) 1 TeV. The various orders of the signal and backgrounds are same as in Table 10.  summarise the boosted decision tree results in Table 15. For heavy Higgs mass ranging between 400 GeV and 1 TeV, we show the upper limit on σ(pp → H → tt) in Fig. 22 Table 15: Respective background yields for the 1 2j2b + / E T channel after the BDT analyses optimised for m H = (a) 400 GeV, (b) 600 GeV, (c) 800 GeV and (d) 1 TeV. The various orders of the signal and backgrounds are same as in Table 10.

The (H/A)bb channel
Finally, we study the process where the resonant (pseudo)scalar is produced in association with a pair of bottom quarks, viz., pp → (bb)H/A. The need to study this process lies in the fact that one can probe and impose strong limits on the lower part in the m A − tan β plane, as will be discussed in section 5. The cross-section of the inclusive (bb)H process receives contribution from both the 5-flavour (5F) and the 4-flavour (4F) processes. At leading order (LO), the 5F scheme is dominated by the process bb → H. However, for scenarios involving b-tagging in the final state, the process where the resonant scalar is produced with b-quark or gluon becomes important, viz. gb → g/b H. There are two processes contributing to the 4F scheme at LO where the heavy Higgs is produced in association with two b-quarks, one is via the gluon fusion process (gg → bbH) and the other is via quarks (qq → bbH). The total inclusive cross-section is obtained by matching the 4F and 5F scheme numbers in which both these cross sections are multiplied by proper weight factors [148]. The matched cross-section is computed as follows.
is the weight factor. The analysis can be subdivided according to the number of b-tagged jets. However, we specifically focus on the category with ≥ 1 b-jets upon following a recent study performed by the ATLAS collaboration [128]. Furthermore, we consider the heavy Higgs decaying to a pair of τ -leptons and we specifically focus on the scenario where both the τ s decay hadronically.
On the analysis front, we generate the signal events (at LO in SM) with MG5 aMC@NLO and shower them via Pythia-8 [149]. We use different parton distribution functions (PDFs) for the sample generations. Specifically, we use the CT10nlo nf4 [150] for the 4F bbH process, MSTW2008nnlo68cl [151] for the 5F bbH process and CT10 [152] for the ggF process. The cross-sections for the 5F scheme as well as for the ggF (at NNLO) are estimated using SusHi [153]. Whereas, for the 4F scheme we use the MadGraph LO cross-sections. The H and A masses are varied between 200 GeV and 1 TeV.
The various backgrounds at play are Z/γ * + jets, multijets, W + jets, V V (V = W ± , Z), tt and single top. The Z/γ * + jets with the Z-boson decaying to a pair of leptons (e, µ and τ ) is the dominant background for the τ h τ category (a category that we will not address in the present work) but also gives significant contribution to the τ h τ h category. We simulate this background merged with three additional partons and some specific generation level cuts which are described in Appendix-A. Similarly, the W + jets is also generated with up to three additional partons and the W -boson is then decayed leptonically. In order to include the dominant multijets background in the τ h τ h category, we generate an exclusive bbjj sample where j includes light quarks and gluon. These light jets can fake hadronically decaying τ s. Finally, we include the top-quark related backgrounds viz., tt and single top. Next, we describe our analysis for the ≥ 1 b-tagged jets category upon closely following Ref. [128].

The τ h τ h Channel : b-tag category
We select events containing at least one b-tagged jet with p T > 20 GeV and two τ -tagged jets with p T > 65 GeV 16 . These two τ -tagged jets must have opposite electric charge (from their track reconstruction). We also veto events having leptons (e, µ) or τ -tagged jets with 1.37 < |η τ | < 1.52, in the final state. The azimuthal angle separation between the two τ -tagged jets has to fulfil the condition, |∆φ(τ, τ )| > 2.7. The b-and the τ -jets must have an angular separation in the η − φ plane, viz., ∆R(b, τ ) > 0.2. Besides, we also impose a minimum bound on the visible invariant mass of the two hadronically decaying τ leptons to be m vis. τ τ > 50 GeV. For the fake bbjj background, we demand the two light jets to satisfy the τ jet configuration during the analysis and we later multiply the event yield with the j → τ fake rate. Similarly, for the W (→ τ ν or ν)+ jets background, we demand at least one extra light jet satisfying the τ jet requirement on top of the τ jet ensuing from W -boson decay. After imposing the aforementioned cuts, we improve our analysis by optimising over some other kinematic variables viz., the transverse momentum of the hardest τ -tagged jet (p T,τ 1 ), sum of the cosine of the azimuthal angle separation between the τ -jets and / E T ( τ 1,2 cos ∆φ) and the transverse mass of the total system which is defined below.
where the symbols have their usual meanings. These four kinematic variables are shown in Fig. 23. The optimised cuts along with the signal efficiencies and the background yields for each benchmark point are shown in Table 16. We perform our analysis upon considering both the 4F and 5F signal samples separately. Finally, we add add them by multiplying these cross-sections with the aforementioned weight factor in order to obtain the upper limit on the matched bbH production cross section. We show the 95% exclusion for σ(pp → bbH) × BR(H → τ h τ h ) in Fig. 24. The 95% upper limit varies between ∼ 15 fb and ∼ 3 fb for m H varying between 300 GeV and 500 GeV. This is close to an order of magnitude improvement over the existing bounds at 13

The future of the pMSSM parameter space
The Higgs sector in the MSSM comprises of two Higgs doublets which give rise to five massive Higgs states after the electroweak symmetry breaking. The Higgs spectrum is thus composed of two CP -even scalars, h and H, one CP -odd scalar, A, and two charged scalars, H ± (detailed studies on the Higgs sector of MSSM can be found in Ref. [36,154]). In addition to the extended Higgs sector, the SUSY particle spectrum boasts a multitude of particles, viz., the sleptons, squarks, gluinos and electroweakinos. A majority of direct searches at the LHC have excluded stops and gluinos below the TeV scale (Refs. [155][156][157][158][159][160][161][162][163] show some such limits in various supersymmetric interpretations). This more or less nullifies the prospect of observing these particles unless the luminosity is enhanced significantly. The electroweakino sector has also been probed in numerous studies [164][165][166][167][168][169] and bounds have been obtained on their masses within simplified scenarios [170][171][172][173][174][175][176][177][178][179]. In many of these studies, the electroweakino masses are excluded from between a few hundred GeVs to about half a TeV and are comparatively weakly coupled compared to the gluinos and stops. Within a generic SUSY parameter space without any correlation between the choice of the electroweakino mass parameters, these bounds can become considerably weaker. The ATLAS and CMS collaborations have also performed several studies to search for resonant Higgs through their decay into SM final states [57,110,128,130,[180][181][182][183][184][185]. However, none of these searches could find any significant excess over the SM expectations and thus only imposed upper limits on the production cross-section of the heavy Higgs bosons times their branching ratio into various SM final states. In this section, we present where σ i represents the MSSM (or any specific model in question) Higgs production cross-section in the i th production mode (i = ggF , V BF , tth or V h) at the LHC and σ SM i denotes the corresponding SM cross-section. BR f corresponds to the branching fraction of the MSSM Higgs into a particular SM final state (f = W W, ZZ, bb, γγ, τ τ ) and BR SM f is the corresponding SM value. We apply all these constraints over our parameter space by demanding that all our signal strengths simultaneously lie within 2σ of their experimental counterparts. The latest Higgs signal strengths (13 TeV, 36 fb −1 ) measured by both the CMS and ATLAS collaborations are listed in Table 17.  In order to evaluate the current allowed region in the parameter space of the phe-nomenological MSSM (pMSSM), we perform a random scan over a wide range of pMSSM input parameters, as described below. The parameters relevant to our study are the pseudoscalar mass variable (m A ), ratio of the vacuum expectation values of the two Higgs doublets (tan β), the third generation soft squark mass parameters (MQ 3 , Mũ 3 , Md 3 ), the trilinear coupling of the stop (A t ) and sbottom (A b ) and the gluino mass parameter (M 3 ). These parameters are varied in the following range.
The particle spectra and the branching fractions of the SM and SUSY particles are obtained using FeynHiggs [194]. We consider only those parameter points which satisfy the light Higgs mass constraint defined above. Furthermore, we allow only those points which lie within 2σ uncertainty of each of signal strength variables listed in Table 17. The parameter space points which are allowed by the aforementioned light Higgs mass constraint and the Higgs signal strength constraints are referred to as the allowed parameter space points in the remainder of this section and are shown in grey in Fig. 26.
The ATLAS and CMS collaborations have also performed numerous searches for the heavy Higgs bosons through their decay into SM final states, however, none of these searches have been able to observe any significant excess over the SM expectation. Consequently, upper limits have been set on the production cross-section of the heavy Higgs boson (σ H/A ) times its branching ratio into SM states. In this analysis, we consider the latest search limits on σ ggH ×Br ( [110,130,184] derived by the CMS and ATLAS collaborations upon using the Run-II dataset with an integrated luminosity of ∼ 36 fb −1 . The gluon fusion channel is undoubtedly the dominant Higgs production mode at the LHC for low values of tan β. However, it gets overrun by the bbH/A production channel at high tan β values. In the current analysis, while evaluating the impact of the existing upper limits on pp → H → hh → 4b, 2b2γ, 2b2τ , only the contributions from the gluon fusion production are taken into account. This choice is motivated by the fact that the the H → hh decay modes gain dominance only in the low and intermediate tan β values where the gluon fusion mode overshadows the bbH/A channel. Although the current search limits on H → hh do not impose any constraints on our parameter space, the future runs have the potential to probe the low m A and low tan β regime. The impact of these future limits are discussed in the later part of this section. The H → ZZ/W W limits also turn out to be ineffective in constraining our parameter space and will require improvements of about three orders of magnitude for making any impact. We would like to mention that the upper limits derived by ATLAS in the H → γγ search channel is on the fiducial crosssection times BR(H → γγ). We compare these upper limits against a combination of the ggF + bbH/A production cross-sections and observe that an improvement of around two orders of magnitude will be required in order to affect our parameter space. Limits from searches in the H/A → τ τ channel impose the strongest constraints on the parameter space. Constraints from σ bbH/A × BR(H/A → τ τ ) yield stronger limits compared to their gluon fusion counterparts and exclude the low m A and high tan β region. The current search limits from ATLAS and CMS furnish roughly equivalent impact and rule out tan β 16 for m A ∼ 1 TeV. Before presenting the results in the m A −tan β plane, we show the current allowed branching fractions, viz., H → hh, H → tt and H → τ + τ − in Fig. 25. The H → hh branching ratio dominates for tan β 8 and for m H ≤ 2m t . All points are allowed by the Higgs mass and Higgs signal strength constraints. However, the grey regions are excluded by the present direct searches of the heavy Higgs. In Fig. 26, we show the impact of the latest direct search limits from σ bbH/A × Br(H/A → τ τ ) in the m A − tan β plane. The parameter space points shown in Fig. 26 (grey and orange) are obtained after implementing the light Higgs mass constraints and the Higgs signal strength measurements. The grey points are excluded upon imposing the aforementioned direct search limits.
Our main concern in this section is to quantify the impact of the projected direct search limits for the HL-LHC which were derived in the previous sections. In this regard, we consider the projected direct search limits for the HL-LHC in the H → hh (Sec. 2), H → tt (Sec. 3) and bbH/A → bbτ h τ h (Sec. 4) channels. Among the various final states of the H → hh channel, the bbγγ final state furnishes the strongest limit in the m A 600 GeV regime, while the 4b final state imposes the strongest upper limits in the m A 700 GeV region. The H → hh decay mode gains dominance in the low tan β region and especially before the tt mass threshold is attained. The same is reflected in the left panel of  where the brown points represent the region excluded at 95% CL by the projected 2σ reach from the H → hh → bbγγ channel. The 4b final state, on the other hand, is rendered ineffective on account of reduced production cross-section at high values of m A . The upper limits derived from searches in the remaining H → hh channels furnish much weaker bounds and will not be able to probe the pMSSM parameter even at the HL-LHC. The couplings of the heavy Higgs bosons with the up-type quarks have an inverse dependence on tan β and thus consequently the H → tt channel has the potential to play an important role in the low tan β regime. The parameter space points excluded at 95% CL by the H → tt HL-LHC search limits derived in Sec. 3 are shown in green in Fig. 27. The strongest future limits are obtained by the bbH → bbτ τ channel (derived in Sec. 4). This will be able to exclude (at 2σ) until ∼ tan β ∼ 8 at m A ∼ 1 TeV as shown in Fig. 27, where the orange points are excluded by the same. The blue points in Fig. 27 denotes the parameter space which will evade the direct searches at the HL-LHC as well. The right panel in Fig. 27, however, shows the discovery potential at 5σ.
The results discussed till now assume that the heavy Higgs bosons underwent decays only into SM final states. The branching fractions of the heavy Higgs bosons into SM final states can, in principle, undergo significant modifications in the presence of light SUSY  [195,196]. In the remainder of this section we study the impact on the pMSSM parameter space in the presence of non-SM decay modes of the heavy Higgs bosons at the HL-LHC. Here we will restrict ourselves to the case of light electroweakinos. These electroweakinos are required to be an admixture of gauginos (bino and wino) and higgsinos in order to have couplings with the Higgs bosons. The pure gauginos and higgsinos do not couple with the Higgs states. LHC searches in the chargino-neutralino pair production mode furnishes the most stringent constraints on the electroweakino sector and excludes degenerate wino-like χ 0 2 and χ ± 1 of mass 450 GeV [179], for an LSP neutralino of mass ∼ 100 GeV. However, such constraints do not apply to scenarios where the LSP and NLSP are almost degenerate in mass. We explore this fact and vary M 2 and µ in such a way that |M 2 − µ| < 10 GeV, M H/A > (χ 0 1 + χ 0 2 ) and M 2 , µ > 200 GeV [179]. Closeness between M 2 and µ ensures that the χ 0 1 , χ 0 2 , χ 0 3 and χ ± 1 have significant admixtures from both winos and higgsinos. The remaining input parameters are randomly varied within the range specified in Eqn. 5.2 except for M 1 which we fix at 1 TeV. In presence of these H/A → ino decay modes, the branching fraction of the heavy Higgs bosons to SM final states undergoes modifications and manifests in weaker limits on the parameter space, as shown in Fig. 28. Correspondingly, the orange and blue regions in Fig. 27 shift upward. The brown and green regions shrink further down. For m A varying between 400 GeV and 700 GeV, tan β as low as 3 is excluded at 95% CL. In presence of these non-SM decay modes, the current (13 TeV, 36 fb −1 ) limits on σ bbH/A × Br(H/A → τ τ ) exclude tan β 22 for m A ∼ 1 TeV.
The HL-LHC reach weakens out till tan β ∼ 11 at m A ∼ 1 TeV. In the current scenario, the projected limits from H → tt lose sensitivity on the pMSSM parameter space under study. The HL-LHC projections from H → hh → bbγγ also imposes weaker constraints and excludes tan β < 6 at m A ∼ 400 GeV.

Discussion on systematic uncertainties
In case of the bbγγ channel, the signal over background ratio, S/B, is large and hence the upper limit is mildly affected by incorporating a systematic uncertainty as large as 10%. Upon its inclusion, the 2σ upper limit changes from [143, 37.04] fb to [179. 87, 40.26] fb for m H ∈ [275, 400] GeV. For bbH, the effects of systematic uncertainties become negligible for m H > 400 GeV. However, for lower masses, the backgrounds are larger and thus the inclusion of such uncertainties weaken the limits. The 2σ upper limit changes dramatically from ∼19 fb to ∼205 fb for m H = 300 GeV. However, at m H = 400 GeV, the limit changes from 15.12 fb to 15.18 fb. Since the H → tt channel has a small S/B ratio, adding systematic uncertainty will drastically change the upper limit on the cross-section.

Summary
In this work, we have studied the prospects of observing or excluding a resonant heavy Higgs or pseudoscalar in the purview of the HL-LHC at 14 TeV. Various experimental observations and theoretical motivations necessitate physics beyond the Standard Model (BSM). Several searches are performed, either in the context of specific models or in a generic model-independent fashion, to gauge the type of new physics. In this work, we specifically focus on neutral heavy Higgs bosons (both CP odd and even). Run-II data at the LHC has already constrained strongly interacting BSM particles like gluinos and stops to O(≥ 1) TeV. However, the LHC still hasn't imposed such strong constraints on extended Higgs sectors. The Standard Model Higgs self-coupling being still unknown, we are yet to fully understand the scalar sector of new physics. Here, we studied three major search channels for such a heavy Higgs (or pseudoscalar). Specific to the CP -even heavy Higgs, we studied the prospects of constraining σ(pp → H → hh) in multifarious channels, viz., bbγγ, bbbb, bbτ + τ − , bbW W * and γγW W * . We took guidance from the present searches and optimised each channel carefully to obtain upper limits on σ(pp → H → hh) from each channel. Corroborating the present searches, we find that the bbγγ and bbbb final states serve as the golden channels for m H in the range [300, 600] GeV and [700, 1000] GeV, respectively. The bbγγ sets a 95% CL upper limit on σ(pp → H → hh) between [104,17] fb in the aforementioned mass range. The 4b channel on the other hand sets a corresponding cross-section limit between [8, 4] fb for m H ∈ [800, 1000] GeV. The limits from the remaining three channels are not so promising. On the other hand, if the mass of the scalar or the pseudoscalar Higgs is above the tt threshold and the Higgs is dominantly produced via gluon fusion, then it can also have a dominant decay into tt. Upon studying the fully leptonic as well as the semi-leptonic final states, we find the strongest limits on σ(pp → H → tt) from the semi-leptonic category. The 95% CL upper limits lie between 210 fb and 40 fb for m H ∈ [400, 1000] GeV. Finally, we studied the bbH/A production in the bbτ τ final state upon demanding at least one b-tagged jet in the final state and demanding two hadronic τ s. The 95% CL upper limit on σ(pp → bbH) varies between 15 fb and 3 fb for m H lying between 300 GeV and 500 GeV. All our searches for the scalars are mostly model independent and can be translated to models with multiple Higgs bosons with narrow widths.
In this work, we considered the specific example of supersymmetry, more specifically the pMSSM. We apply present constraints from the SM-like Higgs boson mass measurement and all its signal strengths into multiple final states. The future limits obtained in this work constrain different regimes of the parameter space. The H → hh search, mostly in the bbγγ channel excludes tan β to as low as 4, at 95% CL, for m A ∼ 2m t GeV. The H → tt has a similar exclusion on tan β for m A varying between [400,800] GeV. The bbH channel in the di-τ + ≥ 1b-tagged jet final state excludes tan β as low as 6.5 for m A = 1 TeV. The blue region in Fig. 27 will not be probed even by direct searches if the heavy Higgs bosons decay only to SM particles. We will require higher energy colliders in order to be able to probe this region. This scenario might change if there are light electroweakinos, below m H/A /2. In that situation, the m A − tan β parameter region changes. Upon considering a scenario where |M 2 − µ| < 10 GeV, M H/A > (χ 0 1 + χ 0 2 ) and M 2 /µ > 200 GeV, one finds that the H → hh and H → tt excludes tan β down to 3 for m A ∈ [400, 700] GeV. The exclusion bound on tan β decreases to 11 for m A = 1 TeV at 95% CL.