Enlarging the scope of resonant di-Higgs searches: hunting for Higgs-to-Higgs cascades in 4b final states at the LHC and future colliders

We extend the coverage of resonant di-Higgs searches in the bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{b} $$\end{document}bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{b} $$\end{document} final state to the process pp → H1→ H2H2→ bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{b} $$\end{document}bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{b} $$\end{document}, where both H1,2 are spin-0 states beyond the Standard Model. Such a process constitutes a joint discovery mode for the new states H1 and H2. We present the first sensitivity study of this channel, using public LHC data to validate our analysis. We also provide a first estimate of the sensitivity of the search for the HL-LHC and future facilities like the HE-LHC and FCC-hh. We analyze the discovery potential of this search for several non-minimal scalar sector scenarios: an extension of the SM with two extra singlet scalar fields, the two-Higgs-doublet model and a two-Higgs doublet model plus a singlet, which captures the scalar potential features of the NMSSM. We find that this channel represents a novel, very powerful probe for extended Higgs sectors, offering complementary sensitivity to existing analyses.


Introduction
While analyses at the Large Hadron Collider (LHC) by ATLAS and CMS show that the properties of the Higgs particle h with mass m h ∼ 125 GeV are, at present, compatible with those of the Standard Model (SM) Higgs boson h SM [1][2][3], the detailed nature of the scalar sector responsible for electroweak (EW) symmetry-breaking remains to be determined. It is particularly important to ascertain whether the scalar sector consists of just one SU(2) L doublet or has a richer structure with additional states. Addressing this question is a key goal of present and future studies at the LHC.
Searching for the existence of additional Higgs bosons at the LHC constitutes the main avenue for probing non-minimal scalar sectors, allowing to directly access the spectrum and properties of the scalars beyond the SM (BSM). Among such direct searches, those targeting decay chains involving several scalar states (henceforth Higgs-to-Higgs) are of particular importance. They depend on the scalar self-couplings and could therefore provide insight into the structure of the scalar potential. Resonant di-Higgs production, pp → H → h SM h SM , is -1 -JHEP02(2020)002 the prime (and simplest) example of a Higgs-to-Higgs process, where a resonantly produced BSM H state decays into a pair of 125 GeV Higgs bosons h SM (see [4,5] for a review). ATLAS and CMS have looked for this process at √ s = 8 TeV and 13 TeV in a wide range of final states, including bbbb [6][7][8][9], bbW + W − [10,11], bbτ + τ − [12][13][14] and bbγγ [15][16][17][18].
Non-minimal Higgs sectors generically feature several BSM states. In such a case, Higgs-to-Higgs decays with both the parent particle and its decay products as BSM states are possible, and may constitute the most promising avenue for their discovery. This has been emphasized in the literature for certain processes within the two-Higgs-doublet model (2HDM) [19][20][21][22][23], the 2HDM plus a scalar singlet [24], the next-to-minimal supersymmetric Standard Model (NMSSM) [25][26][27][28][29] and the SM extended by several singlet scalars [30].
In this work we show that it is possible to enlarge the scope of resonant di-Higgs pp → H → h SM h SM searches to probe more general Higgs-to-Higgs processes involving two BSM states. We present a detailed sensitivity study of the channel pp → H 1 → H 2 H 2 → bbbb where the heavier state H 1 is assumed to be produced via gluon fusion with a subsequent on-shell decay into a pair of H 2 bosons (i.e we consider m H 1 > 2m H 2 ). The potentially dominant H 1 → H 2 H 2 branching fraction for m H 1 m H 2 combined with a large H 2 → bb branching fraction 1 typical of light scalars (which would at the same time make the discovery of H 2 via direct production challenging, see e.g. [31]) make this search channel an important, yet unexplored, probe of non-minimal Higgs sectors. While no ATLAS or CMS analysis of the pp → H 1 → H 2 H 2 → bbbb signature exists at present 2 , we can use its similarity to resonant di-Higgs searches in the bbbb final state to validate our analysis for m H 2 = 125 GeV, before extending it to the 2D mass parameter space (m H 1 , m H 2 ). Specifically, we follow the recent √ s = 13 TeV CMS search for a narrow spin-0 or spin-2 di-Higgs resonance in the bbbb final state with 35.9 fb −1 of integrated luminosity [8]. The detailed public information available for this search allows us to reproduce their reported selection efficiencies and 95% confidence level (C.L.) exclusion sensitivities with our simulation. We then obtain the expected signal efficiencies for the pp → H 1 → H 2 H 2 → bbbb process in the mass plane (m H 1 , m H 2 ) and provide the 13 TeV LHC 95% C.L. exclusion sensitivity on the signal cross section with 35.9 fb −1 , as a function of m H 1 and m H 2 .
In addition to the above, we provide an extrapolation of the current exclusion sensitivity for the High-Luminosity (HL)-LHC with √ s = 14 TeV, as well as to future collider proposals like the High-Energy (HE)-LHC with √ s = 27 TeV and a √ s = 100 TeV protonproton collider (henceforth referred to as FCC-hh). We discuss the impact of multi-jet background systematic uncertainties on the expected sensitivity of our proposed search, as well as the effect of possible future improvements on the b-tagging and trigger efficiencies.
Finally, we assess the reach of the proposed analysis within specific BSM models: a two-singlet extension of the SM, a 2HDM scenario and a 2HDM plus a real scalar or pseudoscalar singlet. This allows us to compare the projected sensitivity of the search to other analyses targeting BSM scalars, studying their complementarity and identifying where the search pp → H 1 → H 2 H 2 → bbbb provides the leading probe of the existence of the H 1 and H 2 states. 1 A concrete example of such a scenario would be a 2HDM with a large mass splitting mH mA between the CP-even (H) and CP-odd (A) neutral BSM scalars. 2 We note there are existing LHC analyses for hSM → H2H2, with mH 2 < 62 GeV, see e.g. [32,33].

JHEP02(2020)002
Our work is organised as follows: in section 2 we reproduce the efficiencies of the CMS √ s = 13 TeV resonant di-Higgs search in the bbbb final state, in order to validate our subsequent analysis. In section 3 we derive the present 95% C.L. exclusion sensitivity prospects for the pp → H 1 → H 2 H 2 → bbbb process, and in section 4 we provide extrapolations to the HL-LHC, the HE-LHC and the FCC-hh. In section 5 we cast these prospects into a twosinglet extension of the SM, a 2HDM scenario and a 2HDM + singlet scalar/pseudoscalar, comparing in each case the sensitivity of our proposed search with other LHC searches for BSM scalars. Finally, we summarize our results in section 6.

Implementation and validation of the CMS search
As stated in section 1, there is no current experimental search at the LHC for BSM spin-0 states H 1,2 (i.e. belonging to an extended Higgs sector) through the process pp → H 1 → H 2 H 2 . However, the process bears a strong similarity to resonant di-Higgs production pp → H 1 → h SM h SM , which has been actively searched for by ATLAS and CMS since LHC Run 1. This similarity can thus be exploited to extend current resonant di-Higgs searches to include processes where both scalars belong to a BSM sector. In analogy with resonant di-Higgs, different analysis strategies can be devised, depending on the decay channels of the H 2 scalar. In this work, we concentrate on the bbbb final state and make use of the latest √ s = 13 TeV CMS search for a narrow di-Higgs resonance in this channel with 35.9 fb −1 of integrated luminosity [8]. 3 In this section, we validate our implementation of the CMS analysis by fixing m H 2 = 125 GeV and reproduce both the signal selection efficiencies and the 95% C.L. cross-section upper limits reported in [8], before extending our analysis to the 2D mass plane (m H 1 , m H 2 ) in section 3.

