Singlet-Doublet dark matter freeze-in: LHC displaced signatures versus cosmology

We study the Singlet-Doublet dark matter model in the regime of feeble couplings, where the dark matter abundance is obtained via the freeze-in mechanism. As a consequence of the small couplings, the heavier particles in the model are long-lived with decay length at typical scales of collider experiments. We analyse the collider signatures of the model, characterised by displaced h and Z bosons plus missing momentum, employing current LHC searches for displaced vertices and missing energy to significantly constrain the parameter space of the model. We also take into account the cosmological bounds relevant for our light dark matter candidate arising from Lyman-α forest constraints. Our analysis emphasises the interplay between displaced signatures at the LHC and cosmology for dark matter candidates whose relic abundance is obtained through the freeze-in mechanism.


Introduction
Compelling observational evidence supports the existence of dark matter (DM), the most abundant form of matter in the Universe [1]. Yet, at present, the nature of the dark sector, possibly including dark matter particle(s) and new mediators driving the interactions between dark matter and the Standard Model (SM), is still unknown. In the last decades most of the attention has been devoted to the weakly interacting massive particle (WIMP) paradigm, in which dark matter is a new type of elementary particle with weak-type interactions with the SM. In WIMP models the abundance of the dark matter is obtained through the freeze-out mechanism. The weak couplings involved typically imply a viable -i.e. giving rise to the observed relic abundance -dark matter particle with mass of the order of the electroweak (EW) scale. This remarkable coincidence -the so-called WIMP miracle -has motivated a large effort in the study of WIMP dark matter scenarios both JHEP09(2018)037 in top-down approaches, triggered by supersymmetry, and bottom-up approaches focusing on effective theories and simplified models. In the latter framework, an important effort has been deployed in the last years to characterise dark matter simplified model parameter space, to search for the dark matter particle and the associated mediators at colliders, and to explore the complementarity between the LHC, the direct and indirect detection experiments, see e.g. [2][3][4]. The recurrent null results in the search for WIMPs (both at the LHC and in direct/indirect detection experiments) provide good motivations though to take a step back and to investigate alternative dark matter paradigms and to study thoroughly their phenomenology.
Here we focus on the so-called freeze-in mechanism for producing a feebly coupled dark matter candidate, i.e. particles that were not in thermal equilibrium with the SM in the early universe, see e.g. [5][6][7]. Because of the very small coupling involved, these models can give rise to displaced signatures at the LHC in terms of long-lived mediators decaying into dark matter plus SM fields, as it has already been underlined in a number of works [5,[8][9][10][11]. Long-lived particles and displaced signatures at colliders in relation to dark matter simplified models have been discussed also in [12][13][14][15][16][17][18]. The requirement of a viable frozen-in dark matter scenario giving rise to displacement at colliders of detector size points directly towards light dark matter candidates, with a mass of the order of the keV [5]. This is typically the mass scale currently tested by cosmology and astrophysics probes in the framework of warm dark matter scenarios, see e.g. [19][20][21], that lead to the suppression of the small scale structure formation and, by the same token, can help to alleviate the small scale crisis [22][23][24][25][26][27] in the ΛCDM (the Standard Cosmological scenario) see e.g. the discussion in [28][29][30][31][32][33][34][35]. On the other hand, due to the feeble coupling involved, direct and indirect detection dark matter searches are challenging, see however e.g. [36][37][38][39].
In this paper we will study how the combination of collider and cosmological constraints can probe a significant portion of the parameter space of frozen-in dark matter models. First, note that the long-lived mediators typically possess sizeable couplings with the SM particles that keep them in thermal equilibrium with the SM bath, and thus they can copiously be produced at the LHC. Second, the reach of displaced vertex signatures at the LHC can actually extend beyond the regime where the mean decay length of the heavy mediator, cτ , is within the tracking detector with typical radius 1 m. This is because of the background-free nature of displaced signatures and because of the exponential decay distribution. 1 In general, the study of displaced signatures with missing energy at the LHC is thus a powerful probe of frozen-in dark matter in the small mass regime (keV to a few MeV), complementing the existing constraints from astro-physics and cosmology on light free-streeming DM candidates.
For concreteness, we focus here on the Singlet-Doublet dark matter model [41][42][43][44][45][46][47][48][49][50][51][52]. In the Singlet-Doublet model, the SM is augmented with a pair of electroweak doublet Weyl fermions and one Majorana singlet, interacting with the SM Higgs via Yukawa couplings. The dark sector contains thus a charged fermion and three neutral fermions, whose lightest 1 The latter is indeed such that a substantial fraction of the mediators produced at LHC would decay within the detector even for cτ > 1 m, see e.g. [40].

JHEP09(2018)037
state constitutes a stable dark matter candidate. Within the freeze-in regime of such model, the Yukawa couplings are feeble and thus the dark matter particle is mostly singlet while the role of the long-lived mediators is played by the components of the electroweak doublets. Two of them are neutral and hence a smoking-gun signature of this scenario at the LHC consists of displaced Higgs bosons or Z-bosons plus missing momentum. We analyse the LHC sensitivity on such interesting final states by recasting an existing ATLAS search for displaced vertices and missing energy [53]. Making use of the large statistics already collected, we show that the proper decay length that can be (will be) tested is actually significantly larger than the detector size, reaching more than 10 (100) meters. In the corresponding parts of the parameter space with the correct dark matter abundance, this currently constrains the dark matter mass to be as large as 500 keV and could reach a few MeV with 300 fb −1 , hence well beyond the warm dark matter regime.
The rest of the paper is organised as follows. In section 2, we discuss the model in the feeble-coupling regime and analyse the decay modes of the dark sector particles. Section 3 is devoted to the calculation of the abundance of dark matter produced through the freezein mechanism. We then study the constraints on light dark matter from cosmology in section 4. In section 5, we analyse in detail the main collider signatures -disappearing charged tracks, as well as displaced Higgs or Z plus missing transverse momentum (E / T ) -of our scenario and we present the recasting of the ATLAS search [53]. In section 6, we combine the results of the previous sections and show the interplay between collider and cosmological signatures in probing our model of freeze-in dark matter. We summarise and conclude in section 7, while we present some technical details in the appendices.

The feebly-coupled Singlet-Doublet DM model
We perform our analysis within the Singlet-Doublet dark matter model [41], which consists in adding to the SM a pair of Weyl doublet fermions, ψ u and ψ d , with opposite hypercharges and one fermionic singlet, ψ s : The subscripts indicate the SU(2) L × U(1) Y quantum numbers and these new fields are assumed to be odd under an unbroken Z 2 symmetry, under which the SM fields are even, so as to guarantee the stability of dark matter. The mass terms and Yukawa interactions of the model read where H is the Higgs doublet (with hypercharge 1/2) and · indicates a contraction of the SU(2) L indices through the antisymmetric tensor ab , see appendix A for more details.
For later convenience, we also define the following alternative parameterisation of the two Yukawa couplings: y u ≡ y sin θ, y d ≡ y cos θ.

JHEP09(2018)037
As is apparent from the field content, the model is a generalisation of the "Bino-Higgsino" system of supersymmetric models with free couplings y u and y d (whereas supersymmetry would relate them to the SM gauge couplings).

The spectrum
Upon EW-symmetry breaking, the Lagrangian in eq. (2.2) leads to mixing among the neutral components (ψ s , ψ 0 d , ψ 0 u ) of the Z 2 -odd fermions. The resulting mass matrix reads where v 246 GeV is the vev of the Higgs field. The above matrix is diagonalised by a rotation matrix U , U MU T =M. The mass eigenstates are then given by We employ the convention |m χ 1 | < |m χ 2 | < |m χ 3 |, thus our dark matter candidate is χ 1 . The model has already been extensively investigated within the framework of the freeze-out mechanism of dark matter production [41][42][43][44][45][46][47][48][49][50][51][52]. The latter requires the couplings y u and y d to be of the order of 10 −2 − 1 for the Yukawa interactions to drive the relic abundance to the observed one. In that case the model is constrained by direct and indirect detection experiments and also features interesting collider signatures, typically resembling the Bino-Higgsino system of supersymmetric models (but with arbitrary Yukawa couplings). In contrast, in this work, we focus on the freeze-in mechanism of dark matter production associated to very feeble Singlet-Doublet interactions. As it will become clear from the discussion in section 3, the typical Yukawa couplings of interest for our analysis range from 10 −9 to 10 −6 and a large mass difference between the singlet and the doublet mass-scales will have to be considered: |m s | |µ|. As a result, the model features suppressed mixing between the singlet and doublet and the singlet fermion ψ s χ 1 is the lightest of the neutral fermions.
In the limit |y u |, |y d | 1, |m s | |µ|, (2.6) we can expand the mass eigenvalues at the first order in y 2 u,d and get 2 From these expressions, we see that, in the feebly-coupled regime, there is one neutral state of mass approximately m s (corresponding to χ 1 ∼ ψ s ), two neutral fermions χ 2,3 JHEP09(2018)037 with mass approximately µ (equal mixture of ψ 0 u and ψ 0 d ), and one charged fermion with mass m ψ = µ, that we denote with ψ ± .
The set of states with tree-level mass µ is further split by quantum corrections at one loop, which increases the mass of the charged state. Using the results of [54], one finds that the splitting between the charged and the neutral states is where α 2 = g 2 /(4π) (with g being the SU(2) L gauge coupling) and θ W is the weak mixing angle. Considering |µ| > 100 GeV, ∆M spans the following range 250 MeV ∆M 350 MeV . (2.9) A sketch of the spectrum of the model and the possible decay modes, described in the next section, is shown in figure 1.

Decay modes and decay lengths
In this subsection we study the decay modes of the fermion mass eigenstates in the feeble coupling regime that control the phenomenology of the model. General expressions for the decay widths through the model's Yukawa interactions can be found in the appendix of ref. [47], while here we specialise to the regime of (2.6) making use of the expression for the mixing matrix reported in appendix A. Throughout this work, we consider |µ| > m W as a doublet mass lower than about 90 GeV is excluded by searches for charged fermions performed at LEP (see [55], for a recent reassessment). In this regime, the heavy charged states ψ ± can decay directly into the lightest mass eigenstate, ψ ± → W ± χ 1 , via the suppressed Singlet-Doublet mixing, or to the heavier neutral states (∼ neutral components of the doublet) and a soft pion, ψ ± → π ± χ 2,3 , via gauge interactions. The latter decay mode occurs via an off-shell W and, for the range of charged-neutral state mass splitting reported in eq. (2.9), it dominates over possible leptonic modes involving ± ν instead of π ± . Half of the decays into pions JHEP09(2018)037 Figure 2. Left: branching ratios of the decay of the charged fermion into χ 2,3 plus a pion. Right: decay length (in cm) of the charged fermion; the shaded region is excluded by the ATLAS search for disappearing tracks [57,58]. The mass of the lightest neutral fermion is fixed to m χ1 = 10 keV. plus neutral states go to π ± χ 2 and the other half to π ± χ 3 , with the partial decay widths given by [56]: where f π 130 MeV, G F is the Fermi constant and θ c the Cabibbo angle. As mentioned above, this decay mode competes with the decay into the lightest fermion eigenstate W ± χ 1 induced by the small Singlet-Doublet mixing. At leading order in the Yukawa couplings, taking tan θ = 1 (i.e. y u = y d = y √ 2 ), the decay width reads (for µ ≷ 0): α is the electromagnetic constant, and s W ≡ sin θ W . It can be shown that Γ[ψ ± → W ± χ 1 ] has a negligible dependence on tan θ in the limit (2.6), and hence the formula (2.11) will suffice to our purposes. Making use of eqs. (2.10) and (2.11), we show in figure 2 the branching ratios for the decay process ψ ± → π ± χ 2,3 (left panel) and the contours for fixed values of the charged fermion decay length cτ ψ (right panel) on the (m ψ , y) plane in the custodial symmetry limit y u = y d for m χ 1 = 10 keV. We can see that for y ≤ 10 −8 the decay mode into π ± χ 2,3 is the dominant one. Moreover, comparing the two panels of figure 2, we can see that, when ψ ± decays preferably into pions plus heavy neutral state, the decay length is about 1 cm

JHEP09(2018)037
that is approximately the minimal length, to which LHC searches for disappearing charged tracks have sensitivity. For illustrative purposes, we thus also show, in the right panel of figure 2, a purple region excluded by a recent ATLAS analysis [57,58]. We will address the collider constraints in more details in section 5. Let us mention that the overall picture depends neither on m χ 1 nor on tan θ in the limit (2.6) that we are interested in here.
The two neutral states χ 2 and χ 3 can decay either into Zχ 1 or into hχ 1 through the Yukawa interactions. At leading order in y u and y d , the decay widths read (for µ ≷ 0): after having redefined the fermionic fields in order for m χ i to be positive. Notice that, in practice, in the regime of (2.6), the above decay widths are not affected by the sign of µ.
The typical decay lengths of χ 2 (solid line) and χ 3 (dashed line) are shown as a function of their mass in figure 3 for two values of tan θ and m s = 10 keV. As can be seen, for a y coupling of the order 10 −7 , the decay length is around 1 cm, while it exceeds 10 m for y ∼ 10 −9 . In addition, the decay lengths of χ 2 and χ 3 appear to be essentially equal for m χ 1,2 300 GeV or for tan θ 1. Let us also emphasize that, in the limit of eq. (2.6), the decays χ 3 → χ 2 are not allowed due to the tiny mass splittings, as it can be verified by inspecting the expressions in eq. (2.7).

Dark matter abundance from Freeze-in
In the feeble Yukawa coupling regime that we are considering, the dark matter candidate χ 1 ∼ ψ s has strongly suppressed interactions with the SM particles. Hence it is not in thermal equilibrium with the SM bath at the time of production. In contrast, the components of the electroweak doublet are in thermal equilibrium because of their unsuppressed gauge interactions. Assuming zero initial abundance of χ 1 , the dominant production mechanism for the dark matter particle is through the decay of the heavy mediators (ψ ± and χ 2,3 ) along the cosmological evolution. 3 This production "freezes-in" when the abundance of JHEP09(2018)037 Decay length of the two heavy neutral fermions, χ 2 (solid lines) and χ 3 (dashed lines) for two different choices of tan θ = y u /y d . The mass of the lightest neutral fermion is set to m χ1 = 10 keV. the heavy mediators is Boltzmann suppressed, that is approximately when the temperature drops below their mass. This is the framework in which we carry out our analysis; see e.g. [5-9, 11, 59, 60] for some previous examples. The DM comoving number density induced through the decay of A → B χ DM simply reduces to [5] where g A counts the spin degrees of freedom of the mother particle A, g * is the number of degrees of freedom at the freeze-in temperature T ∼ m A , and M Pl = 1.22 × 10 19 GeV is the Planck mass. This result is obtained making the following simplifying assumptions: (i) the mother particle A and the daughter particle B are in thermal equilibrium with the SM thermal bath; (ii) A follows a Maxwell-Boltzmann distribution function; (iii) we can neglect the Pauli-blocking/stimulated emission effects associated to B. See also the discussion in [61].
Let us notice that the freeze-in mechanism considered here -taking place through the decay of a mediator that is in thermal equilibrium with the Standard Model sector in the early universe -shares some similarities with the so-called superWIMP mechanism [62]. In the latter case, the late decay of the WIMP mother particle -occurring considerably after its freeze-out -give rise to the dark matter abundance. The main difference between the two mechanisms is the life-time of the mother particle in thermal contact with the Standard Model. In the superWIMP case, the life-time is typically much longer (τ ∼ 10 5 − 10 8 s) and the scenario is subject to constraints from big bang nucleosynthesis [62,63].

Freeze-in: collider and cosmology interplay
In the context of the Singlet-Doublet model, three mass degenerate heavy states can decay into dark matter, namely the charged fermion ψ ± and the two neutral fermions χ 2,3 , giving rise to the DM yield: 2) where g ψ = 2 takes into account the number of degrees of freedom of the charged fermion. Notice that the contributions of the heavy neutral fermions χ 2,3 directly depend on the total decay widths of χ 2,3 , which we will denote Γ χ 2,3 in the following, as the decays into χ 1 are the only available decay modes. On the other hand, for the charged fermion, only the partial decay width into the W + χ 1 final state appears since the ψ ± decays into χ 2,3 are already accounted for by the two first contributions associated to the thermal equilibrium abundances of χ 2,3 . We can now compute the DM relic density in terms of Y χ 1 : where the present entropy density and critical density are respectively s 0 = 2.8912×10 9 m −3 and ρ c = 10.537 h 2 GeV/m 3 . Considering that, in our scenario, we have m χ 2,3 m ψ µ, we obtain as a result where g A i is the number of degrees of freedom of the mother particle Obtaining the dark matter yield on more general grounds, starting from Maxwell-Boltzmann statistics, requires a fully numerical treatment of the evolution equations, which makes the computation and the interpretation of the freeze-in mechanism less straightforward. The authors of ref. [61] have however recently delivered the public code micrOMEGAs5.0 that allows to easily handle such computations. We have explicitly checked that the analytical results presented here are in excellent agreement with the ones obtained with micrOMEGAs5.0 (employing the Singlet-Doublet model files from our implementation in FeynRules [64]) in the Maxwell-Boltzmann limit. 4 Beyond this simplifying assumption, the full numerical treatment of the evolution equations gives rise to a moderate positive correction to Ω χ 1 h 2 with respect to our analytical result (about 25%). In this paper, we choose to discuss the results of our analysis with the Maxwell-Boltzmann approximation, thus neglecting the above small correction, in order to have a fully analytical understanding of the parameter space of the model yielding the observed DM abundance.

