Collider Searches for Dark Matter through the Higgs Lens

Abstract

Despite the fact that dark matter constitutes one of the cornerstones of the standard cosmological paradigm, its existence has so far only been inferred from astronomical observations, and its microscopic nature remains elusive. Theoretical arguments suggest that dark matter might be connected to the symmetry-breaking mechanism of the electroweak interactions or of other symmetries extending the Standard Model of particle physics. The resulting Higgs bosons, including the 125 GeV spin-0 particle discovered recently at the Large Hadron Collider, therefore represent a unique tool to search for dark matter candidates at collider experiments. This article reviews some of the relevant theoretical models as well as the results from the searches for dark matter in signatures that involve a Higgs-like particle at the Large Hadron Collider.

1. Introduction

The concept of dark matter (DM) was originally introduced to reconcile the observations of the high velocity dispersion of galactic clusters [1] and the flat rotational curves of spiral galaxies [2] with the predictions of Newton’s gravity. Since then, precise cosmological and astrophysical observations [3] have strengthened the evidence for the existence of DM, establishing that around 25% of the energy budget of the observable universe consists of a matter component that in the standard DM picture is electrically neutral, weakly interacting and nonrelativistic.

Since the Standard Model (SM) of particle physics does not provide any candidate particle with the above characteristics, the existence of DM constitutes evidence for physics beyond the SM (BSM). The detection and identification of the nature of DM at the microscopic level constitutes therefore one of the major challenges for particle physics, with complementary searches pursued by direct detection (DD), indirect detection (ID) and collider experiments. The former two types of searches rely on the observation of recoils from the elastic scattering of DM particles on the nuclei in the detector material or on the observation of annihilation products of DM pairs such as monochromatic photons. For recent reviews of DD and ID search strategies, see for instance [4,5]. Collider searches for DM on the other hand rely on the production of DM particles in high-energy particle collisions and can be separated in two broad classes according to the experimental signature that they produce: (i) searches for missing transverse momentum ($E_T^{\mathrm{miss}}$) plus X signatures, also known as mono-X, where the $E_T^{\mathrm{miss}}$ resulting from the DM particles leaves the detectors unnoticed and the visible, i.e., detectable, final state X is used for triggering, and (ii) searches containing only visible particles such as pairs of leptons or jets that aim to detect the particles mediating the interactions between the DM and the SM particles through the observation of a new resonance or a modification of the kinematics of the final-state particles. Both types of searches have been pursued at the Large Hadron Collider (LHC) since its very beginning and previously at the Tevatron. See recent reviews [6,7,8] for general overviews on DM collider phenomenology.

The discovery of a SM-like Higgs boson by the ATLAS and CMS collaborations [9,10] has opened up a new avenue in the searches for DM, allowing one to probe the possible connection between the Higgs boson and the dark sector, i.e., a BSM sector that contains the DM particle and is almost decoupled from the SM. In fact, there are both experimental and theoretical arguments that suggest that the Higgs boson partakes in the mediation between the dark and the visible sector. Experimentally, the Higgs sector is compared to the gauge or fermionic sector of the SM far less explored and constrained, while theoretically, the SM Higgs doublet plays a special role, because it is the only SM field that allows one to write down a renormalisable coupling to the dark sector, if DM is uncharged under the SM gauge group. The Higgs sector therefore offers an interesting portal to the hidden sector [11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29] and dark sectors in particular. The goal of this review is to discuss experimental and theoretical aspects of BSM models that feature a DM–Higgs connection.

Given that the microscopic nature of DM is essentially unknown, there is a lot of freedom in the description of the interactions between the dark and the visible sector, leading to a vast array of phenomenological models. The models discussed in this review can be broadly classified into four categories according to the structure of the dark and visible sector and the type of particle(s) mediating between the two sectors:

(i)

Models in which the SM-like Higgs boson itself mediates between the dark sector and the SM. This model class represents the simplest realisation of the Higgs portal idea with a minimal particle content, containing, besides the SM states, only a single DM field that can be of spin-0, spin-$1/2$ or spin-1. The DM–SM interactions can be formulated in terms of composite operators of dimension four and higher in an effective field theory (EFT) framework, which allows one to describe the DD, the ID and the collider phenomenology in a model-independent fashion. In Section 2, we discuss in detail the simplest realisation of these Higgs portal EFTs, namely the case of a real singlet scalar field, considering both a marginal dimension-four and a derivative dimension-six coupling. Special emphasis is thereby put on highlighting the complementary of the different noncollider and collider search strategies in constraining the parameter space of the two models. We also briefly discuss the case of the Higgs portal model with fermionic DM that is a benchmark model used by both ATLAS and CMS to interpret their invisible Higgs decay searches.

(ii)

Models with an extended Higgs sector in which a spin-0 particle that mixes with one of the visible non-SM Higgs bosons mediates between the dark and the visible sector. Our discussion of this class of models is presented in Section 3 and focuses on the 2HDM+a [30,31,32,33] and the 2HDM+s [34,35] models. These two models are the simplest gauge-invariant and renormalisable models in this class that feature a fermionic DM candidate and therefore represent the natural extension of the simplified pseudoscalar and scalar DM models, as defined in [36,37]. The particle content of these models involves four additional BSM spin-0 states, besides the DM candidate. The extra spin-0 states have an important impact on the collider phenomenology of the 2HDM+a and 2HDM+s models, as they allow for resonant mono-Higgs, mono-Z and $tW+E_T^{\mathrm{miss}}$ production, thereby leading to a far richer mono-X phenomenology at the LHC when compared to the spin-0 simplified DM models. We review the state of the art of the phenomenology and experimental constraints on the 2HDM+a and 2HDM+s models, pointing out, in particular, the similarities and differences between the two models.

(iii)

Models with extended Higgs and gauge sectors in which both a spin-0 and a spin-1 portal connect the dark and the visible sector. The resulting theories fall into the class of dark Higgs or dark $Z^{\prime }$ models, which typically have a rich collider and DM phenomenology. In Section 4, we discuss two representative models in more detail: the 2HDM+$Z^{\prime }$ model [38] and the dark Higgs two-mediator DM (2MDM) model [39], which were used to guide and interpret existing mono-Higgs searches by ATLAS and CMS. While the scalar sectors of this models are quite different, both the 2HDM+$Z^{\prime }$ and the 2MDM models contain a spin-1 mediator that couples the DM particles to some of the SM states. This feature allows one to the test the models by searching for resonant production of non-$E_T^{\mathrm{miss}}$ final states such as dijets, $t\overline{t}$ and $Zh$. Using the latest available LHC data, we derive the constraints on the 2HDM+$Z^{\prime }$ and the 2MDM model that arise from the relevant searches for resonant SM final states and the mono-jet signature. In both cases, we show that for the benchmark scenarios considered by ATLAS and CMS, non-$E_T^{\mathrm{miss}}$ searches exclude additional parameter space not probed by the existing mono-Higgs interpretations.

(iv)

Models that feature exotic decays of the $125\mathrm{GeV}$ Higgs boson into hidden sector particles with macroscopic proper decay lengths of $c\tau >{10}^{-2}\mathrm{m}$. Due to the long-lived nature of the hidden particles, their decays into SM particles are displaced from the primary interaction vertex, which represents a striking experimental signature. In Section 5, we discuss three representative models that feature such long-lived particles (LLPs). These models derive either from the idea of neutral naturalness [40,41,42,43,44] or involve a hypercharge portal [45,46] or a hypercharge and a fermion portal [47,48,49,50,51]. While the spin of the particle that connects the hidden and the visible sector differs, all three models share a common feature: they can lead to both prompt and displaced signatures, depending on the strength of the interaction that links the two sectors. ATLAS, CMS and, in some cases, also LHCb have already searched for exotic Higgs decay signals, and we compare the available findings in a set of summary plots. Whenever possible, we include in these plots other experimental limits as well, to emphasise the unique role that searches for LLPs can play in constraining the parameter space of the considered theories.

Each of the four sections mentioned above is structured in a similar fashion. We first present the most relevant theoretical aspects of the considered models and then discuss the experimental constraints that apply in each case, focusing in many but not all cases on the collider bounds. The relevant constraints are combined into state-of-the-art summary plots in various benchmark scenarios of the examined DM models. Whenever several results that address a particular signature for a given model are available, we focus on the most recent measurements that typically provide the highest sensitivity. Hence, most ATLAS and CMS searches presented here are based on the full LHC Run 2 data set of around $140{\mathrm{fb}}^{-1}$ collected at $13\mathrm{TeV}$. For the readers mostly interested in collider phenomenology,

Table 1. The list of LHC signatures most relevant to this review, the models that they constrain, and the section(s) and/or figures(s) where the corresponding discussion can be found in the manuscript. The LHC signatures are grouped by double lines into three classes: (i) processes with a significant amount of $E_T^{\mathrm{miss}}$, (ii) prompt signals involving SM final states and (iii) signatures relevant in the context of LLP searches. In each class, the signatures are ordered as they appear in the text.

Signature	Model	References
$h\to \mathrm{inv}$	EFT Higgs portals	Section 2.2 and Section 2.5, Figure 3 and Figure 4
	2HDM+a	Section 3.1
	Neutral naturalness	Section 5.1, Figure 20
	Dark photons	Section 5.2, Figure 22 and Figure 23
	Vector plus fermion portal	Section 5.3, Figure 26
$\mathrm{VBF}+E_T^{\mathrm{miss}}$	EFT Higgs portals	Section 2.2, Figure 3
$t\overline{t}+E_T^{\mathrm{miss}}$	EFT Higgs portals	Section 2.2, Figure 3
	2MDM	Section 4.2
$tW+E_T^{\mathrm{miss}}$	EFT Higgs portals	Section 2.2, Figure 3
	2HDM+a	Section 3.1, Figure 7 and Figure 8
$h+E_T^{\mathrm{miss}}$	2HDM+a	Section 3.1, Figure 7, Figure 8, Figure 9 and Figure 10
	2HDM+s	Section 3.2, Figure 12
	2HDM+$Z^{\prime }$	Section 4.1, Figure 14
$Z+E_T^{\mathrm{miss}}$	2HDM+a	Section 3.1, Figure 7, Figure 8, Figure 9 and Figure 10
	2HDM+s	Section 3.2, Figure 12
$s+E_T^{\mathrm{miss}}$	2MDM	Section 4.2, Figure 18
$j+E_T^{\mathrm{miss}}$	2MDM	Section 4.2, Figure 18
$H^{\pm }\to tb$	2HDM+a	Section 3.1, Figure 7 and Figure 8
$h\to 4f$	2HDM+a	Section 3.1
Dijets	2HDM+$Z^{\prime }$	Section 4.1, Figure 14 and Figure 15
	2MDM	Section 4.2, Figure 18
$Z^{\prime }\to t\overline{t}$	2HDM+$Z^{\prime }$	Section 4.1, Figure 14 and Figure 15
	2MDM	Section 4.2, Figure 18
$Z^{\prime }\to Zh$	2HDM+$Z^{\prime }$	Section 4.1, Figure 14 and Figure 15
$h\to 4f$	Neutral naturalness	Section 5.1, Figure 20
	Dark photons	Section 5.2, Figure 22 and Figure 23
	Vector plus fermion portal	Section 5.3, Figure 26
$h\to \mathrm{inv},\mathrm{undet}$	Neutral naturalness	Section 5.1, Figure 20
	Dark photons	Section 5.2, Figure 22
Dileptons	Dark photons	Section 5.2, Figure 23
	Vector plus fermion portal	Section 5.3, Figure 26

provides a list of the LHC signatures that are discussed in this review, indicating which model is constrained by a given search and the place(s) where the corresponding discussion can be found. Our whole review is tied together in Section 6, where we present an outlook. We commence without further ado.

2. Higgs Portal Models

One of the special features of the SM Higgs doublet H is that $H^{†}H$ is the only Lorentz and gauge invariant operator with a mass dimension of two. The operator $H^{†}H$ therefore furnishes a portal to the dark or hidden sector [11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29]. In particular, at the level of dimension-four operators, one can write down the couplings of $H^{†}H$ to dark spin-0 and spin-1 fields, while the leading interactions with dark spin-$1/2$ fields are of dimension five. If the resulting EFT is equipped with a suitable symmetry, the dark field becomes stable, giving rise to a scalar, vector and fermionic DM candidate, respectively. Such a symmetry can for instance be a ${\mathbb{Z}}_2$ exchange symmetry.

2.1. Theory

In this section, we consider the simplest possibility of these Higgs portal models, namely the case of a real scalar $\varphi$ that is a singlet under the SM gauge group, but odd under a ${\mathbb{Z}}_2$ symmetry, i.e., $\varphi \to -\varphi$. This guarantees the stability of $\varphi$, making it a suitable DM candidate. The interactions between the dark sector and the SM that we consider are

$${\mathcal{L}}_{\varphi H}=c_m{\varphi }^2\left(H^{†}H\right)+\frac{c_d}{{\Lambda }^2}\left({\partial }_{\mu }{\varphi }^2\right)\left({\partial }^{\mu }\left(H^{†}H\right)\right),$$

(1)

where the first (second) term is the so-called marginal (derivative) Higgs portal, the parameter $c_m$ ($c_d$) denotes the corresponding coupling or Wilson coefficient, and $\Lambda$ is a mass scale that suppresses the derivative Higgs portal that corresponds to a dimension-six operator. Such a derivative coupling with the Higgs field arises in models where DM is a pseudo-Nambu–Goldstone boson (pNGB) [24,25,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68]. In such a case, $\Lambda$ is associated with the scale of global symmetry breaking that gives rise to the appearance of the pNGB(s). Besides the two types of interactions introduced in (1), explicit spin-0 ultraviolate (UV) completions of Higgs portal models can contain additional operators (see [22,25] for a full classification of operators up to dimension six). In order to highlight the complementarity of collider and noncollider bounds on Higgs portal models in a simple fashion, we focus in what follows on the subclass of models in which the leading effects are well captured by the EFT Lagrangian ${\mathcal{L}}_{\varphi H}$. After drawing our general conclusions, we do, however, also briefly discuss the possible impact of other operators not included in (1). Recent detailed phenomenological studies of Higgs portal models with vector and fermionic DM can be found for instance in [26,28]. See also the ATLAS and CMS publications [69,70,71,72,73].

2.2. Collider Constraints

A common feature of Higgs portal models is that they predict Higgs to invisible decays if the DM candidate is kinematically accessible, i.e., $m_{\varphi }<m_h/2$ in the case of the real scalar $\varphi$ with mass $m_{\varphi }$ and $m_h\simeq 125\mathrm{GeV}$ the mass of the SM-like Higgs boson. The most important Higgs production channels for searches of Higgs to invisible decays are displayed in Figure 1. For the effective interactions (1), the relevant partial Higgs decay width reads

$$\Gamma \left(h\to \varphi \varphi \right)=\frac{v^2}{8\pi m_h}{\left(1-\frac{4m_{\varphi }^2}{m_h^2}\right)}^{1/2}{\left(c_m+\frac{m_h^2{c}_d}{{\Lambda }^2}\right)}^2.$$

(2)

where $v\simeq 246\mathrm{GeV}$ denotes the vacuum expectation value (VEV) of the $125\mathrm{GeV}$ Higgs boson. The Formula (2) can be used to translate experimental limits on the Higgs to invisible branching ratio $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ into constraints on the strength of the marginal and derivative Higgs portals. In fact, in the limit $m_{\varphi }\ll m_h$, the best existing 95% confidence level (CL) exclusion LHC bound [73] of

$$\mathrm{BR}\left(h\to \mathrm{inv}\right)<0.11,$$

(3)

leads to

$$|c_m|<5.1·{10}^{-3},\frac{\Lambda }{\sqrt{|c_d|}}>1.7\mathrm{TeV},$$

(4)

when the SM value ${\Gamma }_h^{\mathrm{SM}}\simeq 4.07\mathrm{MeV}$ [3] of the total Higgs decay width is used. Notice that the bound (3) results from a statistical combination of searches for invisible Higgs decays, where the Higgs is produced according to the SM via VBF (see the upper-left Feynman diagram in Figure 1) or in association with a pair of top quarks (see the upper-right Feynman diagram in Figure 1) in final states with zero or two leptons. At the high-luminosity upgrade of the LHC (HL-LHC), it may be possible to set a limit on the Higgs to invisible branching ratio of $\mathrm{BR}\left(h\to \mathrm{inv}\right)<2.5·{10}^{-2}$ [74]. This implies that the bounds (4) may be improved to $2.3·{10}^{-3}$ and $2.6\mathrm{TeV}$ by the end of the LHC era.

If the DM candidate is too heavy to be pair produced as a real particle in the decay of the $125\mathrm{GeV}$ Higgs boson, $E_T^{\mathrm{miss}}$ signatures still arise from off-shell Higgs production. Possible channels to search for signals of this kind are VBF Higgs production in the $jj+E_T^{\mathrm{miss}}$ channel as well as $t\overline{t}+E_T^{\mathrm{miss}}$ and $tW+E_T^{\mathrm{miss}}$ production. Relevant Feynman diagrams are shown in Figure 2. The LHC reach of these channels in the context of (1) has been studied recently in [25,68], and we summarise the main findings of these articles below.

2.3. DM Phenomenology

DM states that couple to the $125\mathrm{GeV}$ Higgs typically lead to spin-independent (SI) DM–nucleon cross–sections (${\sigma }_{\mathrm{SI}}$), which are severely constrained by the existing DM DD experiments such as XENON1T. While the marginal Higgs portal leads to an unsuppressed SI DM–nucleon cross-section, the DM–nucleon interactions that are mediated by the derivative Higgs portal are suppressed by $q^2/{\Lambda }^2\lesssim {\left(100\mathrm{MeV}\right)}^2/{\Lambda }^2$ where $q^2$ characterises the momentum transfer in DM scattering with heavy nuclei. Explicitly one finds

$${\sigma }_{\mathrm{SI}}^{\varphi N}=\frac{c_m^2{m}_N^4{f}_N^2}{\pi m_h^4{\left(m_{\varphi }+m_N\right)}^2},$$

(5)

where $m_N\simeq 939\mathrm{MeV}$ is the average of the nucleon mass and $f_N\simeq 0.31$ [75,76,77,78] parameterises the strength of the Higgs–nucleon interactions. For $m_{\varphi }=100\mathrm{GeV}$ the latest XENON1T 90% CL upper limit on the SI DM–nucleon cross-section reads ${\sigma }_{\mathrm{SI}}<9.12·{10}^{-47}{\mathrm{cm}}^2$ [79]. By means of (5), this bound can be translated into a limit on the marginal Higgs portal coupling:

$$|c_m|<5.0·{10}^{-3}.$$

(6)

By contrast, the derivative Higgs portal coupling $c_d$ remains unconstrained by DD experiments due to the aforementioned momentum suppression. In fact, it turns out that up to dimension-six, the derivative Higgs portal is the only spin-0 DM–Higgs operator that does naturally satisfy the constraints imposed by ${\sigma }_{\mathrm{SI}}$, once radiative corrections are considered [68].

In order to understand the physics of DM ID and thermal-freeze out in Higgs portal models described by (1), let us write the velocity-averaged cross-section for the annihilation of DM into a SM final state X as

$$〈\sigma \left(\varphi \varphi \to X\right)v〉\left(T\right)=a_X+T{b}_X,$$

(7)

where T denotes the DM temperature. Notice that in today’s universe $T_0\simeq 0$, while at freeze-out $T_f\simeq m_{\varphi }/25$. The p-wave coefficient $b_X$ can therefore usually be neglected in the calculation of the ID constraints. However, it can be relevant in the case of the DM relic density (${\Omega }_{\mathrm{DM}}{h}^2$), in particular, if the s-wave coefficient $a_X$ is parametrically suppressed.

For $m_b<m_{\varphi }\lesssim m_W$ with $m_b\simeq 4.2\mathrm{GeV}$ ($m_W\simeq 80.4\mathrm{GeV}$) the bottom-quark (W-boson) mass, DM annihilation into bottom–antibottom quark pairs is the dominant contribution to ${\Omega }_{\mathrm{DM}}{h}^2$. The corresponding s-wave coefficient reads

$$a_{b\overline{b}}=\frac{3m_b^2}{\pi }{\left|\frac{1}{4m_{\varphi }^2-m_h^2+i{m}_h{\Gamma }_h}\left(c_m+\frac{4m_{\varphi }^2{c}_d}{{\Lambda }^2}\right)\right|}^2,$$

(8)

with ${\Gamma }_h$ denoting the total decay width of the $125\mathrm{GeV}$ Higgs boson including contributions from $h\to \varphi \varphi$ (see (2)). In the case $m_{\varphi }\gtrsim m_W$ the $\varphi \varphi \to W^{+}{W}^{-},ZZ,hh,t\overline{t}$ channels dominate DM annihilation. These processes all receive unsuppressed s-wave contributions. For DM masses sufficiently far above v, the relevant coefficients take the following form:

$$a_X=\frac{N_X{m}_{\varphi }^2}{\pi }{\left(\frac{c_m}{4m_{\varphi }^2}+\frac{c_d}{{\Lambda }^2}\right)}^2,a_{t\overline{t}}=\frac{3m_t^2}{\pi }{\left(\frac{c_m}{4m_{\varphi }^2}+\frac{c_d}{{\Lambda }^2}\right)}^2,$$

(9)

with $X=W^{+}{W}^{-},ZZ,hh$ and $N_{W^{+}{W}^{-}}=2$, $N_{ZZ}=N_{hh}=1$. Notice that in the case of $m_{\varphi }\gg v$, DM annihilation to W and Z bosons reduces to three times the contribution from annihilation to the $125\mathrm{GeV}$ Higgs boson. This is an expected feature in the $SU{\left(2\right)}_L\times U{\left(1\right)}_Y$ symmetric limit. In addition to the DM annihilation channels discussed above, DM annihilation into monochromatic photons can also be relevant in the context of (2). The corresponding formulas can be found for instance in [68].

