The domain of thermal dark matter candidates

We consider, in general terms, the possible parameter space of thermal dark matter candidates. We assume that the dark matter particle is fundamental and was in thermal equilibrium in a hidden sector with a temperature $T'$, which may differ from that of the Standard Model temperature, $T$. The candidates lie in a region in the $T'/T$ vs. $m_{\rm dm}$ plane, which is bounded by both model-independent theoretical considerations and observational constraints. The former consists of limits from dark matter candidates that decoupled when relativistic (the relativistic floor) and from those that decoupled when non-relativistic with the largest annihilation cross section allowed by unitarity (the unitarity wall), while the latter concerns big bang nucleosynthesis ($N_{\rm eff}$ ceiling) and free streaming. We present three simplified dark matter scenarios, demonstrating concretely how each fits into the domain.


I. INTRODUCTION
Of all the scenarios for dark matter (DM), one of the most appealing possibilities is that it is made of elementary particles that were in thermal equilibrium in the early universe. If this is the case, then just like the photons of the cosmic microwave background, their abundance today reflects the one they had when they decoupled, or 'froze out', from their thermal environment [1][2][3]. We refer to DM candidates that were once in thermodynamic equilibrium as thermal candidates. One class of thermal candidates consists of particles that decoupled when they were non-relativistic. Because of the exponential dependence of their Maxwell-Boltzmann distribution, the moment of freeze-out is critical: the later the decoupling, the smaller their relic abundance and vice versa. Alternatively, the DM particles may have decoupled when they were relativistic. In that case, their relic abundance depends on the moment of decoupling only through entropy transfer effects (like for the Standard Model (SM) relic neutrinos). The primary question we would like to address in this work is the following: what is the domain of all possible thermal candidates?
The answer is well known if the DM particles were in thermal equilibrium with SM particles. As we will review in section II, thermal candidates which give the correct relic abundance lie somewhere between a few eV (a value which is ruled out by structure formation constraints, as is well known) and about 100 TeV. The lower bound corresponds to DM that decoupled when relativistic [1], and the upper bound to DM that decoupled when non-relativistic with the largest possible annihilation cross section compatible with unitarity [4]. But what if the DM is secluded and lives in a hidden sector (HS) with a temperature, T , that differs from that of the thermal bath of SM particles (the visible sector, or VS), which has temperature T ? The possibility that DM resides in a HS (see e.g. [5][6][7]), possibly with its own temperature, perhaps due to some reheating story after inflation (see e.g. [8,9]), has by now been studied in many scenarios. However the full domain of possibilities has not been systematically discussed in the literature. In this work, we aim to determine it in general. This domain of all DM thermal candidates is two dimensional, as it depends both on the DM mass and on the temperature ratio, ξ = T /T . In particular, we will derive the equivalent of the lower and upper mass bounds on thermal DM candidates as a function of T /T . We will also determine which part of this domain is excluded by various observational constraints.
The plan of this paper is as follows. To lay the groundwork, we will begin in section II by briefly reviewing the standard result for particles that were in thermal equilibrium with the SM. In section III, we will extend this to DM particles which reside in a HS. In section III A, we will discuss the theoretical constraints, in particular from unitarity. Next, in section III B, we will discuss the bounds from limits on the number of relativistic degrees of freedom at the time of big bang nucleosynthesis (BBN) on one hand, and from free-streaming (FS) on the other. The main result of our analysis is summarised in Fig. 2, which depicts the domain (white region) of all possible thermal candidates which survive these constraints in the plane T /T vs. m dm . The temperature ratio here is defined at the moment of decoupling of the DM particles. Three simple but illustrative explicit DM models are considered in section IV. We draw our conclusions in section V. This work is built upon and extends results from [10,11].