Validation of the selection efficiencies for the signal
The CMS collaboration reports the signal efficiencies at various stages of the spin-0 analysis event selection, namely from trigger level up to the signal region (SR) definition. The search defines two kinematic regions that feature different event selection criteria: a lowmass-region (LMR) for masses m H 1 ∈ [250, 650] GeV, and a medium-mass-region (MMR) for masses m H 1 ∈ [550, 1200] GeV. 4 The transition region m H 1 ∼ 580 GeV is determined by the respective sensitivities of the LMR and MMR selection strategies [8].
Events are selected with an online trigger that requires either of the following conditions to be satisfied: i) Four reconstructed jets of p T > 30 GeV and |η| < 2.4 of which two satisfy p T > 90 GeV and at least three b-tagged jets.
ii) Four reconstructed jets of p T > 45 GeV of which at least three are b-tagged. 3 We choose the CMS analysis [8] and not the ATLAS analysis [9] since in the former the public information available and the optimization for a low mass resonance made the validation and the extrapolation of their results easier to perform. Our strategy could nevertheless be applied to the ATLAS analysis leading to similar results. 4 For mH 1 > 1200 GeV, the angular separation between the two b-quarks from a Higgs decay is typically too small to satisfy jet isolation criteria, causing a large drop in the signal selection efficiencies. A different analysis strategy making use of jet-substructure techniques is needed in this regime.

JHEP02(2020)002
The analysis then requires all four selected jets to be b-tagged 5 and be within |η| < 2.4. This initial selection stage is labelled 4b and is common to both the LMR and MMR categories. For the LMR selection the analysis then identifies two 125 GeV Higgs boson candidates by pairing the b-jets and requiring |m bb − 120 GeV| < 40 GeV for each pair, while for the MMR selection the two b-jet pairs must satisfy ∆R bb < 1.5. This selection is named "HH candidate". 6 Finally, the SR is defined in the two dimensional space of the reconstructed masses of the lighter Higgs boson candidates, m H 1 2 and m H 2 2 , as the circular region with χ < 1, where χ is defined as The values of the parameters C and R are set to (C, R) = (120, 20) GeV and (C, R) = (125, 20) GeV for the LMR and MMR category respectively. We note that further improvements in the Higgs boson mass resolution through multivariate regression techniques applied by the CMS analysis are not included in our analysis. These increase the sensitivity of the CMS analysis by 5-20% depending on the mass hypothesis [8], and thus our validation is expected to yield a potential mismatch of at least that order. For our validation we have implemented the relevant interactions for the spin-0 BSM state in the Feynrules package [36] and generated hard-scattering events through the Madgraph5 aMC@NLO platform [37]. These events have been generated at leading order (LO) with fixed widths of 10 GeV and 1 GeV for H 1 and H 2 respectively (in order to ensure the narrow width approximation for the signal, as in the CMS analysis [8]) and up to two additional jets in the matrix element. The matching and merging between hard-scattering and parton shower has been performed via the MLM procedure 7 [38] with PYTHIA8 [39] using the shower-k T scheme. Finally, Delphes [40] is used for a simulation of the CMS detector performance which also makes use of the Fastjet [41] algorithm to cluster anti-k T [42] jets with radius R = 0.4. A crucial ingredient in this last step concerns the 13 TeV CMS b-tagging efficiencies, as well as the c-jet and light-jet mis-tag rates which are functions of the jet p T and η. We have modeled these rates using the information from [34], assuming the performance of the DeepCSV b-tagging algorithm for the same operating point as used in [8] (see appendix A for details).
Our simulated signal efficiencies at the HH and SR stages are shown in figure 1 for both the LMR and the MMR regions, together with the corresponding CMS efficiencies from [8]. Overall, we find the agreement between our validation efficiencies and those reported by the √ s = 13 TeV pp → H 1 → h SM h SM → bbbb CMS analysis to be better than 50% (except for the very low LMR masses, where the agreement is worse and the mismatch at the HH 5 The DeepCSV b-tagging medium working point used yields an average b-tagging efficiency of 68% and respective mistag probabilities for c-jets and light-jets of 12% and 1.1% [34] (see appendix A for details). 6 For both LMR and MMR selection categories, in case of multiple HH candidate combinations in an event, the combination that minimizes χ as defined in eq. (2.1) is chosen. We also note that ∆R bb depends only on the mass ratio mH 1 /m h SM [35], and as such the ratio of signal efficiencies at 4b and HH candidate stage could for the MMR category in principle be approximately extrapolated to a 2D mass plane, modulo acceptance effects that depend on the individual scalar masses. 7 We have set xqcut=qcut=mH 1 /4 GeV. stage can reach 90%). The mostly moderate mismatch at the HH selection level can be understood from our use of a fast detector simulation and our relatively limited information in the modeling of b-tagging efficiencies, as discussed in appendix A. The agreement is nevertheless very good (better than 15%) for the SR selection with m H 1 > 450 GeV, i.e. for the whole mass region of the MMR category and part of the LMR category, 8 as can be seen from figure 1.

Validation of the cross section upper limits
Besides reproducing the selection efficiencies of the CMS experimental analysis, it is crucial to check that our procedure can provide upper limits on the signal cross section consistent with those obtained by the CMS collaboration. CMS provides the inclusive background yield, dominated by QCD multi-jet processes, at the SR selection stage in the (m H 1 2 , m H 2 2 ) plane (recall eq. (2.1)) for the MMR category [8], while such information is not available for the LMR category. This background yield is independent of the value of m H 1 considered and in the region defined by eq. (2.1), approximately 2630 SM background events are found. This results in a 95% C.L. upper limit on the signal event yield of ∼ 105 events considering only the statistical uncertainty on the SM background, using a significance measure of N S / √ N S + N B (with N S and N B respectively the number of signal and SM background events). The derived limits on the inclusive pp → H 1 → h SM h SM → bbbb signal cross section are shown in the right panel of figure 2  presence of SM background systematic uncertainties, our significance measure gets modified to N S / N S + N B + u 2 B N 2 B , with u B the SM background systematic error, and we also derive the corresponding inclusive limits assuming a 3% systematic uncertainty on the SM background, illustrated in the right panel of figure 2 as a dot-dashed blue line. The 3% value chosen for the background systematics is mildly conservative for the MMR category and clearly shows the degrading of the limits due to systematic uncertainties in figure 2. This value is chosen based on a comparison of our analysis with HE-LHC projections for resonant di-Higgs production [43], where we find that a value of 2% reproduces the projected limits quoted therein (see section 4.2 for more details).
The above inclusive limits for the MMR category (we note again that for the LMR category it is not possible to extract an inclusive limit from public CMS data) are a factor ∼ 3 − 4 weaker than the CMS 95% C.L. upper limit on the signal cross section from [8], shown in the right panel of figure 2 as a solid (dashed) black line for the 95% C.L. observed (expected) limit. The reason is that the limit is not computed in an inclusive manner (that is, solely from the signal and SM background event yield after SR selection defined by eq. (2.1)); rather, it is extracted by fitting the SM QCD background distribution after the SR selection as a function of the invariant mass of the four b-jet system m 4b , and considering only SM background events within a certain width around the signal hypothesis -6 -JHEP02(2020)002 m 4b ∼ m H 1 . Using the SM background m 4b distribution after SR selection provided by the CMS analysis [8] for the MMR (provided in [8] for m H 1 > 550 GeV) and LMR categories, and defining a ±2 Γ H 1 mass window 9 around each m H 1 signal hypothesis, we obtain the corresponding fitted 95% C.L. upper limits on the pp → H 1 → h SM h SM → bbbb signal cross section, both without systematic uncertainties and again assuming a 3% systematic uncertainty on the SM background. These are respectively shown in figure 2 for the MMR (right panel) and LMR (left panel) categories as solid red lines (no systematics) and dashdot red lines (3% systematics), showing good agreement with the expected 95% C.L. upper limits reported by the CMS analysis. These results validate our extrapolation of the CMS analysis [8] to search for BSM scalars, which we do in the next section.

