WIMPs Without Weakness: Generalized Mass Window with Entropy Injection

We study general freeze-out scenarios where an arbitrary number of initial and final dark matter particles participate in the number-changing freeze-out interaction. We consider a simple sector with two particle species undergoing such a thermal freeze-out; one of the relics is stable and gives rise to the dark matter today, while the other one decays to the Standard Model, injecting significant entropy into the thermal bath that dilutes the dark matter abundance. We show that this setup can lead to a stable relic population that reproduces the observed dark matter abundance without requiring weak scale masses or couplings. The final dark matter abundance is estimated analytically. We carry out this calculation for arbitrary temperature dependence in the freeze-out process and identify the viable dark matter mass and cross section ranges that explain the observed dark matter abundance. This setup can be used to open parameter space for both heavy (above the unitarity bound) or light (sub-GeV) dark matter candidates. We point out that the best strategy for probing most parts of our parameter space is to look for signatures of an early matter-dominant epoch.


I. INTRODUCTION
The particle nature and non-gravitational interactions of dark matter (DM) are major mysteries of particle physics and cosmology. A central classifying feature of a given dark matter scenario is its production mechanism. While attractive alternative scenarios, such as freezein [1], have been put forward recently, thermal production scenarios via DM freeze-out remain well-motivated. Perhaps the best-studied production mechanism is the classic thermal freeze-out via 2-to-2 annihilations for Weakly-Interacting Massive Particles (WIMPs) [2][3][4]. The dark matter yield is given in this case by Y DM ∼ x f.o. /(m DM M Pl σv ), valid as an order of magnitude estimate, as long as the dark sector was in thermal contact with the SM initially [5]. Given the measured temperature of matter-radiation equality T eq , the annihilation rate yielding the correct late time dark matter abundance can be determined by setting T eq ∼ Y DM m DM , implying that: This result has three profound implications. First, if σv remains close to constant after freeze-out, there is a fixed target cross section for indirect detection searches. Second, perturbative unitarity provides an upper bound on the DM particle mass. Third, given the naive cross section scaling of σv ∼ α 2 /m 2 DM , one expects a relationship between the DM mass scale and the interaction * Email: pasadi@mit.edu; † Email: tslatyer@mit.edu; ‡ Email: juri.smirnov@fysik.su.se; ORCID: 0000-0002-3082-0929 Top (DM freeze-out during a RD epoch): First, the unstable relic φ freezes out in a RD epoch at T = T φ , followed by the DM candidate χ freezing out at T = Tχ. Then φ takes over the energy budget of the universe at T = Ti and MD begins. At T = Tτ φ , φ particles decay to SM, diluting the χ abundance. The decay of φ should take place before TBBN ∼ 3 MeV given the strong bounds [6][7][8][9]. Bottom (DM freezeout during a MD epoch): Similar timeline as above, but now Tχ Ti, i.e. the universe enters the early MD epoch before χ freezes out, which affects its final abundance.
strength of the theory. When this scaling holds, it follows that in this simple scenario: where M Pl is the Planck mass, and α is a coupling constant assumed to be comparable to the coupling of Standard Model (SM) weak interactions. This relation is often referred to as the "WIMP miracle"; it describes the observation that a weak scale DM mass requires an interaction strength similar to the electroweak interaction.
The interaction topology of the freeze-out process can be generalized to the case where DM depletion is domi-nated by number-changing interactions with p DM particles in the initial state and q DM particles in the final state, with p > q; see e.g. Ref. [10] for a 3 → 2 freeze-out scenario known as Strongly Interacting Massive Particle (SIMP) freeze-out. Parametrizing the interaction rate factor σv p−1 ∼ α p /m 3p−4 DM , requiring the correct relic density relates the mass and interaction strength as m χ ∼ α T eq M Pl T eq 1/p ∼ α 10 28/p−9 GeV .
For p > 2 the DM mass is constrained to be in the sub-GeV range; for such a light DM relic, strong astrophysical/cosmological bounds exist due to DM self interaction constraints [11] and the verified predictions of Big Bang Nucleosynthesis (BBN) [6][7][8].
Furthermore, if multiple dark-sector species are produced by similar freeze-out mechanisms, but one or more species decays with a timescale shorter than the age of the universe, then the decays of the unstable components can inject significant entropy into the visible sector and modify the yield of any remaining dark components. We consider a simple two-state version of this scenario, where an already frozen out relic χ is diluted by entropy injection from another relic φ (following the notation of Ref. [12]). Possible timelines for this scenario are shown in Fig. 1. To obtain a non-negligible entropy injection, φ should be both abundant and long-lived enough to dominate the energy budget of the universe, giving rise to an early matter-dominant (MD) epoch starting at temperature T i [13].
While φ always freezes out during a radiationdominant (RD) epoch, for some points of the parameter space χ's freeze-out can happen after φ dominates the energy density of the universe, i.e. during a MD epoch (bottom figure in Fig. 1); see Ref. [14] for a previous work on DM freeze-out during an early MD epoch.
We study the freeze-out of both χ and φ from the thermal bath when they are controlled by general interactions beyond the conventional WIMP and SIMP scenarios. These freeze-outs happen, respectively, at temperatures T χ and T φ . We analytically calculate the final abundance of χ after the entropy injection and will show that such a minimal extension can vastly expand the viable parameter space. We also check our analytic formulas against numerical results for a few benchmark scenarios and find better than 10% agreement in the final relic abundance calculation.
Given the vastly different search strategies required to look for DM in different parts of the parameter space, identifying new dynamics that can expand the viable parameter space can have profound implications for the experimental efforts. We find that our minimal extension can substantially affect the viable mass range of thermal relics. We will argue that the best way to probe the high mass end of our parameter space is to look for signatures of an early MD epoch, as discussed in Refs. [15][16][17].
The effects of entropy injection in the early universe have been considered previously in the literature. It has been shown that in a setup with DM particles freezing out during a RD epoch and with 2-to-2 interactions controlling the freeze-out, the DM mass can be as large as 10 10 GeV [12]. If the baryon asymmetry of the universe (BAU) can be generated after the dilution, even higher masses are permitted. The effect of such a dilution on non-thermal [18] or hot thermal relics [19] has also been studied in the literature. Entropy injection has been studied in setups where the DM and SM have different temperatures and where the 2-to-2 interactions controlling DM freeze-out are χχ → φφ [20]. The decay of inflatons to decoupled SM and DM particles (without freeze-out or freeze-in) has been studied in Ref. [21].
In this work, we extend these studies by analytically studying general number-changing interactions for both φ and χ freeze-outs, including the case with p 4 initial particles. We also consider general temperature dependence in the freeze-out cross sections and the effect of χ freezing out during an early MD epoch. We will use the result of our analytic calculations to demonstrate the following: • The setup naturally opens up parameter space between the freeze-in and freeze-out DM annihilation cross sections.
• Lower interaction rates arise due to small couplings between DM and SM particles for DM masses at or below the GeV scale. We thus find new open parameter space for sub-GeV DM masses, which can be tested in the future by direct detection experiments with sensitivity to very-low-energy recoils, and improved indirect detection searches.
• Even with φ lifetimes as short as the electroweak phase transition timescale, it is possible to have a DM candidate with a mass that exceeds the WIMP unitary bound [22] of ∼ 100 TeV, and with unexpectedly large coupling values α ∼ O(1).
• In scenarios with 3 → 2 freeze-out interactions we find new parameter space with masses far above the GeV scale.
• Finally, we show that in general DM freeze-out via interactions with p 4 are viable even for heavy DM relics.
In Sec. II we summarize the result of our analytic calculation for entropy injection and general freeze-out mechanisms. In Sec. III we use our analytic calculation to identify the viable range of cross sections and masses that give rise to the observed DM abundance today for a few different benchmark freeze-out interactions. In Sec. IV we show the parameter space that entropy injection opens in a sample minimal model with a Z mediator. We conclude in Sec. V. Further details on the entropy injection calculation, freeze-out during a RD epoch, and freeze-out during a MD epoch are provided in Apps. A-C, respectively. We also comment more on a few other benchmark freeze-out scenarios in App. D.

