Search strategies for pair production of heavy Higgs bosons decaying invisibly at the LHC

The search for heavy Higgs bosons at the LHC represents an intense experimental program, carried out by the ATLAS and CMS collaborations, which includes the hunt for invisible Higgs decays and dark matter candidates. No significant deviations from the SM backgrounds have been observed in any of these searches, imposing significant constraints on the parameter space of different new physics models with an extended Higgs sector. Here we discuss an alternative search strategy for heavy Higgs bosons decaying invisibly at the LHC, focusing on the pair production of a heavy scalar $H$ together with a pseudoscalar $A$, through the production mode $q \bar q \to Z^\ast \to HA$. We identify as the most promising signal the final state made up of $4b + E_T^\text{miss}$, coming from the heavy scalar decay mode $H \to hh \to b \bar b b \bar b$, with $h$ being the discovered SM-like Higgs boson with $m_h=125$ GeV, together with the invisible channel of the pseudoscalar. We work within the context of simplified MSSM scenarios that contain quite heavy sfermions of most types with ${\cal O}(10)$ TeV masses, while the stops are heavy enough to reproduce the 125 GeV mass for the lightest SM-like Higgs boson. By contrast, the gauginos/higgsinos and the heavy MSSM Higgs bosons have masses near the EW scale. Our search strategies, for a LHC center-of-mass energy of $\sqrt{s} =$ 14 TeV, allow us to obtain statistical significances of the signal over the SM backgrounds with values up to $\sim$ 1.6$\sigma$ and $\sim$ 3$\sigma$, for total integrated luminosities of 300 fb$^{-1}$ and 1000 fb$^{-1}$, respectively.


Introduction
The discovery at the LHC of a SM-like Higgs boson [1,2], with the most recent measurement of its mass set at m h SM = 125.09 ± 0.21 (stat.) ±0.11 (syst.) GeV [3], has been a major confirmation of the standard model of particle physics (SM). The experiments reveal that the Higgs boson mass value agrees quite well with the range preferred by the analysis of electroweak precision tests [4,5], and the spin, parity, and some of its couplings to the SM particles do not show deviations from the standard expectations [6]. Further studies of the Higgs couplings are required in order to test more precisely its SM nature [7,8], or to find evidence of physics beyond the SM (BSM). In fact, the LHC has already provided important bounds on the scale of new physics [9].
On the other hand, the presence of invisible Higgs decays, which would be another clear signal of new physics, is very well-motivated and predicted in many extensions of the SM, such as the MSSM (for a review, see, e.g., [59]). The searches for these exotic decays at the LHC are centered on the production of an invisibly decaying Higgs boson via gluon fusion [60], vector boson fusion [60][61][62][63], and in association with a vector boson [60,61,[63][64][65][66]. These searches are also carried out considering the production of dark matter candidates in association with a vector boson [66,67], with a Higgs boson [68][69][70][71][72][73], or with heavy-flavor quarks [74,75]. As it happens with the searches for heavy Higgs bosons, no significant excesses have been observed over the SM backgrounds in these invisible decay searches and limits are placed on the parameter spaces of the different models, production cross sections, and invisible branching ratios.
It is worth to stress that none of the searches listed above considers the heavy Higgs boson pair production. With the aim to explore this production mechanism, we analyze in a jointed framework the invisible heavy Higgs decays together with the decays into SM particles as final products. In particular, this combination represents a probe of the coupling between two heavy Higgs bosons and the dark matter particles taking part of the hard interaction. Remarkably, besides the gluon fusion plus jets channel, the heavy Higgs-pair production is unique to study invisible decays of the pseudoscalar. Without the vector boson fusion and bremsstrahlung modes, which are absent for the pseudoscalar, the other possible production mechanisms are the associated with a light Higgs and with a pair of quarks. However, the former is dynamically suppressed within the MSSM whereas the latter has a cross section at the edge of accessibility at HL-LHC. We propose then an alternative search strategy for the heavy neutral Higgs bosons, by means of the production of a pair of Higgs bosons H + A through the tree-level mode qq → Z * → H A, and taking the possibility of invisible Higgs decays into account. In order to make quantitative statements, we work within a particular MSSM scenario, so-called Slim SUSY [76,77], but the conclusions are general for any given scenario with a similar mass spectrum. In the Slim SUSY scenario the only new particles at the electroweak (EW) scale are the heavy Higgs bosons and the charginos and/or the neutralinos. Within the context of Slim SUSY, the R parity is conserved and therefore the lightest supersymmetric particle (LSP), which is assumed to be the lightest neutralino χ , is stable, which allows for invisible Higgs decays H, A → χχ if M H,A > 2M χ . In the low-tanβ regime, these invisible decay channels can be sizeable, and even the dominant decay mode of the pseudoscalar, as we will show below. This low-tanβ regime [78], recently noticed in [76,77] and highlighted in [54][55][56][57][58]79,80], gives rise to a rich Higgs phenomenology since the decay rates of the heavy Higgs bosons are not dominated by the decay channels into bottomquark and τ -lepton pairs, as occurs for moderate and large values of tan β. This situation opens the possibility of heavy Higgs boson searches through interesting channels, as H → W + W − , H → ZZ, H → hh, A → hZ, the decay modes into top-quark pairs H, A → tt , and the invisible decays H, A → χχ.
The paper is organized as follows: in Section 2 we present the MSSM scenarios in the low-tan β regime with invisible Higgs decays that give rise to the mentioned rich Higgs phenomenology. Section 3 is devoted to the heavy Higgs-pair production modes at the LHC, the description of the different potential final states, and the confrontation with the current LHC searches. In Section 4 a dedicated search strategy for the promising channel pp → H A → 4b + E miss T is performed, showing the expected significance of the signal at the LHC and the future prospects at the HL-LHC. Finally, perspectives and conclusions are presented in Section 5.

