Simultaneous Search for Extra Light and Heavy Higgs Bosons via Cascade Decays

Models with extended Higgs sectors can contain several additional Higgs states, heavier or lighter than the SM Higgs boson. The couplings of lighter extra states to SM particles can be strongly reduced, leading to small cross sections for their direct production. Heavier extra states can have larger couplings to SM particles and, moreover, have large branching fractions into lighter extra states, notably into a SM-like Higgs boson accompagnied by another Higgs state which can be lighter or heavier than 125$\sim$GeV. Motivated by corresponding scenarios in the NMSSM we study the prospects for the discovery or exclusion of cascade decays $ggF \to H_3 \to H_2 + H_1$ in the $b\bar{b}b\bar{b}$, $b\bar{b}\tau\tau$ and $b\bar{b}\gamma\gamma$ final states where either $H_1$ or $H_2$ can be SM-like. Significant regions of the NMSSM parameter space can be tested by these searches. These are, however, not confined to models of the NMSSM type.


Introduction
Extended Higgs sectors are frequent properties of models beyond the Standard Model (BSM). Such extra states can have very small couplings to quarks, leptons and SM gauge fields. For instance, for singlets under the SM gauge symmetries such renormalizable couplings are disallowed by gauge invariance. The direct production cross sections for these states are then strongly suppressed in all channels. On the other hand, couplings of singlets to SU (2) Higgs doublets of the SM-or BSM-type are possible and typically present in BSM models. This allows for the discovery of such states in cascade decays of heavy BSM SU(2) Higgs doublets, provided the production cross sections of the latter are large enough.
The final states after BSM-Higgs to BSM-Higgs + SM-Higgs cascades typically correspond to the ones in searches for resonant SM-Higgs (H 125 ) pair production: mainly bbbb, bbτ τ and bbγγ. Corresponding searches have been performed at the LHC by ATLAS [1][2][3][4][5][6] and by CMS [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. However, one of the SM-like Higgs bosons would now be replaced by a lighter or heavier BSM-Higgs boson. One can argue that the cross sections for such processes can be more promising than for resonant SM-Higgs pair production: a) A sizeable gluon-gluon-fusion (ggF ) production cross section of a heavy scalar or pseudoscalar Φ, i.e. a sizeable coupling of Φ to top quarks, requires Φ to possess a sizeable SU(2)-doublet component. However, since H 125 is also a SU(2)-doublet, trilinear couplings Φ − H 125 − H 125 (with Φ a pure doublet) violate the SU(2) symmetry and must be proportional to a SU(2) symmetry breaking vev; the latter is limited from above by the Z/W masses. This limits the possible partial width for Φ → H 125 + H 125 , whereas the concurrent decay Φ → tt is always possible if Φ can be produced in ggF and is heavier than 2 m top . b) In the case ggF → Φ → H 125 + H with Φ a pure doublet, the trilinear coupling Φ − H 125 − H can be SU (2) invariant if H is a singlet. In models with extended Higgs sectors including both an extra doublet and a singlet, such a coupling can thus be much larger than the Z/W masses leading to sizeable Φ → H 125 + H partial widths.
In Two-Higgs-Doublet-Models of type II such as the Minimal Supersymmetric Standard Model (MSSM) the production cross sections for extra CP-even (H) and CP-odd (A) Higgs doublets are not suppressed, and are dominated by ggF for tan β not too large [22][23][24]. H or A can thus play the rôle of Φ above. The Next-to-Minimal Supersymmetric Standard Model (NMSSM) [25,26] contains additional CP-even (H S ) and CP-odd (A S ) singlet-like states with masses typically below M H ∼ M A . One finds that the BR(H → H S + H 125 ) and BR(A → A S + H 125 ) can be up to ∼ 50% [27][28][29][30][31][32], for the reasons given above and detailed in the next section.
In the NMSSM this offers the possibility to produce otherwise practically invisible mostly singlet-like states H S /A S in cascade decays of H/A [27][28][29][30][31][32]. It is the aim of the present paper to study the prospects for discovery or exclusion of, simultaneously, H/A and H S /A S states in ggF → H → H S + H 125 or ggF → A → A S + H 125 in the final states bbbb, bbτ τ and bbγγ. Supersymmetry plays no rôle here, accordingly our results are applicable to any models with similarly extended Higgs sectors; see, e.g., [33].
On the other hand, the analyses presented here are complicated by the fact that the masses M H S /M A S are not known a priori. An important aspect of optimal search strategies are M H S /M A S dependent selection criteria (cuts) on events, hence different analyses should be performed, varying the assumptions on M H S /M A S . Only at the end of each of these analyses a search for a resonance-like bump in the total invariant mass of the H S /A S plus H 125 decay products, which should correspond to M H /M A , is proposed.
In the next section we discuss shortly the Higgs sector of the NMSSM and the couplings relevant for the processes considered here. In section 3 we present features of our signal simulations. In section 4 we discuss the optimal search strategy for the bbbb final state, and compare expected 95% CL upper limits and 5 σ discovery limits on the cross sections times branching fractions to the ones possible in the NMSSM. Sections 5 and 6 are devoted to the bbτ τ and bbγγ final states. All these search strategies and results are identical for ggF → H → H S + H 125 and ggF → A → A S + H 125 , for notational simplicity we will refer to H → H S + H 125 only. In section 7 we conclude with a summary and an outlook.