Searching for new scalars via
Having validated our implementation of the CMS experimental analysis, we can now proceed to extend the search to the (m H 1 , m H 2 ) mass plane. For the event generation we follow the procedure discussed in the previous section, within a 2D mass grid defined as follows: • m H 1 is varied in the range [300, 1000] GeV in steps of 50 GeV.
• For each m H 1 value, m H 2 is varied in the range [65 GeV, m H 1 /2], taking ten equally spaced values.
The various parameters defining the SR selection in eq. (2.1) for LMR and MMR categories need to be modified accordingly. For the MMR category, we maintain R = 20 GeV and set C = m H 2 , while for the LMR category we also keep R = 20 GeV and set instead C = (120/125)×m H 2 . In addition, for the LMR the HH candidate selection on the b-jet pairs needs to be modified to |m bb − (120/125) m H 2 | < Max[20 GeV, (40/125) m H 2 ]. These modifications match the CMS analysis HH and SR selection criteria for m H 2 = 125 GeV, while prodiving a natural generalization of those for other values of m H 2 . 10 For the estimate of the SM QCD multi-jet background, we first extrapolate the SM background event yield in the SR as a function of m H 2 to the region m H 2 > 300 GeV using a smoothly falling exponential fit. 11 Then, we adopt the following procedure: • In the MMR category with m H 1 > 550 GeV, the fitted SM background event yield is computed as described in the previous section, with an additional overall rescaling of the (MMR) SM background m 4b distribution. The factor is determined by the ratio of the inclusive background yield in the SR with C = m H 2 over the SR background yield for C = 125 GeV. This procedure assumes that the SM background m 4b shape 9 Specifically, we adopt [mH interpolate linearly between them. 10 We stress that our selection is parameter-space dependent, and it may be possible to significantly optimise the sensitivity of the search by exploiting Machine Learning algorithms as those used in ref. [44] for generic BSM resonant searches. 11 We observe that such a fit provides a very good description of the measured SM background yield in the SR for 125 GeV < mH 2 < 300 GeV.

JHEP02(2020)002
remains approximately unchanged and only its overall normalization varies when the SR selection from eq. (2.1) is redefined by setting C = m H 2 .
• In the LMR category with m H 1 < 550 GeV, 12 we follow the same strategy as for the MMR category above, performing the aforementioned overall rescaling of the SM background m 4b distribution (now for the LMR category). However, since no inclusive SM background yield after SR selection is provided by the CMS analysis for the LMR category, we use the same rescaling factor (as a function of m H 2 ) as for the MMR category.
• For the MMR category with m H 1 < 550 GeV, it is not possible to apply the above strategy, since the SM background m 4b distribution for the MMR category is not provided by the CMS analysis [8] in this region (and thus the fitted 95% C.L. upper limits on the signal cross section cannot be derived). Instead, we use the total SM QCD background yield in the SR (redefined by C = m H 2 ) to obtain the inclusive 95% C.L. upper limits on the signal cross section.
These procedures for the estimate of the SM background are nevertheless expected to fail both for m H 1 < 300 GeV and m H 1 → 2m H 2 , since in these regions the CMS measured multi-jet data do not follow a smoothly falling distribution, but rather display a kinematic feature near the threshold region, driven by the kinematic selection of the analysis (mainly trigger effects) [8]. This leads to a feature that depends on m H 2 , and peaks around 300 GeV in the CMS analysis. We avoid being near these regions of the (m H 1 , m H 2 ) plane by imposing m H 1 > 300 GeV and m H 1 > 2m H 2 + 25 GeV in our analysis.
The 95% C.L. upper limits on the pp → H 1 → H 2 H 2 → bbbb cross section, in the narrow width approximation, are shown in figure 3 for the MMR category and in figure 4 for the LMR category as color coded heat maps. For the MMR category, figure 3 shows both the fitted limit for m H 1 > 550 GeV assuming a 2% SM background systematic uncertainty and the inclusive limit for m H 1 < 550 GeV with a SM background 0.1% systematic uncertainty. The same 0.1% background systematic uncertainty is assumed for the LMR category. These choices for the SM background systematic uncertainty are motivated in section 4.2, and are conservative given the information from the CMS analysis [8]. The SM background systematic uncertainties are driven by the background modeling, resulting in a smaller error at lower invariant masses given the larger statistics in that region, which explains the difference in systematics between MMR and LMR categories. We also show the 95% C.L. signal cross section upper limits with no SM background systematic uncertainties in appendix B. By comparing the results from figures 3 and 4 with figure 16 (top-left panel) in appendix B, we see that the effect of background systematics is not very important for current signal upper limits, which are at present statistically dominated.
Comparing the inclusive signal 95% C.L. upper limits for the MMR category (m H 1 < 550 GeV) to those for the LMR category from figure 4, we see the latter are much stronger except for the small region m H 2 80 GeV, m H 1 400 GeV (see figure 4). We thus omit 12 In the rest of the paper, we consider the boundary between LMR and MMR categories at mH 1 = 550 GeV.  Figure 3. 95% C.L. expected upper limit on the pp → H 1 → H 2 H 2 → bbbb signal cross section (in fb) in the (m H2 , m H1 ) plane for the MMR category, for m H1 > 550 GeV (fitted limit, assuming a 2% systematic uncertainty) and for m H1 < 550 GeV (inclusive limit, assuming a 0.1% systematic uncertainty), extending the CMS analysis [8]. The various contours correspond to the 95% C.L. upper limits on the sensitivity κ 2 × BR from eq. (3.1), see text for details.
from now on the use of the MMR inclusive results for m H 1 < 550 GeV and use the LMR and MMR fitted 95% C.L. signal cross section upper limits respectively for m H 1 lighter and heavier than 550 GeV. We can then parametrize our BSM cross section as withσ H 1 the inclusive production cross section of a SM-like Higgs boson with mass m H 1 , κ 2 an effective rescaling factor with respect to the SM-like Higgs boson cross section and BR = BR(H 1 → H 2 H 2 → bbbb). Through this parametrization we can translate the derived upper limits on the cross sections into limits on . These limits, shown as isocontours in figures 3 and 4 keeping in mind that they are valid in the narrow width approximation, serve as a reference point for understanding the possible impact of our search in specific BSM models, which are discussed in detail in section 5. From figure 4 we see that in the LMR category, for a fixed value of m H 1 the sensitivity of the search increases with the mass m H 2 . On the other hand, for the MMR category and a fixed m H 1 , increasing m H 2 from 65 GeV results in an increase in sensitivity up to an optimal value of m H 2 ∼ 120 − 180 GeV (depending on the value of m H 1 ), above which the sensitivity of the search drops and the fitted 95% C.L. limit on the signal cross section  quickly becomes very large, as shown in figure 3. This behaviour of the sensitivity for the LMR and MMR categories can be understood from the interplay between the signal acceptance and the SM multi-jet background yield in the SR. The SM background yield in the SR decreases rapidly as m H 2 increases, which explains the behaviour observed for the LMR category, as well as the initial growth in sensitivity for m H 2 > 65 GeV in the MMR category. For the MMR category, the SR acceptance decreases for the BSM signal as m H 2 increases for a fixed m H 1 and eventually overcomes the decrease in the SR background yield, and the sensitivity drops again. This however does not occur for the LMR category, which retains sensitivity to the 2 m H 2 → m H 1 region. Overall, it is interesting to note that values of σ(pp → H 1 → H 2 H 2 → bbbb) in the ballpark of several tens of fb can be accessed for the MMR category in the optimal region m H 2 ∼ 100 − 150 GeV, yielding a sensitivity to κ 2 × BR 0.02 − 0.1 depending on the value of m H 1 . At the same time, the sensitivity to κ 2 × BR for the LMR category reaches values as low as 7 × 10 −3 , with most of the LMR parameter space being constrained to κ 2 ×BR < 0.05 at 95% C.L. This suggests that our proposed search may indeed be sensitive to new scalars in realistic BSM scenarios. We discuss this in more detail in section 5.  the targeted total integrated luminosity. The different future setups that we consider 13 are shown in table 1. We furthermore study the impact of systematic uncertainties, and briefly comment on the role of trigger and b-tagging efficiencies, which are critical for this final state.
4.1 General procedure for the extrapolation to higher √ s The main challenge in computing the reach for hadron colliders with higher centre of mass energy is the scaling of the multi-jet background. This is estimated in a data-driven way by the CMS and ATLAS experimental collaborations and cannot be reliably simulated within our framework. Both the SM background and the signal increase for higher collider energies, while the parton shower is expected to provide a larger number of significantly harder jets, thus possibly changing the kinematic features of the events. Here we perform our analysis by naively assuming that both the signal and the SM background cross sections scale roughly by the same amount, given by the ratio of inclusive gg → H 1 production cross sections at √ s = 13 TeV and √ s = X TeV, namely This assumption is motivated by the fact that the overall partonic centre of mass energy of both the signal and the SM background after the full event selection will be peaked around √ŝ ∼ m H 1 and is valid to the extent that the gg luminosity dominates the multijet rate. To compute the rescaling factor r X we use SusHi v1.7.0 [47,48], which gives NNLO accuracy for the production of a SM-like Higgs boson in the infinite top mass limit [49][50][51][52][53]. These cross sections are reported in figure 5 for the various collider scenarios of table 1. We further assume that the acceptance and selection efficiencies of the proposed search remain approximately constant for the different collider scenarios considered. This is expected to be a good approximation for the case of the HL-LHC (modulo improvements in trigger and b-tagging efficiencies, which we discuss at the end of section 4.2), while it will not be very accurate for the HE-LHC and FCC-hh. Our results should be interpreted as a conservative first estimate of the pp → H 1 → H 2 H 2 → bbbb search channel sensitivity at the HL-LHC and future hadron colliders, bearing in mind that future experimental analyses may improve upon this estimate.
In general, the expected improvement in sensitivity on the κ 2 × BR factor of eq. (3.1) for a future collider with respect to the current √ s = 13 TeV limits from section 3 may be 13 We also note that, as of now, ATLAS and CMS have recorded ∼ 150 fb −1 each, currently under analysis. simply written as is given in the absence of systematic uncertainties simply by As an example, for m H 1 = 1 TeV the ratio σ 95%C.L.
is given, in the absence of systematic uncertainties, by 0.12, 0.13 and 0.36 respectively for HL-LHC, HE-LHC and FCC-hh, which would then result sensitivity improvements of in respective I = 10, 50 and 300 respectively. Nevertheless, the impact of systematic uncertainties on the ratio σ 95%C.L.
13 may be important, and we discuss this in more detail in section 4.2.

Validation & Extrapolation to HL-LHC
Before presenting our results for the HL-LHC extrapolations, we validate our procedure by comparing our findings for the case m H 2 = 125 GeV with the official CMS projections for the HL-LHC in the h SM h SM → bbbb channel [43]. The CMS collaboration reports projected 95% C.L. cross section upper limits of 46, 7.3 and 4.4 fb respectively for m H 1 = 300 GeV, 700 GeV and 1 TeV respectively, whereas we obtain 32.1, 2.4 and 1.4 fb with our extrapolation procedure and without systematic uncertainties. We find that we approximately  reproduce 14 the CMS projected 95% C.L. cross section upper limits by assuming a 2% systematic uncertainty for the MMR benchmarks, m H 1 = 700 GeV, 1 TeV, and a 0.1% systematic uncertainty for the LMR benchmark m H 1 = 300 GeV. Hence in our HL-LHC extrapolations (as well as for the HE-LHC and FCC-hh extrapolations from section 4.3) we will consider a flat systematic error of 2% for the MMR category and of 0.1% for the LMR category (we also show the HL-LHC, HE-LHC and FCC-hh extrapolations in the absence of systematic uncertainties in appendix B), which are also the values adopted in section 3 for the current upper limits. We stress again that these uncertainties are dominated by the data-driven SM background modeling. 15 The projected 95% C.L upper limits on the signal cross section for HL-LHC are shown in figure 6 for both the MMR (m H 1 > 550 GeV) and LMR categories, together with the κ 2 ×BR sensitivity defined in eq. (3.1). By comparing them with the results from section 3, we observe a factor I 5 − 10 improvement in the sensitivity for the HL-LHC with respect to the present reach in the κ 2 × BR factor for the LMR category, becoming larger as m H 2 14 We however stress that we reproduce the CMS projected 95% C.L. cross section upper limits without It is worth studying here in more detail the impact of systematic errors on the projected HL-LHC 95% C.L. signal upper limits, since in this case (as opposed to that of section 3) the upper limits start to become systematics dominated. For the discussion of systematic errors we follow ref. [54]. Considering that we are dealing mostly with b-jets (with a reported 2-6% overall systematic uncertainty), and taking into account additional O(1%) sources of systematics (e.g. integrated luminosity, jet energy scale. . . ) treated as uncorrelated, we may expect a 1-10% range for the overall systematic uncertainty. We nevertheless stress that due to the large statistics for the data-driven SM background, background systematic uncertainties are expected to be small. In figure 7 we show the variation of the projected HL-LHC 95% C.L. upper limit on the signal cross section as a function of the SM background systematic uncertainty, for fixed m H 2 = 125 GeV (left panel) and fixed m H 1 = 1 TeV (right panel). We observe how increasing the SM background systematic error leads to a saturation of the 95% C.L. upper limit on the signal cross section when systematic uncertainties dominate over statistical ones. In the left panel of figure 7 we show that for low m H 2 , where the multi-jet SM background is large, already a few % systematic error on the SM background leads to an important increase of the 95% C.L. signal upper limit, thus stressing the key importance of having systematic effects under control for this study. At the same time, increasing the mass m H 2 leads to a significant decrease of the SM background yield in the SR and thus to a smaller impact of the corresponding background systematics, as can be seen from the right panel of figure 7. This panel also shows that for the MMR category (in this case, for a benchmark m H 1 = 1 TeV), the interplay between the decrease of SM background yield and the decrease of signal efficiency in the SR as m H 2 increases leads to a minimum for the 95% C.L. signal upper limit (as a function of m H 2 ) for m H 2 in the range m H 2 ∼ 140-200 GeV (recall the discussion at the end of section 3).  Table 2. Value of κ 2 ×BR for different m H2 , m H1 (in GeV) benchmarks, for our proposed search with current ( √ s = 13 TeV, 35.9 fb −1 ) LHC data (assuming respectively a 2% and 0.1% systematic uncertainty for the MMR and LMR benchmarks), as well as the extrapolations to HL-LHC, HE-LHC and FCC-hh (assuming the same uncertainties). We also give the sensitivity improvement I for HL-LHC, HE-LHC and FCC-hh.
We close this section with a few remarks on the impact of b-tagging and trigger efficiency. We first note that accounting for potential improvements on b-tagging for the HL-LHC and future colliders in our projections is rather difficult due to: i) the b-tagging efficiency being a phase-space dependent (p T , η) quantity. This is already an important issue for the current analysis, as detailed in appendix A; ii) our lack of knowledge of the b-jet truth content of the SM multi-jet background; iii) noting that improvements would come from a deep-learning algorithm which will most likely not yield a flat rescaling in phase-space; iv) noting that improvements strongly depend on both the detector capabilities and their performance, as well as on the pileup conditions, all of which are not fully known for future colliders like HE-LHC and FCC-hh. Nevertheless, we can obtain a rough idea of potential improvements quoted in the recent literature. On one hand, the effect of changes in b-tagging efficiency on the overall signal strength uncertainty has been evaluated by the CMS collaboration, showing that an improvement of 10% in the b-tagging efficiency leads to a relative improvement in the signal strength uncertainty of up to 6% [55]. On the other hand, the inclusion of timing information (which helps to reduce the number of spurious reconstructed secondary vertices by ∼ 30%) provides an increase in the b-tagging efficiency of about 4-6% depending on the pseudorapidity, evaluated for the same mis-tag rate [56]. Finally, while the challenging data-taking conditions at the HL-LHC could worsen the b-tagging efficiency, the new inner tracker detector as well as novel reconstruction techniques could provide a sizeable improvement. For example, it has been estimated that the upgrades of the inner tracker would lead to an 8% improvement in efficiency [57].