JHEP09(2018)037
The result of eq. (3.4) is rather generic for the freeze-in scenarios, independently of the underlying dark matter model, and indicates the typical order of magnitude of the physical quantities involved, that is the dark matter mass, the mediator(s) mass, and the mediator(s) widths. Note that the decay length of a particle is related to the total decay width through Hence mother particles or equivalently mediators, A, with a total decay width allowing for the dark matter density to be in agreement with the observed abundance Ωh 2 0.12, are in the right ballpark to give rise to macroscopically long displacements at colliders. In order for the mediators to be produced at the LHC, their mass cannot typically exceed the TeV scale. For such mass scale, a DM mass in the heavy range 1 MeV m χ 1 1 GeV corresponds to mediators escaping the detectors, while for light dark matter, 1 keV m χ 1 1 MeV, the signature will be characterised by displaced vertices visible inside the detectors. This highlights the natural interplay among LHC long-lived or displaced signatures, the freezein mechanism, and cosmological or astrophysical probes of light ( keV) dark matter. An important remark is that these considerations and correlations are strictly correct for mediators that can decay into dark matter only, which is the case of χ 2,3 here. We study this complementarity in further detail in sections 5 and 6.

The viable dark matter parameter space
We can now employ the decay widths that we computed in section 2.2 to derive the predicted value of the dark matter relic abundance on the parameter space of the Singlet-Doublet model. Expanding the expressions of the decay widths in the limit of small m χ 1 one finds at leading order: These expressions show that the combinations of the decay widths entering in the computation of the relic abundance do not depend on tan θ at zeroth order in m χ 1 . Plugging these expressions into eq. (3.3), we find the following approximate expression which accounts for the correct relic abundance up to a few percent level error when µ 400 GeV. Eq.  In the left panel of figure 4, we show the dependence of the dark matter abundance through the freeze-in mechanism on the parameters of the model for a fixed DM mass m χ 1 = 20 keV on the (m χ 2,3 , y) plane, or equivalently (µ, y) plane. It appears that, for m χ 1 ∼ few tens of keV, the coupling should be y = O(10 −8 ) in order to reproduce the observed dark matter density. In the right panel of figure 4 instead, we show on the (m χ 2,3 , m χ 1 ) plane contours of the values of the coupling y that yield Ω χ 1 h 2 = 0.12. We can see that, for m χ 1 in the [100 MeV, 1 keV] mass range and a µ scale of relevance for colliders, the required size of the Yukawa coupling is in the range 10 −8 y 10 −10 . (3.10) The largest values of the coupling y ∼ 10 −8 allow for very light dark matter candidates (few keV) to account for all the dark matter while heavier particles of hundreds of MeV requires even more suppressed Yukawa interactions with y ∼ 10 −10 . This observation will be relevant when comparing the reach of collider experiments to the one of cosmology probes on the frozen-in Singlet-Doublet dark matter parameter space.