II. DOMAIN FOR T = T
We start by reviewing a classic result regarding possible thermal candidates that were in thermal equilibrium with the SM thermal bath. If this were the case, the lightest possible thermal candidate would have a mass at the eV scale. This lower bound is obtained assuming the largest possible number of DM particles when it arXiv:2105.01263v1 [hep-ph] 4 May 2021 decouples, that is to say when DM decouples relativistically, in which case Ω dm h 2 = m dm n dm,0 ρ c,0 /h 2 ≡ m dm Y dm,dec s 0 ρ c,0 /h 2 ≈ 0.12 g m dm 6 eV g * s,0 g * s,dec .
The parameter g = c · g dm counts the effective degeneracy of the DM candidate, where c is 1 (3/4) for a boson (fermion) and g dm is the internal DM degrees of freedom, while g * s,0 counts the total number of degrees of freedom contributing to the entropy today, g * s,dec is the same but at DM decoupling and ρ c,0 is the critical energy density today. From this, the DM mass is required to be at least m dm m dm,CM ≈ 6 g g * s,dec g * s,0 eV .
In the context of SM neutrinos, Eq. (2) is the Cowsik-McClelland (CM) bound [1]. 1 It is usually considered as an upper limit on the sum of neutrino masses. In the present context, it is a lower limit since it implies that the DM abundance is too small when m dm < m dm,CM .
A thermal candidate with a larger mass than (2) must decouple when non-relativistic so that its abundance is Boltzmann suppressed. This is the WIMP (weakly interacting massive particle) scenario, in which the abundance is inversely proportional to the annihilation cross section (thermally averaged at the time of freeze-out), which is required to be about σv ≈ 3 · 10 −26 cm 3 /s to match observations of the DM abundance [13]. Although this number depends on the DM mass only logarithmically, the cross-section itself depends on it, albeit in a model-dependent way. In the non-relativistic limit, relevant for WIMP freeze-out, and when the DM mass is the largest relevant mass scale, the annihilation cross section scales like σv ∝ 1/m 2 dm . For fixed m dm , unitarity sets the upper limit on the DM mass by fixing the maximal possible annihilation cross section of a pair of particles with angular momentum J, where v is their relative velocity. This translates into an upper limit on the mass of thermal candidates, known 1 For SM neutrinos, g → g ν = 2 × 3/4 per species and, due to entropy transfer from the electrons and positrons to the photons after SM neutrino decoupling, g * s,0 /g * s,dec = 4/11. For Ωh 2 1, this gives m ν 92 eV. Here, in general, we have in mind that g is of the order of a few, corresponding to O(1) species of DM. See [12] for a recent scenario in which g is considered to be a huge number. as the Griest-Kamionkowski (GK) bound [4], Together, the lower (2) and upper (4) limits define a one-dimensional domain, Below this range, the DM is always under-abundant, while above it is over-abundant. The abundance of the DM increases with the DM mass at both the lower and upper limits in (5) (with Ω dm,CM ∝ m dm and Ω dm,GK ∝ m 2 dm respectively). Thus, assuming continuity of the abundance as a function of the DM mass, Ω dm (m dm ) implies that for any given DM model, one obtains the observed relic abundance for at least one value of m dm within the range of Eq. (5). Where precisely the thermal candidates lie is model-dependent, so the scope of this continuity argument is limited (the DM mass range (5) spans about 13 orders of magnitude), but is robust. Another generic result is that for a given model there is in general an odd number of DM candidates, except for specific and thus fine-tuned choices of parameters (this will be seen more explicitly when we consider specific models in section IV).
These well known facts are illustrated by the case of a massive Dirac neutrino [1,2,14]. Consider Fig. 1, reproduced from [15]. It shows Ω dm h 2 as a function of m dm , i.e. the neutrino mass, the only free parameter in this setup. As the figure shows, there are three possible thermal candidates: one at low mass, corresponding to the CM mass bound; a heavy one, corresponding to the GK mass bound 2 ; and an intermediate one, with a mass in the GeV range. The latter is the Lee-Weinberg (LW) bound [2]. There are thus three possible candidates in this simple model (all of which are excluded by well known constraints).
We want to argue that all models with thermal DM candidates share similar features. What are possible caveats? The most important hypothesis is that DM was in thermal equilibrium in the early universe. We will come back to this condition and its implications in section III A 2. Here, we assume DM in thermal equilibrium and explore the other possible shortcomings, beginning with the GK mass bound. The GK mass bound (4) is quite general but still rests on some assumptions. First, it assumes that the DM particle annihilates with itself or with its antiparticle. Instead, it is possible that the abundance is set by co-annihilation of the DM particle with some almost degenerate companions [16,17]. 2 The behaviour of the relic abundance as a function of the neutrino mass for masses above the Z resonance is peculiar. The fact that it decreases is due to emission of longitudinal gauge bosons with a cross-section that scales as σv ∝ m Provided there is sufficiently fast chemical equilibrium between these companion particles and the DM candidate, the DM abundance is determined by the most efficient annihilation channel(s). The GK mass bound thus applies to the companion particles, whose mass cannot be too much larger than the DM particle itself. Hence m dm m dm,GK in that case too. Another possibility is that DM is complex and that, like for the baryons in our universe, the DM abundance is set by an asymmetry instead of thermal freeze-out, see e.g. [18,19]. This scenario is also subject to a GK bound. Indeed, a key feature of the asymmetric DM scenario is that the thermal, or symmetric, component of the DM abundance must be efficiently depleted, otherwise the asymmetry is hardly relevant. Efficient annihilation of the DM particle-antiparticle pairs (or their companions, as in the co-annihilation scenario) requires a larger cross section than the canonical value 3 · 10 −26 cm 2 /s. This is suggested by the dashed line depicted in Fig. 1. In the absence of an asymmetry, the thermal abundance would be far below observations. Translating this to the maximum possible cross section set by unitarity considerations, the mass of the asymmetric DM must be below the standard GK bound (4), m dm,asym < m dm,GK .
Finally, the GK mass bound only applies if DM particles are fundamental. Indeed, composite DM candidates [4] with radius much larger than their Compton wavelength, r dm 1/m dm , can have a geometric cross section, σ ∼ πr 2 dm (see, for instance, [20] for a recent concrete example). With this possible exception kept in mind, we assume that generic thermal particle DM candidates all fall in the domain (5).

III. DOMAIN FOR T = T
We now extend the established results of section II to scenarios in which the DM was in thermal equilibrium in a hidden sector (HS). This HS is assumed to be feebly interacting with the VS, so that the different sectors may have distinct temperatures, T = T , at the moment of DM decoupling. Then, for a given T /T = 1, we expect a generalised form of the interval (5), with Understanding the dependence of this range on T /T is one of our main goals. The result is summarised in Fig. 2 (white region) in the plane T /T vs. m dm . The boundaries of that regions are explained in sections III A and III B. Illustrations in terms of explicit models are given in section IV. From now on, we use ξ ≡ T /T for the ratio of temperatures of the hidden and visible sectors.

A. Theoretical bounds
In this section, we consider theoretical constraints on the domain of all possible thermal DM candidates. These will be set by 1) the DM abundance from relativistic decoupling (relativistic floor), 2) the condition of thermalisation of the DM within the HS (thermalisation condition), and 3) the DM abundance set from decoupling when non-relativistic with largest possible cross section (unitarity wall). We study the corresponding limits in that order.

Relativistic floor
We first discuss the generalisation of the Cowsik-McClelland bound for ξ = 1. The basic assumption is that the DM particle decoupled when relativistic at a temperature T dec . For given values of m dm and ξ, relativistic decoupling generates the largest possible relic DM abundance. Comparing with the observed DM density gives a relation between ξ (defined at decoupling) and the DM mass, ξ(m dm ). Following [10,11], we call this relation the relativistic floor as for fixed m dm , ξ < ξ(m dm ) gives too little DM. As in the previous section, the degeneracy of the DM species will be denoted by g , including a factor 3/4 for fermionic DM. The effective number of relativistic degrees of freedom in entropy in the HS is denoted with a prime, g * s .
One can distinguish two main scenarios of relativistic decoupling. The first scenario is analogous to the neutrino relativistic decoupling scenario in the SM. In the same way that the interaction rate for SM neu- , we can imagine that the DM interaction rate was ∼ T 5 /Λ 4 below some scale T ∼ Λ, with Λ> T dec m dm . Typically, Λ is the mass of some heavy mediator. The second scenario instead invokes a heavy "cut-off" scale which is not the mass of a heavy mediator but the mass of the particles the DM scatters into. If these particles are heavier than the DM, the annihilation rate became Boltzmann suppressed when T went below their mass. Schematically, this is analogous to non-relativistic DM decoupling in the sense that decoupling was due to a Boltzmann suppression. However, this happens at a higher temperature than m dm , when DM is still relativistic. A difference between these two scenarios is that in the second one, DM is reheated just before its decoupling, as the heavier particle becomes non-relativistic and annihilates into it (similarly to photons which are reheated during the e + e − annihilation catastrophe). Thus, in this case, the initial value of T /T , which we call ξ in (before the heavier particle became non-relativistic), is different from the value of T /T at DM decoupling, ξ dec . They differ by where g * s,in > g * s,dec counts the HS degrees of freedom at the corresponding times. 3 In the first scenario in- 3 As an example, consider a HS consisting of massive dark photons, electrons and positrons. The dark photon is the lightest particle of the HS, m γ < m e , and so is a DM candidate. In that case, g * s,dec counts the dark photon's degrees of freedom while g * s,in also includes the dark electrons and positrons. Also, the initial dark photon temperature corresponds to T in ≈ m e so that the abundance of the dark photons after decoupling is Y γ = (g * s,in /g * s,dec )Y dm (T = m e ). More stead there is no such reheating of the HS and ξ in = ξ dec as (g * s,in /g * s,dec ) = 1. Taking into account this possible effect, we find the general formula where g * s,0 and g * s,dec count the total of degrees of freedom (i.e. from both sectors) at the corresponding times. Inversion of Eq. (8) gives the ξ in required to produce the correct DM abundance, This relation defines the relativistic floor, giving the blue exclusion region in Fig. 2. The small kinks along the floor are due to the evolution of the VS contribution to g * s .