MSSM scenarios in the low-tan β regime with invisible Higgs decays
The rich phenomenology that arises in the MSSM Higgs sector, when the value of tan β is low, was recently remarked in [76,77] and emphasized in [54][55][56][57][58]79,80]. In the low-tan β regime, the Higgs phenomenology is significantly enriched due to the suppression of the heavy neutral MSSM Higgs bosons decays into down-type fermions, which are proportional to tan β. In such a case, these heavy Higgs bosons will have sizeable branching ratios for a variety of interesting decay modes, depending on their masses one could have: H → W + W − , ZZ, hh, tt , and A → hZ, tt . Another possibility, which was not considered in [54][55][56][57][58]79,80], is that the heavy scalar H and the pseudoscalar A decay invisibly into a pair of lightest neutralinos, H, A → χχ (if M χ < M H,A /2), assumed to be the LSP and consequently stable. These invisible decays can also have important branching ratios in the low-tan β regime and they could even be the dominant decay channel for the pseudoscalar.
In order to study this rich Higgs phenomenology, in the low-tan β regime with invisible Higgs decays, we work within the context of Slim SUSY scenarios [76,77], defined with the following assumptions: This class of simplified MSSM scenarios (which we can label as Slim Pheno) can be described with only four parameters at the EW scale: tan β, the pseudoscalar mass M A , the bino mass M 1 , and the higgsino mass μ. 1 We choose μ to be close enough to M 1 in order to obtain an important bino/higgsino admixture, which allows to have sizable Higgs-neutralino-neutralino couplings [84], but large enough to avoid the heavy Higgs decays into a pair of second neutralinos (χ 2 ) or higgsino-like charginos (χ ± ). Therefore, the only new particles present at low energies, relevant for the Higgs phenomenology that we are interested in, are the two heavy neutral Higgs bosons (H and A) and the LSP neutralino, under the condition M χ < M H,A /2. All along this work, we fix tan β = 3 as a reference value, and define three different benchmarks accordingly to the value of M A , which generate the following mass spectra, computed with the SUSY-HIT package [85]: The values of M A chosen here might be in tension with flavor physics observables as B → X s γ , due to the important contributions to the b → sγ transition coming from the charged Higgs bosons [86]. However, in the MSSM the different Higgs, chargino, and gluino contributions to B → X s γ are generically competitive and it is not difficult to obtain cancellation patterns among them [87]. In addition, if the correction factor 1 b to the Higgs-bottom Yukawa coupling is positive, the B → X s γ rate predicted in the MSSM is smaller than in the Type-II 2HDM [88], as considered in [86]. Although a detailed study of the constraints imposed by this class of low energy observables is beyond the scope of this work, we have checked that the predictions for BR(B → X s γ ) in our benchmarks are allowed at the 2σ uncertainty.
We show in Table 1 the values of the branching ratios of the main decay channels of the heavy neutral MSSM Higgs bosons, H and A, in these three simplified scenarios, calculated also with SUSY-HIT [85]. For M A = 200 GeV, the decay modes H, A → tt , H → hh, and A → hZ are kinematically closed, and consequently the dominant decay of the heavy scalar H is into W + W − (59%), followed by the channels into ZZ (23%) and bb (12%), while the dominant decay mode of the pseudoscalar is, interestingly, its invisible channel with a branching ratio of 62%. On the other hand, for M A = 300 GeV, the decay channel H → hh is open and becomes the dominant Table 1 Branching ratios of the main decay modes of the heavy neutral MSSM Higgs bosons in the low-tan β regime, for the three scenarios light-mass, moderate-mass, and heavy-mass (tan β = 3 and M A = 200, 300, and 400 GeV, respectively), computed with SUSY-HIT [85]. "C" and "F" stand for kinematically closed and forbidden, respectively. The values in parentheses correspond to the pseudoscalar Higgs boson A branching ratios.
200 300 400 one for the H Higgs boson, with a branching ratio of 45%. Whilst in this kinematic region the decay A → hZ is also accessible, with a decay rate of 22%, the invisible A decay mode keeps as the dominant one with a branching ratio of 46%. Finally, if the threshold for the production of a top-quark pair is open (which happens for M A = 400 GeV), the dominant decay channel of both H and A is into tt (77% and 91%, respectively). Notice that in this case the decay H → hh stays with a relatively significant branching ratio (11%), and the invisible decay of the pseudoscalar too (4%).