Cosmology probe of light dark matter
Dark matter candidates with non negligible velocity dispersion deep in the radiation dominated era can leave a distinctive imprint in cosmology and astrophysics observations due to their free-streaming that delays the structure growth. Overdensities are suppressed below the comoving free-streaming horizon given by where a is the scale factor and H is the Hubble rate and v is the velocity dispersion of the dark matter ( v is given by the velocity of light for relativistic dark matter). For example,

JHEP09(2018)037
thermal warm dark matter (WDM), which was in thermal equilibrium and relativistic until decoupling at temperature T D , has a free-steaming length of λ f s Mpc (keV/m X ) T D /T ν where T ν is the temperature of active neutrinos and m X is the WDM mass. 5 Such a WDM scenario has served as a benchmark for non-cold dark matter cosmology. In this work we exploit the results of Lyman-α forest studies that have been used to set a lower bound on the thermal WDM mass of at 95% confidence level (CL) [20]. Notice that the above constraint can e.g. relax to m WDM > 3.2 keV at 95% CL when the WDM makes only part (> 80%) of the total DM content [65]. It has also been argued that considering a different thermal history in the treatment of the Lyman-α forest data (especially the ones associated to high redshift quasars), a few keV DM candidate could even provide a good fit to the data, see [19,65,66], a possibility that is strongly challenged by X-ray constraints in the context of sterile neutrinos [65]. Thermal warm dark matter is not the only relic that would suppress small scale structure formation. Other DM candidates with non negligible velocity at the time of production will give rise to similar effect. 6 Among them, one finds (non-)resonantly produced sterile neutrinos [78][79][80][81][82], sterile neutrinos from frozen-in scalars [83,84], mixed dark matter scenarios [85] and -of interest for this paper -other frozen-in particles [5,[86][87][88][89]. The different mechanisms of production involved can typically give rise to distribution functions that can differ from the (thermal) Fermi-Dirac distribution. As a result, the imprint on the linear matter power spectrum should a priori be recomputed making use of the relevant Boltzmann codes. Dedicated hydrodynamical simulations should then be performed so as to extract the non-linear evolution of a baryon+DM population and properly compute the observables relevant to estimate the Lyman-α flux power spectra within a given DM scenario and compare with data. All this procedure is however beyond the scope of this paper.
Here we use the constraints that have been derived in refs. [86,88] on keV dark matter produced through the freeze-in from the decay(s) of some thermalised mother particle A into the DM and another daughter particle B. In ref. [88], the suppression of the linear matter power spectrum in the freeze-in scenario has been computed and compared to the one of thermal WDM with a mass of 4.65 keV. 7 This provided a constraint on the mass of 5 The thermal abundance of WDM is given by ΩX h 2 (TD/Tν ) 3 mX /94 eV where TD/Tν = (g * (Tν )/g * (TD)) 1/3 for entropy conservation with g * (T ) the effective number of relativistic degrees of freedom and g * (Tν ) = 10.75. All together, an injection of large number of relativistic degrees (> 10 4 ) of freedom is needed, compared to the available ones in the SM g * (T > TEW ) = 106.75, so as to be able to get a thermal WDM of a still allowed few keV WDM [19][20][21]. 6 Notice that collisional or Silk damping (in contrast with free streaming ≡ collision-less damping) can also give rise to a similar imprint in small scale structures. This would be typically the case of dark matter interacting with relativistic species, see e.g. [67][68][69][70][71][72][73][74][75][76][77]. 7 In ref. [88], the transfer function of DM (associated to the ratio of cold DM and freeze-in DM linear matter power spectra) produced through freeze-in from the decays of some thermalised mother particle A → B+ DM always appear to have the very same spectral form as the one of thermal WDM. For other references, estimating the range of viable non-cold dark matter candidates based on the derivation of the linear matter power spectrum, see e.g. [71,85,87,[90][91][92].
JHEP09(2018)037 the frozen-in dark matter particle that shows a dependence on the mass splitting between A and B for compressed spectra. Within our framework, the resulting constraint on the dark matter mass reads: where η = 1.9 (as obtained in ref. [88] from a numerical fit), and Γ ij is the decay width 3 and ψ ± denoting the mediators that decay into the dark matter fermion χ 1 and another SM final state B j = Z, h or W ± . The results of ref. [88] imply thus that, in general, frozen-in DM, resulting from the decay of a thermalised mother particle with m DM > 12 keV, evades the constraints from the Lyman-α forest data derived in [20]. Lower DM masses can become allowed when ∆ ij is small, i.e. for small mass splittings between A i and B j (as pointed out in [86]). Our bound is shown in figure 5 as a function of the doublet mass scale µ. As we can see, the lower bound on the DM mass becomes weaker than m χ 1 > 12 keV only for values of the doublet mass µ approaching the mass of the decay products B j = Z, h, W ± .

Signatures at the LHC
In the feebly-coupled regime of the Singlet-Doublet dark matter model, the mediators ψ ± and χ 2,3 are essentially the charged and neutral components of the extra SU(2) L doublets. As a result, they can be produced at the LHC through electroweak processes. These JHEP09(2018)037 production processes are induced by gauge couplings only and thus the cross sections are independent of the couplings y u and y d . They are actually equal to those of a pure Higgsino in supersymmetry (SUSY) that can be computed using public tools such as Prospino2 [101].
Referring to the SUSY nomenclature, the relevant production modes include neutralino pair production, chargino pair production, and associated production of neutralino and chargino: pp → χ 2 χ 3 + X, pp → ψ + ψ − + X, pp → χ 2,3 ψ ± + X. Being substantially decoupled from the SM sector, the singlet dark matter χ 1 can only be produced at the last step of the decay chain, with the possible decay modes as illustrated in figure 1. We report in figure 6 the total production cross section (obtained by summing over all mediator pair and associated production modes) with a continuous red line, and the production cross section of a χ 2 χ 3 pair with a dashed blue line, as a function of the doublet mass scale µ. The cross sections were computed by means of Prospino2 for pp collisions with √ s = 13 TeV at next-to-leading order (NLO).
In order to obtain the collider constraints, we first have to compute the typical decay length of the heavy mediators in the viable dark matter parameter space, i.e. where Ω χ 1 h 2 = 0.12. In figure 7, we present the decay length of the mediators for the model parameters accounting for the whole observed dark matter abundance. In the left panel, we show contours for the decay length of the heavy neutral fermions, cτ χ 2,3 . On general grounds, the results depend on tan θ but, as already noticed in figure 3, the tan θ dependence is negligible as long as µ 300 GeV or tan θ 1. Also, as expected from the discussion in section 3, the figure shows that decay lengths leading to displaced signatures within the volume of LHC detectors correspond to the light dark matter regime, m χ 1 1 MeV. On the right panel of figure 7, the dashed orange contours indicate the branching fraction of the ψ ± decay into pions and χ 2,3 . We see that this decay mode is dominant except in a small corner of the parameter space where m χ 1 = O (1) keV and m ψ = µ is larger JHEP09(2018)037 Figure 7. Left: decay length of the neutral fermions χ 2 (solid lines) and χ 3 (dashed lines). Right: branching ratio (orange dashed lines) and decay length (green solid lines) of ψ ± → χ 2,3 π ± ; the shaded region is excluded by searches for disappearing tracks [57,58]. The coupling y is set in both plots by requiring Ωh 2 = 0.12.
than about 1 TeV. The ψ ± decay mode into W ± χ 1 , driven by the Yukawa interactions and contributing to the dark matter relic abundance, is typically subdominant, due to the feeble couplings involved. As a consequence, the decay length of ψ ± is always of the order of 1 cm in the parameter region relevant for the freeze-in mechanism as shown by the green solid lines in the right panel of figure 7.
We can now discuss the LHC signatures of our freeze-in Singlet-Doublet model.
• Disappearing tracks: independently of the final steps of the decay chain, the charged fermions ψ ± decay with a small displacement (of the order 1 cm at most, cf. figure 7) leading to 'disappearing' charged tracks that can be searched for at the LHC. In fact, a recent ATLAS analysis [57] (see also the similar search [102] from the CMS Collaboration), reinterpreted in [58] in terms of pure Higgsino production (which is exactly our case), excludes the regions shaded in purple in the right panels of figures 2 and 7. This search constrains the mass of the charged fermions to be larger than about 150 GeV in the regime in which ψ ± → π ± χ 2,3 dominates. For future prospects of searches for disappearing tracks and possible strategies to increase their sensitivity, see [103,104], where Higgsino masses up to approximately 400 − 500 GeV are foreseen to be accessible at the future high-luminosity run of the LHC (HL-LHC).
• Displaced h and/or Z + E / T : most of the mediator production modes will eventually produce a pair of heavy neutral fermions (χ 2 χ 2 , χ 2 χ 3 or χ 3 χ 3 ), possibly with extra soft objects that will go undetected. Indeed, as shown above, the relic abundance requirement implies that the charged fermions decay dominantly into the heavy neutral ones, χ 2,3 , plus soft pions. Given the possible decay modes of χ 2,3 , our key collider signature is thus characterised by a final state with displaced ZZ, hh or Zh, plus missing momentum,  Figure 8. Schematic representation of one of the processes leading to displaced Z or h bosons plus missing energy at the LHC. The red lines denote long-lived particles. Similar final states arise from χ 2 or χ 3 pair production. Note that we do not specify the production mechanism of the pair of neutral heavy fermions since it could be produced directly through electroweak processes or through the decay of the charged fermion.
Hence, the precise signal yield in each of the three channels ZZ, hh and hZ is determined by the branching fractions of the two neutral fermions χ 2,3 . In figure 9, we show the branching fractions of the χ 2,3 decays into hχ 1 (dashed line) and Zχ 1 (continuous line) as a function of tan θ and for several benchmark masses. As we can see, for tan θ ≈ 1, i.e. y u ≈ y d , one of the two heavy fermions decays predominantly into Z + χ 1 and the other one into h + χ 1 , independently of their mass. This leads to final states with a balanced sample of hh (25%), ZZ (25%) and hZ (50%). The same is true for tan θ 1 or tan θ 1 and when the mass of the neutral fermions is much larger than the Higgs mass (where effectively one has BR[χ 2,3 → Zχ 1 ] = BR[χ 2,3 → hχ 1 ] = 50%). The only configuration where there is not a balance in h and Z is when the mass of the neutral fermions is close to the Higgs mass. In this latter case, kinematics favor the decays into Z + χ 1 , and hence final states with ZZ + E / T are more probable.
In the next subsection, we will estimate the constraints on the three final states with displaced ZZ + E / T , hh + E / T or hZ + E / T that can be obtained from existing LHC searches at 13 TeV for displaced signatures, and we will subsequently study the impact on the parameter space of our model. Notice that searches performed at the LHC with √ s = 8 TeV can also be sensitive to the main signatures of our model, ZZ + E / T , hh + E / T or hZ + E / T . A number of such searches have been considered in ref. [105] and reinterpreted in terms of supersymmetric models. In particular, our scenario is similar to the case of Higgsinos decaying into gravitino in gauge-mediated SUSY models considered in [105], which is constrained mainly by a search for displaced dileptons [106] and a search for displaced jet pairs [107,108], both performed by CMS. In the next subsection, we show a comparison of the sensitivity of these searches with the 13 TeV analysis we are going to recast. Note that our model and the Higgsino-gravitino scenario considered in [105] JHEP09(2018)037 In contrast, in the Singlet-Doublet freeze-in model, the mass splitting between the heavy neutral components is so small -as shown by the expressions in eq. (2.7) -that the two neutral fermions always decay directly to χ 1 Z or χ 1 h with branching ratios as illustrated in figure 9.
• Searches for prompt decays: for small values of the decay length of the mediators (corresponding to moderate/large values of y), we expect that limits from standard prompt searches can be effective. A combination of recent searches at the LHC with √ s = 13 TeV for production of supersymmetric charginos and neutralinos can be found in ref. [109]. Possible final states are hh + E / T , ZZ + E / T and hZ + E / T , which are typical signatures of searches targeting Higgsino-like neutralinos in gauge mediated supersymmetry breaking [110], or Higgsinos decaying into light Bino, see e.g. [111]. In these final states, limits to the Higgsino mass up to 600-700 GeV were obtained. 8  Note that the configurations of the model giving rise to the observed relic abundance through freeze-in considered in section 3 never give rise to prompt decays, i.e. the decay length is always larger than about 1 mm. We discuss thus these prompt decay searches briefly. For more detailed discussions of the prompt signatures of the Singlet-Doublet model, see refs. [43,47].
• Mono-X searches: in the region with very large decay length, where the neutral fermions escape the detectors, mono-X searches could be the only strategy to look for this model at the LHC (besides the disappearing charged tracks associated to the charged fermion). In this regime the collider signature of our model is very similar to an Higgsino dark matter scenario in which the mass splitting among the Higgsino components is tiny, as already mentioned. There have been several investigations on this scenario and the corresponding mono-X signatures, e.g. in [112][113][114][115][116][117]. Some of these investigations have exploited the soft leptons that would be present for a mass splitting of the order of few GeV, which is however not the case of our model. Instead, the case of a pure mono-jet signal has been shown to be not promising, with an estimated reach on the Higgsino mass of order 200 GeV at HL-LHC with 3000 fb −1 [115]. Hence we decide not to include these signatures in our analysis.

Recasting strategy for displaced h and Z + E / T
The aim of this subsection is to estimate the current LHC limit on displaced ZZ, Zh and hh + E / T using public information from ATLAS and CMS searches. In table 1, we report on the relevant ATLAS and CMS searches for displaced signatures, focussing on the most recent analyses at 13 TeV. Among these searches, we identify the recent ATLAS analysis on displaced vertices (DV) with jets and E / T [53] as the most promising for our scenario. The motivation is manifold: (i) this analysis exploits the largest available dataset among those listed in table 1; (ii) the large hadronic branching fractions of h and Z imply that our model yields a sizeable production cross section in this channel; (iii) our final states contain a relevant source of E / T , and the analysis of ref. [53] is the only one targeting it with a dedicated selection; (iv) and finally, detailed auxiliary material is provided with the information needed for a recasting [122].

JHEP09(2018)037
The ATLAS DV+E / T analysis [53] targets final states with at least one displaced vertex with jets and large missing transverse momentum. The results are interpreted in a model with long-lived gluinos decaying into jets and the lightest neutralino. In the auxiliary material [122], the efficiencies for the missing momentum and the displaced vertex reconstruction are provided. In particular, the efficiencies of the displaced vertex reconstruction are given prior to detector simulation, as a function of the invariant mass of the vertex, of the number of tracks and of the displacement.
In order to estimate the sensitivity of this search on our final states, we have first implemented the model in FeynRules [64], and then simulated the relevant samples with MadGraph5 [123], combined with Pythia8 [124] for the parton showering and the underlying pp collision, and Delphes3 [125] (with the standard ATLAS card) for the detector simulation. The displacement is applied to the simulated events a posteriori, taking into account the four momenta of the long-lived particle in order to properly compute the displacement, 9 which is obtained by sampling an exponential distribution with mean decay length cτ χ 2,3 . In appendix B, we discuss the details of the selection and the validation of our implementation. The latter was performed by reproducing the exclusion limits set by the ATLAS search on the simplified model they considered with long-lived gluino.
After this validation, we can now estimate the efficiency of the ATLAS DV+E / T analysis in the final states we are interested in, which are ZZ + E / T , hh + E / T or hZ + E / T . 10 We do this as a function of the lifetime τ χ 2,3 and of the mass of the long-lived particles, which are simply the two neutral fermions χ 2 and χ 3 produced in pairs (cf. appendix B for plots displaying the resulting efficiencies). The mass is important in order to determine the boost factor in the displacement, as well as to get the correct p T distribution of the displaced tracks. We can now use the obtained selection efficiencies to evaluate the reach of the ATLAS analysis in three simplified models with fixed branching fractions, that serve for illustrative purposes: In order to constrain the above simplified models, we consider the total production cross section of the doublet fermions states, computed at NLO by Prospino2 [101], summing all production modes shown in (5.1), corresponding to the solid red line in figure 6. With no background in the signal region, the parameter configuration of a model is excluded at 95% confidence level (CL) or more if it yields a number of selected events ≥ 3.0. The resulting estimated exclusion is depicted in figure 10 by the three solid lines. As we can see, the difference in the efficiencies among the three simplified models results in only a 9 In this approach we neglected possible distortions of the kinematic distributions of the final state charged tracks due to the displacement. 10 If Z or h decay into bb, (some of) the resulting tracks may have an additional displacement, which makes the reconstruction of the DV more involved, as discussed in more detail in appendix B. In the following, we neglect possible issues related to this for the reasons discussed in the appendix. small impact on the sensitivity. Moreover, the largest doublet mass (about 1.3 TeV) is probed for a decay length around cτ ≈ 5 cm. Also, the exclusion curves are not symmetric in cτ with respect to this maximal reach. This is due to the fact that the exponential distribution determining the displacement is falling very rapidly for a displacement larger than a given cτ , while it goes to zero less steeply for displacement smaller than cτ . This also explains why the reach of the analysis extends to regions with very large decay lengths, up to cτ ≈ 50 m.

JHEP09(2018)037
In figure 10, we also show for comparison the exclusion from the 8 TeV searches, as reported by ref. [105]. This is depicted as a dashed orange line and includes both searches targeting displaced leptons [106] and displaced di-jets [107,108]. The displayed 8 TeV limit has been obtained in ref. [105] in a simplified model with an Higgsino-like neutralino undergoing displaced decays into gravitino plus Z or h in the large tan β regime, which roughly corresponds to our simplified model ii). 11 As we can see, in the region of low doublet mass, the sensitivity of the ATLAS DV+E / T analysis is diminished because the spectrum is compressed and jet/E / T cuts become more severe. This is where the 8 TeV searches, in particular the one targeting displaced dileptons, become instead more efficient,

JHEP09(2018)037
despite the small leptonic branching fractions of the bosons. 12 Finally, we also display on the same plot an estimate of the reach of the prompt searches (as a purple dot-dashed line), considering the BR[χ 2,3 → Zχ 1 ] = BR[χ 2,3 → hχ 1 ] = 50% case reported in [109]. We stress that this limit will never be relevant in the parameter region leading to the correct freeze-in dark matter abundance, but we report it here for illustrative purposes. In order to draw this line, we have compared the cross section limits reported in ref. [109] to the total production cross section of the doublets multiplied by the probability that both produced particles decay promptly given a certain mean decay length cτ χ 2,3 . 13

DV + E / T constraints on the Singlet-Doublet model
We can now use the recasting presented above to provide estimates for the ATLAS exclusion on the parameter space of the Singlet-Doublet freeze-in model. For this purpose, at each point of the parameter space, we sum the production cross sections over the production channels weighted by the appropriate branching fraction in order to determine the signal strength for each of the possible final states. For instance, the signal cross section in hh + E / T is given by and analogous expressions can be written for the hZ + E / T and ZZ + E / T . The production cross section in each channel is hence a function of the parameter space of the model through the branching fraction dependence on (y u , y d , µ, m s ). We multiply these three type of signal cross sections with the corresponding efficiencies (derived in appendix B) to obtain the final estimate on the number of expected events. Each efficiency is also a function of the parameters (y u , y d , µ, m s ), since it depends on the mass of the long-lived particle, which is simply µ, and on the decay lengths that follow from eqs. (2.12)-(2.15). For simplicity, we take the average of the decay length of χ 2 and χ 3 as the mean decay length setting the displacement. As we have discussed above, this is an excellent approximation as long as tan θ 1 or µ 300 GeV (see figure 3). We neglect the extra displacement induced by the decay of the charged fermion. Note that this is indeed typically a small fraction of the overall displacement in the relevant portion of the parameter space, as illustrated by the green contours in the right panel of figure 7. In our estimate we also consider the same efficiency in the case in which the neutral heavy fermions are directly pair produced as in the case in which the neutral fermions are produced through the decay of the charged fermion. We checked this hypothesis on a few benchmark points and it induces an effect of at most 20%, which is largely negligible for the purpose of our recasting. As for the 12 We also remark that, as discussed in detail in the appendix B, our implementation of the ATLAS DV+E / T analysis tends to overestimate the exclusion in the compressed region (for mass splittings 100 GeV). The complementarity with the 8 TeV searches is thus welcome. 13 As a rough estimate, we consider to be prompt the events with a total displacement ≤ 0.5 mm that we compute based on cτ only, without taking into account the boost factor. Under the above assumptions, we can assess the current limits on the Singlet-Doublet parameter space from the ATLAS DV+E / T search. The region excluded according to our recasting is shown with filled cyan colour in figure 11. Its shape follows from combining the excluded regions for the simplified models reported in figure 10 with the iso-contours of the average cτ χ 2,3 (denoted as cτ ). The latter are shown with green continuous lines in figure 11 while the dashed cyan curve gives the estimated reach of an analogous DV+E / T search with a dataset of 300 fb −1 . 14 The red continuous curve shows the (y, µ) combinations that account for all the DM for a 12 keV DM candidate. Going above the red line, i.e. to larger values of the coupling y, induces an overabundant dark matter population, while below the red line it is underabundant; see eq. (3.9). Finally, the dot-dashed line delimits the region excluded by LHC searches for the prompt signature W Z + E / T .

JHEP09(2018)037
It is remarkable that the sensitivity of the ATLAS DV+E / T search extends to heavy electroweak states and to quite large values of the decay lengths. This is related to the almost background-free nature of displaced vertices signatures which renders the search very efficient even for small signal cross section. Note that the largest excluded mediator mass is about 1.2 TeV, somewhat smaller than for the simplified model analysis of figure 10. The reason is as follows. In the high mass region, when the lifetime, or equivalently the coupling y, maximises the experimental sensitivity (y ≈ 10 −8 ) the branching fraction of the process ψ ± → χ 1 W ± is not completely negligible (up to about 10%, see figure 2 left), and hence the signal yield into long-lived neutral fermion pairs is slightly diminished. In the case of the DV+E / T analysis extrapolation to 300 fb −1 , we expect to probe masses of the neutral fermions up to 1.7 TeV and decay lengths as large as 100 m.
Let us add here a remark on the uncertainties on our recasting and their effects on the estimated limits in figure 11. Given the steep fall of the production cross section as a function of the mediators' mass (see figure 6), we expect that even O (1) modifications of our estimated efficiencies would have a small impact on the mass reach (for instance a 50% change in the efficiency would only correspond to a change of around 10% in the mediators' mass limits). 15 Let us stress that the collider bounds presented in figure 11 are expected to be independent of the dark matter mass for m χ 1 below the GeV scale. The only curve that is affected by the m χ 1 parameter is the relic abundance continuous red contour. Considering larger values of the dark matter mass the red line would be shifted to lower values of y. As a result, the dark matter candidates with m χ > 12 keV (i.e. compatible with the Lyman-α bound discussed in section 4) are concerned with the excluded region below the red curve of figure 11. For instance, from the right panel of figure 4, one can deduce that for e.g. m χ 1 ≈ 1 MeV the Ωh 2 = 0.12 contour should appear at y ≈ 10 −9 in figure 11. This corresponds to larger values of cτ where the DV+E / T search loses sensitivity in such a way that no constraint can be set at present.
For completeness, let us briefly discuss the prompt decay constraints. In the upper part of figure 11, the size of the coupling y is such that the mediators decays are prompt. In particular, the charged fermion ψ ± predominantly decays into W ± χ 1 . 16 In order to estimate the corresponding constraint, we have computed the W Z + E / T production cross section in the Singlet-Doublet model multiplied by the probability that both heavy particles decay promptly, using the same approximations as for the prompt exclusion in figure 10. Comparing the latter results with the limits on the cross section given in ref. [109] we exclude the region delimited by the dot dashed purple line of figure 11. As discussed above, such constraint lies however in a zone of the parameter space where the frozen-in dark matter scenarios with masses above the Lyman-α bound give an overabundant dark matter relic density. 15 Note in particular that this applies to the possible issues associated with b-jets, discussed at the end of appendix B, that would at most reduce the signal strength by ≈ 25%. 16 The other prompt signatures discussed above are less sensitive as the production cross section is sensibly lower for χ2, χ3 production only. The latter is indeed almost one order of magnitude smaller than the total doublet production as seen in figure 6. Our estimate of the ATLAS DV + E / T exclusion is shaded in cyan ("DV+MET"), the magenta region is excluded by disappearing tracks ("DT"), the Lyman-α bound is shown in gray ("Ly-α").