In terms of (8) and (9), today’s DM relic density is approximately given by

$$\frac{{\Omega }_{\mathrm{DM}}{h}^2}{0.12}=\frac{3·{10}^{-26}{\mathrm{cm}}^3/\mathrm{s}}{{\langle \sigma v\rangle }_f},{\langle \sigma v\rangle }_f=\sum _X〈\sigma \left(\varphi \varphi \to X\right)v〉\left(T_f\right),$$

(10)

where the sum over X involves all annihilation channels that are kinematically open at a given value of $m_{\varphi }$. While the above formulas represent useful expressions to estimate ${\Omega }_{\mathrm{DM}}{h}^2$, we use micrOMEGAs [80] in our numerical analysis to obtain the constraints on the parameter space of (1) that follow from the PLANCK measurement ${\Omega }_{\mathrm{DM}}{h}^2=0.120\pm 0.001$ [81]. The ID exclusions shown below in Figure 3 are also determined with the help of micrOMEGAs.

2.4. Summary Plots

The upper (lower) panel in Figure 3 summarises the most important constraints on the marginal (derivative) Higgs portal introduced in (1). The solid black contours correspond to the current best limit on $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ as given in (3) while the dashed black lines represent the expected HL-LHC 95% CL limit $\mathrm{BR}\left(h\to \mathrm{inv}\right)<2.5·{10}^{-2}$ [74]. The purple region in the upper plot is disfavoured by the 90% CL bounds of XENON1T [79] on ${\sigma }_{\mathrm{SI}}$. The vertical orange shaded bands indicate the DM mass ranges that are excluded at 95% CL by the $\gamma$-ray observations of dwarf spheroidal galaxies (dSphs) of the Fermi-LAT and DES collaborations reported in [82]. The used experimental bounds assume DM annihilation via $\varphi \varphi \to b\overline{b}$ and that ${\Omega }_{\mathrm{DM}}{h}^2=0.12$. Compared to $\varphi \varphi \to b\overline{b}$, the constraints that follow from the latest Fermi-LAT search for monochromatic photons [83] lead to weaker constraints. These limits are hence not shown in the figure. In the parameter space below (above) the red curve, the marginal (derivative) Higgs portal model predicts ${\Omega }_{\mathrm{DM}}{h}^2>0.12$, i.e., larger values of the DM relic density compared to the PLANCK measurement [81]. The green regions correspond to the 95% CL exclusion limits found in [25] from a study of off-shell invisible Higgs production in the $\mathrm{VBF}+E_T^{\mathrm{miss}}$ channel. Finally, the blue domains represent the 95% CL constraints obtained by the combined $t\overline{t}+E_T^{\mathrm{miss}}$ and $tW+E_T^{\mathrm{miss}}$ ($tX+E_T^{\mathrm{miss}}$) analysis strategy discussed in [68]. The latter two types of collider limits assume an integrated luminosity of $3{\mathrm{ab}}^{-1}$, collected at the HL-LHC.

From the upper panel in Figure 3, it is evident that the constraints on the Wilson coefficient $c_m$ of the marginal Higgs portal from searches for Higgs to invisible decays at the LHC are more stringent than the DD bounds for DM masses $m_{\varphi }\lesssim 5\mathrm{GeV}$, while in the range $5\mathrm{GeV}\lesssim m_{\varphi }<m_h/2$, they are roughly comparable in strength. In the case $m_{\varphi }>m_h/2$, the bounds that follow from ${\sigma }_{\mathrm{SI}}$ are however by more than two orders of magnitude stronger than those that mono-X searches at the HL-LHC are expected to set. Off-shell invisible Higgs production in the VBF channel [25] is likely the best probe of the marginal Higgs portal at the LHC if $m_{\varphi }>m_h/2$. Notice however that the study [25] assumes a systematic uncertainty of $1\%$, while the shown $tX+E_T^{\mathrm{miss}}$ exclusion limits are based on a systematic uncertainty of $15\%$ [68]. Assuming a reduction of background uncertainties in $tX+E_T^{\mathrm{miss}}$ down to $5\%$ would bring the $\mathrm{VBF}+E_T^{\mathrm{miss}}$ and $tX+E_T^{\mathrm{miss}}$ exclusion limits closer together. Combining the two mono-X channels as in the case of the LHC searches for the invisible Higgs boson decays (cf. for instance [69,70,71,72,73]) can be expected to improve the HL-LHC reach. The potential of the high-energy option of the LHC, the future circular hadron-hadron collider, the compact linear collider and a muon collider in constraining the marginal Higgs portal through $\mathrm{VBF}+E_T^{\mathrm{miss}}$ off-shell Higgs production has been studied recently in [25]. For earlier analyses, see also [84,85,86,87,88].

For what concerns the derivative Higgs portal model, the lower panel in Figure 3 shows that in the Higgs on-shell region, i.e., for $m_{\varphi }<m_h/2$, HL-LHC measurements of invisible Higgs decays exclude large parts of the parameter space that lead to ${\Omega }_{\mathrm{DM}}{h}^2=0.12$. Only a narrow corridor around the $125\mathrm{GeV}$ Higgs-boson resonance survives the DM relic density constraint, which is however excluded by DM ID measurements. Given that ${\sigma }_{\mathrm{SI}}$ is momentum suppressed, the stringent limits from DM DD experiments do not put constraints on the derivative Higgs portal model. This highlights the need to test such models with $m_{\varphi }>m_h/2$ using mono-X searches at the HL-LHC; however, only if ${\Omega }_{\mathrm{DM}}{h}^2<0.12$. The best tests of the derivative Higgs portal model in the Higgs off-shell region seem again to be searches for invisible decays of the $125\mathrm{GeV}$ Higgs boson produced in the VBF channel. This conclusion however depends once more on the actual size of systematic uncertainties of the relevant mono-X channels under HL-LHC conditions.

Our discussion so far was phrased within an EFT, but concrete examples of UV complete models where the two Higgs portal interactions (1) dominate in the low-energy limit have been constructed. For instance, the pNGB DM models in [52,54,60] lead to a sizeable marginal Higgs portal coupling $c_m$, while the constructions in [24,66] manage to suppress this coupling, making the derivative Higgs portal coupling $c_d$ the leading interaction between the dark and the visible sector. As shown above, in the latter category of models, only DM production at the LHC is able to directly probe pNGB DM models. If the DM candidate can be produced as a real particle, the searches for invisible Higgs boson decays play a key role in such explorations, while DM masses above the Higgs threshold can be tested by studying mono-X signatures such as $\mathrm{VBF}+E_T^{\mathrm{miss}}$ or $tX+E_T^{\mathrm{miss}}$. Dedicated experimental searches and/or interpretations by ATLAS and CMS of the relevant mono-X signatures in the pNGB DM context do not exist at the present.

2.5. Further Considerations

As promised, we now return to a brief discussion of Higgs portal models with interactions not encoded in (1) such as models with fermionic DM. In Figure 4, we show a comparison between the upper limits at 90% CL from DD experiments [79,89,90,91] on the SI DM–nucleon cross-section and the exclusion limits that are derived from the latest ATLAS measurement of Higgs to invisible decays [73]. At 90% CL, this measurement leads to $\mathrm{BR}\left(h\to \mathrm{inv}\right)<0.09$—the corresponding 95% CL bound is given in (3). The translation of the $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ bound into a limit on ${\sigma }_{\mathrm{SI}}$ relies on an EFT approach under the assumption that the $125\mathrm{GeV}$ Higgs boson decays to a pair of DM particles are kinematically possible and that the DM particle is either a scalar or a Majorana fermion. In the scalar case, the EFT approach boils down to extract limits on the Wilson coefficient $c_m$ from (2) and to insert the obtained values into (5) to derive bounds on ${\sigma }_{\mathrm{SI}}$. Notice that the limit on the Wilson coefficient $c_m$ does only marginal change when going from DM masses of $10\mathrm{GeV}$ down to $1\mathrm{GeV}$ (see the upper panel in Figure 3) while the bound on the SI DM–nucleon cross-section worsens notable. While this is puzzling at first, one has to realise that in the scalar case ${\sigma }_{\mathrm{SI}}$ scales with $1/{(m_{\varphi }+m_N)}^2$(cf. (5)) which implies that for constant $c_m$ and the DM mass $m_{\varphi }$ sufficiently larger than $m_N$, one has ${\sigma }_{\mathrm{SI}}\propto 1/m_{\varphi }^2$. This feature leads to a deterioration of the limit on ${\sigma }_{\mathrm{SI}}$ with decreasing DM masses, although the bound on $c_m$ in fact even slightly improves.

In order to understand the behaviour of the exclusion limit for the Higgs portal model with Majorana DM $\chi$ as shown in Figure 4, we first provide expressions for the relevant effective interactions, the partial Higgs decay width for $h\to \chi \overline{\chi }$ and the SI DM–nucleon cross-section. The Higgs portal takes the following form:

$${\mathcal{L}}_{\chi H}=\frac{c_{\chi }}{\Lambda }\chi \overline{\chi }\left(H^{†}H\right),$$

(11)

where $c_{\chi }$ is a Wilson coefficient assumed to be real and the scale $\Lambda$ suppresses the interaction because the operator is of dimension five. The corresponding partial Higgs decay width reads

$$\Gamma \left(h\to \chi \overline{\chi }\right)={\left(1-\frac{4m_{\chi }^2}{m_h^2}\right)}^{3/2}\frac{c_{\chi }^2{v}^2{m}_h}{4\pi {\Lambda }^2},$$

(12)

while the SI DM–nucleon cross-section is given by

$${\sigma }_{\mathrm{SI}}^{\chi N}=\frac{4c_{\chi }^2{m}_{\chi }^2{m}_N^4{f}_N^2}{\pi m_h^4{\left(m_{\chi }+m_N\right)}^2{\Lambda }^2}.$$

(13)

We add that for a complex $c_{\chi }$, the SI DM–nucleon cross-section can be strongly suppressed (cf. for instance [92]). From (12), it follows that for $m_{\chi }\to 0$, measurements of $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ simply lead to constant bound on the combination $|c_{\chi }|/\Lambda$ of parameters. The SI DM–nucleon cross-section (13) however scales as $m_{\chi }^2/{(m_{\chi }+m_N)}^2$ where compared to (5), the additional factor of $m_{\chi }^2$ appears due to dimensional reasons. It follows that decreasing the DM mass leads to a steady improvement of the limit on ${\sigma }_{\mathrm{SI}}$ from $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ (see Figure 4) even though $|c_{\chi }|/\Lambda$ stays essentially constant.

The above discussion suggests that from the point of view of collider physics, interpreting the limits on $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ in terms of exclusions on ${\sigma }_{\mathrm{SI}}$ is not the optimal choice. A better way to interpret and to compare the $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ results to the bounds on the SI DM–nucleon cross-section obtained by DD experiments would in our opinion consist in showing limits on the effective interaction strength of the Higgs portal models (i.e., $|c_m|$, $\Lambda /|c_d|$, $\Lambda /|c_{\chi }|$, etc.). In the case of the marginal and derivative Higgs portal, this has been performed in Figure 3, and this is also the standard presentation in most of the theoretical literature [11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29] on Higgs portal models.

3. Portals with Extended Higgs Sectors

In the spin-0 simplified models with a fermionic singlet DM candidate $\chi$, which were recommended in the articles [36,37] as benchmarks for DM searches at the LHC, the interactions between the mediator and the SM fermions are not invariant under the $SU{\left(2\right)}_L\times U{\left(1\right)}_Y$ gauge symmetry. This feature leads to unitarity violation at high energies [93,94]—for discussions of unitarity violation in the case of spin-1 simplified DM models, see [95,96,97]. The simplest way to restore gauge invariance in spin-0 simplified DM models is to introduce a singlet scalar that provides the portal to the dark sector and mixes with the $125\mathrm{GeV}$ Higgs [98,99,100,101,102,103,104]. Compared to the spin-0 simplified DM models, this model leads to additional $E_T^{\mathrm{miss}}$ signatures such as Higgs to invisible decays, mono-W, mono-Z and VBF + $E_T^{\mathrm{miss}}$ production at tree level—see [105] for a detailed discussion of the collider phenomenology—but the stringent constraints on the scalar-Higgs mixing that arises from the Higgs measurements at the LHC (cf. [106,107] for the latest global Higgs analyses by ATLAS and CMS) lead to a suppression of all $E_T^{\mathrm{miss}}$ signals. This suppression still allows the LHC to test the model via (3) if $m_{\chi }<m_h/2$. However, in the case $m_{\chi }>m_h/2$, the LHC coverage of the parameter space of the singlet fermionic DM model is very limited, in particular, if one compares the collider limits to the strong bounds that result from DM DD experiments. We therefore do not discuss the singlet fermionic DM model any further.

Restoring gauge invariance in spin-0 simplified DM models and simultaneously satisfying the stringent LHC constraints from the Higgs boson measurements is possible if the SM Higgs sector is extended. Two-Higgs-doublet models (2HDMs) [108,109] that contain two Higgs doublets $H_1$ and $H_2$ instead of just H are one of the natural choices for such extensions. The tree-level 2HDM scalar potential can be written as

$$\begin{array}{cc}\hfill V_H& ={\mu }_1{H}_1^{†}{H}_1+{\mu }_2{H}_2^{†}{H}_2+\left({\mu }_3{H}_1^{†}{H}_2+\mathrm{h}.\mathrm{c}.\right)+{\mathsf{\lambda }}_1{\left(H_1^{†}{H}_1\right)}^2+{\mathsf{\lambda }}_2{\left(H_2^{†}{H}_2\right)}^2\hfill \\ \hfill & \phantom{xx}+{\mathsf{\lambda }}_3\left(H_1^{†}{H}_1\right)\left(H_2^{†}{H}_2\right)+{\mathsf{\lambda }}_4\left(H_1^{†}{H}_2\right)\left(H_2^{†}{H}_1\right)+\left[{\mathsf{\lambda }}_5{\left(H_1^{†}{H}_2\right)}^2+\mathrm{h}.\mathrm{c}.\right],\hfill \end{array}$$

(14)

where we have imposed a ${\mathbb{Z}}_2$ symmetry under which $H_1\to H_1$ and $H_2\to -H_2$ but allowed for this discrete symmetry to be broken softly by ${\mu }_3{H}_1^{†}{H}_2+\mathrm{h}.\mathrm{c}.$ The ${\mathbb{Z}}_2$ symmetry is the minimal condition necessary to guarantee the absence of flavour-changing neutral currents (FCNCs) at tree level [110,111], and such a symmetry is realised in many well-motivated complete UV theories in the form of supersymmetry (SUSY), a $U\left(1\right)$ symmetry or a discrete symmetry acting on the Higgs doublets. In order to avoid possible issues with electric dipole moments, it is commonly assumed that all parameters in the scalar potential (14) are real. In such a case, the CP eigenstates that arise after spontaneous symmetry breaking from $V_H$ can be identified with the mass eigenstates, giving rise to two CP-even scalars h and H, one CP-odd pseudoscalar A and one charged scalar $H^{\pm }$. Besides the masses of the five Higgs bosons, the 2HDM parameter space involves the angles $\alpha$ and $\beta$. The former angle describes the mixing of the two CP-even Higgs bosons while the latter encodes the ratio of the VEVs $v_1$ and $v_2$ of the two Higgs doublets $tan\beta =v_2/v_1$ with $v=\sqrt{v_1^2+v_2^2}$.

In the context of the inert doublet model [112,113,114], the scalar potential (14) alone already allows for interesting DM and collider phenomenology [115,116,117,118,119,120,121,122,123,124,125,126,127,128,129]. The DM particle in the inert doublet model is the lightest neutral component of the second Higgs doublet, making it a spin-0 state. To accommodate the possibility of having a spin-1/2 singlet DM candidate (such as in the case of the simplified DM models [36,37]) considering only (14) is therefore not enough and one generically needs an additional mediator, besides the DM particle. If one insists that this mediator has spin-0 and one does not want to violate CP, the additional mediator can either be a scalar or pseudoscalar. Dedicated studies of the collider aspects and DM properties in the latter type of next-generation simplified DM models have been presented in [30,31,32,33,34,35,130,131,132,133,134,135,136,137,138,139,140,141,142]. Below, we review in detail two of the models discussed in these articles.

3.1. 2HDM+a Model

The 2HDM+a model [30,31,32,33] is the simplest gauge-invariant and renormalisable extension of the simplified pseudoscalar DM model [36,37]. It includes a Dirac fermion $\chi$, which transforms as a singlet under the SM gauge group and therefore provides a DM candidate, and a CP-odd mediator P that furnishes the dominant portal between the dark and the visible sector. Since the DM DD constraints are weaker for models with pseudoscalar mediators compared to models with scalar mediators, the observed DM relic abundance can be reproduced in large regions of parameter space. These features allow for a host of $E_T^{\mathrm{miss}}$ signatures at colliders which can be consistently compared and combined, making the 2HDM+a model one of the main pillars of the LHC DM search programme [136,143,144,145,146,147,148,149,150,151,152].

3.1.1. Theory

If the Dirac DM field $\chi$ and the pseudosalar P are taken to transform under the ${\mathbb{Z}}_2$ symmetry as $\chi \to -\chi$ and $P\to P$, the only renormalisable DM–mediator coupling that is allowed by symmetry is

$${\mathcal{L}}_{\chi }=-i{y}_{\chi }P\overline{\chi }{\gamma }_5\chi .$$

(15)

Here, it is assumed that the dark-sector Yukawa coupling $y_{\chi }$ is real in order not to violate CP. Besides (14), the scalar potential in the 2HDM+a model contains the following terms

$$V_{HP}=P\left(i{b}_P{H}_1^{†}{H}_2+\mathrm{h}.\mathrm{c}.\right)+P^2\left({\mathsf{\lambda }}_{P1}{H}_1^{†}{H}_1+{\mathsf{\lambda }}_{P2}{H}_2^{†}{H}_2\right),$$

(16)

that connect singlets to doublets. The parameters $b_P$, ${\mathsf{\lambda }}_{P1}$ and ${\mathsf{\lambda }}_{P2}$ are taken to be real to not violate CP. Notice that the first term in $V_{HP}$ breaks the ${\mathbb{Z}}_2$ symmetry softly. The singlet potential reads

$$V_P=\frac{1}{2}{m}_P^2{P}^2.$$

(17)

A quartic term $P^4$ is not included in (17) since it does not lead to any substantial modification of the LHC phenomenology—in particular, such an addition would have no relevant effects in any of the $E_T^{\mathrm{miss}}$ observables discussed below.

Including the mass of the DM particle, the Lagrangian of the 2HDM+a model contains 14 free parameters in addition to the SM ones. After rotation to the mass eigenbasis, these 14 parameters can be traded for seven physical masses, three mixing angles and four couplings:

$$\left\{\begin{array}{c}{\mu }_1,{\mu }_2,{\mu }_3,b_P,m_P,m_{\chi },\\ y_{\chi },{\mathsf{\lambda }}_1,{\mathsf{\lambda }}_2,{\mathsf{\lambda }}_3,{\mathsf{\lambda }}_4,{\mathsf{\lambda }}_5,\\ {\mathsf{\lambda }}_{P1},{\mathsf{\lambda }}_{P2}\end{array}\right\}⟺\left\{\begin{array}{c}v,m_h,m_A,m_H,m_{H^{\pm }},m_a,m_{\chi },\\ cos(\beta -\alpha ),tan\beta ,sin\theta ,\\ y_{\chi },{\mathsf{\lambda }}_3,{\mathsf{\lambda }}_{P1},{\mathsf{\lambda }}_{P2}\end{array}\right\}.$$

(18)

Here, $sin\theta$ represents the mixing of the two CP-odd weak spin-0 eigenstates, and the additional CP-odd mediator a is mostly composed of P for $sin\theta \simeq 0$. The parameters appearing on the right-hand side of (18) are used as input in the analyses of the 2HDM+a model. Since the VEV v and the Higgs mass $m_h$ are already fixed by observations, there are 12 input parameters.

Some of the 2HDM+a parameters are constrained by Higgs physics, electroweak (EW) precision observables (EWPOs), vacuum stability considerations, flavour physics and LHC searches for additional spin-0 bosons. The mixing angle $\alpha$ between the CP-even scalars h and H is for instance constrained by Higgs coupling strength measurements [106,107]. For arbitrary values of $tan\beta$, only parameter choices with $cos(\beta -\alpha )\simeq 0$ are experimentally allowed. In order to satisfy the constraints from Higgs physics, the existing experimental 2HDM+a analyses [143,144,145,146,147,148,149,150,151,152] have concentrated on the so-called alignment limit of the 2HDM where $cos(\beta -\alpha )=0$ [153], treating $tan\beta$ as a free parameter.

The measurements of the EWPOs constrain the differences between the masses of the additional scalar and pseudoscalar particles $m_H$, $m_A$, $m_{H^{\pm }}$ and $m_a$, because the exchange of spin-0 states modifies the propagators of the EW bosons starting at the one-loop level. In [33], it was shown that the sum of the potentials (14) and (16) has a custodial symmetry if $cos(\beta -\alpha )=0$ and $m_A=m_H=m_{H^{\pm }}$. For such parameter choices, the EWPOs are satisfied for any value of $sin\theta$ and $m_a$, which renders the choice $m_A=m_H=m_{H^{\pm }}$ a good starting point to explore the 2HDM+a parameter space.

The requirement that the scalar potential $V_H+V_{HP}+V_P$ of the 2HDM+a model is bounded from below restricts the possible choices of the spin-0 boson masses, mixing angles and quartic couplings. Assuming that ${\mathsf{\lambda }}_{P1},{\mathsf{\lambda }}_{P2}>0$ and $m_A=m_H=m_{H^{\pm }}$, one can show [136] that there are two bounded from below conditions:

$${\mathsf{\lambda }}_3>\frac{m_h^2}{v^2}\simeq 0.26,{\mathsf{\lambda }}_3>\frac{m_A^2-m_a^2}{v^2}{sin}^2\theta -\frac{m_h^2}{v^2}{cot}^2\left(2\beta \right).$$