Thermalisation condition
By definition, a thermal DM candidate was in thermodynamic equilibrium. We must make sure that this condition was satisfied. In the following, we refer to thermal equilibrium as both kinetic and chemical equilibrium, so that the DM initial abundance is entirely determined by the HS temperature, T . 4 Consider again the case of SM neutrinos in the early universe. At high temperatures, larger than the electroweak scale, their interaction rate was Γ ∼ α 2 W T while the expansion rate of the universe was H ∼ √ g * T 2 /m pl [3]. Thus, SM neutrinos were in thermal equilibrium at temperatures T ≤ T eq α 2 W m pl / √ g * ∼ 10 14 GeV.
Translating this to relativistic particles in a HS interacting with a rate Γ ∼ α 2 T , gives where ξ eq = (T /T ) eq and α is a HS analogue of α W . So, if ξ < 1 (ξ > 1), particles in a HS entered thermal equilibrium at a later (resp. earlier) time in the history of the universe than those in the VS. This is illustrated in Fig. 3, where we plot Γ/H as a function of generally, the dark photons could have inherited the entropy from all the charged particles that were once in equilibrium in the HS. m dm , Γ ∼ T in all cases. As the temperature of the HS drops, we consider three scenarios for DM decoupling, see text: 1) solid lines, the DM particles become non-relativistic, 2) dot-dashed lines, the DM annihilates into heavier particles of mass m (here chosen to be m /m dm = 3), and 3) dashed lines, the DM interacts through some heavy mediator and decouples while still relativistic. We show these three cases for the same DM mass but for two choices of ξ = T /T (orange and grey lines, respectively). For the orange curves, the intersections with the horizontal line correspond to the temperatures T eq and T dec . Between these two temperatures, the DM is in thermal equilibrium, Γ > H. The grey curves correspond to smaller ξ. In particular, the solid line is such that a DM particle barely reaches thermal equilibrium around T ∼ m dm ∼ T eq = T dec .
x ≡ m dm /T , with thermal equilibrium corresponding to Γ/H > 1.
The way DM decoupled is model-dependent. Here, we consider three generic possibilities. The first is that the DM stayed in thermal equilibrium until it became non-relativistic and its abundance changed from n dm ∝ T 3 to being Maxwell-Boltzmann suppressed at T m dm , with n dm ∝ exp(−m dm /T ). This case is represented by solid lines in Fig. 3. The second and third possibilities are the two relativistic decoupling scenarios already discussed above, of annihilation cut off by a heavy mediator mass or by the mass of the particles that the DM scatters into. The corresponding rates are shown in Fig. 3 as dashed lines and dot-dashed lines respectively. Their behaviour differs from the first case when T becomes smaller than the corresponding cut-off mass scale. Other scenarios are possible, such as combinations of the above three cases. However, in what follows, and in particular in considering explicit models in section IV, we will focus on the three simple possibilities considered here.
These three cases are depicted in Fig. 3 assuming the same DM mass, but for two choices of ξ (orange and grey lines). Candidates along the orange curves were in thermal equilibrium between T eq and T dec . If the ratio ξ increases (decreases), these curves move up (resp. down). Thus, a generic feature which can be read from Fig. 3 is that the temperature range within which the DM was in chemical equilibrium shrinks as ξ decreases. If ξ is too small, all other things being kept the same, then the particle was never in equilibrium. The absolute minimum value of ξ for which there is thermalisation for a given DM model is therefore obtained by setting T eq ≈ T dec (≈ m dm ), see the solid grey curve in Fig. 3. This corresponds to the case of freeze-out of a mildly non-relativistic DM particle, thus it lies close to the relativistic floor.
To be specific, let us take Γ ∼ α 2 T as above when T m dm . This leads to the T temperature range for thermal equilibrium 5 or equivalently, This range shrinks to zero as ξ decreases down to a minimum value found by taking T dec ≈ T eq , This condition of thermalisation of the HS crosses the relativistic floor, Eq. (9), at corresponding to ξ min ∼ 10 −6 α −2/5 .
As expected, ξ can be lower if the DM and its companion particles have stronger interactions. However, for fundamental particles, cross sections are constrained 5 This range assumes that the DM mass is the only relevant mass scale (the first of the three possibilities discussed above). If, instead, the DM is interacting with heavier particles of, say, mass m , then the condition of thermalisation is more constraining as decoupling occurs for T ∼ m . Thus m dm → m in Eqs. (11,12,13) and the thermalisation condition becomes independent of the DM mass for m dm m , We will meet such a situation in section IV.

6
from above by unitarity [4]. Thus, even if the DM abundance does not depend on the cross section in the case of relativistic decoupling, the DM mass along the relativistic floor cannot be arbitrarily large due to the requirement of thermalisation.
To determine the absolute limits on the DM mass along the relativistic floor, we consider the maximal, thermally averaged cross section allowed by unitarity. Assuming annihilation in J = 0 state, it is given by [11] σv = π where I is , where = 1 for fermions, −1 for bosons and 0 in the classical (Maxwell-Boltzmann) case, and k −,max = 1 − 4x 2 /w k 2 + − w. As we will consider regimes in which the DM is relativistic or mildly non-relativistic, it is important to keep track of quantum statistics effects. In the relativistic limit, with n +1 = 5π/12 and n −1 = 15π/16. For a classical Maxwell-Boltzmann distribution, (16) takes the form σv = 4π This gives n 0 = π/2 in the relativistic limit and 1/v = Using Eqs. (16)- (17), we can find the lowest temperature T dec = T eq at which a HS can be in thermal equilibrium, given the maximally allowed cross section of the DM. This condition is depicted as the orange region in Fig. 2. It is set by simply requiring that n dm σv ≈ H at T dec = T eq , which gives what we call the thermalisation condition, where g * ,dec counts the effective number of degrees of freedom contributing to the expansion rate at T dec = T eq . Modulo the dependence through x dec = m dm /T dec , we see that essentially ξ eq ∝ m 1/2 dm as in (13). As above, the thermalisation condition crosses the relativistic floor, Eq. (8), at which gives the lowest possible temperature for a HS with thermal DM that decouples when (barely) relativistic. This temperature corresponds a DM candidate of mass