JHEP09(2018)037
Green contours correspond to the average χ 2,3 decay length. The cyan dashed line is the estimated exclusion of LHC with 300 fb −1 . The coupling y is fixed such that Ω χ1 h 2 = 0.12 everywhere.

Displaced vertices vs cosmology for freeze-in DM
We can now combine the LHC limits and the cosmological bound derived in the previous sections, in order to characterise the experimental sensitivity on the viable parameter space of the freeze-in Singlet-Doublet model. As at the end of section 3, we present our results in the DM mass vs mediator mass plane fixing in each point the coupling y to the value that accounts for the observed relic abundance through the freeze-in mechanism. On the same two dimensional plane, we can show the combination of the existing (and future) constraints on the model. Our summary plot is shown in figure 12. As before, the green lines indicate contours of fixed average decay length of the neutral fermions, which controls the phenomenology at colliders. The magenta shading at low mediator masses represents the region excluded by searches for disappearing charged tracks (DT) [57,58]. It does not depend on the DM mass (or equivalently on the value of the y coupling) since, in this region, the decay length of the charged fermion is independent of m χ 1 , as can be seen in the right panel of figure 7. The cyan region and the dashed cyan line are the estimated exclusion and future prospect of the ATLAS DV+E / T search, discussed in section 5. The gray region, finally, is excluded by the Lyman-α forest data, as discussed in section 4. Figure 12 summarises the findings of this paper, as it nicely shows the interplay between collider searches for displaced signatures and cosmological constraints in our freeze-in dark matter model. On the one side, the observed DM abundance implies a relation among the parameters of the theory, leaving only two free parameters (plus a third one, tan θ, that affects the phenomenology of the model very mildly in our limit, as we discussed in JHEP09(2018)037 the previous sections). On the other side, for the range of decay lengths that are a priori optimal for studying displaced signatures at the LHC (O(10) cm) and µ scales within the reach of the collider (µ O(1) TeV), our dark matter model can leave a testable imprint on small scale structures. In such a region, we have a complementary constraint from the Lyman-α forest observations, which is essentially independent of the mediator mass. In contrast, the reach of LHC searches is intrinsically limited by the production cross section of the mediators, hence by their mass, and by the size of the detector. Due to the very low background of the recast search and the large dataset available, the current LHC limit actually extends to rather large values of the mediators' lifetime and, likewise, it can probe DM masses larger than those tested by cosmology, reaching up to m χ 1 = O(1) MeV.

