Dark and Bright Signatures of Di-Higgs Production

If the Higgs boson decays to a pair of invisible particles, the number of di-Higgs events, where each Higgs decay into Standard Model (SM) particles, are reduced by a factor of two-third taking into account the current LHC bound on invisible decay width of the Higgs boson. We investigate the sensitivity of the upcoming high luminosity run of the LHC to di-Higgs production and subsequent decay to dark matter in the context of the singlet scalar extension of the SM augmented by a fermionic dark matter in the dark and bright channel $\gamma\gamma+\not\!\! E_T$. Once systematic uncertainties on background yields are considered, this dark and bright channel presents competitive limits than $b\bar{b}+\not\!\! E_T$ after a careful tuning of the kinematical cuts that raise the signal over background ratio. We further show that in a multivariate analysis, for an invisible branching fraction as low as $\sim 10$%, we obtain stronger bounds for the Higgs trilinear coupling from the $\gamma\gamma+\not\!\! E_T$ channel compared to the $b\bar{b}\gamma\gamma$ final state. Finally, we demonstrate that the three channels $\gamma\gamma+\not\!\! E_T$, $b\bar{b}+\not\!\! E_T$ and $b\bar{b}\gamma\gamma$, complement each other in the search for di-Higgs production with non-SM trilinear couplings when an invisible decay mode is present.


I. INTRODUCTION
After the discovery of the Higgs boson at the LHC [1,2], and the confirmation that this spin-0 particle plays a role in the electroweak symmetry breaking (EWSB) mechanism [3][4][5], the LHC still has a long road to unravel the details of the Higgs boson self-interactions by measuring its trilinear and quartic couplings. These measurements are key in the understanding of the stability of the Higgs potential and the nature of the electroweak phase transition (EWPT). The Standard Model (SM) Higgs potential, by itself, is metastable and cannot trigger a first order EWPT either. Additional scalar fields present in a beyond the SM model (BSM) can facilitate stabilizing the Higgs potential and also can trigger EWPT with profound implications for the matter-antimatter asymmetry. Thus, the stability of the vacuum and a 1 st order EWPT make the Higgs self-interactions possible targets for signals of BSM physics. At the HL-LHC, such signals might reveal themselves either as deviations of the trilinear Higgs coupling as compared to its SM value or new resonances.
Another physics goal of the upcoming runs of the LHC is to search for dark matter (DM) and its possible connection to the Higgs. The upcoming high luminosity run of the LHC (HL-LHC) has an enormous potential to discover or exclude DM models, including the Higgs portal models [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] where Higgs bosons couple to the DM field. Such interaction is expected if the DM acquires its mass via the EWSB mechanism.
The LHC experiments have searched for invisibly decaying Higgs bosons and the current most stringent bound for its branching ratio into an invisible dark state is provided by the CMS collaboration, which is 19% [21]. In contrast, in all the searches for di-Higgs by CMS and ATLAS, they do not consider any significant beyond the SM decay channel of the * Electronic address: aalves@unifesp.br † Electronic address: tghosh@hawaii.edu ‡ Electronic address: farinaldo.queiroz@iip.ufrn.br Higgs in addition to the SM modes, which is, of course, a reasonable assumption given the state of affairs. However, if the Higgs boson indeed possesses an invisible decay channel with branching ratio, BR inv , then decay rates of the Higgs to all the SM channels will be reduced by (1 − BR inv ). Hence, the di-Higgs production rate with both Higgs bosons decaying to the SM particles will go down by a factor of (1 − BR inv ) 2 , which is around 0.65 taking into account the current LHC bound. Thus, all the discovery prospects of di-Higgs will also drop nearly by the same factor. In this scenario, looking for channels where one of the Higgs bosons decaying invisibly becomes an important task. Recently, the hh → bb+ E T [22] channel was studied in the context of the SM and a singlet scalar extension of it, called xSM [23][24][25]. xSM is an economical model with interesting phenomenological consequences like detectable gravitational waves signals from strongly 1 st order EWPT, and a new massive scalar that decays to hh. Among the hh channels involving E T , this is the final state with the largest expected number of events but also with the largest backgrounds, mainly from bbZ and Zh. In the unlikely case where systematic effects can be neglected, this channel might reach around 4σ of statistical significance in the SM case using a Boosted Decision Trees (BDT) algorithm. However, its signal-tobackground ratio is tiny, and once even very small systematic backgrounds are taken into account, that significance drops below 1σ.
In this paper, we study the mode hh → γγ+ E T , a dark and bright decay of the di-Higgs system, and explore non-SM trilinear couplings in the context of xSM [23][24][25] augmented by a fermionic DM. If the DM is a scalar particle, the Z 2 symmetry required to stabilize the DM will preclude cubic terms in the scalar potential, which facilitate to trigger a 1 st order EWPT. However, in such models 1 st order EWPT can still be induced via zero-temperature loop effects, thermally driven scenarios, and modifications to the scalar potential at tree level by renormalizable or non-renormalizable operators (see [26] and references therein). On the other hand, for a vector DM one needs an extra scalar to generate its mass via spontaneous symmetry breaking [16][17][18][19][20]. A detailed study on these intricacies will be presented elsewhere.
We focus on the part of the xSM parameters space that respects several experimental and theoretical constraints, and where new heavy Higgs bosons are not expected to be observable at the LHC but shift the SM Higgs trilinear coupling away from its SM value. In this case, our results can be useful to constrain other models with non-SM trilinear couplings. We do not restrict ourselves to points that might deliver a strong gravitational wave signal in future space-based interferometers; however, the parameters space of interest in the present study is expected to cover that where a strongly 1 st order EWPT might occur [27].
We also compare this channel's performance with bb+ E T and bbγγ final states. Even though the expected number of signal events is much lower than bb+ E T , the backgrounds are also much smaller, and a more favorable S/B is helpful to tame the systematic uncertainties for similar signal significances. We show that by carefully tuning the kinematic cuts using a Gaussian Process algorithm, the photons plus missing energy channel has a competitive reach both for discovery as well as exclusion compared to the bb+ E T channel once realistic systematic uncertainties are included in the computation. Moreover, using BDT along with tuned cuts, we show that γγ+ E T becomes competitive with bbγγ in terms of signal significance, with invisible branching ratios close to the current LHC bound for non-SM trilinear couplings which make the hh cross sections large.
Besides, we want to emphasize that estimating the reach of γγ+ E T is also important in view of the combination with another channel like bb+ E T . Combining several search channels is often possible in experimental studies which take systematic effects and correlations into account. This is explicitly shown in the case of double Higgs production into SM channels in Ref. [28].
Our paper is organized as follows. In Section II, we introduce our xSM model assumptions. We impose the constraints on the xSM parameter space for our collider analysis in Section III. In Section IV, we lay out our strategy to analyze the γγ+ E T final state. The results of the cut-and-count and multivariate analyses on the γγ+ E T channel are presented in Sections V and VI, respectively. Finally, we conclude in Section VII.

