Collider Searches for Dark Matter through the Higgs Lens

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 \, {\rm 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.


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 non-relativistic.
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 miss T ) plus X signatures, also known as mono-X, where the E miss T 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 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 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 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 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 emphasise is thereby put on highlighting the complementary of the different non-collider 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 besides the DM candidate four additional BSM spin-0 states. 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 miss T 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 similarities and difference 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 models, which typically have a rich collider and DM phenomenology. In Section 4 we discuss two representative models in more detail: the 2HDM+Z 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 and the 2MDM model contain a spin-1 mediator that couples the DM particles to some of the SM states. This feature allows to the test the models by searching for resonant production of non-E miss T final states such as dijets, tt and Zh. Using the latest available LHC data, we derive the constraints on the 2HDM+Z 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 miss T searches exclude additional parameter space not probed by the existing mono-Higgs interpretations. (iv) Models that feature exotic decays of the 125 GeV Higgs boson into hidden sector particles with macroscopic proper decay lengths of cτ > 10 −2 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 Neutral naturalness Section 5.1, Figure 20 Dark photons Section 5.2, Figures 22 and 23 Vector plus fermion portal Section 5.3, Figure 26 VBF + E miss T EFT Higgs portals Section 2.2, Figure 3 tt + E miss      Vector plus fermion portal Section 5.3, Figure 26 h → inv, 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 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 miss T , (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.
in these plots also other experimental limits 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 fb −1 collected at 13 TeV. For the readers mostly interested in collider phenomenology, Table 1 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.

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 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 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 Z 2 exchange symmetry.

Theory
In this section, we will consider the simplest possibility of these Higgs portal models, namely the case of a real scalar φ that is a singlet under the SM gauge group, but odd under a Z 2 symmetry, i.e. φ → −φ. This guarantees the stability of φ making it a suitable DM candidate. The interactions between the dark sector and the SM that we consider are 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 Λ 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 Λ 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 non-collider 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 L φH . After drawing our general conclusions we will 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][74].

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 φ < m h /2 in the case of the real scalar φ with mass m φ and m h 125 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 where v 246 GeV denotes the vacuum expectation value (VEV) of the 125 GeV Higgs boson. The formula (2) can be used to translate experimental limits on the Higgs to invisible branching ratio BR(h → inv) into constraints on the strength of the marginal and derivative Higgs portals. In fact, in the limit m φ m h the best existing 95% confidence level (CL) exclusion LHC bound [74] of when the SM value Γ SM h 4.07 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 BR(h → inv) < 2.5 · 10 −2 [75]. This implies that the bounds (4) may be improved to 2.3 · 10 −3 and 2.6 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 GeV Higgs boson, E miss T 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 miss T channel as well as tt + E miss in Figure 2. The LHC reach of these channels in the context of (1) has been studied recently in [25,68], and we will summarise the main findings of these articles below.

DM phenomenology
DM states that couple to the 125 GeV Higgs typically lead to spin-independent (SI) DM-nucleon cross sections (σ SI ), which are severely constrained by the existing DM DD experiments like XENON1T. While the marginal Higgs portal leads to an unsuppressed SI DMnucleon cross section, the DM-nucleon interactions that are mediated by the derivative Higgs portal are suppressed by q 2 /Λ 2 (100 MeV) 2 /Λ 2 where q 2 characterises the momentum transfer in DM scattering with heavy nuclei. Explicitly one finds where m N 939 MeV is the average of the nucleon mass and f N 0.31 [76][77][78][79] parameterises the strength of the Higgs-nucleon interactions. For m φ = 100 GeV the latest XENON1T 90% CL upper limit on the SI DM-nucleon cross section reads σ SI < 9.12 · 10 −47 cm 2 [80]. By means of (5) this bound can be translated into a limit on the marginal Higgs portal coupling: In 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 dimensionsix the derivative Higgs portal is the only spin-0 DM-Higgs operator that does naturally satisfy the constraints imposed by σ 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 where T denotes the DM temperature. Notice that in today's Universe T 0 0, while at freeze-out T f m φ /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 (Ω DM h 2 ), in particular, if the s-wave coefficient a X is parametrically suppressed.
For m b < m φ m W with m b 4.2 GeV (m W 80.4 GeV) the bottom-quark (W-boson) mass, DM annihilation into bottom-antibottom quark pairs is the dominant contribution to Ω DM h 2 . The corresponding s-wave coefficient reads with Γ h denoting the total decay width of the 125 GeV Higgs boson including contributions from h → φφ see (2) . In the case m φ m W the φφ → W + W − , ZZ, hh, tt 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 with X = W + W − , ZZ, hh and N W + W − = 2, N ZZ = N hh = 1. Notice that in the case of m φ v, DM annihilation to W and Z bosons reduces to three times the contribution from annihilation to the 125 GeV Higgs boson. This is an expected feature in the SU(2) L × U(1) 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 where the sum over X involves all annihilation channels that are kinematically open at a given value of m φ . While the above formulas represent useful expressions to estimate Ω DM h 2 , we will use micrOMEGAs [81] in our numerical analysis to obtain the constraints on the parameter space of (1) that follow from the PLANCK measurement Ω DM h 2 = 0.120 ± 0.001 [82]. The ID exclusions shown below in Figure 3 are also determined with the help of micrOMEGAs.

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 BR(h → inv) as given in (3) while the dashed black lines represent the expected HL-LHC 95% CL limit BR(h → inv) < 2.5 · 10 −2 [75]. The purple region in the upper plot is disfavoured by the 90% CL bounds of XENON1T [80] on σ SI . The vertical orange shaded bands indicate the DM mass ranges that are excluded at 95% CL by the γ-ray observations of dwarf spheroidal galaxies (dSphs) of the Fermi-LAT and DES collaborations reported in [84]. The used experimental bounds assume DM annihilation via φφ → bb and that Ω DM h 2 = 0.12. Compared to φφ → bb, the constraints that follow from the latest Fermi-LAT search for monochromatic photons [85] 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 Ω DM h 2 > 0.12, i.e. larger values of the DM relic density compared to the PLANCK measurement [82]. The green regions correspond to the 95% CL exclusion limits found in [25] from a study of off-shell invisible Higgs production in the VBF + E miss T channel. Finally, the blue domains represent the 95% CL constraints obtained by the combined tt + E miss T and tW + E miss T (tX + E miss T ) analysis strategy discussed in [68].  [75] are indicated by dashed black lines. The purple region in the upper panel is disfavoured by the 90% CL bound on the SI DM-nucleon cross section σ SI set by XENON1T [80]. The vertical orange shaded bands display the range of DM masses that is excluded at 95% CL by Fermi-LAT and DES [84]. The red curves correspond to the value Ω DM h 2 = 0.12 [82]. In the parameter space below (above) the red curve the Universe is overclosed in the case of the marginal (derivative) Higgs portal model. The green regions indicate the 95% CL exclusion limit obtained in [25] from a study of off-shell Higgs production in the VBF + E miss T channel, while the blue regions represent the corresponding exclusion limits derived in [68] from a study of tX + E miss T final states. For further details see the main text.
The latter two types of collider limits assume an integrated luminosity of 3 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 φ 5 GeV, while in the range 5 GeV m φ < m h /2 they are roughly comparable in strength. In the case m φ > m h /2, the bounds that follow from σ 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 φ > m h /2. Notice however that the study [25] assumes a systematic uncertainty of 1%, while the shown tX + E miss T exclusion limits are based on a systematic uncertainty of 15% [68]. Assuming a reduction of background uncertainties in tX + E miss T down to 5% would bring the VBF + E miss T and tX + E miss T exclusion limits closer together. Combining the two mono-X channels as done in the case of the LHC searches for the invisible Higgs boson decays (cf. for instance [69][70][71][72][73][74]) 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 VBF + E miss T off-shell Higgs production has been studied recently in [25]. For earlier analyses see also [86][87][88][89][90].
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 φ < m h /2, HL-LHC measurements of invisible Higgs decays exclude large parts of the parameter space that lead to Ω DM h 2 = 0.12. Only a narrow corridor around the 125 GeV Higgs-boson resonance survives the DM relic density constraint, which is however excluded by DM ID measurements. Given that σ 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 φ > m h /2 using mono-X searches at the HL-LHC, however only if Ω 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 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, 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 VBF + E miss T or tX + E miss T . 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 present.

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 [80,[91][92][93] on the SI DM-nucleon cross section and the exclusion limits that derive from the latest ATLAS measurement of Higgs to invisible decays [74]. At 90% CL this measurement leads to BR(h → inv) < 0.09 -the corresponding 95% CL bound is given in (3). The translation of the BR(h → inv) bound into a limit on σ SI relies on an EFT approach under the assumption that the 125 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)  bounds on σ SI . Notice that the limit on the Wilson coefficient c m does only marginal change when going from DM masses of 10 GeV down to 1 GeV (see the upper panel in Figure 3) while the bound on the SI DM-nucleon cross section worsens notable. While this is puzzeling at first, one has to realise that in the scalar case σ SI scales with 1/(m φ + m N ) 2 cf. (5) which implies that for constant c m and the DM mass m φ sufficiently larger than m N , one has σ SI ∝ 1/m 2 φ . This feature leads to a deterioration of the limit on σ 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 χ as shown in Figure 4 we first provide expressions for the relevant effective interactions, the partial Higgs decay width for h → χχ and the SI DM-nucleon cross section. The Higgs portal takes the following form where c χ is a Wilson coefficient assumed to be real and the scale Λ suppresses the interaction because the operator is of dimension five. The corresponding partial Higgs decay width reads while the SI DM-nucleon cross section is given by We add that for a complex c χ the SI DM-nucleon cross section can be strongly suppressed (cf. for instance [94]). From (12) it follows that for m χ → 0 measurements of BR(h → inv) simply lead to constant bound on the combination |c χ |/Λ of parameters. The SI DM-nucleon cross section (13) however scales as m 2 χ /(m χ + m N ) 2 where compared to (5) the additional factor of m 2 χ appears due to dimensional reasons. It follows that decreasing the DM mass will lead to a steady improvement of the limit on σ SI from BR(h → inv) (see Figure 4) even though |c χ |/Λ stays essentially constant.
The above discussion suggests that from the point of view of collider physics interpreting the limits on BR(h → inv) in terms of exclusions on σ SI is not the optimal choice. A better way to interpret and to compare the BR(h → inv) 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 |, Λ/|c d |, Λ/|c χ |, etc.). In the case of the marginal and derivative Higgs portal this has been done 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.