A. Entropy injection
Consider a minimal DM setup with a single DM candidate χ, with mass m χ , that freezes out. Unless specified otherwise, we will take χ to be self-conjugate (e.g. a Majorana fermion or real scalar). Throughout this work, we assume that kinetic equilibrium is maintained between the dark sector particles and the SM. If the interaction controlling the freeze-out has p χ = 2 DM particles in the initial, q χ = 0 DM particles, and two SM particles in the final state, the relic abundance will be given by [23] Here s 0 and ρ c are today's entropy density and critical energy density, respectively, and the ratio has a numerical value of [24] s x f.o.,χ = m χ /T χ with T χ being the freeze-out temperature, g * (S),χ is number of relativistic degrees of freedom (for entropy) at T = T χ , M Pl ≈ 1.2 × 10 19 GeV is the Planck mass, and we used σv ≡ σ 0,pχ as the cross section for the interaction controlling the freeze-out, i.e. freeze-out interactions do not have any temperature dependence. Using this, it has been argued that the natural mass range for a DM particle freezing out via p χ = 2 interactions is around a TeV, with an upper bound of ∼ 10 2 TeV on its mass owing to unitarity arguments [22]. However, there are a number of ways to break outside this mass window. One simple possibility is when, in addition to DM and its portal to the SM, there is another relic φ with mass m φ that can freeze out and at some point decay to SM particles, injecting a significant amount of entropy into the SM thermal bath. For this to happen, the unstable relic φ should take over the energy budget of the universe after its freeze-out [23], giving rise to an early MD epoch. If such an epoch starts at T = T i , it can be shown that the injected entropy will dilute the DM relic abundance by a factor: where S f (S i ) is the SM bath total entropy after (before) the entropy injection, Γ φ is φ's decay rate, and g 1/3 * is a weighted average over relativistic degrees of freedom throughout the MD epoch; see App. A for further details. We assume φ can only decay to SM bath particles. The effect of this entropy injection on the DM relic abundance with p χ = 2, q χ = 0 (conventional WIMP) has been studied in Ref. [12].
We extend the work of Ref. [12] to the general case of arbitrary p χ and q χ , and the assumption that the abundance of φ itself is set by a freeze-out process. Furthermore, while the dilution factor ξ in Ref. [12] was treated as a free parameter itself, we consider scenarios where φ freezes out via a process with p φ and q φ particles in the initial and final states respectively. We also consider a general temperature dependence for the freeze-out cross section In the conventional case of p φ = 2 interactions, l φ corresponds to the partial wave expansion of the interaction in the non-relativistic regime (l φ = 0 for s-wave, l φ = 1 for p-wave, etc.); for higher p φ interactions the situation is somewhat more complicated. For simplicity, and to further capture the effect of other phenomena such as Sommerfeld enhancement [25][26][27][28][29][30], we remain agnostic about the interpretation of l φ and merely use it to capture the temperature-dependence of the freeze-out interaction. In App. B we show that where x f.o.,φ = m φ /T φ with T φ being φ's freeze-out temperature, g * (S)i is the relativistic degrees of freedom (for entropy) at T = T i , and g * (S),φ is the relativistic degrees of freedom (for entropy) at φ's freeze-out. In our setup, where all the particles are at the same temperature, we can simplify the equation further with g * = g * S at all temperatures. In deriving this equation we also assumed φ is identical to its anti-particle. Next, we should study the freeze-out of χ, distinguishing between the RD and MD epochs in the two following sections.