(19)

These inequalities suggest that ${\mathsf{\lambda }}_3$ has to be sufficiently large, in particular, if the mass splitting of the two pseudoscalars and/or $sin\theta$ are large. Since the relations (19) are modified in models with more general scalar potentials including dimension-six operators [133,154], the bounded-from-below requirements are not directly imposed in the ATLAS and CMS interpretations of the $E_T^{\mathrm{miss}}$ searches [143,144,145,146,147,148,149,150,151,152]. Instead, the choice ${\mathsf{\lambda }}_3=3$ is employed, which generically allows for heavy Higgses above $1\mathrm{TeV}$, while keeping ${\mathsf{\lambda }}_3$ small enough to stay in the perturbative regime.

The quartic couplings ${\mathsf{\lambda }}_3$, ${\mathsf{\lambda }}_{P1}$ and ${\mathsf{\lambda }}_{P2}$ affect the decay pattern of the heavy CP-odd and CP-even 2HDM scalars. In the alignment limit and assuming that $m_A=m_H=m_{H^{\pm }}$, the $Aah$ and $Haa$ couplings take the following form [33]

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill g_{Aah}& =\frac{1}{2m_{A}v}\left[m_A^2-m_a^2+m_h^2-2{\mathsf{\lambda }}_3{v}^2+2\left({\mathsf{\lambda }}_{P1}{cos}^2\beta +{\mathsf{\lambda }}_{P2}{sin}^2\beta \right)v^2\right]sin\left(2\theta \right),\hfill \\ \hfill g_{Haa}& =\frac{1}{m_{H}v}\left[2cot\left(2\beta \right)\left(m_h^2-{\mathsf{\lambda }}_3{v}^2\right){sin}^2\theta +sin\left(2\beta \right)\left({\mathsf{\lambda }}_{P1}-{\mathsf{\lambda }}_{P2}\right)v^2{cos}^2\theta \right].\hfill \end{array}\end{array}$$

(20)

These expressions imply that parameter choices of the form ${\mathsf{\lambda }}_3={\mathsf{\lambda }}_{P1}={\mathsf{\lambda }}_{P2}$ are well suited to keep the total widths ${\Gamma }_A$ and ${\Gamma }_H$ better under control in the limit of heavy Higgs masses, since such choices lead to cancellations in (20).

Indirect constraints on the charged Higgs-boson mass $m_{H^{\pm }}$ arise from $Z\to b\overline{b}$, $B\to X_s\gamma$, $B_s\to {\mu }^{+}{\mu }^{-}$ and B-meson mixing—see [33,136,142] for details and relevant references. These constraints are more model dependent than the bounds that derive from the EWPOs because they depend on the choice of the Yukawa sector of the 2HDM model. For instance, in the case of the 2HDM of type-II, the inclusive $B\to X_s\gamma$ decay sets a 95% CL lower limit $m_{H^{\pm }}>800\mathrm{GeV}$ [155,156] that for $tan\beta \gtrsim 2$ is practically independent of the specific choice of $tan\beta$. In other 2HDM realisations such as a fermiophobic 2HDM model of type-I (see [157,158] for recent detailed discussions), all flavour bounds can however be significantly relaxed, allowing for EW-scale values of $m_{H^{\pm }}$. Furthermore, since the FCNC constraints arise from loop corrections, they can in principle be weakened by the presence of additional particles that are too heavy to be produced at the LHC. This makes the bounds from flavour indicative, and the analyses [143,144,145,146,147,148,149,150,151,152] have not directly imposed them on the parameter space of the 2HDM+a model.

Direct searches for heavy Higgs bosons can also be used to explore and to constrain the parameter space of the 2HDM+a model. Discussions of the existing LHC constraints can be found in [33,136,141]. The most interesting signatures arise if one deviates from the alignment limit and/or gives up the assumption $m_A=m_H=m_{H^{\pm }}$. In particular, searches for final states involving tops and/or W bosons [132,137,141,145,154] provide interesting avenues to test the 2HDM+a model, and future LHC searches should consider such channels and interpret the obtained results in the context of the 2HDM+a model. In Section 3.1.3, we give an example of a non-$E_T^{\mathrm{miss}}$ that already now provides relevant constraints on the 2HDM+a parameter space.

3.1.2. Anatomy of $E_t^{\mathrm{miss}}$ Signatures

The 2HDM+a model provides a wide variety of $E_T^{\mathrm{miss}}$ signatures that can be searched for at the LHC. As pointed out in the works [32,33,132], the most interesting signals are $h+E_T^{\mathrm{miss}}$ (i.e., mono-Higgs), $Z+E_T^{\mathrm{miss}}$ (i.e., mono-Z) and $tW+E_T^{\mathrm{miss}}$ production, because these signals can be resonantly enhanced. Examples of representative diagrams that lead to the discussed $E_T^{\mathrm{miss}}$ signals in ggF production are displayed in Figure 5. Diagrams with bottom-quark loops and graphs in which the internal a is replaced by an A also exist. In the case of the mono-Higgs signature, it is evident from the figure that for $m_A>m_h+m_a$ the triangle graph shown on the left-hand side in the upper row allows for resonant $h+E_T^{\mathrm{miss}}$ production. Similar resonance enhancements arise from the graphs on the left-hand side for the mono-Z (middle row) and $tW+E_T^{\mathrm{miss}}$ channel (lower row) if $m_H>m_Z+m_a$ and $m_{H^{\pm }}>m_W+m_a$, respectively. Notice that resonant mono-Higgs, mono-Z and $tW+E_T^{\mathrm{miss}}$ production is not possible in the simplified pseudoscalar DM model [36,37] because the mediator couples only to fermions at tree level. As a result, only the Feynman diagrams displayed on the right-hand side of Figure 5 are present in this model. Notice that the appearance of the resonance contributions due to an on-shell A, H and $H^{\pm }$ not only enhances the mono-Higgs, mono-Z and $tW+E_T^{\mathrm{miss}}$ compared to the simplified pseudoscalar DM model but also leads to a quite different kinematics in these channels [32,33,132,136]. Besides ggF production, important contributions to the mono-Higgs and mono-Z channels can arise from $b\overline{b}$-fusion production in 2HDM realisations such as type-II models in which the bottom-quark Yukawa coupling is enhanced by $tan\beta$. The corresponding Feynman diagrams are shown in Figure 6. As for the ggF mono-Higgs and mono-Z signals, there are also resonant (left column) and nonresonant contributions (right column) in the case of $b\overline{b}$-fusion.

An additional $E_T^{\mathrm{miss}}$ signature that allows one to put constraints on the parameter space of the 2HDM+a model are invisible Higgs decays [33]. The relevant decay channel is $h\to aa$ followed by $a\to \chi \overline{\chi }$. The decay of the $125\mathrm{GeV}$ Higgs to a pair of pseudoscalars a is proportional to the square of the $haa$ coupling. For $cos(\beta -\alpha )=0$ and $m_A=m_H=m_{H^{\pm }}$, this coupling takes the form

$$\begin{array}{cc}\hfill g_{haa}& =\frac{1}{m_{h}v}[\left(2m_A^2-2m_a^2+m_h^2-2{\mathsf{\lambda }}_3{v}^2\right){sin}^2\theta \hfill \\ \hfill & -2\left({\mathsf{\lambda }}_{P1}{cos}^2\beta +{\mathsf{\lambda }}_{P2}{sin}^2\beta \right)v^2{cos}^2\theta ].\hfill \end{array}$$

(21)

Ignoring fine-tuned solutions for which $|g_{haa}|\ll 1$, one can show that for sufficiently light DM masses, the bound (3) leads to a lower limit of $m_a\gtrsim 100\mathrm{GeV}$. We add that the direct measurements of the total Higgs width [159,160] that impose ${\Gamma }_h<1.1\mathrm{GeV}$ at 95% CL furthermore imply that $m_a\gtrsim m_h/2\simeq 62.5\mathrm{GeV}$. This bound also holds for heavy DM, because the processes $h\to aa$ with $a\to f\overline{f}$, where f denotes all the kinematically accessible SM fermions, always provides a sizeable non-SM contribution to ${\Gamma }_h$ unless (21) is artificially small.

In addition to the $E_T^{\mathrm{miss}}$ signatures discussed so far, the 2HDM+a model also leads to mono-jet, $t\overline{t}+E_T^{\mathrm{miss}}$ and $b\overline{b}+E_T^{\mathrm{miss}}$ signals. In contrast to the case of the simplified pseudoscalar DM model where these mono-X channels provide the leading constraints (see [94,143,161,162,163,164,165,166,167,168] for the latest theoretical and experimental analyses), this is not the case in the 2HDM+a, because the mono-jet, $t\overline{t}+E_T^{\mathrm{miss}}$ and $b\overline{b}+E_T^{\mathrm{miss}}$ channels do not receive resonantly enhanced contributions. Reinterpreting the existing spin-0 simplified DM model searches in the 2HDM+a context is relatively straightforward [33,136].

3.1.3. Summary Plots

In the articles [136,143,144,145,146,147,148,149,150,151,152], various 2HMD+a benchmark scans have been recommended and studied. The following parameters are common to all existing benchmark scans,

$$m_A=m_H=m_{H^{\pm }},cos\left(\beta -\alpha \right)=0,y_{\chi }=1,{\mathsf{\lambda }}_3={\mathsf{\lambda }}_{P1}={\mathsf{\lambda }}_{P2}=3,$$

(22)

and we have describedthese choices in Section 3.1.1. Unless the relevant parameter is scanned over the plots shown below also employing

$$tan\beta =1,m_{\chi }=10\mathrm{GeV},$$

(23)

where the latter choice ensures a sizeable $\mathrm{BR}(a\to \chi \overline{\chi })$ for the pseudoscalar masses of interest, i.e., $m_a\gtrsim 100\mathrm{GeV}$. The Yukawa sector of the 2HDM is chosen to be of type-II. In the alignment limit, the couplings of the neutral non-SM scalars to the top-quark and bottom-quark therefore scale like

$$\begin{array}{cc}\hfill & g_{Ht\overline{t}}\propto y_{t}cot\beta ,g_{Hb\overline{b}}\propto y_{b}tan\beta ,\hfill \\ \hfill & g_{At\overline{t}}\propto y_{t}cot\beta cos\theta ,g_{Ab\overline{b}}\propto y_{b}tan\beta cos\theta ,\hfill \\ \hfill & g_{at\overline{t}}\propto y_{t}cot\beta sin\theta ,\phantom{i}{g}_{ab\overline{b}}\propto y_{b}tan\beta sin\theta ,\hfill \end{array}$$

(24)

where $y_t=\sqrt{2}{m}_t/v\simeq 0.94$ and $y_b=\sqrt{2}{m}_b/v\simeq 0.02$ are the top-quark and bottom-quark Yukawa coupling, respectively. These expressions imply that for the benchmark choice (23) of $tan\beta$ the top-quark loop diagrams in Figure 5 provide the dominant contributions to the mono-Higgs and mono-Z, while for sufficiently large values of $tan\beta$, the $b\overline{b}$-induced graphs of Figure 6 can give a relevant or even the dominant contribution.

The mono-X searches that provide at present the bulk of the sensitivity to the 2HDM+a model are the $h+E_T^{\mathrm{miss}}$ channel with either $h\to b\overline{b}$($h\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$) [144,150] or $h\to \gamma \gamma$ ($h\left(\gamma \gamma \right)+E_T^{\mathrm{miss}}$) [149] decays, the $Z+E_T^{\mathrm{miss}}$ signal followed by $Z\to {\ell }^{+}{\ell }^{-}$($Z\left({\ell }^{+}{\ell }^{-}\right)+E_T^{\mathrm{miss}}$) [147,151] and $tW+E_T^{\mathrm{miss}}$ production [148]. Thus, we focus on these searches below, disregarding other mono-X searches such as mono-jets or $t\overline{t}+E_T^{\mathrm{miss}}$. To evade the limits from invisible Higgs decays (cf. Section 3.1.2), we consider only $m_a$ values satisfying $m_a>100\mathrm{GeV}$ when studying the 2HDM+a parameter space. In addition, we show that the recent search for charged Higgs bosons in the channel $H^{\pm }\to tb$ [169] allows one to set stringent constraints on $m_{H^{\pm }}$, which under the assumption of degenerated 2HDM spin-0 masses (22), translate into bounds on $m_A$ or $m_H$. As we can see, the obtained limits are to first approximation independent of $m_a$, and as a result, the $H^{\pm }\to tb$ search provides complementary constraints with respect to the mono-X signatures, because these searches depend on the precise choice for $m_a$. Other non-$E_T^{\mathrm{miss}}$ searches such as for instance four-top production [136,145] are not considered.

Scans in the $m_A$–$m_a$ Plane

An assortment of 95% CL exclusions in the $m_A$–$m_a$ plane are shown in Figure 7 for $sin\theta =0.35$ (left panel) and $sin\theta =0.7$ (right panel). At present, the $h\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$ and the $Z\left({\ell }^{+}{\ell }^{-}\right)+E_T^{\mathrm{miss}}$ signatures allow one to constrain the largest parts of the $m_A$–$m_a$ plane in both $sin\theta$ scenarios. In fact, in the high-mass region, the $h\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$ is the most sensitive strategy in the case of $sin\theta =0.35$, while it becomes less sensitive than the $Z\left({\ell }^{+}{\ell }^{-}\right)+E_T^{\mathrm{miss}}$ search close to $m_A=m_a+m_h$ (dashed grey line), where the $E_T^{\mathrm{miss}}$ spectrum becomes softer and the acceptance is reduced due to the high threshold of the $E_T^{\mathrm{miss}}$ triggers. The $Z\left({\ell }^{+}{\ell }^{-}\right)+E_T^{\mathrm{miss}}$ search on the other hand relies on lepton triggers, making it more sensitive to the parameter space close to $m_A=m_a+m_h$. One also observes that for $sin\theta =0.7$, which corresponds to close-to-maximal mixing in the pseudoscalar sector, the $Z\left({\ell }^{+}{\ell }^{-}\right)+E_T^{\mathrm{miss}}$ and $tW+E_T^{\mathrm{miss}}$ signatures have increased sensitivity compared to the case of low $sin\theta$. This feature is readily understood by noticing that the couplings relevant for resonant production scale with $sin\theta$. In the alignment limit with $m_A=m_H=m_{H^{\pm }}$, one has explicitly:

$$\begin{array}{cc}\hfill g_{HaZ}& =\frac{1}{m_{H}v}\left[{(m_A^2-m_a^2-m_Z^2)}^2-4m_a^2{m}_Z^2\right]sin\theta ,\hfill \\ \hfill g_{H^{\pm }a{W}^{\pm }}& =\frac{m_W^2}{m_{H^{\pm }}v}\left[m_W^2-4m_A^2\right]sin\theta .\hfill \end{array}$$

(25)

In the case of the mono-Higgs exclusions, two features that are visible in the plot on the left-hand side of Figure 7 deserve an explanation. To understand the dips in the $h+E_T^{\mathrm{miss}}$ constraints, one first has to realise that (20) and (21) imply that for an on-shell A the real part of the sum of the $gg\to A\to ha\to h+\chi \overline{\chi }$ and $gg\to a\to ha\to h+\chi \overline{\chi }$ amplitudes is proportional to

$$\frac{\left(2m_A^2-2m_a^2+m_h^2\right){sin}^2\theta -6v^2}{m_A^2-m_a^2}.$$

(26)

Here, we have employed the parameter choices (22). For fixed values of $m_a$ and $sin\theta$, the numerator in (26) is not sign definite being negative (positive) for sufficiently small (large) $m_A$. Close to the zero point of (26), the resonant contribution to mono-Higgs production is hence strongly suppressed, leading to a weakening of the constraints which in such a case arise almost entirely from the nonresonant contributions—see the right Feynman diagram in the upper row of Figure 5.

The enhanced sensitivity of the mono-Higgs searches at high $m_A$ is due to the quadratic dependence of (21) on the 2HDM pseudoscalar mass. We also note that for large $m_A$–$m_a$ mass splittings, the widths of all non-SM Higgs bosons become large with the effect being more pronounced for increasing $sin\theta$—cf. (20), (21) and (25). While off-shell effects are taken into account in Figure 7, the possible modifications of the line shape of the intermediate Higgs bosons [170,171,172] are ignored. The latter effects have been studied in [173,174], where it has been shown that for a heavy Higgs boson different treatments of its propagator can lead to notable difference in the inclusive production cross-sections compared to the case of a Breit-Wigner with a fixed width, as used here. The hatched grey regions in the two plots of Figure 7 correspond to the parameter space where any of the additional Higgs bosons have a relative width of more than 20%. The mono-X exclusions in these regions carry some (hard to quantify) model dependence related to the precise treatment of the widths of the internal Higgs bosons.

Notice finally that in the case of the benchmark scenario (22) the constraints from $H^{\pm }\to tb$ are to a very good approximation independent of $m_a$, and we find that charged Higgs boson masses $m_{H^{\pm }}\lesssim 700\mathrm{GeV}$ ($m_{H^{\pm }}\lesssim 600\mathrm{GeV}$) are strongly disfavoured for $sin\theta =0.35$ ($sin\theta =0.7$). Searches for $H^{\pm }\to tb$ hence cover an area largely complementary to the results of the mono-X searches. We also add that the obtained limits are almost as strong as the indirect bounds that follow from $B\to X_s\gamma$ (see Section 3.1.1). With more luminosity to be collected at the LHC, one can expect the direct charged Higgs searches to become competitive with the flavour bounds.

Scans in $tan\beta$

The 95% CL limits that the searches [147,148,149,150,151,169] set in the $m_a$–$tan\beta$ and $m_A$–$tan\beta$ planes are displayed in Figure 8. For $tan\beta \gtrsim 5$, the $b\overline{b}$-induced production becomes dominant for both the mono-Higgs and mono-Z signatures, since the couplings (24) of the bottom quark to the neutral Higgs bosons are proportional to $tan\beta$ in the case of a 2HDM model of type-II. In fact, the sensitivity of the $h+E_T^{\mathrm{miss}}$ searches at high $tan\beta$ solely stems from $b\overline{b}$-fusion production, while at low $tan\beta$, the sensitivity comes entirely from ggF production. Production via $b\overline{b}$-fusion leads to final states with more bottom-quarks jets (b-jets); therefore, dedicated selections with additional b-jets can be employed to increase the sensitivity in this region. In particular, the $h\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$ exclusion at high $tan\beta$ and low $sin\theta$ has been shown to come almost exclusively from a selection with additional b-jets [150]. Notice that the $H^{\pm }\to tb$ search is sensitivity to both low and high values of $tan\beta$, since the associated coupling involves both $y_{t}cot\beta$ and $y_{b}tan\beta$ terms. Parameter regions with $tan\beta \lesssim 0.3$ are incompatible with the requirement of having a perturbative top-quark Yukawa coupling [109], and therefore, this region is not displayed in the exclusion plots, while flavour observables disfavour very high values of $tan\beta$—see for instance [33,136,142].

Scans in $sin\theta$

In Figure 9, we show the 95% CL exclusion limits as a function of $sin\theta$ that follow from [149,150,151] for two sets of pseudoscalar masses $m_A$ and $m_a$. The sensitivity of all searches improves as $sin\theta$ increases from zero, since for $sin\theta =0$, the DM mediator a decouples. The shape of the exclusion curves is due to the interplay between the production cross-section and the acceptance, in particular at intermediate values of $sin\theta$ values, the $E_T^{\mathrm{miss}}$ spectrum becomes harder, leading to an increased signal acceptance. The improvement of the sensitivity of the $h\left(\gamma \gamma \right)+E_T^{\mathrm{miss}}$ search at high $m_A$ and $sin\theta$ is due to the scaling of the $g_{haa}$ coupling (21), as discussed previously.

DM Mass Scan

Mono-X searches can also provide constraints on the mass $m_{\chi }$ of the DM candidate. Figure 10 illustrates these constraints at 95% CL for (22), (23), $m_A=600\mathrm{GeV}$, $m_a=250\mathrm{GeV}$ and $sin\theta =0.35$. In the region $m_{\chi }<m_a/2$, where the $a\to \chi \overline{\chi }$ decay is kinematically allowed, a variation of $m_{\chi }$ has a negligible impact on the production cross-sections and the $E_T^{\mathrm{miss}}$ spectrum, and in consequence the sensitivity of the $X+E_T^{\mathrm{miss}}$ searches is independent of $m_{\chi }$. For higher masses, the production cross-sections decrease significantly, and the $E_T^{\mathrm{miss}}$ becomes softer, significantly worsening the sensitivity. For this particular benchmark, the existing searches can exclude DM masses below $m_{\chi }\lesssim 140\mathrm{GeV}$. The DM relic density calculation for this benchmark is also shown in Figure 10. This calculation is performed using MadDM [175] and relies on the simplified assumption that the DM relic density is solely determined by the interactions predicted in the model. This assumption would be violated in the presence of additional hidden degrees of freedom or interactions. As a result, the overproduction or underproduction of DM should not be interpreted as an argument for excluding certain parameter ranges of the model [136,176].

Outlook

The above discussion shows that the 2HDM+a model predicts a rich phenomenology of processes resulting in a diverse range of final-state signatures with and without $E_T^{\mathrm{miss}}$. The searches that are considered in this review to constrain the 2HDM+a parameter space include the mono-Higgs, mono-Z, $tW+E_T^{\mathrm{miss}}$ and $H^{\pm }\to tb$ channels—see also the recent ATLAS note [152]. The HL-LHC prospects of the $tW+E_T^{\mathrm{miss}}$ and four-top channel have been discussed in [137,145]. All existing 2HDM+a collider studies [32,33,132,136,140,141,142,143,144,145,146,147,148,149,150,151,152] have assumed that the Yukawa sector of the model is of type-II. In this class of models, the bounds from FCNC processes and LHC searches for heavy Higgses are strong, pushing the masses of the additional 2HDM Higgses above the $500\mathrm{GeV}$ range. It is well known (cf. for example [157,158]) that in fermiophobic 2HDM models of type-I the constraints on the additional Higgs boson can be significantly relaxed, thereby allowing for new scalars and pseudoscalars with masses of order of the EW scale. The fermiophobic nature of the Higgs bosons leads to unconventional production mechanism and also the decays of the non-SM spin-0 particles can have unfamiliar patterns. The mono-X phenomenology in fermiophobic 2HDM+a models awaits explorations.