Portals with extended Higgs sectors
In the spin-0 simplified models with a fermionic singlet DM candidate χ, 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(2) L × U(1) Y gauge symmetry. This feature leads to unitarity violation at high energies [95,96] -for discussions of unitarity violation in the case of spin-1 simplified DM models see [97][98][99]. 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 GeV Higgs [100][101][102][103][104][105][106]. Compared to the spin-0 simplified DM models, this model leads to additional E miss T signatures such as Higgs to invisible decays, mono-W, mono-Z and VBF + E miss T production at tree-level -see [107] 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. [108,109] for the latest global Higgs analyses by ATLAS and CMS) lead to a suppression of all E miss T signals. This suppression still allows the LHC to test the model via (3) if m χ < m h /2. However, in the case m χ > 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 it extended. Two-Higgs-doublet models (2HDMs) [110,111] 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 where we have imposed a Z 2 symmetry under which H 1 → H 1 and H 2 → −H 2 but allowed for this discrete symmetry to be broken softly by µ 3 H † 1 H 2 + h.c. The Z 2 symmetry is the minimal condition necessary to guarantee the absence of flavour-changing neutral currents (FCNCs) at tree level [112,113] and such a symmetry is realised in many well-motivated complete UV theories in the form of supersymmetry (SUSY), a U(1) 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 as 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 ± . Besides the masses of the five Higgs bosons, the 2HDM parameter space involves the angles α and β. 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 In the context of the inert doublet model [114][115][116] the scalar potential (14) alone already allows for interesting DM and collider phenomenology [117][118][119][120][121][122][123][124][125][126][127][128][129][130][131]. 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 (like in the case of the simplified DM models [36,37]) considering only (14) is therefore not enough and one generically needs besides the DM particle an additional mediator. 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][132][133][134][135][136][137][138][139][140][141][142][143][144]. Below we review in detail two of the models discussed in these articles.

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 χ, 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 miss T 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 [138,[145][146][147][148][149][150][151][152][153][154].

Theory
If the Dirac DM field χ and the pseudosalar P are taken to transform under the Z 2 symmetry as χ → −χ and P → P, the only renormalisable DM-mediator coupling that is allowed by symmetry is L χ = −iy χ Pχγ 5 χ .
Here it is assumed that the dark-sector Yukawa coupling y χ is real in order not to violate CP. Besides (14) the scalar potential in the 2HDM+a model contains the following terms that connect singlets to doublets. The parameters b P , λ P1 and λ P2 are taken to be real to not violate CP. Notice that the first term in V HP breaks the Z 2 symmetry softly. The singlet potential reads 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 miss T 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: Here sin θ 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 θ 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 α between the CP-even scalars h and H is for instance constrained by Higgs coupling strength measurements [108,109]. For arbitrary values of tan β only parameter choices with cos(β − α) 0 are experimentally allowed. In order to satisfy the constraints from Higgs physics, the existing experimental 2HDM+a analyses [145][146][147][148][149][150][151][152][153][154] have concentrated on the so-called alignment limit of the 2HDM where cos(β − α) = 0 [155], treating tan β 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 ± 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 has been shown that the sum of the potentials (14) and (16) has a custodial symmetry if cos(β − α) = 0 and m A = m H = m H ± . For such parameter choices the EWPOs are satisfied for any value of sin θ and m a , which renders the choice m A = m H = m H ± 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 λ P1 , λ P2 > 0 and m A = m H = m H ± one can show [138] that there are two bounded from below conditions: These inequalities suggest that λ 3 has to be sufficiently large, in particular, if the mass splitting of the two pseudoscalars and/or sin θ are large. Since the relations (19) are modified in models with more general scalar potentials including dimension-six operators [135,156] the bounded from below requirements are not directly imposed in the ATLAS and CMS interpretations of the E miss T searches [145][146][147][148][149][150][151][152][153][154]. Instead the choice λ 3 = 3 is employed which generically allows for heavy Higgses above 1 TeV, while keeping λ 3 small enough to stay in the perturbative regime.
The quartic couplings λ 3 , λ P1 and λ 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 ± the Aah and Haa couplings take the following form [33] These expressions imply that parameter choices of the form λ 3 = λ P1 = λ P2 are well suited to keep the total widths Γ A and Γ 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 ± arise from Z → bb, B → X s γ, B s → µ + µ − and B-meson mixing -see [33,138,144] 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 → X s γ decay sets a 95% CL lower limit m H ± > 800 GeV [157,158] that for tan β 2 is practically independent of the specific choice of tan β. In other 2HDM realisations such as a fermiophobic 2HDM model of type-I (see [159,160] for recent detailed discussions) all flavour bounds can however be significantly relaxed, allowing for EW-scale values of m H ± . 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 [145][146][147][148][149][150][151][152][153][154] 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,138,143]. The most interesting signatures arise if one deviates from the alignment limit and/or gives up the assumption m A = m H = m H ± . In particular searches for final states involving tops and/or W bosons [134,139,143,147,156] 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 will give an example of a non-E miss T that already now provides relevant constraints on the 2HDM+a parameter space.