Heavy Higgs-pair production H + A at the LHC
As discussed in the previous sections, a relevant Higgs production mechanism at the LHC to accomplish our proposal is the Higgs-pair production H + A. This process occurs at tree-level and it is mediated by a virtual Z * coming from light-quark annihilation or via gluon fusion at one loop through a box diagram; both have been previously studied in [89] for the different pairs of MSSM Higgs bosons. In particular, it is worth to note a significant advantage for the tree-level channel qq → Z * → H A. Since its amplitude goes as sin(β − α), which is now constrained from the light Higgs search at the LHC to be of order one, there is not a dynamical suppression as in the case of qq → Z * → hA. In what follows, we will present the corresponding cross sections for the three benchmarks defined in the previous section, discussing the different final states at the LHC in each case. Next, we will probe these three scenarios confronting them with the general searches of new physics at the LHC at 8 TeV and 13 TeV.

Cross sections and final states
The cross sections for the heavy Higgs-pair production modes qq → H A and gg → H A at the LHC with √ s = 14 TeV are shown in Table 2, calculated at next-to-leading order (NLO) by means of the code HPAIR [90][91][92], for M A = 200, 300, and 400 GeV, with tan β = 3. The cross sections predicted for the qq mode are sizable, even for the heavy-mass scenario with M A = 400 GeV, of O(1) fb at least. In contrast, the values for the gg production mode are much smaller, between one and two orders of magnitude lower than qq → H A. 2 According to these numbers, it is clear that in the low-tan β regime, for the reference value of tan β = 3, it is a reliable approximation to consider only the mode qq → H A in order to calculate the H + A production cross section.
In what follows we shall discuss the cross sections of the main decay channels for the H + A pair production through the mode qq → H A. In order to compute the event rate for each decay mode, we use the branching ratios of H and A listed in Table 1 along with the H + A production cross sections obtained with HPAIR [90-92] (see Table 2). Thus, the cross sections are given by As we want that one of the two heavy Higgs bosons decays invisibly, we will consider only the invisible channel of the pseudoscalar, A → χχ, since A has the largest invisible rates, being in addition its dominant decay mode in the two first benchmarks. For those decay chains that include H → W + W − we consider only the semileptonic decay mode, W + W − → lν l + jj , with l = e, μ, so that BR(W W → lν l + jj ) = 0.29. In the case of the decay H → ZZ, we take into account the leptonic Z-boson decays, and thus BR(ZZ → 4l) ' 4 × 10 −4 , with 4l = 4e, 2e2μ, 4μ. For the processes with the H → hh decay mode, we will consider BR(h → bb) ' 0.7 and BR(h → τ + τ − ) ' 0.07. The resulting cross sections with √ s = 14 TeV are shown in Tables 3, 4, and 5 for the three scenarios light-mass, moderate-mass, and heavy-mass, respectively. Some comments are required in each case: • In the case of M A = 200 GeV (light-mass scenario), we see from Table 3 that the dominant decay mode is the final state W + W − (→ lν l jj 0 ) + χχ, being the subdominant decay mode bb + χχ. The former has the advantage of a very large cross section, but the presence of E miss T in one of the W decay would result in an involved analysis since it affects the reconstruction of the missing transverse energy profile of the invisible A decay, possibly preventing an efficient discrimination between signal and background. Moreover, the hadronic W decay would need to be confronted with large QCD backgrounds. On the other hand, the latter seems to be more interesting, but a detailed study of this bb + E miss T channel [93] has demonstrated that it will be challenging to probe it, even during the HL-LHC run. Instead of the H decay into a bottom-quark pair, the decay into a τ -lepton pair is more efficient in reducing the QCD background, but it is affected by its relatively small cross section which is one order of magnitude lower than in bb channel. Moreover, the τ -tagging efficiency is not much better than the b-tagging one. Finally, the cleanest channel would be 4l + E miss T , where the four leptons come from the heavy scalar decay H → ZZ with both Z bosons decaying leptonically. Unfortunately, this channel has a very tiny cross section and it would be a difficult challenge to enhance the signal significance over the SM backgrounds. • The final states allowed for the moderate-mass scenario seem to be more promising, as we can see in Table 4, since the H → hh 3 decay is open. Indeed, the dominant channel is 4b + E miss T , in which the four bottom quarks are the decay products of two SM-like Higgs bosons h originated from this decay. This class of final state, with another different decay chain, is also commented in [93] as a very promising channel to probe new physics at the LHC, but was not analyzed in detail. We could also take advantage of the presence of the decay H → hh and consider the channel 2b2τ + E miss T . Nevertheless, the cross section in this case is one order of magnitude lower than the 4b + E miss T final state, and we would have to deal again with large QCD backgrounds and b and τ taggings anyway. In this kinematic region, the channels W + W − + E miss T and bb + E miss T still have important cross sections, but lower than the 4b + E miss T final state. However, these decay channels have the limitations already discussed in the light-mass scenario.
• The kinematic region in which the threshold of tt production in the heavy Higgs decays is open (heavy-mass scenario) is dominated by far by the 4t final state (see Table 5), with a cross section of ∼1 fb. The largest channel with missing transverse energy is tt + E miss T , with a cross section 20 times smaller than the 4t channel. With such a low cross section the QCD background would overwhelm the signal in the case of hadronic top decays; whereas for leptonic top decays, the E miss T coming from the tops would distort the missing transverse energy distribution from the invisible A decay. Therefore, it seems rather complicated to probe the invisible Higgs decays through the heavy Higgs-pair production in this kinematic region. However, one can just study the heavy Higgs-pair production through the four tops channel, which could be detectable at the LHC with the search strategies designed in [99]. For this scenario, the branching ratio of H → hh keeps still a sizable value (11%). Thus, another channel to try to probe the heavy Higgs-pair production could be 4b + 2t. Regrettably, the cross section is more than one order of magnitude lower than the 4t channel and the QCD background is also very large to obtain a substantial significance.
Taking into account the arguments stated above, we shall focus on the moderate-mass scenario, with the channel qq → H A → 4b + E miss T (H → hh → 4b and A → χχ) as the most promising signal to probe invisible Higgs decays through heavy Higgs-pair production at the LHC. This same final state has been recently analyzed by CMS [100] within the context of gauge-mediated SUSY breaking, considering the electroweak production of two higgsinos that decay into the lightest Higgs boson h and the LSP goldstino G . The results reported by CMS are consistent with the SM background predictions and 95% CL exclusion limits are imposed on the higgsino-pair production cross sections, which are much larger than ours. Moreover, this 4b + E miss T final state, originated from this SUSY cascade process, has a very different kinematic behavior from the final state coming from our proposed signal pp → H A → 4b + E miss T . In the next section we will develop a dedicated search strategy for this channel. Before that, we show below that the recent LHC searches do not exclude a signal compatible with the production process pp → Z * → H A in any of the proposed benchmarks.

