Probing dark matter freeze-in with long-lived particle signatures: MATHUSLA, HL-LHC and FCC-hh

Collider searches for long-lived particles yield a promising avenue to probe the freeze-in production of Dark Matter via the decay of a parent particle. We analyze the prospects of probing the parameter space of Dark Matter freeze-in from the decay of neutral parent particles at the LHC and beyond, taking as a case study a freeze-in Dark Matter scenario via the Standard Model Higgs. We obtain the projected sensitivity of the proposed MATHUSLA surface detector (for MATHUSLA100 and MATHUSLA200 configurations) for long-lived particle searches to the freeze-in Dark Matter parameter space, and study its complementarity to searches by ATLAS and CMS at HL-LHC, as well as the interplay with constraints from Cosmology: Big-Bang Nucleosynthesis and Lyman-α forest observations. We then analyze the improvement in sensitivity that would come from a forward detector within a future 100 TeV pp-collider. In addition, we discuss several technical aspects of the present Dark Matter freeze-in scenario: the role of the electroweak phase transition; the inclusion of thermal masses, which have been previously disregarded in freeze-in from decay studies; the impact of 2 → 2 scattering processes on the Dark Matter relic abundance; and the interplay between freeze-in and super-WIMP Dark Matter production mechanisms.


Introduction
While thermal freeze-out has long been considered the paradigm for Dark Matter (DM) production in the early Universe, experimental efforts over the last 30 years have found no conclusive evidence of the existence of this type of DM. At the same time, searches for DM in direct detection experiments [1,2], via indirect detection [3,4] and through collider probes (e.g. at the LHC [5,6]) yield at present strong bounds on the interaction strength between DM and the Standard Model (SM) particles.

