$Z_2$ SIMP Dark Matter

Dark matter with strong self-interactions provides a compelling solution to several small-scale structure puzzles. Under the assumption that the coupling between dark matter and the Standard Model particles is suppressed, such strongly interacting massive particles (SIMPs) allow for a successful thermal freeze-out through N-to-N' processes, where N dark matter particles annihilate to N' of them. In the most common scenarios, where dark matter stability is guaranteed by a $\mathbb{Z}_2$ symmetry, the seemingly leading annihilating channel, i.e. 3-to-2 process, is forbidden, so the 4-to-2 one dominate the production of the dark matter relic density. Moreover, cosmological observations require that the dark matter sector is colder than the thermal bath of Standard Model particles, a condition that can be dynamically generated via a small portal between dark matter and Standard Model particles, \`a la freeze-in. This scenario is exemplified in the context of the Singlet Scalar dark matter model.


Introduction
One of the major tasks in particle physics is to accommodate dark matter (DM) [1] into the Standard Model (SM) and its extensions. The fact that so far only the gravitational interactions of DM have been confirmed makes this challenging. Among the various strategies used to hunt for DM direct interactions, progress in precision cosmology allows, for the first time, to look for possible signatures of non-gravitational self-interactions among DM particles in the sky. The expectation is that strong enough DM self-interactions, which are mostly likely to be short-ranged, would leave observable imprints on small scales during the growth of cosmic structures. In fact, some small-scale puzzles have been reported, such as the so-called 'cusp vs. core' [2][3][4][5] and the 'too-big-to-fail' [6,7] problems, challenging the standard collisionless cold DM paradigm. Both problems can be addressed by strong DM self-interactions [8][9][10][11][12][13][14], although astrophysical solutions also exist [15][16][17][18]. More recently, a non-vanishing offset between DM and star mass distribution in the cluster Abell 3827 has been reported, also pointing towards to the existence of strong DM self-interactions [19,20].
Strong self-interactions among DM particles would have a significant impact on other aspects of DM phenomenology. For instance, frequent self-scatterings among DM particles modify their phase-space distribution in the Local Group, as well as inside astrophysical objects, leading to novel signatures/interpretations in DM direct and indirect searches [21][22][23][24][25]. More importantly, to achieve sufficient DM self-interactions without contradicting cosmological observations, new mechanisms in DM modeling need to be introduced. 1 The DM self-interactions strength required to alleviate the small-scale structure problems is indeed more than ten orders of magnitude larger than the weak interactions. In this context, the existence of strong DM self-interactions also influences the DM production mechanisms. In particular, the DM thermal freeze-out does not have to be induced by the usual annihilations of two DM particles into a couple of lighter states.
Alternatively, one can consider a framework where the freeze-out proceeds via N -to-N number-changing processes, where N DM particles annihilate into N of them (with N > N ≥ 2), first studied in Ref. [26] and recently named the 'SIMP paradigm' [35]. One potential problem of this mechanism is that when the N -to-N annihilations are effective, DM reheats itself, modifying significantly the formation of structures [38]. Several solutions exist: One can assume either kinetic equilibrium between the DM and the visible sector [35,39], or an enlarged dark sector containing new particles that are relativistic at the moment of the freeze-out [40]. A different approach, adopted in this work, is to assume that the two sectors never reach kinetic equilibrium with each other, and that the DM temperature T is much smaller than the SM temperature T [41]. This can be easily achieved through a very weak portal between DM and SM particles, as shown later.
The most studied cases of the N -to-N processes correspond to 3-to-2 [26,35], because it is typically the dominant one. However, 3-to-2 annihilations, necessarily induced by interaction vertices with an odd number of DM particles, are forbidden in the most common models where the DM stability is guaranteed by a Z 2 symmetry (e.g. R-parity in SUSY or K-parity in Kaluza-Klein scenarios). To allow for 3-to-2 annihilations, one has to assume that DM is protected by a different symmetry such as a Z 3 [39,40,42], or consider models where the DM stability emerges as a result of the DM dynamics [41,[43][44][45][46]. If DM is stabilized by a Z 2 symmetry, the 4-to-2 reactions would be the ones giving rise to the DM relic abundance, while the 3-to-2 annihilations are forbidden. In this work we study such SIMP DM candidates, and illustrate the results in the Singlet Scalar DM model [47,48] as an example.
In Section 2 we briefly review the Singlet Scalar DM model. Section 3 is dedicated to the conditions for this model to alleviate the small scale problems with DM self-interactions. Then we discuss in Section 4 the DM production mechanism compatible with the strong selfinteractions: 4-to-2 DM annihilations, with DM particles colder than the SM. A small Higgs portal can naturally generate such a difference in their temperatures, as shown in Section 5. In Section 6 we briefly discuss ways of mitigating the strong constraints on DM self-interactions from cluster observations, that could allow for lighter DM, before concluding in Section 7.

Singlet Scalar Z Dark Matter
The singlet scalar model [47,48] is one of the minimal extensions of the SM that can provide a viable DM candidate. In addition to the SM framework, this model contains a scalar sector which is only composed by a real scalar S. This particle is a singlet under the SM gauge group, but odd under a Z 2 symmetry, which guarantees its stability. The most general renormalizable scalar potential is given by where H is the SM Higgs doublet. We require that the Higgs gets a non-vanishing vacuum expectation value, v H = 246 GeV, while the singlet does not, S = 0, to ensure the stability of the DM candidate. At tree level, the singlet mass is m 2 S = 2 µ 2 S +λ HS v 2 H . The phenomenology of this model is completely determined by three parameters: the DM mass m S , the Higgs portal λ HS and the quartic coupling λ S . Note that the role of the DM self-coupling λ S was typically disregarded in previous studies.
For the phenomenological study in the following sections, the Lagrangian of the model is implemented in SARAH [71,72] and CalcHEP [73].

Dark Matter Self-Interactions
In the singlet scalar model, self-interactions occur via the contact interaction and the s-, tand u-channel exchange of a Higgs boson (Fig. 1). If DM is light (i.e. m S m h ), the first process typically dominates. In this case, the corresponding cross-section at low redshifts depends only mildly on the non-relativistic DM velocity. In the limit of a small Higgs portal and low velocity, the ratio of the DM self-interaction cross section over the DM mass is given by In order to smooth the inner region of DM halos at scales of few kpc, as required by the 'cusp vs. core' problem, cosmological simulations show that the DM self-interaction needs to lay in the range 0.1 σ SS /m S 10 cm 2 /g, while a smaller value is enough to solve the 'toobig-to-fail' problem [10-14, 74, 75]. On the other hand, too strong self-interactions could cause a collapse in the core of DM halos [10]. Concretely, Bullet Cluster measurements point towards DM self-interactions below 1.25 cm 2 /g at 68% CL [76][77][78]. Moreover, recent observational data on cluster collisions have led to a more stringent bound: σ SS /m S < 0.47 cm 2 /g at 95% CL [79].
Given the large uncertainties in cosmological simulations, we will consider 0.1 σ SS /m S 1 cm 2 /g as the parameter range of interest throughout the paper. Having in mind Eq. (3.1) together with the condition of perturbativity λ S < 4π, the singlet scalar model naturally leads to a DM mass in the sub-GeV region: m S 200 MeV. Such a light DM allows the Higgs to decay into a pair of DM particles. Therefore, both the current limits on the invisible Higgs branching ratio (BR inv 20% [80]) and the Higgs total decay width (Γ tot To summarize, the solution to the small scale problems implies a sub-GeV DM with strong self-coupling λ S ∼ O(1). On top of that, to avoid a too large invisible decay of the Higgs in conflict with the LHC data, a suppressed Higgs portal is also required.

Dark Freeze-out with T < T
In the singlet scalar model, the DM relic abundance can not be generated by the usual pair-annihilating freeze-out, if one wants to solve at the same time the small scale structure problems. In fact, for m S 200 MeV, the kinematically-allowed final states consist of only pions and light fermion pairs (e.g. muons, electrons, etc.). The thermally averaged crosssection in the non-relativistic approximation is which is very suppressed because of the small fermion mass m f compared to the 'heavy' Higgs. For instance, λ HS ∼ 10 −2 and m f = m µ give σv SS→ff ∼ 10 −32 cm 3 /s, far below the canonical value of few 10 −26 cm 3 /s [82], and thus leads to an overclosed Universe. Hence, DM can not be a vanilla thermal WIMP annihilating into a pair of SM particles. Another mechanism for generating the DM relic abundance is therefore needed. Under these conditions, a freeze-out via N -to-N reactions emerges as the dominant DM production mechanism. As DM only annihilates within the dark sector, it is also referred to as 'dark freeze-out' [83,84]. 2 In practice, this scenario can give rise to both large DM self-interactions and the proper DM relic abundance, while being consistent with all experimental constraints. In the singlet scalar model, it happens dominantly via 4-to-2 annihilating channels, as shown in Fig. 2. This is because the Z 2 symmetry forbids the 3-to-2 interactions in the model. Higher-order processes, such as the 5-to-3 or the 6-to-2 annihilations, are subdominant and thus not considered hereafter. The 4-to-2 annihilations have been considered in Ref. [35] assuming that the visible and dark sectors share the same temperature. They concluded that the dark freeze-out requires a too light DM, in tension with cosmological observations. Besides, Refs. [69,70] studied the production of self-interacting singlet scalar DM via the freeze-in without taking into account the effect of 4-to-2 DM annihilations, which is actually important as shown in Section 5.  In the case where the Higgs portal is very suppressed, the Boltzmann equation which describes the evolution of the DM number density n ≡ n (T ) is: where H(T ) is the Hubble rate as a function of visible temperature T and n eq (T ) represents the equilibrium number density of DM particles at a dark temperature T . In the nonrelativistic limit, the 4-to-2 DM annihilation cross section is given by In the Boltzmann equation (4.2), a factor of 2, denoting that each 4 → 2 reaction removes two DM particles, has been absorbed into the definition of the cross section above. This differential equation has been analytically solved in Appendix A, under the assumption that the dark freeze-out happens instantaneously at a visible temperature T FO (or equivalently, T FO for the dark sector), and the results are shown in Fig. 3. The solution is also based on the assumption that the dark freeze-out happens non-relativistically ( where Γ 4→2 = (n eq ) 3 σv 3 4→2 . Conversely, the orange lower band in Fig. 3 yields a (semi-) relativistic freeze-out (i.e. when Eq. (4.4) is not satisfied), which will be addressed afterwards. We emphasize that this inequality requires that the DM particles were in chemical equilibrium between themselves before the freeze-out took place. In Fig. 3, the four black thick lines show the coupling λ S needed to reproduce the observed DM relic abundance through 4-to-2 annihilations. In the case where the dark and the visible sectors are in kinetic equilibrium (thick solid black line T /T | FO = 1), the perturbativity of the λ S coupling implies a DM mass lighter than ∼ 500 keV, in agreement with Ref. [35]. However, this case always produces a too large DM self-scattering. This is shown in Fig. 3 by the three straight diagonal lines, which correspond to the ratio σ SS /m S = 0.1, 1 and 10 cm 2 /g (gray, red and blue lines, respectively), typically needed in order to solve the small-scale problems. The hatched area is disfavored by cluster observations (i.e. σ SS /m S > 1 cm 2 /g).
The constraints from cluster observations become much less severe if the dark temperature T of DM particles is lower than the visible temperature T . The reason is that in the case where T /T | FO 1, the equilibrium DM number density is suppressed by a factor (T /T ) 3 compared to the case where T = T . As a result, in order to obtain the observed DM abundance, an earlier freeze-out and thus a smaller annihilation rate are required. This corresponds to a smaller value for λ S , leading to a weaker DM self-scattering more in line with the cluster bounds. Let us recall that the entropies of the dark and the visible sectors are separately conserved after the moment when the two sectors kinetically decoupled from each other. It is then useful to compute the entropy ratio between them, which is also a conserved quantity. In particular, if the freeze-out happens non-relativistically, one has where s and s are the entropy densities of the SM sector and the dark sector, respectively. This expression was first derived in Ref. [26], and here we add a new term, 2.5 T FO /m S , to take into account second-order corrections. Substituting Ω S h 2 ∼ 0.12 [1] and the solution of x FO in Appendix A into the above equation, we estimate that the proper DM abundance requires ξ ∼ 7 · 10 9 m S GeV The entropy ratio is nearly proportional to the DM mass, up to a weak dependence on the self-coupling λ S . Moreover, to match the observed DM relic abundance there is always a ξ 10 2 for each DM mass. The model therefore never leads to observable extra radiation [26], consistent with BBN/CMB bounds.
At last, we comment the opposite case where the freeze-out happens when the DM is (semi-)relativistic, corresponding to the lower band of Fig. 3. Here the entropy ratio is obtained as ξ ∼ m S 13.5 eV · Ω S h 2 , (4.7) using an analogy with the hot freeze-out of SM neutrinos. The comparison between the previous two cases suggests that Eq. (4.5) can also give a good estimation for the (semi-) relativistic freeze-out by taking T FO ∼ m S .

Generating T < T via the Higgs Portal
In the previous section we studied the dark freeze-out of DM particles colder than the SM sector, without assuming any particular mechanism generating such a difference of temperatures. This could happen even in the absence of the Higgs portal, provided that most DM particles were created at a much earlier time, for instance directly from the inflaton decay [86,87], and then reached a chemical equilibrium within the dark sector. Now we explore the possibility that most of the DM particles were instead created from the SM sector through the Higgs portal. For a very small DM self-coupling, λ S 10 −3 , the freeze-in scenario takes place: once created, DM only redshifts with the expansion of the Universe, without interacting with other particles. Therefore, it is similar to the case of the relativistic freeze-out above, except that now DM particles statistically carry an average kinetic energy similar to the one of the SM particles, i.e. of the order of T . For light DM (m S m h /2), the invisible decay of the Higgs dominates the energy transfer [27,67,68,88], so it is straightforward to obtain (5.1) Nevertheless, since the average kinetic energy of DM particles is given by T , Lyman-α forest surveys apply directly, requiring DM to be heavier than several keV [85]. 3 For higher values of λ S (λ S 10 −3 ), the situation is different. One would expect that if both λ S and λ HS are large enough, after a large amount of DM particles are produced from the SM sector, they can chemically thermalize with themselves: the 'dark thermalization' [89]. In this case the freeze-in scenario does not apply any more. But the condition for obtaining the dark thermalization is not obvious. To simplify the question, we assume that DM particles frequently scatter between themselves, reaching the kinetic equilibrium described by a dark temperature T and a non-vanishing chemical potential µ . However, elastic scatterings do not change the DM number density, so one also needs to consider the leading DM number-changing channel, i.e. the 4↔2 processes, in order to reach chemical equilibrium. The condition for the dark thermalization is therefore given by evaluated just before the 2-to-4 processes become Boltzmann-suppressed. Note that contrary to Eq. (4.4), we are taking here n and not n eq . In practice, this condition is slightly looser than the one of Eq. (4.4), and they mostly coincide. This condition for dark thermalization can be further verified by actually solving the full set of Boltzmann equations, including the Higgs portal interactions, for both the DM number density and its energy density. One numerical solution is depicted in Fig. 4 for a specific parameter set (m S = 1 MeV, λ S = 4 · 10 −3 and λ HS = 6 · 10 −11 ), showing the evolution of both the DM abundance Y S (≡ n /s) (left panel) and its dark temperature T (right panel). The actual evolution is given by the two red lines, labeled Y S and T , respectively. It suggests that even for λ S as low as a few 10 −3 , DM particles are able to thermalize with themselves at a quite early time, and then to freeze-out after becoming non-relativistic, as explained in the four steps below. These four periods are labeled from (a) to (d) in the two panels of Fig. 4. Left: Evolution of DM abundance per comoving volume Y ≡ n /s with respect to the inverse photon temperature x = m S /T . Y S denotes the actual evolution of the DM density and Y eq (T ) the DM equilibrium density, while Y obs. corresponds to the measured relic abundance Ω S h 2 ∼ 0.12 [1]. Y freeze-in gives the abundance in the case where the 4-to-2 processes are ignored. Right: Evolution of the DM temperature T with respect to x. The SM temperature T and the DM mass m S (dotted horizontal line) are provided for comparison. The parameters used are m S = 1 MeV, λ S = 4 · 10 −3 and λ HS = 6 · 10 −11 (leading to ξ ∼ 4 · 10 5 ). Three vertical lines stand for the times when DM reaches chemical equilibrium (x ∼ 2 · 10 −5 ), DM becomes non-relativistic (x ∼ 0.04) and DM freezes-out (x ∼ 0.3).
• During the period (a), the first DM particles are created via out-of-equilibrium decays of the Higgs (h → SS). The 2-to-4 processes are originally very suppressed by the small DM number density and the high velocity of the DM particles (initially T and T are of the same order). However these suppressions become less severe due to the increase of the DM population and to the cooling down of the dark sector. At some point the DM production via the 2-to-4 interactions overtakes the one due to the Higgs decay, giving rise to a dramatic increase of the DM abundance at the cost of reducing its temperature. At the end of this period DM reaches chemical equilibrium.
• After reaching chemical equilibrium among themselves (x ∼ 2 · 10 −5 ), in period (b) the actual DM abundance Y (T ) tracks the equilibrium density Y eq (T ). Moreover, as DM is relativistic during this period, T scales as T .
• When DM becomes non-relativistic (T < m S ) in period (c), its number density becomes Boltzmann suppressed. During this period, due to the conservation of the DM entropy, T scales as 1/ log a [26], where a is the scale factor. This gives rise to a relative increase of T compared to T .
• Eventually when x ∼ 0.3 (or equivalently x ∼ 7) DM freezes-out, leading to a constant value for Y (T ), which characterizes period (d). During this era, DM is non-relativistic and out of chemical equilibrium, so that T decreases as T 2 .
In the left panel of Fig. 4, the blue curve labeled as 'freeze-in' corresponds to the case where the DM abundance is solely produced via Higgs decays (hence corresponding to the freeze-in mechanism), because the 4-to-2 processes (and then the dark thermalization) were artificially ignored. In this case the DM relic abundance is largely underestimated. It is worth mentioning that the importance of the 4-to-2 processes is due to the large value for λ S , required for the solution of the small-scale structure problems. Therefore, in the white region of Fig. 5, above the '(Semi-)relativistic freeze-out' band, the 4-to-2 processes have to be taken into account. This is why the freeze-in mechanism alone can not give simultaneously a large DM self-scattering and the correct DM abundance for m S O(1) MeV in a consistent way, in contrast to previous works which neglected the 4-to-2 processes [69,70].
Under the assumption that all the dark sector was originally created by energy-transfer from the visible sector through the Higgs portal, and taking into account that this transfer occurs mainly around T ∼ m h /3 m S [88], we calculate the ratio of energy densities between the two sectors after the Higgs decay: where g SM * (T ) is the effective degrees of freedom of the visible sector and n h eq the equilibrium Higgs number density. This expression leads to ξ ∼ 2 · 10 5 · 10 −10 λ HS which is independent of m S in the relativistic limit. Eq. (5.4) together with Eq. (4.6) allow to estimate the strength of Higgs portal required to provide the observed DM relic density. This is illustrated in Fig. 5 for various values of λ HS . Besides, a DM dark freeze-out is possible before the complete decay of the Higgs population. Roughly, it corresponds to T FO m h /3. In this case, the (late) Higgs decay into DM particles after T FO increases the DM relic abundance, and thus resumes the DM annihilations. This is commonly referred to as the 'reannihilation' mechanism [88,90,91]. However, we verified that this possibility never takes place in the parameter region of interest here.

Experimental Detection
Sub-GeV DM that couples to the SM via a (very) suppressed Higgs portal is very challenging to detect, both for next generation of DM direct detection [92][93][94][95][96] or indirect detection experiments [97,98]. For instance, the experimental signature of such DM particles in direct detection searches using electron recoil events is suppressed by the smallness of both Higgs portal λ HS and electronic Yukawa coupling y e .
In the context of astrophysical observations, although there exist stringent bounds on Yukawa-like couplings between SM fermions and light (meta-)stable new scalars based on stellar energy-loss arguments [99], they do not directly apply to Z 2 -symmetric DM particles, which can only be created in pairs. Moreover, such DM particles may be trapped in stars by strongly interacting with environmental DM particles without causing any noticeable energy loss [21,100,101].

A Caveat Leading to Lighter Dark Matter
From Figs. 3 and 5, it is clear that cluster bounds on self-interactions (together with the requirement of a non-relativistic freeze-out) impose a strong constraint that makes DM mass lie in the MeV range. In this section, we briefly discuss scenarios where cluster bounds become less severe, opening the possibility of lighter (warm) DM freezing-out through the 4-to-2 processes.
One generic class of models which are able to avoid large self-interaction is the so-called excited (or inelastic) DM [102][103][104][105]. Here we take as an example a pseudo-Dirac fermion invariant under a U (1) D gauge symmetry [34,103,105]. At low energies, after the symmetry breaking of the U (1) D , a small Majorana mass splits the Dirac fermion into two nearlydegenerated Majorana states χ 1 and χ 2 . The physical spectrum contains a light state χ 1 , which is the DM candidate (stabilized by an accidental Z 2 symmetry), an excited partner χ 2 slightly heavier and a massive gauge boson A D . Due to the Majorana nature of the χ i , gauge interaction terms like χ 1 χ 1 A D or χ 2 χ 2 A D are not allowed, while χ 1 χ 2 A D does exist. Therefore, provided that U (1) D mixes with the SM hyperchange, χ 2 decays to χ 1 and a photon, via the kinetic mixing with A D . The gauge boson A D is also unstable and decays to χ 1 and χ 2 , or SM particles.
If A D is heavier than both χ 1 and χ 2 , the DM abundance can not be generated via the usual freeze-out into SM particles, because its annihilation cross section is extremely suppressed by the tiny kinetic mixing coupling (typically below 10 −10 for A D mass in the MeV range [106,107]). Nevertheless, if the gauge coupling of U (1) D is large enough, after the decoupling of A D the chemical equilibrium among χ 1 and χ 2 can still be maintained via (co-)annihilation processes, such as 2χ i + 2χ j → 2χ i , 3χ i + χ j → χ i + χ j and 4χ i → 2χ j , with i = j. Therefore, the DM abundance can be set by a 4-to-2 dark freeze-out scenario. It is similar to the one described above for the singlet scalar model, although the strength of these processes gets moderately suppressed by the existence of heavier A D .
The situation of self-scatterings for excited DM, on the other hand, differs substantially from the singlet scalar case, due to the absence of the interaction vertex χ 1 χ 1 A D . As a consequence, DM elastic scatterings χ 1 χ 1 → χ 1 χ 1 can only happen via loops, and are thus strongly suppressed [34]. The inelastic scattering χ 1 χ 1 → χ 2 χ 2 at low temperatures is also suppressed as it is endothermal. Therefore, cluster bounds on DM self-interactions can be easily lifted without overturning previous results, such as the dark freeze-out via 4-to-2 processes. Such an alleviation of DM self-interaction bounds in excited DM models helps to extend the allowed parameter region. Similar conclusions can also be reached in models where DM elastic scattering is p-wave suppressed. Consequently, both phenomenological classes of models allow for light DM candidates, even much below MeV. In fact, for a DM mass around O(30) keV, the damping effect induced by its free-streaming may mimic that of conventional keV-scale warm DM on structure formation.

Conclusions
The collisionless cold dark matter (DM) paradigm has been highly successful in accounting for large scale structure of the Universe. However, it is not clear that it can also fully explain the small scales. In fact, long-standing puzzles from observations of small scale objects (e.g. the 'cusp vs. core' and the 'too-big-to-fail' problems), can be alleviated if large DM selfinteractions at scales below few Mpc exist.
The singlet scalar DM model is a very simple and economical model, described by only three parameters: the DM mass m S , the DM quartic coupling λ S and the Higgs portal λ HS . In this framework, imposing strong DM self-interaction naturally points toward a sub-GeV m S and a self-coupling λ S ∼ O(1). Moreover, non-observation of the invisible decay of Higgs requires a suppressed Higgs portal. Given these conditions, it is challenging to produce the measured DM relic abundance through the standard freeze-out mechanism. Indeed, the crosssection of DM annihilating into SM particles is very suppressed because of the smallness of both the Higgs portal and the Yukawa couplings of the relevant SM fermions.
An alternative production mechanism is to consider a dark freeze-out via N -to-N processes, where N DM particles annihilate into N of them (with N > N ≥ 2), namely the SIMP paradigm. The most studied cases correspond to 3-to-2 reactions. However, for models where the DM stability is guaranteed by a Z 2 symmetry, the dominant annihilation channels are 4-to-2 processes (the 3-to-2 being forbidden). Moreover, to be consistent with structure formation, DM particles colder than their SM counterparts are required before the dark freezeout. In the singlet scalar model, it can be easily achieved by introducing a very suppressed Higgs portal.
With this in mind, we first considered the dark freeze-out of DM particles via 4-to-2 annihilations for a dark sector colder than the SM, without assuming any particular mechanism to generate such a difference of temperatures. By analytically solving of the corresponding Boltzmann equation, we have shown that for each parameter set of (λ S , m S ), there is a certain value of temperature ratio T /T , which generates both the measured DM relic abundance and the self-interactions needed to address the small-scale problems. Typically, a ratio T /T ∼ O (20) at the moment of the freeze-out is needed (cf. Fig. 3).
As a second step, we studied how to generate a colder dark sector under the assumption that most of the DM particles were initially created out of the SM through the Higgs portal, à la freeze-in. We found that for most of the parameter region considered, large values of λ S coupling are able to bring DM particles into both kinetic and chemical equilibrium with themselves thanks to the fast DM elastic scatterings and the 4 ↔ 2 interactions. Therefore, in this setup one naturally obtains T < T due to the tiny portal λ HS . Typical values of the Higgs portal, e.g. λ HS ∼ O(10 −12 − 10 −10 ), needed to yield the measured DM abundance are shown in Fig. 5 for the dark freeze-out scenario.
At last, it is worth pointing out that the essence of our conclusion is not exclusive to the singlet scalar case, but general to a specific pattern of models where DM particles solely interact with the SM via a Higgs portal: In order to provide a possible solution to small-scale problems, such models prefer a sub-GeV DM with strong self-interactions but weakly connected to the SM. Together with relevant experimental constraints, this challenges the usual WIMP freeze-out, but favors the dark freeze-out scenario via 4-to-2 (or 3-to-2 if allowed) annihilations of DM. In this scenario, the measured DM relic abundance is produced by having a DM sector colder than the SM one at the moment of freeze-out, which is in turn provided by the very suppressed Higgs portal.

A Freeze-out Approximation
In this Appendix we closely follow the procedure presented in Ref. [41].
After the dark and visible sectors completely decouple (i.e. the energy density transfer via the Higgs decay h → SS is not efficient anymore), the evolution of the DM number density n ≡ n (T ) is given by the Boltzmann equation: which, in the case where the SM energy density dominates the expansion of the Universe, can be rewritten in terms of the SM temperature T as where x ≡ m S /T and Y (x) ≡ n /s(x).
In the non-relativistic approximation, the cross-section σv 3 4→2 is independent of the temperature, and then the solution to Eq. (A.2) that matches the observed DM relic abundance reads Notice that the temperature that enters in this expression is the one of the SM and not the DM one. In order to estimate x FO , it is necessary to establish when the annihilation rate per particle (n eq ) 3 σv 3 4→2 drops below the expansion rate of the Universe. Using Eq. (A.3) it is found that this happens when the freeze-out temperatures satisfy . (A.4) The solution to this equation is shown in Fig. 6, for different temperature ratios at the moment of the DM freeze-out. For m S = 100 keV and assuming kinetic equilibrium between the two sectors until the freeze-out, x FO = x FO ∼ 13.5 which agrees with the result from Ref. [35].