LHC searches
In order to probe the three scenarios introduced in Table 1, we have confronted them with the general searches of new physics at the LHC at 8 TeV and 13 TeV by using the software Check-MATE 2 [101]. We have simulated the process pp → Z * → H A using MadGraph 5 [102] and decayed the heavy neutral Higgs bosons through the branching ratios computed with the SUSY-HIT package [85]. The cross section of the H A production was computed at NLO by using HPAIR [90][91][92], which includes QCD and SUSY-QCD corrections taken from [103] and [104], respectively. Finally, we have implemented PYTHIA [105] for the parton shower and hadronization, and the detector simulation was carried out using Delphes 3 [106].
For both the 8 TeV and the 13 TeV searches, none of the three scenarios is excluded by the CheckMATE 2 validated analyses. For each analysis implemented in CheckMATE 2, we have also computed the signal significance of the most restrictive signal region. We have used the approximate formula S/ √ B, where S is the number of signal events in a given signal region and B is the corresponding number of background events. The results derived from the LHC searches at 8 TeV included in CheckMATE 2 are listed in Table 6. In the last column, we include the name of the most restrictive signal region as given in the corresponding experimental analyses. For the three scenarios, the values of S/ √ B obtained with the signal regions of the general new physics LHC searches at 8 TeV are low. The results corresponding to the LHC searches at 13 TeV are also shown in Table 6 in parentheses. As in the case of the 8 TeV searches, the values obtained for the significance are also low for the three scenarios.
For M A = 200 GeV, the most restrictive signal regions are S7 and bCbv for the LHC searches at 8 and 13 TeV, respectively. The former is defined in the analysis associated to Table 6 Significances corresponding to the LHC at 8 TeV and 13 TeV (in parentheses) obtained from CheckMATE 2, for the three benchmarks with M A = 200, 300, and 400 GeV. In the last column the name of the most restrictive signal region is included. the ATLAS search for direct top-squark pair production in final states with two leptons of opposite charge using 20.3 fb −1 [107], while the latter corresponds to the ATLAS search for top squarks in final states with one isolated lepton, jets, and missing transverse momentum using 13.2 fb −1 [108]. The signal region S7 is defined through the lepton-based stransverse mass, m T 2 , requiring m T 2 > 120 GeV and also in such a way that the number of jets with p T > 20 GeV must be smaller than 2. On the other hand, the signal region bCbv selects events with two or more jets with p T > (120 GeV, 80 GeV) and no b-tagged jets. The following requirements are also imposed: E miss Finally, for M A = 400 GeV the signal region that exhibits the highest significance is called SR3b for both cases of 8 TeV and 13 TeV searches, although its definition is different for each analysis. In both cases, the search is focused on final states with jets and two same sign (SS) leptons or three leptons. For the 8 TeV analysis [110], this signal region requires two SS leptons or three leptons with at least five jets and at least three b-jets, and also m eff > 350 GeV, while for the 13 TeV analysis [111] the requirements becomes: N lept ≥ 2, at least three b-jets with p T > 20 GeV, E miss T > 125 GeV, and m eff > 650 GeV. The analyses considered above show that even the most restrictive signal regions corresponding to general searches are far from being sensitive to the heavy Higgs-pair production. It is then sensible to design a dedicated search strategy for this process. We present in the next section a detailed discussion on this matter. requiring exactly 4 b-tagged jets. However, since in our working point the maximum b-tagging efficiency is ∼ 75%, we have also studied a second signal region, SR2, in which the above b-tagging requirement is relaxed to 3 or 4 b-tagged jets. In addition, both signal regions require 0 or 1 light-jet to suppress multi-jet backgrounds, and include a lepton veto to disfavor the presence of missing transverse energy due to neutrinos in backgrounds involving top quarks. In what follows, we present first the general features of the signal and a procedure to optimize its potential detection. Later on, we will detail the search strategies for both signal regions.