The neutral Higgs sector of the NMSSM
In this section we discuss briefly some properties of the Higgs sector of the CP-conserving Z 3 -invariant NMSSM. It consists in two SU (2) doublets H u , H d and a complex singlet S. The superpotential of the Higgs sector reads where λ and κ are dimensionless Yukawa couplings, andĤ u ,Ĥ d andŜ denote chiral superfields. Once the real component of the superfieldŜ develops a vacuum expectation value (vev) s, the first term in the superpotential generates an effective µ term µ = λs . (2. 2) The vev v u of H u generates up-type quark masses, the vev v d of H d generates masses for down-type quarks and leptons, and both vevs contribute to the Z and W ± masses. Their ratio is tan β = v u /v d . Decays of a heavy Higgs state into two lighter Higgs states occur in the presence of corresponding trilinear Higgs couplings. Most of the trilinear Higgs couplings in the Z 3 -invariant NMSSM originate from quartic terms in the Higgs potential (see [25,26]) proportional to two powers of λ, κ or the electroweak gauge couplings, once the (neutral) Higgs fields are expanded around their vevs and decomposed into their real and imaginary parts: Hence the trilinear couplings are proportional to the vevs v u , v d or s. where the dimensionful parameter A λ can be much larger than Higgs vevs. In order to obtain its impact on trilinear couplings among Higgs mass eigenstates, the mass matrices have to be diagonalized. In the CP-even sector, where one deals with a 3 × 3 mass matrix, a first step in this direction is a rotation in the SU(2) doublet sector into the so-called Higgs basis where the vev of H is zero, and the vev h = v 2 u + v 2 d is equal to the one of the SM Higgs boson. In fact, in most of the phenomenological acceptable regions of the parameter space of the NMSSM (near the alignment limit [29]) the eigenstates of the full 3 × 3 CP-even mass matrix are not very different from h , H and S r , and will be denoted by H 125 (∼ h , approximately SM-like), H (∼ H , approximately MSSM-like) and H S (∼ S r , approximately singlet-like) in the following.
The corresponding rotation of the imaginary components H 0 u,i and H 0 d,i (with β → −β) diagonalizes their 2 × 2 mass submatrix exactly and generates the Goldstone boson together with the MSSM-like pseudoscalar A . The latter still mixes with the singlet-like S i , but typically both differ little from the mass eigenstates A and A S .
Performing the rotation (2.5) in (2.4) and using the previous approximations in the CPeven and CP-odd sectors, one obtains the trilinear couplings where the omitted terms are suppressed by tan β. Hence, for not too small tan β → 1, trilinear couplings g H 125 HH S and g H 125 AA S are generated which have no analog in the MSSM, and are larger than all other trilinear Higgs couplings if λA λ is large. On the other hand the masses M H/A of the nearly degenerate mostly MSSM-like states H/A are approximatively given by  [29].) For M H S > 250 GeV, decays H S → H 125 +H 125 are possible, leading to double-resonant tri-Higgs production (not considered here). Decays Φ S → τ + + τ − are practically always possible. In the regions in the NMSSM parameter space with all present constraints on the signal rates of H 125 being satisfied the BR(Φ S → γ + γ) is in the 0.1 − 0.3% range, making this decay observable as well. Henceforth we will consider resonant bbbb, bbτ + τ − and bbγγ final states originating from Φ S → bb, Φ S → τ + τ − and Φ S → γγ decays.
Of interest will be the product of cross sections times branching fractions σ(ggF → H/A) × BR(H/A → H 125 + H S /A S → bbbb, bbτ + τ − and bbγγ) for various masses M H , M H S and M A S , for realistic regions in the parameter space of the NMSSM. To this end we have performed scans using the public code NMSSMTools 5.1.0 [34,35] including the radiative corrections from [36]. All phenomenological constraints, including the absence of Landau singularities below the GUT scale and, notably, constraints from Higgs searches in various channels at LEP and LHC are applied.
These include searches for scalar and pseudoscalar Higgs production at LEP (including unconventional Higgs decays), constraints from B-physics, constraints on the mass of H 125 (±3 GeV to account for theoretical uncertainties) and on its signal rates from the combined ATLAS and CMS run I results which disallow too large H 125 − H S mixings, constraints from searches for ggF → H S → γγ for M H S = 65 − 122 GeV, and searches for H/A in the H/A → τ τ channel with H/A produced in association with b-quarks and via ggF . We note that constraints from these latter searches in the M A −tan β plane are weak for tan β ≈ 2−3, typical in the NMSSM, and further alleviated if H/A have large branching fractions into the final states considered here.
The results of these scans for σ(ggF → H/A)×BR(H/A → H 125 +H S /A S → bbbb, bbτ + τ − and bbγγ) will be compared to the sensitivities in different final states in the following sections. The ggF production cross sections for H/A have been obtained from the CERN Yellow Report web page [37] at NNLO+NNLL, after an appropriate rescaling of the H/Agluon-gluon coupling provided by NMSSMTools 5.1.0. Also the branching ratios of H/A and H S /A S are taken from NMSSMTools 5.1.0. In the Figures showing the 95% CL exclusion limits and 5 σ discovery cross sections, viable values for the cross sections times branching fractions in the parameter space of the NMSSM will be indicated as light shaded blue regions for ggF → H → H 125 + H S , and as light shaded red regions for ggF → A → H 125 + A S . For simplicity we will use the notation ggF → H → H 125 + H S for the search strategies; the same search strategies apply to ggF → A → H 125 + A S .