Summary and conclusions
Despite many experimental and theoretical efforts, the nature of dark matter remains a mystery. It is thus timely to look for DM beyond the most popular paradigms. In this work, we considered the case of a dark matter candidate with such a tiny coupling that it never reaches thermal equilibrium with the SM in the early universe. It is well known that, despite such suppressed interactions, the observed dark matter density can be accounted for by the freeze-in mechanism, with the dark matter being produced, for instance, via the decay of thermalised mediators. Within this context, we studied the case of the Singlet-Doublet dark matter model that consists in extending the SM with a pair of Weyl electroweak doublet fermions and a singlet Majorana fermion. The new fermions interact with the SM through gauge interactions and/or the Higgs portal induced by Yukawa interaction terms that couple the doublet and the singlet fermions to the Higgs particle. The Singlet-Doublet model rests thus on 4 free parameters only: 2 new mass scales (the doublet mass scale µ and the singlet mass scale m s ), and two Yukawa couplings (y u and y d ). We have shown that, considering these couplings in the range [10 −8 , 10 −10 ] together with a doublet mass scale µ larger than the Higgs mass, the lightest neutral fermion, which is essentially the singlet Majorana fermion, can account for the whole dark matter abundance via the freeze-in mechanism. In this regime, the DM is light with a mass between a few keV up to hundreds of MeV.
Such a dark matter scenario could seem hopelessly beyond the reach of any dark matter experimental search. We show instead that the range of model parameters required for a successful freeze-in naturally gives rise to long-lived/displaced collider signatures that are already strongly bounded by the present LHC data. In addition, it is well know that thermal warm dark matter candidates of a few keV are also constrained by cosmology due to their free-streaming suppressing the growth of small scale structure. Even though frozen-in dark matter was never in thermal equilibrium in the early universe, Lyman-α bounds turn out to constrain the dark matter to be always heavier than 12 keV. In the low dark matter mass region, the model features thus both exotic LHC signals and a testable imprint on cosmology providing two complementary handles to probe the same scenario.
Concerning the collider searches, the relevant signatures of this model consist of disappearing charged tracks, related to the production of the charged component of the doublet JHEP09(2018)037 ψ ± , and displaced h and/or Z + E / T , associated to the decay of the two neutral fermions χ 2,3 . In the first case, ψ ± decays with a small displacement (about 1 cm) and our scenario is essentially the case of pure Higgsino DM in supersymmetric models. Current searches for disappearing tracks thus constrain the doublet mass scale µ to be larger than 150 GeV. In the second case, the χ 2,3 fermions decay with displacements in a wide range, from centimetres to kilometres, depending on the point of the viable parameter space of interest (i.e. for y u , y d giving rise to the right relic abundance through freeze-in). We have argued that, at present, the most constraining search was provided by ATLAS in ref. [53] and we have reinterpreted its results in the framework of the Singlet-Doublet dark matter model. According to our recasting, this analysis can exclude scenarios with a decay length of the heavy neutral mediators as large as ∼ 50 m, mediator masses as large as 1.2 TeV, and dark matter candidates with masses as large as 500 keV.
In figure 12, we have brought together all the experimental signatures which can probe the viable parameter space where the freeze-in production mechanism gives rise to the correct dark matter abundance. This nicely illustrates the interplay between collider searches and cosmology for frozen-in dark matter.
An interesting extension of our work is to enlarge the experimental reach on the parameter space of the model. The LHC sensitivity, shown in figure 12, could be improved towards large mediator masses, or towards small or large decay lengths. This is possible on all these three fronts by exploiting the presence of a displaced Z or h resonance, both in hadronic and leptonic decay channels, such that some of the event selection requirements that currently limit selection efficiencies can be relaxed, while keeping backgrounds to a negligible level. Also, at higher luminosities a dedicated event selection would help to suppress the increasing backgrounds. As a result of our study, we thus advocate dedicated experimental searches for displaced Z + E / T or h + E / T signatures, potentially in association with an extra identified Z or h boson. On the other hand, it would be interesting to also probe the case with large/moderate dark matter mass and very long-lived mediators (upper part of figure 12). For this purpose, one could for instance estimate the reach of the proposed detector MATHUSLA [126] on this scenario.
Finally, we stress again that interplay between exotic collider signatures and cosmology constraints go beyond the Singlet-Doublet model and apply to a large class of simplified models of freeze-in dark matter where the production occurs through decays of heavy mediators in thermal equilibrium with the SM bath. From the model building perspective, it would be interesting to investigate such complementarity in other models, also including those where the freeze-in is not realised through the decays of heavy mediators, but via scattering processes and/or via non-renormalisable interactions [127]. We leave these interesting possibilities for future works.

