Spotting hidden sectors with Higgs binoculars

TL;DR

This paper proposes new search strategies for detecting hidden sectors at the LHC through events with two Higgs bosons and missing energy, using machine learning to optimize discovery potential.

Contribution

It introduces simplified models for hidden sectors with two event topologies and develops optimized machine learning-based analysis methods for these signatures.

Findings

LHC can test significant parameter space for hidden sectors with these signatures.

Di-Higgs production with missing energy is a competitive search channel.

Machine learning enhances signal-background discrimination in these searches.

Abstract

We explore signals of new physics with two Higgs bosons and large missing transverse energy at the LHC. Such a signature is characteristic of models for dark matter or other secluded particles that couple to the standard model through an extended scalar sector. Our goal is to provide search strategies and an interpretation framework for this new signature that are applicable to a large class of models. To this end, we define simplified models of hidden sectors leading to two different event topologies: symmetric decay, i.e., pair-produced mediators decaying each into a Higgs plus invisible final state; and di-Higgs resonance, i.e., resonant Higgs-pair production recoiling against a pair of invisible particles. For both scenarios, we optimize the discovery potential by performing a multi-variate analysis of final states with four bottom quarks and missing energy, employing state-of-the-art machine learning algorithms for signal-background discrimination. We determine the parameter space that the LHC can test in both scenarios, thus facilitating an interpretation of our results in terms of complete models. Di-Higgs production with missing energy is competitive with other missing energy searches and thus provides a new opportunity to find hidden particles at the LHC.

Tables

Table 1. Table 1: Benchmarks S ​ _ ​ m B ​ _ ​ m A ​ _ ​ m χ 𝑆 _ subscript 𝑚 𝐵 _ subscript 𝑚 𝐴 _ subscript 𝑚 𝜒 S\textunderscore m_{B}\textunderscore m_{A}\textunderscore m_{\chi} (symmetric topology, left) and R ​ _ ​ m B ​ _ ​ m A ​ _ ​ m χ 𝑅 _ subscript 𝑚 𝐵 _ subscript 𝑚 𝐴 _ subscript 𝑚 𝜒 R\textunderscore m_{B}\textunderscore m_{A}\textunderscore m_{\chi} (resonant topology, right) for a di-Higgs plus E̸ T subscript italic-E̸ 𝑇 \not{E}_{T} signature. The second column shows the signal rate σ S , R ​ ( p ​ p → b ​ b ¯ ​ b ​ b ¯ ​ χ ​ χ ) subscript 𝜎 𝑆 𝑅 → 𝑝 𝑝 𝑏 ¯ 𝑏 𝑏 ¯ 𝑏 𝜒 𝜒 \sigma_{S,R}(pp\to b\bar{b}\,b\bar{b}\,\chi\chi) at a centre-of-mass energy of 14 TeV, assuming a UV completion by a vector-like quark with m Q = 1 ​ TeV subscript 𝑚 𝑄 1 TeV m_{Q}=1\,\,\mathrm{TeV} and y Q = 1 subscript 𝑦 𝑄 1 y_{Q}=1 , see Appendix  A .

symmetric benchmark #	${𝝈}_{𝑺}$ [fb]
$S\_1000\_475\_\{340,250,150,50\}$	2.94
$S\_1000\_400\_\{265,175,75\}$	2.99
$S\_1000\_325\_\{190,100,0.1\}$	3.00
$S\_1000\_250\_\{115,25\}$	3.01
$S\_750\_350\_\{215,125,25\}$	15.25
$S\_750\_275\_\{140,50\}$	15.32
$S\_750\_200\_\left\{65\right\}$	15.33
$S\_500\_225\_\{90,50,0.1\}$	106.10
$S\_500\_200\_\{65,25\}$	106.12
$S\_500\_175\_\{40,0.1\}$	106.14
$S\_500\_150\_\left\{15\right\}$	106.15

Equations

$$pp\to B\to AA\to (h\chi )(h\chi ).$$

$$pp\to B\to AA\to (hh)(\chi \chi ).$$

$$\mathcal{L}\supset \frac{C_{Bgg}}{\Lambda }B{G}_{\mu \nu }^a{G}^{\mu \nu a}+\frac{m_{BAA}}{2}BAA.$$

$${\mathcal{L}}_S\supset {\lambda }_{A\chi HH}A\chi H^{†}H$$

$$pp\to B\to AA\to (h\chi )(h\chi )$$

$${\mathcal{L}}_R\supset m_{AHH}A{H}^{†}H+\frac{m_{A\chi \chi }}{2}A\chi \chi .$$

$$\Gamma (A\to WW)=2\Gamma (A\to ZZ)=2\Gamma (A\to hh).$$

$$pp\to B\to AA\to (hh)(\chi \chi ),$$

$$pp\to B\to AA\to (VV)(\chi \chi ),pp\to B\to AA\to (VV)(V^{\prime }{V}^{\prime }).$$

$$\mathcal{L}=\frac{C_{Bgg}}{\Lambda }B{G}_{\mu \nu }^a{G}^{\mu \nu a}+\frac{m_{BAA}}{2}BAA+m_{AHH}A{H}^{†}H+\frac{m_{A\chi \chi }}{2}A\chi \chi +{\lambda }_{A\chi HH}A\chi H^{†}H.$$

$$C_{Bgg}=2.1\cdot 10^{-3},\Lambda =1\text{TeV},m_{BAA}=v=246\mathrm{GeV},m_{A\chi \chi }\approx 2m_{AHH}.$$

$${\sigma }_S(pp\to b\bar{b}b\bar{b}\chi \chi )=\sigma (pp\to B)\times \mathcal{B}(B\to AA)\times {\mathcal{B}}^2(A\to h\chi )\times {\mathcal{B}}^2(h\to b\bar{b}),$$

$${\sigma }_R(pp\to b\bar{b}b\bar{b}\chi \chi )=\sigma (pp\to B)\times \mathcal{B}(B\to AA)\times 2\mathcal{B}(A\to hh)\mathcal{B}(A\to \chi \chi )\times {\mathcal{B}}^2(h\to b\bar{b}).$$

$${\sigma }_S(pp\to b\bar{b}b\bar{b}\chi \chi )=\sigma (pp\to B)\times {\mathcal{B}}^2(h\to b\bar{b})\approx 8{\sigma }_R(pp\to b\bar{b}b\bar{b}\chi \chi ).$$

$${ϵ}_b=0.8\frac{30\tanh (3\cdot 10^{-3}{p}_T)}{1+8.6\cdot 10^{-2}{p}_T},{ϵ}_l=0.002(1+3.65\cdot 10^{-3}{p}_T),{ϵ}_c=0.2\frac{\tanh (0.02p_T)}{1+3.4\cdot 10^{-3}{p}_T}.$$

$${\Delta }_{BA}=m_B-2m_A,{\Delta }_{Ah\chi }=m_A-(m_h+m_{\chi }),$$

$${\Delta }_{BA}=m_B-2m_A,{\Delta }_{Ah\chi }=\text{min}(m_A-2m_h,m_A-2m_{\chi }).$$

$$S\text{textunderscore}{m}_B\text{textunderscore}{m}_A\text{textunderscore}{m}_{\chi },R\text{textunderscore}{m}_B\text{textunderscore}{m}_A\text{textunderscore}{m}_{\chi },$$

$$S\text{textunderscore}750\text{textunderscore}350\text{textunderscore}25,R\text{textunderscore}750\text{textunderscore}350\text{textunderscore}25,S\text{textunderscore}1000\text{textunderscore}250\text{textunderscore}25,R\text{textunderscore}1000\text{textunderscore}275\text{textunderscore}25.$$

$$\Sigma =\frac{S}{\sqrt{S+B+(\beta B{)}^2}}.$$

$$m_{\chi }>m_h:$$

$$2m_{\chi }\approx m_h:$$

$$\mathcal{L}\supset -y_{Q}B\bar{Q}Q.$$

$$\frac{C_{Bgg}}{\Lambda }=\frac{g_s^2{y}_Q}{48{\pi }^2{m}_Q},\frac{C_{B\gamma \gamma }}{\Lambda }=\frac{e^2{y}_Q}{72{\pi }^2{m}_Q},$$