Unitarity wall
If we depart from the point corresponding to (22), going along the orange thermalisation condition, the HS temperature increases as T ∝ m dm , so that the possible DM candidates are less and less relativistic and we enter into a secluded regime of non-relativistic freezeout [6,7]. In this case, the relic abundance depends on the annihilation rate. In the instantaneous freezeout approximation, i.e. fixing x dec = m dm /T dec from the condition Γ/H| T =T dec = 1 and determining Y dm by assuming that the yield after freeze-out is equal to n dm,eq /s| T =T dec , the relic abundance obtained is 6 Ω dm h 2 = 4.7 · 10 8 g dm g 1/2 * ,dec x dec ξ g * s,dec m pl σv GeV (23) with x dec given by 7 x dec ln 0.038 ξ 2 σv m pl m dm g dm √ g * ,dec We first consider the case in which the expansion of the universe is dominated by the VS. This is natural for 6 It can be checked that the relic abundance obtained in this way differs from the one obtained from a proper integration of the Boltzmann equation by less than a factor ∼ 2, regardless of the value of ξ considered in the allowed domain of Fig. 2, and as long as x dec 3. The agreement can be somewhat further improved using the following expression: The power of ξ in the logarithm which gives x dec can be obtained different from the value 2 in Eq. (24), if one uses another prescription for determining x dec , see [7,21] for other prescriptions which lead to a value 5/2 or 3/2 in the ξ dec < 1 case. ξ 1. In this case Taking the maximal annihilation cross section allowed by unitarity, Eqs. (16) and (17), gives an upper bound on ξ dec as a function of m dm , This equation leads to the diagonal part of the red boundary in Fig. 2. For fixed ξ dec , a value of the DM mass beyond this diagonal leads to an excess of DM. Clearly, the relation (27) generalises the standard unitarity bound which is set assuming T = T [4].
Let us make a few remarks. First we note that the condition (27) scales as ξ dec ∝ 1/m 2 dm , so heavier DM candidates are possible for T < T . Second, this expression, which is obtained for the unitarity limited cross section, implies that in general ξ dec ∝ σv . Indeed, for fixed DM mass, decreasing the cross-section by decreasing the coupling leads to a increase of the DM abundance (as for standard non-relativistic FO at T = T ), which is compensated by a decrease of the HS temperature. Lastly, we see from Eq. (24) that x dec decreases as the DM mass increases when σv ∝ m −2 dm . While x dec = O(20) for T ∼ T , the DM is less non-relativistic at decoupling if T T , as expected since in this case the DM particles are less numerous already to start with and thus must be less Boltzmann suppressed to account for the relic density.
The crossing of the unitarity constraint (27) and the condition for thermalisation (20) gives the heaviest possible thermal DM candidate (denoted by the little red dot in Fig. 2), which has mass m dm,max ≈ 52 PeV (28) when DM is a Dirac fermion. This is slightly heavier than (22) and corresponds to a temperature ratio, ξ ≈ 6.9 · 10 −5 (29) to be compared with (21). It may be worth mentioning that taking freeze-out at T dec ∼ m dm,max implies that the DM decoupled around T ∼ 10 11 GeV.
So far, we have assumed that the expansion of the universe was dominated by the VS. If there are many degrees of freedom in the HS or if ξ 1, the entropy and energy densities of the HS can be dominant at the time of the decoupling of the DM. In this case, we can approximate the expansion rate and entropy by by neglecting the VS contributions, where g * is the effective number of HS relativistic degrees of freedom. In this approximation, the DM abundance no longer depends on ξ. Put simply, the VS, and its temperature, are irrelevant at the time of DM decoupling. Thus, because the annihilation cross section is set by unitarity, there is a unique value for the maximum DM mass as long as the HS dominates the universe.
In this scenario it is crucial that after DM freezeout, most of the entropy and energy of the HS, which is shared by the DM companions, is transferred to the VS. There are two main possibilities, depending on whether or not the entropy is conserved in this transfer: -If the transfer occurs while the companions are in thermal equilibrium with the VS, entropy is conserved and Eq. (23) together with Eqs. (30) and (31) apply. This gives the upper bound on m dm corresponding to the vertical part of the red exclusion region in Fig. 2, To draw this limit, we considered a minimal scenario, with O(1) companion particles on top of the DM and only the known, SM degrees of freedom in the VS. It intersects the curve (27) at roughly ξ ≈ 4 g * /g * ≈ 3. If there were more particles in the VS than the O(100) of the SM, the line would shift towards the left but, as it scales as the fourth root of the number of particles, their effect would be mild.
-If transfer happens through a slow, out-of-equilibrium decay of the companion particles, entropy is produced. Energy conservation leads to heating of the VS [22] and dilution of the DM abundance. Such a scenario has been considered in [23]. The effect of entropy production is to slightly shift the vertical line in Fig. 2 to the right. The shift can be approximated by multiplying (32) by a factor of obtained by imposing that ρ tot = ρ V S + ρ HS is conserved and then extracting the evolution of the temperature ratio due to this conservation. As (32) is derived assuming entropy conservation, the factor (33) divides by the contribution coming from entropy conservation (g * s , as defined in Eq. (8)) and multiplies the contribution from energy conservation (g * ). This factor lies between f ∼ 1 and f ∼ 1.6 depending on the VS temperature when the entropy is injected (T < 10 keV and T > TeV respectively) for g * = O(1). We note at this stage also that the DM abundance from relativistic freeze-out given by Eq. (8) also assumed entropy conservation. An out-of-equilibrium processes resulting in entropy production would cause the RHS of Eq. (8) to be a factor of (g * s,dec /g * s,in ) 1/4 smaller. Since g * s,dec ≤ g * s,in , this could slightly raise the relativistic floor.

B. Observational constraints
In this section, we consider two observations that set important restrictions on the domain of thermal DM candidates. These are depicted in Fig. 2, see the green (N eff ceiling) and purple (free streaming) regions. Our aim, as in the previous section, is to be as modelindependent as possible.