JHEP09(2018)037
by the Strategic Research Program High-Energy Physics and the Research Council of the Vrije Universiteit Brussel. AM and SL are supported by FWO under the "Excellence of Science -EOS" -be.h project n.30820817.

A.1 Yukawa interactions
For the contraction of the indeces in the Lagrangian of eq. (2.2), we follow the conventions of refs. [47,50]. In particular, the Yukawa interactions can be explicitly written as where i and j are SU(2) L indices.

A.2 Approximate expression of the rotation matrix
In the limit (2.6), the mass eigenvalues at the first order in the couplings y 2 u,d result as shown in eq. (2.7). We report here the rotation matrix, defined as in eq. (2.5), at leading order in y u and y d : We omit the O(y 2 u,d ) terms that are needed in order to diagonalise correctly the mass matrix M obtaining the eigenvalues shown in eq. (2.7). In fact, in our parameter regime, the above expression suffices to reproduce the rotations resulting from numerical diagonalisation to high accuracy. Hence, we employ it to derive the expressions for the decay widths of the heavy particles reported in section 2.2.

B Recasting of the ATLAS search
In this appendix, we provide details about the recasting of the ATLAS search of ref. [53] that we employed in order to set limits on the Singlet-Doublet model and in particular on displaced neutral bosons + E / T final states.
The signature that we consider both in the validation (gluino-neutralino simplified model) and in the Singlet-Doublet model is constituted by a pair of heavy long-lived particles decaying into charged tracks plus missing energy. For the case of two neutral heavy fermions, the process is depicted in figure 8, with the long-lived particle highlighted in red.
We first review the selection cuts of the search, we then validate our simulation with the simplified model studied in the ATLAS analysis, and then we apply the same recasting to our dark matter model.