II. MODEL ASSUMPTIONS
An invisibly decaying SM Higgs boson might arise in Higgs portal models where the Higgs couples either directly to scalar (S) or vector (V ) dark matter fields through renormalizable interactions like |S| 2 |H| 2 and V µ V µ |H| 2 , respectively. In the case of fermionic DM, renormalizable interactions with the Higgs field arises only by the addition of an extra scalar. When such additional scalar acquires a vacuum expectation value (vev), the DM mass generation through spontaneous symmetry breaking occurs.
For example, consider the simple extension where one real gauge singlet is added to the Higgs potential, the so-called xSM [23][24][25]. Now consider that this extra scalar mediates the interaction with DM as follows where H T = (0, h + v EW )/ √ 2 is the SM Higgs doublet and S = s + v S the new scalar, where v EW and v S are the vev of the SM Higgs and the new singlet scalar, respectively. The Higgs-DM interactions read, after EWSB, where the mixing angle between s and h is θ, h 1 is the 125 GeV Higgs and h 2 the heavier one The fermionic dark matter is assumed to get its mass entirely from the new scalar vev; thus it's mass, and Yukawa coupling is related by m χ = y χ v S . We highlight that throughout we will not impose the relic density constraints because for the couplings adopted here one may reproduce the right relic dark matter density either via the thermal freeze-out or non-thermal production of dark matter [29]. If a dominant non-thermal production is assumed the bounds stemming from direct and indirect detection might be circumvented [30]. Due to the mixing of scalars, the SM Higgs boson can decay to a DM pair but with a sin θ suppression. Otherwise, even if no mixing of the scalars is allowed, the SM Higgs decays to DM at one-loop mediated by the heavy scalar. These decays are possible due to the scalars self-interactions The heavy Higgs decay to dark matter, by its turn, can have a sizable decay branching ratio, suppressing its decay to SM Higgs pairs.
The mixing with the new scalar induces deviations in the SM Higgs trilinear coupling given by where b 3 and b 4 are free parameters, and θ is assumed to be small in view of the current bound | sin θ| 0.20 coming from the 1-loop correction to the W boson mass [31].
From Eq.(5) above, we see that a heavy new scalar shifts λ 111 away from its SM value toward large positive values. In Fig. (1) we display κ λ as function of the h 2 mass for some fixed mixing angles and other parameters. We see that scenarios with large shifts for m h2 > 500 GeV are possible. For small mixings, however, large κ λ will occur for h 2 masses in the region of TeVs. For h 2 masses smaller than ∼ 1 TeV, κ λ is more sensitive to variations in b 3 , b 4 and v S , but for the typical parameters we are interested in, the picture does not change much. A dark and bright signature, pp → γγ+ E T , might receive several contributions in xSM. When θ is non-vanishing, treelevel h 1 → χχ and 1-loop h 2 → γγ are possible. Besides the double Higgs production pp → h 1 h 1 → γγ + χχ, more contributions to the process γγ+ E T arise: all of them involving trilinear couplings, and (4) pp → h i h j h k , pp → h i h j and pp → χχ + h 1,2 involving quartic couplings. The h 1 h 1 production with an enhanced trilinear λ 111 coupling, when θ is small and h 2 is heavy, is the dominant contribution while the other contributions are expected to be negligible. We explain them in details below.
(1) pp → h 2 → h 1 h 1 . This is the resonant contribution where the heavy Higgs boson h 2 decays to a pair of SM Higgs bosons. This contribution might produce hard photons and large missing transverse momentum but it is suppressed by the scalar mixing angle as sin 2 θ in the gluon fusion production of h 2 . The decay into dark matter, by its turn, diminishes the branching to h 1 h 1 . The cross section is low for large h 2 masses.
(2) pp → h 1 h 2 . There are three contributions here, actually, with either h 1 or h 2 off-mass shell, and a box diagram. The pp → h * 1 → h 1 h 2 is expected to have a too low cross section. The other channel, pp → h * 2 → h 1 h 2 , is suppressed by sin 2 θ. Note that, in these cases, DM can be produced directly in the decay of h 2 . The branching ratio of h 2 → γγ is suppressed by the mixing angle just like the effective h 2 gg coupling. The box diagram is also suppressed by sin θ.
(3) pp → h 2 h 2 . This is analogous to pp → h 1 h 1 with a triangle and a box contribution where h 2 couples to the top quark. The triangle contribution is suppressed by sin θ whereas the box contribution by sin 2 θ. Moreover, the photons should be produced by the h 2 decay. There is another contribution where h * 1 → h 2 h 2 , but again, its cross section is expected to be very small.
(4) Triple Higgs production, for example, pp → h 2 /h * 2 → h 1 h 1 h 1 , also produces photons plus missing energy but with a tiny production cross-section and an additional suppression coming from the invisible branching ratio of the third Higgs. Double production involving quartic couplings occurs at the 2-loop level only. Finally, pp → χχ + h 1 , with the quartic interaction λ 2111 is also possible but again at the 2-loop level. All these contributions are expected to be small too.
If the mixing angle is not too small and h 2 is not too heavy, it is expected that the resonant contribution (1) substantially enhances the h 1 h 1 production cross section and helps to disentangle the background events efficiently due to the production of hard photons and missing energy. The other channels might also contribute to this scenario. However, if the mixing angle is too small or m h2 is of O(TeV), or both, even the resonant contribution will be negligible. In this case, the sole contribution is pp → h 1 h 1 → γγ + χχ and, as we discussed before, a significant deviation in the trilinear coupling might be possible. Otherwise, model parameters that lead to large Higgs masses and small shifts will be tough to probe at the HL-LHC. We will adopt a conservative and less modeldependent approach and suppose that contributions from h 2 are negligible.
Our results can thus apply to models where the trilinear coupling is shifted away from its SM value. There exist other models, like the Composite Higgs boson model [32][33][34], and the Two Higgs Doublet Models [35], which also predict non-standard trilinear couplings which might enhance double Higgs production. If the Higgs boson decays to dark matter, our results might apply.
In this paper, we are going to explore the double SMlike Higgs h 1 production and decay to an invisible channel plus two photons where the only possibility for an enhanced cross section is by changing the triple Higgs coupling λ 111 . As is well known, the triangle and the box contribution to gg → h 1 h 1 interfere destructively, but for |λ 111 | sufficiently large, the production cross section grows beyond the SM value. Models with new heavy states contributing sizeably to h 1 h 1 production should be easier to explore once they are more discernible from backgrounds. For this moment, we stick to the hardest scenario with no resonances or new contributions to the double Higgs production.
The motivation to explore a final state with large missing energy is that once the Higgs boson possesses an invisible decay mode with branching ratio BR inv , branching ratios of all other decay channels drop by a factor (1 − BR inv ). In particular, for double Higgs production, any final state containing just the SM particles will be suppressed by a factor (1 − BR inv ) 2 . The current LHC bound on the Higgs decay to invisible particles (19% [21]) means that any SM only final state studied so far, will have 1/3 less number of events. Consequently, the LHC reach for di-Higgs production into SM channels would suffer a substantial depletion. Exploring channels with one of the Higgs boson decaying invisibly is an interesting exercise to evaluate their potential in the search for double Higgs production under the hypothesis of the existence of an invisible channel. In the case where the HL-LHC has no sensitivity to the dark matter decays, it is still essential to determine the potential of the collider to explore the classic di-Higgs channels, like bbγγ, with reduced branching ratios in scenarios with large trilinear couplings.