$$\mathcal{L}\supset \frac{C_{B\gamma \gamma }}{\Lambda }B{F}_{\mu \nu }{F}^{\mu \nu }.$$

$$C_{Bgg}=\frac{g_s^2}{48{\pi }^2}\simeq 2.1\cdot 10^{-3},$$

$$\Gamma (B\to gg)$$

$$\Gamma (B\to \gamma \gamma )$$

$$\Gamma (B\to AA)=\frac{{\left|m_{BAA}\right|}^2}{32\pi m_B}\sqrt{1-\frac{4m_A^2}{m_B^2}}.$$

$$\sigma (pp\to B)={\left(\frac{m_t{C}_{Bgg}}{m_Q{C}_{\hat{S}gg}}\right)}^2\sigma (pp\to \hat{S})\simeq {\left(\frac{y_Q}{m_Q}\frac{\sqrt{2}{m}_t}{y_t}\right)}^2\sigma (pp\to \hat{S})={\left(\frac{y_{Q}v}{m_Q}\right)}^2\sigma (pp\to \hat{S}),$$

Full text

I Introduction

Postulating a hidden sector that interacts primarily with the Higgs boson is tempting for good reasons. Higgs couplings to new scalar standard-model (SM) singlets are renormalizable and secluded from visible matter Patt:2006fw ; OConnell:2006rsp . An extended scalar sector can thus serve as a portal to a hidden sector LopezHonorez:2006gr ; MarchRussell:2008yu ; Englert:2011yb ; Athron:2017kgt . At the LHC, the Higgs interaction with a hidden sector is best probed in signatures with one or two Higgs bosons. Searches for invisibly decaying Higgs bosons or for mono-Higgs production in association with missing transverse energy $\not{E}_{T}$ are well-established parts of the LHC program. Invisible Higgs decays probe hidden sectors with particles significantly lighter than the Higgs boson. Mono-Higgs signals are often predicted in models that can also be probed in other channels, such as mono-jet production, mono-$Z$ production, or signatures with missing energy and several leptons and/or jets. For a review of missing energy searches at the LHC, we refer the reader to Refs. Kahlhoefer:2017dnp ; Abe:2018bpo and references therein.

A signal of two Higgs bosons and missing energy is naturally predicted in the context of supersymmetry (SUSY), for instance from Goldstino production in models with gauge mediated SUSY breaking Giudice:1998bp ; Matchev:1999ft , or from chain decays of superpartners into Higgs bosons and neutralinos in the minimal supersymmetric standard model (MSSM) and its extensions Kang:2015nga ; Bernreuther:2018nat ; Titterton:2018pba . More generally, di-Higgs plus missing energy is a signature of models with extra scalars Arganda:2017wjh , such as a pseudo-scalar portal to a dark sector No:2015xqa , axion-like particles Brivio:2017ije , massive right-handed neutrinos Kang:2015nga ; Kang:2015uoc or in the framework of Little Higgs scenarios Kang:2015nga ; Chen:2018dyq . Experimental searches for di-Higgs plus $\not{E}_{T}$ production at the LHC have been performed for a signal of four bottom quarks and missing energy in the context of SUSY Sirunyan:2017obz ; Aaboud:2018htj . This search targets a signature of Higgsino pair production, followed by a decay chain with Higgs bosons and Goldstinos in the final state Matchev:1999ft . Since the analysis is optimized for very light Goldstinos produced via this specific decay chain, its reinterpretation in other scenarios is limited. A systematic exploration of the di-Higgs plus $\not{E}_{T}$ channel at the LHC is still lacking.

Our goal is to provide a minimal, simple framework to exploit the full potential of the LHC to search for new hidden sectors with a di-Higgs plus $\not{E}_{T}$ signature. As a matter of fact, the search strategy for this signature strongly depends on the masses and decays of the relevant particles. Based on two main decay topologies, we define simplified models for $pp\to hh\chi\chi$ production, where $h$ is the SM Higgs boson and $\chi$ is invisible and stable at detector scales. Each model involves two scalar mediators $B$ and $A$, where $B$ couples to gluons and is heavier than two $A$ scalars. The first model, referred to as symmetric topology, is inspired by electroweakino production in the MSSM. A pair of on-shell scalars $A$ is produced from the decay of $B$. Each $A$ subsequently decays into a Higgs boson and an invisible scalar $\chi$. The di-Higgs signature is thus generated by

[TABLE]

In the second model, referred to as resonant topology, each of the pair-produced scalars $A$ decays into either two Higgs bosons or invisibly. The corresponding production chain is

[TABLE]

Such a topology is typical in scalar portal models. Since the definition of the two simplified models is based solely on the kinematic properties of the final state, LHC searches for these simple topologies can easily be recasted in terms of concrete models.

Our analysis focuses on the Higgs decay into bottom quarks, $h\to b\bar{b}$, which maximizes the event rates. The signal thus consists of four $b$-jets and a large amount of missing transverse energy. To reconstruct the two Higgs bosons from the four $b$-jets, we will make ample use of the mature analysis techniques for di-Higgs searches without associated missing energy. Due to its sensitivity to the Higgs self-interaction, Higgs pair production in the SM (see Refs. ATLAS:2018otd ; CMS:2018rig for the latest experimental prospects) and beyond (for a review see e.g. Ref. Zurita:2017sfg and references therein) is a key target of the LHC program and proposed future colliders Contino:2016spe . The prospects to observe a signal of Higgs pairs has evolved from “seemingly impossible” Baur:2003gp to a detailed investigation of the final states $b\bar{b}\,\tau^{+}\tau^{-}$ Dolan:2012rv , $b\bar{b}\,WW^{*}$Papaefstathiou:2012qe and $b\bar{b}\,b\bar{b}$ deLima:2014dta . This tremendous progress was triggered by exploiting novel techniques such as jet substructure and shower deconstruction Butterworth:2008iy ; Soper:2011cr ; Soper:2014rya . Today these techniques are applied by the ATLAS and CMS collaborations in experimental analyses of Higgs pair production Sirunyan:2017isc ; Aaboud:2018knk . In our search for $hh\chi\chi$ production with four $b$-jets and missing energy, we will combine jet substructure techniques with a state-of-the-art multivariate analysis to optimize the sensitivity to our signal. For Higgs pair production in the SM, the channel with four $b$-jets does not have the best performance, due to an immense multi-jet background. In contrast, due to the presence of large missing energy in our signal, the largest background arises from electroweak gauge bosons plus jets, which lies a few orders of magnitude below the QCD multi-jet processes that appear in SM di-Higgs searches. We therefore focus on the four-bottom final state, leaving other decay channels for future exploration.

Our article is organized as follows. In Section II, we introduce simplified models for $hh\chi\chi$ production. In Section III, we discuss the main features of the di-Higgs plus $\not{E}_{T}$ signature in our simplified models and the challenges we face in reconstructing the four-bottom final state, as well as triggering and backgrounds. We attempt a cut-based analysis and investigate its discovery prospects at the High Luminosity LHC (HL-LHC). In Section IV, we explain the details of our multi-variate analysis and demonstrate a large gain in sensitivity compared to the cut-based analysis. We stress that in the context of our simplified model, new physics could be discovered first in the di-Higgs plus $\not{E}_{T}$ channel, while being consistent with all existing (and future HL-)LHC searches. This highlights the importance of carrying out the proposed analysis. In Section V, we explore the validity of the dark matter interpretation of our model. We defer our conclusions to Section VI.

II Simplified models of a hidden scalar sector

In this section, we provide details on the two simplified models that give rise to the di-Higgs plus $\not{E}_{T}$ signature, but with different final-state topologies. In the symmetric topology, the two Higgs bosons stem from a chain decay within the hidden sector, while in the resonant topology they form a di-Higgs resonance. Since the two setups have some structural similarities, we first discuss their common features. We then move on to describe the specific ingredients of the models that lead to the different topologies.

Both models rely on an extended scalar sector with three new real scalar particles $A$, $B$, and $\chi$ that are singlets under the SM gauge group, with a mass hierarchy $m_{B}\gg m_{A}\gg m_{\chi}$. The models also feature a discrete $\mathbb{Z}_{2}$ symmetry under which particles belonging to the hidden sector are odd, while new particles in the visible sector as well as the SM particles are even. For the sake of simplicity, we assume that none of the new scalars develops a vacuum expectation value (VEV).