JHEP03(2020)022
In freeze-in scenarios, DM particles are very feebly coupled to the thermal bath in the early Universe and never achieve thermal equilibrium, yet the coupling between DM and the thermal bath particles allows DM to be produced in decays and/or scatterings of bath particles. 1 Through these processes, the DM abundance slowly increases towards equilibrium, however without ever reaching it. For renormalizable interactions between DM and the thermal bath, the production of DM is most efficient for temperatures T around max(m A , m DM ), m A being the mass of the particle A whose scatterings/decays produce DM. 2 These processes shut-off (DM "freezes-in") soon after T drops below m A or m DM .
Given the feeble interactions between DM and the visible sector, freeze-in DM candidates are naturally compatible with the current null results from DM experimental searches, and probing these scenarios may be more challenging than that of thermal DM candidates. 3 Nevertheless, when freeze-in DM production in the early Universe proceeds via the decay of thermal bath particles [8], the corresponding feeble coupling would make these particles long-lived. Searches for long-lived particles (LLPs) at the LHC and beyond (see [27,28]) have received renewed attention in recent years from the lack of conclusive signals in prompt LHC searches for physics beyond the SM, and provide a promising avenue for probing freeze-in DM scenarios [29][30][31][32][33][34].
In this work we analyze the prospects of probing the parameter space of DM freeze-in from the decay of long-lived neutral parent particles at the LHC and future colliders. As opposed to electromagnetic/QCD charged LLPs, which leave tracks in the LHC detectors and could be discovered via searches for heavy stable charged particles (see [33] for an analysis of DM freeze-in LHC signatures in these scenarios), neutral LLPs leave no trace in the ATLAS and CMS detectors until their decay. They constitute an ideal search objective for the proposed MATHUSLA surface detector [28,[35][36][37] and other recent proposals for LLP detectors [38][39][40][41].
As a concrete model, we study a freeze-in DM scenario via the SM Higgs (a simple version of a supersymmetric Higgsino-Axino [29,42,43] or Higgsino-singlino system, very similar to the so-called singlet-doublet DM scenario from [32,44]). We perform an analysis of the projected sensitivity of MATHUSLA in its 100 m × 100 m (MATHUSLA100) and 200 m × 200 m (MATHUSLA200) design proposals [37], to the freeze-in DM parameter space, and study its complementarity to LLP searches with the ATLAS and CMS detectors which use displaced vertices and missing transverse energy E miss T [32]. We then analyze the improvement in sensitivity than would come from a forward detector within a future 100 TeV pp -collider (FCC-hh).
In addition, we discuss in detail several technical aspects of the present DM freeze-in scenario (the latter three being in fact more general): the role of the breaking of electroweak (EW) symmetry in the early Universe (the EW phase transition) on the freeze-in production of DM; the inclusion of thermal masses on the computation of the DM relic density, previously disregarded in DM freeze-in studies (there have been however very recent JHEP03(2020)022 efforts to include these effects as part of fully-consistent DM freeze-in analyses [17][18][19]); the potential impact of 2 → 2 scattering processes on setting the DM relic abundance and the interplay between freeze-in and super-WIMP DM production mechanisms (see also [15] for a detailed analysis on the latter). Finally, we analyze the constraints from cosmology, namely from Big-Bang Nucleosynthesis (BBN) and from limits on the washout of small-scale structure via Lyman-α forest observations.
Our work is organized as follows: we give a review of the freeze-in DM production mechanism in section 2. In section 3 we discuss DM freeze-in via the SM Higgs, introduce the model and analyze the DM relic density including the impact of the EW phase transition, of including thermal masses in the calculation and of the super-WIMP DM production mechanism. In section 4 we discuss the constraints from BBN and from Lyman-α forest observations. In section 5 we obtain the projected sensitivity of MATHUSLA to the DM freeze-in parameter space, and in section 6 we analyze its complementarity with ongoing LHC searches: mono-jet, disappearing tracks and displaced vertices + E miss T . In section 7 we derive the future prospects of an FCC-hh 100 TeV collider with a forward detector, and finally we conclude in section 8. Various technical details are confined to the appendices: appendix A discusses 1 → 2 vs 2 → 2 processes regarding the DM relic density. Appendix B presents details of the recasting procedure for the ATLAS displaced vertices + E miss T search [45]. Appendix C explicitly presents the NLO+NLL Higgsino pair production cross sections obtained with Resummino-2.0.1 [46,47] and used for our sensitivity estimates throughout the paper.

Dark Matter through the freeze-in mechanism: review
We review here the main aspects of freeze-in production of DM through the decay of a parent particle A in thermal equilibrium with the plasma [8] A → B SM χ , (2.1) where B SM is a SM state and χ is the DM candidate. Initially (after reheating at the end of inflation) the DM abundance is assumed to be negligibly small and subsequently increases continuously as the parent particle A within the thermal bath decays into DM during the radiation-dominated era. 4 This DM production process is effective as long as the parent particle is relativistic. As the temperature of the radiation bath drops below the mass of A, the abundance of the parent particle becomes exponentially suppressed and the DM production process ceases to be effective. DM is then said to freeze-in, with most of the χ particles produced at temperatures T ∼ m A /3. The DM number density n χ evolves according to the Boltzmann equation [8] dn χ dt

JHEP03(2020)022
where Γ A is the decay rate producing DM in eq. (2.1) and K 1,2 (x) are the first and second modified Bessel functions of the 2 nd kind. The parent particle equilibrium number density n eq A can be well-approximated by where g A counts the spin d.o.f. of A and ξ = 2 if the decaying particle is not self-conjugate (otherwise ξ = 1). At high temperatures (T m A ), n eq A ∼ T 3 characteristic of a relativistic species, while it features the Maxwell-Boltzmann exponential suppression for T m A . The Boltzmann equation (2.2) is general, with the details of the cosmological history entering through the Hubble parameter H and the time vs temperature relation. Conservation of entropy s = 2π 2 45 g s * T 3 (with g s * being the number of effective degrees of freedom contributing to the entropy density) in general makes it convenient to track the evolution of the DM population by defining the comoving number density (yield), Y χ ≡ n χ /s, which is dimensionless, and on the other hand allows us to relate time and temperature as dt = −dT /(HT ), where (see e.g. [49])H For a radiation-dominated universe the Hubble parameter is H 1.66 with Y eq A = n eq A /s. Assuming that at very high temperatures T R m A the abundance of χ is vanishing, this equation can be integrated to give with x min = m A /T R 1 and x max → ∞. Far from phase transitions where the number of d.o.f. in the plasma could potentially change abruptly, the DM yield is to a good approximation Requiring then the relic abundance of χ (with ρ c /s 0 = 3.6 × 10 −9 GeV the critical energy density over the entropy density today) to match the observed DM relic abundance Ω DM h = 0.12, we can obtain an estimate of JHEP03(2020)022 the corresponding decay length cτ A of the parent particle: where we have assumed for simplicity g s * (x) = g * (x) at T = m A /3 (approximately the freeze-in temperature, see [49]). We note the macroscopic lifetime for the parent particle A, which according to eq. (2.9) increases linearly with the DM mass m χ as a result of requiring Ω χ h = 0.12.

Freeze-in (from decay) via the Standard Model Higgs
In the present work we concentrate on freeze-in DM scenarios where the parent particle A in eq. (2.1) is neutral, which are significantly more challenging to probe experimentally than those with a charged parent particle, 5 particularly if the decay length of A is large, cτ A 1 m. In the following we consider B SM in eq. (2.1) to be the SM Higgs [29,32,43].

The model
A simple model which features the SM Higgs as B SM consists on adding to the SM a Dirac fermion χ, singlet under the SM gauge group, and an SU(2) L Dirac fermion doublet Ψ with hypercharge 1/2 such that the Lagrangian reads: This model is very similar to the singlet-doublet DM model considered in [32,44], and could be regarded as a simplified version of the Higgsino-Axino system studied in [10,29,42,43], or an analog of a would-be feebly interacting Higgsino-Bino or Higgsino-Singlino system, bearing in mind that here χ and Ψ are Dirac fermions. The neutral particles χ and ψ 0 mix after EW symmetry breaking due to the presence of the y χ coupling, giving rise to mass eigenstates χ 1 (mostly singlet-like) and χ 2 (mostly doublet-like), their respective masses being m 1 and m 2 . We consider m 2 > m 1 and the coupling y χ 1, as needed for the freeze-in mechanism to yield the correct DM relic density. In this limit we have m 1 m s and m 2 m D in eq. (3.2) (m 2 > m 1 is then achieved by setting m D > m s ). The singlet-doublet mixing is simply given by directly from the Yukawa term in eq. (3.2), and Z − χ 2 − χ 1 , W ± − ψ ∓ − χ 1 from the singlet-doublet mixing induced by the Yukawa term after EW symmetry breaking. When kinematically possible, the decay widths for χ 2 → hχ 1 , χ 2 → Zχ 1 and ψ ± → W ± χ 1 are given by with λ(x, y, z) = x 4 + y 4 + z 4 − 2x 2 y 2 − 2x 2 z 2 − 2y 2 z 2 . In figure 1 (left) we show the branching fractions of χ 2 as a function of m 2 for a benchmark 6 m 1 = 10 MeV, highlighting that for m 2 400 GeV the decay χ 2 → Zχ 1 starts to dominate over χ 2 → hχ 1 . For m 2 −m 1 → m h the branching fraction into hχ 1 becomes negligible, while for m 2 −m 1 m h both branching fractions approach 50%.
In addition, loops of EW gauge bosons induce a radiative mass splitting between the charged state ψ ± and the (mostly doublet-like) neutral state χ 2 [50,51]. The value of this mass splitting δm is given by with g the SU(2) L gauge coupling, θ W the Weinberg angle and f (m Z /m 2 ) given by (3.6) 6 Compared to the results in [32], ours correspond to the choice tan θ = 1 in [32].

JHEP03(2020)022
For m 2 ∈ [100, 1000] GeV the range of mass splitting values is δm ∈ [260, 340] MeV. For such splittings the decay mode ψ ± → π ± χ 2 completely dominates over the decay involving charged leptons ψ ± → ± νχ 2 [50], with the width Γ(ψ ± → π ± χ 2 ) given by with f π 130 MeV and θ c the Cabibbo angle. For the very small values of the singletdoublet mixing relevant for DM freeze-in, the decay width from eq. (3.7) is generally much larger than Γ(ψ ± → W ± χ 1 ) and yields the dominant decay mode for ψ ± . This can be seen explicitly in figure 1 (right), where we show the branching ratio BR(ψ ± → W ± χ 1 ) as a function of m 2 m ψ ± and the χ 2 decay length cτ χ 2 (which can be directly related to Γ(ψ ± → W ± χ 1 ) by means of eq. (3.4)) for a benchmark DM mass m 1 = 12 KeV. We also show in figure 1 (right) the DM relic density curve for m 1 = 12 KeV as computed with micrOMEGAs5.0 [49] (see section 3.2 for details); since roughly all parameter space below this curve is excluded by the Lyman-α bound from Cosmology (see section 4), figure 1 (right) highlights that, in order to obtain the observed DM relic abundance, BR(ψ ± → W ± χ 1 ) 0.1 for m 2 < 1 TeV is needed in the present freeze-in scenario.

DM relic density: the role of EW symmetry breaking & thermal masses
As discussed in section 2, the DM relic abundance is obtained via slow χ 1 production from the decays of the parent particles χ 2 , ψ ± present in the thermal bath during the radiation-dominated era. For the computation of the DM relic density, we have used micrOMEGAs5.0 [49] with the additional implementation of several important features, which we detail in the following. Prior to EW symmetry breaking in the early Universe (the EW phase transition), the interactions of DM with the EW gauge bosons are absent, 7 and the states ψ 0 , ψ + in (3.1) decay via the Yukawa interaction y χ into the SM Higgs doublet field components and the DM candidate, Ψ → Hχ. The correct evaluation of DM production in the EW symmetric phase requires the inclusion of the temperature-dependent (thermal) mass Π(T ) of the Higgs doublet, given by [52,53] (see also [54]) with g, g the gauge couplings for SU(2) L and U(1) Y , y t the top quark Yukawa coupling and λ the Higgs quartic coupling. On the other hand, since the parent particle is in equilibrium with the SM, it will also acquire a thermal mass (see e.g. [55,56]) 8 Π Ψ (T ): 9) 7 We are indebted to Thomas Konstandin for discussions on this point. 8 Here we use the thermal mass of a Higgsino (or analogously, that of a lepton SU(2)L doublet with hypercharge Y = 1/2). We note that the DM coupling of Ψ may be safely neglected for this purpose.