III. MODEL CONSTRAINTS
A number of theoretical and phenomenological constraints apply to this model. First, the scalar potential need to be bounded from below, which leads to the following conditions Next, we impose that the potential be stable (at T = 0) by solving the equations We also require that all the 2 → 2 scattering amplitudes involving W, Z, h 1 and h 2 bosons respect perturbative unitarity at the high energy limit. Details of these calculations can be found in Ref. [27].
From the experimental side, we take the constraints for the Higgs mixing angle [36], the Higgs branching into invisible [21], the current bound in the trilinear coupling shifts κ λ [28], and W mass corrections [31,37] | sin θ| < 0.33, For | sin θ| < 0.33, h 2 masses up to a few TeV also respect constraints from oblique parameters S, T and U [6]. The W mass bounds imply | sin θ| < 0.2, effectively.
The branching fraction of the Higgs boson into fermionic DM (invisible) is given, in this model, by where Γ SM = 4.07 MeV. The branching ratio of h 1 to any other SM particles is multiplied by 1 − BR(h 1 → χχ). In our model, the DM-nucleus scattering is spinindependent and nearly isospin conserving. Thus, one can easily compute the dark matter-nucleon scattering cross section and compare with most stringent limits up-to-date that stem from the XENON1T experiment. In cm 2 we find the DM-nucleon spin-independent scattering cross section to be [6,38,39], where µ χN is the reduced mass of the proton-DM system. We impose the latest XENON1T limits [40] on the model. We point out that there are other competitive bounds in the literature stemming from other experiments, but we will adopt the XENON1T for being the most restrictive [41,42].
It is known that fermionic DM Higgs portal models typically lead to an over-abundance of DM in the universe unless the DM mass lies very close to the threshold of resonant annihilations via one of the Higgs bosons [7]. In our case, the SM Higgs and the new Higgs can play the role of mediators between the DM and the SM particles. We assume that the new Higgs is heavy. Thus, DM annihilations through h 2 are not efficient, unless we are near its resonance, but this will not happen here because the dark matter mass is sufficiently small. Because we want to explore the LHC prospects for Higgs into invisible states, we restrict the dark matter mass to m χ < 60 GeV (14) and keep y χ in the perturbative regime. As aforementioned we will assume that the correct relic dark matter density is somehow achieved invoking non-thermal effects that could bring down the abundance to its correct value as measured by PLANCK experiment [43][44][45].
With the constraints of Eqs. (6)(7)(8)(9)(10)(11)(12)(13)(14) at hand, we scan over the relevant portion of the parameters (θ, b 3 , b 4 , v S , m h2 , m χ ) of the xSM augmented with the fermionic DM interactions. In Fig. (2), we show the number of models as function of κ λ which respect all the constraints in six scenarios from θ ≤ 0.03 to θ ≤ 0.20. As θ gets large, the number of models with large trilinear shifts increases, however, the majority of the models lie in the 1 < κ λ < 2 region.