JHEP09(2018)037
Selection criteria of the ATLAS DV+E / T search [53]. The ATLAS analysis [53] is explained in detail in the auxiliary material in [122]. The search targets displaced vertices and missing transverse momentum.
The displaced vertices are identified by analysing the associated displaced tracks. First, a selected displaced track should satisfy the following requirements: • The track is associated to a stable particle; • The particle has a transverse momenta p T > 1 GeV; • The transverse impact parameter d 0 ≡ R decay sin ∆φ > 2 mm, where R decay is the transverse decay length and ∆φ is the azimuthal angle between the heavy decaying particle momentum and the track momentum.
With the following selected tracks, one can construct a candidate displaced vertex which should satisfy the following criteria: • The transverse displacement R decay should be within 4 and 300 mm; • The longitudinal displacement should be smaller than |z decay | < 300 mm; • The number of associated charged tracks should be n tracks ≥ 5; • The invariant mass of the vertex should be larger than 10 GeV (a pion mass for the tracks is assumed).
Given the previous strategy to select displaced tracks and displaced vertices, events are hence required to satisfy the following conditions. These jets should satisfy the requirement that the scalar sum of the p T of the charged particles that are not displaced (according to the previous selection) should not exceed 5 GeV.
3. The events must contain at least one displaced vertex which has passed the selection. JHEP09(2018)037 Figure 13.
Comparison between the result of our simulation and the excluded cross section reported in the ATLAS paper [53] for a long-lived gluino simplified model. The ATLAS results are shown as red lines, while our analysis corresponds to the blue bands. In order to draw our bands, we considered variation of the efficiency of ±50%. The upper plots show the excluded cross section as a function of cτ for a fixed neutralino mass m χ = 100 GeV and two benchmark values for the gluino mass, mg = 1400, 2000 GeV. The lower plots show the excluded cross section as a function of m χ fixing τ = 1 ns and the same two benchmarks for mg.
Recasting and validation. On the selected events, one can then apply the efficiency as reported in the auxiliary material. Indeed, the ATLAS collaboration provides the efficiency for reconstructing the displaced vertices as a function of the number of displaced tracks and of their invariant mass. They also provides the efficiency tables as a function of the missing energy at truth level.
In order to recast this analysis, we have simulated LO events of the new physics process with MadGraph5 + Pythia8 + Delphes3 with standard minimal cuts, default parameters for MC and detector simulator (we used the default ATLAS card), assuming prompt decays of the heavy pair-produced particles. We employed generator level information in order to extract the momenta of the two heavy particles and of their associated tracks and in order to introduce the displacement by hand, including the boost factor of the heavy decaying particle. The decay time was generated through an exponential distribution with a mean lifetime τ . With this information we derived the impact parameter of each track and the other relevant geometrical properties. Then we processed the output following the ATLAS selection cuts strategy, including the reconstruction efficiencies. We first applied this procedure to a simplified model analogous to the one considered in the ATLAS paper, by considering gluino pair production followed by displaced decays into qq plus neutralino. Figure 14. Efficiencies for the simplified final states with hh + E / T , hZ + E / T , ZZ + E / T . Left: efficiencies as a function of cτ for m χ2 = m χ3 = 300, 500, 1000 GeV. Right: efficiencies as a function of the mass m χ2 = m χ3 for cτ = 1, 10, 100 cm.