Simulation of Signal Samples
Signal events for the production of H in ggF are generated by MadGraph5 aMC@NLO [38] with matrix elements at NLO taken from aMC SusHi [39][40][41] using the NNPDF2.3NLO PDF set [42]  For each pair (M H , M H S ) we generated 150k unweighted events, more than the expected number of events at 3000 fb −1 . Accordingly the statistical fluctuations from the Monte Carlo (MC) generation are small compared to the expected statistical fluctuations from the data; the latter will be taken into account.
The output is given to the detector simulation Delphes 3 [44]. Jets are clustered with FastJet 3.0.1 [45] using the anti-k T algorithm with ∆R = 0.4. For b-tagging the ATLAS card is used in Delphes 3.
The p T dependence of the b-tagging and mistagging efficiencies is chosen in the ATLAS card following the parametrizations given in [46]. The default value of the parameters in Delphes 3.3.2 correspond to a working point ε b = 70% as, e.g., in the ATLAS search for Higgs pair production in the 4b final state at 13 TeV in [6]. We will employ the same settings except for the bbγγ final state.

Search Strategy for the bbbb Final State
Searches for resonant SM Higgs pair production in the bbbb final state have been performed before by ATLAS at 8 TeV [1,3,4] and at 13 TeV [6], and by CMS at 8 TeV [8,10,12] and at 13 TeV [14,15].
Searches for ggF → H → H S + H 125 → bbbb are complicated by the presence of two unknown masses of H and H S . A naive approach would be to require one bb pair with a mass near 125 GeV, and to look for simultaneous excesses in the plane of invariant masses of the other bb pair and the total 4b invariant mass. However, this approach does not allow to optimize cuts as function of different masses of H and H S . An at least ∼ 20% gain in efficiency can be obtained as follows: a) Choose a tentative value for M H S , and optimise the cuts and the choice of bb pairs as function of this value; b) Search subsequently for an excess in the total 4b invariant mass (suitably corrected, see below).
Subsequently we describe first our simulation of signal samples and the strategies for the analysis. In the following subsections we discuss the background simulation and validation, and finally the results for the expected 95% CL upper limits and 5 σ discovery limits on the cross sections times branching fractions as function of M H and M H S . The latter are compared to possible production cross sections times branching fractions in the NMSSM.