Anatomy of E miss
T signatures The 2HDM+a model provides a wide variety of E miss T signatures that can be searched for at the LHC. As pointed out in the works [32,33,134] the most interesting signals are h + E miss T (i.e. mono-Higgs), Z + E miss T (i.e. mono-Z) and tW + E miss T production, because these signals can be resonantly enhanced. Examples of representative diagrams that lead to the discussed E miss T signals in ggF production are displayed in Figure 5. Diagrams with bottomquark 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 miss T production. Similar resonance enhancements arise from the graphs on the left-hand side for the mono-Z (middle row) and tW + E miss T channel (lower row) if m H > m Z + m a and m H ± > m W + m a , respectively. Notice that resonant mono-Higgs, mono-Z and tW + E miss T 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 ± not only enhances the mono-Higgs, mono-Z and tW + E miss T compared to the simplified pseudoscalar DM model, but also leads to a quite different kinematics in these channels [32,33,134,138]. Besides ggF production, important contributions to the mono-Higgs and mono-Z channels can arise from bb-fusion production in 2HDM realisations such as type-II models in which the bottom-quark Yukawa coupling is enhanced by tan β. 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 non-resonant contributions (right column) in the case of bb-fusion. An additional E miss T signature that allows to put constraints on the parameter space of the 2HDM+a model are invisible Higgs decays [33]. The relevant decay channel is h → aa followed by a → χχ. The decay of the 125 GeV Higgs to a pair of pseudoscalars a is proportional to the square of the haa coupling. For cos(β − α) = 0 and m A = m H = m H ± this coupling takes the form Ignoring fine-tuned solutions for which |g haa | 1, one can show that for sufficiently light DM masses the bound (3) leads to a lower limit of m a 100 GeV. We add that the direct measurements of the total Higgs width [161,162] that impose Γ h < 1.1 GeV at 95% CL furthermore imply that m a m h /2 62.5 GeV. This bound also holds for heavy DM, because the processes h → aa with a → ff , where f denotes all the kinematically accessible SM fermions, always provides a sizeable non-SM contribution to Γ h unless (21) is artificially small.
In addition to the E miss T signatures discussed so far, the 2HDM+a model also leads to mono-jet, tt + E miss T and bb + E miss T signals. In contrast to the case of the simplified pseudoscalar DM model where these mono-X channels provide the leading constraints (see [96,145,[163][164][165][166][167][168][169][170] for the latest theoretical and experimental analyses) this is not the case in the 2HDM+a, because the mono-jet, tt + E miss resonantly enhanced contributions. Reinterpreting the existing spin-0 simplified DM model searches in the 2HDM+a context is relatively straightforward [33,138].

Summary plots
In the articles [138,[145][146][147][148][149][150][151][152][153][154] various 2HMD+a benchmark scans have been recommended and studied. The following parameters are common to all existing benchmark scans and we have motivated these choices in Section 3.1.1. Unless the relevant parameter is scanned over the plots shown below also employ where the latter choice ensures a sizeable BR(a → χχ) for the pseudoscalar masses of interest, i.e. m a 100 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 where y t = √ 2m t /v 0.94 and y b = √ 2m b /v 0.02 are the top-quark and bottom-quark Yukawa coupling, respectively. These expressions imply that for the benchmark choice (23) of tan β 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 β the bb-induced graphs of Figure 6 can give a relevant or even the dominant contribution.    (22) and (23). The coloured exclusions arise from [149][150][151][152][153]171]. The hatched grey regions correspond to the parameter space where any of the additional Higgs bosons has a relative width of more than 20%. This is indicated by the label Γ/m > 20%. The dashed grey line corresponds to the equality m A = m a + m h . See text for further details.
The mono-X searches that provide at present the bulk of the sensitivity to the 2HDM+a [149,153] and tW + E miss T production [150]. Thus we will focus on these searches below, disregarding other mono-X searches such as mono-jets or tt + E miss T . To evade the limits from invisible Higgs decays (cf. Section 3.1.2), we consider only m a values satisfying m a > 100 GeV when studying the 2HDM+a parameter space. In addition, we will show that the recent search for charged Higgs bosons in the channel H ± → tb [171] allows to set stringent constraints on m H ± , which under the assumption of degenerated 2HDM spin-0 masses (22), translate into bounds on m A or m H . As we will see, the obtained limits are to first approximation independent of m a , and as a result the H ± → 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 miss T searches such as for instance four-top production [138,147] are not considered. 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 θ = 0.7, which corresponds to close-to-maximal mixing in the pseudoscalar sector, the Z ( + − ) + E miss understood by noticing that the couplings relevant for resonant production scale with sin θ. In the alignment limit with m A = m H = m H ± one has explicitly:

Scans in the m
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 miss T 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 → A → ha → h + χχ and gg → a → ha → h + χχ amplitudes is proportional to Here we have employed the parameter choices (22). For fixed values of m a and sin θ, 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 non-resonant contributions -see the right Feynman diagram in the upper row of Figure 5.  (22) and (23). In the plot on the right-hand side the parameter space where any of the non-SM Higgs bosons has a relative width of more than 20% is indicated by the hatched grey region and the label Γ/m > 20%.
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 θ -cf. (20), (21) and (25). While off-shell effects are taken into account in Figure 7, possible modifications of the line shape of the intermediate Higgs bosons [172][173][174] are ignored. The latter effects have been studied in [175,176], 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 has 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 ± → tb are to a very good approximation independent of m a , and we find that charged Higgs boson masses m H ± 700 GeV (m H ± 600 GeV) are strongly disfavoured for sin θ = 0.35 (sin θ = 0.7). Searches for H ± → 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 → X s γ (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 β
The 95% CL limits that the searches [149][150][151][152][153]171] set in the m a -tan β and m A -tan β planes are displayed in Figure 8. For tan β 5, the bb-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 β in the case of a 2HDM model of type-II. In fact, the sensitivity of the h + E miss T searches at high tan β solely stems from bb-fusion production, while at low tan β the sensitivity comes entirely from ggF production. Production via bb-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  (22) and (23) are m A = 600 GeV, m a = 250 GeV and sin θ = 0.35. The solid lines correspond to the limits arising from the mono-X searches [151][152][153], while the dashed blue line corresponds to the calculated relic density for the studied 2HDM+a benchmark model. this region. In particular, the h(bb) + E miss T exclusion at high tan β and low sin θ has been shown to come almost exclusively from a selection with additional b-jets [152]. Notice that the H ± → tb search is sensitivity to both low and high values of tan β, since the associated coupling involves both y t cot β and y b tan β tems. Parameter regions with tan β 0.3 are incompatible with the requirement of having a perturbative top-quark Yukawa coupling [111] and therefore this region is not displayed in the exclusion plots, while flavour observables disfavour very high values of tan β -see for instance [33,138,144].

Scans in sin θ
In Figure 9 we show the 95% CL exclusion limits as a function of sin θ that follow from [151][152][153] for two sets of pseudoscalar masses m A and m a . The sensitivity of all searches improves as sin θ increases from zero, since for sin θ = 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 θ values, the E miss T spectrum becomes harder, leading to an increased signal acceptance. The improvement of the sensitivity of the h(γγ) + E miss T search at high m A and sin θ 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 χ of the DM candidate. Figure 10 illustrates these constraints at 95% CL for (22), (23), m A = 600 GeV, m a = 250 GeV and sin θ = 0.35. In the region m χ < m a /2, where the a → χχ decay is kinematically allowed, a variation of m χ has a negligible impact on the production cross sections and the E miss T spectrum, and in consequence the sensitivity of the X + E miss T searches is independent of m χ . For higher masses, the production cross sections decrease significantly and the E miss T becomes softer, significantly worsening the sensitivity. For this particular benchmark the existing searches can exclude DM masses below m χ 140 GeV. The DM relic density calculation for this benchmark is also shown in Figure 10 . This calculation is performed using MadDM [83] 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, overproduction or underproduction of DM should not be interpreted as an argument for excluding certain parameter ranges of the model [138,177].

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 miss T . The searches that are considered in this review to constrain the 2HDM+a parameter space include the mono-Higgs, mono-Z, tW + E miss T and H ± → tb channels -see also the recent ATLAS note [154]. The HL-LHC prospects of the tW + E miss T and four-top channel have been discussed in [139,147]. All existing 2HDM+a collider studies [32,33,134,138,[142][143][144][145][146][147][148][149][150][151][152][153][154] 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 GeV range. It is well-known (cf. for example [159,160]) 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.

2HDM+s model
Instead of mixing an additional CP-odd singlet P with the 2HDM pseudoscalar A, as done 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]. Like 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 miss T 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 will give an example of a 2HDM+s model realisation to which the DM DD searches are blind. Detailed discussions of the 2HDM+s model have been presented in [34,35,142] and we will summarise the most important findings of these articles below.

Theory
The construction of the DM-mediator interactions and the full scalar potential of the 2HDM+s model proceeds like in the case of the 2HDM+a model with minor modifications. First, the mediator S couples to the dark scalar current χχ and not to the pseudoscalar current χγ 5 χ. Second, the coupling of the S to the term (H † 1 H 2 + h.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,142] the mixing angle θ between H and S has been defined such that for sin θ = 0 the state S 2 is a pure S state. We will 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 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 The SI DM-nucleon cross section takes the general form where c N is the Wilson coefficient of the dimension-six nucleon operator O N = χχNN 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,142]  0.89 is the effective gluon-nucleon coupling. Besides the obvious possibilities to suppress the Wilson coefficient c N , i.e. y χ → 0, θ → 0 or m S 1 → m S 2 , in the case of the 2HDM of type-II there is a fourth possibility to achieve that c N 0. The trick is to choose tan β such that the square bracket in (30) vanishes. Numerically, this happens for tan β 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 GeV Higgs to provide the DM portal [100][101][102][103][104][105][106] 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 σ SI . See [35] for a detailed discussion of this point.

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 [142] the ATLAS and CMS searches for the h(bb) + E miss T [146,178] and the Z ( + − ) + E miss T [179,180] final states were examined, and it was found that like 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  Figure 11. Feynman diagrams that give rise to a resonant mono-Higgs signal (left) and a mono-Z signature (right) in ggF production in the 2HDM+s model.
in Figure 11. Notice that resonant contributions to the tW + E miss T 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 − aW − vertex by a H − S 2 W − vertex. The tW + E miss T signature has so far not been studied in the context to the 2HDM+s model.
The analysis performed in [142] employs the generalised alignment limit and uses the following benchmark parameter choices where λ S1 and λ S2 are the analogues of λ P1 and λ 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 λ 3 and the vanishing quartic couplings λ S 1 and λ 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 λ S 1 cos 2 β + λ S 2 sin 2 β, while the AS 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 β plane (right panel). The shown exclusions have been obtained in [142] by recasting the ATLAS searches for h(bb) + E miss T [178] and Z ( + − ) + E miss T [179]. 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 twodimensional parameter planes. For comparison, also the exclusions in the 2HDM+a model using the same input are 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 → A → S 2 Z and gg → H → aZ. However, for a scalar and pseudoscalar of the same mass and with the same couplings to quarks one has approximately σ(gg → A)/σ(gg → H) 2 and BR(A → S 2 Z)/BR(H → aZ) 1/2 [142], which leads to the same signal strength in both cases. For the mono-Higgs signal, the relevant resonant processes are instead gg → S 1 → S 2 h and gg → A → ah. In this case one has however that σ(gg → S 1 )/σ(gg → A) 1/2 and BR(S 1 → S 2 h)/BR(A → ah) 1 [142], 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    (31) and are taken from the publication [142]. The solid blue contour in the right plot corresponds to the parameter space that is excluded at 90% CL by the XENON1T upper limit on the SI DM-nucleon cross section [80]. For further details consult the main text. 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 [80] do not lead to any constraint because of the parameter choice tan β = 1. The situation is different for the m S 2 -tan β plane (right panel) where tan β 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 become 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 β = O(1).
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 [135] 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 [138] 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,142] give a detailed account of the relevant E miss T signatures, the DM DD and ID bounds and the relic density calculation, one would hope that combined analyses such as [154] will be done in the future in the case of the 2HDM+s model as well.

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 involve extensions of both the Higgs and gauge sector of the SM. The resulting theories fall into the class of dark Z or dark Higgs models, which typically have a rich collider and DM phenomenology. In the following we will 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.