B. Freeze-out during a RD epoch
We consider interactions where p χ DM particles in the initial state convert to q χ DM particles (p χ q χ + 1) and an arbitrary number of SM particles in the final state. We also assume for both χ and φ particles and anti-particles are identical; this affects the symmetry factors used in the equations.
The Boltzmann equation for such relics will be where p = p χ , p φ (q = q χ , q φ ) refers to the initial (final) number of χ or φ particles involved in their respective freeze-out processes, s is the SM entropy density, Y is the relic's yield, and σv p−1 is the relevant cross section. Notice that the number of SM particles involved in the interaction does not enter the equation. Solving this equation for relic χ, we find (see App. B for further details): where the interaction cross section is defined as with σ 0,pχ a prefactor, x χ = m χ /T , and l χ capturing the temperature dependence of the interaction. (A slightly more general result can be found in App. B, see Eq. (B13).) We can easily check that in the limit of p χ = 2, q χ = 0, Eq. (10) reduces to Eq. (4). As a cross check, one can verify that this expression has the right dimensions; in doing so, it should be noted that σ 0,pχ has mass dimension −3p χ + 4. A similar calculation can be repeated for the freeze-out of φ and its abundance before its decay. We checked Eq. (10) explicitly against numerical calculations for a few benchmark scenarios and found better than 10% agreement.
Combining Eq. (10) with Eqs. (6)-(8), we find an expression for the χ relic abundance after the entropy injection The following assumptions were used to arrive at this equation: 1. both φ and χ freeze out during a RD epoch, 2. the entropy injection happens after χ has frozen out from the thermal bath, and 3. we neglected the factor of 1 in Eq. (6) for ξ, i.e. we assumed ξ 1.
Furthermore, as argued in Ref. [12], if the BAU is generated before φ decays, we should check that ξ > ∼ η b with η b = 6 × 10 −10 [24] denoting the observed baryon-tophoton ratio; otherwise any pre-existing asymmetry will be washed out to values too small to explain today's BAU. However, if a baryogenesis mechanism is provisioned for after the entropy injection, there is no obstruction to going to smaller values of ξ.
Equation (12) shows that Ω χ is most sensitive to p φ,χ (since they appear in powers of various quantities), masses, freeze-out cross sections, and φ's decay rate to SM Γ φ ; other parameters mostly affect the relic abundance by O(1) factors. Unlike the case of p χ = 2 we find that the final abundance does explicitly depend on m χ for p χ 3. We also find that the final abundance always increases as σ 0,pχ decreases. We can also calculate x f.o.,χ as a function of other parameters, see App. B.
Notice that Γ φ is a free parameter of the calculation at this level. By varying this parameter, it is possible to go to different regions of the parameter space with different DM masses. In particular, the entropy injection dilutes the DM abundance, which can be used to open up new parameter space for DM masses beyond the unitarity bound [22]. We will discuss this in more detail in the upcoming section.