Analyses of Signal Samples
After the simulation of signal samples as described in section 3, at least four b-tagged jets with p T > 30 GeV and |η| < 2.5 are required. Four b-tagged jets can be paired in six different ways. The two invariant masses of two bb pairs are tentatively denoted by M bb (H 125 ) and M bb (H S ). The subsequent procedure depends on the chosen value for M H S , and has to be repeated for each choice.
An event is kept only if a pairing exists for which M bb (H 125 ) is sufficiently close to 125 GeV, and M bb (H S ) is sufficiently close to M H S . In practice, the measured invariant masses of bb pairs are often somewhat smaller than the mass of the decaying Higgs boson. Therefore "sufficiently close to" should better be replaced by "slightly below" as in [6,14]. Generalizing the conditions applied in [6,14], an event is kept only if a pairing exists for which In the case of heavy resonances with masses above ∼ 1 TeV, the two b-jets from a single Higgs boson tend to merge into a single fat jet. Accordingly a "boosted" analysis based on single ∆R = 1.0 jets per Higgs boson was applied in [6] for searches for such heavy resonances. We found that for such heavy H states the production cross sections become too small for reasonable sensitivities, and limit ourselves to M H ≤ 1 TeV subsequently. We had tried to vary the jet reconstruction parameter ∆R without observing, however, a major impact on the possible sensitivities. The next step in the event selection are cuts on p T of the bb pairs associated with H 125 and H S , respectively.
These cuts can be optimised with respect to M 4b (defined by the event) and the tentative value for M H S . To this end we considered many samples of M H ∼ M 4b and M H S . In each case we studied the dependence of the sensitivity s / √ B ( s and B denote the efficiencies of the signal and background, respectively) on the cuts on p T . Maximizing the sensitivities, we obtain different values for the optimal cuts on p T for each sample of M 4b and M H S . These different values for the cuts on p T are well approximated by the following functions of M 4b and M H S : GeV The efficiencies of the three cuts a) 4b with p T > 30 GeV, b) signal region (4.1) and c) dijet p T (4.2) are shown in Figs Here the two b-jets from H S are boosted and hardly get resolved by the standard jet clustering algorithm; dedicated boosted analyses as in [6] could be envoked for this configuration.
For the search for a resonance H we found it useful to replace M 4b by M X , with M X defined such that uncertainties in the measurements of M bb relative to M H 125 and M H S (typically due to radiation out of the jet cones) are corrected: The empirical variable M X was already used by CMS in resonant double Higgs production search in the bbγγ channel [14]. We found that replacing M X by a full Lorentz covariant expression in terms of M H 125 , M H S and the two values of M bb did not improve the sensitivities.
In Fig. 4 we show the distributions of the reconstructed masses M 4b and M X for the signal samples (M H , M H S ) = (350, 195), (500, 310), (750, 450), (1000, 600) GeV. We clearly see a sharpening of the peaks using M X in each case.
In order to obtain the required production cross sections times branching fractions for 5 σ discovery or 95% CL exclusion we have to obtain the background distribution of M X .