Extrapolation to HE-LHC and FCC-hh
Using the results from section 4.1, here we provide an extrapolation of the 95% C.L. σ(pp → H 1 → H 2 H 2 → bbbb) upper limits to a √ s = 27 TeV HE-LHC collider and a √ s = 100 TeV FCC-hh collider, see table 1. As for the extrapolation to HL-LHC performed in the previous section, we assume here a background systematic error of 0.1% and 2% for the LMR and MMR categories, respectively (the results for HE-LHC and FCC-hh without background systematic uncertainties are shown in appendix B). Our results for HE-LHC and FCC-hh are shown respectively in figure 8 and figure 9. By comparing these with the HL-LHC results from figure 6, we can readily see that SM background systematic uncertainties significantly hinder the potential improvement in sensitivity of HE-LHC and FCC-hh for the MMR category. The 2% SM background systematics result in the signal sensitivity being dominantly driven by this error, such that the much larger integrated luminosities of

as a probe of extended Higgs sectors
We now analyze the LHC sensitivity of the proposed search pp → H 1 → H 2 H 2 → bbbb in the context of specific extensions of the SM. Our aim here is two-fold. First, we will assess the reach of this search within the parameter space of several well-studied BSM models. Second, we will compare the projected sensitivity of the search to other analyses for BSM scalars, studying their complementarity and identifying where pp → H 1 → H 2 H 2 → bbbb provides the leading probe of the existence of these new scalars. For our analysis we consider three benchmark models: • A simplified model with two real scalar singlets added to the SM in section 5.1, which allows for a direct mapping of the sensitivities derived in sections 3 and 4 to many BSM scenarios.
• A Type-I 2HDM scenario in section 5.2.
• A 2HDM (of Type-I) with the addition of a real (pseudo)scalar singlet, which captures the features of more complicated scalar sectors as e.g. the one of the NMSSM in section 5.3.
We stress that the use of our CMS analysis validation (recall sections 2 and 3) to derive limits for these extensions of the SM assumes the narrow width approximation for the signal, which we have nevertheless verified to hold in a large portion of the parameter space of the models.