IV. THE SEARCH CHANNEL γγ+ ET
The di-Higgs channel involving E T , which produces the most significant number of events is, naturally, bb+ E T that was studied in Ref. [22]. In that work, a cut-and-count and a machine learning (ML) analysis were performed for both nonresonant and resonant scenarios within the xSM. The nonresonant study was performed for the SM trilinear coupling only, and ∼ 2.3σ and ∼ 4.3σ signal significances were estimated from the cut-and-count and the ML analysis, respec- tively. As we mentioned before, the study does not take systematic uncertainties into account and assume that errors are statistically dominated. The caveat of those results is that the signal-to-background ratio, S/B, is tiny, amounting to just 0.03 in the most powerful computation using BDT. The backgrounds to this channel include bbZ, bbW , Zh 1 and tt. Once even a minimal uncertainty is taken into account in the number of background events, the significance drops below 1σ. The situation in the resonant case is most promising being less sensitive to systematic effects.
The signal cross section of γγ+ E T suffers from the low branching ratio to photons but, the backgrounds to this channel are expected to be much lower compared to bb+ E T . The dominant ones are: (1) continuum Zγγ, (2) qq → Zh 1 , and (3) gg → Zh 1 at 1-loop level, where the Z boson decays to neutrinos. The reducible contributions W γγ and W h 1 can be neglected after vetoing a hard charged lepton and imposing the full selection cuts in all the subsequent analysis. In fact, we are going to show that despite starting with a much smaller signal cross section due to low h → γγ branching ratio, we can achieve a signal-to-background ratio that is an order of magnitude higher than that of bb+ E T .
We simulate signal and background events at the 14 TeV LHC using MadGraph5 aMC@NLO v2.3.6 [46] at leading order. Higher order QCD corrections are included in total cross sections via appropriate K-factors given in Table (I). The dark matter mass is fixed at 50 GeV in the simulations. The distributions depend weakly on the DM mass once the missing transverse momentum corresponds to the Higgs p T up to detector effects. We use Pythia8 [47] for hadronization and showering, and the detector simulation is performed by Delphes3 [48]. The number of SM signal and backgrounds events with QCD corrections, assuming 3 ab −1 , are shown in Table (I) with the following basic selection criteria  where the minimum transverse momentum within the fiducial volume of the electromagnetic calorimeters were required for the two photon hardest photons of the event. Here, p T (γ), η(γ), ∆R γγ are the transverse momentum, rapidity and the distance between the photons in the η × φ plane, respectively, while E T is the missing transverse momentum of the event.