JHEP03(2020)022
such that before the EW phase transition we consider the mass of Ψ to be The temperature dependence from both Π H (T ) and Π Ψ (T ) will then be inherited by the parent decay width in the EW symmetric phase 9 with T EW the EW phase transition temperature, corresponding to T EW ≈ 160 GeV in the SM [57,58]. Since the thermal mass of the particle Ψ is smaller than the thermal mass of the Higgs, we have that for T m D ( m 2 ) the decay Ψ → Hχ will not be kinematically open. When this is the case, it is instead the Higgs decay H → Ψχ which becomes kinematically open and sources DM production, its rate being Γ(H → Ψχ) Γ(Ψ → Hχ) × (m Ψ /Π H ) 3 (with m Ψ and Π H given respectively by (3.10) and (3.8)). The temperature at which this change in DM production process approximately occurs is T 2.96 m 2 . For lower temperatures (yet generally above T EW ), the decay Ψ → Hχ sources DM production and yields the (by far) leading contribution to the DM relic density prior to the EW phase transition.
After the EW phase transition the decays χ 2 → Zχ 1 and ψ ± → W ± χ 1 become possible and source DM production together with the decay χ 2 → hχ 1 . The decay width Γ A of the parent particle(s) responsible for DM freeze-in is in this case given by where for simplicity we do not consider the temperature dependence of the Higgs vev after EW symmetry breaking and also neglect thermal corrections for temperatures below the EW phase transition. The temperature dependence of the freeze-in decay processes from eqs. (3.11) and (3.12) has been implemented in a modified version of micrOMEGAs5.0, and in the remainder of this work we consider the range between T EW = 160 GeV and T EW = 50 GeV (the latter corresponding to a generic benchmark for a strongly super-cooled first order EW phase transition) when discussing/computing the freeze-in DM relic density.
Both the EW phase transition and the inclusion of thermal masses in the calculation can then have a notable influence on the DM freeze-in process. 10 In particular, the DM production prior to the EW phase transition can be strongly affected by thermal corrections. The production of DM after the EW phase transition is however not expected to be significantly affected, and as discussed above we are for simplicity neglecting here the thermal effects below T EW . It becomes clear that thermal effects will be more important JHEP03(2020)022 when the contribution to the total DM relic density from DM production for T > T EW is sizable. In figure 2 we show the DM relic density with the inclusion of thermal corrections Ωh 2 TC and the corresponding DM relic density when these are disregarded 11 Ωh 2 NoTC , as well as the partial contributions before and after the EW phase transition, as a function of m 2 for m 1 = 1 GeV and T EW = 50 GeV (left panel), T EW = 100 GeV (middle panel) and T EW = 160 GeV (right panel).
As seen in figure 2, for m 2 m h + m 1 , T EW the inclusion of thermal masses yields a very mild suppression Ωh 2 TC /Ωh 2 NoTC ∼ 0.9. In this case, DM freezes-in while the EW symmetry is restored, and the broken phase decays χ 2 → Zχ 1 and ψ ± → W ± χ 1 do not contribute to DM production in the early Universe. In this regime we find the DM relic abundance to be rather insensitive to the value of T EW . On the other hand, figure 2 shows that the DM production in the EW broken phase starts becoming appreciable for m 2 6 T EW , and for yet smaller m 2 the DM production in the broken phase quickly becomes dominant when T EW m h . In contrast, for T EW m h the DM production prior to the EW phase transition is still very important for smaller m 2 masses, and it is in this region where the inclusion of the thermal corrections yields the largest effect, providing a significant enhancement to the DM relic density compared to the value neglecting thermal 11 Here thermal corrections are disregarded by setting ΓΨ→Hχ 4 × Γ(χ2 → hχ1)| m h =125 GeV . JHEP03(2020)022  corrections. The difference in behaviour of the two regimes (T EW much above or below the Higgs mass) is due to the fact that the latter case allows for the decay process prior to the EW phase transition to be kinematically open for a longer period. All this may be seen explicitly in figure 3 where Ωh 2 TC /Ωh 2 NoTC is shown as a function of m 2 for several values of m 1 and T EW . Note that below m 2 m h the enhancement to the DM relic abundance coming from thermal corrections is significantly larger for m DM = 10 GeV than for m DM = 1 GeV or less. This is due to phase space arguments, the parent decay being open or closed in the EW symmetric phase depending on the inclusion of thermal corrections or not, mostly for the lowest T EW considered.
Finally, we stress that in the above discussion, and throughout the general discussion of freeze-in from decay in [8], outlined in section 2, it has been implicit that (1 → 2) decays are the dominant processes producing DM. We nevertheless note that there are contributions to DM production at the same order in the feeble coupling y χ from 2 → 2 scattering processes such as q χ 2 → q χ 1 (with q a SM quark) mediated by a Higgs boson in the t-channel. In a similar manner, for T < T EW analogous q χ 2 → q χ 1 and q χ ± → q χ 1 scattering processes with a Z and W boson in the t-channel become possible. These processes are discussed in detail 12 in appendix A, where we find that decays A → B SM χ 12 Appendix A does however not analyze in detail the impact of other 2 → 2 processes in the present freeze-in scenario involving SM gauge bosons VSM, such as HΨ → χVSM and VSMΨ → χH. While such processes may be quantitatively more important than scattering processes q Ψ → q χ [60], they are still expected to be highly subdominant to 1 → 2 decays. A precise assessment of the impact of these processes is left for future work.

JHEP03(2020)022
completely dominate over scatterings X SM A → Y SM χ (with X SM and Y SM SM particles, and B SM in the t-channel of the scattering process) except when the former are strongly phase-space suppressed, i.e. m B SM + m χ → m A . In addition, other contributions like soft gauge interaction processes (in 1 → 2 and 2 → 2 DM production processes) before the EW phase transition (see [60] for a discussion) are not expected to be significant in the present case (contrary to e.g. the Type-I see-saw scenario [60]), due to the specific mass spectrum of our freeze-in scenario, with m 2 > 100 GeV ( m ) and m 2 m χ . 13 Finally, scatterings do become dominant in the region m B SM + m χ > m A > m χ where 1 → 2 decays are kinematically forbidden and instead A decays 14 via a 3-body process. In our model, the latter situation occurs for m 2 − m 1 < m W , when all the decays of χ 2 and ψ ± that produce DM are 3-body (e.g. χ 2 → bb χ 1 through an off-shell Higgs boson). Since we do not consider these parameter space regions in the present work due to experimental constraints (e.g. the region m ψ ± < m W is strongly constrained by LEP [61,62]), 2 → 2 scattering processes can be safely disregarded in our scenario.