N eff ceiling
We know that the universe was dominated by radiation at the time of big bang nucleosynthesis (BBN) and until about the recombination epoch.. A crude constraint is set by considering the possible contribution of the HS particles to the number of relativistic degrees of freedom at the time of BBN and recombination. This is expressed in terms of N eff , the effective number of neutrinos. Too large a value of ∆N eff at temperatures around T ∼ O(1) MeV will increase the Hubble rate and thereby impact the abundances of light nuclei, which are rather well measured. The latest CMB measurement by Planck [13] gives N eff = 2.99 ± 0.17 (the results from BBN are similar [24]), which is to be compared with the SM prediction, N eff = 3.04, see e.g. [25]. Hence, we impose the constraint that ∆N eff ≤ 0.29 at 2σ. We will distinguish two contributions to ∆N eff , one from the DM degrees of freedom and one from other HS particles. The latter is much more model-dependent than the former.
We first look at the contribution of the DM degrees of freedom. The shift in N eff due to the dark matter at T = 1 MeV can be computed by taking the ratio of the DM and neutrino energy densities. Using the fermionic (+) or bosonic (-) equilibrium number densities, we can express ∆N eff as where z = E dm /T with T BBN = ξ T BBN Any other relativistic degrees of freedom in the HS will give an additional contribution to ∆N eff , so we are being very conservative here. The constraint it leads to is depicted by the green region of Fig. 2. Its two key features may easily be understood. Firstly, it is clear from Eq.
(34) that the contribution of a DM particle to ∆N eff is only sizeable when it is relativistic at the time of BBN, T BBN ∼ 1 MeV, i.e. for ξ m dm /MeV. This basic condition sets the diagonal part of the green boundary for m dm 10 MeV. Next, suppose that DM is indeed relativistic at BBN, i.e. m dm T BBN . For fermionic DM, we have ∆N eff g dm (11/4) 4/3 ξ 4 /2, meaning that a Dirac fermion (g dm = 4) counts as two families of neutrinos, while a Majorana DM particle (g dm = 2) counts as one. For bosonic DM instead, ∆N eff 4g dm (11/4) 4/3 ξ 4 /7 (with g dm = 1 for a real scalar). Collectively, we write this as ∆N eff = 4g dm,eff (11/4) 4/3 ξ 4 /7 with g dm,eff = B g dm,B + 7/8 F g dm,F . The bound from Planck then corresponds to which excludes any value of ξ larger than ∼ 1 when m dm T BBN . As soon as T is below T , however, the ξ 4 suppression allows a large number of relativistic degrees of freedom [6]. Notably, the constraint (35) is independent of the DM mass. This limit corresponds to the horizontal part of the green region. To draw it, we assumed g dm,eff = 3.5 (i.e. a Dirac fermion), but taking instead g dm,eff = 10, for instance, barely changes the figure.
Next, we discuss the possible additional contribution of companions in the HS. This is much more modeldependent as it is related to the number of particles which are left after DM freeze-out and how their number changes afterwards. Consider first the typical situation in which DM annihilates into lighter particles which decouple relativistically. In this case, if they disappear (say, by decaying to SM particles) before the BBN epoch they are harmless for BBN. However, since these particles will not disappear until they are nonrelativistic, this requires their mass to be larger than T BBN . In this case, if the mass of the companion particle is about the DM mass the constraint remains as given by the green area in Fig. 2. If the companion particle is lighter than the DM by a fixed ratio, the diagonal part instead moves to the right by a factor of m dm /m comp . Alternatively, for a fixed value of m comp , one gets an horizontal line, i.e. a fixed upper bound on the temperature ratio: ξ < m comp /T BBN . In the extreme case, where the companion is much lighter than the DM particle and/or it decays after BBN, the bound (35) becomes ξ 0.60/(g dm,light ) 1/4 where g dm,light now includes all HS degrees of freedom that are still abundant at the BBN time (see also [6]). This corresponds to a horizontal exclusion line which for low DM mass is somewhat below the horizontal boundary of the green region in Fig. 2 Fig. 2. If this is not the case, values of ξ larger than ∼ 1 are excluded. The above discussion was deliberately general. For specific models, the constraints can be expected to be stronger than the generic BBN bound we discussed here. The most delicate situation is when the DM and/or its companion particles annihilate or decay at the time of BBN, as the production of light elements can be affected by changes to the expansion rate and, more importantly, by energy transfer into the VS, see e.g. [26][27][28]. Determining this, however, necessarily demands specifying a HS scenario.

Free-streaming
Thermal DM cannot be too light, otherwise it remains relativistic for too long and does not permit the formation of large-scale structures that we observe. For a thermal relic with ξ = 1, the strongest bound is m dm > 5.3 keV [29], obtained from Lyman-α forest data. This bound can be generalised to account for different values of T /T by converting it to a limit on the DM free-streaming horizon, the average distance a DM particle travels after production.
The average momentum of a population of particles in thermal equilibrium with temperature T is This average persists after the population goes out of thermal equilibrium. The particle species can then be said to be non-relativistic when p ≤ m, i.e. below T nr = m dm /3.15 for a fermion, with the corresponding time being t nr = 1/[2H(T nr /ξ)]. If t nr < t eq = 1.9 × 10 11 s, the time of matter-radiation equality, then the free-streaming horizon can be estimated as [3] t eq t nr a eq 5 + log t eq t nr g * s,0 + ξ 3 g * s,0 g * s,dec + ξ 3 g * s,dec , where a eq = 8.3 × 10 −5 is the scale factor at t eq . On the other hand, if t nr > t eq , we have λ FS 6(t 2 eq t nr ) 1/3 − t eq a eq g * s,0 + ξ 3 g * s,0 g * s,dec + ξ 3 g * s,dec since for t > t nr , v(t) a nr /a (t nr /t) 2/3 . Note that Eqs. (38) and (39) coincide when t nr = t eq .
An early decoupled fermionic thermal relic of mass 5.3 keV has a free-streaming horizon λ F S 0.066 Mpc. Imposing this upper limit on the free-streaming horizon leads to the bound given by the purple region in Fig. 2. The behaviour is quite different depending on whether ξ is smaller or larger than 1. When T T , then t nr ∝ ξ 2 /m 2 dm , hence λ FS ∝ ξ/m dm , up to the logarithm in Eq. (38) (see also [11]) . Conversely, when T T , the hidden sector dominates, with H(T nr ) ∼ m 2 dm /m pl . Then if t nr < t eq , as is the case for m dm 10 eV, we have λ FS ∝ 1/(ξm dm ), neglecting the logarithm. As a result of all that, the absolute minimum value of m dm lies at the intersection of the free-streaming constraint and the relativistic floor, which is at m dm 1.0 keV and ξ 0.16. Thus we obtain see also [11].

