Interplay of super-WIMP and freeze-in production of dark matter

Non-thermalized dark matter is a cosmologically valid alternative to the paradigm of weakly interacting massive particles. For dark matter belonging to a $Z_2$-odd sector that contains in addition a thermalized mediator particle, dark matter production proceeds in general via both the freeze-in and superWIMP mechanism. We highlight their interplay and emphasize the connection to long-lived particles at colliders. For the explicit example of a colored t-channel mediator model we map out the entire accessible parameter space, cornered by bounds from the LHC, big bang nucleosynthesis and Lyman-alpha forest observations, respectively. We discuss prospects for the HL- and HE-LHC.


I. INTRODUCTION
The evidence for dark matter (DM) in our Universe provides a strong motivation for extending the standard model (SM) of particle physics by a dark sector containing a thermally or nonthermally produced relic. While the hypothesis of a thermalized and frozen-out DM candidate-such as a weakly interacting massive particle (WIMP)-is an attractive and thus widely studied possibility, it is by far not the only viable explanation. In particular, in view of many null-results from WIMP searches, an exploration of alternative scenarios is vital to exploit the current experimental capabilities and identify the nature of DM.
One such scenario is feebly interacting DM that never reaches thermal equilibrium with the SM throughout the cosmological history. In this case, DM production may proceed via occasional scatterings or decays of particles in the thermal bath [1][2][3][4][5], so-called freeze-in [6]. Another possibility is the out-of-equilibirum decay of a thermally decoupled mother particle, i.e., through the super-WIMP mechanism [7,8]. The latter is realized in models where the mother particle belongs to a Z 2 -odd dark sector, that forbids its decay into SM particles, while it may have sizeable couplings to the SMİn this case, the mother particle freezes out similarly to a WIMP while DM is produced through its decay, that typically becomes efficient much later in cosmic history. In addition, in general, a contribution to DM production from freeze-in is also present within this setup, as long as the mediator decay is possible [6].
In this article, we highlight the phenomenological implications of the interplay of super-WIMP and freezein production of DM and provide up-to-date experimental constraints and prospects. We consider a Z 2 -odd dark sector comprising a feebly interacting DM particle and a mediator that transforms nontrivially under the SM gauge groups, such that its gauge interactions drive it towards thermal equilibrium in the early Universe. In contrast, the feeble DM interactions prevent it from thermalizing. For concreteness, we focus on a Majorana fermion DM candidate and a colored t-channel mediator, mapping out the entire accessible parameter space.
The DM density constraint imposes a fairly general relation between the involved masses and the lifetime of the mediator [6]. In a wide range of the cosmologically allowed parameter space, the mediator has macroscopic decay length allowing for experimental tests for long-lived particles at colliders as well as effects on big bang nucleosynthesis (BBN) through late decaying mediators. They constrain the parameter space towards small mediator masses and small mass splittings, respectively. We discuss current constraints from searches for detector-stable R-hadrons as well as future projections for stable and metastable mediators at the HL-and HE-LHC. For large mass splittings and a significant super-WIMP contribution to DM production, large deviations of the DM momentum distribution from the thermal one can arise. This leads to a large free-streaming length, suppressing the amplitude of the matter power spectrum on small scales, which can be probed via Lyman-α forest observations [9][10][11]. This constrains the parameter space towards large mass splittings. The same observation constrains the parameter space towards very small DM masses (few keV) where freeze-in dominates. The parameter space is hence cornered from all sides.
The remainder of this work is structured as follows. We first introduce the model under consideration in Sec. II and refer to possible embeddings and variations. In Sec. III, we detail the computation of the DM density and provide some model-independent phenomenological considerations. Finally, Sec. IV provides results for the cosmologically viable parameter space, experimental constraints and future projections. We conclude in Sec. V.