Super-WIMP contribution to the DM relic abundance
For decay widths of the parent particle significantly smaller than its inverse freeze-out timescale which in the present scenario roughly corresponds to cτ χ 2 (100 GeV/m 2 ) 2 × 3 meters, the thermal freeze-out of χ 2 followed by its decay into χ 1 will also contribute to the DM relic density. The contribution to the DM relic abundance from this mechanism is given by Ω 2 h 2 × m 1 /m 2 , where Ω 2 h 2 is the freeze-out abundance of χ 2 . Assuming that χ 2 undergoes freeze-out before decaying as implied by eq. (3.13), the freeze-out abundance of χ 2 corresponds to that of thermal Higgsino DM, which for m 2 m W is simply given by [63] (3.14) When Ω 2 h 2 ×m 1 /m 2 yields the dominant contribution to the DM relic density, it is usually referred to as the super-WIMP [64] mechanism (see also [15] for a recent work on the interplay between the freeze-in and super-WIMP scenarios). Combining the super-WIMP and freeze-in contributions to the DM relic density, we get the total DM relic abundance The above discussion allows to quantify the relative importance of the super-WIMP and freeze-in contributions to the DM relic abundance, with their ratio (SW/FI) given by (m 1 Ω 2 h 2 )/(m 2 Ω 1 h 2 ). In figure 4 we show the SW/FI ratio in the (m 2 , cτ χ 2 ) plane, fixing here for simplicity T EW = 160 GeV and two benchmark values for the DM mass, m 1 = 13 A detailed analysis of these processes in freeze-in scenarios of the type discussed in this paper is however necessary to quantitatively assess this statement, and still lacking in the literature. 14 We assume that BSM is itself unstable, otherwise for mB SM + mχ > mA the state A becomes stable.
JHEP03(2020)022  1 GeV and m 1 = 300 GeV (we nevertheless note that the SW/FI ratio is approximately independent of m 1 , except for m 2 − m 1 → m h , as can be seen in figure 4). We also show in figure 4 the corresponding DM relic abundance condition Ω DM h 2 = 0.12 using micrOMEGAs5.0 and eq. (3.15) for various choices of the DM mass m 1 . As is apparent from eq. (3.14), in the present scenario the super-WIMP mechanism could only account for a significant fraction of the observed DM relic density for m 2 > 1.1 TeV. In addition large DM masses m 1 300 GeV are required, as shown in figure 4 and expected in order to partially overcome the otherwise large suppression factor m 1 /m 2 in the super-WIMP contribution from (3.15). Such large DM masses imply large values of cτ χ 2 for DM not to be overproduced by the freeze-in mechanism, which in turn may become constrained from Big Bang Nucleosynthesis (BBN). We discuss these constraints in the next section.

Constraints on Dark Matter freeze-in from cosmology
As previously mentioned, there are several important cosmological constraints on this class of freeze-in DM scenarios (for a discussion of such constraints in similar scenarios, see also [32,33]). The first constraint we need to consider comes from Big Bang Nucleosynthesis (BBN), which accurately explains the measured primordial abundances of light elements in the Universe 15 (see e.g. [66,67]). If the state χ 2 lives long enough to decay during BBN, JHEP03(2020)022 its visible decay products may induce several processes that alter the predictions of BBN 16 In order to obtain the corresponding bounds on the parameter space of m 1 , m 2 and cτ χ 2 , we use the recent re-analysis of BBN constraints on long-lived decaying particles from [69]. Concerning the hadronic final statesqq (with q either a first/second generation quark or a b quark), ref. [69] derives bounds on the amount of energy injected in the radiation bath in the form of hadrons from the decay of a long-lived particle X, as a function of the particle lifetime τ X and under the assumption that X decays dominantly intoqq, 17 meaning that the hadronic energy injection per particle decay is Eq q ∼ m X . In the present case the hadronic energy injection from the decay of χ 2 corresponds only to the energy fraction carried by the visible decay products of χ 2 , either h or Z (both decays χ 2 → hχ 1 and χ 2 → Zχ 1 are present during the BBN epoch, which occurs for temperatures much below the EW phase transition temperature T EW ). This is given by From [69] we then obtain the limits on the total hadronic energy injection onto the thermal bath from χ 2 decays, given by with Ω 2 h 2 from (3.14). Assuming m 1 m 2 , m h , the corresponding limits on the (m 2 , cτ χ 2 ) parameter space are shown in figure 5 (left). We note that leptonic decays of χ 2 also affect BBN, yet these constraints are found to be significantly weaker than the hadronic ones (see [69]) and are therefore not relevant in the present case.
A second important set of limits originate from constraints on the washout of smallscale structure by DM with a non-negligible velocity dispersion (partially relativistic, or warm). The leading such constraint comes from Lyman-α forest observations. For thermal warm dark matter (WDM) these observations place a lower limit on the DM mass in the range m DM 4.09 − 5.3 keV [70][71][72]. In two recent works [73,74] (see also [75,76]) the Lyman-α limit for the case of freeze-in DM produced via two-body decays of a parent particle in thermal equilibrium with the plasma (which is precisely our scenario) has been estimated by comparing the suppression in the linear matter power spectrum from a thermal WDM scenario (with m DM given by the WDM Lyman-α limit, taken to be 4.65 keV [71,77]) to that of the freeze-in scenario. The Lyman-α bound in this case reads 3) 16 For even longer lifetimes τ 10 6 seconds, constraints from the Cosmic Microwave Background (CMB) spectral distorsion become very strong (see [68,69]). Such long lifetimes are however not relevant for our study. 17 The corresponding limits for thebb andūu final states found in [69] are very similar, and thus approx- where the sum runs over all decay channels of the parent particle(s) that contribute to DM production, Γ ij are the corresponding decay widths Γ(A i → B j χ 1 ) (in our case given by eq. (3.4)), g i are the number of degrees of freedom of the parent particle A i , the parameter ∆ ij yields the mass splitting between the parent particle and the visible decay product in each channel (e.g. ∆ ij = 1 − m 2 Z /m 2 2 for the decay χ 2 → Zχ 1 ) and η 1.9 [74]. Assuming that freeze-in DM saturates the observed DM relic density, figure 5 (right) shows the Lyman-α bound on m 1 in the present scenario. This bound can in turn be translated into a lower limit for cτ χ 2 , shown in figure 5 (left).
The combination of BBN and Lyman-α bounds yields an allowed range of parent particle decay lengths cτ χ 2 for viable DM freeze-in, cτ χ 2 ∼ 10 −1 −10 7 meters, which provide a clear target for LLP searches at colliders. We investigate the current and projected sensitivity of such searches to the Higgs-assisted freeze-in from decay DM scenario analyzed in this work in sections 5-7. Such long lifetimes make it challenging to search for χ 2 with the LHC's ATLAS and CMS detectors (we will nevertheless discuss the current ATLAS and CMS bounds, as well as future prospects in section 6), except maybe for rather small values of m 1 close to the Lyman-α bound [32]. Here we analyze the prospects for probing the parameter space of the present scenario with the proposed MATHUSLA surface detector for long-lived JHEP03(2020)022 particles [28,[35][36][37]. The large detector volume of MATHUSLA together with its capability to effectively operate as a background-free environment for LLP searches (see [28,37] for details) at the LHC, yields an ideal setup to probe freeze-in DM scenarios. At the LHC, the states χ 2 and ψ ± can be produced via the Drell-Yan processes pp → χ 2 χ 2 , pp → χ 2 ψ ± , pp → ψ + ψ − . For our study, we implement the Lagrangian from eq. (3.2) in FeynRules [78] and simulate the various Drell-Yan production processes for χ 2 and ψ ± at √ s = 13 TeV LHC in Madgraph aMC@NLO [79]. We then normalize the respective cross sections to the corresponding √ s = 13 TeV LHC NLO+NLL charged/neutral Higgsino production cross sections σ LHC 13 computed with Resummino-2.0.1 [46,47] 18 (with the PDF set MSTW2008nlo90cl from LHAPDF [81]). The Higgsino pair production cross sections for the various Drell-Yan processes are shown explicitly in appendix C. In the following, we also consider that the decays ψ ± → χ 2 + X (X = π ± , ± ν being very soft), which occur within a few cm of the interaction point, effectively convert all the charged states ψ ± into neutral χ 2 ones with approximately identical kinematics.
The MATHUSLA detector would be located at a distance ∼ O(100 m) from the ATLAS or CMS interaction point. The probability for a χ 2 particle to decay inside the MATHUSLA detector volume is given by with P decay (β cτ χ 2 , L a , L b ) given by Here L a and L b are the distances from the interaction point for which χ 2 respectively enters and leaves the MATHUSLA detector volume, and β is the χ 2 boost factor | p χ 2 |/m 2 . In this work we consider two possible MATHUSLA design configurations: MATHUSLA100 and MATHUSLA200. In the former case, the MATHUSLA detector volume for LLP decays as measured from the interaction point is [37] x ∈ [−50, 50] m, y ∈ [100, 120] m, z ∈ [100,200] m, while in the latter case (corresponding to the original MATHUSLA proposal [35,37]) the detector volume as measured from the interaction point is x ∈ [−100, 100] m, y ∈ [100, 120] m, z ∈ [100, 300] m. The MATHUSLA geometric acceptance geometric in eq. (5.2) then measures the fraction of generated events contained within the MATHUSLA solid angle as seen from the interaction point, which we extract from our event simulation. 19 Specifically, for the computation of geometric we consider events for which the χ 2 trajectory intersects the MATHUSLA tracking layers (see [35,36] for details on the tracking layout within the MATHUSLA detector) and assume for simplicity that the visible χ 2 decay products would then also hit the tracking layers if the decay happens inside MATHUSLA (for m 2 m h + m 1 this occurs automatically, since the decay products of χ 2 will be fairly collinear to the χ 2 trajectory, while for m 2 → m h + m 1 addressing the possible modification of our acceptances requires a more detailed event and detector simulation beyond the scope of this work). We also assume perfect MATHUSLA detection efficiency for χ 2 decays inside the detector volume (see nevertheless [82] for a discussion on this point).
The number of expected signal events at MATHUSLA is then given by where L is the integrated luminosity and the integral over P MATH decay denotes a generic integration over phase-space and decay length for the signal event sample. Considering the MATHUSLA detector as an essentially background-free environment (see [37] for a detailed discussion on this point), the 95% C.L. exclusion sensitivity corresponds to N events = 3.
In figure 6 we show the MATHUSLA100 and MATHUSLA200 95% C.L. sensitivities in the (m 2 , cτ χ 2 ) plane for an integrated luminosity L = 300 fb −1 (solid) and L = 3000 fb −1 (dashed). We also show the lines yielding the observed DM relic abundance for m 1 = 1 GeV (black), m 1 = 10 MeV (dark green) and m 1 = 100 KeV (light green) obtained with micrOMEGAs5.0, assuming the temperature of the EW phase transition to be in the range T EW ∈ [50 GeV, 160 GeV]. The results of figure 6 show that MATHUSLA could cover a wide region of the viable freeze-in parameter space between the Lyman-α and BBN constraints, and probe DM masses up to m 1 ∼ 1−10 GeV. In the next section, we compare these sensitivity projections with those of searches by ATLAS and CMS, to investigate the complementarity between MATHUSLA and the existing LHC detectors.