V. CUT ANALYSIS AND COMPARISON WITH bb+ ET AND bbγγ CHANNELS
First of all, we want to show that once systematic uncertainties are taken into account, γγ+ E T presents better prospects than bb+ E T at the HL-LHC. We are going to compare the signal significance of both channels for the same level of systematic effect despite this is not likely to be realistic. Jet energy calibration, b-tagging efficiency and a higher number of backgrounds with QCD radiation are some features which are expected to impact more severely the bb+ E T channel [51] while the much more precise photon identification of the detectors and non-QCD backgrounds favors γγ+ E T .
In order to maximize the signal significance in the photons plus missing energy mode, we tuned the kinematic cuts using CutOptmize [52], a Python package aimed to maximize the statistical significance in particle physics searches at colliders in cut-and-count and ML analysis using a Gaussian Process algorithm as implemented in Hyperopt [53]. Especially, CutOptmize learns to tame the systematic uncertainties in the number of background events by raising S/B. We did not attempt to optimize the bottoms channel but we believe that optimizing cuts can also raise the S/B ratio making it less sensitive to a common systematic uncertainty, ε sys . The number of signal and background events for bb+ E T was taken from Ref. [22] and the significance metrics as well, which is given by .
The kinematic variables used for cuts and the BDT classification are the following: (1) the transverse momentum of the two leading photons: the distance between the photons in the η × φ plane: the number of jets in the event: N j with p T > 20 GeV and |η| < 2.5.
In Fig. (3), we show the signal and background distributions of the variables E T and the distance between the photons, ∆R γγ at the upper and lower panel, respectively. The √ŝ min (0) spectrum displays a behavior which is observed in other kinematic variables with energy dimensions as well. They are: (1) the Higgsstrahlung qq → Zh 1 and the continuum Zγγ backgrounds are softer than the signals for trilinear couplings around the SM value, (2) the harder spectra occur for λ 111 closer to λ SM , the values with the strongest destructive interference and, consequently, with smaller cross sections, (3) as |λ 111 | increases, the signal distributions get softer, (4) the gg → Zh 1 background resembles the SM signal distributions more closely. The photons distance in η × φ plane, ∆R γγ , presents a similar behavior: the Higgsstrahlung qq → Zh 1 and the continuum Zγγ backgrounds are more easily distinguishable from the signals with λ 111 ∼ λ SM , however, gg → Zh 1 is, again, similar to the SM.
A note of caution is necessary at this point. For photons produced centrally in the detector, |η γ | < 1.5, and with energies up to ∼ 200 GeV, the Delphes3 parametrization of the photon energy resolution leads to a full width at half of maximum (FWHM) of the photons pair mass of 5 GeV which is roughly twice as large as the CMS and ATLAS resolutions [28]. This difference can be attributed in part to the simplified approach to the photons simulation from Delphes3 which neglects e + e − conversions in the electromagnetic calorimeter. It is not clear, however, if the discrepancy has other sources. Our typical cut requirements are hard enough to select events with high energy though. For hard photons, the energy resolution is not an issue and we can trust the distributions generated by Delphes3. Moreover, we do not select events with too narrow γγ widths, limiting the distance to the Higgs mass to 5 GeV at least, large enough to accommodate the FWHM from Delphes3.
The production cross section of Higgs pairs depends quadratically on λ 111 where the minimum occurs approximately at ∼ 2.5λ SM , the same point where the signal kinematic distributions get maximally discernible from the backgrounds as we see in Fig. (3). Optimizing the cuts raises the cut efficiency and the background rejection for each λ 111 and makes the total cross section a more determinant factor across the analysis. FIG. 4: The 1, 2, and 5σ reaches of the γγ+ ET and bb+ ET channels for three systematics scenarios in a cut analysis. The cuts of γγ+ ET are optimized for each κ λ . In the lower panel, we also show the ATLAS [56] prospects for bbγγ in the 5 and 10% systematics cases. The color code denotes the density of xSM models that pass the criteria from Eqs. (6)(7)(8)(9)(10)(11)(12)(13)(14).
In order to get a straightforward comparison with the results of bb+ E T [22], we fixed the cut efficiency for the SM trilinear coupling and just rescaled the results for other λ 111 by µ hh = σ(pp → h 1 h 1 ) BSM /σ(pp → h 1 h 1 ) SM . This approximation was adopted in Ref. [22] to estimate the 5σ reach and the 95% CL exclusion limit in the µ hh versus BR(h 1 → invisible) plane. We borrow the number of signal (for the SM case) and background bb+ E T final state events after all cuts, which are S = 298 and B = 11231 respectively, from Ref. [22]. For each µ hh we calculate the corresponding κ λ .
We show, in Fig. (4), the 1, 2, and 5σ reaches after 3 ab −1 at the 14 TeV LHC. First, we reproduced the results of Ref. [22] using the information available, they are represented by the dashed lines. Note that in the absence of systematic errors, bb+ E T can exclude the SM trilinear coupling down to BR(h 1 → invisible) ∼ 6% and all κ λ down to BR inv ∼ 15% at 68% CL. Assuming the current LHC limit of BR inv = 0.19, trilinear couplings in the interval |κ λ − 2.5| 1 can be excluded at 95% CL and it is capable to discover di-Higgs production down to |κ λ − 2.5| 2, as we see in the upper panel of fig. (4). In the statistics dominated scenario, bottoms plus missing energy is more promising than γγ+ E T whose limits and discovery prospects are shown in solid blue lines. We emphasize that these limits take into account just the total number of expected events, in excess of backgrounds, for a given coupling and invisible branching ratio. More stringent limits might be possible in a shape-based analysis of a suitable kinematic distribution just like in those cases where it is possible to reconstruct the di-Higgs mass.
Once we include systematic uncertainties, however, the γγ+ E T channel becomes competitive. Contrary to bb+ E T , whose cuts are fixed, we tuned the cut thresholds of the events of photons plus missing energy aiming the maximization of the signal significance of Eq. (16) for each integer κ λ from −5 up to 12. The optimized cuts for the SM case are quoted in Eq. (17) with no systematic errors included In this case, and for other non-SM λ 111 as well, the cuts found in the optimization process confirm our intuition about the signal rich region in some kinematic variable distributions but, not always. The SM is a good example. The cut of ∆R γγ < 1.5 seems obvious if we take a look at the bottom panel of Fig. (3) but not a loose cut on E T . The cut in the missing transverse energy turned out to be redundant once a hard cut on √ŝ min (0) is imposed as we see in Eq. (17) above. For the SM trilinear coupling, the number of signal events surviving these criteria is 17, and 159 for the total background, corresponding to 1.4σ, and S/B = 0.11 with no systematics included. The signal-to-background ratio is four times larger than that found in the bb+ E T study but the signal significance drops by half. As we discussed, this much larger S/B ratio makes γγ+ E T more promising than bb+ E T when we include systematic effects. In all the subsequent analysis, we tuned the cuts in order to maximize Eq. (16) with ε sys = 0; however, it would be possible to tune the cuts for each systematics level to raise S/B even more.
In the middle panel of Fig. (4), we display the results for 5% systematic errors in both signal and background normal-izations. The limits of the photons channel move by less than one unit in κ λ compared to the 0% systematic error case, approximately, while the bottoms channel gets much less constraining compared to the statistics dominated scenario. Yet, bb+ E T is still more constraining by a small margin. Raising the systematic errors to 10% now flips the situation: the γγ+ E T becomes a better channel to probe di-Higgs production with dark matter decays as we see in the lower panel of Fig. (4) where we also display the reach of the bbγγ channel using the results of the ATLAS Collaboration [56] projecting the prospects of the 14 TeV LHC and 3 ab −1 in a cut-andcount study. To calculate the corresponding κ λ of a given invisible branching ratio, we multiplied the number of SM signal events quoted by ATLAS by (1 − BR inv ) 2 and then obtained the κ λ that would enhance the significance to 1, 2, and 5σ. As a final remark, note that bbγγ probes more significant regions in κ λ than the other two channels for all the relevant invisible branching ratios in the cut analysis.