II. THE MODEL
As a simple example of a Z 2 -odd new physics sector we consider a top-philic, colored scalar t-channel mediatort and a Majorana DM fermion χ interacting with the SM through the Lagrangian where D μ is the covariant derivative, t is the top quark Dirac field and λ χ is the new physics coupling. Thet particle is a SUð2Þ L singlet and has hypercharge identical to t R , similar to a right-handed squark field in supersymmetry. The model introduces the three parameters m χ , mt and λ χ . In this work, we focus on the regime of sufficiently small couplings λ χ , such that the χ particle never reaches thermal equilibrium with the SM bath. This means that neithert − χ conversions, such as (inverse) decayst ↔ tχ, nor annihilations, such as χχ ↔ tt, occur at rates comparable to the Hubble expansion rate throughout cosmic history. Details on further possible interactions and its phenomenology in the case of thermalized DM can be found in [12,13].
The simplified model, and variants with different spinand gauge-quantum numbers, can be part of generic extensions of the SM of particle physics. For example, the model possesses a natural embedding in supersymmetric models. In this case, the Z 2 -symmetry can be identified with R-parity, and the mediator with the lightest superpartner of the SM particles (being the right-handed stop for the specific model from above). The feebly interacting DM particle can be realized in the context of a hidden sector, that features an unbroken hidden Uð1Þ gauge symmetry. After supersymmetry breaking, a small kinetic mixing with the SM Uð1Þ Y hypercharge leads to a small bino-admixture of the hidden gaugino, providing a small coupling λ χ of the form introduced above [14,15].
Arguably, also supersymmetric models featuring gravitino DM and a long-lived next-to lightest supersymmetric particle (NLSP) share similarities with the type of models studied here if R-parity is conserved, but also exhibit differences due to the nonrenormalizeable interactions [1,2] (see [16] for a recent analysis of stop NLSP, and references therein for other possibilities).
A variant of the model considered here, but without Z 2 -symmetry, has been studied in [17,18].

III. DARK MATTER PRODUCTION
For small enough values of the coupling λ χ that connects DM to the SM, the DM particle χ is never in equilibrium with the SM thermal bath. In this case, any process throughout the cosmic history leading to the production of χ particles contributes to an accumulated χ population. Immediately after the end of inflation, χ particles may be produced during the reheating process. In this work, we assume that this process leads to a negligible contribution to the abundance of χ particles, and adopt the common assumption that reheating produces a thermal bath of SM particles, with maximal temperature given by T R . Furthermore, we assume T R ≫ mt, such that the mediatort thermalizes due to its gauge interactions. 1 In this case, within the simple model considered here, there are two distinct sources of χ particle production. First, the freeze-in mechanism that is most efficient for T ∼ mt, and second, the super-WIMP mechanism, corresponding to the late decay of the frozen-out population oft. In the following, we discuss both sources in turn.