General signal features and cut optimization procedure
Let us start remembering that the proposed signal final state is made up of 4 b-jets and it has large E miss T , having a total cross section of 0.453 fb, and stating that the coming discussion on the procedure used to maximize the sensitivity of the signal will be centered on the irreducible backgrounds, postponing our treatment of the reducible contributions to Sections 4.2 and 4.3. The motivation is that the former expectedly follow the signal more closely than the latter. Thus, an optimization of the signal significance based on the irreducible backgrounds may be sufficient to achieve an efficient discrimination among the signal and all background sources. In other words, we first develop the search strategy by taking into account the signal and only the irreducible backgrounds, and subsequently, we estimate the reducible contribution. The main irreducible SM backgrounds for this final state, in both signal regions SR1 and SR2, are ttbb, Z(→ νν)bbbb, and Z(→ νν)bbbbj . Both the signal and the main irreducible backgrounds have been generated with MadGraph_aMC@NLO [102] and showered with PYTHIA [105], simulating the detector response with Delphes 3 [106]. The events have been generated within the four-flavor scheme with the factorization and renormalization scales set to μ R = μ F = H T /2 and using the PDF MSTW2008nlo68cl_nf4. Besides, the following basic cuts have been imposed at generator level 4 For the signal and each one of the main irreducible backgrounds, the number of generated events has been much larger than the expected amounts for an integrated luminosity of ∼1500 fb −1 and a center-of-mass energy of √ s = 14 TeV. 5 Let us remember that the signal cross section has been obtained at NLO by means of the codes HPAIR and SUSY-HIT (see Section 3), whilst for the backgrounds we have used the NLO cross section calculated with MadGraph_aMC@NLO only in the case of ttbb. The cross sections obtained with MadGraph_aMC@NLO for the three main irreducible backgrounds are given by σ (ttbb) = 1633 fb, σ (Z(→ νν)bbbb) = 3.27 fb, σ (Z(→ νν)bbbbj ) = 2.45 fb . 4 These cuts are consistent with the usual selection performed by the trigger system and also by the reconstruction algorithms in the experimental analyses. Since we have checked that the distortions on these kinematic variables due to the parton shower are statistically negligible, we apply these cuts at the generator level, which makes the simulation process more efficient. Finally, note that the same cuts are applied both to light and b-jets, which prevents the impact of potential misidentification at detector level on the corresponding kinematic distributions. 5 Due to the rather large cross section, the number of generated events for the ttbb background has been five times larger than the amount of events expected at 300 fb −1 . After the event generation with the basic cuts, the next step on the analysis will be to require that both the signal and background events pass a set of selection cuts at detector level, whose definition depends on the signal region: The jet multiplicity of the backgrounds is expected to be much larger than the signal one, as it can be observed in Fig. 1, in which we display the distribution of the number of light jets after the selection cuts, without the light-jet requirement, for the signal and the main irreducible backgrounds in the signal regions SR1 (left panel) and SR2 (right panel). From these two plots it is clear that a cut in the number of light jets, demanding 0 or 1 light jet at most, will keep a large part of the signal while removing a considerable portion of the backgrounds rates, especially for ttbb, which has the largest cross section by far.
Another interesting feature of the signal is the expected very large values of E miss T as compared with the backgrounds. Therefore, a potential powerful discriminator between the signal and backgrounds is the variable E miss q H nj T > 10 GeV 1/2 for the former, which points out that a cut requiring values of the E miss T significance larger than 10 GeV 1/2 would largely increase the signal-over-background ratio. Besides, in contrast to the E miss T alone, the E miss T significance leads to a better discrimination even among the different backgrounds allowing in turn a higher significance of the signal. In particular, a cut on E miss T significance will have a considerable impact on the ttbb background which after imposing the selection cuts is characterized by events with instrumental