C. Freeze-out during a MD epoch
When there is a large hierarchy between the masses of the two relics, φ takes over the energy budget of the universe and the universe enters an early MD epoch before χ freezes out, i.e. T χ T i , see the bottom row of Fig. 1.
Here φ's freeze-out goes forward as before. We can check from the calculation of App. B that its asymptotic yield Y φ ∞ will be given by and its energy density can be written as where s(T ) is the SM entropy at temperature T . In the second equality we neglect the change in g * S after the φ freeze-out and use Eq. (8). This was done to simplify the temperature-dependence of ρ φ (T ) so as to simplify the upcoming integrations. We will see eventually that as long as we do not introduce any new degrees of freedom this will only introduce O(1) change in the final DM abundance expression.
In App. C we provide more details on the Boltzmann equations governing χ freeze-out during this MD epoch. The final result for the χ abundance before the entropy injection is Comparing this with Eq. (10) we observe changes in the power of x f.o. and DM mass, as well as the appearance of T i in the calculation.
We checked this equation against the numerical result as well and found very good agreement. Using Eq. (6) and Eq. (8) for the entropy injection, we find an expression for the final DM abundance today: We should reiterate that in deriving this equation we neglected the evolution of g * S between T φ and T χ for analytic clarity; this will merely introduce O(1) effects on the DM abundance calculation. This dependence should be restored in a fully numerical treatment. It should also be noted that T i introduces some further dependence on m φ and σ 0,φ , see Eq. (8).

III. VIABLE DM PARAMETER SPACE
We can use Eqs. (12) and (16) to study the available parameter space in models with different freeze-out dynamics. In this section, we consider a few benchmark scenarios and comment on the new parameter space opened thanks to the entropy injection.
A. Broadening the thermal cross section range Let us start by considering the WIMP-like freeze-out case of p χ,φ = 2, q χ,φ = 0. In this setup, if the freezeouts happen during a RD epoch the final relic abundance does not explicitly depend on the masses. Without any entropy injection the observed DM abundance would require In our scenario, however, no absolute cross section or mass scale is singled out. The final abundance is instead determined by the ratio of the χ and φ annihilation cross sections. After the entropy injection, the right relic abundance is obtained even with cross section values below the standard WIMP thermal cross section.
In Fig. 2 we show the annihilation cross section values for χ and φ that yield the correct DM abundance, at different lifetimes for φ. We use l χ,φ = 0 (s-wave interactions) in generating this plot, and assume both relics freeze out during a RD epoch. We assume a geometric s-wave cross section, This allows us to relate each σ χ v to a maximal DM mass. This should be viewed as an approximate upper limit for the s-wave cross section during freezeout, for a given m χ , since the unitarity bound on the non-relativistic s-wave cross section is σv < 4π/k 2 mχ < 4  10 4 TeV mχ < 10 5 TeV mχ < 4  10 6 TeV mχ < 10 8 TeV mχ < 4  10 9 TeV The viable range of cross sections for a 2 → 0 freezeout process for both φ and χ. Using the unitarity bound the cross sections can be associated with an upper bound on the mass of the particle, denoted by dashed gray lines. On the red lines (below the red lines) we can get the observed relic abundance (overclose the universe) by assuming the denoted lifetime for φ. In the upper left corner of the parameter space the dilution factors will always be too small (for Tτ φ Tχ) and we underclose the universe. For the points in the lower gray region, the dilution factor will be small enough that the BAU is completely washed away at or after BBN. For τ φ ∼ 0.1 s there is no time left to reproduce this asymmetry before BBN; yet, for smaller lifetimes the BAU can be generated after φ's decay. and freeze-out occurs when the DM is only mildly nonrelativistic (i.e. k is parametetrically similar to m χ ). The m χ values obtained by this relation can thus also be viewed as approximate upper bounds on the DM mass for a fixed cross section.
An additional factor to consider is the baryon to photon ratio of the universe. If the dilution factor ξ becomes smaller than around η b ≈ 6 × 10 −10 , any previouslyproduced baryon asymmetry dilutes to values smaller than today's observed asymmetry, unless φ decays generate a SM baryon asymmetry. Thus, for this part of the parameter space, a baryogenesis mechanism should be provisioned that is active after the φ entropy injection.
In general we will consider lifetimes τ φ 0.1 s; this choice ensures the bulk of the energy density stored in φ has decayed by BBN. A more careful treatment of BBN constraints may open a small parameter space for larger lifetimes, see Refs. [6,7]; here we use τ φ 0.1 as a conservative bound and leave a more rigorous study of the BBN bounds for future works.
For the case of τ φ ∼ 0.1 second, the region with ξ η b is ruled out since in that region the observed BAU can not be produced before BBN. However, proposals exist for baryogenesis at temperatures as low as O(MeV) [31][32][33][34]; using these models the parameter space with ξ η b and τ φ 0.1 second can still be viable.
The lower bound on the DM cross sections in Fig. 2 for which our analysis is viable comes from the condition that the dark sector has to be thermally populated at early times, assuming a high reheating temperature. As we go to lower cross sections, eventually the two sectors will not start from equilibrium and we transition to a freeze-in [1] scenario instead. Our analytic analysis can be repeated for these cross sections; we leave this study for future work.
The existence of entropy injection and DM dilution allows DM annihilation cross sections to SM to be below the usual thermal WIMP cross section of Eq. (17). This allows us to have viable models with annihilation cross sections between the conventional freeze-out cross section and the freeze-in cross section. In the heavy DM mass end, this will open up parameter space for DM masses beyond the unitarity bound. At the lower end of the DM mass spectrum, the lower annihilation rates to SM particles lead to open parameter space that is accessible to upcoming experiments, see Sec. IV for an example of such models.