T (mono-X) signatures
Due to the long lifetime of χ 2 , standard "mono-X" searches for DM at the LHC may be sensitive to the freeze-in parameter space when the decay of χ 2 happens outside the ATLAS/CMS detector. Here we analyze the sensitivity of such "traditional" DM searches, focusing on the mono-jet signature which generally yields stronger constraints as compared e.g. to mono-γ. Considering the state χ 2 to be invisible (decaying outside the relevant detector volume), jet + E miss T signatures can be obtained in the current scenario via Drell-Yan production pp →χ 2 χ 2 , pp → χ 2 ψ ± , pp → ψ + ψ − in association with an initial-stateradiation (ISR) jet, see figure 7 (left). In the latter two processes ψ ± decays into χ 2 and a very soft π ± . The particle ψ ± will in principle leave a short track in the detector's inner tracker before decaying, but since the decay length of ψ ± is O(1 cm) (see eq. (3.7) and section 6.2), such a short track typically does not affect the mono-jet event selection. 20 We first focus on the Drell-Yan pp →χ 2 χ 2 process with an ISR jet, which is mediated by an off-shell Z-boson. This scenario can be reinterpreted as a simplified model for DM (with χ 2 as the invisible particle) with a vector mediator [83,84] with m med = m Z , and with the mediator couplings to SM quarks (g q ) and to χ 2 (g χ ) being fixed by the respective gauge quantum numbers of the quarks and χ 2 . This allows us to use the existing constraints on simplified DM models with vector mediators from the CMS mono-jet analysis at √ s = 13 TeV with an integrated luminosity L = 35.9 fb −1 [6]. Moreover, since the p T distribution of the ISR jet is similar for all three processes pp →χ 2 χ 2 , pp → χ 2 ψ ± , pp → ψ + ψ − , as shown in figure 7 (right), we can to a good approximation rescale the cross 20 The mono-jet analysis features only indirect track vetoes, in the form of well-identified electrons, muons and hadronic τ leptons, all of which require further activity in the calorimeters or muon stations. That way, the processes pp → χ2ψ ± , pp → ψ + ψ − with an ISR jet may also contribute to the mono-jet signal. We thank Steven Lowette for clarifications on this point.  section for pp →χ 2 χ 2 + j by the total cross section including all three processes (noting also that there is no interference among the different processes). We then use the publicly provided 95% C.L. exclusion limit on the mono-jet signal strength µ CMS (mono-jet cross section divided by the theoretical cross section for a DM simplified model with a vector mediator and g q = 0.25, g χ = 1) as a function of the DM mass from the CMS analysis [6] for m med = m Z , weighted by the decay probability of χ 2 inside the CMS detector. This probability is approximately given by exp[−L CMS /cτ χ 2 ], withL CMS a mean CMS relevant detector size taken here to be the average hadronic calorimeter radius,L CMS ∼ 2.5 meters. 21 The corresponding CMS 95% C.L. exclusion limit on the present scenario is then given by equating where µ CMS (m 2 ) is the CMS exclusion limit given in [6] for m DM = m 2 and m med = m Z , g Z L and g Z R are the couplings of the Z-boson to the left and right-handed SM quarks 22 and r σ (m 2 ) is the ratio between the cross section for pp →χ 2 χ 2 + j and the total cross section including all other Drell-Yan processes (which also involve ψ ± production), which we 21 More accurately, the detector geometry should be convoluted with the rapidity distribution of our signal events and the probability of vetoing a signal event in the mono-jet analysis if the decay happens in different parts of the CMS detector. This is however beyond the scope of this work. The exponential fall-off of the mono-jet sensitivity for small values of cτχ 2 shown in figure 8 should then only be taken as approximate. 22 Here we take the up-type quark couplings, which yield the dominant contribution, to compute the exclusion limit. In practice, both the up and down-type quark contributions should be taken into account, bearing in mind that they are slightly different and this difference can be taken into account by re-weighting the different PDF contributions to the Drell-Yan process as a function of the momentum transfer. We however do not perform such re-weighting here.
Using a simple √ L luminosity rescaling, we also provide a naive projection of the limit to L = 300 fb −1 with √ s = 13 TeV, shown in figure 8, yielding m 2 210 GeV in the limit cτ χ 2 L CMS . We however stress that since this projection has been derived assuming the SM background uncertainties in the mono-jet search are dominantly statistical, and in reality the mono-jet search is systematics dominated and so the true mono-jet sensitivity for L = 300 fb −1 should be significantly weaker than the above value of m 2 . On the other hand, future improvements in the systematic errors would make our (overly aggressive) limit more realistic.