E miss
T and, precisely, the E miss T significance is highly efficient to discriminate fake from genuine E miss T . Apart from the cuts described above, which are common for both signal regions SR1 and SR2, we have also considered a wide set of angular variables (more than 40 altogether), related to the azimuthal angle φ and the pseudorapidity η. These angular variables have been introduced with the aim to amplify the differences between the signal and background distributions, which still have similar kinematics after the imposition of: the selection cuts, the restrictions over the E miss T significance, and over the invariant mass variables (only in the case of SR1). As we will show in the next section, the use of some of the angular variables turns out to be useful to improve the significance of the signal.
The procedure to select the optimal values of the cuts for the different kinematic variables has been developed by means of a sequential optimization of the statistical significance, S, calculated according to the expression [112]: where S and B denote the number of signal and background events, respectively. The ultimate goal of the cut optimization is to determine a search strategy that maximizes the statistical significance in each signal region. The procedure starts with the selection of events that pass the discrete cuts on the number of b-jets, light-jets, and leptons. After that, the distributions of the kinematic variables (transverse momenta of the final state particles, E miss T significance, invariant masses, angular variables) are then scrutinized. Once the variable that shows a better discriminating power between signal and background is chosen, the optimization algorithm is executed just for this variable and the optimal value of the cut is applied to the events. This procedure is then repeated by picking and optimizing the most discriminating variable at each step until the maximum statistical significance is achieved, keeping at least one signal event. Finally, the acceptances for the signal and backgrounds, defined as the proportion of events at the detector level that pass the cuts implemented at each step with respect to the amount of events generated with the basic cuts, are used to apply correction factors to the expected number of events at a certain luminosity.

SR1: 4 b-jets signal region
In this case, besides the main irreducible backgrounds ttbb, Z(→ νν)bbbb, and Z(→ νν)bbbbj mentioned above, we have to deal with the major reducible backgrounds that arise from tt , tt + jets and Z + jets. Specifically, we have considered 6 tt, ttj, ttjj , Zbbjj, Zjjjj , Zbbjjj, Zjjjjj .  More precisely, we restrict the search to the region 250 GeV < m 4b < 340 GeV, where m 4b is the invariant mass of the four b-jets, and also make use of the variable χ hh [29] defined as: requiring that over all possible combinations of b 1 -b 4 at least one of them gives a χ hh value below 2.7. This last cut on χ hh is intended to select events in which at least two disjoint pairs of b-jets are most likely to arise from the decay h → bb. (d) We perform a cut on the azimuthal angular separation between a b-jet and the missing transverse energy. Among all the possible values, we require the second minimum between any b-jet and E miss T in the event (1φ 2 nd min bE miss T ) to be above 1.6. 6 There is an additional source of reducible background that we have not included in our analysis arising from multi-jet production. The cross section for this background is huge and a data driven approach is necessary. However, since in this case there is not a genuine source of missing transverse energy but just a fake contribution as a consequence of imperfect reconstruction of objets and energy resolution, an E miss T significance ∼1 GeV 1/2 is expected for this background, much below the value of the cut applied to this variable, as we will see below. Table 7 Signal and background event cut flow corresponding to the first signal region (SR1) for √ s = 14 TeV and a total integrated luminosity of L = 300 fb −1 .

Process
Signal  Table 8 Signal and background event cut flow corresponding to the first signal region (SR1) for √ s = 14 TeV and a total integrated luminosity of L = 1500 fb −1 .

Process
Signal The cut flow for the signal and background events obtained by applying the steps described above is shown in Table 7 for a total integrated luminosity of L = 300 fb −1 . The cut flow is the result of the sequential optimization of the statistical significance, S, as introduced in Section 4.1. At each step we display the events corresponding to the signal and to the irreducible backgrounds along with the corresponding significance. Besides, we display in Table 8 the cut flow corresponding to a total integrated luminosity of 1500 fb −1 . In this case no new optimization procedure has been applied because the acceptances can be safely rescaled, since the amount of generated events is consistent with the required luminosity. In addition, since the Z + jets reducible backgrounds are not attainable through our simulations, we show in the sixth column an estimation of the number of events corresponding to these reducible backgrounds that remain after imposing the cuts involved in each step. For this estimation, we have used conservative misidentification rates for light jets, ² j = 5 ×10 −3 , and c-jets, ² c = 0.15, and a nominal b-tagging efficiency of 75%. For the steps (b) to (d), we have also assumed Z + jets to have the acceptances of Zbbbb or Zbbbbj (according to the number of jets in the final state), since similar kinematic distributions of events are potentially expected for this type of backgrounds. Finally, the tt + jets reducible backgrounds are also beyond our computational resources to simulate them. However, the acceptances of the tt + jets are smaller at each step of the cut flow than the acceptances of the corresponding irreducible background ttbb. 7 Therefore, it seems fairly sensible to assume that no tt + jets events would survive at the end of this search strategy. Table 9 Signal and background event cut flow corresponding to the second signal region (SR2) for √ s = 14 TeV and a total integrated luminosity of L = 300 fb −1 .