B. Generalized freeze-out mass window
The entropy injection shifts the required interaction rates for the DM candidates. Using a cross section parametrization, as well as unitarity arguments, this can be translated into a generalized viable DM mass window.
In Fig. 3, we show the available parameter space on the plane of m χ − m φ for different benchmark scenarios. We again use geometric cross sections as estimated upper bounds for the cross sections during freeze-out. Specifically, we take the non-relativistic generalized unitarity bound computed in Eq. (10) of Ref. [35] and replace the particle momentum with the mass, k → m: As in the s-wave case discussed above, this is justified by the observation that the DM is only mildly nonrelativistic during freeze-out; however, if the unitarity bound is truly saturated then the cross sections could be slightly larger than the values given here (and of course much smaller cross sections are always allowed).
In App. D we show variations of these plots under rescaling of the assumed maximal cross sections. Since the geometric cross section for χ annihilation should be close to the maximal value allowed by unitarity, using it will allow us to estimate the maximum viable DM mass in each scenario; this is the meaning of the m χ axis in Fig. 3. If we go to low enough DM masses, the right relic abundance can be obtained without an entropy injection (and with smaller freeze-out cross sections). This is indicated by the gray region in the plots. In all these plots we set the φ lifetime for every point in the param-eter space such that the final Ω χ matches the observed value [36].
In each plot, the part with τ φ lifetime longer than ∼ 0.1 second is ruled out by the BBN bounds [6,7]. For large enough hierarchies between m φ and m χ , even if the φ particles decay right after χ freezes out, they still inject too much entropy into the SM bath and dilute the DM abundance to below the observed value; thus, in this part of the parameter space, our model can account only for a fraction of DM today.
Clearly, depending on the type of interaction controlling the freeze-outs, the viable mass range shifts significantly. In particular, when p χ = 2 (top plot in Fig. 3), without any entropy injection we find the right DM abundance for m χ ∼ 100 TeV (note that the cross sections used in this plot saturate the unitarity bound, thus this agrees with the results of Ref. [22]). With entropy injection, we can have viable parameter space with m χ > ∼ 100 TeV. For the case of p χ = 3 this bound on the mass in the case of zero entropy injection is shifted to m χ ∼ 0.1 GeV (consistent with Ref. [10]). For m φ 10 18 GeV, in the p χ = 2 (p χ = 3) case the DM mass can go up to even ∼ 10 14 GeV (∼ 10 12 GeV), provided the unstable relic φ is long-lived enough.
It is worth reiterating that the final parameter space is not especially sensitive to q χ,φ as these quantities merely give rise to O(1) changes in the final Ω χ . The plots in Fig. 3 are made assuming l χ,φ = 0 interactions; we find that changing these parameters does not affect the viable mass window perceptibly either. Results of Sec. II can be used to remake these plots with smaller cross sections or different freeze-out interactions for χ or φ as well. For completeness, in App. D we include some plots on the viable mass range in a handful of other freeze-out scenarios.
The bottom plot in Fig. 3 shows that the entropy injection allows heavy DM masses even for p χ 4 freeze-out interactions. Given this enormous viable mass range, it will be interesting to find actual models in which the freeze-out is controlled by p χ 4 interactions. For instance, this can be achieved if we assume the process controlling the freeze-out respects a large discrete symmetry, e.g. Z N with large enough N , between DM particles.
We also assumed there are unspecified 2-to-2 elastic interactions between the SM and each of the relics φ and χ that keep them in kinetic equilibrium with the thermal bath. The crossing symmetry can relate this elastic process to a number changing process with p χ,φ = 2 that can affect the freeze-out dynamics. In Ref. [10], however, it was shown that for masses above 1 MeV there can always exist a range of couplings for which kinetic equilibrium can be maintained during the freeze-out via a 2-to-2 elastic scattering while the associated p χ = 2 annihilation process is suppressed.
All in all, Fig. 3 shows a vast mass range that can explain the observed DM abundance which, above all parameters, is most sensitive to p χ,φ . In particular, in such a simple setup the unstable particle φ can decay even be- χ��→�� ϕ��→�� � χ�ϕ →�� σ→σ ����   FIG. 3: The viable mass range of φ and χ assuming the interactions controlling the χ freeze-out are 2 → 0 (top), 3 → 2 (middle), and 4 → 3 (bottom), with φ freezing out via 2 → 0 interactions, and no temperature dependence in the freeze-out cross sections. The freeze-out cross sections of both χ and φ are set to the maximum value allowed by the unitarity bound for each mass point. We set the lifetime of φ such that we get the right relic abundance for every point in the parameter space. Contours of a constant lifetime (in seconds) are denoted by dashed green lines. Every point in the orange region requires τ φ > ∼ 0.1 second to give rise to the observed DM abundance and is ruled out by the BBN bounds. In the gray region, the entropy injection is negligible; the observed DM abundance can be obtained here without any entropy injection and smaller cross sections. On the black line, we get the correct relic abundance with cross sections that saturate the unitarity bound [35] and without any entropy injection. In the purple region, if the entropy injection happens after the χ freeze-out, we will always underproduce DM. For every point to the right of the red line the BAU has to be generated after the φ entropy injection since in that region ξ η b ∼ 6×10 −10 .
fore the electroweak phase transition, i.e. τ φ 10 −11 s, and still inject a significant amount of entropy into the SM bath. With such a short lifetime, even if the dilution factor is small enough to wash away any pre-existing baryon asymmetry (ξ η b ), we can regenerate BAU via electroweak baryogenesis proposals [37].
Our result shows that even in minimal scenarios including an entropy injection, there is essentially no preferred mass window for thermal relic DM. Such a broad mass window became viable thanks to the entropy injection, which in turn takes place if and only if we have an early MD epoch in the universe. Looking for other signatures of such an epoch could be the most model-independent way to probe the high mass end of our viable parameter space, see for example Refs. [15][16][17].