VI. BDT ANALYSIS
Although it improves the reach of the LHC for double Higgs production with one invisibly decaying Higgs boson, the γγ+ E T channel is still not competitive with bbγγ in constraining the Higgs trilinear coupling in the presence of an invisible decay mode as we discussed in the last section. The prospects change when we perform a multivariate analysis.
We used BDTs, as implemented in XGBoost [57], to better classify our signal and background events. Just like the cutand-count analysis, we tuned the cuts but also the BDT hyperparameters in a joint optimization of the signal significance of Eq. (16) with no systematics which were included just in the final computation of the statistical significance. The optimization of a dedicated classifier was performed for all the trilinear couplings corresponding to an integer κ λ from −5 to 12 taking into account the changes in the fifteen kinematic distributions described in Section V.
We split our dataset, of around 100 thousand events for each one of the three background classes and the signal class (totaling 400k events), in two equal parts, one for testing and the other for training. We used CutOptimize to perform 200 iterations in the joint parameters space to maximize the significance. In each iteration, we randomly split the data 5 times in training and test sets of equal parts and calculated the BDT outputs of each class, estimating their distributions. These distributions were then used to place a final cut in order to better separate the signal and the backgrounds. The median of the five significance signals gives a final significance which constituted the optimization objective of the maximization algorithm. Once the best cuts and BDT hyperparameters were found, we perform a final 10-fold cross validation by randomly splitting the data set in train/test sets of equal size to estimate the final significance and its variance. In all the final results, the variance of the signal significance was small indicating the robustness of the best parameters found.
A feature importance analysis using the SHAP [58][59][60] package, shows that E T , M γγ , √ŝ min (0), M T B and N jets are the most important variables for the BDT classification for the majority of trilinear couplings. We took the number of signal and backgrounds events for bb+ E T from Ref. [22] after their BDT analysis using TMVA, namely, S = 593 and B = 19466 but again fixing the cut efficiency for other κ λ values and just scaling the signal significance by the total cross sections of the signal. It must be emphasized that this approach is only an approximation since the cut efficiency of the signal varies with κ λ .
We also took the prospects for the bbγγ channel from a recent CMS study [28] of the di-Higgs production at the HL- The 68, 95% CL limits, and 5σ reaches of the γγ+ ET and bb+ ET channels for 5% systematics scenarios in the BDT analysis. The cuts of the dark and bright channel are optimized for each λ. In all panels, we also show the CMS [28] prospects in the bbγγ.
LHC using BDTs to purify the samples and using a parametric maximum likelihood fit to photons, m γγ , and bottoms, m bb , masses distributions. In order to estimate the couplings that can be probed with reduced branching ratios into bb and γγ due to the presence of the dark matter decay, we just rescaled the signals significance of the CMS work by (1 − BR inv ) 2 and looked for the κ λ with that new, smaller, significances.
In the upper panel of Fig. (5), we display the 68% CL limits of the three search channels in the κ λ × BR inv plane for 5% systematics in both signal and background normalizations. First of all, we note that the optimization process of γγ+ E T is able to find parameters for which the results in terms of signal significance vary smoothly as a function of κ λ and BR inv . This behavior indicates that the optimization algorithm finds the path of the points of the maximum of the objective function in the multidimensional parameters space. The γγ+ E T and bb+ E T channels perform nearly the same but, again, the bottoms channel is not optimized. The prospects using bbγγ events are the better at this confidence level for all the relevant invisible Higgs branching ratios but, as expected, these limits soften as BR inv gets larger. The color code of the heatmaps superimposed on the plots of Fig. (5) shows the density of models respecting all constraints of Eqs. (6)(7)(8)(9)(10)(11)(12)(13)(14) in the κ λ versus BR inv plane. The blue region concentrates the majority of models. Part of that region can be probed in bbγγ and γγ+ E T at 68% CL, as we see in the upper panel. The blank regions do not contain viable points of xSM, but we show them once they might constrain other interesting models with negative shifts in the trilinear coupling.
The middle panel of Fig. (5) presents the 95% CL limits. The dark matter channels now probe a region moved approximately one unit in |κ λ | upwards compared to the 68% CL case. However, the bbγγ limits, as obtained by the CMS Collaboration, soften more intensely. Interestingly, a complementarity now arises where bbγγ probes small invisibly Higgs decays up to ∼ 8%, while the other two channels constrain higher branching ratios. Only trilinear couplings either larger than 5.5λ SM or smaller than 0.5λ SM can be observed at this confidence level.
In the last panel of Fig. (5), we show the discovery prospects assuming, again, a 5% systematics. Now, γγ+ E T performs better than bb+ E T for all invisible branching ratios. This channel is also better than bbγγ for BR inv 8%. We estimate that it is possible to discover di-Higgs production and decay to bright and dark states for λ 111 7.7λ SM in xSM, and λ 111 −1.4λ SM in models with negative trilinear shifts, for BR inv < 19%.
Finally, we raise the systematics to 10% and estimate the prospects of the HL-LHC at the 95% CL using BDTs. In Fig. (6), we see that κ λ down to 6 can be probed for BR inv = 19% with γγ+ E T events and a complementarity with bbγγ occurs around BR inv = 10%. Compared to these topologies, raising the systematics level has a more deleterious effect on the reach of bb+ E T channel. The 95% CL limits of the γγ+ ET and bb+ ET channels for 10% systematics scenarios in the BDT analysis. The cuts of the dark and bright channel are optimized for each κ λ . We also show the CMS [28] prospects in the bbγγ channel.