Two singlet scalar extension of the SM
One of the simplest possibilities is to consider that both the H 1 and H 2 states come from singlet scalar fields S 1,2 (see e.g. [30] for a recent phenomenological analysis of this scenario). This is a simplified framework to which more complicated models could be mapped. The most general scalar potential for the SM Higgs H and two singlet scalars S 1 , S 2 has the following form with 2 ≤ a + b + 2c ≤ 4 (we disregard tadpole terms for S 1,2 ). While the most general potential from eq. (5.1) has 17 free parameters (once the SM Higgs vev v and the Higgs mass m h are fixed), in practice most of them are phenomenologically unimportant and may be safely ignored in the present analysis. In particular, considering the process 16 pp → S 1 → S 2 S 2 → bbbb we only care about the ggS 1 , S 1 S 2 S 2 and S 2 bb interactions. Regarding the former (effective) coupling between S 1 and the gluons, this would naively come from the first singlet mixing with the SM Higgs, the mixing given by sin α 1 . Yet, there are other possibilities, e.g. the additional presence of vector like quarks at/above the TeV scale which couple to S 1 (see e.g. [58]). These latter interactions allow for having 16 With some abuse of notation, we label for the rest of this section the singlet-like scalar mass eigenstates as S1 and S2 (even if they do not correspond exactly to the singlet scalar states from (5.1) due to singletdoublet mixing). L. sensitivity to the pp → S 1 → S 2 S 2 cross section (in fb) for m S1 < 550 GeV (LMR category) and m S1 > 550 GeV (MMR category). Also shown are contours of 95% C.L. current upper limits (left panel) and HL-LHC projections (right panel) on sin 2 α 1 assuming BR(S 1 → S 2 S 2 ) 1. The narrow width approximation for the signal has been assumed. the ggS 1 coupling C 1 g as a free parameter in our setup, which can then be traded for the production cross section σ(pp → S 1 ). Nevertheless, we also discuss below the interplay of our analysis with other LHC searches for S 1 when its production at the LHC comes purely from the singlet-doublet mixing sin α 1 .
The interaction S 1 S 2 S 2 is also a free parameter, dominantly controlled by the coupling λ 1,2,0 in (5.1) for small singlet-doublet mixing. We note that in this small mixing limit it is possible for the branching fraction BR(S 1 → S 2 S 2 ) to approach unity even with moderate values of λ 1,2,0 , strongly suppressing the sensitivity of other search channels for S 1 in this case. Finally, the coupling of S 2 to the SM fermions is generated via mixing 17 between the second singlet and the SM Higgs, the mixing given by sin α 2 . This results in branching fractions of S 2 to the SM states equal to those of a SM Higgs boson with mass m S 2 , independently of the value of sin α 2 (for m S 2 < 250 GeV; above this, the possible decay of S 2 into two 125 GeV Higgs bosons would modify this behaviour). Incidentally, this has 17 In certain scenarios it would also be possible to generate the interaction between S2 and the SM fermions via the dimension 5 effective operator (cqyq/Λ) S2QLHqR (see e.g. [59,60]). We however do not consider this possibility here.