IV. VIABLE PARAMETER SPACE IN A BENCHMARK MODEL
As a concrete example of how entropy injection expands the parameter space of simple thermal relic DM models, in this section, we study a gauged B − L symmetry with some entropy injection at around τ φ ≈ 0.1 seconds. We do not consider the freeze-out of φ for simplicity and treat the dilution factor ξ as a free parameter.
We assume the Z mediator of the U (1) B−L gauge gets a mass via the Stuckelberg mechanism [38]. A new vector-like fermion χ with a B − L charge n χ = 1 is then accidentally stable, and provides a minimal DM candidate [39,40]. The dark matter particle χ freezes out via usual p χ = 2 interactions. The relevant interactions are where n χ (n f ) is the DM (SM fermions) charge under U (1) B−L , and g B−L is the gauge coupling of the U (1) B−L gauge group. The charges of the SM particles are −1 for leptons and 1/3 for quarks, such that nucleons couple with a B − L charge n N = 1. The DM charge will be taken as n χ = 3 for concreteness. In Fig. 4 and Fig. 5, we show the plane of dark matter mass vs. spin-independent elastic cross section. We study both limits of heavy (above the ∼ 100 TeV unitarity bound) and light (sub-GeV) DM masses. The U (1) B−L gauge coupling is uniquely determined for any point on this plot. We then choose the right dilution factor ξ for each point such that the right DM abundance today is obtained; contours of required ξ are shown.
In Fig. 4, we show the heavy DM scenario with a fixed mediator mass. The indirect detection limit from the FermiLAT satellite [41] is only relevant for the case that no dilution is present, i.e. ξ = 1. However, as discussed in Refs. [15][16][17], the early MD epoch can lead to increased production of high-density DM mini-halos. This effect could lead to an enhanced DM annihilation signal, if those halos are not disrupted at later times. In particular in the case of a low reheating temperature, this effect could improve the indirect detection prospects by several orders of magnitude. The LEP collider bound of m Z /g B−L 6.9 TeV [42,43] provides an upper bound on the spin-independent cross section σ DM−SM 5 · 10 −44 n 2 χ cm 2 . The current Xenon1T limit [44] excludes the no-dilution scenario in a wide range of masses as well.
For simplicity, in making this plot we assumed χ freezes out during a RD epoch, i.e. T i T χ . Since the entropy injection should occur before t ∼ 0.1 second, i.e. Γ φ (0.1 second) −1 , we can calculate a lower bound on ξ using Eq. (6), see the orange region in Fig. 4. Combining this limit with the collider searches mentioned above puts an upper bound of m χ 10 10 GeV on DM mass in this benchmark model.
With this example we demonstrate that for ξ 1 a wide parameter space remains open for DM masses well above the perturbative unitarity bound, and significant gauge interactions g B−L ∼ O(1). Such scenarios without the weakness of the WIMP couplings are within the reach of upcoming large volume detectors, such as XENONnT [45], and provide excellent experimental targets. We also find viable parameter space below the neutrino floor [46][47][48] which motivates developing new search techniques, see for example Refs. [49][50][51][52][53][54][55][56].
In Fig. 5, we show contours of required ξ to get the right DM abundance today in the light DM regime. We assume the DM freeze-out happens during a MD epoch  Fig. 4 is used for the constant ξ contours. We show the current Xenon1T [44] bounds (green), expected reach of low threshold superfluid helium searches [58] (red dot-dashed contour), the collider bounds on the B − L mediator from BaBar [57] (see also [43]) (blue), indirect DM annihilation limits from Planck observations [59] (purple), BBN bounds (gray), thermalization threshold (orange), and expected neutrino floor [48] (dashed magenta). The entropy injection opens new viable parameter space for such light DM candidates; see the text for further details. and we fix τ φ = 0.1 second, such that the entropy injection happens right before the BBN epoch. Note that in most of the viable parameter space ξ < η b ; in those parameter regions the BAU should be generated either after the entropy injection, e.g. using models from Refs. [31][32][33][34], or by φ decay itself.
Superimposed on this figure are the limits from searches for the B − L gauge boson in the sub-GeV mass range at BaBar [57] (see also Ref. [43]), as well as the limit on the s-wave annihilation cross section from the CMB [59], and current Xenon1T bounds [44]. At very small coupling values the system never enters thermal equilibrium with the SM and our treatment is not applicable. The lower m χ end of the plot is ruled out by BBN constraints.
Overall, while the no-dilution scenario is ruled out in this DM mass range, a large parameter space with ξ 10 −8 is widely open, motivating new search strategies that can cover the entire viable region. We find that part of the remaining parameter space is within the reach of upcoming low threshold direct detection techniques with superfluid helium [58,60] with one kg-year exposure; larger low threshold detectors can fully test the low mass end of this thermal DM scenario within the next decade [60].