VII. CONCLUSIONS
The production of two Higgs bosons is of prime importance to understand the scalar sector of the Standard Model. Determining the self-interactions of the Higgs boson might reveal a more profound structure of the scalar potential with clues to important open questions as the stability of the vacuum and an electroweak phase transition that has driven the baryon matter asymmetry, for example.
Along with the nature of the Higgs sector, one of the most obvious targets of colliders is dark matter. If just like the other known particles, the dark matter has its mass generated through the Higgs mechanism, it is plausible to observe a Higgs-DM interaction. Especially, if the Higgs boson decays to DM, new possibilities to search for double Higgs production may open up. On the other hand, the number of events from decay channels involving only the SM particles are expected to be reduced by a factor of (1 − BR inv ) 2 , which is around 2/3 given the current LHC bound on the invisible decay width of the Higgs. All this makes the study of double Higgs production and decay to DM an interesting and timely task.
In this work, we investigated the di-Higgs production in three complementary channels: γγ+ E T , bb+ E T and bbγγ. The prospects of the LHC in the dark and bright channel γγ+ E T were estimated for the first time in both a cut-andcount and multivariate analysis using a dedicated algorithm to optimize their signal significances. The estimates for bb+ E T and bbγγ were taken from the literature adapting to the inclusion of the Higgs decays to DM when necessary.
Our approach aimed at the singlet scalar extension of the SM, the xSM, which is the simplest extension of the SM that leads to first order EWPT. Under the assumptions that we made, namely, a massive new Higgs boson with small mixing with the SM Higgs boson, a more or less modelindependent estimate follows and an effective field theory ap-proach is reliable. In order to keep new ingredients to the SM to a minimum, we augmented xSM with a fermionic DM coupling to the new scalar. Concerning the DM sector, we demanded points are respecting both the latest XENON1T bounds. Moreover, for DM masses not too close to Higgs mediator thresholds the relics abundance bound can also be evaded.
For the xSM, we found that portions of the parameters space that will give us significant deviation of the trilinear Higgs coupling relative to the SM value can be probed by, at least, one of the decay channels studied here. The γγ+ E T channel, in particular, becomes a better option than bb+ E T once systematic uncertainties are taken into account. Moreover, the dark and bright channel presents better prospects for exclusion and discovery than bbγγ for BR inv as low as ∼ 10% in the multivariate analysis using BDTs. Overall, all three channels complement each other in the region 0 < BR inv < 19% and −5 ≤ κ λ ≤ 12, which is still allowed by the data. Our results also show that those parameters of xSM with a trilinear coupling close to the SM and an invisible Higgs branching ratio smaller than ∼ 15% will be tough to be probed at the LHC relying only on shifts of the SM Higgs selfinteractions. Scenarios with not so heavy new Higgs bosons are potentially more easily accessed in colliders and a complementary study across many search channels like this one is due.