Disappearing track signatures
Due to the very small mass splitting δm ∼ 300 MeV between the states ψ ± and χ 2 (recall eq. (3.5)), the state ψ ± is relatively long-lived, with its decay length in the range cτ ψ ± ∈ [0.7, 2.1] cm for m ψ ± between 90 GeV and 1 TeV. The dominant decay ψ → χ 2 π ± can then lead to a disappearing track signature, due to the softness of the produced pion. Both CMS [85] and ATLAS [86] have performed searches for disappearing tracks with LHC 13 TeV data and an integrated luminosity of 36.1 fb −1 (38.4 fb −1 ) in the case of ATLAS (CMS). We note that ATLAS can reconstruct tracks as short as ∼ 12 cm, while for the CMS tracker the minimum reconstructed track length is ∼ 25 − 30 cm (see e.g. the discussion in [33,87]). As a result, the ATLAS search can constrain smaller lifetimes than the CMS study. In fact, the CMS search [85] only constrains chargino lifetimes τ > 0.07 nanoseconds (ns) corresponding to cτ > 2.1 cm, already at the edge of the maximum value of cτ ψ ± possible in the present scenario. Using the exclusion from the τ ∈ [0.07, 0.1] ns, m 2 ∈ [50, 150] GeV bin in the CMS analysis [85] (yielding σ × B < 17.60 pb) and using the chargino production cross sections provided in [80], we derive the tentative current bound m ψ ± > 84 GeV, which is in fact weaker than existing LEP constraints on the state ψ ± [61,62]. In contrast, the ATLAS collaboration recently reinterpreted their analysis [86] in terms of pure Higgsino 23 production [88], obtaining a present 95% C.L. exclusion of m ψ ± 145 GeV. In addition, the ATLAS collaboration has performed a preliminary study of their future sensitivity to a disappearing track signature in the pure Higgsino scenario [89], which would yield a 95% C.L. sensitivity of m ψ ± ∼ 260 GeV at the HL-LHC with √ s = 14 TeV and L = 3000 fb −1 . We show the corresponding limits in figure 9.

Displaced Vertex (DV) signatures
Neutral long-lived particles may be searched for via displaced vertex signatures in ATLAS, CMS and LHCb, and a wide variety of analyses already exist for LHC √ s = 7 , 8 and JHEP03(2020)022 • CMS/ATLAS searches looking for displaced lepton pairs [90], which would be sensitive to decays χ 2 → Z χ 1 (Z → ), and to a lesser extent (due to the much smaller leptonic branching fraction) to χ 2 → h χ 1 (h → ν ν ).
• ATLAS searches looking for activity in the muon spectrometer [92,96], which could in principle be sensitive to very large decay lengths cτ χ 2 1 m.
• ATLAS searches for displaced vertices + E miss T [45]. These would take advantage of the large hadronic branching fraction of both Z and h, as well as of the missing momentum in the decays χ 2 → Z/h χ 1 .
In this work, we concentrate on the latter DV + E miss T search by ATLAS [45] at √ s = 13 TeV with 32.8 fb −1 , leaving an exploration of the other two searches highlighted above for future work. The ATLAS DV + E miss T search [45] targets events with large missing transverse momentum and one or more displaced vertices with large track multiplicity (thus likely to correspond to hadronic decays). The analysis provides very detailed documentation 24 for validation and recasting in its auxiliary material [97], and we give details on the analysis selection in appendix B. In order to validate our approach, we have applied the analysis to the model used by ATLAS for interpretation of their results, corresponding to pair-production of long-lived gluinos, with the gluino eventually decaying to a pair of SM quarks and a neutralino. The results of our validation are shown in appendix B, which allow us to perform a recast of the ATLAS search in terms of the Higgs-assisted freeze-in DM scenario explored in this work (an analogous recasting analysis for this scenario has already been performed in [32], with similar results).
As in section 5, we use the Lagrangian implementation of the freeze-in DM model in FeynRules [78]. We simulate Drell-Yan production pp → χ 2χ2 with the decays of the χ 2χ2 pair to Zχ 1 Zχ 1 , Zχ 1 hχ 1 or hχ 1 hχ 1 in Madgraph aMC@NLO [79]. Since we find the kinematics of χ 2χ2 , χ 2 ψ ± , and ψ + ψ − production to be very similar (both without and with the production of an extra hard jet, cf. figure 7), we normalize our signal cross section to the sum of respective LHC √ s = 13 charged/neutral Higgsino production cross sections at NLO+NLL obtained from Resummino-2.0.1 [46,47]. After parton level event generation (we use the time of flight option within Madgraph aMC@NLO to introduce a non-zero displacement in the decay of χ 2 ), our events are passed to Pythia8 [98] which decays the Z and h bosons, and performs showering and hadronisation. Finally we cluster jets using the FastJet [99] implementation within Delphes 3 [100]. Further details of our analysis are given in appendix B. The derived LHC 13 TeV 95% C.L. limits from the ATLAS DV + E miss T search with 32.8 fb −1 on the (m 2 , cτ χ 2 ) parameter space of our present scenario are shown in figure 9, corresponding to a number of signal events N events = 3 (since the expected number of 24 We note that CMS analyses searching for a events with an electron and a muon with large impact parameters [95] are very well-documented and sensitive to the present scenario via decays χ2 → hχ1 (h → eνe µνµ). The very small Higgs leptonic branching fraction (∼ 0.4 %) results in a much weaker sensitivity w.r.t. DV + E miss T searches (see e.g. [94]). background events in the signal region is 0.02 +0.02 −0.01 , the current search may be considered as background-free). We also show in figure 9 the 95% C.L. sensitivity projection for the DV + E miss T search with 300 fb −1 and 3000 fb −1 . In the former case, the search can still be regarded as background-free and we use N events = 3, whereas in the latter we expect ∼ 2 background events, and the corresponding required number of signal events is found via the CLs method [101] to be N events = 6. Finally, we also include in figure 9 the present and projected 95 % C.L. bounds from disappearing track searches (see section 6.2) and the projected 95 % C.L. sensitivity from MATHUSLA100 and MATHUSLA200 (see section 5) to highlight the complementarity among the various searches. While for decay legths cτ χ 2 100 m the LHC DV + E miss T search yields the most sensitive probe of these scenario, for larger decay lengths MATHUSLA provides a major increase in sensitivity w.r.t. the main LHC detectors. face detector proposal to search for very long-lived particles. In this section we analyze the sensitivity to the Higgs-mediated freeze-in DM scenario explored in this work of a √ s = 100 TeV hadron collider (from now on referred to as FCC-hh) and a forward LLP detector with volume given by z ∈ [20,40] m, ρ = x 2 + y 2 ∈ [5, 30] m as measured from the FCC-hh interaction point.
We simulate Drell-Yan production pp → χ 2χ2 , pp →χ 2 ψ − , pp → χ 2 ψ + , pp → ψ + ψ − at √ s = 100 TeV in Madgraph aMC@NLO, normalizing the respective cross sections to the corresponding NLO+NLL pure Higgsino production cross sections computed with Resummino-2.0.1, and shown in appendix C. The decays of the charged states into χ 2 andχ 2 are again considered to yield neutral states with approximately identical kinematics to that of the decaying charged states. We follow the same procedure used in section 5 for MATHUSLA to estimate the event yield at the forward detector, computing from simulation the number of signal events N events from the probability for a particle χ 2 to decay within the forward detector volume (we assume perfect detector performance as well as a background-free environment), given by (5.3). The forward detector geometric acceptance JHEP03(2020)022 for our signal is found from simulation to be geometric ∼ 0.5, which is roughly ten times larger than the one of MATHUSLA due to the rather forward nature of our Drell-Yan events for √ s = 100 TeV proton-proton collisions.
In figure 10 we show our FCC-hh forward detector 95% C.L. sensitivity projections, given by N events = 3, for two choices of FCC-hh integrated luminosity, L = 3 ab −1 and L = 30 ab −1 . The mass reach shows a dramatic increase w.r.t. MATHUSLA and HL-LHC, 25 being able to probe LLPs with cτ χ 2 ∼ 50 m up to m 2 ∼ 10 TeV. figure 10 also shows that the forward detector would be able to probe lifetimes up to the BBN bound for m 2 600 GeV, and approach the boundary of the super-WIMP parameter space region (solid black line in figure 10).