2HDM+Z model
The dark Z model that has received by far the most attention at the LHC (see for example [146,151,152,178,[181][182][183][184]) is the so-called 2HDM+Z models introduced in [38]. In this model only the right-handed up-type quarks are assumed to be charged under the U(1) Z symmetry, while all the other SM fermion fields are chosen as SM singlets. This choice allows LHC production of the Z boson in uūand cc-fusion, but avoids the stringent constraints from searches for dilepton resonances since the charged leptons carry no U(1) Z charge. The Higgs sector of the 2HDM+Z 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 χ is a spin-1/2 or a spin-0 particle have been sketched in [38]. In the first case χ is a Majorana fermion that arises from singlet-doublet DM [185][186][187][188], while in the second case χ 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 χ need to be tuned, but after tuning it is possible to obtain large branching ratios of A to DM for χ masses in the ballpark of 100 GeV in both cases. For concreteness we will assume hereafter that the DM candidate in the 2HDM+Z model is a Majorana fermion that couples to the A with the coupling strength g χ . The same assumption is made in the ATLAS and CMS analyses [146,151,152,178,[181][182][183][184]. In addition we will focus our discussion on the LHC phenomenology of the 2HDM+Z model, ignoring the DM phenomenology because it is more model dependent.

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 model. In this limit the partial decay widths of the Z boson that are relevant for the discussion hereafter take the following form Here g Z denotes the U(1) Z gauge coupling, we have used the abbreviations x i/j = m 2 i /m 2 j and have neglected contributions that are suppressed by the mass ratio m 2 Z /m 2 Z of the Z-boson and the Z -boson mass. Such terms arise from Z-Z mixing, but are numerically subdominant in the partial decay rates of the Z boson. Notice that the Z boson can only decay to up-quark and charm-quark pairs as well as to top-quark pairs if m Z > 2m t with m t 173 GeV denoting the top-quark mass, but not to the other SM fermions because these matter fields are assumed to be U(1) Z singlets. If the A is sufficiently heavy, the coupling g χ is sizeable and tan β 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 Notice that the partial decay width A → Zh vanishes in the alignment limit, but depending on the precise value of tan β and/or the masses of the heavy 2HDM spin-0 states the decays A → bb, A → τ + τ − , A → HZ and A → H ± W ∓ can be relevant in this limit. From the expressions for the partial decay widths (32) and (33) it follows that if a heavy Z boson is produced on-shell in the LHC collisions this gives rise to a dijet or a tt final state, a mono-Higgs signature or resonant Zh production. Examples of the corresponding Feynman diagrams are shown in Figure 13. In our numerical study of the 2HDM+Z model we will consider all four collider signatures that arise from the exchange of Z bosons, and study only 2HDM+Z benchmarks in which the dominant A decay modes are described by the expressions (33).
Other laboratory constraints on the 2HDM+Z 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 model, one obtains at tree level a strictly positive correction to T of the form where sin θ w denotes the sine of the weak mixing angle θ w with sin 2 θ w 0.23, respectively, and α 1/128 is the electromagnetic coupling constant in the MS scheme at the mass of the Z boson. A simultaneous fit to the EWPO gives [3] T = 0.03 ± 0.12 ,  (14) is custodially invariant (see the related discussion in Section 3.1). Below we will only consider this limit. As explained in Section 3.1 indirect constraints on the charged Higgs-boson mass arise from flavour physics with B → X s γ providing a particularly strong 95% CL lower limit of m H ± > 800 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 model but turn out to be less restrictive than the Z -boson search strategies discussed below.

Summary plots
The existing LHC analyses [146,151,152,178,[181][182][183][184] that consider the 2HDM-Z model employ the alignment limit and adopt the following parameter choices: This benchmark has been established in [37] and we will also focus on it in what follows. The phenomenological impact that modifications of (36) have will however also be 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 Notice that this bound can be relaxed by decreasing the values of g Z and/or tan β 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 model in the parameter space that is testable at the LHC. Dedicated searches for the 2HDM-Z model have been carried out by ATLAS and CMS in mono-Higgs final states with the Higgs decaying to bottom-quark [146,152,178,181,183], photon [151,182] or tau pairs [182] as well as a combination of all relevant Higgs decay channels [184]. The dominant contribution to the mono-Higgs signal arises in the 2HDM+Z model from the graph shown in the middle of Figure 13 with subleading effects stemming from cc-fusion and pp → Z → Zh followed by Z → νν.
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 model. That this is the case can simply be seen by evaluating the partial decay width of the Z boson. Using (36) together with m Z = 2 TeV and m A = 800 GeV and employing exact formulas for all possible Z decay channels gives in good agreement with the approximations (32). One furthermore has BR(A → χχ) = 41.8% , BR(A → tt) = 57.7% .
The numbers (38) and (39) for the branching ratios of the Z and the A indicate that besides the mono-Higgs channel, also LHC searches for dijet final states (see [189,190] for the latest ATLAS and CMS results), tt resonances (see for example [191,192]) and resonant Zh production [193][194][195][196][197][198] can be expected to lead to relevant constraints in the m Z -m A plane for benchmark parameter choices like (36 (36). The solid blue, solid green and dashed green domains are taken from [151], [152] and [183], respectively. The solid orange, dashed cyan and solid red exclusions have instead been obtained from a recast of the dijet search [189], the tt search [191] and the search for resonant Zh production [198]. For comparison also the indirect constraints arising from the T parameter (dashed brown line) and B → X s γ [157,158] (dashed grey line) are shown. See text for further explanations. existing dijet, tt and Zh resonance searches in the 2HDM-Z model do not exist, but recasts are straightforward as we will show below.
In Figure 14 we summarise the relevant 95% CL constraints on the 2HDM-Z model in the m Z -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 [152] and CMS [183]. The searches h(bb) + E miss T provide stronger constraints than h(γγ) + E miss T because BR h → bb /BR(h → γγ) 250. The sensitivity decreases with increasing m A due to other decay channels like A → tt becoming kinematically accessible, thereby reducing BR(A → χχ). 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 bb + E miss T system, which is better suited for resonant signatures. As a consequence the CMS search is competitive with the ATLAS one despite using not 139 fb −1 but only 36 fb −1 of integrated luminosity collected at 13 TeV.
The solid orange, dashed cyan and solid red vertical lines indicate the 95% CL limits on m Z that follow from our recast of the dijet search [189], the tt search [191] and the resonant Zh search [198], respectively. These bounds are essentially independent of the precise value of the mass of the A, and for m A m Z read m Z > 3690 GeV , m Z > 3520 GeV , m Z > 2820 GeV .
Using [190] and [192] would lead to a slightly weaker dijet and tt 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  [189], the tt resonance search [191] and the search for resonant Zh production [198], respectively. The parameter space excluded by the T parameter is also displayed as a dashed brown line. Further details can be found in the main text.
obtained dijet and tt constraint are so strong that they exceed the maximal mass reach of the existing mono-Higgs searches that reads 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 ± in conflict with B → X s γ if the Yukawa sector of the 2HDM-Z model is taken to be of type II as done in [38] and all existing LHC analyses [146,151,152,178,[181][182][183][184]. Notice that the recent ATLAS Zh search [198] considers the combination of the + − bb and ννbb final states to set limits. In our recast we have only included contributions with a Z Zh vertex that lead to both of these final states (cf. the right diagram in Figure 13). Contributions with a Z Ah vertex where the pseudoscalar decays via A → χχ have instead been neglected (cf. the middle diagram in Figure 13). While such diagrams mimic the ννbb final state they lead to a different E miss T distribution than the former contributions with Z → νν 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, tt and Zh searches in probing the 2HDM-Z parameter space we show in Figure 15 the 95% CL exclusions in the m Z -g Z plane that follow from our recasts of [189], [191] and [198]. We employ the parameters given in (36) with m A = 1.5 TeV while keeping the U(1) Z coupling g Z as a free parameter. One observes that the used dijet (tt) search allows to set better constraints on g Z than the used Zh search in the mass region m Z 1.5 TeV (m Z 2 TeV) for lower Z -boson masses the situation is reversed. For m Z 1 TeV low-mass dijet searches such as [199][200][201][202][203][204] allow to set further relevant constraints on the 2HDM-Z 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 model at the LHC. In fact, as we have argued stringent constraints arise in general from searches that look for resonant Z -boson production in final states such as dijets, tt or Zh. Similar to what has been shown in [148] for the case spin-1 single-mediator DM simplified models, each Z -boson resonance search is sensitive to complementary regions of the m Z -g Z parameter space, and only by combining the whole suite of searches one is able to exploit the full LHC potential. ATLAS and CMS interpretations of dijet, tt and Zh resonance searches in the 2HDM-Z model do no exist but could be easily added to the exotics search canon.

2MDM model
While simplified models with a single spin-1 boson mediating the DM-SM interactions may lead to unitarity violation [97][98][99], it has been shown that the introduction of an extra scalar mediator can cure this problem [98]. 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 [205] and opening up new channels for DM annihilation, which can allow for the relic density constraints to be met [39,206]. A concrete model implementation of this general idea is discussed in what follows.

Theory
The 2MDM model introduced in [39] includes an additional U(1) Z gauge symmetry with a new gauge boson Z and a complex scalar field S that is a singlet under the SM gauge group. DM is taken to be a Majorana fermion χ 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 Here g is the U(1) Z gauge coupling, y χ is a Yukawa coupling, q χ and q q denote the U(1) Z quantum numbers of the DM and the quark fields, respectively, and P L and P R project onto leftand right-handed fields. The U(1) Z charge of S is fixed to be q S = −2q χ 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(1) Z charge to avoid mass mixing between the Z and the Z boson which is severely constrained by the measurements of the EWPOs. In such a case there is also no Z Zh vertex at tree level, and as a result searches for resonant Zh production do not provide any constraint. Note that the Z 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 [207], and realistic models that implement this mechanism have been discussed for example in [208][209][210][211][212][213]. Models of this type generically have additional fermionic degrees of freedom to cancel gauge anomalies (cf. [214][215][216] 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 χ.
The scalar S is assumed to acquire a VEV S = 1/ √ 2 (0, w) T , spontaneously breaking the U(1) Z and thereby giving masses to the Z boson and the DM particles: Introducing now the coupling strengths g χ = g q χ and g q = g q q , the Lagrangian (42) can be rewritten as 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(1) Y and the U(1) Z 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 GeV Higgs can mix with the amount of mixing conventionally parameterised by sin θ.

Experimental constraints
The LHC searches [217][218][219] have provided interpretations in the 2MDM model for presenting exclusion limits in the m Z -m s plane. The used value of sin θ clearly satisfies the aforementioned limit on sin θ imposed by Higgs physics. Furthermore, the decay channels h → χχ and h → Z Z are kinematically closed. Since the mass m s is scanned over in the works [217-219], on-shell Higgs decays of the form h → ss are in principle allowed if m s < m h /2. In the limit sin θ 1, the corresponding partial decay width can be approximated by For the parameters (45) and taking m Z = 1 TeV and m s m Z , one finds numerically Γ(h → ss) = 1.9 · 10 −3 Γ SM h . Such a small modification easily passes the bounds from the direct measurements of the total Higgs width [161,162] as well as the numerous limits that stem from searches for h → 4 f (see for example [220] 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 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, tt or j + E miss T searches. Since all these signatures result from s-channel Z -boson exchange, an important ingredient that will determine the signal strength q q q Z q χ q q q Z χ g Figure 16. Representative Feynman diagrams that lead to dijet or heavy-quark pair production (left) and a mono-jet signal (right) in the 2MDM model.
in a certain channel is the branching ratio of the Z boson into the corresponding final state. The necessary partial decay widths take the following form Using the input parameters given in (45) together with m Z = 1 TeV, these expressions lead to BR(Z → χχ) = 25.5%, BR(Z → qq) = BR(Z → tt) = 12.4% and a relative width of Γ Z /m Z = 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 χ = 1, g q = 0.25 and g = 0. As a result, one expects to find the same hierarchy of sensitivities as in the spin-1 simplified DM case -see [221,222] for the latest DM summary plots for s-channel mediators by ATLAS and CMS -with the non-E miss T searches providing better bounds on m Z than the E miss T searches. We will 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 boson the displayed diagrams contain a dark Higgs s that can either be emitted in a Z Z s vertex (left graph) or a χχs vertex (right graph). The dark Higgs can decay into SM and DM particles. Assuming that the decay channel s → χχ is kinematically inaccessible, the relevant partial decay widths are given in the limit sin θ 1 by 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 [223]) which is expected because the dark Higgs obtains its SM coupling solely by mixing with the 125 GeV Higgs. In consequence, the decay pattern of the mediator s resemble to first approximation those of a SM Higgs boson with mass m s . In the mass range 50 GeV m s 140 GeV the dom- 50% at the upper end of the considered mass range. The latter feature is expected from the SU(2) L × U(1) Y symmetric limit. For m s values above the top-quark threshold also the decay s → tt is relevant. In the following numerical analysis, we restrict ourselves to the mass range 50 GeV < m s < 350 GeV, and thus the dominant final states that arise from the graphs in Figure 17 are bb + E miss T , VV + E miss T with VV = W + W − , ZZ and hh + E miss T . 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 [189], tt [191] and mono-jet search [169] the following 95% CL limits on the Z -boson mass: The dijet [190], tt [192] and mono-jet [170] 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 miss T events with ISR jets (see the right diagram in Figure 16) while we have not considered contributions that follow from s (VV) + E miss T production with the EW gauge bosons decaying fully hadronically (see Figure 17). We will 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. [221,222]).
For the s + E miss T signatures, it is important to realise that the dark Higgs is produced with high momentum due to the large m Z -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 140 GeV. This decay mode is targeted by a reinterpretation [217] of the ATLAS h(bb) + E miss T search that considers the invariant mass distribution of the bb pair down to 50 GeV, allowing to probe m s > 50 GeV. For m s 140 GeV, the s → 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 [218]. The challenging multi-prong decay s → V (qq)V (qq) is reconstructed with a novel technique [224] 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 miss  [217]. The solid orange, dashed cyan and solid purple vertical lines follow from a reinterpretation of the dijet [189], tt [191] and mono-jet search [169], respectively. The pink hatching indicates the direction for which the DM relic density is larger than observed. Further details are discussed in the main text.
two-dimensional fit to the dilepton invariant mass and the transverse mass of the trailing lepton plus E miss T system. The constraints on the m Z -m s plane from the s + E miss T analyses by the ATLAS collaboration [217,218] are presented in Figure 18. The 95% CL exclusion limits (49) are also shown in the figure for comparison. Focusing on the E miss T searches first one observes that the novel s(VV) + E miss T and s(bb) + E miss T search strategies allow to exclude additional parameter space of the 2MDM model with respect to our mono-jet recast of the ATLAS search [169]. In this context one however has to remember that the mono-jet limit given in (49) does not include the contributions of the s(VV) + E miss T and s(bb) + E miss T processes. The obtained mono-jet bound therefore provides a conservative lower limit on the actual sensitivity of the j + E miss T search [169] for the 2MDM benchmark (45) -we estimate that including s + E miss T contributions would strengthen the mono-jet limits by around 15% (5%) at m s = 100 GeV (m s = 200 GeV). It is also evident from the figure that our recast of the dijet [189] and tt search [191] has a higher mass reach in m Z than the mono-X searches. We note that the 2MDM model can also give rise to hh + E miss T and tt + E miss T signatures. The former signal can be constrained for instance by using the results of [225,226] that studies final states with at least three b-jets and E miss T (see also [227] for an exploration of the hh + E miss T signature in simplified models of hidden sectors). Estimating the sensitivity of multi b-jet plus E miss T searches to the hh + E miss T signature in the 2MDM model is however beyond the scope of this review. The latter signal can be targeted by standard tt + E miss T searches (cf. for example [145,166,168]). Employing [168] we expect that for (45) the HL-LHC with 3 ab −1 may be able to set the bounds M Z 500 GeV for m s 750 GeV and M Z 1.5 TeV for m s 400 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 [228].

Models with exotic Higgs decays involving LLPs
BSM scenarios with hidden sectors that are connected to the SM sector through the 125 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 non-trivial and often requires dedicated triggering and reconstruction algorithms -see [229,230] 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 [231]. Below we will consider three hidden sector models in which the LLPs arise from decays of the 125 GeV Higgs boson. The collider phenomenology of axion-like particles (see for example [232,233] for recent comprehensive studies) that might also lead to LLP signatures from Higgs decays will not be discussed.

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 like 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 [234]. 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α s that is comparable in strength to its SM counterpart α s . Connecting DM with neutral naturalness is possible (see for instance [235][236][237][238][239][240][241][242][243][244][245][246]) but more model-dependent than the LLP phenomenology on which we will focus in the following.

Theory
In models of neutral naturalness the coupling between the 125 GeV Higgs boson and the top partners gives rise to an effective coupling hĝĝ between the 125 GeV Higgs h and the hidden gluonsĝ 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ĝĝ interactions can be parameterised by whereĜ a µν denotes the field strength tensor of the hidden SU(3) gauge group andζ is a model-dependent mixing angle. One generically has |ζ| = 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 |ζ| 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 [247,248].
The lightest mirror glueball is a scalar state with quantum numbers J CP = 0 ++ and its massm 0 is entirely determined by the RG running of the strong couplingα s in the hidden sector. In fact,m 0 is related to the mass m 0 1.7 GeV of the 0 ++ glueball in QCD [249][250][251] by the following simple rescaling [44] with Λ Landau 150 MeV the Landau pole in QCD, i.e. the scale where 1/α s (Λ Landau ) = 0. Without knowing the exact mass spectrum of states in the hidden sector it is not possible to give a precise value ofΛ Landau and thereforem 0 . One can however estimate the typical range ofm 0 values by studying the following toy model. Let M * be the scale where the discrete symmetryα s (M * ) = α s (M * ) 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β whereN f denotes the number of flavours that are active for renormalisation scales µ in the range M < µ < M * . Solving the one-loop RG equation ofα s it is easy to show that the Landau pole in our toy model occurs at where 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 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] BR h → 0 where BR(h → gg) SM 8.2% is the Higgs to digluon branching ratio in the SM. The parameter κ(m 0 ) encodes our ignorance about the hadronisation of the lightest mirror glueball and the mixing effects of excited hidden glueball states with the 125 GeV Higgs. In the article [44] it has been argued that κ(15 GeV) 0.1 and κ(50 GeV) 1. Using these numbers as well as the estimateζ v 2 /M 2 , one obtains for the two benchmarks that led to (55) the following range of branching ratios: 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 → 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 ++ → h * → YY. The corresponding partial decay widths have been calculated in the work [248] and in the case of (50) take the form wheref 0 is the annihilation matrix element of the lightest mirror glueball through (50) and Γ(h * → YY) SM denotes the partial decay width for a SM Higgs boson with massm 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 bb, cc and τ + τ − final states, respectively. Using again ζ v 2 /M 2 as well asα s (m 0 )f 0 0.18m 3 0 [44,251] together with (58) the proper decay length of the lightest mirror glueball can be approximated by For typical values ofm 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 bothm 0 and M suggests that cτ 0 ++ can in practice be treated as an almost free parameter in the framework of neutral naturalness. Notice that hidden valley models [252][253][254] 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 π h . Given its pseudoscalar nature the hidden pion preferentially decays to heavy SM flavours for example π h → bb. As argued in the articles [252][253][254], the typical π 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 GeV Higgs boson it is possible to obtain h → π h π h branching ratios that are observable at the LHC [252][253][254]. The phenomenology of the π 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 BR(h → ss) BR(h → aa) and its proper decay length cτ s (cτ 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 will hereafter refer to the LLP produced in exotic Higgs decays as an a particle.

Experimental constraints
The first searches for pair-produced neutral LLPs in the context of Higgs portal models have been performed by the CDF [255] and DØ [256] 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 have been 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 → aa → 4b optimised for prompt decays or small proper decay lengths cτ a 6 · 10 −3 m [257,258]. 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 → h → aa with a decaying to hadronic jets [259], 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 m to 1 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.  Figure 19. Representative diagrams that can give rise to LLP signatures at the LHC. The left and right graph displays associated Zh production followed by Z → + − and h → aa and ggF Higgs production with h → aa, respectively. The a dominantly decays to bottom-quark pairs leading in both case to a four-bottom final state. Most existing searches have targeted this final state but analyses that look for multi-jet events have also been performed.
In order to significantly increase the efficiency for reconstructing displaced vertices, dedicated reconstruction techniques, so-called large-radius tracking (LRT), were developed at ATLAS [260]. Properties of displaced vertices reconstructed from LRT tracks like 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 → aa → 4b and Z → + − . The leptons from the Z → + − decays were used for triggering, and backgrounds were essentially eliminated by requiring two displaced vertices in candidate events [261].
CMS has recently searched for h → aa → 4b and Z → + − using the same process [262]. A trigger and selections based on dilepton Z-boson decays provide sensitivity to light LLPs with masses of 15 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 → aa → 4b process targeting cτ a between 10 −3 m and 3 m by requiring displaced vertices consistent with LLP decays to be reconstructed in the inner tracker [263,264]. Unlike the ATLAS analysis [261], the CMS search [264] 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 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 [264] provides currently the best sensitivity for proper decay lengths between 10 −3 m and 10 −1 m. CMS has also searched for displaced leptons arising from ggF Higgs production followed by h → aa → 4 [265]. Since this search assumes that the LLP has equal probability to decay to two muons and two electrons, whereas (58) implies that BR(a → µ + µ − ) = O(10 −4 ) and BR(a → e + e − ) = O(10 −8 ), an interpretation of [265] in the context of neutral naturalness/hidden valleys leads to no meaningful constraint on the h → aa branching ratio.
For mean decay lengths in the range of 1 m to 10 2 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.
ATLAS searched for h → aa → 4 f decays in the CM targeting event topologies compatible with ggF Higgs production [266]. The decay of a LLP inside the CM is typically reconstructed as a single jet with striking characteristics, namely a narrow width in pseudorapidityazimuthal-angle (η -φ) 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 [267] and when  [257,261,266], CMS [262,264] and LHCb [259] publications, respectively. For comparison the 95% CL limits on the branching ratio of undetected and invisible Higgs decays (60) are also displayed. The boundary between the constraints from undetected and invisible Higgs decays is taken to be cτ a = 3.5 GeV corresponding to m a = 35 GeV. For further explanations consult the main text.
performing the offline analysis employing an artificial neural network. ATLAS recently performed a search for the same process using the MS [268]. 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 [267] based on the overall activity in the MS. The backgrounds are essentially eliminated by requiring two displaced vertices from candidate h → aa → 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 has been applied in a similar analysis [269] and the results were statistically combined with the CM search [266]. Finally, ATLAS carried out a search using the combination of the ID and the MS, targeting short and long mean decay lengths [270].
CMS performed a LLP search using the endcap of the MS [271], targeting cτ a between 10 −3 m and 10 2 m and decay chains h → aa with a → bb, cc, τ + τ − . While many features are similar to the corresponding ATLAS searches [268,269], 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 miss T vector was used as the final discriminant. ATLAS and CMS also performed other searches for neutral LLPs decaying to jets targeting SUSY models [272][273][274][275][276] that are not explicitly optimised for h → aa → 4 f signatures, for example due to very high trigger thresholds, and therefore will not be discussed below.
An indirect way to constrain models where the 125 GeV Higgs boson decays into LLPs is provided by the combination of precision measurements of the Higgs couplings in visible and invisible final states [74,108,109]. In particular, under the assumption that the coupling modifiers κ V = g hVV /g SM hVV of the 125 GeV Higgs boson to EW gauge bosons satisfy κ V ≤ 1, which generically holds in models with Higgs mixing, the following 95% CL constraints can be placed [108]. Here, the undetected category includes events with undetected BSM particles that do not provide a significant E miss T contribution such as the ones typically selected by the aforementioned LLP analyses with short mean decay lengths below O(1 m). In contrast, LLPs that decay outside the tracker and calorimeters, i.e. having mean decay lengths larger than O(1 m), are not registered by standard reconstruction algorithms, resulting in a E miss T contribution. Hence, for moderate boosts the above bound on BR(h → undet) and BR(h → inv) can be interpreted as an indirect limit on BR(h → aa) for proper decays lengths of cτ a 1 m and cτ a 1 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τ a 0.1 m (m a /GeV) [277,278].
The 95% CL limits of the above searches on the branching ratio BR(h → aa) as a function of the proper decay length cτ 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τ a of a few meters the LHC searches are able to set a limit BR(h → aa) 10 −3 . In view of (57) and (59) this is an interesting finding, as it allows to test certain model realisations of neutral naturalness. Notice that the limits from BR(h → undet) and BR(h → inv) provide currently the strongest LHC bound on BR(h → aa) for cτ a 10 −4 m and cτ a 10 3 m, respectively. Besides these two constraints, proper decay lengths cτ a 10 3 m are currently unexplored by collider measurements. Dedicated detectors like MATHUSLA [231], CODEX-b [279] and ANUBIS [280] may address such long mean decay lengths at the HL-LHC.

Dark photons
In the case of neutral naturalness the LLP is a composite spin-0 particle. However, hidden sector models with a spin-1 LLP also exist. The model considered in this subsection is based on an extra U(1) 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 GeV Higgs boson to a pair of dark photons Z d followed by Z d → + − . 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 will not discuss the dark photon searches [281][282][283] which are motivated by the theoretical works [284][285][286][287]. 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 [288][289][290].

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 where B µν = ∂ µ B ν − ∂ ν B µ and X µν = ∂ µ X ν − ∂ ν X µ is the U(1) Y and U(1) X gauge field strength tensor, respectively, denotes the hypercharge mixing parameter, while κ is the Higgs portal coupling. After EW symmetry breaking H = (0, v/ √ 2) T and spontaneous symmetry breaking of the U(1) X symmetry by S = v S / √ 2, the mass spectrum of the hidden sector model contains two heavy neutral gauge bosons and two neutral Higgs bosons. We will denote these states by Z, Z d , h and h d . For small (κ) 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 O( 2 ), while the dark photon vector couplings receive corrections already at O( ). Explicitly, one has where e = √ 4πα is the electromagnetic coupling and Q f (Y f ) is the electric charge (hypercharge) of the relevant fermion. Notice that for m Z d m Z the Z d -boson coupling to fermions is photon-like, while for m Z d m Z the Z d boson couples to fermions like 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 will however follow this established naming convention. At O( ) the hypercharge portal also leads to a coupling between the 125 GeV Higgs, a Z and a Z d boson. To this order the coupling takes the form For O(10 −4 ) a dark photon with mass m Z d 1 GeV decays promptly [45,46]. At the LHC, such dark photons can be searched for in Drell-Yan (DY) dimuon production pp → Z d → µ + µ − and in four-lepton final states that arise from the process pp → h → ZZ d → 4 . 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 GeV Higgs and two dark photons. To lowest order in κ this coupling can be written as For m Z d < m h /2 this coupling allows for Higgs decays of the form h → Z d Z d -a representative Feynman diagram is shown on the right in Figure 21. At the same order in κ also interactions between the 125 GeV Higgs and two dark Higgses and a single dark Higgs and two dark photons exist. In the following, we will assume that the dark Higgs is heavy, i.e. m h d m h /2. In such a case the exotic Higgs decay h → h d h d is kinematically forbidden and the production cross section for pp → h d → Z d Z d is always smaller than pp → h → Z d Z d . As a result, for sufficiently light dark photons the best probe of the Higgs portal parameter κ is the exotic Higgs decay h → Z d Z d . The corresponding partial decay width is given to leading order in κ by the following expression [45,46] 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 κ 1. For example, taking κ = 0.07, m Z d = 30 GeV and m h d = 300 GeV leads to BR(h → Z d Z d ) = 0.15, which is close to the model-independent limit on BR(h → undet) reported in (60).
The best way to constrain the Higgs portal coupling now depends on the size of kinetic mixing. For O(10 −4 ) the dark photon decays promptly and one can again study fourlepton final states [291,292] to constrain the h → Z d Z d branching ratio and thereby κ. For O(10 −4 ) the decay length of the dark photon starts to become macroscopic. In fact, for Figure 21. Example contributions in the dark photon model (61) to DY dimuon production (left diagram) and four-lepton final state production arising from the exotic Higgs decays h → ZZ d (middle diagram) and h → Z d Z d (right diagram), respectively. In the latter two cases the production mechanism of the 125 GeV Higgs boson is not shown.
O(10 −8 ) O(10 −4 ) 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.

Experimental constraints
The ATLAS, CMS and LHCb collaborations have carried out searches for prompt dark photon decays targeting hypercharge mixing parameters 10 −4 in both the DY channel pp → Z d → µ + µ − [293][294][295] and in four-lepton production associated with h → ZZ d and h → Z d Z d [291,292]. ATLAS has also searched for LLPs in the h → ZZ d channel [296]. In addition, smaller values have been tested by dedicated searches for long-lived dark photons in final states with displaced dimuon vertices arising from h → Z d Z d . Building on earlier LHC Run 1 analyses [297,298], the latest ATLAS and CMS results of this kind can be found in [299] and [300], respectively. In the following we review only LHC searches for LLP signatures involving dark photons.
From all possible dark photon decay modes, the Z d → µ + µ − 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 → µ + µ − decay mode targeting proper decay lengths of 10 −3 m < cτ Z d < 10 3 m by requiring displaced vertices consistent with LLP decays to be reconstructed using the muon system [299]. A combination of single muon, dimuon and trimuon triggers with progressively lower thresholds down to p T > 6 GeV, as well as a E miss T 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 m. Finally, displaced vertices of oppositely charged muons with an invariant mass above 15 GeV are selected within the fiducial volume.
CMS also performed a search for dark photons in µ + µ − final states [300]. 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 m < cτ Z d < 10 −1 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 [295,301], 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 GeV. Finally, candidate muon pairs from the Z d → µ + µ − decay are selected in events with two muons targeting DY production, and four muons with an invariant mass consistent with the 125 GeV Higgs boson aiming for pp → h → ZZ d and pp → h → 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 → µ + µ − decays away from the known SM resonances like the J/ψ and Υ mesons.
The region with cτ Z d < 10 −3 m is covered by the CMS search for displaced Z d → e + e − , µ + µ − decays performed at 8 TeV [297], which targets the proper decay length range between 10 −4 m and 10 2 m. This analysis is similar to the 13 TeV CMS search [295], 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 BR(h → Z d Z d ) as a function of the proper decay length cτ Z d of the dark photon. The shown ATLAS limits [299] correspond to the different dark photon masses m Z d = 20 GeV, 60 GeV, while in the case of CMS we display bounds for m Z d = 20 GeV, 50 GeV from [297] and for m Z d = 0.65 GeV, 2 GeV, 15 GeV, 30 GeV, 50 GeV from [300]. Notice that for dark photons with masses of a few tens of GeV, the strongest bound on BR(h → Z d Z d ) typically arises for the lightest choice of the dark photon mass. For m Z d 5 GeV the limits on BR(h → Z d Z d ) however become significantly weaker with decreasing dark photon mass even if data scouting techniques are employed. Specifically, for dark photons with masses above 15 GeV the ATLAS (CMS) measurements exclude BR(h → Z d Z d ) 10 −4 (10 −5 ) for proper decay lengths of around 10 −1 m (10 −2 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 m to 10 2 m the bound on BR(h → Z d Z d ) that derives from the LLP searches [299,300] Figure 23. 95% CL excluded regions in the m Z d -plane following from the ATLAS [299] and CMS [300] measurement. In the case of the ATLAS (CMS) results the shown solid red (dashed blue) bound corresponds to BR(h → Z d Z d ) = 10 −3 . The 90% CL upper limit on the hypercharge mixing parameter that derives from DY dimuon production by LHCb [294] and CMS [295] is also displayed for comparison as a solid green and a dashed yellow line, respectively. The 95% CL limit on from the ATLAS search [291] for pp → h → ZZ d → 4 is finally indicated by a solid orange line. The latter three bounds hold for any value of the h → Z d Z d branching ratio. Finally, the indirect bound that follows from the limit (60) on invisible Higgs decays is shown as a dashed black line. For further details consult the text. stronger than the model-independent limits quoted in (60). For shorter or longer proper decay lengths the ATLAS Higgs coupling measurement [108] however provides currently the best constraint on the branching ratio of h → 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 done in Figure 23 where we show the constraints that derive from the ATLAS [299] and CMS [300] measurement assuming BR(h → Z d Z d ) = 10 −3 . For the same h → Z d Z d branching ratio, the exclusion that follows from [297] almost exactly resembles that of [300] for dark photon masses between 20 GeV and 50 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 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 [294,295] as well as the ATLAS search [291] that focuses on the pp → h → ZZ d → 4 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, searches for pp → h → Z d Z d → 4 [291,292] 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 [74,108]. Imposing BR(h → inv) = 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τ Z d > 0.1m (m Z d /GeV), resulting in a E miss T contribution. Another assumption is that the kinematic distributions such as the E miss T spectrum are the same for pp → h + X → Z d Z d + X and for the SM Higgs production Figure 24. Left diagram: Higgs decay to a pair of the lightest MSSM neutralinos. Middle and right diagram: Possible decays of the MSSM bino to the hidden sector, which follow from gaugino kinetic mixing (66). See text for further details. channels that go into the BR(h → inv) bounds given in (3) and (60). While [278] 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.

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 GeV Higgs boson, but through a cascade involving SUSY and hidden sector particles with masses of O(10 GeV) or below. This implies that the hidden photons have to be even lighter than the other particles in the decay chain, with masses of O(1 GeV) 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 miss T depending on the precise pattern of the Higgs decay.

Theory
A minimal model that realises the general idea proposed in [47][48][49][50][51] contains a massive hidden photon γ 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 (B) and the hidden gaugino (γ d ).
Removing the kinetic mixing between theB and theγ d then gives rise to interactions between all hidden fields charged under U(1) d and the visible neutralinos that are proportional to the hypercharge mixing parameter . In particular, one obtains an interaction term of the following form where g d is the U(1) d gauge coupling, h i d are the hidden scalar fields,h i d are the hidden gauginos and q i is the relevant U(1) 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 Figure 25. Benchmark topologies that have been studied in the context of the FRVZ model at the LHC. In the first process (left diagram) both neutralinos produced in h →Ñ 1Ñ1 decay throughÑ 1  neutralinos (Ñ 1 ), as indicated by the Feynman diagram on the left-hand side of Figure 24. In the pure MSSM without the hypercharge portal, theÑ 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 thẽ 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Ñ 1 decays to the lightest hidden sector gaugino (Ñ d ) and a hidden photon, while in the second case the final-stateÑ d is accompanied by a hidden Higgs boson (h d ). The lightest hidden sector neutralinoÑ d is stable and therefore a DM candidate, but the γ d and the h d decay further. If the hidden photon is the lightest hidden state then γ d → + − and h d → γ d γ d → 4 are possible decay chains that will 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 [302][303][304][305][306][307]) also possess many of the features discussed above. Such models typically contain two fermionic states χ 1 and χ 2 with a small mass splitting ∆m = m χ 2 − m χ 1 > 0. The role ofÑ 1 (Ñ d ) is played by χ 2 (χ 1 ) in inelastic DM models, and the DM candidate χ 1 can be excited to its heavier twin χ 2 by absorbing a massive dark photon γ d . The simplest realisation of such a scenario consists in assuming a hypercharge portal and postulating a U(1) 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 γ dÑ1Ñd coupling (66), the γ d χ 1 χ 2 coupling is not proportional to the hypercharge mixing parameter , but only involves the U(1) 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 γ d χ 1 χ 2 , γ d ff , h d χ 1 χ 2 and hχ 2 χ 2 couplings are effectively treated as free parameters. As a concrete example, the ATLAS collaboration has used the following parameters in their interpretations g γ d χ 1 χ 2 = 0.31 , g h d χ 1 χ 2 = 0.1 , where the choice of g γ d χ 1 χ 2 corresponds to g γ d χ 1 χ 2 = e = √ 4πα with α the electromagnetic fine structure constant at the EW scale. The hypercharge mixing parameter and the coupling g hχ 2 χ 2 are not directly used as external parameters, but expressed through cτ γ d and BR(h → χ 2 χ 2 ), respectively, which then serve as input. Notice that for the parameter choices (67) and sufficiently large BR(h → χ 2 χ 2 ) values, the decays h → χ 2 χ 2 and χ 2 → χ 1 γ 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 BR(h → χ 2 χ 2 ) is identified with BR(h → Z d Z d ).

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 [308,309] and displaced [310,311] 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 → γ d γ d + X topology (cf. left diagram in Figure 25), where the hidden photon is assumed to decay into charged leptons. The ATLAS search [308] included γ d decays to e + e − and µ + µ − , while the CMS search [309] 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 GeV to 2 GeV for ATLAS and 0.25 GeV and 2 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 ∆R = (∆η) 2 + (∆φ) 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 [310] has probed the h → 4γ 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 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 [311] explored the h → γ d γ d + X topology (cf. left diagram in Figure 25), considering only γ d → µ + µ − decays. Hidden photon masses in the range 0.25 GeV and 8.5 GeV were explored, with the upper limit set by the requirement m µµ < 9 GeV in order to sufficiently suppress the background from DY and Υ → µ + µ − production. Proper decay lengths below cτ γ d < 0.1 m were targeted, explicitly  Figure 26. 90% CL upper limits that follow from the ATLAS search [310] and the CMS search [311] for the exotic Higgs boson decay h → γ d γ d + X in the generalised FRVZ model. Both exclusions assume BR(h → γ d γ d + X) = 0.1. For comparison also constraints in the m γ d -plane from LHCb [294], BaBar [312], KLOE-2 [313] and ν-CAL I [314] are shown. The latter bounds hold for any value of the h → γ d γ d + X branching ratio. Finally, the indirect upper limit that follows from (60) is depicted.
Consult the text for further explanations.
including both prompt and displaced signatures. Candidate events were required to have exactly two γ d → µ + µ − 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 γ d -plane in the generalised FRVZ model for hidden photon masses in the range of 0.1 GeV to 10 GeV. The constraints following from the ATLAS [310] and CMS [311] search for the exotic Higgs decay h → γ d γ d + X apply in the generalised FRVZ model and assume BR(h → γ d γ d + X) = 0.1. The ATLAS analysis employs (67) and focuses on the mass range 0.2 GeV m γ d 3.6 GeV and small kinetic mixings that lead to LLP signatures. Depending on m γ d it excludes hypercharge mixing parameters within 1 · 10 −5 3 · 10 −7 . The CMS search probes hidden photon masses in the range 0.25 GeV ≤ m γ d ≤ 8.5 GeV and targets proper decay lengths of cτ γ d ≤ 0.1 m. This analysis is able to exclude hypercharge mixing parameters 3 · 10 −6 ( 7 · 10 −8 ) at low (high) m γ d . The displayed LHCb [294], BaBar [312], KLOE-2 [313] and ν-CAL I [314] bounds have been taken from the Darkcast package developed in [315]. They hold irrespectively of the value of the h → γ d γ 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 BR(h → inv) = 0.1. The limit arising from the recent ATLAS analysis [278] of mono-jet signatures is not shown in the figure because it is not sensitive to BR(h → γ d γ d + X) = 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 like MATHUSLA, CODEX-b, ANUBIS, FASER [316] and FASER2 [317] located at the LHC, in combination with Belle II [306,307], are expected to set additional stringent constraints on the m γ d -plane in the generalised FRVZ or inelastic DM frameworks.

Outlook
In this article we have reviewed the status of the LHC constraints on DM scenarios where the 125 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 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 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 miss T searches for off-shell Higgs production in the VBF, tt and tW, and potentially other channels provide sensitivity to the parameter space where m 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 σ SI bounds do not become constant in the limit m DM 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 BR(h → inv) as well as the limits on σ 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 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 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 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 → bb, γγ in the former case and on the mono-s channel with s → bb, 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 boson in dijet or tt 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 and 2MDM models, searches for X + E miss T signals are important given their distinct experimental signature. We have furthermore pointed out that the 2HDM-Z model can also be probed by a search for Z → Zh resonances. Similarly, we advocate the exploration of the hh + E miss T and tt + E miss T 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.2 and 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 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 GeV Higgs properties in visible and invisible final states. For instance, in the case of neutral naturalness we find that the limits from BR(h → undet) and BR(h → inv) provide currently the strongest LHC bound on BR(h → aa) for cτ a 10 −4 m and cτ a 10 3 m, respectively. To further emphasise the possibility to test LLP scenarios via E miss T searches, we have in the case of the two models with a hypercharge portal performed recasts of the limits imposed by the latest h → 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.