The heaviest of the three new scalars, $B$, is produced via gluon fusion at the LHC and predominantly decays to $AA$ pairs. The relevant interaction terms are

[TABLE]

Here we introduce an effective dimension-five interaction of the scalar $B$ with gluons such that it is resonantly produced via $gg\to B$, in analogy to the dominant Higgs production channel in the SM. We discuss a renormalizable UV completion for this interaction in Appendix A.

The triple scalar coupling $m_{BAA}$ induces the decay $B\to AA$ with a branching ratio near 100%, unless $m_{BAA}$ is very small compared with the Higgs VEV $v$. For values of $C_{Bgg}$ originating from perturbative physics around the TeV scale, the decay into dijets via $B\to gg$ then occurs only at the percent level. We also suppress the decays $B\to\chi\chi$, $B\to A\chi$, and $B\to hh$ by assuming the relevant couplings to be small. Note that a $Bhh$ coupling would induce $B$-$h$ mixing, which is severely bounded by measurements of the Higgs coupling strength Biekotter:2018rhp . As $B$ is produced in gluon fusion, it necessarily belongs to the visible sector, i. e., it is even under the $\mathbb{Z}_{2}$ symmetry. The scalar $\chi$ belongs to the hidden sector and is thus taken to be $\mathbb{Z}_{2}$-odd. As it is the lightest hidden particle, it is stable and appears as missing energy in the LHC detectors.

Depending on the $\mathbb{Z}_{2}$ parity of $A$, two different event topologies for di-Higgs plus $\not{E}_{T}$ can be distinguished,

• •

Symmetric topology. If $A$ is part of the hidden sector, i.e., $\mathbb{Z}_{2}$-odd, it decays via $A\to h\chi$. Di-Higgs plus $\not{E}_{T}$ arises from a symmetric event topology with chain decay in the hidden sector.

• •

Resonant topology. If $A$ instead belongs to the visible sector, i.e., if it is $\mathbb{Z}_{2}$-even, it can decay via $A\to hh$ and $A\to\chi\chi$. The di-Higgs plus $\not{E}_{T}$ signature then arises from an asymmetric event topology and features a di-Higgs resonance.

Figure 1 shows the event topologies for these two cases, which we now discuss in more detail.

II.1 Symmetric topology

In this model, both $A$ and $\chi$ are odd under the discrete $\mathbb{Z}_{2}$ symmetry, i. e. they belong to the hidden sector. The interaction 111Here $H$ denotes the SM Higgs doublet and gauge-invariant field contractions are assumed.

[TABLE]

induces the decay $A\to h\chi$ with a branching ratio of $100\%$. Due to their different $\mathbb{Z}_{2}$ parity, neither $A$ nor $\chi$ can mix with the Higgs boson and thus remain pure singlets under the SM gauge group. In particular, they do not couple to electroweak gauge bosons. The coupling $\lambda_{A\chi HH}$ induces mixing of $A$ and $\chi$ after electroweak symmetry breaking. However, also the mass term $m_{A\chi}^{2}A\chi$ contributes to this mixing. We conclude that the $A$-$\chi$ mixing and the decay $A\to h\chi$ are governed by two independent parameters, so that the mixing can be set to zero without affecting the $A\to h\chi$ decay in our signature. Since $B$ predominantly decays into $AA$, and $A$ decays exclusively into $h\chi$, the process

[TABLE]

is the main discovery channel for the symmetric topology at the LHC.

This model encompasses the well studied case of electroweakino production in the MSSM, with $A$ and $\chi$ corresponding to fermions (for instance, Higgsino production as in Ref. Matchev:1999ft ). A renormalizable coupling connecting $A$, $\chi$ and the SM Higgs implies that at least one of the fermions is non-trivially charged under the SU(2) electroweak group. Hence new states with electric charges appear, which often provide the leading collider signatures for these scenarios: multi-leptons plus $\not{E}_{T}$ (see e.g. Ref. Bramante:2015una ), mono-jet plus soft leptons (see e.g. Refs. Giudice:2010wb ; Schwaller:2013baa ; Low:2014cba ; Barducci:2015ffa ), disappearing tracks (for recent work see e.g. Refs. Mahbubani:2017gjh ; Fukuda:2017jmk ; Filimonova:2018qdc ; Saito:2019rtg ).

II.2 Resonant topology

In this model, both $A$ and $B$ are even under the $\mathbb{Z}_{2}$ parity, while $\chi$ is odd. The coupling $\lambda_{A\chi HH}$ is therefore forbidden and the decay $A\to h\chi$ is absent. Instead we introduce the couplings

[TABLE]

Both of these couplings were forbidden with the symmetry assignment leading to the symmetric topology. The coupling term $m_{AHH}$ induces a mixing of $A$ with $h$ after electroweak symmetry breaking. Unlike $A$-$\chi$ mixing in the symmetric model, $A$-$h$ mixing is unavoidable here: it is induced by the same parameter as the decay $A\to hh$ that is part of the di-Higgs plus $\not{E}_{T}$ signature. As a consequence, $A$ inherits the couplings of $h$, inducing $A\to WW$ and $A\to ZZ$ as relevant decay modes. In the limit $m_{A}\gg m_{h},m_{W},m_{Z}$, the decay rates fulfill the simple relation

[TABLE]

The decay $A\to\chi\chi$ is instead induced by the coupling $m_{A\chi\chi}$, so that its decay rate can be treated as an independent parameter. To maximize the significance of the signature

[TABLE]

we assume a branching ratio of $\mathcal{B}(A\to\chi\chi)=0.5$. In addition to the di-Higgs plus $\not{E}_{T}$ signature, this model also gives rise to signatures with di-boson resonances and signatures with four electroweak bosons $V,V^{\prime}=W,Z,h$ forming two di-boson resonances,

[TABLE]

These signatures and di-Higgs plus $\not{E}_{T}$ typically occur at similar rates. They complement each other in the search for scalar hidden sectors of this kind.

In summary, the two simplified models can be conveniently described by the interaction Lagrangian

[TABLE]

In the symmetric topology, the couplings $m_{AHH}$ and $m_{A\chi\chi}$ are equal to zero, due to the $\mathbb{Z}_{2}$ parity under which both $A$ and $\chi$ are odd. In the resonant topology, on the other hand, $\lambda_{A\chi HH}=0$, as only $\chi$ belongs to the hidden sector. Based on the field content and symmetries of both models, additional terms could be added to the respective Lagrangians. Along this work, we will only consider those that are relevant for the collider phenomenology of the di-Higgs plus $\not{E}_{T}$ signature  222 In general, the UV completion for $C_{Bgg}$ can induce an additional dimension-5 operator $C_{B\gamma\gamma}$, which gives rise to a di-photon resonance in the final state. Whether or not $C_{B\gamma\gamma}$ is correlated with $C_{Bgg}$ depends on the specific UV completion. Colored and electrically neutral particles contribute to $C_{Bgg}$, but not to $C_{B\gamma\gamma}$. Additional vector-like leptons affect $C_{B\gamma\gamma}$, but not $C_{Bgg}$. In the UV completion described in Appendix A, $C_{B\gamma\gamma}$ gives rise to a branching ratio $\mathcal{B}(B\to\gamma\gamma)$ of a few permille for all our benchmark points. We thus neglect $C_{B\gamma\gamma}$ in our analysis, as well as operators such as $C_{Agg}$ and $C_{A\gamma\gamma}$ (which would appear in the resonant model, but not in the symmetric one). . Using Eq. (10), we have implemented both Lagrangians into FeynRules Alloul:2013bka and used the Universal FeynRules Output (UFO) Degrande:2011ua for event generation.

Throughout our analysis, we fix the parameters in the hidden sector as follows,

[TABLE]

The value of $C_{Bgg}$ is motivated by a UV completion with a vector-like quark with mass $m_{Q}=\Lambda=1\,\text{TeV}$ and scalar coupling $y_{Q}=1$, as discussed in Appendix A. In the resonant model, the relation between $m_{AHH}$ and $m_{A\chi\chi}$ ensures a $50\%$ decay of $A\to\chi\chi$ in the limit $m_{A}\gg m_{h},\,m_{\chi}$. The magnitude of these couplings does not affect the signal rate. In the symmetric model, the branching ratio of $A\to h\chi$ is 100%, regardless of the size of the coupling $\lambda_{A\chi HH}$.