A. Freeze-in
Freeze-in production relies on the occasional production of χ particles within a thermal bath. For the model considered here, due to the Z 2 -symmetry in the dark sector, production processes have to involvet in the initial or final state. Since the abundance oft becomes strongly suppressed for T ≪ mt, the relevant temperature range for freeze-in is T ≳ mt. At these temperatures, gauge interactions keept close to thermal equilibrium, i.e., we may assume nt ≃ n eq t ¼ gt We consider both the 1 → 2 processt → χt as well as all allowed 2 → 2 processes ab → χc, includingtt → χg, tg → χt, gt → χt Ã . The Boltzmann equation for the number density n χ reads [6] 1 For T R ≲ m˜t the relic density becomes dependent on T R and the production via freeze-in may dominantly proceed via DM pair production whose rate is suppressed by heavy mediator propagators arising in the t-channel or in loops. Hence, significantly larger couplings are expected to saturate the relic density constraint than found in this work. Here where q =ðE a E b Þ, g a are the internal degrees of freedom (d.o.f.) of species a, and we neglected the loss term, as appropriate for freeze-in, as well as statistical factors ð1 AE f i Þ (see below). The factor 2 in (2) takes into account charge conjugated processes, which contribute equally due to CP symmetry and the Majorana nature of χ. Since n χ appears only on the lefthand side, the Boltzmann equation can be solved by direct integration for the yield Y χ ¼ n χ =s, where s ¼ π 2 45 g ÃS T 3 is the entropy density. The final yield can be split into contributions from 1 → 2 and 2 → 2 processes. The former is, for example, given by where x ¼ mt=T, we assumed dg ÃS =dx ¼ 0 during freezein and introduced where K i denote modified Bessel functions. The integral saturates for x 0 ≳ 1 due to Boltzmann suppression of Y eq t . The contribution to the DM density from freeze-in is given by ðΩh 2 Þ fi ¼ m χ Y fi χ sðT 0 Þh 2 =ρ crit , where sðT 0 Þ and ρ crit denote the entropy-and critical energy density today, respectively.
We compute the freeze-in contribution with MICROMEGAS 5.0.4 [19] which assumes the mediator to follow the equilibrium density. As discussed above we expect this to be a good approximation for the setup considered here. We used the default approximate phase space integration as well as the optional full vegas integration routine, that includes also quantum statistical factors, and found deviations below 5%.
For m˜t > m t þ m χ , the two-body decayt → χt is kinematically allowed. Since 2 → 2 processes are formally suppressed by two powers of a SM coupling constant, one may expect them to give a subdominant contribution in that case. Nevertheless, we find them to contribute at the same level as the decay. This has several reasons. For processes such astt → χg,tg → χt, gt → χt Ã the relevant SM coupling is α s , which is sizeable for most of the parameter space. In addition, 2 → 2 processes are favored kinematically over 1 → 2 for T ≳ mt. Lastly, there is a large number of possible 2 → 2 processes that add up, while only a single 1 → 2 channel exists.
When both 1 → 2 and 2 → 2 processes are kinematically allowed, unphysical divergences may occur related to nearly on-shell propagators (the default routine within MICROMEGAS excludes 2 → 2 processes in that case). We checked that no such effects occur at a sizeable level, except fortg → χt andtγ → χt. For these processes the cross-section σðsÞ becomes enhanced close to the threshold ffiffi ffi s p ≳ mt. These processes feature an s-channel contribution, involving a propagator 1=ðs − m 2 t Þ that gives a large contribution close to threshold. A similar effect occurs for tZ → χt when mt ≫ m Z . The enhancement can be understood as soft initial state radiation that contributes to the next-to leading order correction to the decayt → χt, in an expansion in the SM couplings. We expect it to be regulated for the χ production rate when consistently including all real and virtual corrections. In addition, one may argue that for the processes above the enhancement is cut off when including a thermal mass for the soft initial state particle. Here we do not attempt to provide a full next-to leading order result. Instead, we implement a cut-off ffiffi ffi We checked that as long as R is close to unity the final yield depends only very weakly on the precise choice. We used R ¼ 1.2 in our numerical results. The total value of the final yield is affected at most at the 10% level for mt ≲ 10 5 GeV when choosing R ¼ 1.1 instead. In addition, we checked that when omitting the processestg → χt,tγ → χt andtZ → χt in the abundance calculation, all results remain qualitatively unchanged (see Sec. IVA for details).

B. Super-WIMP
The super-WIMP mechanism relies on the thermal freeze-out of the mediator, that subsequently decays into the DM particle [8]. Within the model considered here, freeze-out oft annihilation into SM particles yields an abundance Y fõ t at temperatures T ≪ mt=25, analogous to WIMP freeze-out, that can be converted into the density parameter ðΩh 2 Þ˜t ¼ mtY fõ t sðT 0 Þh 2 =ρ crit . At a much later time t ≃ 1=Γt →χt , the mediator decays, which yields a contribution to the DM density given by We compute the freeze-out abundance of the mediator in the absence of DM as a function of mt using MICROMEGAS 5.0.4 [19]. We take into account Sommerfeld enhancement of the mediator annihilation cross section as detailed in Appendix B of [20]. The mediator density may also be affected by bound-state effects [21][22][23], which we leave for future work.