JHEP02(2020)002
the consequence that in this model, our proposed search would be most suited for singlet scalar masses m S 2 150 GeV, since for m S 2 150 GeV it ceases to be efficient due to the sharp drop in BR( should be most sensitive to the existence of S 1 and S 2 in this region. To summarize the above discussion, the present model has as relevant free parameters C 1 g , m S 1 , m S 2 , BR(S 1 → S 2 S 2 ) and the mixing sin α 2 . The parameter C 1 g can be traded by either the mixing sin α 1 or directly the production cross section σ(pp → S 1 ). At the same time, the value of sin α 2 is only important regarding the complementarity with other collider probes of S 2 , since BR(S 2 → bb) is solely determined by m S 2 for m S 2 < 250 GeV. Using the results from sections 3 and 4 we show in figure 10 the current (left plot) and HL-LHC (right plot) 95% C.L. exclusion sensitivity to σ(pp → S 1 ) × BR(S 1 → S 2 S 2 ) in the mass plane (m S 2 , m S 1 ). In both cases we have assumed a 2% (0.1%) SM background systematic uncertainty for m S 1 > 550 GeV (m S 1 < 550 GeV), as in sections 3 and 4.2. Then, assuming the gluon fusion production of S 1 to come exclusively from the Higgssinglet mixing and setting also BR(S 1 → S 2 S 2 ) 1, the current/projected 95% C.L. upper limit on σ(pp → S 1 → S 2 S 2 ) from figure 10 can be rephrased as an upper limit on sin 2 α 1 , 18 shown as coloured contours in figure 10. We stress that these assume implicitly the narrow width approximation, which breaks down for large values of BR(S 1 → S 2 S 2 ) and/or sizable values of sin 2 α 1 . More concretely, assuming that S 1 → S 2 S 2 and S 1 → V V are the dominant decays of S 1 , the total width of S 1 can be approximately written as with Γ SM (m S 1 ) the total width of a SM Higgs of mass m S 1 . In order for the narrow width approximation to hold, we require Γ S 1 to be smaller than the value 2 Γ H 1 defining the mass window of the CMS search (recall footnote 9 in section 2.2). When the production of S 1 at the LHC is due exclusively to the Higgs-singlet mixing sin α 1 , it is possible to explore the interplay between the pp → S 1 → S 2 S 2 → bbbb search analyzed in this work, other direct searches for S 1 (and S 2 ) and LHC measurements of the Higgs signal strengths. The latter yield the present limit sin 2 α 1 < 0.073 at 95% C.L. (see [30]) under the assumption |sin α 2 | |sin α 1 |. In figure 11 we fix sin 2 α 1 = 0.07 and show the present sensitivity of our proposed analysis together with the sensitivity of direct BSM Higgs searches in ZZ final states 19 from the latest ATLAS analysis [61] in the plane (m S 1 , BR(S 1 → S 2 S 2 )) for m S 2 = 80, 100, 120, 140 GeV. 20 We also show the region where the narrow width approximation for S 1 ceases to be valid as obtained from eq. (5.2), and one expects a degrading of the various search limits. Properly quantifying this effect is however beyond the scope of this work. Interestingly, we see that our analysis nicely complements existing searches for BSM scalars, providing a new avenue to probe the singlet-like scalar S 1 . 18 We note that sin α1 = κ as defined in eq. (3.1). 19 Other BSM Higgs searches, e.g. those in W W , τ τ or γγ final states, are significantly less sensitive for the model considered here. 20 We assume for simplicity BR(S1 → hSMhSM) = 0 and BR(S1 → S2hSM) = 0. If these branching fractions are sizable, they will weaken the constraints from figure 11. MMR LMR Present 95% C.L. Limits, sin 2 α 1 = 0.07 S1 → S2S2 (m S2 = 80 GeV) S1 → S2S2 (m S2 = 100 GeV) S1 → S2S2 (m S2 = 120 GeV) S1 → S2S2 (m S2 = 140 GeV) S1 → ZZ Figure 11. 95% C.L. limits on the branching fraction BR(S 1 → S 2 S 2 ) as a function of m S1 from S 1 → ZZ searches by ATLAS (observed limit) with 36.1 fb −1 [61] (yellow region) and from our search pp → S 1 → S 2 S 2 → bbbb (expected limit) assuming m S2 = 80 GeV (solid red), m S2 = 100 GeV (dashed red), m S2 = 120 GeV (solid blue) and m S2 = 140 GeV (dashed blue). The Higgs-singlet mixing has been set to sin 2 α 1 = 0.07 (satisfying LHC measurements of Higgs signal strengths), and we assume |sin In the grey region the narrow width approximation for the signal as derived from (5.2) no longer holds and the limits shown are expected to degrade.

2HDM
The 2HDM represents the simplest scenario where the two BSM states H 1 and H 2 are contained in one field, a second Higgs doublet. As opposed to the model discussed in section 5.1, for the 2HDM the interactions of H 1,2 with the SM fermions and gauge bosons are not governed by their mixing with the SM Higgs. The scalar potential for a theory with two Higgs doublets Φ 1,2 (with a softly-broken Z 2 -symmetry and no CP-violation) is given by with all scalar potential parameters being real. The breaking of EW symmetry is shared between the two doublets, whose vacuum expectation values (vevs) are given by v 1,2 (with v 2 1 + v 2 2 = v = 246 GeV, v 2 /v 1 ≡ tan β). In addition to the 125 GeV Higgs state h, the 2HDM scalar sector contains another neutral CP-even scalar H, a neutral CP-odd scalar A and a charged scalar H ± . In the following we identify H and A with our neutral BSM states H 1 and H 2 , respectively.
The couplings of the 2HDM scalar states to SM gauge bosons and fermions are controlled by tan β and by a mixing angle α in the CP-even neutral sector. The limit of a JHEP02(2020)002 SM-like 125 GeV Higgs h = h SM (the so-called "alignment" limit of the 2HDM [62]) corresponds to cos (β − α) = 0. In addition, there are various possible choices for the couplings of the two doublets Φ 1,2 to fermions (see [63] for a detailed discussion of fermion couplings in 2HDMs), and in this work we choose for definiteness a so-called Type-I 2HDM.
Given a set of values for m H , m A , m H ± , cos (β − α) and tan β, theoretical constraints dictate the allowed range for µ 2 in (5.3). These constraints are the boundedness from below of the 2HDM scalar potential as well as the stability of the EW minimum, and the requirements of unitarity [64] and perturbativity on the quartic couplings λ i (see e.g. [22,23] for more details). It is possible that certain choices for m H , m A , m H ± , cos (β − α), tan β yield no allowed range for µ 2 , making these choices not physically viable. When a viable range for µ 2 exists, its specific value has an impact on the 2HDM scalar self-couplings (as shown below).
In the following we consider 21 m H ± m H > 2 m A , such that the decay H → AA is open and the analysis from sections 2-4 may be applied to the 2HDM for H 1 ≡ H and H 2 ≡ A. The coupling g HAA is given by [22,23] with M 2 ≡ µ 2 /(s β c β ) and we use the notation c φ ≡ cos φ, s φ ≡ sin φ. We note that in the alignment limit c β−α = 0, g HAA vanishes for tan β = 1 and/or m 2 H = M 2 . For m H m A the decay H → ZA will be present together with H → AA, such that both decay modes may compete to be the dominant one (for c β−α = 0 the decay modes H → W + W − , H → ZZ and H → hh could also be important). The process pp → H → ZA (Z → , A → bb) has been searched for by CMS at  [47,48] to obtain their LHC production cross section. After imposing 2HDM theoretical constraints and ensuring compatibility with the ATLAS/CMS H → ZA experimental searches, we obtain the maximum possible cross section 23 for the process pp → H → AA → bbbb in the 21 Measurements of EW precision observables require H ± to be close in mass to either A or H, to avoid a large breaking of custodial symmetry. At the same time, bounds from flavour physics constrain m H ± to be above a certain value at 95% C.L. (which depends on the 2HDM Type), see e.g. [65]. These motivate our choice m H ± ∼ mH . 22 The 13 TeV search by ATLAS [67] interprets its results in terms of an A → ZH decay, but these are equally applicable to H → ZA. 23 The procedure we follow is, for a given set of values mH , m H ± , mA, c β−α and tan β within our scan, to maximize the BR(H → AA) within the range of µ 2 allowed by theoretical constraints (which amounts to choosing the value of M 2 for which gHAA in (5.4) is maximal in that range), and subsequently discarding the points which are ruled out by the H → ZA ATLAS/CMS searches. . Also included are the latest constraints from √ s = 13 TeV LHC pp → H/A → τ τ CMS searches [71] and ATLAS di-boson (pp → H → ZZ) searches [61] (the latter only relevant for c β−α = 0.1). In the white region of figure 12, it is not possible to simultaneously satisfy the various theoretical and experimental constraints, while for the coloured region there exist points in our scan that satisfy all constraints, from which we extract a maximum allowed value for σ(pp → H → AA → bbbb) as a function of m A and m H .
As figure 12 highlights, the combination of 2HDM theoretical constraints and current experimental limits from BSM Higgs searches rule out a sizable fraction of the (m H , m A ) parameter space for m H > 2 m A . For c β−α = 0.1 the search H → ZZ yields stringent constraints on the 2HDM parameter space, which translate into allowed values σ(pp → H → AA → bbbb) 10 fb, too low to be probed by our proposed search. In contrast, in the alignment limit c β−α = 0, the regions that survive the combination of 2HDM theoretical constraints and limits from ATLAS/CMS 13 TeV H → ZA searches yield a large cross section for the process pp → H → AA → bbbb, particularly for m A < 130 GeV. Figure 12 shows (red contours) the (m H , m A ) region for which the current sensitivity of our proposed search allows to probe currently unconstrained 2HDM parameter space, making this search highly complementary to other BSM scalar searches in the 2HDM.

2HDM + singlet scalar/pseudoscalar
We now consider the addition of a real scalar/pseudoscalar singlet field S to the above 2HDM (for a detailed discussion of the 2HDM + complex singlet, see [24,29]), assuming for simplicity CP conservation in the scalar potential. This scenario captures the features JHEP02(2020)002 of richer, more complicated scalar sectors than those analyzed in sections 5.1 and 5.2. At the same time, it is well-motivated as a simplified version of the NMSSM and as a portal to a dark matter sector 24 [73][74][75][76]. As emphasized in [24,29], in this class of extended Higgs sectors the Higgs-to-Higgs cascade decays we discuss in this work are ubiquitous.
Besides eq. (5.3), the scalar potential for the 2HDM + real singlet field S contains the following terms We note that if CP is conserved, the last three terms in eq. (5.5) are absent for a real pseudoscalar field S, and µ S has to be purely imaginary. In the following, we discuss separately the phenomenology for scalar and pseudoscalar S, highlighting the differences between them and analyzing the sensitivity of our proposed search in each case.
Pseudoscalar S: for µ S = 0 in eq. (5.5), the pseudoscalar S will mix with the 2HDM CP-odd state after EW symmetry breaking, yielding two CP-odd mass eigenstates, which we label A, a (with m A > m a ). The lighter mass eigenstate a is considered here to be mostly singlet-like. In addition to these, the model contains two CP-even scalars H and h (this last one identified with the 125 GeV Higgs boson) and a charged scalar H ± , as in the 2HDM scenario studied in the previous section. Then, assuming for simplicity the 2HDM alignment limit c β−α = 0, we can analyze the sensitivity of the various LHC probes of the states H and a, as well as the interplay among them. We then perform a comparison of direct searches of a via pp → a → τ τ and pp → a → γγ with searches for pp → H → Za (a → bb) and our proposed search pp → H → aa → bbbb, as a function of m H , m a , tan β, the singlet-doublet mixing s θ and the branching fraction BR(H → aa). The various partial decay widths of a are given in [74], together with the partial decay widths of H in the alignment limit c β−α = 0. The LHC production cross sections for H and a are computed at NNLO in QCD with SusHi [47,48]. Then, assuming the other 2HDM states (H ± and A) are sufficiently heavy to not play a phenomenological role in the following, we show in figure 13 [71]). We find pp → a → γγ searches (e.g. [77]) are currently not sensitive to the parameter space of the model.
As highlighted in figure 13, direct searches for a barely have sensitivity to this model, with pp → a → τ τ only constraining tan β 1, and pp → a → γγ being even less sensitive 24 see also [72] for a study on the connection between dark matter and di-Higgs signatures. at present. In contrast, searches for cascade scalar decays probe a sizable region of the parameter space, with a strong interplay between H → aa and H → Za searches: for large singlet-doublet mixing (e.g. left panel of figure 13) H → Za decays typically drive the sensitivity for m a > 130 GeV, with H → aa providing the strongest sensitivity for m a < 130 GeV (we note the ATLAS search [67] does not go below 130 GeV for the mass of the lighter BSM scalar); however, as the singlet-doublet mixing diminishes (right panel of figure 13) the H → Za sensitivity weakens significantly and the H → aa decay mode (which does not necessarily vanish in the limit s θ → 0 [74]) may become the leading probe of the parameter space of the model.
The independence of the H → aa decay on the zero mixing limit s θ → 0 is a strong point of this particular search mode. This is further emphasised by the fact that EW precision constraints from the ρ parameter directly constrain s θ as a function of the scalar sector masses [74]. For the 2HDM, the well-known BSM contributions to the ρ parameter (see e.g. [78]) vanish for either m H = m ± H or m A = m ± H in the alignment limit. For non-zero s θ the 2HDM pseudoscalar contribution is shared between the two CP-odd mass eigenstates, which yields an additional source of custodial symmetry breaking. In the decoupling regime for A and H ± considered here, the ρ parameter drives the model towards either small singlet-doublet mixing or tuned regions of parameter space. Finally, let us stress that if the states A and H ± are not decoupled (contrary to what has been considered so far), other BSM scalar cascade decays not considered in this work such as A → ha and H ± → W ± a could allow to probe the 2HDM + singlet pseudoscalar scenario.