C. Domain of thermal DM candidates
We summarise here the constraints discussed in sections III A (theory) and III B (observations). Fig. 2 depicts the bounds obtained from the ξ = T /T relativistic floor (Eq. 9, boundary of blue region), the thermalisation constraint (imposing unitarity, Eq. 20, orange) and the unitarity wall (Eqs. 27 and 32, red). Along the ξ = 1 horizontal dashed black line lies the domain of DM candidates that were in equilibrium with the SM bath when they decoupled, Eq. (5). The observational constraints come from the N eff ceiling (green region, Eq. 34, assuming the products of DM annihilation are gone by the time of BBN) and from free-streaming (purple region, Eqs. 38 and 39). All these constraints define the possible domain of thermal DM, as given by the white region. In specific cases, this region can be further reduced. For instance, the dashed green line gives the approximate (somewhat model-dependent) upper bound on ξ which holds when the (or some of the) companion(s) of the DM are still around during BBN. The vertical part of the red region corresponds to Eq. (32) rather than Eq. (33) (i.e. for ξ > 1, conservation of entropy rather than of energy when the HS particles disappear into SM particles later on) but the difference is small.
In the blue region, bounded by the relativistic floor, the thermal DM particles would be under-abundant, while in the red region, bounded by the unitarity wall, they would be over-abundant. Only in between can viable thermal DM candidates exist. To make this clear, consider Fig. 4, in which we show the DM parameter density against the DM mass for different choices of ξ = T /T . The horizontal dashed black line corresponds to the observed density, Ω dm h 2 ≈ 0.12. The solid curves below this line correspond to particles that decoupled when relativistic for different choices of ξ (from ξ = 1 to ξ = 10 −5 ). Thus they give the maximum value of Ω dm one can obtain for the values of m dm and ξ considered. The intersections of these curves with the horizontal line define the relativistic floor. The solid curves above the dashed line correspond to particles that decoupled when non-relativistic, assuming the unitarity bound for their annihilation cross section (thus giving the minimum value of Ω dm one can obtain for the values of m dm and ξ considered). The intersections of these curves with the dashed line define the unitarity wall. For fixed ξ, the line of given ξ that overlaps the horizontal dashed line defines the corresponding thermal DM mass range within which it may be possible to have the observed relic density; this is illustrated by the horizontal black solid line in the case of ξ = 1. As explained above, and as can also be seen on this figure, as ξ decreases the DM mass range shrinks and also shifts toward higher DM masses. Around ξ ∼ 10 −5 , the non-relativistic floor and unitarity wall merge (see blue curves), corresponding to DM candidates around the PeV scale. All together, they form the domain of thermal DM candidates, as depicted in the plane ξ vs. m dm in Fig. 2.
How specific models fit into this picture is the subject of the next section but, as an illustration, in Fig. 4 we reproduce the case of a Dirac neutrino for ξ = 1, see the grey line and Fig. 1. As discussed in section II, this theory has three thermal DM candidates. The one that lies between the floor and the wall (the Lee-Weinberg candidate around m dm ∼ GeV) will move around if, forgetting for the sake of the argument the relation with electroweak interactions, we change its interactions, keeping ξ = 1. For instance, increasing α W will make it moves toward lighter DM while increasing the mass of the gauge boson(s) moves it toward heavier  DM, etc. One may further convince oneself that, playing with the parameters of this model, the horizontal black solid interval may be filled by thermal candidates.

IV. EXPLICIT MODELS
In this section, we illustrate how some concrete models of hidden DM fit in the domain of thermal DM candidates of Fig. 2. We will do so using three simplified DM scenarios. In the first one, fermionic DM particles annihilate into a pair of vector bosons. We dub this the t-channel scenario after the topology of the treelevel annihilation process. In the second scenario, the DM annihilates into HS fermions through an s-channel process. Finally, we consider scalar DM annihilation into HS scalars via a contact interaction. We will skim over the fate of the DM companions, except in the last scenario.

A. Scenario 1 : t-channel
For this first model, the DM consists of a Dirac fermion, χ, charged under a local U (1) . The dark photon has mass m γ ; the origin of its mass is not important here, but m γ can be smaller or larger than the mass of the χ, called m dm as above. We consider separately the cases ξ < 1 and ξ > 1.
1) If ξ < 1, the DM and the dark photons have little impact on the expansion rate. For fixed DM abundance, HS coupling α and dark photon mass, the temperature ratio ξ can be expressed as a function of the DM mass, ξ = ξ(m dm ). An example of such a relation is the black solid curve in Fig. 5. Its salient features are the following: -For m dm < m γ , the DM particles can efficiently annihilate into the γ only as long as T m γ . The annihilation becomes Boltzmann suppressed for T m γ , so the DM particles decouple relativistically around T ∼ m γ > m dm . Thus, the only way to account for the relic density is if DM lies on the relativistic floor, see Fig. 5. This can further be seen in Fig. 6 (a modelspecific version of Fig. 4), which displays the DM abundance as a function of the DM mass for several different values of ξ. For m dm < m γ , the DM density parameter is ∝ m dm , a signature of relativistic decoupling, so the candidates lie along the relativistic floor.
-At the threshold m dm ∼ m γ , the annihilation rate into dark photons, which was strongly Boltzmann suppressed at T m dm for m dm < m γ (proportionally to e −2m γ /m dm ), increases very quickly. The resulting efficient annihilation into dark photons leads to a drop in the abundance for fixed ξ, as can be also seen in Fig. 6. Thus, to account for the relic density one needs a larger value of ξ, as displayed in Fig. 5. This marks a sharp transition from the relativistic freeze-out regime to the non-relativistic freeze-out one, i.e. from the relativistic floor value of ξ to the value of ξ needed when the annihilation rate is no longer Boltzmann suppressed by the higher value of m γ . Note that such a threshold feature is not automatically present. Indeed, if the coupling decreases, it may no longer be possible to deplete the abundance through non-relativistic freeze-out. In this case, the curves of the DM candidates lie along the relativistic floor up to a value of m dm (below m γ ), above which DM can no longer be in thermal equilibrium, see below.
-Well above the threshold, m dm > m γ , the DM lies in the secluded, non-relativistic freeze-out regime, Ω dm ∝ m 2 dm /α 2 , see Fig. 6, and the relationship between ξ and m dm evolves similarly to Eq. (27), with ξ dec ∝ α 2 /m 2 dm . The dependence on α stems from the discussion below Eq. (27), so the impact of changing α is manifest. For instance, for ξ < 1, compare this solid line with the part of the black dashed line which lies in the ξ < 1 region: it is obtained for a larger value of α and consequently lies at larger values of m dm .
-Still along these solid and dashed curves, as the DM mass increases and so ξ decreases, FO occurs for larger T dec and smaller x dec = m dm /T dec , see Eq. (24). Eventually, FO takes place when the DM is only mildly non-relativistic, m dm T . This mildly non-relativistic FO is the reason for the upturn of the solid and dashed curves as they go close to the relativistic floor. However, around such masses, thermalisation of the HS is no longer guaranteed, see III A 2 and in particular Fig. 3. Lack of thermal equilibrium implies that the curves stop at some DM mass, see the black dots on the solid and dashed black curves and Fig. 6 of [10]. This dot corresponds to the intersection of the solid curve with the thermalisation diagonal which can be drawn for the given value of α , Eq. (13) (this is not drawn in Fig.  6 but is simply parallel to the orange thermalisation line, given by Eq. (20)).
-The condition for thermal equilibrium scales like ξ min ∝ m dm /α 2 , see Eq. (14). Thus, as α changes, the endpoint runs parallel and close to the relativistic floor. It ends in the lower right corner of the domain (the little red dot) for maximal coupling, where the unitary wall and the thermalisation condition meet. If, on the other hand, we decrease α , the thermalisation endpoint moves up, eventually reaching m dm ≈ m γ . At this point, the condition for thermal equilibrium becomes ξ min ∝ m γ /α 2 (see footnote 5) so that for smaller α the thermalisation endpoints run along the relativistic floor.
-We pointed out in section II that, varying the DM mass and all other things being kept constant (here ξ, m γ and α ), we can in general expect to have an odd number of DM candidates, except at some fined tuned points (corresponding here to the regime m dm ≈ m γ ). This can be seen from the fact that most contours in Fig. 6 cross the Ω dm h 2 = 0.12 line either once or thrice.
2) So far, we assumed that the model parameters are such that the DM candidates lie in the region ξ < 1. Large values α brings in DM candidates for which ξ > 1, see the dashed black curve in Fig. 5. We recognise on this curve several patterns that are similar to the case of smaller α . Starting from the left of the figure, with low mass DM, the abundance along the solid line is set by relativistic decoupling until m dm ≈ m γ , at which point the annihilation channel opens up, leading to a Boltzmann suppression of the relic density that must be compensated by a sharp rise of ξ. Also, as for smaller values of α , at large values of m dm the candidates follow a diagonal line, ξ dec ∝ α 2 /m 2 dm (parallel to the unitarity wall), corresponding to a non-relativistic secluded freeze-out regime. The line is closer of the unitarity wall than the solid line because α is larger.
There is nevertheless a clear difference between the solid and dashed lines for intermediate values of m dm corresponding to the ξ > 1 region. For large α , the line extends itself to this region because as soon as the annihilation channel opens up at m dm ≈ m γ , the large annihilation cross section leads to a large Boltzmann suppression factor which can only be compensated by a value of ξ 1. In this case, however, the DM particle and the dark photons dominate the expansion rate, so the DM abundance no longer depends on the temperature of the visible sector, see section III A 3. The DM mass then depends only on the dark gauge coupling α : Note that the dark photons can generically be made to decay before BBN, so that the domain below the green shaded area but above the green dashed line is potentially allowed. For the two sets of parameters considered here, where m γ = 10 GeV and m γ = 30 MeV, the constraint associated with the γ , ξ m γ /T BBN , is ξ 10 4 and ξ 30 respectively (or a more stringent constraint if the decay of the γ occurs after BBN has started, see e.g. [26][27][28]). The gap in the mass range for candidates for which ξ 1 (see also Fig. 6) could be filled in several ways for this theory. Changing the coupling and the mass of the dark photon is one way. Another way, but a less effective one, is to play with the ratio of degrees of freedom between the VS and HS. We thus conclude that, with an appropriate choice of parameters (possibly all the way to the maximum cross sections allowed by unitarity), the DM candidates of this simple model can fill the whole domain of thermal DM candidates.