3.2. 2HDM+s Model

Instead of mixing an additional CP-odd singlet P with the 2HDM pseudoscalar A, as performed in (16), it is also possible to mix a scalar singlet S with the CP-even spin-0 states h and H. This gives rise to the class of 2HDM+s models [34,35], which are the natural, gauge invariant extension of the simplified DM models with a single scalar mediator [36,37]. As in the case of the 2HDM+a model, the presence of non-SM Higgs bosons in the 2HDM+s model can lead to novel $E_T^{\mathrm{miss}}$ that are not captured by a DM simplified model with just a single scalar mediator. While the CP nature of the mediator in the 2HDM+s model typically leads to large SI DM–nucleon interactions that are problematic in view of the stringent DM DD constraints, we give an example of a 2HDM+s model realisation to which the DM DD searches are blind. Detailed discussions of the 2HDM+s model are presented in [34,35,140], and we summarise the most important findings of these articles below.

3.2.1. Theory

The construction of the DM–mediator interactions and the full scalar potential of the 2HDM+s model proceed as in the case of the 2HDM+a model with minor modifications. First, the mediator S couples to the dark scalar current $\chi \overline{\chi }$ and not to the pseudoscalar current $\overline{\chi }{\gamma }_5\chi$. Second, the coupling of the S to the term $(H_1^{†}{H}_2+\mathrm{h}.\mathrm{c}.)$ that appears in (16) is taken to be purely real so that only CP-even components can mix. Alternatively, one can assume that the singlet S develops a VEV and, in this way, obtain a mixing between S and the CP-even scalars of the 2HDM doublets [34,35]. In view of the stringent constraints from the measurement of the Higgs signal strengths at the LHC, it is furthermore natural to work in a generalised alignment limit [34] where only the weak eigenstates H and S mix, giving rise to the mass eigenstates $S_1$ and $S_2$. In the works [34,35,140], the mixing angle $\theta$ between H and S is defined such that for $sin\theta =0$ the state $S_2$ is a pure S state. We employ this definition hereafter as well.

Before discussing the constraints on the 2HDM+s model that arise from mono-X searches, let us also spend some words on the restrictions on the parameter space that stem from DM DD experiments. Expressed in terms of mass eigenstates, the DM–mediator interactions in the 2HDM+s model take the form

$${\mathcal{L}}_{\chi }=-y_{\chi }\left(sin\theta S_1+cos\theta S_2\right)\chi \overline{\chi }.$$

(27)

The couplings between $S_1$ and $S_2$ and the SM quarks depend on the choice of the 2HDM Yukawa sector. In the case of a type-II model, one has for instance

$${\mathcal{L}}_Y=-\frac{cos\theta S_1-sin\theta S_2}{\sqrt{2}}\left[\sum _{q=u,c,t}{y}_{q}cot\beta \overline{q}q-\sum _{q=d,s,b}{y}_{q}tan\beta \overline{q}q\right].$$

(28)

The SI DM–nucleon cross-section takes the general form

$${\sigma }_{\mathrm{SI}}={\left(\frac{m_N{m}_{\chi }}{m_N+m_{\chi }}\right)}^2\frac{c_N^2}{\pi },$$

(29)

where $c_N$ is the Wilson coefficient of the dimension-six nucleon operator $O_N=\chi \overline{\chi }\overline{N}N$ that can be found by integrating out the mediators $S_1$ and $S_2$. In the case of (27) and (28), this coefficient reads [34,35,140]

$$\begin{array}{cc}\hfill c_N& =\frac{m_N}{v}\frac{y_{\chi }sin\left(2\theta \right)}{2}\left(\frac{1}{m_{S_1}^2}-\frac{1}{m_{S_2}^2}\right)\hfill \\ \hfill & \phantom{xx}\times \left[cot\beta f_{T_u}^N-tan\beta \sum _{q=d,s}{f}_{T_q}^N+\frac{4cot\beta -2tan\beta }{27}{f}_{T_G}^N\right],\hfill \end{array}$$

(30)

where $f_{T_u}^N\simeq 0.019$, $f_{T_d}^N\simeq 0.045$ and $f_{T_s}^N\simeq 0.043$ [75,76,77,78] are the quark–nucleon matrix elements and $f_{T_G}^N=1-{\sum }_{q=u,d,s}{f}_{T_q}^N\simeq 0.89$ is the effective gluon–nucleon coupling. Besides the obvious possibilities to suppress the Wilson coefficient $c_N$, i.e., $y_{\chi }\to 0$, $\theta \to 0$ or $m_{S_1}\to m_{S_2}$, in the case of the 2HDM of type-II, there is a fourth possibility to achieve that $c_N\simeq 0$. The trick is to choose $tan\beta$ such that the square bracket in (30) vanishes. Numerically, this happens for $tan\beta \simeq 1$, and this cancellation is only possible because the up- and down-type contributions to $c_N$ interfere destructively in the case of the 2HDM of type-II. Note that in a scenario where a singlet scalar mixes with the $125\mathrm{GeV}$ Higgs to provide the DM portal [98,99,100,101,102,103,104] or a 2HDM model of type-I, interference effects between different quarks are not possible as all fermion couplings to the mediator are proportional to the SM Yukawa couplings with the same proportionality coefficient. The above discussion shows that the Yukawa freedom introduced by the 2HDM is an important feature of the 2HDM+s model, because it can be used to mitigate the stringent DM DD constraints. As a result, it is possible to obtain the correct relic density without having excessive contributions to ${\sigma }_{\mathrm{SI}}$. See [35] for a detailed discussion of this point.

3.2.2. Lhc Constraints

While ATLAS and CMS have not published interpretations of the DM searches in the context of the 2HDM+s model, theoretical reinterpretations of LHC searches have been performed. In the article [140], the ATLAS and CMS searches for the $h\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$ [144,177] and the $Z\left({\ell }^{+}{\ell }^{-}\right)+E_T^{\mathrm{miss}}$ [178,179] final states were examined, and it was found that, as in the case of the 2HDM+a model, these signatures provide the best coverage of the parameter space of the model. This is again related to the fact that the mono-Higgs and mono-Z signatures receive resonant contributions. The corresponding ggF production channels are displayed in Figure 11. Notice that resonant contributions to the $tW+E_T^{\mathrm{miss}}$ signal also exist in the 2HDM+s model. The relevant graph is obtained from the left one in the lower row in Figure 5 by replacing the $H^{-}a{W}^{-}$ vertex by a $H^{-}{S}_2{W}^{-}$ vertex. The $tW+E_T^{\mathrm{miss}}$ signature has so far not been studied in the context of the 2HDM+s model.

The analysis performed in [140] employs the generalised alignment limit and uses the following benchmark parameter choices:

$$\begin{array}{cc}\hfill & m_A=m_{S_1}=m_{H^{\pm }},m_{\chi }=10\mathrm{GeV},sin\theta =0.3,\hfill \\ \hfill & \phantom{xx}{y}_{\chi }=1,{\mathsf{\lambda }}_3=\frac{m_h^2}{v^2}\simeq 0.26,{\mathsf{\lambda }}_{S1}={\mathsf{\lambda }}_{S2}=0,\hfill \end{array}$$

(31)

where ${\mathsf{\lambda }}_{S1}$ and ${\mathsf{\lambda }}_{S2}$ are the analogues of ${\mathsf{\lambda }}_{P1}$ and ${\mathsf{\lambda }}_{P2}$ introduced in (16). The Yukawa sector of the 2HDM is taken to be of type-II. While these parameter choices have common features with the 2HDM+a benchmark (22), the low value of ${\mathsf{\lambda }}_3$ and the vanishing quartic couplings ${\mathsf{\lambda }}_{S_1}$ and ${\mathsf{\lambda }}_{S_2}$ lead to a modified hierachy of the mono-X signatures with respect to the 2HDM+a model. This happens because the $S_1{S}_{2}h$ coupling like $g_{Aah}$ in (20) involves a term proportional to ${\mathsf{\lambda }}_{S_1}{cos}^2\beta +{\mathsf{\lambda }}_{S_2}{sin}^2\beta$, while the $A{S}_{2}Z$ coupling similar to $g_{HaZ}$ in (25) does not depend on the quartic couplings. The parameter choices (31) thus favours a mono-Z signal over a mono-Higgs signature, while in the case of the 2HDM+a model, the opposite is the case because of the sizeable quartic couplings employed in the benchmark (22).

In Figure 12, we summarise the most relevant constraints on the 2HDM+s model in the $m_{S_2}$–$m_A$ plane (left panel) and the $m_{S_2}$–$tan\beta$ plane (right panel). The shown exclusions have been obtained in [140] by recasting the ATLAS searches for $h\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$ [177] and $Z\left({\ell }^{+}{\ell }^{-}\right)+E_T^{\mathrm{miss}}$ [178]. As a result of the choice of the quartic couplings in (31), the mono-Z search provides compared to the mono-Higgs search far better constraints in both two-dimensional parameter planes. For comparison, the exclusions in the 2HDM+a model using the same input are also shown in Figure 12. One observes that the mono-Z constraints are similar in shape and reach in both models. This feature is readily understood by noting that the relevant processes that lead to the mono-Z signal in the two models are $gg\to A\to S_{2}Z$ and $gg\to H\to aZ$. However, for a scalar and pseudoscalar of the same mass and with the same couplings to quarks, one has approximately $\sigma \left(gg\to A\right)/\sigma \left(gg\to H\right)\simeq 2$ and $\mathrm{BR}\left(A\to S_{2}Z\right)/\mathrm{BR}\left(H\to aZ\right)\simeq 1/2$ [140], which leads to the same signal strength in both cases. For the mono-Higgs signal, the relevant resonant processes are instead $gg\to S_1\to S_{2}h$ and $gg\to A\to ah$. In this case, one has however that $\sigma \left(gg\to S_1\right)/\sigma \left(gg\to A\right)\simeq 1/2$ and $\mathrm{BR}\left(S_1\to S_{2}h\right)/\mathrm{BR}\left(A\to ah\right)\simeq 1$ [140], and as a result for the same parameters, the 2HDM+a model predicts a higher mono-Higgs cross-section than the 2HDM+s model. This results in stronger exclusion bounds in the former than in the latter case, as evident from Figure 12.

Since the 2HDM+s interactions generically lead to large SI DM–nucleon cross-sections, we show in the two panels of Figure 12 also the restrictions on the parameter space that stem from DM DD experiments. The shown results employ (29), (30) and (31). In the case of the scan in the $m_{S_2}$–$m_A$ plane (left panel), the latest XENON1T results [79] do not lead to any constraint because of the parameter choice $tan\beta =1$. The situation is different for the $m_{S_2}$–$tan\beta$ plane (right panel) where $tan\beta$ values sufficiently different from one are excluded at 90% CL by the XENON1T measurement, as indicated by the solid blue line. Notice that the DD constraint becomes weaker for increasing mass $m_{S_2}$. All these features can be understood from the structure of the coefficient $c_N$ as given in (30). The scans in Figure 12 demonstrate that, in contrast to naive expectation, in the 2HDM+s model of type-II, it is possible to have interesting LHC signatures that a compatible with the severe limits imposed by DM DD experiments if $tan\beta =\mathcal{O}\left(1\right)$.

The above comparison of the mono-X phenomenology in the 2HDM+s and 2HDM+a models furthermore stresses the complementarity of mono-Higgs and mono-Z searches in exploring DM two-mediator models (see also [133] for a related discussion). Unfortunately, the 2HDM+s model has largely escaped the attention of both the theoretical and experimental community. Neither recommendations a la [136] exist for the 2HDM+s model nor is there a single 2HDM+s interpretation by ATLAS and CMS of the plethora of LHC DM searches. In view of the fact that the articles [34,35,140] give a detailed account of the relevant $E_T^{\mathrm{miss}}$ signatures, the DM DD and ID bounds and the relic density calculation, one would hope that combined analyses such as [152] are conducted in the future in the case of the 2HDM+s model as well.

4. Portals with Extended Higgs and Gauge Sectors

In Section 3, we have discussed two-mediator models where both portal particles have spin-0. Constructing renormalisable interactions between the dark and the visible sector that involve a spin-0 as well as a spin-1 state is however also possible. Given that the interactions between DM and SM particles are severely constrained experimentally, such construction typically involves extensions of both the Higgs and gauge sector of the SM. The resulting theories fall into the class of dark $Z^{\prime }$ or dark Higgs models, which typically have a rich collider and DM phenomenology. In the following section, we discuss two of these models in more detail. The choice of models is motivated by the fact that ATLAS and CMS have already searched for the discussed DM models.

4.1. 2HDM+$Z^{\prime }$ Model

The dark $Z^{\prime }$ model that has received by far the most attention at the LHC (see for example [144,149,150,177,180,181,182,183]) is the so-called 2HDM+$Z^{\prime }$ models introduced in [38]. In this model, only the right-handed up-type quarks are assumed to be charged under the $U{\left(1\right)}_{Z^{\prime }}$ symmetry, while all the other SM fermion fields are chosen as SM singlets. This choice allows LHC production of the $Z^{\prime }$ boson in $u\overline{u}$- and $c\overline{c}$-fusion but avoids the stringent constraints from searches for dilepton resonances since the charged leptons carry no $U{\left(1\right)}_{Z^{\prime }}$ charge. The Higgs sector of the 2HDM+$Z^{\prime }$ model is taken to be a 2HDM of type-II. To incorporate DM interactions, it is assumed that the heavy CP-odd pseudoscalar A that arises from the 2HDM possesses a large coupling to DM particles, such that the corresponding branching ratio is close to one. Models in which the DM candidate $\chi$ is a spin-1/2 or a spin-0 particle have been sketched in [38]. In the first case, $\chi$ is a Majorana fermion that arises from singlet–doublet DM [184,185,186,187], while in the second case, $\chi$ is the lightest component of a complex scalar field. To avoid the stringent constraints from invisible Higgs boson decays and/or DM direct detection experiments, the fundamental parameters in the underlying models that give rise to $\chi$ need to be tuned, but after tuning, it is possible to obtain large branching ratios of A to DM for $\chi$ masses in the ballpark of $100\mathrm{GeV}$ in both cases. For concreteness, we assume hereafter that the DM candidate in the 2HDM+$Z^{\prime }$ model is a Majorana fermion that couples to the A with the coupling strength $g_{\chi }$. The same assumption is made in the ATLAS and CMS analyses [144,149,150,177,180,181,182,183]. In addition, we focus our discussion on the LHC phenomenology of the 2HDM+$Z^{\prime }$ model, ignoring the DM phenomenology because it is more model dependent.

4.1.1. Theory

The tight constraints from the Higgs coupling measurements at the LHC are avoided most easily by working in the alignment limit of the 2HDM+$Z^{\prime }$ model. In this limit, the partial decay widths of the $Z^{\prime }$ boson that are relevant for the discussion hereafter take the following form

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \Gamma \left(Z^{\prime }\to q\overline{q}\right)& =\frac{g_{Z^{\prime }}^2{m}_{Z^{\prime }}}{32\pi }{\left(1-4x_{q/Z^{\prime }}\right)}^{1/2}\left(1-x_{q/Z^{\prime }}\right),\hfill \\ \hfill \Gamma \left(Z^{\prime }\to Ah\right)& =\frac{g_{Z^{\prime }}^2{sin}^2\left(2\beta \right)m_{Z^{\prime }}}{768\pi }{\left(x_{A/Z^{\prime }}^2+{(1-x_{h/Z^{\prime }})}^2-2x_{A/Z^{\prime }}(1+x_{h/Z^{\prime }})\right)}^{3/2},\hfill \\ \hfill \Gamma \left(Z^{\prime }\to Zh\right)& =\frac{g_{Z^{\prime }}^2{sin}^4\beta m_{Z^{\prime }}}{192\pi }{\left(x_{Z/Z^{\prime }}^2+{(1-x_{h/Z^{\prime }})}^2-2x_{Z/Z^{\prime }}(1+x_{h/Z^{\prime }})\right)}^{1/2}\hfill \\ & \phantom{xx}\times \left(x_{Z/Z^{\prime }}^2+{(1-x_{h/Z^{\prime }})}^2-2x_{Z/Z^{\prime }}(1+x_{h/Z^{\prime }})+12x_{Z/Z^{\prime }}\right).\hfill \end{array}\end{array}$$

(32)

Here, $g_{Z^{\prime }}$ denotes the $U{\left(1\right)}_{Z^{\prime }}$ gauge coupling, and we have used the abbreviations $x_{i/j}=m_i^2/m_j^2$ and have neglected contributions that are suppressed by the mass ratio $m_Z^2/m_{Z^{\prime }}^2$ of the Z-boson and the $Z^{\prime }$-boson mass. Such terms arise from Z–$Z^{\prime }$ mixing but are numerically subdominant in the partial decay rates of the $Z^{\prime }$ boson. Notice that the $Z^{\prime }$ boson can only decay to up-quark and charm-quark pairs as well as to top-quark pairs if $m_{Z^{\prime }}>2m_t$ with $m_t\simeq 173\mathrm{GeV}$ denoting the top-quark mass, but not to the other SM fermions because these matter fields are assumed to be $U{\left(1\right)}_{Z^{\prime }}$ singlets.

If the A is sufficiently heavy, the coupling $g_{\chi }$ is sizable and $tan\beta$ is small, important decays of the heavy pseudoscalar are to DM and top-quark pairs. In the alignment limit, the corresponding partial decay rates are given by

$$\begin{array}{cc}\hfill \Gamma \left(A\to \chi \overline{\chi }\right)& =\frac{g_{\chi }^2{m}_A}{8\pi }{\left(1-4x_{\chi /A}\right)}^{1/2},\hfill \\ \hfill \Gamma \left(A\to t\overline{t}\right)& =\frac{3y_t^2{cot}^2\beta m_A}{16\pi }{\left(1-4x_{t/A}\right)}^{1/2}.\hfill \end{array}$$

(33)

Notice that the partial decay width $A\to Zh$ vanishes in the alignment limit, but depending on the precise value of $tan\beta$ and/or the masses of the heavy 2HDM spin-0 states the decays $A\to b\overline{b}$, $A\to {\tau }^{+}{\tau }^{-}$, $A\to HZ$ and $A\to H^{\pm }{W}^{\mp }$ can be relevant in this limit.

From the expressions for the partial decay widths (32) and (33), it follows that if a heavy $Z^{\prime }$ boson is produced on-shell in the LHC collisions this gives rise to a dijet or a $t\overline{t}$ final state, a mono-Higgs signature or resonant $Zh$ production. The examples of the corresponding Feynman diagrams are shown in Figure 13. In our numerical study of the 2HDM+$Z^{\prime }$ model, we consider all four collider signatures that arise from the exchange of $Z^{\prime }$ bosons and study only 2HDM+$Z^{\prime }$ benchmarks in which the dominant A decay modes are described by the expressions (33).

Other laboratory constraints on the 2HDM+$Z^{\prime }$ model arise from the measurements of the EWPO and flavour physics. In the context of the EWPO, the most important constraint turns out to be the bound on the amount of custodial symmetry breaking as measured by the Peskin–Takeuchi parameter T. In the 2HDM+$Z^{\prime }$ model, one obtains at tree level a strictly positive correction to T of the form

$$T=\frac{{sin}^2\left(2{\theta }_w\right)g_{Z^{\prime }}^2{sin}^4\beta }{16\pi {\alpha }^2}\frac{x_{Z/Z^{\prime }}}{1-x_{Z/Z^{\prime }}},$$

(34)

where $sin{\theta }_w$ denotes the sine of the weak mixing angle ${\theta }_w$ with ${sin}^2{\theta }_w\simeq 0.23$, respectively, and $\alpha \simeq 1/128$ is the electromagnetic coupling constant in the $\overline{\mathrm{MS}}$ scheme at the mass of the Z boson. A simultaneous fit to the EWPO gives [3]

$$T=0.03\pm 0.12,$$

(35)

and this constraint leads to a lower bound on the mass $m_{Z^{\prime }}$ of the $Z^{\prime }$ boson for fixed values of $g_{Z^{\prime }}$ and $tan\beta$. Constraints on the masses of the extra spin-0 states of the 2HDM in general arise from the T parameter once loop corrections are included. However, such effects are absent if all non-SM Higgs bosons are taken to be mass degenerate, i.e., if $m_A=m_H=m_{H^{\pm }}$ is assumed, because in this case, the 2HDM Higgs potential (14) is custodially invariant (see the related discussion in Section 3.1). Below, we only consider this limit.

As explained in Section 3.1 indirect constraints on the charged Higgs-boson mass arise from flavour physics with $B\to X_s\gamma$ providing a particularly strong 95% CL lower limit of $m_{H^{\pm }}>800\mathrm{GeV}$ in the case of the 2HDM of type-II. Direct searches for heavy Higgs bosons can also be used to explore the parameter space of the 2HDM+$Z^{\prime }$ model but turn out to be less restrictive than the $Z^{\prime }$-boson search strategies discussed below.

4.1.2. Summary Plots

The existing LHC analyses [144,149,150,177,180,181,182,183] that consider the 2HDM-$Z^{\prime }$ model employ the alignment limit and adopt the following parameter choices:

$$g_{\chi }=1,g_{Z^{\prime }}=0.8,tan\beta =1,m_{\chi }=100\mathrm{GeV},m_A=m_H=m_{H^{\pm }}.$$

(36)