JHEP02(2020)002
Scalar S: in this case the last three terms of (5.5) may be present for a CP-conserving potential. After EW symmetry breaking, the singlet state S (we consider for simplicity that the singlet field does not get a vev; the discussion when it does is however analogous) mixes with the CP-even states H and h from the 2HDM, yielding three CP-even mass eigenstates, one of which is the 125 GeV Higgs boson. The other two states we label H 1 and H 2 , with m H 1 > m H 2 , and we consider H 2 to be the singlet-like state. Note that, contrary to the pseudoscalar case, the ρ parameter does not necessarily imply constraints on the singlet-doublet mixing (e.g. in the 2HDM alignment limit these can be avoided by having m A m ± H and the cancellation is no longer spoiled by a−A mixing). For simplicity, we focus in the following on the limit where the 125 GeV Higgs boson is not mixed with the other 2HDM CP-even state (corresponding to the alignment limit of the 2HDM, see section 5.2), and may only have a small singlet admixture as a departure from SM-like properties. 25 In addition, we consider the mixing between the singlet S and the heavy 2HDM state H (see section 5.2) to be sizable, such that the decays of H 2 into SM particles are controlled by this mixing, which we parametrise by sin θ.
Again, we consider the case where the 2HDM states A and H ± are decoupled. 26 The leading decays of the state H 1 are H 1 → H 2 H 2 , decays into SM fermions and the decay into two 125 GeV Higgs bosons H 1 → h 125 h 125 . The latter vanishes if h 125 is the SM Higgs boson, but is non-zero if h 125 has a small singlet admixture. At the same time, the LHC production cross section for H 1 is suppressed by cos 2 θ compared to the production of the 2HDM state H in the alignment limit c β−α = 0. In figure 14 we demonstrate the interplay between resonant di-Higgs searches pp → H 1 → h 125 h 125 and our proposed search pp → H 1 → H 2 H 2 → bbbb. Direct searches for the BSM state H 2 (e.g. pp → H 2 → τ τ and pp → H 2 → γγ) are currently only sensitive to the tan β < 1 and large mixing sin θ → 1 region of the 2HDM + scalar singlet scenario (similarly to what we already found above for the 2HDM + pseudoscalar singlet scenario). Figure 14 then shows, for fixed sin θ = 0.7 and BR(H 1 → h 125 h 125 ) = BR(H 1 → H 2 H 2 ) = 0.3 the 95% C.L. exclusion sensitivity in the (m H 1 , tan β) plane from present ATLAS and CMS resonant di-Higgs searches in the bbbb [8,9], bbγγ [17,18] and bbτ τ [13,14]   states. A rather precise estimate of the LHC sensitivity of such a search is possible given its similarity with CMS and ATLAS resonant di-Higgs searches, which we have used to validate our analysis, specifically choosing for this purpose the latest √ s = 13 TeV CMS resonant di-Higgs search in the bbbb final state [8]. The recasting procedure obtained here has the advantage of being model-independent, as it relies solely on the masses of the two BSM scalar particles. With present data from [8], the pp → H 1 → H 2 H 2 → bbbb search yields sensitivity to production cross sections times branching fractions ranging from the picobarn to tens of femtobarns depending on the BSM scalar masses, showing the power of this simple generalisation of an existing LHC search. We also stress that a dedicated experimental search is likely to yield appreciable improvements in sensitivity with respect to the one obtained in this work.
We have briefly discussed the impact of several experimental features that could affect the analysis, including b-tagging and the role of systematic errors. We have also devised a simplified procedure to obtain an estimate of the sensitivity for future collider machines, appropriately scaling the present LHC results to other center-of-mass energies and total integrated luminosities. We have analyzed the specific examples of the High Luminosity -27 -JHEP02(2020)002 LHC (HL-LHC), the High Energy LHC (HE-LHC) and of the Future Circular Collider in its proton-proton incarnation (FCC-hh).
We have then applied our analysis to three specific scalar extensions of the SM: i) a Higgs sector with two additional singlet scalars; ii) a 2HDM scenario; iii) the 2HDM plus a singlet scalar/pseudoscalar, which is currently a "de-facto" benchmark for dark matter searches at the LHC [76]. Comparing the reach of our proposed study to present LHC searches constraining these models, we explicitly show the parameter space regions that our search renders accessible, stressing its complementarity to existing searches.
Our study shows promising prospects for this yet unexplored probe of heavy Higgs bosons, and highlights explicitly how extending the coverage of current LHC searches for BSM scalars can yield new avenues to probe non-minimal Higgs sectors. Finally, we note that our study represents the first, minimal step in probing scenarios with Higgs-to-Higgs decay topology in which all states come from the BSM sector. As the paradigm for extended scalar sectors evolves into increasingly non-minimal territory, it is essential that the experimental programme continues to extend its searches to probe uncharted model space through Higgs-to-Higgs cascades.
JHEP02(2020)002   [34]. The upper and right panel show the integrated efficiencies in p T and η, obtained by convoluting the 2D efficiency map with a double-differential distribution of b-jets from tt (blue triangles). These are shown alongside the 1D efficiency maps reported in [34] (pink squares).
in the CMS search, we implemented a Delphes card that parametrised the tagging efficiency using the information reported in [34]. In that study, b-tagging efficiencies and c-and lightjet mis-tag rates are determined using a high purity tt sample and quoted in bins of either p T or η, but not both simultaneously. 1D parametrisations as a function of p T are reported in several p T bins for the three working points. The medium working point has an inclusive b-tagging efficiency of 68%, and inclusive mis-tag rates of 12% and 1.1% for c-and lightjets, respectively. In order to have a better modelling of the b-tagging efficiency over the jet kinematics, we extrapolate the reported efficiencies into a 2D function of jet p T and η.
Our b-tagging efficiency map, shown in figure 15, is determined by a fit to these results using the reported 1D efficiencies as boundary conditions, taking into account the kinematical distributions of b-jets from tt. A tt sample was generated at NLO in QCD with Madgraph5 aMC@NLO, to obtain a 2-dimensional distribution in p T and η with binnings that matched the reported efficiency curves of [34]. The 'unfolded' 2D efficiency map must obey the constraint that it reproduces the reported 1D maps when integrated along either axis as well as the inclusive efficiency when fully integrated over. The integration procedure is a convolution of the binned efficiency map with the double-differential distribution in p T and -29 -JHEP02(2020)002 η. As shown in the upper and right panels of figure 15, the solution is able to satisfy the constraints (the last 10 p T bins were combined due to lack of MC statistics in the tt sample). However, this is clearly an under-constrained problem, and the solutions are found to be somewhat sensitive to the initial conditions of the least-squares minimisation procedure used. The map shown in figure 15 used the average of the corresponding efficiencies in the 1D p T and η maps as a starting point. Taking randomised initial conditions leads to noisy solutions that oscillate around the former. This was verified by averaging over a stochastic sample of solutions with random initial conditions between 0 and 1, observing that the resulting map was within 10-20% of that obtained from the average initial conditions. The variance of the obtained efficiency in the regions where most of the tt sample resided was found to be around 10-20%. This region has a dominant impact on the integrated efficiencies. In bins poorly populated by tt, the results fluctuated more, with a standard deviation of order 50%. These bins, however, do not have a big impact of the overall efficiency.
For the c-and light-jet mis-tag rates, a simpler, more approximate procedure was employed. The 1D efficiencies in p T and η were taken as independent and used to construct a 2D map, for each p T bin where the polynomial parametrisations were provided by the CMS collaboration. In each bin, the efficiency was taken as the product of the p T -dependent parametrisation and a fit to the (inclusive) η efficiency distribution of the medium working point, divided by the integral of said η distribution, such that the integral of the new map matched the integral of the 1D p T parametrisation. These rates are not expected to have any significant impact on our analysis as the multi-jet background is determined by data-driven methods. They therefore did not warrant a high-statistics MC simulation of the non b-jet composition/kinematics in tt that would be required to repeat the procedure employed for the b-tagging.
B H 1 → H 2 H 2 → bbbb limits/projections without systematic uncertainties Here we provide our current ( √ s = 13 TeV, 35.9 fb −1 ) LHC estimates and future HL-LHC, HE-LHC and FCC-hh projections for the 95% C.L. cross section sensitivity for the process pp → H 1 → H 2 H 2 → bbbb assuming no systematic uncertainties, illustrated in figure 16. In addition, we give the values of κ 2 × BR and the improvement in sensitivity I (for HL-LHC, HE-LHC and FCC-hh with respect to the current sensitivity) for the (m H 2 , m H 1 ) plane benchmarks defined previously in table 2, in the absence of SM background systematic uncertainties, in table 3.