Simultaneous Search for Extra Light and Heavy Higgs Bosons via Cascade Decays
Ulrich Ellwanger Matias Rodriguez-Vazquez
hep-ph hep-ex

=1














 









10pt
-0.0cm
0.6cm
0cm 
16.5cm
22.5cm
6mm
myheadings













addtoresetequationsection




 0.3ex<-0.75em-1.1ex∼ 
 0.3ex>-0.75em-1.1ex∼ 

    [ κ̨łλβ̱θ
        
        definitiondefinitionDefinition[section]
        
        
        
        LPT Orsay 17-33 
        
        
        
         Simultaneous Search for Extra Light and Heavy Higgs Bosons
        via Cascade Decays
        Ulrich Ellwanger^a and Matías Rodríguez-Vázquez^a
        
        ^a Laboratoire de Physique Théorique, UMR 8627, CNRS,
        Université de Paris-Sud, Université
        Paris-Saclay, 91405 Orsay, France
        
        
        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 GeV.
        Motivated by corresponding scenarios in the NMSSM we study the prospects for
        the discovery or exclusion of cascade decays ggF → H_3 → H_2 + H_1 in the
        bb̅bb̅, bb̅ττ and bb̅γγ 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 
        bb̅bb̅, bb̅ττ and bb̅γγ.
        Corresponding searches have been performed at the LHC by ATLAS
        <cit.> and by CMS 
        <cit.>. 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
        <cit.>.
        H or A can thus play the rôle of Φ above.
        The Next-to-Minimal Supersymmetric Standard Model (NMSSM)
        <cit.> 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%<cit.>, 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<cit.>.
        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 bb̅bb̅,
        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., <cit.>.
        
        We will adopt various strategies from the searches for resonant SM
        Higgs pair production by ATLAS
        <cit.> and by CMS 
        <cit.>. Moreover,
         for M_H_S near 125 GeV we can
        compare our backgrounds and expected 95% CL upper limits on the
        cross sections times branching fractions to the ones obtained in these
        publications.
        
        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 bb̅bb̅ 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 ℤ_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
        
        
        W_Higgs=λŜĤ_u·Ĥ_d+κ^3/3Ŝ^3
        
        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  .
        
        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 corres­ponding trilinear Higgs couplings. Most of the trilinear Higgs couplings in the
        ℤ_3-invariant NMSSM originate from quartic terms in the Higgs
        potential (see <cit.>) 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:
        
        H^0_u= v_u + 1/√(2)(H^0_u,r + i H^0_u,i),   
        H^0_d= v_d + 1/√(2)(H^0_d,r + i H^0_d,i),   
        S=s+1/√(2)(S_r +i S_i)  .
        
        Hence the trilinear couplings are proportional to the vevs v_u, v_d or s.
        Whereas v_u, v_d are limited from above by M_Z^2 = g_1^2+g_2^2/2(v_u^2+v_d^2),
        a large vev s can generate
        a trilinear coupling ∼ H_u· H_d S. Another source for such a
        coupling is a trilinear Higgs-dependent soft SUSY breaking term
        λ A_λ H_u· H_d S + h.c.
        
        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
        
        H^0_u,r = sinβ h' - cosβ H',     H^0_d,r = cosβ h' + sinβ H'
        
        where the vev of H' is zero, and the vev <h'>=√(v_u^2+v_d^2)
        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 <cit.>) 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 (<ref>) in (<ref>) and using the
        previous approximations in the CP-even and CP-odd sectors, one obtains the trilinear couplings
        λ A_λ/√(2)(tan^2β-1/tan^2β+1 
        H_125(H H_S - A A_S) + …)
        
        where the omitted terms are suppressed by tanβ. Hence, for not too small tanβ→ 1,
        trilinear couplings g_H_125 H H_S and g_H_125 A A_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
        
        M_H/A^2 ∼μ(A_λ+κ/λμ)
        1+tan^2β/tanβ ,
        
        which limits A_λ from above for fixed M_H/A, tanβ,
        small |κ/λ| and |μ|  100 GeV (as
        required by the lower LEP bound on higgsino-like charginos). Accordingly the
        trilinear couplings g_H_125 H H_S and g_H_125 A A_S can be
        larger for larger M_H/A^2.
        
        The production cross sections for the mostly 
        MSSM-like states H/A is dominated by ggF<cit.>;
        at √(s)=13-14 TeV and for tanβ∼ 2-3 (typical in the NMSSM) they are
        O(1 pb) up to M_H/A∼ 600 GeV. The trilinear couplings 
        g_H_125 H H_S and g_H_125 A A_S
        induce the decays H → H_125+ H_S and A → H_125+ A_S if kinematically allowed. The branching fractions
        BR(H/A → H_125+ H_S/A_S) can be as large as ∼ 50%, in contrast to the
        decay H → H_125 +H_125.
        
        The singlet-like states Φ_S = H_S/A_S
         have small couplings to quarks, leptons and gauge fields induced by mixings
        with h', H' and A'. Hence the production cross sections for Φ_S
        are typically small, and their discovery may have to rely on
        H/A → H_125+ Φ_S decays. Via the couplings induced by mixing, Φ_S can decay
        into the same channels as H_125 and H/A. For M_Φ_S > 2m_top, decays
        into tt̅ are dominant, whereas decays Φ_S→ bb̅ dominate typically for
        M_Φ_S < 2m_top. (The branching ratio for H_S→ W^+ W^- can also be
        sizeable <cit.>.)
        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 bb̅bb̅, 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
         → bb̅bb̅,  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
        <cit.> including
        the radiative corrections from <cit.>. 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
         → bb̅bb̅,  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 <cit.>
        at NNLO+NNLL, after an appropriate rescaling of the H/A-gluon-gluon
        coupling provided by . Also the branching ratios
        of H/A and H_S/A_S are taken from .
        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 <cit.> with matrix elements at NLO taken from <cit.> using the NNPDF2.3NLO PDF
        set <cit.> in the 4-flavour scheme. Renormalization and factorization
        scales are chosen as H_T/2 on an event-wise basis. <cit.> is used for the H → H_125+ H_S decays, the
        H_125 and H_S decays and the subsequent showering and hadronization.
        The total widths of H are below M_H/50 in all cases, below M_H/100 for
        M_H < 500 GeV, hence the narrow width approximation is well satisfied.
        
        Separate signal events have been generated for each pair (M_H, M_H_S).
        For M_H we chose M_H=425, 500, 625, 750 and 1000 GeV.
        We varied M_H_S in steps of 10 GeV in the range 25-225 GeV,
        and in steps of 20 GeV above 225 GeV up to the kinematic boundary
        M_H-125 GeV (except for M_H=1000 GeV where M_H_S was varied in
        steps of 25/50 GeV).
        
        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 <cit.>. Jets are clustered with <cit.> using the anti-k_T algorithm with Δ R = 0.4.
        For b-tagging the ATLAS card is used in . 
        
        The p_T dependence of the b-tagging and mistagging efficiencies is
        chosen in the ATLAS card following the parametrizations given in
        <cit.>. The default value of the parameters in
         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 <cit.>. We will
        employ the same settings except for the bb̅γγ final state.
        
        
        
        § SEARCH STRATEGY FOR THE BB̅BB̅ FINAL STATE
        
        
        Searches for resonant SM Higgs pair production in the bb̅bb̅ final state
        have been performed before by ATLAS at 8 TeV <cit.> and at 13 TeV <cit.>, and
        by CMS at 8 TeV <cit.> and at 13 TeV
        <cit.>.
        
        Searches for ggF → H → H_S + H_125→ bb̅bb̅ 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 <ref>,
        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 <cit.>. Generalizing the conditions applied in
        <cit.>, an event is kept only if a pairing exists for which 
        χ = √((
        M_bb̅(H_125)-115 GeV/0.1  M_bb̅(H_125))^2 +
        (M_bb̅(H_S)-0.85  M_H_S/0.1  M_bb̅(H_S))^2)
        < 2  .
        
        If different pairings within a given event satisfy (<ref>), the combination
        that minimizes χ is chosen.
        In Figs. <ref> we show, for M_H=500 GeV, the distributions of the dijet masses
        M_bb̅(H_125) and M_bb̅(H_S) for the pairing minimizing χ,
        for two different benchmark points M_H_S=85 GeV and M_H_S=340 GeV in the case where
        M_H_S for the analysis was chosen correctly. The black
        contours indicate the signal regions defined by χ < 2.
        [An alternative pairing algorithm based on the angular distances of the
        constituent b-jets was studied. We found that the resulting sensitivities are
        mostly similar to those obtained using the mass-based algorithm described above when the
        jets are sufficiently collimated. However, as M_H_S increases, the jets become
        too back-to-back, making the angular-pairing contraproductive. For this reason, we used the
        mass-based algorithm for the reconstruction of the Higgs candidates.]
        
        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 <cit.> 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:
        
        p_T(bb̅(H_125))    >    1.6 GeV + 0.4  M_4b-0.13  M_H_S - M_H_S/M_4b· 160 GeV
        
        p_T(bb̅(H_S))    >    12 GeV + 0.4  M_4b-0.15  M_H_S - M_H_S/M_4b· 166 GeV
        
        The efficiencies of the three cuts a) 4b with p_T > 30 GeV, b) signal region
        (<ref>) and c) dijet p_T (<ref>) are shown in Figs. <ref>
        for M_H_S=75 GeV and M_H_S=325 GeV as function of M_H, and in
        Figs. <ref> for M_H=625 GeV and M_H=1000 GeV as function of M_H_S.
        One observes a decrease of the efficiency of the condition of four b-tagged jets
        for M_H_S 80 GeV (for M_H=625 GeV) and M_H_S 100 GeV (for M_H=1000 GeV).
        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 <cit.>
        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:
        
        M_X = M_4b+125 GeV-M_bb̅(H_125) + M_H_S-M_bb̅(H_S)  .
        
        The empirical variable M_X was already used by CMS in resonant double Higgs
        production search in the bb̅γγ channel <cit.>.
        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. <ref> 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 bb̅bb̅ final state are QCD multijet processes
        including jet misidentifications, and tt̅<cit.>.
        The QCD multijet background is difficult to obtain from Monte Carlo simulations
        alone, and estimated from sidebands in <cit.>.
        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 bb̅bb̅ final state
        coincides with the one searched for in <cit.>.
        
        We simulated bb̅bb̅, bb̅cc̅, bb̅jj (with
        j≠ c/c̅) and tt̅ processes as we did for the signal samples
        in section <ref>.
        After applying the cuts of the previous subsection the relative contributions are
        ∼ 85-88% from bb̅bb̅, ∼ 7-8% from bb̅cc̅, 
        ∼ 4-8% from tt̅ (depending on M_H_S and M_4b).
        The bb̅jj contribution is only 0.8%
        and will subsequently be neglected. We checked that bb̅bb̅+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 <cit.> into account.
        
        In <cit.> 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 <cit.>, 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 <cit.>, the previous cuts are slightly
        modified: The bb̅ pairs are ordered in p_T according to M_bb̅^lead
        and M_bb̅^subl and, as in <cit.>, the signal region
        is defined by
        χ = √((
        M_bb̅^lead-120 GeV/0.1  M_bb̅^lead)^2 +
        (M_bb̅^subl-0.85  M_H_S/0.1  M_bb̅^subl)^2)
        < 1.6  .
        
        
        Still (and expectedly) our simulated background falls below the measured data
        given in  <cit.>.
        On the left hand side of Fig. <ref> we show the M_4b
        distribution measured by <cit.> 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 <cit.>, 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 <cit.>
        is shown on the right hand side of Fig. <ref>.
        
        
        
        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 (<ref>)
        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. (<ref>) with M_H_S dependent fit parameters.
        (The bb̅bb̅, bb̅cc̅, and tt̅ background contributions
        to the M_X distributions were fitted separately.)
        In Figs. <ref> we show the sum of these fits for M_H_S=85 GeV
        and M_H_S=350 GeV.
        
        
        
        Of course the remaining statistical fluctuations of the background can still be evaluated and
        combined with the systematic uncertainty (not shown in Figs. <ref>).
        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 <cit.>
        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 bb̅bb̅ 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 <cit.>, 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. <ref> 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 choose four representative values for M_H=425, 500, 750 and 1000 GeV,
        and show the 95% CL exclusion limits and 5 σ discovery cross sections
        as function of M_H_S in each case. For 300 fb^-1 integrated luminosity
        these are shown in Figs. <ref>, for 3000 fb^-1 integrated
        luminosity in Figs. <ref>.
        
        The expected limits become weaker for M_H_S 50 GeV (for M_X=425-500 GeV)
        and M_H_S 100 GeV (for M_H=1000 GeV). As stated in subsection <ref>
        here the bb̅ pair from H_S becomes too boosted and is no longer resolved
        by the standard jet clustering algorithm.
        
        The shaded blue regions in Figs. <ref> and <ref>
        indicate viable values for the cross sections times branching fractions
        for σ(ggF→ H→ H_125+H_S→ bb̅bb̅) in
        the parameter space of the NMSSM, see section <ref>.
        Typically the viable values for σ(ggF→ A→ H_125+A_S→ bb̅bb̅)
        are smaller; if not we show them as shaded red regions.
        In the region of the NMSSM parameter space corresponding to M_H  500 GeV,
        the partial width for H_S→ H_125+H_125 becomes relatively large (≈ 10 MeV)
        if kinematically allowed. As a consequence the branching fractions of H_S into
        bb̅ (and the other channels considered in this paper) decrease, leading
        to a decrease of the possible production cross sections times branching fractions
        for M_H  500 GeV, M_H_S 250 GeV.
        
        
        
        
        The following conclusions can be drawn from Figs. <ref>
        and <ref>:
        For M_H 500 GeV wide ranges of M_H_S
        in the NMSSM parameter space can be discovered
        or, at least, excluded. For larger M_H testable regions
        in the NMSSM parameter space exist, but for M_H∼ 1 TeV
        only for 3000 fb^-1 integrated luminosity.
        In fig. <ref> 
        we summarize these results showing the 95% C.L. expected upper limits in the M_H vs M_H_S plane.
        
        
        
        
        We recall, however, that the sensitivities to cross sections in
        Figs. <ref> and <ref> 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 <cit.>, and
        by CMS at 13 TeV in <cit.>.
        Following these searches we concentrate on the τ_hτ_h, τ_hτ_e
        and τ_hτ_μ modes. As in the case of the bb̅bb̅ 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
        
        
        For the simulation of signal samples the same series of codes as for the
        bb̅bb̅ final state was used, see section <ref>.
        Events are required to have exactly two b-tagged jets with p_T(b) > 30 GeV
        and |η|<2.5. For the b-tagging efficiency a working point with
        ε_b = 70% is chosen.
        If the event has exactly two hadronic taus τ_h, both are required to have
        p_T(τ_h) > 45 GeV. Events with one hadronic tau are required to have
        exactly one additional isolated lepton ℓ = e,μ of opposite charge and with 
        p_T(ℓ) > 20 GeV. For the transverse mass m_T^ℓ of leptons we
        require
        
        m_T^ℓ≡√(2p_T(ℓ) E_T^miss(1-cos(ϕ(E_T^miss)
        -ϕ(ℓ)))) < 40 GeV  .
        
        All objects are required to have |η|<2.47.
        In the case of τ_h τ_ℓ final states the invariant mass M_ττ
        was reconstructed using the collinear mass, i.e. assuming that the neutrino
        from the τ_ℓ decay is emitted collinear to the lepton and responsible
        for all E_T^miss. 
        
        Considering first the case where the ττ pair originates from H_S
        with an assumed mass M_H_S (and hence that M_bb̅ should be close to
        125 GeV), a M_H_S dependent signal region is defined
        by
        χ = √((
        M_bb̅-110 GeV/0.35· M_bb̅)^2 +
        (M_ττ-0.92· M_H_S/Max(0.35· M_ττ, 30 GeV))^2)
        < 1  .
        
        
        If different pairings within a given event satisfy (<ref>), the combination
        that minimizes χ is chosen.
        In Figs. <ref> we show, for M_H=500 GeV, the distributions of the dijet mass
        M_bb̅(H_125) and M_ττ(H_S) for the pairing minimizing χ,
        for two different benchmark points M_H_S=85 GeV and M_H_S=340 GeV in the case where
        M_H_S for the analysis was chosen correctly. The black
        contours indicate the signal regions defined by χ < 1.
        
        
        
        Cuts on the transverse momenta of Higgs candidates are chosen as
        
        p_T(bb̅)    >    52 GeV + 0.14  M_X-0.2  M_H_S - M_H_S/M_X·
        202 GeV
        
        p_T(ττ)    >    24 GeV + 0.19  M_X-0.02  M_H_S - M_H_S/M_X·
        128 GeV
        with
        
        M_X = M_bb̅ττ + 125 GeV - M_bb̅  .
        
        The numerical coefficients in eqs. (<ref>)–(<ref>) were
        obtained by optimizing the relative signal to background efficiency.
        
        In the case where the bb̅ pair originates from H_S
        with an assumed mass M_H_S (and hence that M_ττ should be close to
        125 GeV), the signal region (<ref>) is replaced by
        χ = √((
        M_ττ-120 GeV/0.35· 120 GeV)^2 +
        (M_bb̅-0.85· M_H_S/Max(0.35· M_H_S, 35 GeV))^2)
        < 1
        
        and the cuts on the transverse momenta of Higgs candidates are
        
        p_T(ττ)    >    118 GeV + 0.02  M_X-0.55  M_H_S - M_H_S/M_X·
        380 GeV
        
        p_T(bb̅)    >    16 GeV + 0.19  M_X-0.02  M_H_S - M_H_S/M_X·
        137 GeV
        with
        
        M_X = M_bb̅ττ + M_H_S - M_bb̅  .
        
        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.
        
        
        
         §.§ Background Estimation
        
        
        Backgrounds originate from tt̅ (and single top) and QCD+electroweak
        bb̅ττ production. Contributions from one or more jets misidentified as
        τ_h are seen to become small after the cuts on p_T(ττ).
        We have generated 1.5· 10^7tt̅ events using  <cit.>; the LO cross section was rescaled by a (NNLO+NNLL) K-factor 1.7
        obtained from <cit.>.
         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 bb̅eτ_h, bb̅μτ_h and bb̅τ_hτ_h
        channels, can be found in Fig. 1 in <cit.>.
        
        We have reproduced the cuts in <cit.> using our background
        samples. In Figs. <ref> we show the measured total invariant
        mass distribution from Fig. 1 in <cit.> 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 (<ref>)
        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 <cit.>. Both functions are defined in the appendix A.
        In Figs. <ref> we show these fits for M_H_S=85 GeV and
        M_H_S=350 GeV where the ττ pair originates from H_S.
        
        
        
         §.§ 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 we can, as before,
        obtain values for 95% CL exclusion and 5 σ discovery for cross sections
        times branching fractions into the H_125→ bb̅, H_S→ττ and
        H_125→ττ, H_S→ bb̅ final states as function
        of the integrated luminosity, M_H and M_H_S.
        We choose four representative values for M_H=425, 500, 750 and 1000 GeV,
        and show the 95% CL exclusion limits and 5 σ discovery cross sections
        as function of M_H_S in each case. For H_125→ bb̅, H_S→ττ
        at 300 fb^-1 integrated luminosity
        these are shown in Figs. <ref>, for 3000 fb^-1 integrated
        luminosity in Figs. <ref>. For H_125→ττ,
        H_S→ bb̅ at 300 fb^-1 integrated luminosity
        these are shown in Figs. <ref>, for 3000 fb^-1 integrated
        luminosity in Figs. <ref>.
        The uncertainties include
        statistical uncertainties and, added linearly, ± 0.24 considered as an estimate of the
        systematic uncertainty originating from the normalisation of the background.
        
        
        
        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.
        
        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. <ref>–<ref>,
        are ∼ 1/9 of the ones in
        Figs. <ref>–<ref> for the
        bb̅bb̅ 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. <ref> and
        <ref>.)
        
        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 bb̅bb̅ final 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 <cit.>
        and at 13 TeV <cit.>,
        by CMS at 8 TeV in <cit.> and at 13 TeV in <cit.>.
        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 <cit.>.
        Events are required to have exactly two b-tagged jets with p_T(b) > 40 GeV
        and |η|<2.5. Following ATLAS <cit.>
        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
        ∑_i p_T,i/p_T,γ < 0.1
        
        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
        
        E_T^lead/M_γγ>0.35,      E_T^subl/M_γγ>0.25,     
        |η|<2.37  .
        
        
        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  .
        
        As in the previous cases it is useful to define a corrected invariant
        mass M_X for the bb̅γγ system:
        
        M_X=M_bb̅γγ+125 GeV - M_bb̅  .
        
        
        Given an assumed value for M_H_S it turned out to optimize the signal
        to background
        efficiency applying a M_X and M_H_S dependent cut on the diphoton
        invariant mass M_γγ, since the measured distribution
        |M_γγ - M_H_S| broadens somewhat with M_X:
        
        |M_γγ - M_H_S| < 4.3 GeV + 0.016 M_X  .
        
        Also the cuts on the H_125 and H_S candidates depend on M_X:
        
        p_T(bb̅)    >    17 GeV + 0.18 M_X  ,
        
        E_T(γγ)    =    68.3 GeV + 0.25 M_X  .
        
        (E_T(γγ) instead of p_T(γγ) allows for M_H_S
        independent cuts.)
        
        Assuming H_125→γγ and H_S→ bb̅ the previous cuts
        are modified as follows: First, for M_X we take
        
        M_X=M_bb̅γγ+ M_H_S - M_bb̅  ,
        
        and for the diphoton invariant mass we require
        
        |M_γγ - 125 GeV| < 2 GeV + 0.02 M_X  .
        
        The mass window for M_bb̅ is now
        
        0.9 M_H_S - 30 GeV < M_bb̅ < 0.9 M_H_S + 20 GeV  .
        
        The cuts on the H_125 and H_S candidates are:
        
        p_T(bb̅)    >    -5.7 GeV + 0.29 M_X  ,
        
        E_T(γγ)    =    7.45 GeV + 0.33 M_X  .
        
        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̅), bb̅jγ and tt̅H. We
        simulated these backgrounds again using .
        Due to the relatively large b-tagging efficiency ϵ=0.85 mistagging
        rates are relatively large. bb̅jγ 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 bb̅jγ 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 tt̅H is always
        3%.
        
        For bb̅γγ, cc̅γγ and jjγγ
        an extra jet was allowed in the final state, bb̅jγ was multiplied
        by a NLO K factor 1.58 from <cit.>, and tt̅H (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 <cit.>
        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
         <cit.>.
        
        (In the case of M_γγ we could check that the efficiency coincides.)
        Extrapolating from sidebands with less b-jets, ATLAS <cit.>
        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 <cit.> by a two parameter Landau
        distribution given in appendix A. We found that the Landau distribution fits the
        M_X distribution from eq. (<ref>) 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.
        
        
        
        
         §.§ 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 we can, as before,
        obtain values for 95% CL exclusion and 5 σ discovery for cross sections
        times branching fractions into the H_125→ bb̅, H_S→γγ and
        H_125→γγ, H_S→ bb̅ final states as function
        of the integrated luminosity, M_H and M_H_S.
        
        After completing our analysis the CMS search <cit.> 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 <cit.>
        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. <ref>. 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.
        
        
        
        Our expected 95% CL exclusion limits and 5 σ discovery cross sections
        for H_125→ bb̅, H_S→γγ as function of M_H_S
        at 300 fb^-1 integrated luminosity
        are shown in Figs. <ref> for
        four representative values for M_H=425, 500, 625 and 750 GeV,
        and for 3000 fb^-1 integrated
        luminosity in Figs. <ref>. For H_125→γγ,
        H_S→ bb̅ at 300 fb^-1 integrated luminosity
        these are shown in Figs. <ref>, for 3000 fb^-1 integrated
        luminosity in Figs. <ref>.
        As in the bb̅bb̅ case, we also present the results for the
        expected 95% C.L. upper limits in the M_H_125 vs M_H_S plane in
        fig. <ref>.
        
        
        
        
        
        
        The uncertainties include
        statistical uncertainties and, added linearly, a factor ± 0.35 considered as
        systematic uncertainty originating from the normalisation of the background.
        
         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
        bb̅bb̅ 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. <ref> and
        <ref>, 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
        bb̅bb̅, 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 bb̅bb̅ are typically in the region
        λ∼ 0.50 - 0.70 (the conservative upper bound from the absence of
        a Landau singularity below M_GUT), κ∼ 0.4 - 0.7, tanβ∼ 2 - 3.5,
        μ_eff∼ 180 GeV - 290 GeV,
        A_λ∼ 150 GeV - 550 GeV,
        A_κ∼ -830 GeV - -210 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 bb̅bb̅.
        
        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 <cit.> in the bb̅γγ
        final state (and actually also from a comparison to the recent CMS search
        <cit.> in the bb̅ττ final state). Thus we hope
        that such promising searches will be performed in the future at the LHC.
        
        
        
        § ACKNOWLEDGEMENTS
        
        
        This project has received support from the European Union's Horizon 2020 research and
        innovation programs 
        ITN Higgs­Tools (PITN-GA-2012-316704),
        ITN Elusives (Marie Sklodowska-Curie grant agreement No 674896),
        RISE InvisiblesPlus (Marie Sklodowska-Curie grant agreement No 690575),
        RISE NonMinimalHiggs (Marie Sklodowska-Curie grant agreement No 645722),
        the ERC advanced grant Higgs@LHC,
        and the Défi InPhyNiTi project N2P2M-SF.
        M.R.V. thanks Marius Wiesemann, Fady Bishara, Dirk E. Zerwas,
        Nikola Makovec, Sophie Henrot,  Alberto Escalante and Davide
        Napoletano for fruitful discussions and help, and the hospitality received
        at the IPPP at the University of Durham, where part of this work has been
        done.
        
        
        
        
        § APPENDIX A: FUNCTIONS TO FIT BACKGROUND DISTRIBUTIONS
        
        
        In this appendix we define the functions used to parametrize the
        total invariant mass (or M_X) distributions of the backgrounds to
        the various final states. The best choice among the functions and the values of the
        corresponding parameters have been obtained by maximum likelihood estimates.
        
        
        
         §.§ Four Parameter Gamma distribution:
        
        f(M_X;α,β,γ,μ)=   γ/βΓ(α)
        e^-((M_X-μ)/β)^γ(M_X-μ/β)^αγ-1    for M_X≥μ ,
        
        f(M_X;α,β,γ,μ)=   0    for   M_X ≤μ ,
        
        
         §.§ Frechet distribution:
        
        f(M_X;α,β,μ)=   α/β
        e^-((M_X-μ)/β)^-α(M_X-μ/β)^-α-1    for M_X≥μ ,
        
        f(M_X;α,β,μ)=   0    for   M_X ≤μ ,
        
        
         §.§ GaussExp function:
        
        
        A Gaussian with an exponential tail:
        
        
        f(M_X;μ,σ,k)=   
        e^-(M_X-μ)^2/2σ^2    for M_X-μ/σ≤ k ,
        
        f(M_X;μ,σ,k)=   
        e^k^2/2-k(M_X-μ)/σ    for M_X-μ/σ > k .
        
        
         §.§ Landau Distribution
        
        f(M_X;μ,σ)=
        ∫_0^∞sin(2t)
        e^-t(M_X-μ)/σ-2t/πlog(t)dt
        
        
        § APPENDIX B: 5 Σ DISCOVERY AND 95% CL EXCLUSION LIMITS
        
        
        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 <cit.>.
        
        As shown in Fig. <ref> 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 ϵ:
        
        s_i ∼ L_int·σ_s·<A·ϵ>_s
        
        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
        
        b_i =  L_int·σ_b·<A·ϵ>_b·∫_bin i f(M_X) dM_X  . 
        
        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
        
        L(σ_a,σ_b) = ∏_i(b_i+σ_a s_i)^b_i + σ_b s_i/(b_i + σ_b s_i)! e^-(b_i+σ_a s_i)
        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 = -2lnL(0,σ_sig)/L(σ_sig,σ_sig)  .
        
        Following <cit.> the significance Z_disc is then
        
        Z_disc=√(t_disc)  ,
        
        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
        
        t_excl = -2lnL(σ_sig,0)/L(0,0)  .
        
        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.
        
        
        
        
        99
        ATL-CONF-2014-005
          [ATLAS Collaboration],
          “A search for resonant Higgs-pair production in the
            bb̅bb̅ final state in pp collisions at √(s)=8 TeV”, 
             ATL-CONF-2014-005
        
        Aad:2014yja
          G. Aad  et al. [ATLAS Collaboration],
          “Search For Higgs Boson Pair Production in the γγ bb̅
           Final State using pp Collision Data at √(s)=8 TeV from the ATLAS Detector,”
          Phys. Rev. Lett.  114 (2015) no.8,  081802
         
          [arXiv:1406.5053 [hep-ex]].
        
        Aad:2015uka
          G. Aad  et al. [ATLAS Collaboration],
          “Search for Higgs boson pair production in the bb̅bb̅ final state from pp collisions
            at √(s) = 8 TeV with the ATLAS detector,”
          Eur. Phys. J. C  75 (2015) no.9,  412
          
          [arXiv:1506.00285 [hep-ex]].
        
        Aad:2015xja
          G. Aad  et al. [ATLAS Collaboration],
          “Searches for Higgs boson pair production in the hh→ bbττ, γγ WW^*,
           γγ bb, bbbb channels with the ATLAS detector,”
          Phys. Rev. D  92 (2015) 092004
          
          [arXiv:1509.04670 [hep-ex]].
        
        
        TheATLAScollaboration:2016ibb 
          The ATLAS collaboration,
          “Search for Higgs boson pair production in the bb̅γγ final state
            using pp collision data at √(s)=13 TeV with the ATLAS detector,”
          ATLAS-CONF-2016-004.
        
        ATL-CONF-2016-049
          [ATLAS Collaboration],
          “Search for pair production of Higgs bosons in the bb̅bb̅ final state using 
           proton–proton
            collisions at √(s) = 13 TeV with the ATLAS detector
             ATL-CONF-2016-049
        
        
        Khachatryan:2016sey 
          V. Khachatryan  et al. [CMS Collaboration],
          “Search for two Higgs bosons in final states containing two photons and
           two bottom quarks in proton-proton collisions at 8 TeV,”
          Phys. Rev. D  94, no. 5, 052012 (2016)
        
          [arXiv:1603.06896 [hep-ex]].
        
        CMS:2014eda
          CMS Collaboration [CMS Collaboration],
          “Search for di-Higgs resonances decaying to 4 bottom quarks,”
          CMS-PAS-HIG-14-013.
        
        Khachatryan:2014jya
          V. Khachatryan  et al. [CMS Collaboration],
          “Searches for heavy Higgs bosons in two-Higgs-doublet models and for t→ ch
           decay using multilepton and diphoton final states in pp collisions at 8 TeV,”
          Phys. Rev. D  90 (2014) 112013
          
          [arXiv:1410.2751 [hep-ex]].
        
        Khachatryan:2015yea
          V. Khachatryan  et al. [CMS Collaboration],
          “Search for resonant pair production of Higgs bosons decaying to two bottom
           quark-antiquark pairs in proton-proton collisions at 8 TeV,”
          Phys. Lett. B  749 (2015) 560
         
          [arXiv:1503.04114 [hep-ex]].
        
        Khachatryan:2015tha
          V. Khachatryan  et al. [CMS Collaboration],
          “Searches for a heavy scalar boson H decaying to a pair of 125 GeV Higgs 
           bosons hh or for a heavy pseudoscalar boson A decaying to Zh,
           in the final states with h →ττ,”
          Phys. Lett. B  755 (2016) 217
        
          [arXiv:1510.01181 [hep-ex]].
        
        Khachatryan:2016cfa
          V. Khachatryan  et al. [CMS Collaboration],
          “Search for heavy resonances decaying to two Higgs bosons in final states containing four b quarks,”
          Eur. Phys. J. C  76 (2016) no.7,  371
          
          [arXiv:1602.08762 [hep-ex]].
        
        Sirunyan:2017tqo
          A. M. Sirunyan  et al. [CMS Collaboration],
          “A search for Higgs boson pair production in the bbττ final state in proton-proton collisions at √(s) = 8 TeV,”
          arXiv:1707.00350 [hep-ex].
        
        
        
        
        CMS:2016tlj
          CMS Collaboration [CMS Collaboration],
          “Search for resonant pair production of Higgs bosons decaying to two bottom quark-antiquark
          pairs in proton-proton collisions at 13 TeV,”
          CMS-PAS-HIG-16-002.
        
        CMS:2016pwo 
          CMS Collaboration [CMS Collaboration],
          “Search for heavy resonances decaying to a pair of Higgs bosons in four b quark final state
           in proton-proton collisions at √(s) = 13 TeV,”
          CMS-PAS-B2G-16-008.
        
        
        CMS-PAS-HIG-16-013
          CMS Collaboration [CMS Collaboration],
          “Search for resonant Higgs boson pair production in the
                               bbτ^+τ^- final state,”
          CMS-PAS-HIG-16-013.
        
        CMS-PAS-HIG-16-029
          CMS Collaboration [CMS Collaboration],
          “Search for resonant Higgs boson pair production in the bb̅τ^+τ^-
            final state using 2016 data ,”
          CMS-PAS-HIG-16-029.
        
        CMS-PAS-HIG-16-032
          CMS Collaboration [CMS Collaboration],
          “Search for H(bb̅)H(γγ) decays at √(s)= 13 TeV,”
          CMS-PAS-HIG-16-032.
        
        Sirunyan:2017djm
          A. M. Sirunyan  et al. [CMS Collaboration],
          “Search for Higgs boson pair production in events with two bottom quarks and two τ
            leptons in proton-proton collisions at √(s)= 13 TeV,”
          arXiv:1707.02909 [hep-ex].
        
        CMS-PAS-HIG-17-006
          CMS Collaboration [CMS Collaboration],
          “Search for resonant and non-resonant Higgs boson pair production in the bb̅ℓνℓν
            final state at √(s)= 13 TeV,”
          CMS-PAS-HIG-17-006.
        
        CMS-PAS-HIG-17-008
          CMS Collaboration [CMS Collaboration],
          “Search for Higgs boson pair production in the final state containingtwo photons and two bottom quarks in proton-proton collisions at √(s)= 13 TeV,”
          CMS-PAS-HIG-17-008.
        
        Djouadi:2005gj
          A. Djouadi,
          
          Phys. Rept.  459 (2008) 1
        
          [hep-ph/0503173].
        
        Djouadi:2013vqa 
          A. Djouadi and J. Quevillon,
          
          
          JHEP  1310, 028 (2013)
         
          [arXiv:1304.1787 [hep-ph]].
        
        Djouadi:2015jea
          A. Djouadi, L. Maiani, A. Polosa, J. Quevillon and V. Riquer,
          
          JHEP  1506 (2015) 168
          [arXiv:1502.05653 [hep-ph]].
        
        Maniatis:2009re
          M. Maniatis,
          Int. J. Mod. Phys.  A  25 (2010) 3505
          [arXiv:0906.0777 [hep-ph]].
        
        Ellwanger:2009dp
          U. Ellwanger, C. Hugonie and A. M. Teixeira,
          
          Phys. Rept.  496 (2010) 1
          [arXiv:0910.1785 [hep-ph]]. 
        
        Kang:2013rj
          Z. Kang, J. Li, T. Li, D. Liu and J. Shu,
          
          Phys. Rev. D  88 (2013) no.1,  015006
        
          [arXiv:1301.0453 [hep-ph]].
        
        King:2014xwa
          S. F. King, M. Mühlleitner, R. Nevzorov and K. Walz,
          
          Phys. Rev. D  90 (2014) no.9,  095014
        
          [arXiv:1408.1120 [hep-ph]].
        
        Carena:2015moc
          M. Carena, H. E. Haber, I. Low, N. R. Shah and C. E. M. Wagner,
          
          Phys. Rev. D  93 (2016) no.3,  035013
        
          [arXiv:1510.09137 [hep-ph]].
        
        Ellwanger:2015uaz
          U. Ellwanger and M. Rodriguez-Vazquez,
          
          JHEP  1602 (2016) 096
        
          [arXiv:1512.04281 [hep-ph]].
        
        Costa:2015llh
          R. Costa, M. Mühlleitner, M. O. P. Sampaio and R. Santos,
          
          JHEP  1606 (2016) 034
        
          [arXiv:1512.05355 [hep-ph]].
        
        Baum:2017gbj
          S. Baum, K. Freese, N. R. Shah and B. Shakya,
          
          Phys. Rev. D  95 (2017) no.11,  115036
        
          [arXiv:1703.07800 [hep-ph]].
        
        vonBuddenbrock:2016rmr
          S. von Buddenbrock  et al.,
          
          Eur. Phys. J. C  76 (2016) no.10,  580
          doi:10.1140/epjc/s10052-016-4435-8
          [arXiv:1606.01674 [hep-ph]].
        
        Ellwanger:2004xm
          U. Ellwanger, J. F. Gunion and C. Hugonie,
          
          
          JHEP  0502 (2005) 066
          [arXiv:hep-ph/0406215].
        
        Ellwanger:2005dv
          U. Ellwanger and C. Hugonie,
          
          
          Comput. Phys. Commun.  175 (2006) 290
          [arXiv:hep-ph/0508022].
        
        Degrassi:2009yq
          G. Degrassi and P. Slavich,
          
          
          Nucl. Phys.  B  825 (2010) 119
          [arXiv:0907.4682 [hep-ph]].
        
        BSMAt13TeV
        Alwall:2014hca
          J. Alwall  et al.,
          
           
          JHEP  1407 (2014) 079
          
          [arXiv:1405.0301 [hep-ph]].
        
        Harlander:2012pb
          R. V. Harlander, S. Liebler and H. Mantler,
          
           
          Comput. Phys. Commun.  184 (2013) 1605
          
          [arXiv:1212.3249 [hep-ph]].
        
        
        Liebler:2015bka
          S. Liebler,
          
          Eur. Phys. J. C  75 (2015) no.5,  210
          
          [arXiv:1502.07972 [hep-ph]].
        
        Mantler:2015vba
          H. Mantler and M. Wiesemann,
          
          Eur. Phys. J. C  75 (2015) no.6,  257
          
          [arXiv:1504.06625 [hep-ph]].
        
        Ball:2013hta
          R. D. Ball  et al. [NNPDF Collaboration],
          
          Nucl. Phys. B  877 (2013) 290
        
          [arXiv:1308.0598 [hep-ph]].
        
        Sjostrand:2006za
          T. Sjostrand, S. Mrenna and P. Z. Skands,
          
          JHEP  0605 (2006) 026
          
          [hep-ph/0603175].
        
        deFavereau:2013fsa
          J. de Favereau  et al. [DELPHES 3 Collaboration],
          
          JHEP  1402 (2014) 057
          [arXiv:1307.6346 [hep-ex]].
        
        hep-ph/0512210
          M. Cacciari and G. P. Salam,
          
          Phys. Lett. B  641 (2006) 57
          [hep-ph/0512210].
        
        ATL-PHYS-PUB-2015-022 ATLAS Collaboration,
          “Expected performance of the ATLAS b-tagging algorithms in Run-2”
        Cacciari:2011hy
          M. Cacciari, M. Czakon, M. Mangano, A. Mitov and P. Nason,
          
          Phys. Lett. B  710 (2012) 612
        
          [arXiv:1111.5869 [hep-ph]].
        Czakon:2011xx
          M. Czakon and A. Mitov,
          
          Comput. Phys. Commun.  185 (2014) 2930
         
          [arXiv:1112.5675 [hep-ph]].
        Baernreuther:2012ws
          P. Baernreuther, M. Czakon and A. Mitov,
          
          Phys. Rev. Lett.  109 (2012) 132001
         
          [arXiv:1204.5201 [hep-ph]].
        Czakon:2012zr
          M. Czakon and A. Mitov,
          
          JHEP  1212 (2012) 054
        
          [arXiv:1207.0236 [hep-ph]].
        Czakon:2012pz
          M. Czakon and A. Mitov,
          
          JHEP  1301 (2013) 080
        
          [arXiv:1210.6832 [hep-ph]].
        Czakon:2013goa
          M. Czakon, P. Fiedler and A. Mitov,
          
          Phys. Rev. Lett.  110 (2013) 252004
          [arXiv:1303.6254 [hep-ph]].
        
        Cowan:2010js
          G. Cowan, K. Cranmer, E. Gross and O. Vitells,
          
          Eur. Phys. J. C  71 (2011) 1554
           Erratum: [Eur. Phys. J. C  73 (2013) 2501]
        
          [arXiv:1007.1727 [physics.data-an]]. ]