This benchmark has been established in [37], and we also focus on it in what follows. The phenomenological impact that modifications of (36) have, however, is also discussed briefly. Before analysing the relevant LHC constraints, we remark that a combination of (34) and (35) leads for (36) to the 95% CL bound

$$m_{Z^{\prime }}>1080\mathrm{GeV}.$$

(37)

Notice that this bound can be relaxed by decreasing the values of $g_{Z^{\prime }}$ and/or $tan\beta$ compared to (36), but such parameter choices generically also reduce the signal strengths of the LHC signatures. The T parameter therefore always provides a relevant constraint in the 2HDM-$Z^{\prime }$ model in the parameter space that is testable at the LHC.

Dedicated searches for the 2HDM-$Z^{\prime }$ model have been carried out by ATLAS and CMS in mono-Higgs final states with the Higgs decaying to bottom quark [144,150,177,180,182], photon [149,181] or tau pairs [181] as well as a combination of all relevant Higgs decay channels [183]. The dominant contribution to the mono-Higgs signal arises in the 2HDM+$Z^{\prime }$ model from the graph shown in the middle of Figure 13 with subleading effects stemming from $c\overline{c}$-fusion and $pp\to Z^{\prime }\to Zh$ followed by $Z\to \nu \overline{\nu }$.

Mono-Higgs searches are, however, not the only relevant collider constraints that need to be considered if one wants to obtain a global picture of the allowed parameter space in the 2HDM-$Z^{\prime }$ model. That this is the case can simply be seen by evaluating the partial decay width of the $Z^{\prime }$ boson. Using (36) together with $m_{Z^{\prime }}=2\mathrm{TeV}$ and $m_A=800\mathrm{GeV}$ and employing exact formulas for all possible $Z^{\prime }$ decay channels give

$$\begin{array}{cc}\hfill & \sum _{q=u,c}\mathrm{BR}\left(Z^{\prime }\to q\overline{q}\right)=62.8\%,\mathrm{BR}\left(Z^{\prime }\to t\overline{t}\right)=30.7\%,\hfill \\ \hfill & \mathrm{BR}\left(Z^{\prime }\to Ah\right)=0.8\%,\mathrm{BR}\left(Z^{\prime }\to Zh\right)=1.3\%,\hfill \end{array}$$

(38)

in good agreement with the approximations (32). One furthermore has

$$\mathrm{BR}\left(A\to \chi \overline{\chi }\right)=41.8\%,\mathrm{BR}\left(A\to t\overline{t}\right)=57.7\%.$$

(39)

The numbers (38) and (39) for the branching ratios of the $Z^{\prime }$ and the A indicate that besides the mono-Higgs channel, LHC searches for dijet final states (see [188,189] for the latest ATLAS and CMS results), $t\overline{t}$ resonances (see for example [190,191]) and resonant $Zh$ production [192,193,194,195,196,197] can also be expected to lead to relevant constraints in the $m_{Z^{\prime }}$–$m_A$ plane for benchmark parameter choices such as (36). To our knowledge, ATLAS and CMS interpretations of the existing dijet, $t\overline{t}$ and $Zh$ resonance searches in the 2HDM-$Z^{\prime }$ model do not exist, but recasts are straightforward as we show below.

In Figure 14, we summarise the relevant 95% CL constraints on the 2HDM-$Z^{\prime }$ model in the $m_{Z^{\prime }}$–$m_A$ plane. The shown exclusions correspond to the benchmark parameter choices (36). The blue and green exclusions represent the latest mono-Higgs constraints by ATLAS [150] and CMS [182]. The searches $h\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$ provide stronger constraints than $h\left(\gamma \gamma \right)+E_T^{\mathrm{miss}}$ because $\mathrm{BR}\left(h\to b\overline{b}\right)/\mathrm{BR}\left(h\to \gamma \gamma \right)\simeq 250$. The sensitivity decreases with increasing $m_A$ due to other decay channels such as $A\to t\overline{t}$ becoming kinematically accessible, thereby reducing BR$(A\to \chi \overline{\chi })$. It should be noted that while the ATLAS search (solid green line) uses the Higgs invariant mass as a discriminant to separate signal from background events, the CMS search (dashed green line) employs the transverse mass of the $b\overline{b}+E_T^{\mathrm{miss}}$ system, which is better suited for resonant signatures. As a consequence, the CMS search is competitive with the ATLAS one despite using not $139{\mathrm{fb}}^{-1}$ but only $36{\mathrm{fb}}^{-1}$ of integrated luminosity collected at $13\mathrm{TeV}$.

The solid orange, dashed cyan and solid red vertical lines indicate the 95% CL limits on $m_{Z^{\prime }}$ that follow from our recast of the dijet search [188], the $t\overline{t}$ search [190] and the resonant $Zh$ search [197], respectively. These bounds are essentially independent of the precise value of the mass of the A, and for $m_A\ll m_{Z^{\prime }}$ read

$$m_{Z^{\prime }}>3690\mathrm{GeV},m_{Z^{\prime }}>3520\mathrm{GeV},m_{Z^{\prime }}>2820\mathrm{GeV}.$$

(40)

Using [189,191] would lead to a slightly weaker dijet and $t\overline{t}$ limit, respectively. Notice that the exclusions (40) are significantly more stringent than the limit (37), which follows from the T parameter. The latter limit is indicated by a dashed brown line in Figure 14. In fact, the obtained dijet and $t\overline{t}$ constraint are so strong that they exceed the maximal mass reach of the existing mono-Higgs searches that reads

$$m_{Z^{\prime }}>3190\mathrm{GeV},$$

(41)

and is obtained for light pseudoscalar masses. As can be seen from the dashed grey horizontal line, such small values of $m_A$ are under the assumption $m_A=m_H=m_{H^{\pm }}$ in conflict with $B\to X_s\gamma$ if the Yukawa sector of the 2HDM-$Z^{\prime }$ model is taken to be of type II as conducted in [38] and all existing LHC analyses [144,149,150,177,180,181,182,183]. Notice that the recent ATLAS $Zh$ search [197] considers the combination of the ${\ell }^{+}{\ell }^{-}b\overline{b}$ and $\nu \overline{\nu }b\overline{b}$ final states to set limits. In our recast, we have only included contributions with a $Z^{\prime }Zh$ vertex that lead to both of these final states (cf. the right diagram in Figure 13). Contributions with a $Z^{\prime }Ah$ vertex where the pseudoscalar decays via $A\to \chi \overline{\chi }$ have instead been neglected (cf. the middle diagram in Figure 13). While such diagrams mimic the $\nu \overline{\nu }b\overline{b}$ final state, they lead to a different $E_T^{\mathrm{miss}}$ distribution than the former contributions with $Z\to \nu \overline{\nu }$ and thus to a different experimental acceptance. Neglecting diagrams with an internal A in the reinterpretation is however always a conservative approach, because this contribution necessarily leads to a larger signal strength.

In order to highlight the complementary of dijet, $t\overline{t}$ and $Zh$ searches in probing the 2HDM-$Z^{\prime }$ parameter space, we show in Figure 15 the 95% CL exclusions in the $m_{Z^{\prime }}$–$g_{Z^{\prime }}$ plane that follow from our recasts of [188,190,197]. We employ the parameters given in (36) with $m_A=1.5\mathrm{TeV}$ while keeping the $U{\left(1\right)}_{Z^{\prime }}$ coupling $g_{Z^{\prime }}$ as a free parameter. One observes that the used dijet ($t\overline{t}$) search allows one to set better constraints on $g_{Z^{\prime }}$ than the used $Zh$ search in the mass region $m_{Z^{\prime }}\gtrsim 1.5\mathrm{TeV}$ ($m_{Z^{\prime }}\gtrsim 2\mathrm{TeV}$) for lower $Z^{\prime }$-boson masses the situation is reversed. For $m_{Z^{\prime }}\lesssim 1\mathrm{TeV}$, low-mass dijet searches such as [198,199,200,201,202,203] allow one to set further relevant constraints on the 2HDM-$Z^{\prime }$ model. These constraints are however not included in the figure. Notice that the shown LHC constraints are stronger than the indirect constraint from the T parameter (dashed brown line) in the entire mass range shown. We finally add that for the parameter choices employed in Figure 15, the existing mono-Higgs searches do not provide any constraint because the A is too heavy.

The above discussion should have made clear that mono-Higgs searches are not the only way to probe the parameter space of the 2HDM-$Z^{\prime }$ model at the LHC. In fact, as we have argued stringent constraints arise in general from searches that look for resonant $Z^{\prime }$-boson production in final states such as dijets, $t\overline{t}$ or $Zh$. Similar to what has been shown in [146] for the case spin-1 single-mediator DM simplified models, each $Z^{\prime }$-boson resonance search is sensitive to complementary regions of the $m_{Z^{\prime }}$–$g_{Z^{\prime }}$ parameter space, and only by combining the whole suite of searches one is one able to exploit the full LHC potential. ATLAS and CMS interpretations of dijet, $t\overline{t}$ and $Zh$ resonance searches in the 2HDM-$Z^{\prime }$ model, do not exist but could be easily added to the exotics search canon.

4.2. 2MDM Model

While simplified models with a single spin-1 boson mediating the DM–SM interactions may lead to unitarity violation [95,96,97], it has been shown that the introduction of an extra scalar mediator can cure this problem [96]. Furthermore, when this scalar field acquires a VEV, it gives rise to a dark Higgs mechanism, generating masses for all the particles in the dark sector [204] and opening up new channels for DM annihilation, which can allow for the relic density constraints to be met [39,205]. A concrete model implementation of this general idea is discussed in what follows.

4.2.1. Theory

The 2MDM model introduced in [39] includes an additional $U{\left(1\right)}_{Z^{\prime }}$ gauge symmetry with a new gauge boson $Z^{\prime }$ and a complex scalar field S that is a singlet under the SM gauge group. DM is taken to be a Majorana fermion $\chi$ in order to evade the constraints from DM DD experiments. The terms of the Lagrangian that contain the interactions relevant for the discussion below can be written as

$$\begin{array}{cc}\hfill \mathcal{L}& =-\frac{g^{\prime }}{2}{q}_{\chi }{Z}^{\prime \mu }\overline{\chi }{\gamma }_{\mu }{\gamma }_5\chi -\frac{y_{\chi }}{2}\overline{\chi }\left(P_{L}S+P_R{S}^{\ast }\right)\chi \hfill \\ \hfill & \phantom{xx}+4g^{\prime 2}{q}_{\chi }^2{Z}^{\prime \mu }{Z}_{\mu }^{\prime }{S}^{†}S-g^{\prime }{q}_q{Z}^{\prime \mu }\sum _q\overline{q}{\gamma }_{\mu }q.\hfill \end{array}$$

(42)

Here, $g^{\prime }$ is the $U{\left(1\right)}_{Z^{\prime }}$ gauge coupling, $y_{\chi }$ is a Yukawa coupling, $q_{\chi }$ and $q_q$ denote the $U{\left(1\right)}_{Z^{\prime }}$ quantum numbers of the DM and the quark fields, respectively, and $P_L$ and $P_R$ project onto left- and right-handed fields. The $U{\left(1\right)}_{Z^{\prime }}$ charge of S is fixed to be $q_S=-2q_{\chi }$ so that the Yukawa term does not break the symmetry explicitly. In addition, it has been assumed that the SM Higgs doublet carries no $U{\left(1\right)}_{Z^{\prime }}$ charge to avoid mass mixing between the Z and the $Z^{\prime }$ boson which is severely constrained by the measurements of the EWPOs. In such a case, there is also no $Z^{\prime }Zh$ vertex at tree level, and as a result, the searches for resonant $Zh$ production do not provide any constraint. Note that the $Z^{\prime }$ boson has flavour-diagonal and flavour-universal vector couplings to quarks, while it does not couple to leptons. The former feature automatically avoids FCNCs, while the latter feature is crucial in view of the stringent bounds from dilepton searches. Such a charge assignment can be obtained by gauging baryon number [206], and realistic models that implement this mechanism have been discussed for example in [207,208,209,210,211,212]. Models of this type generically have additional fermionic degrees of freedom to cancel gauge anomalies (cf. [213,214,215] for recent discussions), but it is tacitly assumed that the only new fermion that plays an important role at LHC energies is the DM candidate $\chi$.

The scalar S is assumed to acquire a VEV $\langle S\rangle =1/\sqrt{2}{(0,w)}^T$, spontaneously breaking the $U{\left(1\right)}_{Z^{\prime }}$ and thereby giving masses to the $Z^{\prime }$ boson and the DM particles:

$$m_{Z^{\prime }}=2g^{\prime }{q}_{\chi }w,m_{\chi }=\frac{y_{\chi }w}{\sqrt{2}}.$$

(43)

Introducing now the coupling strengths $g_{\chi }=g^{\prime }{q}_{\chi }$ and $g_q=g^{\prime }{q}_q$, the Lagrangian (42) can be rewritten as

$$\begin{array}{cc}\hfill \mathcal{L}& =-\frac{g_{\chi }}{2}{Z}^{\prime \mu }\overline{\chi }{\gamma }_{\mu }{\gamma }_5\chi -\frac{g_{\chi }{m}_{\chi }}{m_{Z^{\prime }}}s\chi \overline{\chi }+2g_{\chi }^2{Z}^{\prime \mu }{Z}_{\mu }^{\prime }\left(\frac{m_{Z^{\prime }}}{g_{\chi }}s+s^2\right)-g_q{Z}^{\prime \mu }\sum _q\overline{q}{\gamma }_{\mu }q,\hfill \end{array}$$

(44)

where s is the scalar field excitation around the VEV of S. Besides the interaction terms given in (44), the full 2MDM model contains further interactions. First, there is the possibility of kinetic mixing between the $U{\left(1\right)}_Y$ and the $U{\left(1\right)}_{Z^{\prime }}$ gauge bosons. The kinetic mixing is constrained to be very small by the measurements of EWPOs (see for example [39]) and therefore ignored in the following. Second, the dark Higgs and the $125\mathrm{GeV}$ Higgs can mix with the amount of mixing conventionally parameterised by $sin\theta$. The good agreement between the Higgs signal strength measurements at the LHC and the SM expectations (cf. [106,107]) requires $|sin\theta |<0.25$ at 95% CL if no additional Higgs decay channels are open. For the values of $sin\theta$ considered below, this bound plays no role, and we also show that additional Higgs decay channels are not an issue.

4.2.2. Experimental Constraints

The LHC searches [216,217,218] have provided interpretations in the 2MDM model for

$$g_{\chi }=1,g_q=0.25,sin\theta =0.01,m_{\chi }=200\mathrm{GeV},$$

(45)

presenting exclusion limits in the $m_{Z^{\prime }}$–$m_s$ plane. The used value of $sin\theta$ clearly satisfies the aforementioned limit on $sin\theta$ imposed by Higgs physics. Furthermore, the decay channels $h\to \chi \overline{\chi }$ and $h\to Z^{\prime }{Z}^{\prime }$ are kinematically closed. Since the mass $m_s$ is scanned over in the works [216,217,218], on-shell Higgs decays of the form $h\to ss$ are in principle allowed if $m_s<m_h/2$. In the limit $sin\theta \ll 1$, the corresponding partial decay width can be approximated by

$$\Gamma \left(h\to ss\right)=\frac{g_{\chi }^2{sin}^2\theta m_h^3}{8\pi m_{Z^{\prime }}^2}{\left(1-4x_{s/h}\right)}^{1/2}{\left(1+2x_{s/h}\right)}^2.$$

(46)

For the parameters (45) and taking $m_{Z^{\prime }}=1\mathrm{TeV}$ and $m_s\ll m_{Z^{\prime }}$, one finds numerically $\Gamma \left(h\to ss\right)=1.9·{10}^{-3}{\Gamma }_h^{\mathrm{SM}}$. Such a small modification easily passes the bounds from the direct measurements of the total Higgs width [159,160] as well as the numerous limits that stem from searches for $h\to 4f$ (see for example [219] for a recent detailed discussion of these bounds). It follows that Higgs physics leads to no restrictions on the 2MDM model for the parameter choices (45).

In Figure 16, we show two example diagrams that give rise to relevant LHC signatures in the 2MDM model. These diagrams only involve the exchange of the $Z^{\prime }$ boson and lead to a signature with two light- or heavy-quark jets in the final state or to a classic mono-jet signature where the jet arises as initial-state radiation (ISR). These topologies can hence be probed for instance via dijet, $t\overline{t}$ or $j+E_T^{\mathrm{miss}}$ searches. Since all these signatures result from s-channel $Z^{\prime }$-boson exchange, an important ingredient that determines the signal strength in a certain channel is the branching ratio of the $Z^{\prime }$ boson into the corresponding final state. The necessary partial decay widths take the following form:

$$\begin{array}{cc}\hfill \Gamma \left(Z^{\prime }\to \chi \overline{\chi }\right)& =\frac{g_{\chi }^2{m}_{Z^{\prime }}}{24\pi }{\left(1-4x_{\chi /Z^{\prime }}\right)}^{3/2},\hfill \\ \hfill \Gamma \left(Z^{\prime }\to q\overline{q}\right)& =\frac{g_q^2{m}_{Z^{\prime }}}{4\pi }{\left(1-4x_{q/Z^{\prime }}\right)}^{1/2}\left(1+2x_{q/Z^{\prime }}\right).\hfill \end{array}$$

(47)

Using the input parameters given in (45) together with $m_{Z^{\prime }}=1\mathrm{TeV}$, these expressions lead to $\mathrm{BR}\left(Z^{\prime }\to \chi \overline{\chi }\right)=25.5\%$, $\mathrm{BR}\left(Z^{\prime }\to q\overline{q}\right)=\mathrm{BR}\left(Z^{\prime }\to t\overline{t}\right)=12.4\%$ and a relative width of ${\Gamma }_{Z^{\prime }}/m_{Z^{\prime }}=4.0\%$. These branching ratios are rather similar to the values of the branching ratios that are predicted in the spin-1 simplified DM benchmark models with $g_{\chi }=1$, $g_q=0.25$ and $g_{\ell }=0$. As a result, one expects to find the same hierarchy of sensitivities as in the spin-1 simplified DM case—see [220,221] for the latest DM summary plots for s-channel mediators by ATLAS and CMS—with the non-$E_T^{\mathrm{miss}}$ searches providing better bounds on $m_{Z^{\prime }}$ than the $E_T^{\mathrm{miss}}$ searches. We see below that this naive expectation is in fact correct.

The 2MDM model however also gives rise to signatures not present in the spin-1 simplified DM models. These novel signatures are illustrated in  Figure 17. Besides the $Z^{\prime }$ boson, the displayed diagrams contain a dark Higgs s that can either be emitted in a $Z^{\prime }{Z}^{\prime }s$ vertex (left graph) or a $\chi \overline{\chi }s$ vertex (right graph). The dark Higgs can decay into SM and DM particles. Assuming that the decay channel $s\to \chi \overline{\chi }$ is kinematically inaccessible, the relevant partial decay widths are given in the limit $sin\theta \ll 1$ by

$$\begin{array}{cc}\hfill \Gamma \left(s\to b\overline{b}\right)& =\frac{3y_b^2{sin}^2\theta m_s}{16\pi }{\left(1-4x_{b/s}\right)}^{3/2},\hfill \\ \hfill \Gamma \left(s\to W^{+}{W}^{-}\right)& =\frac{{sin}^2\theta m_s^3}{16\pi v^2}{\left(1-4x_{W/s}\right)}^{1/2}\left(1-4x_{W/s}+12x_{W/s}^2\right),\hfill \\ \hfill \Gamma \left(s\to ZZ\right)& =\frac{{sin}^2\theta m_s^3}{32\pi v^2}{\left(1-4x_{Z/s}\right)}^{1/2}\left(1-4x_{Z/s}+12x_{Z/s}^2\right),\hfill \\ \hfill \Gamma \left(s\to hh\right)& =\frac{{sin}^2\theta m_s^3}{32\pi v^2}{\left(1-4x_{h/s}\right)}^{1/2}{\left(1+2x_{h/s}\right)}^2.\hfill \end{array}$$

(48)

Notice that the first three expressions have the same functional form than the corresponding partial decay widths of the SM Higgs boson (see for instance [222]), which is expected because the dark Higgs obtains its SM coupling solely by mixing with the $125\mathrm{GeV}$ Higgs. In consequence, the decay pattern of the mediator s resembles, in the first approximation, those of a SM Higgs boson with mass $m_s$. In the mass range $50\mathrm{GeV}\lesssim m_s\lesssim 140\mathrm{GeV}$ the dominant decay mode is $s\to b\overline{b}$ with branching ratios of around 80–40%. In the mass range $140\mathrm{GeV}\lesssim m_s\lesssim 300\mathrm{GeV}$ the dominant decay modes are to $W^{+}{W}^{-}$, $ZZ$ and $hh$ with $\mathrm{BR}\left(s\to W^{+}{W}^{-}\right)\simeq 100\%$ at $m_s\simeq 170\mathrm{GeV}$ and $\mathrm{BR}\left(s\to W^{+}{W}^{-}\right)\simeq 2\mathrm{BR}\left(s\to ZZ\right)\simeq 2\mathrm{BR}\left(s\to hh\right)\simeq 50\%$ at the upper end of the considered mass range. The latter feature is expected from the $SU{\left(2\right)}_L\times U{\left(1\right)}_Y$ symmetric limit. For $m_s$ values above the top-quark threshold, the decay $s\to t\overline{t}$ is also relevant. In the following numerical analysis, we restrict ourselves to the mass range $50\mathrm{GeV}<m_s<350\mathrm{GeV}$, and thus, the dominant final states that arise from the graphs in Figure 17 are $b\overline{b}+E_T^{\mathrm{miss}}$, $VV+E_T^{\mathrm{miss}}$ with $VV=W^{+}{W}^{-},ZZ$ and $hh+E_T^{\mathrm{miss}}$.