III Higgs-pair production with missing energy at the LHC

In order to develop a search strategy for di-Higgs plus $\not{E}_{T}$ at the LHC, we first analyze the kinematic features of this signature and their parameter dependence in each simplified model. We assume that $B$ is resonantly produced. The signal rate is then well approximated by

[TABLE]

The couplings of $B$ do not affect the decay kinematics. Unless the $AA$ pair is produced near threshold, the heavy scalar $B$ decays almost fully via $B\to AA$, with a branching ratio of $\mathcal{B}(B\to AA)\approx 1$. Away from the threshold, the production rate for $pp\to B\to AA$ thus depends only on the mass $m_{B}$ and the coupling $C_{Bgg}$. For fixed $m_{B}$, $m_{A}$, and $C_{Bgg}$, the number of produced $AA$ pairs is the same in both simplified models. Concerning the decays of $A$, in the symmetric topology $\mathcal{B}(A\to h\chi)=1$, while in the resonant topology the maximal decay rate into $hh\,\chi\chi$ is obtained for $\mathcal{B}(A\to\chi\chi)=0.5$. Taking into account the possible decays of $A$ into pairs of gauge bosons (see Section II), this corresponds to $\mathcal{B}(A\to hh)\approx 0.125$. The signal rate in the symmetric model is thus about 8 times higher than in the resonant model,

[TABLE]

The production cross section $\sigma(pp\to B)$ can be obtained by rescaling the SM Higgs production cross section, as described in Appendix A.

For our numerical analysis, we generate signal events at the parton level with MG5_aMC@NLO 2.6.1 Alwall:2014hca , using the NNPDF30_lo_as_0118_nf_4 Ball:2014uwa parton distribution functions implemented in the LHAPDF 6.1.6 Buckley:2014ana interface. We employ Pythia 8.2 Sjostrand:2014zea for parton showering and hadronization and Delphes 3.3.3 deFavereau:2013fsa for a basic detector simulation, using the default implementation of the ATLAS detector. Crucial inputs for our analysis are the $b$-tagging efficiency, $\epsilon_{b}$, and the light and charm jet mistag rates, $\epsilon_{l}$ and $\epsilon_{c}$, which are given by the following $p_{T}$-dependent functions

[TABLE]

Hence, for a jet transverse momentum of $p_{T}(j)=50\,(250)\,\,\mathrm{GeV}$, we find $\epsilon_{b},\epsilon_{l},\epsilon_{c}=67,0.24,0.26\,(73,0.3,1.0)\%$.

III.1 Kinematics and benchmarks

The kinematics of our signature is driven by the available phase space in the $B$ and $A$ decays. We parametrize this in either model in terms of the mass differences of the involved particles,

[TABLE]

We fix the scalar mass $m_{B}=\{1000,750,500\}\,\,\mathrm{GeV}$ and scan $\Delta_{BA}$ in steps of $75\,\,\mathrm{GeV}$ over the kinematically accessible region. In the symmetric model, $\Delta_{Ah\chi}$ is scanned in steps of $100\,\,\mathrm{GeV}$. We require $\Delta_{Ah\chi}>10\,\,\mathrm{GeV}$ to prevent too strong a phase-space suppression in $A$ decays. In the resonant model, we fix $m_{\chi}=25\,\,\mathrm{GeV}$ to satisfy the different kinematic boundaries.333As long as long as the $A\to\chi\chi$ final state is kinematically accessible, the mass $m_{\chi}$ has no impact on the di-Higgs plus $\not{E}_{T}$ phenomenology in the resonant model. The so-obtained benchmark scenarios for the symmetric ($S$) and resonant ($R$) model are labelled as

[TABLE]

where the masses of the scalars are given in units of GeV. They are shown for both models in Table 1, where we also give the corresponding signal rates, $\sigma_{S}$ and $\sigma_{R}$, as defined in Eq. (12). We have verified that $\mathcal{B}(B\to AA)\gtrsim 97\,\%$ in both topologies for all benchmarks.

In addition to the di-Higgs plus $\not{E}_{T}$ signature, our models induce processes like $pp\to B\to jj$ and $pp\to B\to\gamma\gamma$, as well as mono-jet, di-jet plus $\not{E}_{T}$, and mono-Higgs signatures. The highest sensitivity to our scenarios is expected in searches for top squarks, bottom squarks and gluinos Aad:2016eki ; Aaboud:2017ayj ; Aaboud:2017vwy ; ATLAS:2017dnw , which focus on large $\not{E}_{T}$ together with $b$-jets and light jets. Using CheckMATE2 Dercks:2016npn , we have verified that all benchmarks evade existing LHC searches for these and similar processes at $\sqrt{s}=8$ and $13\,\,\mathrm{TeV}$. Current searches for $4b+\not{E}_{T}$ lack sensitivity to our models, because they focus on phase-space regions that are sparsely populated by our signal.

We have also verified that future searches at the HL-LHC with $3\,\text{ab}^{-1}$ of data will not be competitive with our final state. Since the $14\,$TeV analyses available in CheckMATE2 are only a handful, we have estimated the reach of the remaining analyses by naively rescaling the expected reach of the current $13\,\,\mathrm{TeV}$ studies with the square root of the luminosity. We conclude that our benchmarks are not only viable today, but also they will not be probed by future LHC searches. Di-Higgs plus $\not{E}_{T}$ can thus be considered the discovery channel for these scenarios. From the resonant topology, we predict additional signatures with di-boson resonances and missing energy, like $WW+\not{E}_{T}$ and $ZZ+\not{E}_{T}$, as well as signatures with four electroweak bosons forming two di-boson resonances, cf. Section II. Since these signatures are expected to occur with similar rates as di-Higgs plus $\not{E}_{T}$ production, they can serve as complementary discovery channels for the resonant topology.

To illustrate the phenomenology of our simplified models, we use the specific benchmarks

[TABLE]

The first two benchmarks correspond to scenarios with little phase space for the $B\to AA$ decay (compressed spectrum), the second two with large phase space (split spectrum). In our benchmarks, a compressed (split) spectrum is parametrized by a small (large) $\Delta_{BA}$, which determines the boost of $A$. The kinematic differences between the two models originate in the respective $A$ decays, which are imprinted on the $\not{E}_{T}$ distribution. In the left panel of Figure 2, we show the $\not{E}_{T}$ distribution for the four benchmarks, illustrating the effects of a compressed and split spectrum in each topology. We also present the transverse momentum distributions of the Higgs bosons at parton level, $p_{T}(h)$, for a split spectrum in the symmetric topology (center panel) and the resonant topology (right panel).

In the symmetric topology, the two $A$ particles are produced back-to-back in the center-of-mass frame and split their transverse momentum into $\chi$ and $h$. The vector sum of their transverse momenta is thus subject to cancellations. The peak position of the $\not{E}_{T}$ distribution depends on the available phase space in both the $B$ and $A$ decays, $\Delta_{BA}=m_{B}-2m_{A}$ and $\Delta_{Ah\chi}=m_{A}-(m_{h}+m_{\chi})$. For larger values of $\not{E}_{T}$, the distribution drops fast. In the resonant topology, the $\not{E}_{T}$ distribution is equal to the transverse momentum distribution of $A$. The peak of the $\not{E}_{T}$ distribution depends now only on $\Delta_{BA}=m_{B}-2m_{A}$, and the spectrum is harder at large $\not{E}_{T}$ than for the symmetric topology. A trigger on missing energy will thus favor one or the other topology, depending on the position of the $\not{E}_{T}$ cut. The transverse momentum distribution of the Higgs bosons peaks at lower momenta than the $\not{E}_{T}$ distribution. As can be seen by inspecting the transverse momenta of the Higgs bosons, the $b$-jets from Higgs decays are less likely to pass the trigger requirements, which typically imply strong cuts on the jet transverse momentum Aaboud:2016leb . Triggering on missing transverse energy rather than on the $b$-jets yields a more efficient signal selection.

III.2 Jet substructure technique