Background Estimation
Dominant backgrounds for the bbbb final state are QCD multijet processes including jet misidentifications, and tt [6,14]. The QCD multijet background is difficult to obtain from Monte Carlo simulations alone, and estimated from sidebands in [6,14]. Such data is not available for the different values of M H S studied here, with one exception: For M H S ∼ 125 GeV, the bbbb final state coincides with the one searched for in [6,14].
We simulated bbbb, bbcc, bbjj (with j = c/c) and tt processes as we did for the signal samples in section 3. After applying the cuts of the previous subsection the relative contributions are ∼ 85 − 88% from bbbb, ∼ 7 − 8% from bbcc, ∼ 4 − 8% from tt (depending on M H S and M 4b ). The bbjj contribution is only 0.8% and will subsequently be neglected. We checked that bbbb + jets processes have little impact on the event shapes. Hence we did not simulate them separately, but took a NLO K-factor of 1.7 [38] into account.
In [6] the M 4b distribution of the multijet background has been obtained from a signal free sideband (with two b-tags only, and M bb outside the search window) and appropriate rescaling using data with four b-tags, subtracting the tt contribution. The measured M 4b distribution is available in Fig. 5 in [6], where it is compared to the estimated background. This allows us to proceed in a similar fashion: In order to compare to the data in [6], the previous cuts are slightly modified: The bb pairs are ordered in p T according to M lead bb and M subl bb and, as in [6], the signal region is defined by Still (and expectedly) our simulated background falls below the measured data given in [6]. On the left hand side of Fig. 5 we show the M 4b distribution measured by [6] using 10.1 fb −1 of integrated luminosity, and our MC result with statistical errors. The number of generated MC events corresponds to an equivalent integrated luminosity of ∼13 fb −1 . The statistical error per bin is thus obtained from the number of MC events per bin rescaled by 10/13. The lower panels show the ratio MC/data bin by bin, with the uncertainties from the data and from our MC combined. At least for the interesting region M 4b > ∼ 350 GeV an overall rescaling of our multijet background, like it was performed in [6], seems appropriate. The rescaling factor is obtained taking the average of the data/MC ratio of all bins weighted by the corresponding uncertainties. We obtain a rescaling factor of 1.55 ± 0.27. (The tt background is left untouched, and remains at ∼ 3 − 6%.) The comparison of the M 4b distribution of our background after rescaling to the data from [6] is shown on the right hand side of Fig. 5.
Clearly it is somewhat optimistic to assume that the rescaling of the multijet background by 1.55 ± 0.27 remains valid for M H S = 125 GeV. In the absence of data from sidebands this is, however, the best we can do. Subsequently ± 0.27 will be used as an estimation of the systematic uncertainty of our background for all M H S , a number to be considered as indicative.   For forecasts at 300 or 3000 fb −1 integrated luminosity the statistical uncertainties of the background are much smaller. It is then convenient to fit the shape of the M X background distributions (4.3) after cuts, which will be used in the following, by continuous functions. We found that the best fits are provided by a four parameter Gamma distribution defined in eq.
Expected ± 2σ Figure 7: Expected 95% CL upper limits from ATLAS [6] in blue, their ±2 σ incertainty bands, the expected 95% CL upper limits from our Monte Carlo in black and, for completeness, the 95% CL upper limits obtained from the data.
and combined with the systematic uncertainty (not shown in Figs. 6). It turns out, however, that for forecasts at 300 or 3000 fb −1 the statistical fluctuations are negligibly small relative to the systematic uncertainty from the rescaling by 1.55 ± 0.27.

Future 95% CL Exclusion Limits and 5 σ Discovery Cross Sections
Given the M X distribution of the background for various hypothetical values of M H S and the M X distributions of signals one can, following the statistical methods from [53] and described in the appendix B, obtain values for 95% CL exclusion limits and 5 σ discovery limits for cross sections times branching fractions into the bbbb final state as function of the integrated luminosity, M H and M H S . In the case of an integrated luminosity of 13.3 fb −1 at 13 TeV we can compare the expected 95% CL exclusion limits on cross sections times branching fractions to the ones given by ATLAS in Fig. 11 in [6], for M X = 300 . . . 1000 GeV and M H S ∼ 125 GeV. (This ATLAS search was actually dedicated to spin 2 resonances decaying to SM Higgs pairs, but the differences to spin-0 resonances are expected to be small.) In Fig. 7 we show the expected 95% CL upper limits from ATLAS, their ±2 σ incertainty bands, the expected 95% CL upper limits from our Monte Carlo and, for completeness, the 95% CL upper limits obtained from the data. We see that our expected 95% CL upper limits coincide well with the ones expected by ATLAS.
Since the background was fitted to data at 13 TeV c.m. energy we will show our results also for 13 TeV, for 300 and 3000 fb −1 integrated luminosity.   We recall, however, that the sensitivities to cross sections in Figs. 8 and 9 are model independent and valid for arbitrary (e.g. non-supersymmetric) extensions of the Higgs sector.

Search strategies for the bbτ τ final state
Searches for resonant H 125 pair production in the bbτ τ final state have been performed by ATLAS at 8 TeV [4], and by CMS at 13 TeV in [16,17,19]. Following these searches we concentrate on the τ h τ h , τ h τ e and τ h τ µ modes. As in the case of the bbbb final state we optimise the cuts as function of a tentative value for M H S . A priori the τ τ pair can originate from H S or H 125 ; both cases will be studied below. For the analysis we will make no assumptions on the relative branching ratios BR(H S → bb) and BR(H S → τ τ ). The aim is to obtain separate 95% CL exclusion limits and 5 σ discovery cross sections for the processes ggF → H → H S (→ bb) + H 125 (→ τ τ ), and ggF → H → H S (→ τ τ ) + H 125 (→ bb).