We now turn to the numerical analysis of the constraints that the various LHC searches impose on the parameter space of the 2MDM model. We start by considering the exclusions that can derived from the signatures that result from the diagrams displayed in Figure 16. For the benchmark parameter choices (45), we find by reinterpreting the dijet [188], $t\overline{t}$ [190] and mono-jet search [167] the following 95% CL limits on the $Z^{\prime }$-boson mass:

$$m_{Z^{\prime }}>3680\mathrm{GeV},m_{Z^{\prime }}>3190\mathrm{GeV},m_{Z^{\prime }}>1520\mathrm{GeV}.$$

(49)

The dijet [189], $t\overline{t}$ [191] and mono-jet [168] searches have comparable sensitivities to the searches leading to (49). We emphasise that our recast of the ATLAS mono-jet search includes only $j+E_T^{\mathrm{miss}}$ events with ISR jets (see the right diagram in Figure 16), while we have not considered contributions that follow from $s\left(VV\right)+E_T^{\mathrm{miss}}$ production with the EW gauge bosons decaying fully hadronically (see Figure 17). We comment on this simplification below. As anticipated, the limits (49) have the same hierarchy as the bounds obtained in the standard benchmark considered in spin-1 simplified DM models (cf. [220,221]).

For the $s+E_T^{\mathrm{miss}}$ signatures, it is important to realise that the dark Higgs is produced with high momentum due to the large $m_{Z^{\prime }}$–$m_s$ mass splitting. As a result, the dark Higgs decay products are strongly collimated. According to (48), the dark Higgs decays predominantly to a pair of b-jets for $m_s\lesssim 140\mathrm{GeV}$. This decay mode is targeted by a reinterpretation [216] of the ATLAS $h\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$ search that considers the invariant mass distribution of the $b\overline{b}$ pair down to $50\mathrm{GeV}$, allowing one to probe $m_s>50\mathrm{GeV}$. For $m_s\gtrsim 140\mathrm{GeV}$, the $s\to VV$ decay mode becomes important. ATLAS targeted this decay mode by exploring fully hadronic final states where both V bosons decay into a quark pair each [217]. The challenging multi-prong decay $s\to V\left(q\overline{q}\right)V\left(q\overline{q}\right)$ was reconstructed with a novel technique [223] aimed at resolving the dense topology from boosted $VV$ pairs using reclustered jets in the calorimeter and tracking information. Recently, CMS has probed the same dark Higgs mass range in the $W^{+}{W}^{-}+E_T^{\mathrm{miss}}$ final state where each W boson decays leptonically [218]. The signal was extracted from a two-dimensional fit to the dilepton invariant mass and the transverse mass of the trailing lepton plus $E_T^{\mathrm{miss}}$ system.

The constraints on the $m_{Z^{\prime }}$–$m_s$ plane from the $s+E_T^{\mathrm{miss}}$ analyses by the ATLAS collaboration [216,217] are presented in Figure 18. The 95% CL exclusion limits (49) are also shown in the figure for comparison. Focusing on the $E_T^{\mathrm{miss}}$ searches, first one observes that the novel $s\left(VV\right)+E_T^{\mathrm{miss}}$ and $s\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$ search strategies allow one to exclude additional parameter space of the 2MDM model with respect to our mono-jet recast of the ATLAS search [167]. In this context, one however has to remember that the mono-jet limit given in (49) does not include the contributions of the $s\left(VV\right)+E_T^{\mathrm{miss}}$ and $s\left(b\overline{b}\right)+E_T^{\mathrm{miss}}$ processes. The obtained mono-jet bound therefore provides a conservative lower limit on the actual sensitivity of the $j+E_T^{\mathrm{miss}}$ search [167] for the 2MDM benchmark (45)—we estimate that including $s+E_T^{\mathrm{miss}}$ contributions would strengthen the mono-jet limits by around 15% (5%) at $m_s=100\mathrm{GeV}$ ($m_s=200\mathrm{GeV}$). It is also evident from the figure that our recast of the dijet [188] and $t\overline{t}$ search [190] has a higher mass reach in $m_{Z^{\prime }}$ than the mono-X searches. We note that the 2MDM model can also give rise to $hh+E_T^{\mathrm{miss}}$ and $t\overline{t}+E_T^{\mathrm{miss}}$ signatures. The former signal can be constrained for instance by using the results of [224,225] that studies final states with at least three b-jets and $E_T^{\mathrm{miss}}$ (see also [226] for an exploration of the $hh+E_T^{\mathrm{miss}}$ signature in simplified models of hidden sectors). Estimating the sensitivity of multi b-jet plus $E_T^{\mathrm{miss}}$ searches to the $hh+E_T^{\mathrm{miss}}$ signature in the 2MDM model is however beyond the scope of this review. The latter signal can be targeted by standard $t\overline{t}+E_T^{\mathrm{miss}}$ searches (cf. for example [143,164,166]). Employing [166], we expect that for (45) the HL-LHC with $3{\mathrm{ab}}^{-1}$ may be able to set the bounds $M_{Z^{\prime }}\gtrsim 500\mathrm{GeV}$ for $m_s\simeq 750\mathrm{GeV}$ and $M_{Z^{\prime }}\gtrsim 1.5\mathrm{TeV}$ for $m_s\simeq 400\mathrm{GeV}$. We finally add that the 2HDM model can lead to interesting LLP signatures if the mixing between the SM and the dark Higgs is switched off, because in such a case, the dark Higgs can only decay through EW gauge boson loops. For details, see the recent publication [227].

5. Models with Exotic Higgs Decays Involving LLPs

BSM scenarios with hidden sectors that are connected to the SM sector through the $125\mathrm{GeV}$ Higgs boson are being actively explored at colliders. Such scenarios are often characterised by new electrically neutral LLPs. These LLPs typically decay into SM particles, leaving a displaced vertex signature in the detector. The identification of such signatures is nontrivial and often requires dedicated triggering and reconstruction algorithms—see [228,229] for a detailed review of experimental aspects of LLPs at the LHC. Generally speaking, LLP signatures can arise in a multitude of BSM models ranging from SUSY, theories of neutral naturalness, hidden valley models and composite Higgs theories, to just name a few examples. Depending on the context, the appearance of BSM LLPs can be theoretically motivated for instance by the naturalness problem of the Higgs mass, the DM puzzle, baryogenesis or the smallness of neutrino masses—a comprehensive collection of BSM theories with LLPs can be found in [230]. Below, we consider three hidden sector models in which the LLPs arise from decays of the $125\mathrm{GeV}$ Higgs boson. The collider phenomenology of axion-like particles (see for example [231,232] for recent comprehensive studies) that might also lead to LLP signatures from Higgs decays are not discussed.

5.1. Neutral Naturalness

The discovery of a light, seemingly elementary Higgs boson has escalated the seriousness of the EW hierarchy problem, while the steadily increasing LHC limits on coloured BSM states exclude more and more of the natural parameter space of the standard solution to the EW hierarchy problem such as SUSY or compositeness. Models of neutral naturalness such as twin Higgs [40], folded SUSY [41], quirky little Higgs [42] and orbifold Higgs [43] provide compelling alternative solutions to the EW hierarchy problem. In these theories, the large radiative corrections to the Higgs mass associated with the top quark, are cancelled by top partners that carry no colour, thereby relaxing the most stringent LHC constraints that follow from strong production. The cancellation is achieved with discrete symmetries that must be nearly exact in the top sector but may be approximate for the other partner or mirror states [233]. In fact, since the QCD coupling drives the renormalisation group (RG) running of the top-quark Yukawa coupling, for the near-exact discrete symmetry in the top sector to be preserved, viable theories of neutral naturalness contain at least one new QCD-like hidden gauge group with a coupling ${\widehat{\alpha }}_s$ that is comparable in strength to its SM counterpart ${\alpha }_s$. Connecting DM with neutral naturalness is possible (see for instance [234,235,236,237,238,239,240,241,242,243,244,245]) but more model dependent than the LLP phenomenology on which we focus in the following section.

5.1.1. Theory

In models of neutral naturalness the coupling between the $125\mathrm{GeV}$ Higgs boson and the top partners gives rise to an effective coupling $h\widehat{g}\widehat{g}$ between the $125\mathrm{GeV}$ Higgs h and the hidden gluons $\widehat{g}$ at the one-loop level. This is in full analogy to the SM where top-quark loops provide the dominant contribution to the effective $hgg$ coupling. The effective $h\widehat{g}\widehat{g}$ interactions can be parameterised by

$$\mathcal{L}=\frac{{\widehat{\alpha }}_s}{12\pi }\frac{h}{v}{\widehat{G}}_{\mu \nu }^a{\widehat{G}}^{a,\mu \nu }\widehat{\zeta },$$

(50)

where ${\widehat{G}}_{\mu \nu }^a$ denotes the field strength tensor of the hidden $SU\left(3\right)$ gauge group and $\widehat{\zeta }$ is a model-dependent mixing angle. One generically has $|\widehat{\zeta }|=\mathcal{O}(v^2/M^2)$ where depending on the model M is either the scale of spontaneous global symmetry breaking or the mass of the top partners. Fine-tuning arguments bound the scale M to lie at or below the TeV scale, implying that the relevant mixing angles fall into the range $0.1\lesssim |\widehat{\zeta }|\lesssim 1$.

The coupling (50) provides a portal for production of states in the hidden QCD sector. Once produced, states in the hidden sector cascade down to the lightest accessible mirror state, which is typically a bound state of hidden QCD. This state can then decay back to SM particles with the associated lifetime depending on the exact nature of the hidden sector. A canonical signal that can appear in models of neutral naturalness without light mirror matter, such as folded SUSY and quirky little Higgs and some realisations of twin Higgs, is the formation of mirror glueballs [246,247].

The lightest mirror glueball is a scalar state with quantum numbers $J^{CP}=0^{++}$, and its mass ${\widehat{m}}_0$ is entirely determined by the RG running of the strong coupling ${\widehat{\alpha }}_s$ in the hidden sector. In fact, ${\widehat{m}}_0$ is related to the mass $m_0\simeq 1.7\mathrm{GeV}$ of the $0^{++}$ glueball in QCD [248,249,250] by the following simple rescaling [44]

$${\widehat{m}}_0=\frac{{\widehat{\Lambda }}_{\mathrm{Landau}}}{{\Lambda }_{\mathrm{Landau}}}{m}_0,$$

(51)

with ${\Lambda }_{\mathrm{Landau}}\simeq 150\mathrm{MeV}$ the Landau pole in QCD, i.e., the scale where $1/{\alpha }_s\left({\Lambda }_{\mathrm{Landau}}\right)=0$. Without knowing the exact mass spectrum of states in the hidden sector, it is not possible to give a precise value of ${\widehat{\Lambda }}_{\mathrm{Landau}}$ and therefore ${\widehat{m}}_0$. One can however estimate the typical range of ${\widehat{m}}_0$ values by studying the following toy model. Let $M_{\ast }$ be the scale where the discrete symmetry ${\widehat{\alpha }}_s\left(M_{\ast }\right)={\alpha }_s\left(M_{\ast }\right)$ between the two strong couplings is broken and denote with M the scale of the lightest top partner. The one-loop beta function in the hidden sector takes the form