Conclusions
Freeze-in constitutes a well-motivated and appealing mechanism for DM production in the early Universe, yet challenging to probe experimentally given the feeble interactions between DM and the visible sector. In this sense, collider searches for LLPs have been recently regarded as possible probes of DM freeze-in production via the decay of a parent particle which is accessible at colliders. In this work we have studied the prospects of probing DM freeze-in from the decay of neutral parent particles which belong to the thermal bath, through LLP signatures. Using as a case-study a model featuring Higgs-associated LLP decays, we have obtained the experimental sensitivity of the proposed MATHUSLA surface detector and ATLAS/CMS at HL-LHC to the cosmologically relevant freeze-in DM parameter space. Our results show the high degree of complementarity between MATH-USLA and the main LHC detectors, together being able to probe DM masses from the Lyman-α bound (of a few keV) up to masses of few GeV, and neutral parent particle masses up to ∼ 2 TeV. Besides, for parent particle masses of O(100) GeV, MATHUSLA could probe lifetimes close to the BBN limit.
We have also analyzed the improvement in sensitivity that would come from a forward LLP detector within a future 100 TeV hadron collider (FCC-hh), finding that such an improvement would be quite dramatic and would allow to probe parent LLPs with masses up to 10 TeV for DM masses in the MeV range, as well as reach the BBN limit for parent LLP masses below 600 GeV.
In addition to the above results, we have discussed in some detail several technical aspects of freeze-in DM scenarios: the first aspect, being particular to our Higgs-assisted freeze-in DM model, concerns the EW phase transition, which switches-on some DM production processes when it takes place and plays a key role in regularising certain 2 → 2 scattering processes which would otherwise make DM production UV-sensitive. At the same time, our calculation of the DM production prior to the EW phase transition includes 25 Obviously, a DV + E miss T search at FCC-hh would significantly increase the mass reach of the corresponding search at HL-LHC. Still, already for cτχ 2 1 m the forward detector is expected to be more sensitive than the main FCC-hh detector(s), as the former will have much less background and the geometric acceptances of both are comparable (as opposed to HL-LHC vs MATHUSLA, case in which the smaller geometric acceptance geometric of MATHUSLA results in MATHUSLA being more sensitive than the LHC DV + E miss T search only for cτχ 2 100 m).