JHEP09(2018)037
In figure 13 we show our estimated cross section exclusion compared to the ATLAS results. We find a good agreement in the region of un-compressed spectrum, while our simulation overestimate the exclusion power in the compressed spectrum case. We argue that this is due to our implementation of the jet cut (number 2. in the list above), for which we have only a limited amount of information provided by the ATLAS documentation. The plots of figure 13 shows that, on the other hand, our simulation consistently reproduces the ATLAS analysis for a mass difference between the heavy particle and its decay products larger than ≈ 100 GeV, taking into account an uncertainty of ±50% on the efficiency.
In figure 14 we display the result of our recasting: the efficiency curves for the final states characterising the Singlet-Doublet model. The samples have been generated at LO with MadGraph5 + Pythia3 + Delphes3 after implementing the model in FeynRules. We JHEP09(2018)037 simulated the following three cases pp → χ 2 χ 3 → ZZχ 1 χ 1 , pp → χ 2 χ 3 → hhχ 1 χ 1 , pp → χ 2 χ 3 → hZχ 1 χ 1 , where the decay of the bosons is performed in Pythia. We then processed the output with the selection procedures explained above to extract the efficiencies as a function of the mean lifetime cτ and the mass of the decaying heavy particles m χ 2 = m χ 3 . As a final remark, let us notice that the efficiencies displayed in figure 14 have been obtained by treating displaced heavy flavour jets like light-flavour ones. However, the case of Z or h decaying into bb pairs requires in principle additional care, since the recast DV+E / T search associates tracks to a displaced vertex based on track-vertex compatibility requirements, and merges displaced vertices if within 1 mm. For displaced b jets, these requirements are difficult to recast. As we mentioned, we choose to neglect this possible issue (an interesting discussion of which can be found in ref. [128]). We argue that this simplification has a limited impact on our estimated exclusions (shown in figures 11 and 12) for the following reasons: (i) Only one DV is enough to satisfy the analysis' requirements, thus there is no loss of sensitivity if at least one of the pair-produced heavy particles decays into a Z decaying into light flavours (and the DV is reconstructed); (ii) Due to the gluons radiated by the b quarks, a DV can still be formed on the tracks not coming from the b-decay vertex; (iii) Part of the b decays will still happen within the required 1 mm.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.