Depending on the model and the mediator spectrum, the $b$-jets from Higgs decays can be produced with a large boost. The $b$-jets are thus collimated and cannot be resolved as individual jets. To reconstruct the boosted $h\to b\bar{b}$ decays, we crucially rely on jet substructure techniques. The current substructure module in Delphes, SoftDrop, is a modified version of the BDRS algorithm Butterworth:2008iy that includes $b$-tagging and flavor tagging for fat jets. To make the tool applicable for our purposes, we have extended these functionalities to subjets. Based on SoftDrop, we have developed two new modules called JetFlavorAssociationSubjets and BTaggingSubjets. 444The corresponding code can be obtained from the authors upon request. These modules allow us to access the four-momenta and $b$-tags of each fat jet in the event, and also of each subjet associated to it in the Delphes output. We will speak of “$x$-$y$ $b$-tags” to describe an event selection where one fat jet contains at least $x$ $b$-tagged subjets and another fat jet contains at least $y$ $b$-tagged subjets. The performance of our tagging technique depends on the fat-jet radius, $R$. We use $R=1.2$ for the symmetric and $R=0.6$ for the resonant topology. 555The choice of these values and their impact on the analysis will be explained at the end of this subsection.

Due to the limited $b$-tagging efficiency and rejection efficiency of light (i.e., non-$b$-tagged) jets, as well as the jet rapidity cut of $|\eta_{b}|<2.5$, not all of the four $b$-subjets in our signal will be tagged. To quantify this statement, we show in Figure 3 the number of $b$-tagged subjets, $N_{b}$, versus the number of light jets, $N_{j}$, for an exemplary benchmark of the resonant topology, $R\textunderscore 1000\textunderscore 275\textunderscore 25$. Other resonant benchmarks show a similar pattern, and the behavior is similar for the symmetric topology.

It is apparent that most of the time only three or fewer $b$-subjets are reconstructed. Note, however, that the amount of missing $b$-jets is larger than the naive estimation from the plain $b$-tagging efficiency. This loss is due to a significant number of reconstructed jets that are either soft, i.e., carry $p_{T}(j)<20\,\,\mathrm{GeV}$, and/or collinear, i.e., with angular separation of $\Delta R_{bb}<0.4$. 666While the detector can resolve jets and subjets with smaller angular separation, regular jets are clustered in ATLAS with a $k_{T}$ algorithm with $R=0.4$. In Figure 4, we show the $p_{T}$ distribution of the softest $b$-quark (left panel) and the minimum distance between any pair of $b$-quarks ($\Delta R_{bb}^{\rm min}$) in an event for the symmetric benchmark $S\textunderscore 1000\textunderscore 250\textunderscore 25$ (center panel), at the parton level (red) and after showering (blue). In the right panel of Figure 4, we show $\Delta R_{bb}^{\rm min}$ for the resonant benchmark $R\textunderscore 1000\textunderscore 275\textunderscore 25$.

From the left panel, we clearly see that the parton shower reduces the average transverse momentum of the $b$-quarks below the detector threshold of $20\,\,\mathrm{GeV}$. A tight event selection with $2$-$2$ $b$-tags would cut away a large amount of signal. We therefore apply looser requirements on the $b$-tags in our analysis (see Section III.3).

Now we turn our attention to the $\Delta R_{bb}^{\rm min}$ distribution. As expected, the parton shower barely changes the collinearity of the $b$-quarks. Therefore, if $b$-jets are not reconstructed as such, this is due to the characteristic mass spectrum of the model, rather than parton showering. The $\Delta R_{bb}^{\rm min}$ distribution depends both on the boost of the Higgs bosons and on the event topology. In a split spectrum, the Higgs bosons are more boosted, so that $b$-quarks from the same Higgs decay are closer to each other. In a compressed spectrum, the Higgs bosons are softer and the $b$-quarks are emitted with a larger angular separation. Naively one might thus think that smaller values of $\Delta R_{bb}^{\rm min}$ are preferred in the symmetric model, where the Higgs bosons carry larger transverse momenta. In Figure 4, however, we observe the opposite behaviour. This is due to the different event topology: In the resonant model, the two Higgs bosons stem from the decay of one $A$ boson and are thus much closer than in the symmetric model, where they originate from opposite sides of the event. Consequently in the resonant model all four $b$-jets tend to be collimated, while in the symmetric model the two pairs of $b$-jets are well separated. We hence conclude that only in the symmetric model $\Delta R_{bb}^{\rm min}$ is a direct measure of the Higgs transverse momentum. In the resonant model, on the other hand, the closest $b$-jets do not always stem from the same Higgs decay, so that $\Delta R_{bb}^{\rm min}$ is also sensitive to the boost of $A$. In the symmetric model, the parton level requirement $\Delta R_{bb}>0.4$ only cuts away a few percent of the signal events. In contrast, in the resonant model this cut has an important impact on the signal. This loss of events, together with the lower total event rates discussed in Section II, suggests that the resonant topology is harder to find that the symmetric one for a given particle spectrum. The fact that the maximum of the $\Delta R_{bb}^{\rm min}$ distribution in the resonant model lies at lower values than in the symmetric model motivates different choices of fat-jet radii. We use $R=1.2$ for the symmetric and $R=0.6$ for the resonant topology.

III.3 Backgrounds and cutflow analysis

The main SM backgrounds to our signal of $4b+\not{E}_{T}$ are due to $W+\text{jets}$, $Z+\text{jets}$, as well as top-antitop production with one leptonic and one hadronic decay. All backgrounds have been simulated using Sherpa 2.2.1 Gleisberg:2008ta at leading order (LO) in QCD, including parton shower and hadronization effects. We use the same setup as for the signal generation, as described at the end of Section III.1. Our analysis has been performed with ROOT Brun:1997pa .

The cutflow analysis is summarized for the symmetric benchmark $S\textunderscore 750\textunderscore 350\textunderscore 25$ in Figure 5 and for the resonant benchmark $R\textunderscore 750\textunderscore 350\textunderscore 25$ in Figure 6. In the right panel of each figure, we list the cross section for the dominant background processes for $\not{E}_{T}>200\,\,\mathrm{GeV}$ and after applying a lepton veto and requiring at least two $b$-tagged subjets from fat jets with radius $R=1.2$ (symmetric topology) and $R=0.6$ (resonant topology), respectively. We have checked that contributions from di-boson plus jets production are smaller. The latter will be neglected in our analysis.

To discriminate between our signal and the backgrounds, we apply a cut-and-count procedure. Throughout our analysis, we apply an initial cut of $\not{E}_{T}>200\,\,\mathrm{GeV}$ and a lepton veto. To study the impact of our $b$-tagging technique, we request various $x$-$y$ $b$-tags with $0\leq x,y\leq 2$, one by one.777The upper limit on $x,y$ is due to the fact that we only consider the two hardest subjets within a given fat-jet. We furthermore apply an optional Higgs Mass Window (HMW) by requesting that the mass of each of the two identified fat jets lies within the window $75\,\,\mathrm{GeV}<m_{J}<175\,\,\mathrm{GeV}$. In addition, we allow for a variable lower cut on $\not{E}_{T}$ (in steps of 50 GeV).888For the symmetric topology, a looser $\not{E}_{T}$ cut is preferred to optimize the significance. For the resonant topology, due to the harder $\not{E}_{T}$ spectrum and the low number of signal events, the preferred cut lies at higher $\not{E}_{T}$. In order to establish a fair comparison of the remaining selection criteria in the two topologies, we discuss the cutflow for a fixed cut of $\not{E}_{T}>200\,\,\mathrm{GeV}$. To optimize the choice of the jet radius, we have carried out our analysis for four different fat-jet radii $R=0.6,\,0.8,\,1.0,\,1.2$. For the symmetric topology, the significance is maximized for $R=1.2$, while the resonant model favors a smaller fat-jet radius of $R=0.6$.999For the resonant model, the sensitivity with $R=0.6$ is a factor of 2 higher than for $R=1.2$.

The impact of the various cuts on signal and background is shown for the symmetric model in the left panel of Figure 5. We see that applying $2$-$2$ $b$-tags plus a Higgs mass window leaves us with about 40 signal events, while the sum of backgrounds ranges around 700 events. In the cut-and-count analysis, we therefore do not achieve a signal-to-background ratio of ${\cal O}(1)$, so that the significance depends critically on systematic uncertainties that can affect the analysis. In the center panel, we illustrate this dependence by showing the significance defined as