JHEP03(2020)022
the corrections due to the thermal mass of the Higgs doublet and the parent particle. The inclusion of thermal masses in DM relic density computations is a rather generic issue that has so far been disregarded for freeze-in via decay scenarios 26 and can have a significant quantitative impact on the results. The other technical aspects analyzed also apply more generally to models of freeze-in DM production from decay, and correspond to the impact of 2 → 2 scattering processes on the DM relic abundance and the interplay between freeze-in and super-WIMP Dark Matter production mechanisms.
A 1 → 2 vs. 2 → 2 freeze-in processes As outlined in section 3.2, whenever there exists a decay A → B SM χ contributing to freezein DM production, there will be related 2 → 2 scattering processes also contributing to DM freeze-in at the same order in the freeze-in coupling y χ , where the SM particle B SM appears in the t-channel of the scattering process X SM + A → Y SM + χ, with both X SM and Y SM SM particles. This is shown in figure 11, with the SM coupling among X SM , Y SM and B SM denoted generically by g SM .
We define a measure of the relative importance of these 2 → 2 scattering processes w.r.t. the 1 → 2 decays A → B SM χ, simply given by the ratio of DM yields from 1 → 2 and 2 → 2 processes (with thermal corrections neglected in this discussion): Figure 11. Feynman diagrams for the processes A → B SM χ (left) and X SM + A → Y SM + χ (right).
with Γ A the decay width of the process A → B SM χ, σ(s) the cross section for the process X SM + A → Y SM + χ, m X and m A the mass of the initial states A and X SM and K 1 the first modified Bessel function of the 2 nd kind. In the last step of (A.1) we have performed explicitly the integration in x = m A /T (only possible when thermal corrections are neglected). As discussed above, Γ A , σ(s) ∝ y 2 χ and so the ratio (A.1) is independent of y χ and might turn out to be sizable. In the following, we compute the above ratio in our freeze-in scenario for B SM = h and X SM = Y SM = q, a SM quark with Yukawa coupling y q (other possibilities for X SM and Y SM such as SM gauge bosons W , Z would yield similar results in this case, being the nature of the t-channel mediator the relevant ingredient here). The cross section σ(s) has a simple analytic form in the limit m q , m χ → 0 Inserting (A.2) into (A.1) and with Γ A given by Γ(χ 2 → hχ 1 ) in eq. (3.4), we obtain the ratio Y 2→2 /(y 2 q Y 1→2 ) as a function of m A , shown in figure 12. We nevertheless note that for q being the top quark, which corresponds to the leading 2 → 2 scattering process (since σ(s) ∝ y 2 t in that case), the top mass m t cannot be neglected. In this case we obtain the cross section σ(s) numerically and show the ratio Y 2→2 /(y 2 t Y 1→2 ) as a function of m A in figure 12, keeping for simplicity m χ → 0. It becomes clear that the decay processes are largely dominant except in the limit m A → m h , and for m A < m h when 2-body decays are forbidden and the ratio (A.1) should instead be understood as Y 2→2 /Y 1→3 (the relevant decays of A are 3-body). In figure 12 we also show the region m A < m h down to m A = 100 GeV, where the ratio Y 2→2 /(y 2 q Y 1→3 ) is very large and 2 → 2 scatterings clearly dominate freeze-in DM production. We stress that 2 → 2 scattering processes with X SM and Y SM being SM states other than quarks (e.g. X SM = V SM (a SM gauge boson) and Y SM = h) may be more important, yet their contribution is expected to be qualitatively of the same order of magnitude as the former (see e.g. figure 5 of [60] in the context of the Type-I see-saw scenario). A complete, detailed analysis of the impact of 2 → 2 scattering processes in freeze-in scenarios is left for future work.
The case where B SM is an EW gauge boson nevertheless behaves in a different fashion to the above, and is worth discussing separately. Focusing on B SM = Z in the EW broken phase (the analysis for B SM = W ± is analogous), we note that the cross section for the JHEP03(2020)022 2 → 2 scattering process X SM +A → Y SM +χ (with X SM = Y SM = q, a SM quark) mediated by a t−channel Z boson is given in the limit m q , m χ → 0 by and scales as m −2 A (independent of s) for s m 2 A , m 2 Z , as opposed to (A.2) which scales as s −2 in that limit. Naively, this would lead to a UV divergent result for (A.1) (without considering thermal corrections, which are crucial in this case), signaling a UV dominated freeze-in production mechanism from 2 → 2 scattering processes mediated by a t−channel Z boson. We note however that the interaction among A, χ and the EW gauge bosons vanishes when the EW symmetry is restored above the EW phase transition temperature (recall the discussion in section 3.2), which regulates the integral ds in (A.1) and again yields a subdominant contribution from 2 → 2 scattering processes to the freeze-in DM abundance.
Finally, it is worth stressing that in a general scenario (such as ours) with various 1 → 2 and 2 → 2 processes contributing to DM freeze-in, the accurate expression analogous to (A.1) may not be as simple, with different contributions becoming relevant at different temperatures (particularly due to EW symmetry breaking in our scenario).

search recast
We describe here our recast of the ATLAS search for displaced vertices plus missing transverse energy [45] at √ s = 13 TeV with 32.8 fb −1 . The auxiliary material from the ATLAS search [97] provides efficiencies that can be applied to simulated truth level event samples. The search defines tracks (which will be associated to displaced vertices DV below) as: • The particle associated to the track is charged and stable.
• The particle has p T > 1 GeV.
• The particle has a transverse impact parameter d 0 > 2 mm.
Displaced vertices are constructed from these tracks. The DV must satisfy: • 4 mm < R < 300 mm, where R = x 2 + y 2 is the transverse distance to the interaction point, and |z| < 300 mm.
• The DV must have 5 or more tracks (as defined above) associated to it.
• The DV mass m DV must be at least 10 GeV, calculated from the three-momenta of the tracks (as defined above) associated to it and assuming a charged pion mass for each of the tracks.
Finally, the whole event is required to have: • E miss T > 200 GeV.
• 75% of the events should have at least one jet 27 with p T > 70 GeV or at least two jets with p T > 25 GeV.
• At least one displaced vertex (as defined above).
After our signal event simulation with Madgraph aMC@NLO [79] and Pythia8 [98], followed by jet-clustering using the FastJet [99] implementation within Delphes 3 [100], we perform the track, vertex, and event selection using the ROOT [102] output of Delphes 3. For those DV in events that pass all the cuts, we apply the vertex-level and event-level efficiencies provided within the auxiliary material of the ATLAS search, obtaining overall event selection efficiencies for ZZ + E miss T , Zh + E miss T and hh + E miss T signal decay channels, respectively. These overall efficiencies for the three separate channels are a function of the χ 2 mass and lifetime, and we therefore generate events for the three different different decay channels and different χ 2 masses and decay lengths 28 cτ χ 2 . 27 The ATLAS analysis requires trackless jets, where a trackless jet is one where the scalar sum of the charged particle pT is less than 5 GeV for those particles with small impact parameter with respect to the primary vertex. However, since the small impact parameter is not defined by ATLAS, we are unable to impose the trackless requirement in our recast. 28 In practice, we have used a fixed cτχ 2 = 0.1 m for event generation, noting that since each displacement variable is proportional to cτχ 2 , we can trivially rescale the position space coordinates of each particle to different values of cτχ 2 , before applying the cuts listed above. Our recast ATLAS result Figure 13. Excluded cross sections at 95% C.L. as a function of the gluino lifetime, for a long lived gluino decaying to a quark pair and a neutralino, compared to the published result from ATLAS [45]. The gluino mass is fixed to 2 TeV (top) or 1.4 TeV (bottom). The neutralino mass is fixed to 100 GeV.

JHEP03(2020)022
In order to validate the above analysis, we have applied our recast to a long-lived gluino model, where the gluino decays to a quark pair and a neutralino. This is the model used by ATLAS to interpret the results from their DV + E miss T search. We show our derived 95% C.L. exclusion limits, as compared to those shown by ATLAS [45], in figure 13 as a function of the gluino lifetime and in figure 14 as a function of the neutralino mass. We find overall good agreement, except in the compressed region (neutralino mass approaching the gluino mass) where our derived limits are stronger than those of ATLAS, signaling that our recasting procedure is less trustworthy for a compressed scenario. In our freeze-in DM model, this happens when the mass of χ 2 approaches the SM Higgs mass (the DM candidate χ 1 is in general approximately massless on LHC scales), such that our derived limits in that region will be stronger than the actual, would-be ATLAS limits. Our recast ATLAS result Figure 14. Excluded cross sections at 95% C.L. as a function of the neutralino mass, for a long lived gluino decaying to a quark pair and a neutralino, compared to the published result from ATLAS [45]. The gluino mass is fixed to 2 TeV (top) or 1.4 TeV (bottom). The gluino lifetime is fixed to 1 ns.
C Higgsino pair production cross sections at √ s = 13, 100 TeV In this appendix, we show the NLO+NLL pure Higgsino pair production cross sections (assuming decoupled squarks) used in this work, corresponding to pp →H 0