Analyses of Signal Samples
Hence, for each tentative value of M H S two different analyses using different cuts are to be performed, resulting in two (slightly) different distributions of M X .
MadGraph5 aMC@NLO was also used to generate QCD+electroweak bbτ τ events; the LO cross section was rescaled by a NLO K-factor 2.9. After applying the cuts of the previous subsection, the relative contributions of the SM backgrounds depend on M H S and on whether the τ τ pair originates from H S or H 125 , although tt is always dominant: For H S → τ τ the tt contribution increases from ∼ 60% for M H S ∼ 50 GeV to ∼ 100% for M H S > ∼ 350 GeV, the remaining background stems from QCD+electroweak bbτ τ production. For H 125 → τ τ the tt contribution is always ∼ 90%.
In order to validate the background contribution to the M X distribution after cuts we use again a search for H 125 pair production at 13 TeV, now in the bbτ τ channel. Measurements of distributions of the (slightly corrected) total invariant mass by CMS, separately in the bbeτ h , bbµτ h and bbτ h τ h channels, can be found in Fig. 1 in [17].
We have reproduced the cuts in [17] using our background samples. In Figs. 12 we show the measured total invariant mass distribution from Fig. 1 in [17] in black, and our MC results including the statistical uncertainties corresponding to 12.9 fb −1 of integrated luminosity, in the three channels, in orange. Due to the smaller number of events the statistical uncertainties are now larger than in the 4b case.
Still we can ask which rescaling of our simulated background, independent of the total invariant mass and common to all three channels (to improve the statistics), provides a best fit to the data. We find a factor 1.01 ± 0.24, and will subsequently use ±0.24 as an estimate of the systematic uncertainty of the background normalisation. For forecasts at 300 or 3000 fb −1 integrated luminosity the shape of the M X background distributions (5.4) will again be parametrized by continuous functions with M H S dependent parameters: For the tt background the Frechet distribution, and for the bbτ τ background (all channels combined) the GaussExp function already used in [14]. Both functions are defined in the appendix A. In Figs. 13 we show these fits for M H S = 85 GeV and M H S = 350 GeV where the τ τ pair originates from H S .  The following observations can be made: First, the expected sensitivities on cross sections times branching ratios differ hardly among the cases H 125 → bb and H S → τ τ versus H 125 → τ τ and H S → bb; if at all, the analyses aiming at H 125 → bb and H S → τ τ are typically somewhat more sensitive.

Future 95% CL Exclusion Limits and 5 σ Discovery Cross Sections
Second, in Two-Higgs-Doublet models of type II as well as in the NMSSM the branching fractions into bb and τ τ of both H 125 and H S are always related by a factor ∼ 9 : 1. Accordingly the possible cross sections times branching fractions in the NMSSM parameter space for both σ(ggF → H → H 125 + H S → bbτ τ ) and σ(ggF → H → H 125 + H S → τ τ bb), indicated in blue in Figs. 14 -17, are ∼ 1/9 of the ones in Figs. 8 -9 for the bbbb final state. (The same reasoning applies to σ(ggF → A → H 125 + A S → bbτ τ ) and σ(ggF → A → H 125 + A S → τ τ bb); the viable NMSSM points correspond to the ones in Figs. 8 and 9.) Then one can ask, for a given point in parameter space, which of the analyses considered up to now is the most sensitive. According to our results this is the search in the bbbb final mH= 1000GeV Expected ± 2σ state which allows to test somewhat larger regions in parameter space.

Search strategies for the bbγγ final state
Searches for resonant H 125 pair production in the bbγγ final state have been performed by ATLAS at 8 TeV [2,4] and at 13 TeV [5], by CMS at 8 TeV in [7] and at 13 TeV in [18,21]. A priori the diphotons can originate from H S or H 125 ; both cases will be studied below. As in case of the previous final states we optimise the cuts as function of a tentative value for M H S .

Analyses of Signal Samples
For the simulation of signal samples we used again MadGraph5 aMC@NLO [38]. Events are required to have exactly two b-tagged jets with p T (b) > 40 GeV and |η| < 2.5. Following ATLAS [5] a working point with b = 0.85 was chosen for the b-tagging efficiency in order to increase the statistics. At least two photons are required in each event which have to satisfy the isolation criteria  where the sum over i includes all tracks within a cone ∆R = 0.4 around the photon. The two leading photons are required to satisfy Additional cuts depend on whether the diphoton pair is assumed to originate from H 125 or H S , and the assumed value of M H S . First we consider the case H S → γγ. Then the bb pair is assumed to originate from H 125 , and we require 100 GeV < M bb < 150 GeV . mH= 1000GeV Expected ± 2σ p T (bb) > −5.7 GeV + 0.29M X , (6.11) E T (γγ) = 7.45 GeV + 0.33M X . (6.12) All the numerical values above have been obtained by optimising signal-to-background ratios using the MC events.