[TABLE]

We assume a systematic uncertainty of $\beta=0,1,5,10\%$, respectively. For a small uncertainty $\beta=1\%$, a maximum significance of about $\Sigma=2$ can be reached for the considered benchmark. In the resonant model, shown in Figure 6, the signal rate is significantly lower than in the symmetric model. With our basic cut-and-count analysis, we therefore do not achieve a noticeable sensitivity to our signal. The fact that the significance of our signal depends on a combination of various kinematic variables suggests to perform a multi-variate analysis to optimize the sensitivity.

IV Multi-variate analysis and results

In this section, we describe the strategy pursued in our multi-variate analysis (MVA) and present our results. Depending on the respective phase-space point, discriminating between the di-Higgs plus $\not{E}_{T}$ signal and the backgrounds can be very challenging. In order to maximize the sensitivity for each of our benchmarks in the two models, we perform a multi-variate analysis. In this study, we use the scikit-learn Pedregosa:2012toh implementation of AdaBoost Freund:1997xna , employing the SAMME.R algorithm to perform a Boosted Decision Tree (BDT) classification. As our best setup, we choose to train 70 trees with a maximal depth of 3, a learning rate of 0.5 and a minimum node size of 0.025 of the total weights.

Before running the BDT, we place basic kinematic selection cuts on the missing transverse energy and the jets. As in the cutflow analysis, we apply a $\not{E}_{T}>200\,\,\mathrm{GeV}$ cut and veto events containing isolated leptons.101010We require electrons (muons) to have $p_{T\text{lep}}>10\,(10)$ GeV and $|\eta_{\text{lep}}|<2.5(2.7)$, and consider them isolated if $\frac{p_{T\text{lep}}}{p_{T\text{had}}}<0.12\,(0.25)$ within $R=0.5\,(0.5)$ of the electron (muon) momentum. Looser lepton selection criteria might result in a better rejection of the $W$+jets background and therefore improve the significance of our analysis. However, in this work we consistently use our conservative lepton definition. The jets are defined as Cambridge-Aachen fat jets $J$ with a jet radius of $R=1.2$ ($R=0.6$) for the symmetric (resonant) topology and transverse momentum $p_{T}(J)>20\,\,\mathrm{GeV}$. A fat jet $J_{i}$ is accepted if it contains two subjets $j_{i}^{k}$, $i,k=\{1,2\}$, where at least one of them is $b$-tagged.

In our multi-variate analysis, we use kinematic information on the two hardest $b$-tagged fat jets, $J_{1}$ and $J_{2}$ and their corresponding subjets $j_{1}^{k}$ and $j_{2}^{k}$. The complete set of variables used for our analysis can be classified in four categories:

• •

Global variables: missing transverse energy, $\not{E}_{T}$; $H_{T}$ (computed using the fat jets); total number of fat jets, $N_{J}$; total number of $b$-tagged fat jets, $N_{Jb}$; total number of $b$-tagged subjets within all fat jets, $N_{jb}$;

• •

Single fat-jet variables: transverse momentum $p_{T}(J_{i})$; pseudo-rapidity $\eta(J_{i})$; jet mass $m_{J_{i}}$; azimuthal angular separation between fat jet and missing momenta, $\Delta\phi(J_{i},\not{E}_{T})$; ratio of transverse momenta $p_{T}(J_{i})/\not{E}_{T}$;

• •

Two fat-jet variables: distance between two fat jets, $\Delta R(J_{1},J_{2})$; invariant mass of two fat jets, $m_{J_{1}J_{2}}$; maximum jet mass ratio, max($m_{J_{1}}/m_{J_{2}},m_{J_{2}}/m_{J_{1}}$);

• •

Subjet variables: transverse momentum $p_{T}(j^{k}_{i})$; pseudo-rapidity $\eta(j^{k}_{i})$; distance between subjets, $\Delta R(j^{k}_{i},j^{l}_{i})$.

We employ $80\%$ of our events for training and $20\%$ for evaluation purposes. The different backgrounds are weighted according to their relative cross section after applying the basic selection cuts. The BDT thus focuses on the dominant backgrounds when trained to avoid misidentification of the respective backgrounds as a signal. To make sure that the BDT will take equal effort in correctly classifying the overall number of signal events and background events, we scale the total weight of all signal events to match the total weight of all background events. This is especially important, since in the training we involve more Monte Carlo events for the background processes than for the signal. The BDT will assign a score (or threshold in machine learning (ML) terminology) to each event, which reflects the likelihood of it being signal.

In Figure 7, we show a typical BDT result for the symmetric benchmark $S\textunderscore 750\textunderscore 350\textunderscore 25$.

The three plots to the left show the training results, while the three plots to the right display the outcome of the evaluation. In the upper left plot of each panel, we present the Receiver Operating Characteristic (ROC) curve. Shown is the true positive rate (TPR, also referred to as recall in ML language) against the false positive rate (FPR, also referred to as fall-out). 111111The TPR is the probability that a signal event gets tagged as signal, while the FPR is the probability that a background event gets tagged as signal. The most relevant information is the Area Under the Curve (AUC), which quantifies the BDT capability to discriminate between signal and background. By construction, the AUC ranges between 0.5 and 1. For all our benchmarks, the AUC is at or above $0.9$, which proves that our signal/background classifier has an impressive performance.

The upper right plot of each panel displays the precision (or positive predictive value) as a function of the recall. The precision is defined as the fraction of true signal events among those events the BDT classified as signal. This curve illustrates how reliable a classification as signal is, depending on which fraction of signal events are classified correctly. In the displayed curve, we used the adjusted signal weights, ensuring that signal and background are on equal footing (see the discussion above). Hence the minimum value for the precision is 0.5.

The lower plot in each panel shows the significance $\Sigma$ defined in Eq. (18) as a function of the expected number of signal events $S$ that would be left after cutting on a given score. The prediction is made for the HL-LHC with an integrated luminosity of $L=3\,\text{ab}^{-1}$. We present the significance for this benchmark for three different assumptions of systematic uncertainties. It is apparent that the use of a BDT enhances the significance by about an order of magnitude compared to our basic cut-and-count analysis.121212The small spikes in the significance curve in the evaluation sample are due to a lack of Monte Carlo statistics in the $Wjjjj$ background for high BDT scores. We have simulated $10^{8}$ Monte Carlo events for this background. Owing to the large number of colored final states, however, the event generation is computationally intense. Moreover, the lepton veto cannot be reliably implemented at parton level, which enhances the number of required events. Generating a larger sample of $Wjjjj$ events would soften the spikes, but is beyond the scope of our analysis.

In order to determine the sensitivity to a given benchmark scenario, we cut the BDT score (in our case, the number of signal events) at the peak of the evaluation significance. In our example $S\textunderscore 750\textunderscore 350\textunderscore 25$ and assuming $5\%$ systematics, this corresponds to a significance of $\Sigma=12$ and $S=325$ and $B=238$ remaining signal and background events, compared to $S=42$ and $B=713$ in the basic cut-and-count approach. The increase in sensitivity is due to a better selection of signal events, while at the same time having a similar improvement in background rejection. All of our benchmarks in the symmetric model feature a signal-to-background ratio close to unity, which suggests good discovery prospects in the presence of systematic uncertainties. To take into account the statistical uncertainty on our $Wjjjj$ background simulation (see footnote 12), in what follows we take a conservative approach and claim a “discovery” at a significance of $\Sigma=7$ (corresponding to 7 standard deviations from the standard-model hypothesis for Gaussian statistics), instead of the common $5\sigma$ threshold.

In Figure 8, we present our results for the resonant benchmark $R\textunderscore 750\textunderscore 350\textunderscore 25$.

Again, we see the excellent performance of our BDT classifier. However, in the resonant model the lower signal rate severely limits the sensitivity. While the BDT improves the sensitivity, we cannot reach the discovery level for our models. Still, a significance of $\Sigma\approx 2$ can be achieved in benchmarks with a lighter scalar $B$, corresponding to a larger signal rate. The HL-LHC can thus test parts of the parameter space, but a more refined strategy (or a combination of multiple channels) would be required to reach a higher sensitivity. In summary, compared with the cut-and-count analysis the BDT enhances the sensitivity to the resonant topology by about a factor of 10 in most benchmarks. Still, the sensitivity is lower than for the symmetric topology, mostly due to the reduced signal rate.