Process
Signal

SR2: 3 or 4 b-jets signal region
The main backgrounds in this signal region are the same as those involved in SR1, except for the addition of two new reducible backgrounds, Z(→ νν)bbj and Z(→ νν)jjj , with cross sections of 2536 fb and 7564 fb, respectively. By using the same set of simulated events as in the case of SR1, the search strategy for SR2 involves the following steps: (a) We select events with 3 or 4 b-jets, 0 leptons, 0 or 1 light-jet, and also require the transverse momentum of the third leading b-jet to be above 30 GeV.
T sums over the three leading jets. (c) We apply two additional cuts: 1φ max bb < 1.90 and 1η max bb < 1.42, where 1φ max bb is the maximum azimuthal angular separation between two b-jets and 1η max bb is the maximum difference in pseudorapidities between the b-jets in the event.
The cut flow of the signal and background events obtained by applying the steps (a), (b), and (c) are shown in Tables 9 and 10, again for total integrated luminosities of L = 300 fb −1 and 1500 fb −1 , respectively. Following the same approach as in SR1, Table 9 has been obtained through the sequential optimization of the statistical significance whereas Table 10 is a projection of the results at 300 fb −1 . We have dealt with the reducible backgrounds in a similar way as in SR1. They have been introduced into the cut flow through an estimation of their impact using the same misidentification rates for light jets (² j = 5 × 10 −3 ), c-jets (² c = 0.15), and the nominal b-tagging efficiency (75%). For the steps (b) to (d), we have also assumed the Z + jets backgrounds to have the same acceptances corresponding to Zbbbb or Zbbbbj . Finally, the tt + jets backgrounds have received the same treatment as in SR1 since analogous arguments are valid in this signal region as well. Table 11 Number of signal and background events along with the significances obtained without taking into account systematic uncertainties (see Eq. (1)) and assuming that this source of uncertainty amounts to 30% (see Eq. (2)). The results corresponding to luminosities of 300, 1000, and 1500 fb −1 are shown.

Discussion
From the comparison between the Tables 7 and 9 (300 fb −1 ), and also between the Tables 8 and 10 (1500 fb −1 ), we can see that the search strategy for the SR1 signal region tends to give a higher significance. However, notice that the signal region SR2 does not involve any cut in invariant mass, what makes this search strategy more model-independent than the SR1 one. Even so, the m 4b invariant mass window considered in SR1 is broad enough to cover a large region of Higgs masses and it contains the moderate mass scenario which is in fact the scenario that we attempt to probe with the search strategies presented here. In addition, the signal region SR2 has the advantage of retaining more signal events; indeed, we expect 4 signal events for L = 300 fb −1 , while just one signal event is expected when applying the SR1 strategy for the same luminosity. This increase in the signal events comes at the expense of keeping at the same time more background events and adding two more reducible backgrounds. Moreover, in SR2 the backgrounds are dominated by their reducible contribution. The presence of more background events is risky when the possible systematic uncertainties are taken into account; in fact, when such a source of uncertainties is included, the statistical significance defined in Eq. (1) is modified as follows [113]: where σ B = (1B)B, with 1B being the relative systematic uncertainty. From the comparison between Eqs. (1) and (2), it can be seen that the bigger the number of background events, the higher the degradation of the signal significance due to the systematic uncertainties.
In Table 11 we summarize the results for the signal regions SR1 and SR2 including the significance S sys , obtained assuming 30% of systematic uncertainties. Besides presenting the results corresponding to L = 300 and 1500 fb −1 , we also display the results for an intermediate luminosity of L = 1000 fb −1 . From the Table 11, it is clear that the systematic uncertainties have a meaningful impact on the significances corresponding to SR2, while the values of the significance are slightly reduced in the case of SR1. With an integrated luminosity of 1500 fb −1 , the SR1 retains 5 signal events with a significance of more than 3σ , even when a 30% of systematic uncertainties is considered. The results for this luminosity in the case of SR2 are similar, except that now the significance is degraded to ∼ 1.8σ due to the effect of the systematic uncertainties. The results for the intermediate luminosity (1000 fb −1 ) are also promising, achieving evidence of the new physics signal in both signal regions. Again, note that with the conservative value chosen for the systematic uncertainties (30%), the significance for SR2 decreases to ∼ 1.7. It is important to mention here that our choice of the signal cross sections is rather conservative and if we had considered an enhancement of the final cross section, we could have obtained better significances even for the lowest luminosity of 300 fb −1 . For example, an increase of a factor 2 in BR(H → hh) × BR(A → χχ), which is not harebrained, 8 would mean a proportional enhancement on the signal events, resulting in statistical significances of 2.73 and 2.96 at 300 fb −1 for SR1 and SR2, respectively. These results would not be so affected by the systematic uncertainties, with values for S sys of 2.65 and 2.27, respectively, which are also very close to the evidence threshold.
Finally, although we have not generated enough number of background events to provide a prediction of the significances at 3000 fb −1 , we have naively assumed that the background rates scale as the cross section of the signal [57], which is a conservative standpoint. Under this assumption, we obtain that the statistical significances S for SR1 and SR2 at 3000 fb −1 are 5.14 and 5.04, respectively, which reach the discovery level. On the other hand, if we assume conservatively that the 30% systematic uncertainties will persist at these luminosities, the significances are reduced to 4.40 for SR1 and 1.89 for SR2. It is clear that, even under these conservative assumptions, the search strategies developed along this work offer a chance to show a first hint at the HL-LHC of this class of new physics signals of heavy Higgs bosons decaying invisibly.