V. CONCLUSION AND OUTLOOK
We studied the freeze-out abundance of thermal relic DM (denoted χ) in a next-to-minimal freeze-out scenario featuring an unstable relic φ that decays to the SM after its freeze-out. We analytically calculated the final χ abundance for general freeze-out interactions with p χ,φ 2 initial particles and q χ,φ 0 final particles and for arbitrary temperature dependence of the freeze-out cross sections.
We found that the final χ abundance is very sensitive to m χ,φ , their cross sections, and especially the initial number of particles in the interaction (p χ,φ ). The viable DM mass window moves as p χ,φ changes, even in the presence of entropy injection. However, even for the conventional WIMP scenario ((p χ , q χ ) = (2, 0)) the entropy injection means there is no particular mass value singled out -unlike the traditional WIMP models in which the weak scale emerges as the natural DM mass scale.
Using our analytic calculation we identified viable parameter space between the freeze-in and freeze-out cross section values (Fig. 2). We also pointed out that for the WIMP scenario, if the BAU is generated after (before) the entropy injection, the DM can be as heavy as ∼ 10 14 GeV (∼ 10 8 GeV) for m φ M Pl as shown in Fig. 3. We also used our results to argue that for freezeout interactions with p χ 3 there is viable parameter space with DM masses as high as m χ ∼ 10 12 GeV.
As a proof of principle, we studied the parameter space of a simple B − L gauge boson portal and showed that the entropy injection can open up parameter space for light, and heavy DM masses that were previously considered ruled out owing to direct detection bounds for the former and CMB bounds for the latter ( Figs. [4][5].
Our calculation and analysis can be extended in a few interesting ways. It is intriguing to consider the case of either the χ or φ relic density being set by a freezein instead of freeze-out. In particular, this way we can extend the range of Fig. 2 to lower cross sections in each direction. Furthermore, given the large viable masses that will not be accessible to conventional DM searches in foreseeable future, the best way to probe our setup is likely by looking for signatures of an early MD epoch. There are already proposals in the literature for searching for an early MD epoch [15][16][17]. Such an epoch can also affect the evolution of the Hubble parameter and thus the spectrum of gravitational waves from various sources in the early universe [61][62][63][64][65].
There are also intriguing studies of UV-complete scenarios which we leave for future explorations. In particular, now that we have a large parameter space for freeze-out via p χ → q χ interactions with p χ 4, it is interesting to think about natural particle physics models that could give rise to such freeze-out behavior.
Another interesting direction will be studying the correlation between DM and SM abundances in more detail. One of the best arguments for the existence of a portal between the DM and the SM is the closeness of their re-spective abundances. Since the SM abundance today is set by the BAU, a complete study should include a way to relate this asymmetry to the DM abundance today. Since the DM abundance in our setup is a function of the φ lifetime, a natural direction for exploration is UVcomplete models where the BAU is generated by decays of φ itself.
The relevant calculation for the entropy injection is carried out for a single heavy relic decaying to the SM in Ref. [23] section 5.3. (Similar calculation is carried out elsewhere in the literature as well; see for instance Refs. [20,66].) For completeness we repeat that calculation here; we also generalize that to the case where the dark sector has a different temperature than SM. While this is not the case in our study, this can be useful for future works.
We want to calculate the amount of entropy injected into SM when φ decays into SM bath. First, we should keep in mind that the energy density of φ is controlled by with ρ φ being φ's energy density, Γ φ their decay rate, and H the Hubble constant. The solution to this equation is where the subscript 0 here refers to some arbitrary origin of time (and not today) and t is the time passed since then. This equation and solution is valid for both RD and MD portion of the evolution, as long as φ is nonrelativistic.
As the φ particles decay to SM, they inject heat and thus entropy into it. The entropy injection can be calculated by: where in the second equality we used Eq. (A1). Notice that dS here is the infinitesimal change in the entropy of SM. The entropy of the SM can be written as where T and g * S refer to SM temperature and number of relativistic degrees of freedom (#dofs) for entropy. Without any entropy injection from decoupled particles, this total entropy of SM remains constant. Putting the last two equations together, we havė This equation can be solved as where now i, f are just two labels for the beginning and end of any interval we are studying. Combining this with Eq. (A2), we find where We can simplify this quantity as below From this point on, again for simplicity we assume i labels a point right after the universe becomes MD, i.e. ρ φ,i ≈ ρ R,i with ρ R denoting the radiation bath energy density. When the universe becomes MD, and before φ decays (t τ φ ), we have In the MD universe, we have R ∝ t 2/3(1+w) with w = 0, i.e. R ∝ t 2/3 . Hence, in the MD epoch, we have Combining these equations, we can rewrite Eq. (A10) as For t i τ φ , we can simply use t i = 0. However, we should be careful about the upper bound of the integral. As mentioned above, we are considering a MD universe. So, the equations are really only reliable as long as most of the φs have not decayed. Yet, let us for the moment assume the universe actually does remain MD even after the φ decays. This introduces some error since the R/R i scaling changes after φ decays, but with this assumption, we can take t f → ∞ and use the definitions of the Gamma functions: which is off by only around 10% from Eq. (5.72) of [23], where the integral is evaluated numerically with the right scaling of R in the RD epoch after our MD era. We can now use this in Eq. (A8) to find Note should be taken that g * and S i,f are referring to SM quantities. The entropy and the energy density of SM at t i can be written as Finally, again using ρ φ,i ≈ ρ R,i = π 2 /30 × g * ,i T 4 i , we find The only place where the temperature difference of the two sectors could enter the calculation is in #dofs g * (S) . Notice that all these #dofs factors entered from energy or entropy density of SM and did not care about #dofs in the dark sector. Thus, even if the two sectors had different temperatures, the calculation would have remained unchanged. Furthermore, since all particles are at the same T , g * = g * S , which is the final formula we use in Sec. III. The g * average is defined in Eq. (A10) and given its small power, it will always be an O(1) number.
Again we neglected the variation in #dofs in the integration above. (This can introduce some error for DM mass in the O(0.1)-O(10) GeV ballpark; for this mass range, the calculation should be carried out numerically for better precision.) Now we should replace H MD (m χ ) from Eq. (C4) and β from Eq. (B5) (for χ). The final abundance of χ becomes This is the equivalent of Eq. (B12) in a MD epoch freeze-out. Specializing to the simpler form of the cross section σv pχ−1 ≡ σ 0,pχ x −lχ This equation is in agreement with the results of Ref. [14]. To turn this yield into Ω χ h 2 for a massive relic, we use: By replacing the value of γ we find the final expression used in Sec. II C.
We can also use Eq. (C7) to find an expression for x f.o. during a MD epoch as well. In this case, we find (C14)

Appendix D: Additional Scenarios
In the body of the paper we reported the viable mass window and Γ φ values for a few (p χ,φ , q χ,φ , l χ,φ ) scenarios and assuming geometric cross sections. Our general formulas in Eqs. (12) and (16) can be used for other values of these quantities as well. We include a few more mass window plots in Figs. 6-9; comparing these figures to Fig. 3 can provide us with intuition on the effect of various parameters on the viable DM mass range.
We show the viable mass window for the case of p χ = q χ + 1 = 5 in Fig. 6; we find a large viable parameter space similar to the case of p χ = q χ + 1 = 4, see the bottom plot of Fig. 3.
In Fig. 8 we keep (p χ,φ , q χ,φ ) = (2, 0) and use the geometric cross sections for each freeze-out but vary the value of l χ,φ . We find that increasing l χ (l φ ) slightly shifts the viable parameter space to lower m χ (m φ ) values. However, the effect of l χ,φ on the viable mass range is clearly not as strong as that of p χ,φ .
These plots can be compared to those of Fig. 9 where we keep l χ,φ = 0 and instead change p χ,φ . We again observe a shift to lower m χ or m φ values, respectively, when p χ or p φ increases; the effect of changing p χ,φ is clearly stronger than changing l χ,φ in Fig. 8.
χ��→�� ϕ��→�� � ϕ →�� σ→σ ���� FIG. 8: Similar to Fig. 3 but with different lχ and l φ values. We find that increasing each of these quantities slightly shifts the parameter space to lower values of mχ or m φ , respectively. We note that the effect of changing l χ,φ is much smaller than changing p χ,φ , via comparison to the top plot of Fig. 3.