Background Estimation
SM backgrounds originate from bbγγ, ccγγ, jjγγ (j = c/c), bbjγ and ttH. We simulated these backgrounds again using MadGraph5 aMC@NLO. Due to the relatively large b-tagging efficiency = 0.85 mistagging rates are relatively large. bbjγ contribute if another fake photon appears. After applying the cuts of the previous subsection the relative contributions of the SM backgrounds depend strongly on M H S and M H , and on whether the γγ pair originates from H S or H 125 . For a light H S ∼ 85 GeV the background from bbjγ is always important, ∼ 60% for H S → γγ and ∼ 30% for H S → bb. However, all other γγ + X SM backgrounds can also contribute several 10% individually. The contribution from ttH is always < ∼ 3%.
For bbγγ, ccγγ and jjγγ an extra jet was allowed in the final state, bbjγ was multiplied by a NLO K factor 1.58 from [38], and ttH (with SM couplings) was simulated at NLO.
Still we cannot expect that the total SM background cross section is completely captured by the MC simulation; in searches by ATLAS and CMS in [5,18,21] sidebands are used for this purpose. Thus we proceed as before and correct the total background cross section using data driven methods in the particular case M H S ∼ 125 GeV equivalent to (resonant) SM Higgs pair production. To this end we modify slightly the cuts on p T (b) (> 55 GeV and > 35 GeV for the leading and next-to-leading b-jets), M bb and M γγ such that they coincide with the ones in the ATLAS search [5]. (In the case of M γγ we could check that the efficiency coincides.) Extrapolating from sidebands with less b-jets, ATLAS [5] obtained 1.63 ± 0.3 expected background events in the signal region, whereas we found 1.06 ± 0.14 events from summing all MC events. Since all our backgrounds were simulated to similar order (NLO) in the QCD coupling and separate higher order K factors are not available, we multiply their sum by 1.54 ± 0.35 where the latter uncertainty will again be treated as an estimate of the systematic uncertainty contributing to our final results.
The dependence of the background on the total invariant mass M bbγγ was parametrized in [5] by a two parameter Landau distribution given in appendix A. We found that the Landau distribution fits the M X distribution from eq. (6.8) as well (with M H S dependent parameters), and used it for the expected M X distributions of the various background contributions for our forecasts at 300 fb −1 and 3000 fb −1 .  [21] for resonant Higgs pair production in the channel H 125 +H 125 → bbγγ based on 35.9 fb −1 appeared. The expected 95% CL exclusion limits given in [21] can be compared to ours for M H S = 125 GeV for the same integrated luminosity; this comparison as function of M H is shown in Fig. 18. The expected limits coincide within 1 σ for M H > ∼ 500 GeV, and within 2 σ everywhere. Our expected limits are systematically more conservative; we note that the CMS analysis employs a trained boosted decision tree in order to separate the signal from the backgrounds which is not available here.

Future 95% CL Exclusion Limits and 5 σ Discovery Cross Sections
Our expected 95% CL exclusion limits and 5 σ discovery cross sections for H 125 → bb, As before viable NMSSM regions for scalar production are shown in shaded blue, for pseudoscalar production in shaded red in case they potentially exceed the ones for scalar production. Again a sizeable region in the NMSSM parameter space can be tested in this final state provided M H is not too large. It is, however, not the same region potentially visible in the bbbb final state: The branching fraction of H S into γγ can vary in the 0.2% ± 0.1% range, and is anticorrelated with its branching fraction into bb. Moreover, as it is visible from the shaded red regions in Figs. 19 and 20, the signal rates σ(ggF → A → H 125 + A S → bbγγ) can be relatively large. These correspond to very singlet-like pseudoscalars A S with very suppressed couplings to quarks and leptons, but sizeable coupling ∼ λ to higgsinos. Then the charged higgsino-loop induced coupling to diphotons can dominate, leading to a large BR(A S → γγ). The coupling A − A S − H 125 is not suppressed in this case, leading to potentially large signal rates.
On the experimental side, the comparison of the upper limits on H 125 → bb and H S /A S → γγ versus H 125 → γγ and H S /A S → bb has a simple answer depending on M H S /A S : For M H S /A S < 125 GeV the search for H S /A S → γγ is sensitive to smaller signal rates, whereas for M H S /A S > 125 GeV the search for H 125 → γγ, H S /A S → bb is typically sensitive to smaller signal rates. However, different regions in the parameter space of underlying models are tested by these searches.