B. Scenario 2: s-channel
We next consider a model in which the DM can annihilate into a companion particle through a mediator in the s-channel. It is potentially richer than scenario 1 as it contains an extra coupling and more particles. Nevertheless, most of the features discussed in the case of scenario 1 are similar, so we will be brief. For defi-  niteness, we consider the following Lagrangian, Here the DM is χ, a Dirac fermion, and its companion particles are a scalar, φ, and another Dirac fermion, ψ.
If the ψ is substantially heavier than the DM (and has a subleading contribution to the DM relic density or decays), we have essentially scenario 1, albeit with a spin zero particle instead of a dark photon. The new aspect is that the annihilation of χ into ψ, which is mediated by φ, can be resonant and depends on α x = y χ y ψ /4π. Consider Fig. 7, in particular, the solid curve, which mostly lies in the ξ < 1 domain and for which m φ = 100 GeV, m ψ = 5 GeV and α x = 0.001. The features of this curve are similar to those of the solid curve in Fig. 5 for the t-channel scenario. At low masses m dm < m ψ , the DM candidates lie along the relativistic floor, then there is the threshold effect at m dm ∼ m ψ . The most notable new feature is due to resonant annihilation at m dm ∼ m φ which occurs here in the non-relativistic freeze-out regime. It peaks because the sharp drop of the DM abundance around the resonance (resulting from a sharp rise of the annihilation cross section) must be compensated by an increase of ξ. To the left of the resonance but after threshold, σv ∼ α 2 x m 2 dm /m 4 φ so the curves grow as ξ dec ∝ α 2 x m 2 dm /m 4 φ . To the right of the resonance, the mediator mass becomes less and less relevant and we recover the behaviour already observed in the t-channel scenario, ξ dec ∝ α 2 x /m 2 dm . Again, the curve stops when the DM cannot thermalise. The fea-  tures of the dashed line, obtained for a large coupling value, are similar to those of the dashed line in Fig. 5. It shows a mass gap for the same reasons as discussed for that model.