To show the dependence of the signal sensitivity on the respective model parameters, we present our results in terms of the scalar masses $m_{A}$ and $m_{\chi}$. In Figure 9, we display the luminosity required to discover the symmetric benchmark scenarios from Table 1 at the HL-LHC, assuming a systematic uncertainty of $5\%$. The mass of the heavy scalar is fixed to $m_{B}=500\,\,\mathrm{GeV}$ (left panel) and $m_{B}=750\,\,\mathrm{GeV}$ (right panel).

Apart from one benchmark, all scenarios are well within the reach of the HL-LHC. We also see that the significance is particularly high for benchmarks with a compressed spectrum. As anticipated from Figure 2, for $m_{B}\gtrsim 2m_{A}$, the $\not{E}_{T}$ spectrum is harder and the cut on missing energy is thus more efficient in rejecting the background. In Figure 9, this effect can be seen by looking at fixed values of $\Delta_{Ah\chi}$ and increasing $m_{A}$ (along the diagonal). For a more compressed spectrum, the required luminosity is drastically reduced by a factor of 10-20, depending on the actual value of $m_{B}$. A similar but milder effect occurs if the decay $A\to h\chi$ proceeds close to threshold. In this case, the $h\chi$ pair follows the direction of $A$, resulting again in a harder $\not{E}_{T}$ distribution. In the plot, this corresponds to fixed values of $m_{A}$ and increasing $m_{\chi}$ (along the vertical direction), resulting in a reduction of the required luminosity by a factor of about 2-3 at most. In summary, a di-Higgs signal could be discovered in an early phase of the HL-LHC, provided that the scalar resonance $B$ is produced at a sizeable rate.

We present our results in a second way, which is particularly convenient for recasting purposes. In Figure 10, we report the cross section that can be probed at the discovery level for a luminosity of 3 ab*-1*. Again, we fix the heavy scalar mass at $m_{B}=500\,\,\mathrm{GeV}$ (left panel) and $m_{B}=750\,\,\mathrm{GeV}$ (right panel), respectively.

We see that in the most difficult benchmark topology, namely in the benchmark that requires the largest discovery luminosity at the HL-LHC, $S\textunderscore 750\textunderscore 275\textunderscore 50$, a cross section of about 15 fb would be required to claim a discovery. It is interesting to compare this value with the latest di-Higgs predictions for the SM that require a total rate for $hh\to 4b$ of about 13 fb for discovery. In the standard model, the Higgs pair is not produced through a resonance, and furthermore the final state of four $b$-quarks without $\not{E}_{T}$ is difficult to identify. In contrast, in our scenarios the scalar $B$ is resonantly produced and the large $\not{E}_{T}$ in the final state facilitates an efficient discrimination against the backgrounds. The fact that we can probe cross sections of a few to several femtobarns is an important result, which motivates a study of the di-Higgs plus $\not{E}_{T}$ signature in the context of complete models. For the benchmarks with $m_{B}=1\,\,\mathrm{TeV}$, the signal rates are too low to claim a discovery at the HL-LHC with $3$ ab*-1*. However, already with slightly more luminosity (or larger cross section) a discovery of these heavy scalar scenarios is possible.

In the resonant model, the planned HL-LHC luminosity is not sufficient to discover any of the benchmark scenarios, due to the lower production rates. We therefore confine ourselves to presenting the cross sections required to discover a particular resonant benchmark in Figure 11. As explained in Section III, the mass of the lightest scalar, $m_{\chi}$, does not affect the sensitivity, since the boost of $A$ does not depend on $m_{\chi}$ or $m_{h}$. We therefore present our results in terms of the heavier scalar masses $m_{A}$ and $m_{B}$.

From the figure, we see that we can only test cross sections in the femtobarn and sub-femtobarn regime. As in the symmetric model, the significance increases when the spectrum is compressed. Increasing $m_{A}$ for fixed $m_{B}=1000\,\,\mathrm{GeV}$ lowers the testable cross section by a factor of $3$. Using CheckMATE2 we have verified that even with the largest possible cross section displayed here, the search for di-Higgs plus $\not{E}_{T}$ is still the most sensitive channel for the resonant model.

V Dark matter

In our models, the lightest scalar $\chi$ is automatically stable, due to the imposed $\mathbb{Z}_{2}$ symmetry. Here we explore the hypothesis that $\chi$ is a dark matter candidate. We discuss the dark matter phenomenology for our simplified models, focusing on dark matter-nucleon scattering in direct detection experiments and thermal freeze-out in the early universe.

A contribution to spin-dependent nucleon scattering arises from the portal operator $H^{\dagger}H\,\chi\chi$, which induces a $h\chi\chi$ interaction after electroweak symmetry breaking. The effect of this operator has been well studied elsewhere (see e.g. Ref. Athron:2017kgt for a recent study). Since it is not relevant for our di-Higgs plus $\not{E}_{T}$ final state, we will assume that it is absent. We follow the same philosophy for other operators, namely we only study the implications of operators that play a role in the collider phenomenology.

In the symmetric model, $A$ and $\chi$ can mix through the coupling $\lambda_{A\chi HH}$ from Eq. (4) upon electroweak symmetry breaking. This mixing induces a Higgs coupling to the lightest scalar, $\lambda_{A\chi HH}\sin\theta_{A\chi}$, where $\theta_{A\chi}$ is the $A$-$\chi$ mixing angle. Direct detection experiments set a strong bound on this coupling. However, as we explained in Section II, neither $\lambda_{A\chi HH}$ nor the mixing affects the signal strength of di-Higgs plus $\not{E}_{T}$. We have therefore set $A$-$\chi$ mixing to zero in our analysis, $\theta_{A\chi}=0$. The lightest scalar $\chi$ can be a viable dark matter candidate that leaves a di-Higgs plus $\not{E}_{T}$ signature at the LHC while agreeing with the null results from direct detection experiments.

In the resonant model, $A$ and $h$ mix through the operator $m_{AHH}AH^{\dagger}H$ after electroweak symmetry breaking. For fixed $m_{A}$ and $m_{\chi}$, direct detection experiments set an upper bound on the product of the $m_{A\chi\chi}$ and $m_{AHH}$ couplings from mixing-induced nucleon scattering. Since the signal strength of di-Higgs plus $\not{E}_{T}$ production depends on the relative size of $m_{AHH}$ and $m_{A\chi\chi}$, but not on their overall magnitude (see Section II), we conclude that in the resonant model the lightest scalar $\chi$ is not ruled out as a dark matter candidate by direct detection experiments if $m_{A\chi\chi}m_{AHH}$ is sufficiently small.

Assuming that our dark matter candidate is a relic from thermal freeze-out in the early universe sets additional constraints on the parameter space of our models. In the symmetric model, dark matter annihilation can be efficient in either of the following scenarios,131313Additional annihilation processes such as $\chi\chi\to B\to gg$ via the $C_{Bgg}$ coefficient are possible. However, such channels would also require a $B\chi\chi$ interaction, which is irrelevant for the di-Higgs plus $\not{E}_{T}$ signature. As in the case of direct detection, we will neglect these interactions.

[TABLE]

If dark matter is heavier than the Higgs boson, it can annihilate by $t$-channel $A$ exchange, scaling as $\lambda_{A\chi HH}^{4}$. The observed dark matter abundance of $\Omega_{\chi}h^{2}=0.1199$ Ade:2015xua can be obtained with moderate couplings and mediator masses $m_{A}\lesssim 1\,\,\mathrm{TeV}$. If dark matter is lighter than the Higgs boson, it can only annihilate through $A$-$\chi$ mixing, which is strongly suppressed by the null results from direct detection experiments. It is still possible to satisfy the observed relic abundance for dark matter pair masses near the Higgs mass. In this case, $s$-channel annihilation through the Higgs boson occurs resonantly, which compensates for the coupling suppression.

In the resonant model, the direct detection bounds on $m_{A\chi\chi}\sin\theta_{Ah}$ suppress all interactions of dark matter. The observed relic abundance can only be obtained for dark matter pairs in the Higgs resonance region,

[TABLE]