C. Some model-independent considerations
Before discussing the parameter space and phenomenology of the specific model considered here in more detail, we highlight several properties that apply more generally to models containing a Z 2 -odd mediator with sizeable interactions with the SM, along with a DM particle that is very weakly coupled. These type of models are not constrained by WIMP searches for direct-and indirect detection. However, the presence of the mediator leads to testable signatures. It is in thermal equilibrium in the Early Universe, and can potentially be produced in laboratory experiments. The mediator needs to be heavier than DM, such that it can decay into the stable DM state. Due to the weak coupling, the mediator decay rate is suppressed, leading generically to long-lived particles with special implications for phenomenology.
In this case, both freeze-in and super-WIMP contributions are present in general. The former depends on the production rate, which in turn depends on the small DM coupling, in our case ðΩh 2 Þ fi χ ∝ λ 2 χ . Since, for small enough λ χ , this coupling plays no role for the mediator freezeout, the super-WIMP contribution is independent of the DM coupling, i.e., ðΩh 2 Þ sW χ ∝ λ 0 χ . Therefore, if we require that the total abundance matches the observed value, ðΩh 2 Þ fi χ ðλ χ Þ þ ðΩh 2 Þ sW χ ¼ 0.12, solutions can only exist for points in parameter space for which ðΩh 2 Þ sW χ ≤ 0.12. In that case, the condition above can be used to determine the value of λ χ to explain the measured DM density. We therefore expect in general that within a large portion of the parameter space freeze-in dominates ("bulk"). The viable region in parameter space is then bounded by a hypersurface on which the super-WIMP mechanism saturates the DM density constraint ("boundary"). We will see below that this expectation is borne out in the model considered here (cf. [24]).
Provided the process that corresponds to mediator decay gives a sizeable contribution to the freeze-in abundance (in the model considered here this is the case for m˜t > m t þ m χ , such that two-body decay is kinematically allowed), we can estimate the freeze-in abundance by generalizing (5) whereH ¼ H=m 2 med , and c ¼ 1ð2Þ for a neutral (charged) mediator. Within the model considered here, Γ med → Γt →χt , m med → mt, g med → gt, c → 2.
As long as the temperature of mediator freeze-out is well above the electroweak scale, the number of relativistic d.o.f. is approximately constant such that Y eq med ðxÞ andHðxÞ are functions of x only, without reference to m med . Furthermore, the number of internal d.o.f. of the mediator cancels out inside the integrand in Eq. (8). Hence, the integral in Eq. (8) is a constant. Consequently, within the "bulk" region of parameter space (for which ðΩh 2 Þ sW χ ≤ 0.12), where we negected the 2 → 2 contribution in order to obtain the parametric estimate above (cf. [6]). For a given decay rate Γ med , or equivalently mediator lifetime, this imposes a correlation m med ∝ m 1=2 χ between the mediator and DM mass. We expect this finding to be applicable to the general class of models discussed above, see e.g., [25].
In addition, within the "bulk" region of parameter space, for which freeze-in dominates, Eq. (9) can be used to estimate the time t dec ≃ Γ −1 med when the (subdominant) population of frozen-out mediator particles decays. In terms of temperature, and well above the electroweak scale, T dec ≃ 6 × 10 8 GeV × ðΓ med =GeVÞ 1=2 and hence x dec ¼ m med =T dec is given by That is, there is a relation m χ ∝ x 2 dec , which we again expect to apply to the class of models discussed above. This also shows that the super-WIMP production via mediator decay is well separated in time from the freeze-in regime x ≃ Oð1Þ for m χ ≫ 0.1 keV.
Note that Eq. (9) furthermore implies a modelindependent statement about the region in the mediator-DM mass plane that provides long-lived particles at the LHC. For proper decay length in the range [1 m; 1 mm] we find where the lower edge of the mass range corresponds to the upper edge of the decay length and vice versa while smaller masses provide mostly detector-stable mediators. Note that in case of additional contributions to DM production the lifetime becomes larger, shifting the respective mediator mass range to larger values. Within freeze-in scenarios long-lived particle signatures at the LHC were studied in [6,[26][27][28][29][30][31].

A. Parameter space and DM density
Out of the three free parameters m χ , mt; λ χ , one may be fixed by the condition that the sum of freeze-in and super-WIMP contributions to the χ density equals the observed DM abundance. We choose to fix λ χ by this condition, and use the DM mass m χ and the mass difference to describe the remaining two-dimensional parameter space. Dark matter stability requires Δm > 0.
In Fig. 1, we show the resulting coupling λ χ (green lines) as function of the DM mass, for fixed Δm ¼ 2.5 TeV (left) and Δm ¼ 100 TeV (right). The freeze-in contribution dominates for m χ ≪ m crit χ ≃ 1.6 TeV and 40 GeV in the two cases, respectively, corresponding to the "bulk" region discussed before. The coupling λ χ required to obtain the measured relic density increases towards lower DM masses. This can be understood in the following way: as discussed above, for Δm ≫ m χ , the freeze-in yield is approximately independent of m χ , such that When increasing m χ for fixed Δm, the super-WIMP contribution becomes larger. Since it is independent of λ χ , its value saturates the constraint ðΩh 2 Þ sW → 0.12 for some finite value of m χ → m crit χ . At this point λ χ → 0, and no solutions providing the measured DM abundance exist for larger values of m χ . In order to quantify the uncertainty due to the approximate treatment of 2 → 2 contributions with threshold enhancement, we show an error band in Fig. 1 around the green line. The lower boundary corresponds to the result obtained when using a cut parameter R ¼ 1.1 (see Sec. III A). For the upper boundary we omit the enhanced 2 → 2 processes when computing the freeze-in abundance.
The full two-dimensional parameter space is shown in Fig. 2, covering the entire accessible region of parameter space (left), going up to very large mediator masses, and the patch for Δm < 50 TeV, with m χ on a linear scale (right), respectively. The value of λ χ obtained from imposing the condition ðΩh 2 Þ fi þ ðΩh 2 Þ sW ¼ 0.12 is indicated by the green contour lines. For the left panel in Fig. 2, the contours show decades in log 10 λ χ , and for the right panel the value of λ χ normalized to 10 −12 . In this two-dimensional parameter space, the "boundary" corresponds to the thick black line, for which ðΩh 2 Þ sW → 0.12. Values of ðm χ ; ΔmÞ above that line are excluded due to DM overproduction. The region below the black line corresponds to the "bulk" as discussed before. Well inside the "bulk," freeze-in production dominates. The black dashed (right panel only) and dotted curves show the contour of constant relative super-WIMP contribution ðΩh 2 Þ sW χ =ðΩh 2 Þ tot χ ¼ 50% and 10%, respectively. Below the dotted curve the super-WIMP contribution is subdominant.
In the following, we discuss observational signatures that can probe different regions of the parameter space, including long-lived colored particles at the LHC (R-hadrons) and during BBN, as well as Lyman-α forest observations.

B. Collider constraints and projections
As discussed in Sec. III C the relic density constraint implies the existence of long-lived particles in a large part of the parameter space within the class of models considered here. In Figs. 1 and 2, we indicate the proper decay length by the cyan dotted contours for cτ ¼ 1 mm and 1 m, corresponding to the range in which decays typically take place inside the detector. Below the latter curve a significant fraction of mediators decay outside the detector. The colored mediatort can be copiously produced at hadron colliders. For large mediator decay lengths, cτ ≳ ðdetector sizeÞ, searches for detector-stable R-hadrons provide a promising discovery channel at the LHC.
Here we constrain the model by current searches at the 13 TeV LHC and estimate projections for the HL-and HE-LHC.
Current searches for detector-stable top-squarks with the CMS detector exclude masses up to m˜t ¼ 1250 GeV at 95% CL [32]. This limit is directly applicable to our model in the region where cτ ≫ 1 m. For intermediate lifetimes, cτ ≲ 1 m, relevant for DM masses m χ ≲ 100 keV, the limit is weakened due to the exponential suppression of the fraction of decays outside the detector. We use the reinterpretation of the above limit for finite lifetimes provided in [20] considering the 'generic model' for hadronization. The resulting 95% CL exclusion is shown in Fig. 2 (blue shaded region). For large DM masses the limit lies entirely in the detector-stable regime and its drop is simply caused by the chosen presentation in terms of m χ and Δm. Towards small masses Δm ≃ mt and the drop in the limit is due to the exponential suppression of the detector-stable fraction. Still, R-hadron searches constrain the parameter space towards small mediator masses down to the smallest m χ consistent with Lyman-α bounds, see Sec. IV D.
In order to illustrate the future sensitivity to the model, we consider R-hadron searches using 3 ab −1 at 14 TeV (HL-LHC) and 10 ab −1 at 27 TeV (HE-LHC). We compute the signal cross sections at the 14 and 27 TeV with NLLFAST [33] and PROSPINO [34], respectively. As the search is based on anomalous ionization loss and timeof-flight the signal efficiencies depend crucially on the velocity distribution of the produced mediators. To first approximation the velocity distribution stays unchanged for constant mt= ffiffi ffi s p . We therefore estimate the signal efficiencies by rescaling the ones from [32] (and [20] for finite lifetimes):  and analogous for 27 TeV. 2 We estimated the background by rescaling the one reported in [32] by the cross section ratio σ 14 TeV =σ 13 TeV (σ 27 TeV =σ 13 TeV for 27 TeV) computed with MADGRAPH5_AMC@NLO [37] for the leading background to heavy stable charged particles which is Drell-Yan production of muons.
In Fig. 2, we draw the corresponding projected 95% CL s -limits for the HL-(blue dashed) and HE-LHC (blue dotted). They reach mediator masses up to 2000 and 4050 GeV, respectively (see also the left panel of Fig. 2).  Note that a naive use of the recasting of the 8 or 13 TeV search for heavy stable charged particles [35,36] does not resemble this behavior up to ffiffi ffi s p ¼ 27 TeV, but results in significantly smaller signal efficiencies. The reason for this is the decreasing efficiency towards large p T for the current CMS detector. We hence implicitly assume an improved performance towards high p T for the HE-LHC. The latter can probe the entire DM mass range up to the boundary (where Δm ≃ 1500 GeV).
In addition to searches for detector stable objects, in the region m χ ≲ 100 keV, signatures of mediators decaying inside the tracker may provide further sensitivity. We expect searches for disappearing R-hadron tracks and displaced tops to be further promising discovery channels at the HL-and HE-LHC.

C. BBN bounds
The presence of a metastable colored mediator during the epoch of BBN affects the predictions for the primordial abundances of light elements through the energy release from its decay [38,39] as well as through bound-state formation with baryonic matter [40,41]. Due to strong hadro-dissociation processes the former effect is dominant for a hadronically decaying mediator. We estimate the respective constraints on the parameter space by applying the results from [38] for a hadronic branching ratio of 1 using the mediator freeze-out abundance Yt and lifetime as computed by MICROMEGAS 5.0.4. The slight dependence on the mediator mass is approximated by linearly interpolating (and extrapolating) the results for 100 GeV and 1 TeV in log-log space.
The resulting constraints are shown as the red shaded regions in Figs. 1 and 2. For fixed m χ and Δm, BBN imposes an upper bound on the lifetime, which translates into a lower bound on the coupling λ χ . For m χ ≪ Δm, both the lifetime and the mediator abundance Yt become independent of m χ , explaining the almost horizontal exclusion contour in Fig. 1. For Δm ¼ 2.5ð100Þ TeV we find λ χ ≥ 1ð6Þ × 10 −14 (see Fig. 1). The value of the coupling required for Ωh 2 ¼ 0.12 is consistent with BBN for most of the parameter space, except for a small strip close to the "boundary", at which λ χ → 0 (not resolved in Fig. 2), as well as the region Δm ≲ m t , for which the two-body decayt → χt is kinematically forbidden, such that the mediator lifetime is increased (see Fig. 2, right panel). BBN bounds are stronger than Lyman-α constraints (see below) for m χ ≫ keV and Δm ≲ 10 TeV.

D. Lyman-α forest bounds
In this section, we consider constraints on the DM and mediator mass from free-streaming of DM particles, that leads to a suppression of the amplitude of the matter power spectrum on length scales smaller than the free-streaming scale Here vðzÞ is the typical velocity of DM particles, and z prod is the redshift at which DM is (dominantly) produced.
The power on small scales can be probed by observations of absorption features in the spectra of distant light sources (quasars) imprinted by intervening clouds of neutral hydrogen, known as Lyman-α forest. The interpretation of these data depends on various properties of the intergalactic medium (including its redshift-dependent temperature and adiabatic index) as well as bias parameters that relate the hydrogen distribution to the underlying DM density field. These astrophysical effects are often described by a number of "nuisance" parameters, that need to be varied together with the cosmological parameters in order to obtain constraints from comparing theoretical predictions based on hydrodynamical simulations with observations [42,43].
Since a dedicated analysis is beyond the scope of this work, and in view of astrophysical uncertainties, we estimate Lyman-α constraints on the model considered here by computing the free-streaming length, and comparing to the maximally allowed value taken from [44]. More specifically, we translate the 2σ limits on the warm DM mass, as function of the warm DM fraction 0 ≤ f ≤ 1, into an f-dependent upper bound λ max fs ðfÞ, and then apply the latter to the model considered here (see below). For f ¼ 1ð0.2Þ the analysis of [44] yields m WDM ≥ 4.0ð1.5Þ keV, corresponding to λ max fs ð1Þ ¼ 0.10 Mpc and λ max fs ð0.2Þ ¼ 0.21 Mpc (we use cosmic parameters from Planck [45] for the conversion).
Even though, within the model considered here, DM is composed of a single particle species, the two populations produced via freeze-in and the super-WIMP mechanism, respectively, feature a different momentum distribution, and therefore different free-steaming lengths, denoted by λ fi fs and λ sW fs (see below). We denote the corresponding fractions of the DM density by f fi ¼ ðΩ χ h 2 Þ fi =0.12 and f sW ¼ ðΩ χ h 2 Þ fi =0.12, such that f fi þ f sW ¼ 1 in the cosmologically allowed parameter region. If for example λ fi fs approaches the maximal allowed value, we find that λ sW fs is negligibly small within the majority of the accessible parameter space, such that the fraction of "warm" DM is in this case f fi , while the rest behaves as cold DM on the relevant scales. The same is true vice versa. Therefore, we impose the bound as for both a ¼ fi, sW. This procedure is expected to fail when both production mechanisms produce a comparable freestreaming length, of the order of the maximally allowed values. In this case, the corresponding matter power spectrum can have a more complicated scale-dependence, which requires a dedicated analysis (see e.g., [46][47][48][49] for related discussions) which is beyond the scope of this work.
Let us now estimate the free-streaming lengths for both production mechanisms. We assume production is dominated at a redshift interval around some redshift z prod with a typical momentum p prod of the DM particles. Due to cosmic expansion, the momentum redshifts according to and the typical velocity entering (14) is given by Together with the standard expression this can be used to compute λ fs as function of z prod and p prod . We note that the integral over z is dominated by redshifts z ≳ 10 3 such that the value computed with lower integration boundary at z ¼ 0 and at the redshifts z ∼ 2-5 relevant for Lyman-α observations is practically identical. For the super-WIMP mechanism, i.e., late decayst → tχ of the mediator, we assume for the time of production where Γ decay is thet decay rate. The decay time can be converted to the redshift using and then solving T ¼ ðg ÃS ðTÞ=g ÃS ðT 0 ÞÞ 1=3 T 0 ð1 þ zÞ for the temperature, where g ÃðSÞ are the usual relativistic d.o.f. At the decay time, the kinetic energy of the mediator can be neglected, such that it decays at rest. The momentum is therefore given by The above equations amount to vγðz prod Þ ≃ mt=ð2m χ Þ for large mt where γ ¼ ð1 − v 2 Þ −1=2 . Note that for very long mediator lifetimes comparable to the time of recombination, the interaction of the mediator with baryonic matter can lead to further suppression of the power spectrum [50]. The mediator lifetimes we consider are much smaller than this, even though there may be a tiny region in parameter space very close to the "boundary", for which λ χ → 0, where this becomes relevant. We do not consider this possibility further here. For the freeze-in contribution the dominant part arises before the onset of a (strong) Boltzmann suppression of the mediator abundance, that is at T ∼ mt. We approximate the production redshift by where x fi corresponds to the value of x ¼ mt=T for which the production rate dY χ =dx is maximal. For simplicity, we consider production via decay only for which we find x fi ≃ 2.4. Production via scatterings peaks at a value of comparable size, and we therefore expect the choice adopted here to give a subleading contribution to the error budget. Note that scatterings are taken into account for the abundance computation. In order to estimate the typical momentum, we compute the average energy of DM particles produced from a thermal distribution oft. As before we focus on production via decay, yielding hE χ iðx fi Þ ≃ 0.92mt for Δm ≫ m t . The resulting bounds on the parameter space are shown by the purple-shaded regions in Figs. 1 and 2. The two distinct exclusion regions in Fig. 1 correspond to the cases λ fi fs > λ max fs ðf fi Þ and λ sW fs > λ max fs ðf sW Þ, respectively. In the former case, the free-streaming length becomes large due to the small DM mass, and in the latter case due to an interplay of the large mass splitting, a long mediator lifetime and a large super-WIMP fraction.
The resulting exclusion region on ðm χ ; ΔmÞ when imposing Ωh 2 ¼ 0.12 is shown in Fig. 2. For small DM masses and mt ≪ 10 7 TeV, the production is dominated by freeze-in (i.e., f fi ≃ 1). In this case, Lyman-α data impose a lower bound on the DM mass of the order of m χ ≳ 6 keV (left part of purple-shaded area in Fig. 2). This bound is comparable to warm DM bounds, but slightly stronger, because freeze-in produces a nonthermal spectrum with slightly higher momentum than in the corresponding thermal case. However, we conservatively attribute an uncertainty of a factor of two to the precise value, due to the approximate estimate based on the free-streaming length, as well as astrophysical uncertainties (see above).
For the region in parameter space for which the super-WIMP mechanism gives a sizeable contribution, Lyman-α data exclude large mediator masses, and require e.g., m˜t ≲ 10 6 TeV for m χ ≃ MeV (right part of purple-shaded area in Fig. 2). The decay of mediators above this bound would lead to a too large free-streaming length, because the heavy mediator converts its rest mass into kinetic energy of the DM particles. For mt ≳ 10 7 TeV this excludes points in parameter space with a fraction of DM produced via the super-WIMP mechanism down to f sW ≳ 10%.
The exclusion contour in Fig. 2 close to the crossing point of the two regions discussed above should be regarded as conservative, because both populations have sizeable free-streaming length in that part of parameter space. The approach outlined above tends to underestimate the suppression of the power spectrum in this case. As mentioned above, a proper treatment of this region goes beyond the scope of this work.

V. CONCLUSION
In this work, we studied a class of DM models comprising a Z 2 -odd dark sector that contains a feebly interacting DM particle along with a mediator that transforms nontrivially under the SM gauge group. The DM particle never thermalizes and is produced via a combination of freeze-in on the one hand, and late decays of frozenout mediator particles, known as super-WIMP mechanism, on the other hand. Despite the fact that DM interactions with the SM are tiny-well in agreement with null-searches in direct and indirect detection experiments-the presence of the mediator leads to characteristic signatures that can be probed by searches for long-lived particles at colliders, via Lyman-α forest observations, as well as the determination of primordial element abundances.
The interplay of freeze-in and super-WIMP production leads to a finite region in parameter space for which the observed DM abundance can be explained. Taking a Majorana DM particle χ and a colored top-philic scalar mediatort as an example, the DM mass is bounded by m χ ≤ m max χ ≃ 2700 GeV. In addition, there is a maximal possible mediator mass, depending on the value of m χ . For example, mt ≤ m max med ≃ ð6; 3 × 10 3 ; 2 × 10 6 Þ TeV for m χ ¼ ð10 3 ; 1; 10 −3 Þ GeV. In addition, we highlight a simple parametric relation between the mediator lifetime and the masses. For DM mass in the MeV range and mediator decay length on detector scales (mm-m), the mediator mass is in the (multi-)TeV region.
We explore experimental probes within the entire accessible parameter space of DM and mediator masses, and provide up-to-date exclusion limits. We find that the parameter space is constrained from all sides. Mediator masses around the TeV-scale and below are excluded from R-hadron searches at the 13 TeV LHC for m χ ≲ 1 TeV, and by BBN for m χ ≳ 1 TeV. On the other hand, very heavy mediators with mass above 10 6 -10 7 TeV are in conflict with recent Lyman-α forest observations, as are DM masses below about 6 keV. In addition, a tiny region close to the boundary of the accessible parameter space, with mt ≲ m max t , is excluded by BBN and Lyman-α observations for mt below and above 10 TeV, respectively (not fully resolved in Fig. 2). In this region, a large fraction of the DM abundance is produced via the super-WIMP mechanism, and the mediator lifetime becomes particularly large. We also provide projections for high luminosity and 27 TeV high energy upgrades of the LHC, that are sensitive to mediator masses up to around 2 and 4 TeV, respectively.
While we focus our phenomenological analysis in this work on a specific simplified model, we emphasize that the qualitative features remain the same in general. For example, for a lepto-philic mediator m max χ and m max med would be smaller due to its larger freeze-out abundance, while the collider bound would be weaker due to the reduced production cross section. In the future, it would be interesting to perform a dedicated analysis of Lyman-α forest bounds for combined super-WIMP and freeze-in production.