C. Scenario 3: contact interaction
Finally, we consider a model with no mediator, wherein the DM and companion particle interact via a contact interaction. To that end, we introduce two real scalar fields, Φ and η, which are charged as (+1, −1) and (−1, +1) respectively under a Z 2 × Z 2 symmetry. The potential in the HS is then For concreteness, we take the Φ as the DM candidate (thus the Φ field has no vev). Being scalars, they both could have portal couplings to the Higgs, |H| 2 Φ 2 and |H| 2 η 2 . As before, we assume that these are small enough that the HS does not thermalise with the VS. The features of the ξ(m dm ) contours are similar to the t-channel scenario, see Fig. 8. First, there is a regime of relativistic freeze-out, then a sharp rise at the threshold m dm ≈ m η for ΦΦ → ηη annihilation, and finally the ξ dec ∝ λ 2 Φη /m 2 dm behaviour parallel to the unitarity wall. We depict these for the cases m η = 1 GeV, λ Φη = 1 and m η = 100 GeV, λ Φη = 0.01.
Despite these similarities, this simple model also allows us to consider the following question. So far we looked only at 2 → 2 annihilation processes. Could we go outside the domain by considering processes involving more particles? Although we will not look at this question in full generality, this model illustrates the fact that in general one can expect that the answer to this question is no. Besides the 2 → 2 processes, which can put η and Φ particles in thermal equilibrium, (with rate per unit volume Γ 2↔2 ∝ λ 2 i n 2 i /m 2 dm with i = Φ, η), the DM abundance can also be changed by 4 ↔ 2 processes, with Γ 4↔2 ∝ λ 4 i n 4 i /m 8 dm (we follow the notation of [30], see also [31]). Such processes are slower than the usual 2 ↔ 2 ones because of extra couplings, phase-space and, in the NR regime, Boltzmann factors. However, they can nevertheless play an important role.
First, we can imagine that the companion is simply absent (or heavier than the DM) in which case the HS consists only of the DM, Φ. This is the scenario that was studied in detail in [30]. It is clear that FO cannot occur when the DM is relativistic, since Γ 4↔2 ∝ T for T m dm , rather FO occurs as the rate becomes Boltzmann suppressed for m dm T , see Fig. 3. As there are no DM candidates along the relativistic floor, and no threshold from the companion, all DM candidates, for a given choice of self-coupling, lie (roughly) on a line that runs to the unitarity wall, but this refers to the unitarity of the 4 ↔ 2 process, not the DM freeze-out from a companion. The freeze-out of 4 ↔ 2 processes is distinct from that of the 2 ↔ 2 case. According to the analysis of [30], the relic abundance is determined by Γ 4→2 ∼ n i H. This condition leads to ξ ∝ (λ Φ /m dm ) 4/7 , as opposed to ξ ∝ m −2 dm in the case of 2 ↔ 2 processes, see Eq. (27). The slope is actually much closer to that of the relativistic floor, ξ ∝ m −1/3 dm , see Eq. (9). This stems from entropy conservation in the HS, which leads to reheating of the HS at the same time as the scalar particles deplete their number via self-interactions. The net effect is still a diminution of the particle abundance, but a much less drastic one than in the case of standard non-relativistic freeze-out. Still, the analysis of [30] reveals that that cosmic DM abundance can be reached for a broad range of DM masses, with a maximum possible mass m dm ∼ 10 5 GeV (see Fig. 8 in [30]). Thus, DM candidates of this minimal scenario are well within the domain of thermal DM candidates. They do not reach the unitarity wall (based on 2 ↔ 2 processes), being qualitatively closer to the case of freeze-out in the relativistic regime.
A further interesting feature of this type of scenarios which is worth to point out is that a similar mechanism could lead to the depletion of the companion particles themselves. If m dm m η , such that the DM decouples before the η self-interactions go out of equilibrium, then the situation is precisely the single scalar scenario studied by [30], except in this case the η scalar is not DM and so its abundance should be subdominant. Considering their Fig. 8 and fixing m η = 10 MeV for concreteness, we can see that relatively small ξ are required in order for the η abundance not to be too large, and moreover there is a mild dependence on the quartic coupling. Taking λ η = 10 (about the largest allowed by unitarity), we find that Ω η < Ω dm for ξ 0.1, while taking λ η = 10 −3 (about the smallest that allows η thermalisation), we have Ω η < Ω dm for ξ 0.03. Note that if m dm m η , the results of [30] do not apply, indeed φ − η interactions could further suppress the η abundance, thus in principle allowing larger ξ. These considerations only apply to ξ 1 since the entropy is conserved in the HS. If ξ 1, the DM must have companions (as otherwise its abundance is too large), and they must decay back to the VS.

V. CONCLUSIONS
Thermal relic dark matter candidates may come in many forms, and there is a vast literature concerning this class of models. For scenarios in which the candidate was in equilibrium with the SM bath, the allowed range of dark matter masses is well known, as reviewed in II. It is, however, also plausible that dark matter thermalised within some hidden sector with a temperature, T , different to the temperature of the SM bath, T . This would happen for any hidden sector which involves relatively large interactions between the particles it contains but is connected to the SM thermal bath via significantly weaker interactions. In this paper we studied this general scenario and identified the allowed domain of thermal dark matter candidates in terms of the DM mass and the ratio of temperatures, ξ = T /T . This domain is given in Fig. 2. While parts of this result is implicit in many works, see [7,10,11,23,30], it provides a unifying, and to a large extent modelindependent, picture. In section III A we explored the theoretical bounds which lead to the exclusion of the blue (relativistic floor), red (unitarity wall) and orange (no thermalisation) regions depicted in the figure. Moreover, trying to maintain the generality of our discussion, we placed two observational bounds on the DM, as discussed in section III B. This led to the green (∆N eff ceiling) and purple (free streaming) exclusion regions in Fig. 2.
Putting everything together, we have identified the largest and smallest allowed mass and temperature of thermal DM candidates. The DM mass range, when it is a Dirac fermion, is The possible temperature ratio range is ξ ∈ [1.4 · 10 −5 , 6.9 · 10 5 ] .
In particular, the lower right corner of the domain depicted in Fig. 2 corresponds to a candidate that decoupled while being mildly non-relativistic, m dm ∼ T fo . Hence, the corresponding temperature of the VS at the time of decoupling is bounded from above by T fo ∼ m dm /ξ ∼ 10 11 GeV. Several other features, although rather obvious in ret-rospect, are made clear from our analysis. Firstly, for ξ < 1, the permitted window of DM masses shrinks, and shifts to larger values of m dm , as ξ decreases. Secondly, all other factors being kept constant, the function ξ(m dm ) has in general an odd number of DM candidates, except for very fine-tuned instances. These features are illustrated by three simple models, discussed in section IV. The results, plotted in Figs. 5, 6, 7 and 8, give specific examples of the general findings of Fig. 2.
The different types of HS interactions considered-tchannel, s-channel and contact interaction-illustrate the applicability of our model-independent conclusions. The last model also includes scenarios in which the DM abundance is set by 4 → 2 processes. Several possible developments could be of interest. First, we treated thermal decoupling with a broad brush. Although we do not expect our results to change significantly, it could be interesting to study more carefully and precisely how the DM abundance evolves if it is barely in thermal equilibrium, possibly in the vein of [32]. Second, the observational constraint based on N eff is very conservative and more stringent constraints, especially for candidates which decoupled (or have a companion that decayed back to the VS) around T BBN , could and should be derived. Another possibility is that some candidates have self-interactions that are constrained by e.g. the Bullet cluster [23]. All this is, however, model-dependent and beyond our scope.
Finally, we assumed that the portal interactions between the HS and the VS played little role in determining the relic abundance of the DM. Yet, they could be necessary to get rid of DM companions if their own abundance becomes a nuisance. This is particularly true for ξ > 1 scenarios. We briefly mentioned the impact of such a connection if the mediator decays back to the SM, as for instance studied in [23]. One could also question how our picture changes if a portal interaction leads to a reannihilation regime [7]. In this case, the HS thermalises while DM is being produced from the VS, and DM becomes non-relativistic when the production is still operative. However, one can check 8 that DM candidates produced through reannihilation lie within the domain of thermal DM candidates.

Acknowledgments
This work is supported by the F.R.S./FNRS under the Excellence of Science (EoS) project No. 30820817 8 Reannihilation regimes occur when interactions within the HS thermalise while the energy transfer from the VS to the HS is relevant but comparatively slow. Thus ξ < 1 for reannihilation candidates. For unitarity limited interactions in the HS, from [7,33] one can see that, for fixed ξ < 1, the unitarity wall is (slightly) shifted towards smaller m dm and so the candidates lie within the thermal domain. Concretely, the relic density differs from the one obtained in the non-relativistic secluded freeze-out regime through a shift of x f by ln([ σv portal / σv HS ] 1/2 /ξ 2 ).
-be.h "The H boson gateway to physics beyond the Standard Model", by the "Probing dark matter with neutrinos" ULB-ARC convention, by the FRIA, and by