Conclusions
In this work we have developed a search strategy for invisible Higgs decays through heavy Higgs-pair production at the LHC, which makes plausible its detection. Our proposed alternative to the search strategy for heavy neutral Higgs bosons focuses on the production of a pair of Higgs bosons H + A through the tree-level mode qq → Z * → H A, and considers the elusive possibility of detecting invisible pseudoscalar decays. In order to make quantitative statements, we have worked within a particular MSSM scenario, so-called Slim SUSY, but the conclusions hold in general for any given scenario with a similar mass spectrum. The only new particles at the EW scale in these MSSM scenarios are the heavy Higgs bosons and the charginos and/or the neutralinos. This fact allows for invisible Higgs decays H, A → χχ if M H,A > 2M χ . Within the moderate-mass scenario (M A = 300 GeV) the channel qq → H A → 4b + E miss T (H → hh → 4b and A → χχ) is found to be the most promising signature to probe invisible Higgs decays through heavy Higgs-pair production at the LHC. We have shown that the invisible channel is the dominant decay mode of the pseudoscalar both for the light and moderate mass scenarios, so that the cross sections of those channels that involve invisible particles in the final state are higher than those that lead to visible final states. This, together with the possibility of using powerful discriminator variables based on E miss T significance, makes the search for heavy neutral scalars through their invisible decays more auspicious than considering visible final states. It is important to note that this is not the case for the large mass scenario, which is clearly dominated by the four tops channel, where both heavy scalars decay visibly (see Table 5).
We have defined two different signal regions, namely, SR1 with exactly 4 b-tagged jets and SR2, in which the b-tagging requirement is relaxed to 3 or 4 b-tagged jets. A detailed search strategy has been performed in each signal region, based on a sequential optimization of the corresponding cuts that maximizes the statistical significance in each case. As a result, we have obtained prospects of evidence (∼3σ ) at 1000 fb −1 and a conceivably possibility of discovery (∼5σ ) at 3000 fb −1 in both of the two signal regions. If we take into account systematic uncertainties of 30%, the evidence significance is degraded to 2.80 (1.67) for the signal region SR1 (SR2), whilst the discovery significance reduces to 4.40 (1.89). We can conclude that the search strategies developed along this work offer the opportunity to discover this class of new physics signals of heavy Higgs bosons decaying invisibly at the HL-LHC.
On the other hand, although a thoughtful analysis of the reducible backgrounds is beyond the scope of this work, our results are conservative. In principle, it should not be difficult to improve them. A better cut optimization through multivariate analysis (MVA) with the boosted decision tree (BDT) algorithm might ameliorate the signal/background ratio. Our estimation of the systematic uncertainties is very rough and one could expect a significant reduction in the future. In addition, although a combination of significances of both signal regions (SR1 & SR2) would imply dealing with large correlations and thus making it highly nontrivial, it might result in an increase of the total significance, with potential interest even for the lowest luminosity of 300 fb −1 considered here. Furthermore, if it were possible to keep the same value of the signal cross section for larger values of the heavy Higgs bosons masses (increasing, for instance, the rates of the decay channels H → hh and A → χχ), the E miss T in the signal events would be greater than in the case of the moderate-mass scenario. Therefore, a more stringent cut on the E miss T significance would lead to an improvement on the significances. Last but not least, larger signal cross sections could be generated within the MSSM scenarios considered in this work or in other BSM models that drive us to the same final state topology. For example, it is not harebrained to obtain an increase of a factor of 2 in BR(H → hh) × BR(A → χχ), which would mean a proportional enhancement on the signal events, resulting in statistical significances of 2.73 and 2.96 at 300 fb −1 for SR1 and SR2, respectively. Considering systematic uncertainties of 30%, these significances would be 2.65 and 2.27, respectively, which are also close to the evidence level.