$$\widehat{\beta }=\left\{\begin{array}{cc}{\displaystyle 11-\frac{2}{3}{\widehat{N}}_f,}& M<\mu <M_{\ast },\\ 11& \mu <M,\end{array}\right.$$

(52)

where ${\widehat{N}}_f$ denotes the number of flavours that are active for renormalisation scales $\mu$ in the range $M<\mu <M_{\ast }$. Solving the one-loop RG equation of ${\widehat{\alpha }}_s$, it is easy to show that the Landau pole in our toy model occurs at

$${\widehat{\Lambda }}_{\mathrm{Landau}}=M_{\ast }{\left(\frac{M_{\ast }}{M}\right)}^{-\frac{2{\widehat{N}}_f}{33}}{e}^{-\frac{2\pi }{11{\alpha }_s\left(M_{\ast }\right)}},$$

(53)

where

$${\alpha }_s\left(M_{\ast }\right)={\alpha }_s\left(m_t\right){\left[1-\frac{7{\alpha }_s\left(m_t\right)}{2\pi }ln\left(\frac{m_t}{M_{\ast }}\right)\right]}^{-1},$$

(54)

with ${\alpha }_s\left(m_t\right)\simeq 0.11$. Assuming now that there is only one light hidden state of mass M, while all other mirror states have a mass $M_{\ast }$, and taking $M=v$ and $M_{\ast }=1\mathrm{TeV}$ or $M=2\mathrm{TeV}$ and $M_{\ast }=20\mathrm{TeV}$, one finds using (51), (53) and (54) the following approximate range of mirror glueball masses:

$$15\mathrm{GeV}\lesssim {\widehat{m}}_0\lesssim 50\mathrm{GeV}.$$

(55)

While more sophisticated calculations (see for instance [44]) lead to slightly larger mass ranges, they do not change the conclusion that a representative mirror sector gives rise to glueballs that can be pair-produced in the decays of the $125\mathrm{GeV}$ Higgs. Glueball production is therefore a smoking-gun signature in many theories of neutral naturalness.

Under the assumption that the $0^{++}$ mirror glueballs are dominantly produced in symmetric two-body Higgs decays, one can estimate the corresponding exclusive Higgs branching ratio. One finds [44]

$$\mathrm{BR}\left(h\to 0^{++}{0}^{++}\right)={\left(1-\frac{4{\widehat{m}}_0^2}{m_h^2}\right)}^{1/2}{\left(\frac{{\widehat{\alpha }}_s\left(m_h\right)}{{\alpha }_s\left(m_h\right)}\widehat{\zeta }\right)}^2\kappa \left({\widehat{m}}_0\right)\mathrm{BR}{\left(h\to gg\right)}_{\mathrm{SM}},$$

(56)

where $\mathrm{BR}{\left(h\to gg\right)}_{\mathrm{SM}}\simeq 8.2\%$ is the Higgs to digluon branching ratio in the SM. The parameter $\kappa \left({\widehat{m}}_0\right)$ encodes our ignorance about the hadronisation of the lightest mirror glueball and the mixing effects of excited hidden glueball states with the $125\mathrm{GeV}$ Higgs. In the article [44], it was argued that $\kappa \left(15\mathrm{GeV}\right)\simeq 0.1$ and $\kappa \left(50\mathrm{GeV}\right)\simeq 1$. Using these numbers as well as the estimate $\widehat{\zeta }\simeq v^2/M^2$, one obtains for the two benchmarks that led to (55) the following range of branching ratios:

$$2·{10}^{-5}\lesssim \mathrm{BR}\left(h\to 0^{++}{0}^{++}\right)\lesssim 1·{10}^{-2}.$$

(57)

In view of all the approximations and estimates that went into (57), the given range should only be taken as an indication of the typical values of the $h\to 0^{++}{0}^{++}$ branching ratio that arise in theories of neutral naturalness.

Since the $0^{++}$ mirror glueball has the same quantum numbers as the SM Higgs, both states mix by virtue of (50). Once produced, the $0^{++}$ mirror glueballs can hence decay to all kinematically available SM particles Y via an off-shell Higgs, i.e., through the process $0^{++}\to h^{\ast }\to YY$. The corresponding partial decay widths have been calculated in the work [247] and in the case of (50) take the form of

$$\Gamma \left(0^{++}\to YY\right)={\left(\frac{{\widehat{\alpha }}_s\left({\widehat{m}}_0\right){\widehat{f}}_0\widehat{\zeta }}{6\pi v\left(m_h^2-{\widehat{m}}_0^2\right)}\right)}^2\Gamma {\left(h^{\ast }\to YY\right)}_{\mathrm{SM}},$$

(58)

where ${\widehat{f}}_0$ is the annihilation matrix element of the lightest mirror glueball through (50) and $\Gamma {\left(h^{\ast }\to YY\right)}_{\mathrm{SM}}$ denotes the partial decay width for a SM Higgs boson with mass ${\widehat{m}}_0$. It follows that the decay pattern of the lightest scalar mirror glueball resembles that of a SM Higgs boson of appropriate mass. For mirror glueball masses in the range (55), the $0^{++}$ decays around 80%, 10% and 10% of the time to $b\overline{b}$, $c\overline{c}$ and ${\tau }^{+}{\tau }^{-}$ final states, respectively. Using again $\widehat{\zeta }\simeq v^2/M^2$ as well as ${\widehat{\alpha }}_s\left({\widehat{m}}_0\right){\widehat{f}}_0\simeq 0.18{\widehat{m}}_0^3$ [44,250] together with (58), the proper decay length of the lightest mirror glueball can be approximated by

$$c{\tau }_{0^{++}}\simeq 2\mathrm{m}·{\left(\frac{15\mathrm{GeV}}{{\widehat{m}}_0}\right)}^7{\left(\frac{M}{1\mathrm{TeV}}\right)}^4.$$

(59)

For typical values of ${\widehat{m}}_0$ and M realised in theories neutral naturalness, the proper decay length of the $0^{++}$ mirror glueball ranges from microns to kilometers. In fact, the strong scaling of (59) with both ${\widehat{m}}_0$ and M suggests that $c{\tau }_{0^{++}}$ can in practice be treated as an almost free parameter in the framework of neutral naturalness.

Notice that hidden valley models [251,252,253] share many of the phenomenological features discussed above. Like in theories of neutral naturalness also in hidden valley models, a new confining gauge group is added to the SM. However, the confining gauge group in hidden valley models makes, in full analogy to QCD, hidden hadrons out of hidden quarks. If the hidden sector comprises two light flavours of quarks, the spectrum of hidden hadrons contains a hidden pion ${\pi }_h$. Given its pseudoscalar nature, the hidden pion preferentially decays to heavy SM flavours, for example ${\pi }_h\to b\overline{b}$. As argued in the articles [251,252,253], the typical ${\pi }_h$ masses and proper decay lengths fall into the ballpark of (55) and (59), respectively. Hidden valley models may also contain hidden Higgses which partake in the mass generation of the hidden quarks. If one of these hidden Higgs fields mixes with the $125\mathrm{GeV}$ Higgs boson, it is possible to obtain $h\to {\pi }_h{\pi }_h$ branching ratios that are observable at the LHC [251,252,253]. The phenomenology of the ${\pi }_h$ is therefore very similar to that of the $0^{++}$ with the most obvious difference that the hidden pion is a pseudoscalar, whereas the lightest mirror glueball is a scalar.

While the above considerations broadly motivate searches for displaced Higgs decays at the LHC, the models presently employed in the interpretation of such searches by the experimental collaborations are more generic than the theories of neutral naturalness or the hidden valley models discussed above. The used simplified models assume SM production of a Higgs boson and its subsequent decay to a pair of scalar (s) or pseudoscalar ($a)$ particles. The s (a) is assumed to decay like the SM Higgs, while its mass $m_s$ ($m_a$), the relevant Higgs branching ratio $\mathrm{BR}\left(h\to ss\right)$($\mathrm{BR}\left(h\to aa\right)$) and its proper decay length $c{\tau }_s$ ($c{\tau }_a$) are treated as free parameters. For what concerns models of neutral naturalness, this approach is motivated by (55), (57) and (59). The scalar (pseudoscalar) case can be thought to cover theories of neutral naturalness (hidden valley models) with the lightest mirror glueball (the hidden pion) being the LLP. For definiteness, we hereafter refer to the LLP produced in exotic Higgs decays as an a particle.

5.1.2. Experimental Constraints

The first searches for pair-produced neutral LLPs in the context of Higgs portal models were performed by the CDF [254] and D [255] collaborations at the Tevatron. Both searches looked for displaced vertices in their tracking system only, thereby setting limits on LLP mean decay lengths of the order of a few centimetres. At the LHC, searches for Higgs decays into LLP were carried out by the ATLAS, CMS and LHCb collaborations in different final states, covering a wide range of mean decay lengths. The LLP mean decay length determines the search strategies and reconstruction techniques that are employed. Below, we discuss the relevant LHC searches, starting with the shortest mean decay lengths considered.

ATLAS has performed searches for the decay $h\to aa\to 4b$ optimised for prompt decays or small proper decay lengths $c{\tau }_a\lesssim 6·{10}^{-3}\mathrm{m}$ [256,257]. The searches select events corresponding to associated $Vh$ production with decays of the EW bosons into leptons, as displayed on the left-hand side of Figure 19. Since the targeted LLP mean decay lengths were small, standard track/vertex reconstruction and b-jet identification techniques were used. LHCb has also performed a search for $gg\to h\to aa$ with a decaying to hadronic jets [258], which is sensitive to small mean decay lengths in the ballpark of a few millimetres. The corresponding signal process is illustrated on the right in Figure 19. At least one displaced vertex was required in the event due to the limited acceptance of the LHCb vertex detector. The data were recorded by requiring the presence of an energetic charged lepton or hadron in the event in the hardware trigger, and either one track with high transverse momentum ($p_T)$ and a large impact parameter, or a displaced vertex of at least two tracks in the software trigger. Further requirements on the displaced vertex properties were applied in the last stage of the software trigger and in the data analysis.

For mean decay lengths in the range of ${10}^{-3}\mathrm{m}$ to $1\mathrm{m}$, a substantial fraction of the LLPs is expected to decay inside the inner detector (ID) of the LHC experiments. This allows for a direct reconstruction of the displaced decay vertex and hence dramatically reduces the SM background rate, which becomes dominated by long-lived hadrons and instrumental backgrounds. At the same time, the macroscopic mean decay lengths significantly decrease the efficiency of standard track reconstruction algorithms and thereby also the identification efficiency of b-jets from LLP decays. Searches relying on displaced vertex signatures in the inner tracker systems of ATLAS and CMS are described in the following section.

In order to significantly increase the efficiency for reconstructing displaced vertices, dedicated reconstruction techniques, so-called large-radius tracking (LRT), were developed at ATLAS [259]. The properties of displaced vertices reconstructed from LRT tracks such as the invariant mass or the number of associated tracks can then be used to discriminate signal from background. Since the LRT reconstruction algorithms are computationally intensive, they were not employed in the ATLAS trigger. Hence, the ATLAS search employing LRT techniques targeted the associated $Zh$ production channel with $h\to aa\to 4b$ and $Z\to {\ell }^{+}{\ell }^{-}$. The leptons from the $Z\to {\ell }^{+}{\ell }^{-}$ decays were used for triggering, and the backgrounds were essentially eliminated by requiring two displaced vertices in candidate events [260].

CMS has recently searched for $h\to aa\to 4b$ and $Z\to {\ell }^{+}{\ell }^{-}$ using the same process [261]. A trigger and selections based on dilepton Z-boson decays provide sensitivity to light LLPs with masses of $15\mathrm{GeV}$ or less. The decays of the LLPs are selected by requiring the presence of displaced jets which are identified using information from the tracking system. CMS in addition searched for the $h\to aa\to 4b$ process targeting $c{\tau }_a$ between ${10}^{-3}\mathrm{m}$ and $3\mathrm{m}$ by requiring displaced vertices consistent with LLP decays to be reconstructed in the inner tracker [262,263]. Unlike the ATLAS analysis [260], the CMS search [263] considered the ggF topology and relied on jets for triggering. Two triggers were used, both requiring a large scalar sum of transverse jet energies in the event and at least two energetic jets consistent with a LLP decay, i.e., at most two associated prompt tracks with $p_T>1\mathrm{GeV}$ and at least one track consistent with originating from a displaced vertex. Benefiting from the large ggF Higgs production cross-section, the CMS search [263] currently provides the best sensitivity for proper decay lengths between ${10}^{-3}\mathrm{m}$ and ${10}^{-1}\mathrm{m}$. CMS has also searched for displaced leptons arising from ggF Higgs production followed by $h\to aa\to 4\ell$ [264]. Since this search assumes that the LLP has equal probability to decay to two muons and two electrons, whereas (58) implies that $\mathrm{BR}\left(a\to {\mu }^{+}{\mu }^{-}\right)=\mathcal{O}\left({10}^{-4}\right)$ and $\mathrm{BR}\left(a\to e^{+}{e}^{-}\right)=\mathcal{O}\left({10}^{-8}\right)$, an interpretation of [264] in the context of neutral naturalness/hidden valleys leads to no meaningful constraint on the $h\to aa$ branching ratio.

For mean decay lengths in the range of $1\mathrm{m}$ to ${10}^2\mathrm{m}$, a significant fraction of the LLPs is expected to decay in the outermost layers of the detector, i.e., inside the calorimeter (CM) or the muon spectrometer (MS). Searches in this regime are described in the following section.

ATLAS searched for $h\to aa\to 4f$ decays in the CM targeting event topologies compatible with ggF Higgs production [265]. The decay of a LLP inside the CM is typically reconstructed as a single jet with striking characteristics, namely a narrow width in the pseudorapidity-azimuthal-angle ($\eta$–$\varphi$) space, a high ratio of energy deposited in the hadronic CM to that registered in the electromagnetic CM, and no or only a few low-momentum tracks associated with it. The search uses these characteristics both in dedicated triggers [266] and when performing the offline analysis employing an artificial neural network. ATLAS recently performed a search for the same process using the MS [267]. The tracking capabilities of the MS allow for an explicit reconstruction of a displaced vertex, which dramatically reduces SM and instrumental backgrounds, and allows for a dedicated triggering strategy [266] based on the overall activity in the MS. The backgrounds are essentially eliminated by requiring two displaced vertices from candidate $h\to aa\to 4b$ decays within the MS. The sensitivity for even longer decay times can be extended by requiring only one displaced vertex to be detected in the MS. This strategy was applied in a similar analysis [268], and the results were statistically combined with the CM search [265]. Finally, ATLAS carried out a search using the combination of the ID and the MS, targeting short and long mean decay lengths [269].

CMS performed a LLP search using the endcap of the MS [270], targeting $c{\tau }_a$ between ${10}^{-3}\mathrm{m}$ and ${10}^2\mathrm{m}$ and decay chains $h\to aa$ with $a\to b\overline{b},c\overline{c},{\tau }^{+}{\tau }^{-}$. While many features are similar to the corresponding ATLAS searches [267,268], the search philosophy differs in one key aspect: the magnetic field return yoke, interleaved with the tracking layers of the MS was used as a sampling CM. Only the endcap region of the CMS MS is considered, as it provides a greater depth of up to 30 nuclear interaction lengths. A dedicated algorithm was applied to cluster the hits in the muon system, and the hit multiplicity in an azimuthal slice close to the $E_T^{\mathrm{miss}}$ vector was used as the final discriminant.

ATLAS and CMS also performed other searches for neutral LLPs decaying to jets targeting SUSY models [271,272,273,274,275] that are not explicitly optimised for $h\to aa\to 4f$ signatures, for example, due to very high trigger thresholds, and therefore are not discussed below.

An indirect way to constrain models where the $125\mathrm{GeV}$ Higgs boson decays into LLPs is provided by the combination of precision measurements of the Higgs couplings in visible and invisible final states [73,106,107]. In particular, under the assumption that the coupling modifiers ${\kappa }_V=g_{hVV}/g_{hVV}^{\mathrm{SM}}$ of the $125\mathrm{GeV}$ Higgs boson to EW gauge bosons satisfy ${\kappa }_V\le 1$, which generically holds in models with Higgs mixing, the following 95% CL constraints

$$\mathrm{BR}\left(h\to \mathrm{undet}\right)<0.19,\mathrm{BR}\left(h\to \mathrm{inv}\right)<0.09,$$

(60)

can be placed [106]. Here, the undetected category includes events with undetected BSM particles that do not provide a significant $E_T^{\mathrm{miss}}$ contribution such as the ones typically selected by the aforementioned LLP analyses with short mean decay lengths below $\mathcal{O}\left(1\mathrm{m}\right)$. By contrast, LLPs that decay outside the tracker and calorimeters, i.e., having mean decay lengths larger than $\mathcal{O}\left(1\mathrm{m}\right)$, are not registered by standard reconstruction algorithms, resulting in a $E_T^{\mathrm{miss}}$ contribution. Hence, for moderate boosts, the above bound on $\mathrm{BR}\left(h\to \mathrm{undet}\right)$ and $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ can be interpreted as an indirect limit on $\mathrm{BR}\left(h\to aa\right)$ for proper decays lengths of $c{\tau }_a\lesssim 1\mathrm{m}$ and $c{\tau }_a\gtrsim 1\mathrm{m}$, respectively. In the case of larger boosts, the boundary between the constraints from undetected and invisible Higgs decays is instead approximately given by $c{\tau }_a\simeq 0.1\mathrm{m}(m_a/\mathrm{GeV})$ [276,277].

The 95% CL limits of the above searches on the branching ratio $\mathrm{BR}\left(h\to aa\right)$ as a function of the proper decay length $c{\tau }_a$ of the LLP are shown in Figure 20. The displayed masses of the LLP a all fall into the central region of $0^{++}$ mirror glueball masses as predicted by neutral naturalness (55). For proper decay lengths $c{\tau }_a$ of a few meters, the LHC searches are able to set a limit $\mathrm{BR}\left(h\to aa\right)\lesssim {10}^{-3}$. In view of (57) and (59), this is an interesting finding, as it allows one to test certain model realisations of neutral naturalness. Notice that the limits from $\mathrm{BR}\left(h\to \mathrm{undet}\right)$ and $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ provide currently the strongest LHC bound on $\mathrm{BR}\left(h\to aa\right)$ for $c{\tau }_a\lesssim {10}^{-4}\mathrm{m}$ and $c{\tau }_a\gtrsim {10}^3\mathrm{m}$, respectively. Besides these two constraints, proper decay lengths $c{\tau }_a\gtrsim {10}^3\mathrm{m}$ are currently unexplored by collider measurements. Dedicated detectors such as MATHUSLA [230], CODEX-b [278] and ANUBIS [279] may address such long mean decay lengths at the HL-LHC.

5.2. Dark Photons

In the case of neutral naturalness, the LLP is a composite spin-0 particle. However, the hidden sector models with a spin-1 LLP also exist. The model considered in this subsection is based on an extra $U{\left(1\right)}_X$ symmetry in the hidden sector, where the associated vector field X is given a mass via a dark Higgs mechanism involving the singlet scalar field S. As explained below, in certain regions of parameter space, this model predicts displaced dilepton vertex signatures that arise from the exotic decays of the $125\mathrm{GeV}$ Higgs boson to a pair of dark photons $Z_d$ followed by $Z_d\to {\ell }^{+}{\ell }^{-}$. This feature makes the discussed hidden sector model experimentally distinct from theories of neutral naturalness where the LLP decay products consist primarily of hadrons. We do not discuss the dark photon searches [280,281,282] which are motivated by the theoretical works [283,284,285,286]. We emphasise that the dark photon models that are discussed in this subsection do not have a DM candidate. For comprehensive discussions of dark photon models with an additional DM candidate, see for example the reviews [287,288,289].

5.2.1. Theory

The interaction terms of the hidden sector model that we consider include both a hypercharge and a Higgs portal. See [45,46] for details and further relevant literature. We write these terms in the following way:

$$\mathcal{L}=\frac{ϵ}{2cos{\theta }_w}{B}_{\mu \nu }{X}^{\mu \nu }-\kappa \left(H^{†}H\right)S^2,$$

(61)

where $B^{\mu \nu }={\partial }^{\mu }{B}^{\nu }-{\partial }^{\nu }{B}^{\mu }$ and $X^{\mu \nu }={\partial }^{\mu }{X}^{\nu }-{\partial }^{\nu }{X}^{\mu }$ is the $U{\left(1\right)}_Y$ and $U{\left(1\right)}_X$ gauge field strength tensor, respectively. $ϵ$ denotes the hypercharge mixing parameter, while $\kappa$ is the Higgs portal coupling. After EW symmetry breaking $\langle H\rangle ={(0,v/\sqrt{2})}^T$ and spontaneous symmetry breaking of the $U{\left(1\right)}_X$ symmetry by $\langle S\rangle =v_S/\sqrt{2}$, the mass spectrum of the hidden sector model contains two heavy neutral gauge bosons and two neutral Higgs bosons. We denote these states by Z, $Z_d$, h and $h_d$. For small $ϵ$ ($\kappa$), the Z (h) is essentially the SM Z boson (Higgs boson), while the dark photon $Z_d$ (dark Higgs $h_d$) is mostly X-like (S-like).

The hypercharge portal leads to a modification of the couplings of the neutral gauge bosons to fermions. In the case of the Z-boson couplings, the corrections start at $\mathcal{O}\left({ϵ}^2\right)$, while the dark photon vector couplings receive corrections already at $\mathcal{O}\left(ϵ\right)$. Explicitly, one has

$$g_{Z_{d}f\overline{f}}=ϵe{Q}_f+ϵe\left(Q_f-\frac{Y_f}{{cos}^2{\theta }_w}\right)\frac{m_{Z_d}^2}{m_Z^2},$$

(62)

where $e=\sqrt{4\pi \alpha }$ is the electromagnetic coupling and $Q_f$ ($Y_f$) is the electric charge (hypercharge) of the relevant fermion. Notice that for $m_{Z_d}\ll m_Z$, the $Z_d$-boson coupling to fermions is photon-like, while for $m_{Z_d}\simeq m_Z$, the $Z_d$ boson couples to fermions such as the SM Z boson. Referring to the $Z_d$ boson as dark photon is hence a bit of a misnomer, because a massive $Z_d$ boson always couples to neutrinos. We however follow this established naming convention. At $\mathcal{O}\left(ϵ\right)$, the hypercharge portal also leads to a coupling between the $125\mathrm{GeV}$ Higgs, a Z and a $Z_d$ boson. To this order, the coupling takes the form

$$g_{hZ{Z}_d}=\frac{2ϵtan{\theta }_w}{v}\frac{m_{Z_d}^2{m}_Z^2}{m_Z^2-m_{Z_d}^2}.$$

(63)

For $ϵ\gtrsim \mathcal{O}\left({10}^{-4}\right)$, a dark photon with mass $m_{Z_d}\gtrsim 1\mathrm{GeV}$ decays promptly [45,46]. At the LHC, such dark photons can be searched for in Drell–Yan (DY) dimuon production $pp\to Z_d\to {\mu }^{+}{\mu }^{-}$ and in four-lepton final states that arise from the process $pp\to h\to Z{Z}_d\to 4\ell$. The corresponding Feynman diagrams are shown on the left-hand side and in the middle of Figure 21.

The Higgs portal gives rise to a coupling between the $125\mathrm{GeV}$ Higgs and two dark photons. To lowest order in $\kappa$, this coupling can be written as

$$g_{h{Z}_d{Z}_d}=\frac{2\kappa v{m}_{Z_d}^2}{m_{h_d}^2-m_h^2}.$$

(64)

For $m_{Z_d}<m_h/2$, this coupling allows for Higgs decays of the form $h\to Z_d{Z}_d$—a representative Feynman diagram is shown on the right in Figure 21. At the same order in $\kappa$, interactions between the $125\mathrm{GeV}$ Higgs and two dark Higgses and a single dark Higgs and two dark photons also exist. In the following, we assume that the dark Higgs is heavy, i.e., $m_{h_d}\gg m_h/2$. In such a case, the exotic Higgs decay $h\to h_d{h}_d$ is kinematically forbidden, and the production cross-section for $pp\to h_d\to Z_d{Z}_d$ is always smaller than $pp\to h\to Z_d{Z}_d$. As a result, for sufficiently light dark photons, the best probe of the Higgs portal parameter $\kappa$ is the exotic Higgs decay $h\to Z_d{Z}_d$. The corresponding partial decay width is given to leading order in $\kappa$ by the following expression [45,46]

$$\Gamma \left(h\to Z_d{Z}_d\right)=\frac{{\kappa }^2{v}^2}{32\pi m_h}\frac{{(1-4x_{Z_d/h})}^{1/2}}{{\left(1-x_{h_d/h}\right)}^2}\left(1-4x_{Z_d/h}+12x_{Z_d/h}^2\right).$$

(65)

Notice that given the small total width of the SM Higgs boson, the branching ratio following from (65) can easily reach the few percent level for values of $\kappa \ll 1$. For example, taking $\kappa =0.07$, $m_{Z_d}=30\mathrm{GeV}$ and $m_{h_d}=300\mathrm{GeV}$ leads to $\mathrm{BR}(h\to Z_d{Z}_d)=0.15$, which is close to the model-independent limit on $\mathrm{BR}(h\to \mathrm{undet})$ reported in (60).

The best way to constrain the Higgs portal coupling now depends on the size of kinetic mixing. For $ϵ\gtrsim \mathcal{O}\left({10}^{-4}\right)$, the dark photon decays promptly, and one can again study four-lepton final states [290,291] to constrain the $h\to Z_d{Z}_d$ branching ratio and thereby $\kappa$. For $ϵ\lesssim \mathcal{O}\left({10}^{-4}\right)$, the decay length of the dark photon starts to become macroscopic. In fact, for $\mathcal{O}\left({10}^{-8}\right)\lesssim ϵ\lesssim \mathcal{O}\left({10}^{-4}\right)$ the decays of the $Z_d$ bosons are displaced with a large fraction of events ending up in the LHC detectors [45,46]. This leads to the exciting opportunity to probe very small values of $ϵ$ that are inaccessible by other means, provided some Higgs mixing is present in the dark photon model. Below we summarise the LHC searches that have considered the case of LLPs in dark photon models with both a hypercharge and a Higgs portal.

5.2.2. Experimental Constraints

The ATLAS, CMS and LHCb collaborations carried out searches for prompt dark photon decays targeting hypercharge mixing parameters $ϵ\gtrsim {10}^{-4}$ in both the DY channel $pp\to Z_d\to {\mu }^{+}{\mu }^{-}$ [292,293,294] and in four-lepton production associated with $h\to Z{Z}_d$ and $h\to Z_d{Z}_d$ [290,291]. ATLAS also searched for LLPs in the $h\to Z{Z}_d$ channel [295]. In addition, smaller $ϵ$ values were tested by dedicated searches for long-lived dark photons in final states with displaced dimuon vertices arising from $h\to Z_d{Z}_d$. Building on earlier LHC Run 1 analyses [296,297], the latest ATLAS and CMS results of this kind can be found in [298,299], respectively. In the following section, we review only LHC searches for LLP signatures involving dark photons.

From all possible dark photon decay modes, the $Z_d\to {\mu }^{+}{\mu }^{-}$ process is experimentally the most accessible one. First, the LHC experiments can trigger on muons with a high efficiency of up to 90%. Second, the rate of SM background processes is moderate, allowing for relatively low $p_T$ thresholds in the trigger. Third, the sophisticated MSs have a large acceptance for muon tracks with large impact parameters coming from displaced vertices. As a result, a wide range of proper lifetimes can be covered by the combination of the inner tracker employing LRT techniques and the MS.

ATLAS searched for dark photons in the $Z_d\to {\mu }^{+}{\mu }^{-}$ decay mode targeting proper decay lengths of ${10}^{-3}\mathrm{m}<c{\tau }_{Z_d}<{10}^3\mathrm{m}$ by requiring displaced vertices consistent with LLP decays to be reconstructed using the muon system [298]. A combination of single muon, dimuon and trimuon triggers with progressively lower thresholds down to $p_T>6\mathrm{GeV}$, as well as a $E_T^{\mathrm{miss}}$ trigger have been employed. The reconstruction efficiency for muon tracks from displaced vertices ranges from 70% for small impact parameters to 10% at the acceptance limit of $4\mathrm{m}$. Finally, the displaced vertices of oppositely charged muons with an invariant mass above $15\mathrm{GeV}$ are selected within the fiducial volume.

CMS also performed a search for dark photons in ${\mu }^{+}{\mu }^{-}$ final states [299]. In contrast to the aforementioned ATLAS search, this search analysed muon tracks reconstructed in the inner tracker, targeting small proper decay lengths of ${10}^{-3}\mathrm{m}<c{\tau }_{Z_d}<{10}^{-1}\mathrm{m}$. The novelty of this analysis is that it searches for pairs of oppositely charged muons with an invariant mass down to the dimuon production threshold. This is achieved using the so-called data scouting technique [294,300], where only partial event information, as reconstructed by the high-level trigger system, is recorded. This dramatically reduces the event size, thereby allowing to store orders of magnitude more candidate events for analysis, and lower the muon $p_T$ thresholds down to $4\mathrm{GeV}$. Finally, candidate muon pairs from the $Z_d\to {\mu }^{+}{\mu }^{-}$ decay are selected in events with two muons targeting DY production, and four muons with an invariant mass consistent with the $125\mathrm{GeV}$ Higgs boson aiming for $pp\to h\to Z{Z}_d$ and $pp\to h\to Z_d{Z}_d$. A sliding window fit is then applied to the distributions of dimuon invariant masses in several signal regions to look for a potential excess from $Z_d\to {\mu }^{+}{\mu }^{-}$ decays away from the known SM resonances such as the $J/\psi$ and ${\rm Y}$ mesons.

The region with $c{\tau }_{Z_d}<{10}^{-3}\mathrm{m}$ is covered by the CMS search for displaced $Z_d\to e^{+}{e}^{-},{\mu }^{+}{\mu }^{-}$ decays performed at $8\mathrm{TeV}$ [296], which targets the proper decay length range between ${10}^{-4}\mathrm{m}$ and ${10}^2\mathrm{m}$. This analysis is similar to the $13\mathrm{TeV}$ CMS search [294], except that it does not apply data scouting techniques and considers both the dielectron and the dimuon channel. Both search channels rely on the dilepton triggers.

In Figure 22, we summarise the existing 95% CL bounds on $\mathrm{BR}\left(h\to Z_d{Z}_d\right)$ as a function of the proper decay length $c{\tau }_{Z_d}$ of the dark photon. The shown ATLAS limits [298] correspond to the different dark photon masses $m_{Z_d}=20\mathrm{GeV},60\mathrm{GeV}$, while in the case of CMS, we display bounds for $m_{Z_d}=20\mathrm{GeV},50\mathrm{GeV}$ from [296] and for $m_{Z_d}=0.65\mathrm{GeV},2\mathrm{GeV},15\mathrm{GeV},30\mathrm{GeV}$, $50\mathrm{GeV}$ from [299]. Notice that for dark photons with masses of a few tens of GeV, the strongest bound on $\mathrm{BR}\left(h\to Z_d{Z}_d\right)$ typically arises for the lightest choice of the dark photon mass. For $m_{Z_d}\lesssim 5\mathrm{GeV}$ the limits on $\mathrm{BR}\left(h\to Z_d{Z}_d\right)$ however become significantly weaker with decreasing dark photon mass even if data scouting techniques are employed. Specifically, for dark photons with masses above $15\mathrm{GeV}$ the ATLAS (CMS) measurements exclude $\mathrm{BR}\left(h\to Z_d{Z}_d\right)\gtrsim {10}^{-4}$ (${10}^{-5}$) for proper decay lengths of around ${10}^{-1}\mathrm{m}$ (${10}^{-2}\mathrm{m}$), while for GeV-scale dark photons the corresponding CMS exclusions are weaker by about two orders of magnitude. One also observes that for proper decay lengths in the range of around ${10}^{-3}\mathrm{m}$ to ${10}^2\mathrm{m}$ the bound on $\mathrm{BR}\left(h\to Z_d{Z}_d\right)$ that derives from the LLP searches [298,299] is stronger than the model-independent limits quoted in (60). For shorter or longer proper decay lengths, the ATLAS Higgs coupling measurement [106] however provides currently the best constraint on the branching ratio of $h\to Z_d{Z}_d$.

The limits presented in Figure 22 can also be translated into 95% CL exclusions in the $m_{Z_d}$–$ϵ$ plane. This is achieved in Figure 23 where we show the constraints that derive from the ATLAS [298] and CMS [299] measurement assuming $\mathrm{BR}\left(h\to Z_d{Z}_d\right)={10}^{-3}$. For the same $h\to Z_d{Z}_d$ branching ratio, the exclusion that follows from [296] almost exactly resembles that of [299] for dark photon masses between $20\mathrm{GeV}$ and $50\mathrm{GeV}$. This limit is therefore not displayed in the figure. One observes that $ϵ$ values of order ${10}^{-8}$ and ${10}^{-7}$ are excluded by these searches. Notice that the CMS search leads to weaker limits in terms of $ϵ$ than the ATLAS search, because CMS does not consider decay lengths beyond ${10}^{-1}\mathrm{m}$. For comparison we also depict in Figure 23 the 90% CL bounds on the hypercharge mixing parameter that follow from the searches by LHCb and CMS in DY dimuon production [293,294] as well as the ATLAS search [290] that focuses on the $pp\to h\to Z{Z}_d\to 4\ell$ channel. These prompt limits are a few orders of magnitude weaker than the bounds that arise from the LLP searches, but they make no assumption about the amount of Higgs mixing in the dark photon model. In addition, the searches for $pp\to h\to Z_d{Z}_d\to 4\ell$ [290,291] also provide some sensitivity to $ϵ$, but the resulting upper limits are significantly weaker than those that stem from the DY searches and hence not reported in the figure.

An indirect constraint on the $m_{Z_d}$–$ϵ$ plane can also be derived from the available limits on the invisible Higgs branching ratio [73,106]. Imposing $\mathrm{BR}\left(h\to \mathrm{inv}\right)=0.1$, we obtain the dashed black contour shown in Figure 23 under the assumption that the dark photon is not registered by standard reconstruction algorithms for $c{\tau }_{Z_d}>0.1m(m_{Z_d}/\mathrm{GeV})$, resulting in a $E_T^{\mathrm{miss}}$ contribution. Another assumption is that the kinematic distributions such as the $E_T^{\mathrm{miss}}$ spectrum are the same for $pp\to h+X\to Z_d{Z}_d+X$ and for the SM Higgs production channels that go into the $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ bounds given in (3) and (60). While [277] suggests that these are good assumptions, we cannot quantify the associated systematic uncertainties. The shown indirect bound from invisible Higgs decays has therefore only an indicative character. However, it should be straightforward for ATLAS and CMS to directly reinterpret their Higgs to invisible searches in the context of long-lived dark photons, thereby improving on our naive estimate.

5.3. Models with a Vector and a Fermion Portal

Hidden sector models with both a vector and a fermion portal such as the Falkowski–Ruderman–Volansky–Zupan (FRVZ) model [50,51] also include a hidden photon as a possible LLP and can lead to signatures with displaced charged leptons. However, in models of this type, the hidden photons are not directly produced in the exotic decays of the $125\mathrm{GeV}$ Higgs boson but through a cascade involving SUSY and hidden sector particles with masses of $\mathcal{O}\left(10\mathrm{GeV}\right)$ or below. This implies that the hidden photons have to be even lighter than the other particles in the decay chain, with masses of $\mathcal{O}\left(1\mathrm{GeV}\right)$ in order to be kinematically accessible. Due to their small mass, the hidden photons are preferentially produced with large boosts at the LHC, resulting in so-called lepton-jets [47,48,49], i.e., collimated groups of leptons in a jet-like structure. These lepton-jet events are accompanied by varying amounts of $E_T^{\mathrm{miss}}$ depending on the precise pattern of the Higgs decay.

5.3.1. Theory

A minimal model that realises the general idea proposed in [47,48,49,50,51] contains a massive hidden photon ${\gamma }_d$ that communicates with the visible sector through mixing with the hypercharge field—see the first term in (61). In the FRVZ model, the particle content of the visible sector is that of the minimal supersymmetric SM (MSSM). Supersymmetrising the hypercharge portal leads to a mixing of the visible bino ($\tilde{B}$) and the hidden gaugino (${\tilde{\gamma }}_d$). Removing the kinetic mixing between the $\tilde{B}$ and the ${\tilde{\gamma }}_d$ then gives rise to interactions between all hidden fields charged under $U{\left(1\right)}_d$ and the visible neutralinos that are proportional to the hypercharge mixing parameter $ϵ$. In particular, one obtains an interaction term of the following form:

$$\mathcal{L}=-\frac{ϵg_d}{cos{\theta }_w}\tilde{B}\sum _i{q}_i{h}_d^{i†}{\tilde{h}}_d^i,$$

(66)

where $g_d$ is the $U{\left(1\right)}_d$ gauge coupling, $h_d^i$ are the hidden scalar fields, ${\tilde{h}}_d^i$ are the hidden gauginos and $q_i$ is the relevant $U{\left(1\right)}_d$ charge.

The interactions (66) lead to vertices involving a visible neutralino, a hidden neutralino and a hidden photon or a hidden Higgs boson. An exotic Higgs decay signal can therefore arise in the FRVZ model as follows: initially, the Higgs decays into a pair of the lightest visible neutralinos (${\tilde{N}}_1$), as indicated by the Feynman diagram on the left-hand side of Figure 24. In the pure MSSM without the hypercharge portal, the ${\tilde{N}}_1$ could be a DM candidate, i.e., the lightest SUSY particle or LSP, but in the FRVZ model, the presence of the term (66) allows the ${\tilde{N}}_1$ to decay into hidden sector states. Example diagrams are shown in the middle and on the right in Figure 24. In the first case, the ${\tilde{N}}_1$ decays to the lightest hidden sector gaugino (${\tilde{N}}_d$) and a hidden photon, while in the second case the final-state ${\tilde{N}}_d$ is accompanied by a hidden Higgs boson ($h_d$). The lightest hidden sector neutralino ${\tilde{N}}_d$ is stable and therefore a DM candidate, but the ${\gamma }_d$ and the $h_d$ decay further. If the hidden photon is the lightest hidden state, then ${\gamma }_d\to {\ell }^{+}{\ell }^{-}$ and $h_d\to {\gamma }_d{\gamma }_d\to 4\ell$ are possible decay chains that lead to the aforementioned lepton-jet signatures. Since both the hypercharge portal in (61) as well as the gaugino kinetic mixing (66) are proportional to $ϵ$, depending on the magnitude of this parameter, the lepton-jets can be either prompt or displaced.

Notice that models of inelastic DM (see for instance [301,302,303,304,305,306]) also possess many of the features discussed above. Such models typically contain two fermionic states ${\chi }_1$ and ${\chi }_2$ with a small mass splitting $\Delta m=m_{{\chi }_2}-m_{{\chi }_1}>0$. The role of ${\tilde{N}}_1$ (${\tilde{N}}_d$) is played by ${\chi }_2$ (${\chi }_1$) in inelastic DM models, and the DM candidate ${\chi }_1$ can be excited to its heavier twin ${\chi }_2$ by absorbing a massive dark photon ${\gamma }_d$. The simplest realisation of such a scenario consists in assuming a hypercharge portal and postulating a $U{\left(1\right)}_X$ symmetry that is spontaneously broken by a dark Higgs $h_d$. This dark Higgs couples to a pair of dark photons. The particle content and coupling structure of inelastic DM models therefore resemble quite closely those of the simplest FRVZ models. The difference between the two types of models is that, unlike the ${\gamma }_d{\tilde{N}}_1{\tilde{N}}_d$ coupling (66), the ${\gamma }_d{\chi }_1{\chi }_2$ coupling is not proportional to the hypercharge mixing parameter $ϵ$ but only involves the $U{\left(1\right)}_X$ coupling constant. The simplest inelastic DM models are hence in some sense a generalisation of the FRVZ idea. In fact, all the existing FRVZ interpretations of LHC searches correspond to such a generalisation, where the ${\gamma }_d{\chi }_1{\chi }_2$, ${\gamma }_{d}f\overline{f}$, $h_d{\chi }_1{\chi }_2$ and $h{\chi }_2{\chi }_2$ couplings are effectively treated as free parameters. As a concrete example, the ATLAS collaboration has used the following parameters in their interpretations

$$\begin{array}{cc}\hfill & g_{{\gamma }_d{\chi }_1{\chi }_2}=0.31,g_{h_d{\chi }_1{\chi }_2}=0.1,\hfill \\ \hfill & m_{{\chi }_1}=2\mathrm{GeV},m_{{\chi }_2}=5\mathrm{GeV},m_{h_d}=2\mathrm{GeV},\hfill \end{array}$$

(67)

where the choice of $g_{{\gamma }_d{\chi }_1{\chi }_2}$ corresponds to $g_{{\gamma }_d{\chi }_1{\chi }_2}=e=\sqrt{4\pi \alpha }$ with $\alpha$ the electromagnetic fine structure constant at the EW scale. The hypercharge mixing parameter $ϵ$ and the coupling $g_{h{\chi }_2{\chi }_2}$ are not directly used as external parameters but expressed through $c{\tau }_{{\gamma }_d}$ and $\mathrm{BR}\left(h\to {\chi }_2{\chi }_2\right)$, respectively, which then serve as input. Notice that for the parameter choices (67) and sufficiently large $\mathrm{BR}\left(h\to {\chi }_2{\chi }_2\right)$ values, the decays $h\to {\chi }_2{\chi }_2$ and ${\chi }_2\to {\chi }_1{\gamma }_d$ are necessarily prompt. From a theoretical point of view, the generalised FRVZ model is therefore quite similar to the dark photon model discussed in Section 5.2 if $\mathrm{BR}\left(h\to {\chi }_2{\chi }_2\right)$ is identified with $\mathrm{BR}\left(h\to Z_d{Z}_d\right)$.

5.3.2. Experimental Constraints

The ATLAS and CMS collaborations have searched for collimated groups of charged leptons or light hadrons in a jet-like structure to constrain exotic Higgs decays by exploring both prompt [307,308] and displaced [309,310] signatures. The results were interpreted in the generalised FRVZ framework described above, as well as in the context of other portal models. In the case of the generalised FRVZ model, the published LHC searches have focused on the two benchmark processes illustrated in Figure 25, with the Higgs boson produced in the ggF topology.

The searches for prompt Higgs decays targeted the $h\to {\gamma }_d{\gamma }_d+X$ topology (cf. left diagram in Figure 25), where the hidden photon is assumed to decay into charged leptons. The ATLAS search [307] included ${\gamma }_d$ decays to $e^{+}{e}^{-}$ and ${\mu }^{+}{\mu }^{-}$, while the CMS search [308] only considered dimuons. Hadronic decays of the hidden photon were not included, since they cannot be easily separated from the QCD multijet background. These searches probe signal hypotheses with hidden photon masses in the ranges of $0.1\mathrm{GeV}$ to $2\mathrm{GeV}$ for ATLAS and $0.25\mathrm{GeV}$ and $2\mathrm{GeV}$ for CMS. In the ATLAS analysis, events are required to have at least two lepton-jets, which are reconstructed by clustering tracks from the primary vertex within a radius $\Delta R=\sqrt{{(\Delta \eta )}^2+{(\Delta \varphi )}^2}=0.5$ of the highest $p_T$ track, and subsequently matching them to electron or muon candidates. In the CMS analysis, events are required to have at least two dimuon pairs, which are reconstructed by combining muon-candidate tracks into a common vertex.

In the domain of searches targeting LLP signatures, the ATLAS analysis [309] has probed the $h\to 4{\gamma }_d+X$ topology (cf. right diagram in Figure 25) considering both leptonic and hadronic decays of the hidden photons. The mass of the hidden photon was set to $0.4\mathrm{GeV}$, resulting in a 10% branching ratio into pions, and the rest of the branching ratio divided equally between electrons and muons [50]. The dark photon decays into muon pairs were targeted by reconstructing jets of muons that were reconstructed using the MS only, while the decays into electron pairs and hadrons were reconstructed as calorimeter jets with a high fraction of energy deposited in the hadronic calorimeters. In both channels, the jets were required to be isolated from any ID activity, as expected for LLPs, and multivariate analysis techniques using timing and topological information of the jet were employed to discriminate signal from background. The CMS search [310] explored the $h\to {\gamma }_d{\gamma }_d+X$ topology (cf. left diagram in Figure 25), considering only ${\gamma }_d\to {\mu }^{+}{\mu }^{-}$ decays. Hidden photon masses in the range $0.25\mathrm{GeV}$ and $8.5\mathrm{GeV}$ were explored, with the upper limit set by the requirement $m_{\mu \mu }<9\mathrm{GeV}$ in order to sufficiently suppress the background from DY and ${\rm Y}\to {\mu }^{+}{\mu }^{-}$ production. Proper decay lengths below $c{\tau }_{{\gamma }_d}<0.1\mathrm{m}$ were targeted, explicitly including both prompt and displaced signatures. Candidate events were required to have exactly two ${\gamma }_d\to {\mu }^{+}{\mu }^{-}$ candidates that are isolated from significant activity in the tracking system and have an invariant mass consistent with each other.

Figure 26 contains a summary of 90% CL exclusions in the $m_{{\gamma }_d}$–$ϵ$ plane in the generalised FRVZ model for hidden photon masses in the range of $0.1\mathrm{GeV}$ to $10\mathrm{GeV}$. The constraints following from the ATLAS [309] and CMS [310] search for the exotic Higgs decay $h\to {\gamma }_d{\gamma }_d+X$ apply in the generalised FRVZ model and assume $\mathrm{BR}\left(h\to {\gamma }_d{\gamma }_d+X\right)=0.1$. The ATLAS analysis employs (67) and focuses on the mass range $0.2\mathrm{GeV}\lesssim m_{{\gamma }_d}\lesssim 3.6\mathrm{GeV}$ and small kinetic mixings that lead to LLP signatures. Depending on $m_{{\gamma }_d}$, it excludes hypercharge mixing parameters within $1·{10}^{-5}\lesssim ϵ\lesssim 3·{10}^{-7}$. The CMS search probes hidden photon masses in the range $0.25\mathrm{GeV}\le m_{{\gamma }_d}\le 8.5\mathrm{GeV}$ and targets proper decay lengths of $c{\tau }_{{\gamma }_d}\le 0.1\mathrm{m}$. This analysis is able to exclude hypercharge mixing parameters $ϵ\lesssim 3·{10}^{-6}$ ($ϵ\lesssim 7·{10}^{-8}$) at low (high) $m_{{\gamma }_d}$. The displayed LHCb [293], BaBar [311], KLOE-2 [312] and $\nu$-CAL I [313] bounds have been taken from the Darkcast package developed in [314]. They hold irrespectively of the value of the $h\to {\gamma }_d{\gamma }_d+X$ branching ratio. The exclusion that follows from the limit (60) is obtained from our recast as described in Section 5.2 and assumes $\mathrm{BR}\left(h\to \mathrm{inv}\right)=0.1$. The limit arising from the recent ATLAS analysis [277] of mono-jet signatures is not shown in the figure, because it is not sensitive to $\mathrm{BR}\left(h\to {\gamma }_d{\gamma }_d+X\right)=0.1$. Figure 26 clearly shows that in the context of the generalised FRVZ model, LHC searches for exotic Higgs decays provide an opportunity to probe values of $ϵ$ that are at present inaccessible by other means. Future LLP experiments such as MATHUSLA, CODEX-b, ANUBIS, FASER [315] and FASER2 [316] located at the LHC, in combination with Belle II [305,306], are expected to set additional stringent constraints on the $m_{{\gamma }_d}$–$ϵ$ plane in the generalised FRVZ or inelastic DM frameworks.

6. Outlook

In this article, we have reviewed the status of the LHC constraints on DM scenarios where the $125\mathrm{GeV}$ Higgs boson plays a prominent role. Specifically, we have covered the case of SM Higgs portals (Section 2), models with extended Higgs and gauge sectors with one mediator (Section 3) or two mediators (Section 4), as well as theories that can lead to LLP signatures induced by the decay of the $125\mathrm{GeV}$ Higgs (Section 5). In all cases, we have collected the existing constraints from the different LHC experiments to obtain state-of-the-art summary plots that show the current LHC sensitivities to the relevant model parameters in various benchmark scenarios. In some cases, we have also indicated in these summary plots the restrictions that other measurements impose, including those of the properties of the $125\mathrm{GeV}$ Higgs boson, flavour physics, EW precision measurements, DM DD and ID experiments, and the relic density, to just name a few.

While our review mostly focuses on “known knowns”, i.e., existing experimental constraints and their interpretations, we have also tried to discuss aspects of the considered models that are in principle known but that might be unknown to at least some of the readers, i.e., “unknown knowns”. For instance, in the case of the Higgs portal models, we have emphasised in Section 2.4 that $E_T^{\mathrm{miss}}$ searches for off-shell Higgs production in the VBF, $t\overline{t}$ and $tW$, and potentially other channels provide sensitivity to the parameter space where $m_{\mathrm{DM}}>m_h/2$. In fact, in the case of the kinetic Higgs portal, such searches are the only known way to test DM masses above the Higgs threshold, making dedicated experimental searches and/or interpretations of the relevant mono-X signatures in our opinion an important goal for future LHC runs. In the context of the Higgs portal models, we have furthermore argued in Section 2.5 that interpreting LHC limits on the invisible Higgs branching ratio in terms of a SI DM–nucleon cross-section leads to somewhat unintuitive results, in the sense that the indirectly obtained ${\sigma }_{\mathrm{SI}}$ bounds do not become constant in the limit $m_{\mathrm{DM}}\ll m_h/2$, as one would naively expect. We have explained the reason for this unexpected feature. To avoid potential confusion, we suggest to present the LHC limits on $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ as well as the limits on ${\sigma }_{\mathrm{SI}}$ obtained by DD experiments in terms of the effective interaction strength of the different Higgs portal models, for example, in terms of $|c_m|$ in the case of the marginal Higgs portal. This is the standard presentation in the theoretical literature on EFT Higgs portals, and adding it to the interpretation of results from the LHC and/or DM DD experiments would further promote a fruitful exchange between these communities as well as the theory community.

We have stressed at the end of Section 3.1 that all existing collider studies of the 2HDM+a model have assumed a Yukawa sector of type-II. For this choice, the bounds from FCNC processes and LHC searches for heavy Higgs bosons are strong, pushing the masses of the additional 2HDM Higgs bosons above the $500\mathrm{GeV}$ range. In fermiophobic 2HDM models of type-I, the constraints on the additional Higgs bosons can be significantly relaxed, thereby allowing for new scalars and pseudoscalars with EW-scale masses. The fermiophobic nature of the Higgs bosons can lead to unconventional production mechanisms and decay patterns of the non-SM spin-0 particles. In our opinion, the mono-X phenomenology in a fermiophobic 2HDM+a deserves dedicated studies. The Yukawa sector of the model also plays an important role for the 2HDM+s model considered in Section 3.2. In particular, we have shown that a suitable choice of the Yukawa sector allows one to tame excessive 2HDM+s corrections to the SI DM–nucleon cross-section, opening up the possibility to test model realisations at the LHC which lead to a viable DM phenomenology that is consistent with the observed relic density.

The existing LHC interpretations of the 2HDM-$Z^{\prime }$ model and the 2MDM model discussed in Section 4.1 and Section 4.2, respectively, have all focused on the mono-Higgs channel with $h\to b\overline{b},\gamma \gamma$ in the former case and on the mono-s channel with $s\to b\overline{b},VV$ in the latter case. In order to present a more global picture of the LHC sensitivity to these models, we have reinterpreted existing LHC searches for heavy spin-1 resonances decaying to visible particles. Our studies show that mono-X searches are not the only way to probe these two-mediator models at the LHC. In fact, both discussed DM models are in general more tightly constrained by resonance searches for the $Z^{\prime }$ boson in dijet or $t\overline{t}$ final states than by mono-X searches. We believe this to be a rather generic feature of DM models with a heavy spin-1 mediator. Despite the dominant sensitivity of resonance searches in the context of the 2HDM-$Z^{\prime }$ and 2MDM models, the searches for $X+E_T^{\mathrm{miss}}$ signals are important given their distinct experimental signature. We have furthermore pointed out that the 2HDM-$Z^{\prime }$ model can also be probed by a search for $Z^{\prime }\to Zh$ resonances. Similarly, we advocate the exploration of the $hh+E_T^{\mathrm{miss}}$ and $tt+E_T^{\mathrm{miss}}$ topologies in the context of the 2MDM model. We are convinced that only through a combined exploration of the whole suite of searches in various channels can the full LHC potential be exploited.

Searches for BSM LLPs have gained a lot of momentum in LHC Run 2. In Section 5.1–Section 5.3, we have presented state-of-the-art summary plots for models of neutral naturalness/hidden valleys, dark photon models and BSM theories with both a hypercharge and a fermion portal, respectively, that can lead to LLP signatures in the decay of the $125\mathrm{GeV}$ Higgs boson. We have stressed that an indirect way to constrain models of this type is provided by the precision measurements of the $125\mathrm{GeV}$ Higgs properties in visible and invisible final states. For instance, in the case of neutral naturalness, we find that the limits from $\mathrm{BR}\left(h\to \mathrm{undet}\right)$ and $\mathrm{BR}\left(h\to \mathrm{inv}\right)$ provide currently the strongest LHC bound on $\mathrm{BR}\left(h\to aa\right)$ for $c{\tau }_a\lesssim {10}^{-4}\mathrm{m}$ and $c{\tau }_a\gtrsim {10}^3\mathrm{m}$, respectively. To further emphasise the possibility to test LLP scenarios via $E_T^{\mathrm{miss}}$ searches, we have in the case of the two models with a hypercharge portal performed recasts of the limits imposed by the latest $h\to \mathrm{inv}$ results to derive upper limits on the hypercharge mixing parameter $ϵ$ as a function of the dark photon mass.

The “known knowns” and the “unknown knowns” discussed in this review are clearly only a snapshot of the broad landscape of collider searches for DM through the Higgs lens. We tried our best to represent a complete picture of experimental searches and the relevant theoretical works and apologise for any potential omission. We hope that our review can trigger the experimental exploration of new search strategies and possibly even groundbreaking new ideas or discoveries in the years leading to the HL-LHC era.