Obtaining the observed relic abundance away from the Higgs resonance requires additional dark matter annihilation channels beyond what is predicted in our simplified model. In any case, the dark matter hypothesis should not constrain the search strategy for di-Higgs plus $\not{E}_{T}$ at the LHC. For instance, a di-Higgs plus $\not{E}_{T}$ signature could also arise in models with hidden sectors, where $\chi$ decays visibly at a later time and outside the detector, so that its decay products could be caught by dedicated detectors such as FASER Feng:2017uoz ; Ariga:2018uku or MATHUSLA Evans:2017lvd ; Curtin:2018mvb ; Alpigiani:2018fgd .

VI Conclusions

In this work we have developed a search strategy for $hh\chi\chi$ production at the LHC using the final state with four $b$-jets and missing transverse energy. For the purpose of our and future studies, we have built two simplified models that give rise to this final state in different kinematic topologies. Both models feature a hidden sector of three new scalar singlets $A$, $B$ and $\chi$, the latter being stable due to an imposed parity symmetry. Since the scalars $A$ are pair-produced via a resonant scalar $B$, event rates in both models are sizeable at the HL-LHC. The decay of $A$ depends on the properties of the particles in the hidden sector and determines the event topology. In our symmetric model, both $A$ and $\chi$ belong to the dark sector, so that $A\to h\chi$ decays occur on both sides of the decay chain. In the resonant model, only $\chi$ belongs to the hidden sector and $A$ decays via $A\to hh$ or $A\to\chi\chi$ with similar branching ratios. We stress that these simplified models can be embedded in more complete models featuring an enlarged scalar sector or other particles in the hidden sector.

To demonstrate the LHC potential to discover hidden sectors with the di-Higgs plus $\not{E}_{T}$ signature, we have performed a full-fledged numerical analysis of the multi-$b+\not{E}_{T}$ final state, including a detailed study of SM backgrounds and detector effects. Dominant backgrounds are due to $t\bar{t}$, as well as $Wjjjj$ and $Zjjjj$ production. Employing the inclusive $\not{E}_{T}$ trigger for event pre-selection, we have first carried out a cut-and-count analysis, followed by a multi-variate analysis based on a boosted decision tree. Both analyses rely on the use of jet-substructure techniques in a modified version of the BDRS algorithm. In our cut-and-count analysis, we employ missing energy, jet and subjet variables, as well as flavor tags to efficiently discriminate between signal and background. With this approach, we have obtained significances of at most $2\sigma$ for the symmetric model, while in the resonant model our analysis turned out to be insensitive. It is thus necessary to optimize the signal-background discrimination by optimally exploiting all kinematic features using a multivariate analysis. In particular we have shown that the use of machine-learning tools is well suited for our analysis.

In the symmetric model, our BDT analysis predicts a significance well above $5\sigma$ for most of our scanned points, thus opening the possibility of an (early) discovery at the HL-LHC. In the resonant model, the signal rates are significantly lower (in the sub-femtobarn regime), which reduces the sensitivity. In order to claim a discovery, an increase of our benchmark cross section by a factor of 2-10 would be needed. In any case, the enhancement from our simple cut-and-count to the multivariate analysis shows that the sensitivity to the di-Higgs plus $\not{E}_{T}$ signal relies on a variety of kinematic features in both signal and background. The BDT is thus the appropriate approach to search for such a many-body final state in an environment with large SM backgrounds. While we have focused on the $4b+\not{E}_{T}$ channel with the largest event rates, additional final states like $b\bar{b}\,\gamma\gamma+\not{E}_{T}$, $b\bar{b}\,WW^{*}+\not{E}_{T}$, $b\bar{b}\,\tau^{+}\tau^{-}+\not{E}_{T}$ can contribute significantly to enlarge the search potential of di-Higgs plus $\not{E}_{T}$.

In the context of dark matter, it is interesting to investigate the interplay of this collider signature with searches at direct and indirect detection experiments in complete models where the relic abundance is satisfied. Within our simplified models, parts of the parameter region could provide us with a viable thermal dark matter relic. In general, potential bounds on viable dark matter models should not limit the scope in di-Higgs plus $\not{E}_{T}$ searches at the LHC.

So far, the LHC collaborations have searched for the di-Higgs plus $\not{E}_{T}$ signature in a specific scenario of supersymmetry with very light invisible particles in the final state. We encourage our experimental colleagues to use our simplified models and a search strategy similar to ours to fully exploit the discovery potential of di-Higgs plus $\not{E}_{T}$. The use of machine learning techniques is crucial to achieve a high significance, in our case an enhancement by an order of magnitude over a cut-and-count analysis. Searching for di-Higgs plus $\not{E}_{T}$ links the efforts at the dark matter frontier with those on the di-Higgs frontier, which in the last few years have seen a spectacular development in both theory and experiment. The sensitivity to this and similar signatures will greatly benefit from merging the techniques developed for Higgs pair measurements in and beyond the standard model with missing energy searches.

VII Acknowledgements

We thank Martin Bauer, Ulrich Haisch, José Miguel No and Tilman Plehn for helpful discussions. This project has been realized within the Heidelberg Karlsruhe Strategic Partnership HEiKA. We acknowledge HEiKA funding through the research bridge “Particle Physics, Astroparticle Physics and Cosmology (PAC)” that facilitated our collaboration. The research of MB and SW is supported by the DFG Collaborative Research Center TRR 257 “Particle Physics Phenomenology after the Higgs Discovery”. SK acknowledges the support of the DFG-funded Doctoral School “Karlsruhe School of Elementary and Astroparticle Physics: Science and Technology” (KSETA). SW acknowledges funding by the Carl Zeiss foundation through an endowed junior professorship (Junior-Stiftungsprofessur).

Appendix A UV completing the $Sgg$ coupling

In this appendix, we introduce a minimal perturbative UV completion for the effective coupling $C_{Bgg}$ in Eq. (10). We add a heavy quark $Q$ with mass $m_{Q}$, odd $\mathbb{Z}_{2}$ parity, and vector-like weak interactions to our model. We furthermore assume that $Q$ couples to the heavy scalar $B$ via

[TABLE]

For the sake of simplicity, we take $Q$ to be an $SU(2)$ singlet with hypercharge $-1/3$. Assuming that $Q$ dominantly decays via $Q\to b\chi$, its mass is constrained by sbottom searches at the LHC Aaboud:2017wqg ; Sirunyan:2017kiw , as well as by more inclusive searches for jets plus $\not{E}_{T}$ Aaboud:2017vwy ; Sirunyan:2017cwe . We estimate the current bound to be at the level of $1\,\,\mathrm{TeV}$, but leave a more thorough investigation for future work.

Following Ref. Falkowski:2015swt , integrating out $Q$ generates the following effective couplings to gluons and photons, 141414If we had instead introduced $Q$ as a vector-like top partner, the coupling $C_{B\gamma\gamma}$ would be larger by a factor of four.

[TABLE]

where $g_{s}$ is the coupling constant of QCD and $C_{Bgg}$ is defined in Eq. (10). The effective coupling to photons is defined analogously by

[TABLE]

With $m_{Q}=1\,$TeV and $y_{Q}=1$, we find151515Notice that the naive dimensional analysis (NDA) estimate for $C_{Bgg}$, assuming a loop-induced perturbative UV completion at the scale $\Lambda=1\,$TeV, is larger by about a factor of three, $C_{Bgg}\sim 1/(16\pi^{2})$.

[TABLE]

and $C_{B\gamma\gamma}$ smaller by more than two orders of magnitude. Both $C_{Bgg}/\Lambda$ and $C_{B\gamma\gamma}/\Lambda$ scale as $y_{Q}/m_{Q}$, so that the effective couplings do not change when simultaneously increasing both $m_{Q}$ and $y_{Q}$. The partial decay widths mediated by these couplings are Falkowski:2015swt

[TABLE]

This should be compared with the decay width into $AA$ pairs,

[TABLE]

The production cross section $\sigma(pp\to B)$ can be estimated by making use of the results of the LHC Higgs cross-section working group deFlorian:2016spz , which provides the contribution of gluon-gluon fusion to the production cross section of a heavy scalar $\hat{S}$ with Higgs-like couplings to quarks. These numbers can then simply be rescaled to obtain the production cross section of $B$,

[TABLE]

where $m_{t}$ and $y_{t}$ are the top-quark mass and Yukawa coupling.