Conclusions and Outlook
Searches for resonant SM Higgs pair production are performed with considerable effort by ATLAS and CMS. As explained in the introduction searches for ggF → Φ → H 1 + H 2 can be more promising where either H 1 or H 2 can be SM-like, and the other state being possibly CP-odd (which does not affect the search methods). This scenario is manifest in the NMSSM where the rôle of Φ is played by the MSSM-like heavy doublet, but the argument is more general. In the present paper we have studied the prospects for corresponding searches in the bbbb, bbτ τ and bbγγ final states, including SM backgrounds. We found that significant regions in the NMSSM parameter space can be tested by these searches: The NMSSM specific parameters testable by bbbb are typically in the region λ ∼ 0.50 − 0.70 (the conservative upper bound from the absence of a Landau singularity below M GU T ), GeV. An exception is the case M H S ∼ 85 GeV − 110 GeV where LEP constraints on the coupling of H S to the Z boson are somewhat weaker; here values of λ down to 0.16 and µ down to 100 GeV (together with tan β up to 4.5) can lead to testable points. These ranges of potentially testable parameters depend little on the total integrated luminosity, but higher luminosity increases of course the number testable parameters within these ranges. Most of the parameters testable by bbγγ are in the same region except for κ ∼ 0.08 − 0.3, A κ ∼ −50 GeV − 10 GeV which indicates that this region is non-overlapping  with the one testable by bbbb.
We are convinced that the here proposed search methods can still be refined, and that the estimated sensitivities to cross sections times branching fractions presented here are conservative. This becomes clear from a comparison to the recent CMS search for resonant SM-Higgs pair production [21] in the bbγγ final state (and actually also from a comparison to the recent CMS search [19] in the bbτ τ final state). Thus we hope that such promising searches will be performed in the future at the LHC.     In this appendix we sketch the computations of 5 σ discovery and 95% CL exclusion limits based on the M X distributions of the background, and the different M X signal distributions (depending on M H S ), following [53]. As shown in Fig. 4 the M X distribution of the signal, after event selection and cuts, corresponds to a certain number s i of expected signal events per bin. (The bin size in M X is 20 GeV, and we have checked that the final results do not vary with this size.) From the MC simulation we know how this event number depends on the total signal cross section σ sig , the integrated luminosity L int and the acceptance A times efficiency : where ∼ indicates a bin-dependent proportionality factor < 1. The backgrounds after event selection and cuts have been fitted by continuous functions of f (M X ) normalized to 1. Thus the number b i of expected background events per bin per integrated luminosity is Due to the large number of simulated events the statistical uncertainties around the median values A · s and A · b are negligibly small, < ∼ 1% in practically all cases. Subsequently we use likelihood functions with (b i + σ b s i )! interpolated by the Γ function for non integer (b i + σ b s i ). For 5 σ discovery limits on σ sig we look for the value of σ sig for which the background only hypothesis is rejected at the 5 σ level. The observed number of events per bin would be b i + σ sig s i , and the likelihood function corresponding to the background only hypothesis is L(0, σ sig ). As function of the number of events per bin it has its maximum at L(σ sig , σ sig ). Hence the test statistics t disc for discovery is t disc = −2 ln L(0, σ sig ) L(σ sig , σ sig ) .

(B.4)
Following [53] the significance Z disc is then and for a 5 σ discovery we determine σ sig such that Z disc = 5. For 95% CL exclusion limits on σ sig we look for the value of σ sig for which the signal hypothesis is rejected at 95% CL. The observed number of events per bin would be b i , and the likelihood function corresponding to the signal hypothesis is L(σ sig , 0). As function of the number of events per bin it has its maximum at L(0, 0). Hence the test statistics t excl for exclusion is For exclusion at 95% CL we determine σ sig such that √ t excl = 1.64 since we consider only positive signal contributions to the number of events.
Uncertainties from the background are estimated as follows: ±(1 − 2) σ statistical uncertainties can be obtained bin by bin. To these we add linearly (to be conservative) the estimated (1−2) σ systematic uncertainties from the overall normalisation of the background. Then the above likelihoods are recomputed with correspondingly larger and smaller values for all b i from which we deduce the ±(1 − 2) σ uncertainties of σ sig for the 5 σ discovery limits and 95% CL exclusion limits.