Phenomenology of self-interacting dark matter in a matter-dominated universe

Abstract

We study production of self-interacting dark matter (DM) during an early matter-dominated phase. As a benchmark scenario, we consider a model where the DM consists of singlet scalar particles coupled to the visible Standard Model (SM) sector via the Higgs portal. We consider scenarios where the initial DM abundance is set by either the usual thermal freeze-out or an alternative freeze-in mechanism, where DM was never in thermal equilibrium with the SM sector. For the first time, we take the effect of self-interactions within the hidden sector into account in determining the DM abundance, reminiscent to the Strongly Interacting Massive Particle (SIMP) scenario. In all cases, the number density of DM may change considerably compared to the standard radiation-dominated case, having important observational and experimental ramifications.

1 Introduction

The existence of dark matter (DM) seems indisputable. From the Cosmic Microwave Background radiation (CMB), large scale structure of the Universe and different physics at galactic scales, one can infer that there must be a long-lived, dynamically non-hot, non-baryonic matter component, whose abundance exceeds the amount of ordinary ‘baryonic’ matter roughly by a factor of five [1,2,3,4] and which has been there from the hot Big Bang era until the present day. However, the non-gravitational nature of the DM component remains a mystery.

For a long time, Weakly Interacting Massive Particles (WIMPs) have been among the best-motivated DM candidates. The increasingly strong observational constraints on DM (see e.g. Ref. [5]) are, however, not only puzzling as such but are now forcing one to ask: is the standard WIMP paradigm just waning, or is it already dead? If so, what alternative explanations for the production and properties of DM do we have?

A simple alternative for the standard WIMPs is provided by relaxing the usual assumption that DM is a thermal relic, produced by the freeze-out (FO) mechanism in the early Universe, and assuming that it never entered in thermal equilibrium with the particles within the Standard Model of particle physics (SM). If that was the case, then the present DM abundance could have been produced by the so-called freeze-in (FI) mechanism, where the abundance results from decays and annihilations of SM particles into DM [6,7,8,9,10]. Assuming that DM never entered into thermal equilibrium with the particles in the visible SM sector typically amounts to choosing a very small coupling between the two sectors. A good thing about this is that then these so-called Feebly Interacting Massive Particles (FIMPs) easily evade the increasingly stringent observational constraints, yet an obvious hindrance is that this also makes the scenario inherently very difficult to test. For a recent review of FIMP DM models and observational constraints presented in the literature, see Ref. [11].

Another way to evade the experimental constraints is to consider non-standard cosmological histories [12]. We know that the Universe was effectively radiation-dominated (RD) at the time of Big Bang Nucleosynthesis (BBN) and one usually assumes that this was the case also at the time the DM component was produced, was it at the time of electroweak cross-over or at higher energy scales. However, there are no obvious reasons for limiting the DM studies on such cosmological expansion histories,¹ as alternatives not only can lead to interesting observational ramifications but are also well-motivated. For example, an early matter-dominated (MD) phase can be caused by late-time reheating [18], massive meta-stable particles governing the energy density of the Universe (see Refs. [19,20,21] for recent works), moduli fields [22,23,24], and so on. The effect on the resulting DM yield can then be outstanding, as recently studied in detail in e.g. Refs. [19,20,21, 25,26,27,28,29,30,31,32].

Indeed, when the expansion rate of the Universe differs from the usual RD case, it tends to effectively dilute the DM abundance when the era of non-standard expansion ends and the visible sector gets reheated (see also Refs. [27, 31] for DM production in fast-expanding universes and Refs. [26, 30] for co-decaying DM). This means, for example, that when the expansion was faster than in the RD case and the DM particles were initially in thermal equilibrium with the visible sector, they generically have to undergo freeze-out earlier than in the usual RD case, thus resulting in larger DM abundance to match the observed one. In case the DM particles interacted so feebly that they were never part of the equilibrium heat bath, the coupling between DM and the visible sector typically has to be orders of magnitude larger than in the usual freeze-in case to compensate the larger expansion rate. Production of DM during a non-standard expansion phase may thus result to important experimental and observational ramifications. Studying the effect non-standard cosmological histories have on different particle physics scenarios is thus not only of academic interest and also not limited to the final DM abundance, as different possibilities to test for example an early MD phase include formation of ultracompact substructures such as microhalos [33] or primordial black holes [34,35,36], as well as cosmological phase transitions with observational gravitational wave signatures [37] (see also Ref. [38]).

In this paper we will consider DM production during such an early MD phase. We will study DM production by both the freeze-out and freeze-in mechanisms, taking for the first time into account the effect that non-vanishing DM self-interactions can have. Instead of performing an intensive full-parameter scan, in this paper we will perform an analytical study of the different representative cases previously mentioned, which allows us to capture the essence of each scenario. Results of an exhaustive scan over the full parameter space in the usual freeze-out and freeze-in cases are presented in a companion paper [39], where we also discuss the effect of other non-standard cosmological histories. However, as we will show, already with the best-motivated non-standard case, an early phase of matter-domination, the DM phenomenology is very rich when the effect of DM self-interactions is taken into account, which is one of the reasons why we devote a separate paper for the analysis of this scenario only. Another important difference to Ref. [39] is that in this paper we will we make the usual assumption that the eventual decay of the energy density component responsible for the early matter-domination is instantaneous, whereas in Ref. [39] the duration of decay is taken to be finite. In this way, the two studies complement each other.

As we will show, the observational limits on DM self-interactions do not only rule out part of the parameter space for the model we will consider in this paper, but taking the detailed effect of DM self-interactions into account is crucial for determination of the final DM abundance, reminiscent to the so-called Strongly Interacting Massive Particle (SIMP) or cannibal DM scenarios [40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63]. We will also discuss other prospects for detection of DM including collider, direct, and indirect detection experiments.

The paper is organized as follows: In Sect. 2, we will present a simple benchmark model where the DM particle is a real singlet scalar odd under a discrete $\mathbb {Z}_2$ symmetry, and discuss what are the requirements for having an early MD phase prior to BBN. In Sect. 3, we turn into the DM production, discussing production by the usual freeze-out mechanism in Sect. 3.1 and by the freeze-in mechanism in Sect. 3.2. In Sect. 4, we discuss the experimental and observational ramifications, and present not only what part of the parameter space is already ruled out but also what part of it can be probed in the near future. Finally, we conclude with an outlook in Sect. 5.

2 The model

We study an extension of the SM where on top of the SM matter field content we assume a simple hidden sector consisting of a real singlet scalar s. The only interaction between this hidden singlet sector and the visible SM sector is via the Higgs portal coupling $\lambda _{hs}|\varPhi |^2s^2$, where $\varPhi$ is the SM Higgs field. The scalar potential is

$$\begin{aligned} V(\varPhi ,s) = \mu _h^2|\varPhi |^2{+}\lambda _h|\varPhi |^4+\frac{\mu _s^2}{2} s^2+\frac{\lambda _s}{4}s^4{+}\frac{\lambda _{hs}}{2}|\varPhi |^2 s^2 ,\nonumber \\ \end{aligned}$$

(1)

where $\sqrt{2}\varPhi ^\mathrm{T}=(0,v+h)$ is the SM SU(2) gauge doublet in the unitary gauge and $v=246$ GeV is the vacuum expectation value of the SM Higgs field. A discrete $\mathbb {Z}_2$ symmetry, under which the DM is odd and the whole SM is even, has been assumed to stabilize the singlet scalar and make it a possible DM candidate. We assume $\lambda _s>0$ and $\mu _s>0$, so that the minimum of the potential in the s direction is at $s=0$ and $m_s^2 \equiv \mu _s^2 + \lambda _{hs}\,v^2/2$ is the physical mass of s after the spontaneous symmetry breaking in the SM sector. This implies $\lambda _{hs}<2\,m_s^2/v^2$.

2.1 An early matter-dominated period

We assume that the Universe was MD for the whole duration of DM production down to $T \gtrsim 4$ MeV, where the lower limit is given by BBN [64,65,66,67]. By this time, the matter-dominance must have ended, the SM sector must have become the dominant energy density component and the usual Hot Big Bang era must have begun. We assume that when DM was produced, both the SM and the singlet sector were energetically subdominant, so that

$$\begin{aligned} 3\,H^2 M_\mathrm{P}^2 = \rho _\mathrm{total} \simeq \rho _\mathrm{M} \gg \rho _\mathrm{SM},\,\rho _s , \end{aligned}$$

(2)

where H is the Hubble scale, $M_\mathrm{P}$ is the reduced Planck mass, and $\rho _\mathrm{M}$ is the energy density of the matter-like component that is assumed to dominate over the SM energy density $\rho _\mathrm{SM}$ and the singlet scalar energy density $\rho _s$. We also assume that the SM was in thermal equilibrium for the whole duration of the early MD phase, so that

$$\begin{aligned} \rho _\mathrm{SM} = \frac{\pi ^2}{30}\,g_* \,T^4, \end{aligned}$$

(3)

where $g_*$ is the usual effective number of relativistic degrees of freedom² and T is the SM bath temperature.

The magnitude of the Hubble expansion rate can be understood by first discussing the dynamics in the usual RD case where the SM is the dominant energy density component. In that case, the Friedmann equation (2) gives at $T=m_h$ the result

$$\begin{aligned} \frac{H_\mathrm{EW}^\mathrm{rad}}{m_h} = \sqrt{\frac{\pi ^2g_*(m_h)}{90}}\frac{m_h}{M_\mathrm{P}} \simeq 1.76\times 10^{-16} , \end{aligned}$$

(4)

where we used $g_*(m_h)=106.75$ and denoted $H_\mathrm{EW}\equiv H(T=m_h)$. However, in a MD Universe at $T=m_h$ we have

$$\begin{aligned} 3\,H_\mathrm{EW}^2\,M_\mathrm{P}^2 = \left. \left( \rho _\mathrm{M} + \rho _\mathrm{SM}\right) \right| _{T=m_h} \simeq \left. \rho _\mathrm{M}\right| _{T=m_h}, \end{aligned}$$

(5)

so that in this case $H_\mathrm{EW}/m_h\gg H_\mathrm{EW}^\mathrm{rad}/m_h$, i.e. the Universe expands much faster than in the standard RD case. Determining the ratio $H_\mathrm{EW}/m_h$ more accurately than this is not possible without specifying the underlying dynamics causing the early MD, so in the remaining of this paper we simply take it to be a free parameter for generality.

2.2 Constraints on the scenario

In all cases, both the model parameters in Eq. (1) and the cosmological parameters are subject to constraints that come from observational data. In this paper, we make the usual assumption that the matter component governing the total energy density decays instantaneously into the SM radiation. The first condition then is that the SM temperature after the matter-like component has decayed into SM particles, $T_\mathrm{end}'$, must be larger than the BBN temperature $T_\mathrm{BBN}=4$ MeV. Second, the temperature has to be smaller than either the final freeze-out temperature or smaller than $m_h$ in the freeze-in case in order not to re-trigger the DM yield after the decay of the matter-like component. As shown in the end of Appendix A, this amounts to requiring

$$\begin{aligned}&5\times 10^{-7}\left( \frac{H_{\mathrm{EW}}/m_h}{10^{-16}}\right) ^{-2/3} \nonumber \\&\quad \lesssim \frac{T_{\mathrm{end}}}{m_h} \lesssim {\left\{ \begin{array}{ll} 2\times 10^{-3}\left( \frac{H_{\mathrm{EW}}/m_h}{10^{-16}}\right) ^{-2/3}\left( \frac{m_s}{\mathrm{GeV}}\right) ^{4/3}x_{\mathrm{FO}}^{-4/3} \quad \text {freeze-out}, \\ \left( \frac{H_{\mathrm{EW}}/m_h}{10^{-16}}\right) ^{-2/3} \quad \text {freeze-in}, \end{array}\right. }\nonumber \\ \end{aligned}$$

(6)

where $T_\mathrm{end}$ is the SM temperature just before the end of matter-domination and $x_\mathrm{FO}\equiv m_s/T_\text {FO}$, with $T_\text {FO}$ being the DM freeze-out temperature. In the following, we will take the above ratio $T_\mathrm{end}/m_h$ to be a free parameter, so that together with $H_\mathrm{EW}$ it constitutes the set of our cosmological parameters, characterizing the duration of the early MD phase. The total parameter space is thus five-dimensional, consisting of the particle physics parameters $\lambda _s$, $\lambda _{hs}$ and $m_s$, in addition to the cosmological parameters $H_\text {EW}/m_h$ and $T_\text {end}/m_h$.

Third, we require that DM freeze-out always occurs while the s particles are non-relativistic, $x_\mathrm{FO}>3$, as otherwise the scenario is subject to relativistic corrections that we are not taking into account in the present paper. Fourth, as discussed above, in a MD Universe $H_{\mathrm {EW}}/m_h\gg 10^{-16}$. Fifth, as discussed below Eq. (1), the portal coupling has to satisfy $\lambda _{hs}<2\,m_s^2/v^2$. Finally, the portal coupling has a further constraint when requiring or avoiding the thermalization of the two sectors, for the case of freeze-out and freeze-in, respectively. Depending on the strength of the portal coupling $\lambda _{hs}$, the singlet scalar particles may or may not have been part of the equilibrium in the SM sector at the time the initial DM density was produced. The threshold value for $\lambda _{hs}$ above which the DM sector equilibrates with the SM is

$$\begin{aligned} \lambda _{hs}^\text {eq} \simeq \sqrt{\frac{128\pi ^3}{\zeta (3)}\frac{H_\mathrm{EW}}{m_h}}. \end{aligned}$$

(7)

This results from requiring that the SM particles do not populate the hidden sector so that they would start to annihilate back to the SM in large amounts, $\langle \sigma _{hh \rightarrow ss}v\rangle n_h/H \simeq \lambda _{hs}^2\,\zeta (3)\, m_h/(128\pi ^3H_\text {EW}) < 1$ [68,69,70], where $\langle \sigma _{hh\rightarrow ss}v\rangle$ is the thermally averaged cross-section for the process $hh\rightarrow ss$ and $\zeta (3)\simeq 1.20$ is the Riemann zeta function. For the freeze-out case we demand $\lambda _{hs}\gg \lambda _{hs}^\text {eq}$ whereas for the freeze-in $\lambda _{hs}\ll \lambda _{hs}^\text {eq}$.

Before concluding this section let us note that the fact that now $H_\mathrm{EW}\gg H_\mathrm{EW}^\mathrm{rad}$ means that in the freeze-out case the value of the portal coupling required to produce the observed DM abundance must be smaller than in the usual RD case, as the DM has to decouple earlier from the thermal bath in order to retain the required abundance. However, the faster expansion rate also means that now the threshold value for thermalization, Eq. (7), can be orders of magnitude larger than the corresponding value $\lambda _{hs}\simeq 10^{-7}$ in the usual RD case. This makes the freeze-in scenario particularly interesting, as it might lead to important experimental ramifications, as we will discuss in Sect. 4.

3 Dark matter production

We start by reviewing the DM production within this model, briefly discussing two fundamental mechanisms that account for it: the freeze-out and the freeze-in scenarios.

Assuming that there is only one DM particle, s, its number density evolution is described by the Boltzmann equation:

$$\begin{aligned}&\frac{dn_s}{dt}+3\,H\,n_s= -\int d\varPi _{s}d\varPi _{\mathrm{a_{1}}}d\varPi _{\mathrm{a_{2}}}\cdots d\varPi _{b_{1}}d\varPi _{\mathrm{b_{2}}}\cdots \nonumber \\&\quad \times \left( 2\pi \right) ^{4}\delta ^{4}\left( p_{s}+p_{\mathrm{a_{1}}}+p_{\mathrm{a_{2}}}...-p_{\mathrm{b_{1}}}-p_{\mathrm{b2}}\cdots \right) \nonumber \\&\quad \times \left[ \left| \mathcal {M}\right| _{\mathrm {s+a_{1}+a_{2}\cdots \rightarrow b_{1}+b_{2}\cdots }}^{2}\,f_{s}\,f_{\mathrm {a_{1}}}\cdots \left( 1\pm f_{\mathrm {b_{1}}}\right) \right. \nonumber \\&\quad \left. \times \left( 1\pm f_{\mathrm {b_{2}}}\right) \cdots -\left| \mathcal {M}\right| _{\mathrm {b_{1}+b_{2}\cdots \rightarrow s+a_{1}+a_{2}\cdots }\,}^{2}f_{\mathrm {b_{1}}}\right. \nonumber \\&\quad \left. \times \,f_{\mathrm {b_{2}}}\cdots \left( 1\pm f_{s}\right) \left( 1\pm f_{\mathrm {a_{1}}}\right) \cdots \right] , \end{aligned}$$

(8)

considering the process $s+a_{1}+a_{2}+\cdots +a_{k}\rightarrow b_{1}+b_{2}+\cdots +b_{j}$, where $a_i, b_j$ are particles in the heat bath. Here $n_s$ is the DM number density, $p_i$ is the momentum of the particle i, $\left| \mathcal {M}\right| ^{2}$ is the squared transition amplitude averaged over both initial and final states, $f_i$ is the phase space density, $+$ applies to bosons and − to fermions and

$$\begin{aligned} d\varPi _i\equiv \frac{g_i}{\left( 2\pi \right) ^{3}}\,\frac{d^{3}p_i}{2E_i} \end{aligned}$$

(9)

is the phase space measure, where $g_i$ is the number of intrinsic degrees of freedom and $E_i$ the energy of the particle i. In the following, we will solve the relevant Boltzmann equations analytically in the regions of interest where different processes dominate at a time. A full parameter scan is performed in the pure freeze-out and freeze-in cases in Ref. [39].

In the freeze-out mechanism, DM was initially in thermal equilibrium with the SM sector. As soon as the interactions between the DM and the SM particles were no longer able to keep up with the Hubble expansion, the system departed from thermal equilibrium and the comoving DM abundance became constant. We will study the case of the DM freeze-out in an early MD era in Sect. 3.1.1 and then consider how a so-called cannibalism phase affects the DM yield in Sect. 3.1.2.

In the freeze-in scenario, the DM was never in thermal equilibrium with the visible sector, due to the very feeble interactions between them. The particles produced by this mechanism are known as FIMPs and their initial number density is, in the simplest case, negligible. The DM abundance is produced by the SM particle decays and annihilations, lasting until the number density of the SM particles becomes Boltzmann-suppressed. At this point, the comoving number density of DM particles becomes constant and the comoving DM abundance is said to ‘freeze in’. The evolution of the initial s number density can be tracked by the Boltzmann Eq. (8) as well. We discuss the DM freeze-in in an early MD era without cannibalism in Sect. 3.2.1 and with it in Sect. 3.2.2.

3.1 The freeze-out case

To study the effects of MD and DM self-interactions in a simple yet accurate way, in this section we assume the mass hierarchy $m_b<m_s < 50$ GeV, where $m_b$ is the mass of the b-quark and the upper limit is chosen to avoid complications with the Higgs resonance in our analytical calculations. Therefore, in this subsection, we will consider DM produced only by $b\bar{b}$ annihilations and present the more general analysis in Ref. [39] for the pure freeze-out case without cannibalism.

3.1.1 Freeze-out without cannibalism

In this scenario, we assume that the DM was initially in thermal equilibrium with the SM particles. In the most simple case that we are considering here, only the annihilation and inverse annihilation processes $ss\leftrightarrow b\bar{b}$ are taken into account for the abundance, and the equation governing the evolution of the DM number density, (8), becomes

$$\begin{aligned} \frac{dn_{s}}{dt}+3\,H\,n_{s}=-\left\langle \sigma _{ss\rightarrow b\bar{b}} v\right\rangle \left[ n_{s}^{2}-\left( n_{s}^\text {eq}\right) ^{2}\right] , \end{aligned}$$

(10)

where $\langle \sigma _{ss\rightarrow b\bar{b}} v\rangle$ is the thermally-averaged DM annihilation cross-section times velocity and $n_{s}^\text {eq}$ corresponds to the DM equilibrium number density.

When the interactions between the DM and the visible sector cannot keep up against the expansion of the Universe any more, the DM decouples and its comoving number density freezes to a constant value. This occurs at $T=T_\text {FO}$ defined by

$$\begin{aligned} \left. \frac{\left\langle \sigma _{ss\rightarrow b\bar{b}}v\right\rangle \,n_{s}}{H}\right| _{T=T_\text {FO}}=1. \end{aligned}$$

(11)

Assuming that DM is non-relativistic when interactions freeze-out, we have

$$\begin{aligned} n_{s}(T)=\left( \frac{m_{s}\,T}{2\pi }\right) ^{\frac{3}{2}}\,e^{-\frac{m_{s}}{T}}, \end{aligned}$$

(12)

whereas the Hubble parameter is given by

$$\begin{aligned} H(T)=H_{\mathrm {EW}}\,\left( \frac{T}{m_h}\right) ^{\frac{3}{2}}\,\left( \frac{g_{*}\left( T\right) }{g_{*}\left( m_h\right) }\right) ^{\frac{1}{2}}. \end{aligned}$$

(13)

Substituting then Eqs. (12) and (13) into (11), the freeze-out condition can be written as

$$\begin{aligned} x_{\mathrm {FO}}= & {} \ln \left[ \frac{\lambda _{hs}^{2}}{2^{9/2}\,\pi ^{5/2}}\,\left( \frac{g_{*}\left( m_h\right) }{g_{*}\left( T_{\mathrm {FO}}\right) }\right) ^{1/2}\,\right. \nonumber \\&\quad \left. \times \left( \frac{H_{\mathrm {EW}}}{m_h}\right) ^{-1}\,\frac{m_b^{2}\,m_{s}^{3/2}}{m_h^{7/2}}\right] , \end{aligned}$$

(14)

where we used $\langle \sigma _{ss\rightarrow b\bar{b}}v \rangle \simeq \lambda _{hs}^2m_b^2/(8\pi \, m_h^4)$ [43, 44] and $x_{\mathrm {FO}}\equiv m_s/T_\mathrm{FO}$ corresponds to the time when DM annihilation into b-quarks becomes smaller than the Hubble parameter. The DM abundance can then be calculated by taking into account the non-conservation of entropy (see Appendix A), yielding:

$$\begin{aligned} \frac{\varOmega _{s}\,h^{2}}{0.12}&\simeq \,3\times 10^{-7}\,x_{\mathrm{FO}}^{3/2}\,e^{-x_{\mathrm{FO}}}\,\nonumber \\&\quad \times \left( \frac{T_{\mathrm{end}}}{m_h}\right) ^{3/4}\,\left( \frac{H_\mathrm{EW}/m_h}{10^{-16}}\right) ^{-3/2}\,\left( \frac{m_{s}}{\mathrm{GeV}}\right) , \end{aligned}$$

(15)

where $x_{\mathrm{FO}}$ is given by Eq. (14). Let us note that in this case, production without cannibalism, the parameter $\lambda _s$ is small ($\lambda _s\lesssim 10^{-3}$) and plays no role in the WIMP DM phenomenology. In the next Subsection we will, however, consider the opposite case where large self-interactions do change the resulting DM abundance.

Figure 1 shows slices of the parameter space that give rise to the observed DM relic abundance. On the upper panel the cosmological parameters are fixed, $H_\text {EW}/m_h=10^{-16}$ (black lines) and $10^{-15}$ (blue lines), and $T_\text {end}/m_h=10^{-6}$ (dashed lines) and $10^{-4}$ (solid lines) while we scan over the relevant particle physics parameters ($\lambda _{hs}$ and $m_s$). The upper left corner in red, corresponding to $\lambda _{hs}>2\,m_s^2/v^2$, is excluded by the requirement discussed below Eq. (1). The figure shows that an increase in the dilution factor due to either an enhancement of the Hubble expansion rate $H_\text {EW}$ or a decrease in the temperature $T_\text {end}$ when the MD era ends has to be compensated with a higher DM abundance at the freeze-out. That, in turn, requires a smaller annihilation cross-section and hence a small $\lambda _{hs}$. The dependence on the DM mass $m_s$ is very mild.

The same conclusion can be extracted from the lower panel of Fig. 1, where the particle physics parameters are fixed, $m_s=20$ GeV (dashed lines) and 50 GeV (solid lines), and $\lambda _{hs}=10^{-3}$ (blue lines) and $10^{-2}$ (black lines) while we scan over the cosmological parameters. The left red band corresponds to a scenario which is not MD($H_\text {EW}/m_h<10^{-16}$), whereas the lower left corner corresponds to a case where the resulting SM temperature after the MD era ends is too small for successful BBN. Both cases are excluded from our analysis. Here the requirement of a non-relativistic freeze-out ($x_\text {FO}>3$) is also taken into account. Other observational constraints on the scenario will be discussed in Sect. 4.

Fig. 1

DM freeze-out without cannibalism. Parameter space giving rise to the observed DM relic abundance. The red regions correspond to the constraints discussed in Sect. 2.2. Other observational constraints are discussed in Sect. 4 and shown in Fig. 8

3.1.2 Freeze-out with cannibalism

The DM and visible sectors seize to be in chemical equilibrium with each other when $\langle \sigma _{ss\rightarrow bb}v \rangle n_s/H=1$. However, the s particles can maintain chemical equilibrium among themselves if number-changing interactions (namely, 4-to-2 annihilations with only DM particles both in the initial and final states, see Fig. 2) are still active. The condition for this so-called cannibalism is given by

$$\begin{aligned} \left. \frac{\langle \sigma _{ss\rightarrow b\bar{b}}v \rangle n_s}{\langle \sigma _{4\rightarrow 2}v^3 \rangle n_s^3}\right| _{x_\mathrm{FO}} \simeq \frac{\pi ^2}{81\sqrt{3}}\frac{\lambda _{hs}^2}{\lambda _s^4}\,x_\mathrm{FO}^3\,e^{2x_\mathrm{FO}} < 1, \end{aligned}$$

(16)

where we used

$$\begin{aligned} \langle \sigma _{4\rightarrow 2}v^3\rangle \simeq \frac{81\sqrt{3}}{32\pi }\,\frac{\lambda _s^4}{m_s^8} , \end{aligned}$$

(17)

in the non-relativistic approximation [20], and where $x_\mathrm{FO}$ is given by Eq. (14). In this case, the DM abundance is driven by the 4-to-2 annihilations and not anymore by the subdominant annihilations into SM particles. The Boltzmann equation governing the DM number density, Eq. (8), becomes

$$\begin{aligned} \frac{dn_{s}}{dt}+3\,H\,n_{s}=-\left\langle \sigma _{4\rightarrow 2} v^3\right\rangle \left[ n_{s}^{4}-n_s^2\left( n_{s}^\text {eq}\right) ^{2}\right] . \end{aligned}$$

(18)

Fig. 2

Examples of Feynman diagrams for the $4\rightarrow 2$ scalar self-annihilation process

If Eq. (16) was satisfied, the DM freeze-out is given by the decoupling of the 4-to-2 annihilations, defined by

$$\begin{aligned} \left. \frac{\langle \sigma _{4\rightarrow 2}v^3 \rangle n_s^3}{H}\right| _{T=T_\text {FO}^\text {c}} = 1, \end{aligned}$$

(19)

as can be inferred from Eq. (18). The time of freeze-out then is

$$\begin{aligned} x^{\mathrm{c}}_{\mathrm{FO}} \equiv \frac{m_s}{T_\text {FO}^\text {c}}= W\left[ 0.2\,\lambda _s^{4/3}\left( \frac{H_{\mathrm{EW}}}{m_h}\right) ^{-1/3}\left( \frac{m_s}{\mathrm{GeV}}\right) ^{-1/6} \right] ,\nonumber \\ \end{aligned}$$

(20)

where $W=W[\lambda _s,\,m_s,\,H_\text {EW}]$ is the 0-branch of the Lambert W function. The DM abundance then becomes (see again Appendix A)

$$\begin{aligned} \frac{\varOmega _{s}\,h^{2}}{0.12}&\simeq \,3\times 10^{-7}\, (x^{\mathrm{c}}_{\mathrm{FO}})^{3/2}\,e^{-x^{\mathrm{c}}_{\mathrm{FO}}}\nonumber \\&\quad \times \left( \frac{T_{\mathrm{end}}}{m_h}\right) ^{3/4}\left( \frac{H_{\mathrm{EW}}/m_h}{10^{-16}}\right) ^{-3/2}\,\left( \frac{m_s}{\mathrm{GeV}}\right) \,. \end{aligned}$$

(21)

When cannibalism is active, the 4-to-2 annihilations tend to increase the DM temperature with respect to the one of the SM bath [41]. However, we have checked that in all cases the DM and SM particles were still in kinetic equilibrium at the time of DM freeze-out, so that temperature of the s particle heat bath was the same as the SM temperature T. The condition for this is $\left. \langle \sigma _{sb\rightarrow sb}v \rangle n_b/H\right| _{x^\mathrm{c}_\mathrm{FO}} > 1$, where we have taken for simplicity $\langle \sigma _{sb\rightarrow sb}v \rangle \simeq \langle \sigma _{ss\rightarrow b\bar{b}}v \rangle$ and $n_b$ is the b-quark number density.

Fig. 3

DM freeze-out with cannibalism. Parameter space giving rise to the observed DM relic abundance, for $\lambda _{hs}=10^{-3}$. The red region corresponds to $\lambda _s>10$

Similar to Fig. 1, Fig. 3 also shows slices of the parameter space that give rise to the observed DM relic abundance. Here the cosmological parameters are fixed, $H_\text {EW}/m_h=10^{-16}$ (black lines) and $10^{-15}$ (blue lines), and $T_\text {end}/m_h=10^{-7}$ (dashed lines) and $10^{-5}$ (solid lines), while we scan over the particle physics parameters $\lambda _s$ and $m_s$ for a fixed $\lambda _{hs}=10^{-3}$. The upper band in red, corresponding to $\lambda _s>10$, is not considered. As in the previous case without cannibalism, an increase in the dilution factor has to be compensated with a higher DM abundance at the freeze-out. In this case with cannibalism, this requires a smaller annihilation cross-section and hence a small $\lambda _s$ or a heavier DM. The behavior with respect to $\lambda _{hs}$ and the cosmological parameters is very similar to the case without cannibalism (see Fig. 1) and is therefore not presented in this figure.

Fig. 4

DM freeze-out without (left column) and with (right column) cannibalism. Parameter space giving rise to the observed DM relic abundance. The red regions correspond to the constraints discussed in Sect. 2.2: the SM temperature after the matter-like component has decayed into SM particles must be larger than the BBN temperature and small enough not to not re-trigger DM production, Eq. (6); the DM freeze-out occurs while the s particles are non-relativistic, $x_\mathrm{FO}>3$; in a MD Universe $H_{\mathrm {EW}}/m_h >1.76\times 10^{-16}$; the portal coupling has to satisfy $\lambda _{hs}<2\,m_s^2/v^2$ and $\lambda _{hs}\ge \lambda _{hs}^\text {eq}$ with $\lambda _{hs}^\text {eq}$ given by Eq. (7). Other observational constraints are shown in Fig. 8

Before closing this subsection, we present the results of an extensive scan over the parameter space for the DM freeze-out without (left column) and with (right column) cannibalism in Fig. 4. The blue regions produce the observed DM relic abundance, whereas the red regions correspond to the constraints discussed in Sect. 2.2. The plots generalize the results of Figs. 1 and 3. First, let us note that the usual RD scenario can be recovered by taking $H_\text {EW}/m_h=H_\text {EW}^\text {rad}/m_h\simeq 1.76\times 10^{-16}$ and $T_\text {end}/m_h=1$. This corresponds to $\lambda _{hs}\simeq 10^{-1}$, in the case where DM mainly annihilates into b-quarks ($m_b<m_s\lesssim 50$ GeV) and does not undergo a cannibalism phase. In the MD scenario the Higgs portal coupling $\lambda _{hs}$ can reach much smaller values down to $\mathcal {O}(10^{-4})$. Such small values naturally need large dilution factors, characterized by large expansion rates $H_\text {EW}/m_h$ up to $\mathcal {O}(10^{-13})$ and low temperatures for the end of the MD era, $T_\text {end}/m_h$ down to $\mathcal {O}(10^{-8})$. In the case with cannibalism, $\lambda _{hs}\lesssim 10^{-2}$ while $\lambda _s\gtrsim 10^{-2}$ due to the fact that the DM annihilation into SM particles must decouple earlier than the 4-to-2 annihilations. Finally, we note that in the scenario where freeze-out occurs during a standard RD phase, cannibalism would generically require non-perturbative values of $\lambda _s$. As shown above, in the MD case the detailed effect of non-vanishing self-interactions can easily be taken into account, as the required values for $\lambda _s$ can be much smaller. This result, along with its observational consequences that we will present in Sect. 4, are among the most important novelties of this work.

3.2 The freeze-in case

In this subsection we assume, for simplicity, the mass hierarchy $m_s < m_h/2$, as we take the Higgs decay into two s to be the dominant production mechanism for DM. A more general analysis is again presented in Ref. [39] for the pure freeze-in case without cannibalism.

3.2.1 Freeze-in without cannibalism

The DM number density can again be computed using the Boltzmann equation (8), which in the absence of DM self-interactions is

$$\begin{aligned} \frac{dn_s}{dt}+3\,H\,n_s=2\,\frac{K_{1}(\frac{m_h}{T})}{K_{2}(\frac{m_h}{T})}\,\varGamma _{h\rightarrow ss}\,n_h^{\mathrm{eq}}, \end{aligned}$$

(22)

where $\varGamma _{h\rightarrow ss}$ is the partial decay width of the Higgs into two s-particles and $n_h^\text {eq}$ is its equilibrium number density. These quantities are given by

$$\begin{aligned} \varGamma _{h\rightarrow ss}= & {} \frac{\lambda _{hs}^2\,m_h}{64\pi \lambda _h}\sqrt{1-\left( \frac{2m_s}{m_h}\right) ^2} , \end{aligned}$$

(23)

$$\begin{aligned} n_h^\mathrm{eq}(T)= & {} \left( \frac{m_h\,T}{2\pi } \right) ^{3/2} e^{-\frac{m_h}{T}}\,. \end{aligned}$$

(24)

By then performing a change of variables, $\chi _s = n_s\,a^3$, where $\chi _s$ is the comoving s number density and a is the scale factor, we get the comoving DM number density at infinity³

$$\begin{aligned} \chi _s^\infty= & {} 2\, \varGamma _{h\rightarrow ss}\int _0^\infty \mathrm{dln}a\left( \frac{m_h\,T}{2\pi }\right) ^{3/2}\nonumber \\&\quad \times e^{-m_h/T}\frac{a^3}{H(a)}\frac{K_1\left( \frac{m_h}{T}\right) }{K_2\left( \frac{m_h}{T}\right) }\nonumber \\\simeq & {} 6.3\,\frac{\varGamma _{h\rightarrow ss}}{H_\text {EW}}\,n_h^\mathrm{eq}(m_h), \end{aligned}$$

(25)

where we have normalized the scale factor so that $a(T=m_h)\equiv a_\mathrm{EW}=1$. The numerical value of the above integral is not sensitive to the upper limit of integration, and we have set it for convenience to $a\rightarrow \infty$. As shown in the Appendix A, the DM abundance today can then be expressed as

$$\begin{aligned} \frac{\varOmega _{s}\,h^{2}}{0.12}\simeq&\, 2\times 10^{22}\,g_*(m_h)^{-1/4}\,\lambda _{hs}^2\nonumber \\&\times \left( \frac{H_{\mathrm {EW}}/m_h}{10^{-16}}\right) ^{-5/2}\,\left( \frac{T_\mathrm{end}}{m_h}\right) ^{3/4}\left( \frac{m_{s}}{\mathrm {GeV}}\right) , \end{aligned}$$

(26)

where we assumed $m_s\ll m_h/2$.

Let us emphasize that the result in Eq. (26) only applies to a scenario where the Universe was effectively MD during the DM yield, and therefore it is not, as such, applicable to other scenarios. To retain the usual RD case, one must set $T_\mathrm{end}=m_h$, use the result of Eq. (4) for $H_\mathrm{EW}$, and use the newly calculated prefactor 11.4 in Eq. (25) instead of 6.3 which we obtained above. These account for the facts that in our case not only there was entropy production at the end of the early MD phase but also that the expansion rate of the Universe at the time of DM freeze-in was different from that in the usual RD case.

Fig. 5

DM freeze-in without cannibalism. Parameter space giving rise to the observed DM relic abundance. The black dotted line shows the parameters yielding the correct DM abundance in the usual RD scenario. The red regions correspond to the constraints discussed in Sect. 2.2. Other observational constraints are shown in Fig. 8

Figure 5 shows slices of the parameter space that give rise to the observed DM relic abundance. On the upper panel the cosmological parameters are fixed, $H_\text {EW}/m_h=10^{-16}$ (black lines) and $10^{-13}$ (blue lines), and $T_\text {end}/m_h=10^{-5}$ (solid lines) and $10^{-1}$ (dashed lines), while we scan over the relevant particle physics parameters ($\lambda _{hs}$ and $m_s$). The upper left corner in red, corresponding to, $\lambda _{hs}>2\,m_s^2/v^2$, is excluded. The figure shows again that an increase in the dilution factor due to either an enhancement of the Hubble expansion rate $H_\text {EW}$ or a decrease in the temperature $T_\text {end}$ when the MD era ends has to be compensated with a higher DM abundance at the freeze-out. This requires an increase in either $m_s$ or the DM production via the Higgs decay (i.e. a bigger $\lambda _{hs}$). The thick dotted black line corresponds to the DM production in the usual RD scenario, characterized by $T_\text {end}/m_h=1$ and $H_\text {EW}/m_h=10^{-16}$. We note that, as expected, in the MD scenario the values for the required values for Higgs portal are always higher than in the RD case.

The same conclusion can be drawn from the lower panel of Fig. 5, where the particle physics parameters are fixed, $m_s=0.1$ GeV (solid lines) and 10 GeV (dashed lines), and $\lambda _{hs}=10^{-9}$ (blue lines) and $10^{-5}$ (black lines), while we scan over the cosmological parameters. The left band corresponds to a scenario which is not MD ($H_\text {EW}/m_h<10^{-16}$). The lower left and the upper right corners correspond to scenarios where the resulting SM temperature after the MD era ends is either too small for successful BBN or so large that it re-triggers the DM yield, respectively. All three cases are excluded from our analysis. Observational constraints on the scenario will be discussed in Sect. 4.

As in the case of freeze-out, the result of Eq. (26) is the final DM abundance only if number-changing DM self-interactions do not become active and the s particles do not reach chemical equilibrium with themselves. This is the scenario we will now turn into.

3.2.2 Freeze-in with cannibalism

Let us now calculate the final DM abundance following the thermalization and consequent cannibalism phase of the s particles. In this case, the Boltzmann equation (8) is

$$\begin{aligned} \frac{dn_s}{dt}+3\,H\,n_s=&\, 2\,\frac{K_{1}(\frac{m_h}{T})}{K_{2}(\frac{m_h}{T})}\,\varGamma _{h\rightarrow ss}\,n_h^{\mathrm{eq}}\nonumber \\&- \left\langle \sigma _{4\rightarrow 2} v^3\right\rangle \left[ n_{s}^{4}-n_s^2\left( n_{s}^\text {eq}\right) ^{2}\right] , \end{aligned}$$

(27)

where $\varGamma _{h\rightarrow ss}$ and $n_h^{\mathrm{eq}}$ are again given by Eqs. (23) and (24), respectively, and $\left\langle \sigma _{4\rightarrow 2} v^3\right\rangle$ by Eq. (17).

For values of the portal coupling required by non-thermalization of the hidden sector with the SM sector $\lambda _{hs}\lesssim (H_\mathrm{EW}/m_h)^{1/2}$, Eq. (7), the initial s particle number density in the hidden sector produced by Higgs decays is always smaller than the corresponding equilibrium number density. Thus, if the self-interactions are sufficiently strong (see below), the s particles can reach chemical equilibrium with themselves by first increasing their number density via 2-to-4 annihilations, and then undergo cannibalism when they become non-relativistic, as discussed in e.g. Refs. [46, 48, 54]. A possible caveat to this is the case where $m_s$ is close to $m_h$, as then the eventual dark freeze-out would occur before the yield from the SM sector has ended. In that case, the production mechanism is dubbed as reannihilation [74, 75]. Because in that case the s particles would not, in general, be in thermal equilibrium at the time of their freeze-out, finding the correct DM abundance requires solving the Boltzmann equation for the DM distribution function instead of number density, which is beyond the scope of this work. In this paper we therefore choose an approach where we solve the Boltzmann equation for DM number density but highlight the regime in our results where reannihilations could potentially alter our conclusions, and leave solving the Boltzmann equation for DM distribution function for future work. Because the freeze-in yield has ended by $T\sim 0.1 m_h$ [75], we take this regime to be determined by $m_s\gtrsim 10$ GeV. As we will show, this is only a small part of the observationally interesting parameter space, especially for DM self-interactions.

In the following, we will solve Eq. (27) in the limit where the self-interactions of s are large, to complement the usual freeze-in scenario discussed above. Note that the 2-to-2 scalar self-annihilations do not have a net effect on the final DM abundance and are therefore not included in Eq. (27).

The number-changing s self-interactions in Eq. (27) become active if

$$\begin{aligned} \left. \frac{\langle \sigma _{4\rightarrow 2}v^3\rangle \left( n_s^\mathrm{init}\right) ^3}{H}\right| _{a_\mathrm{nrel}} > 1 , \end{aligned}$$

(28)

where $n_s^\mathrm{init}(a_\mathrm{nrel}) = \chi _s^\infty (a_\text {EW}/a_\mathrm{nrel})^3$ is the initial s particle abundance produced by Higgs decays, where $\chi _s^\infty$ is given by Eq. (25), and we have invoked the principle of detailed balance. The scale factor $a_\text {nrel}$ when the s particles become non-relativistic can be solved from

$$\begin{aligned} \frac{p_s}{m_s} \simeq \frac{m_h}{2m_s}\frac{a_\mathrm{EW}}{a_\mathrm{nrel}} \simeq 1 , \end{aligned}$$

(29)

so that $a_\mathrm{nrel} \simeq m_h/(2m_s)$ (recall that $a_\mathrm{EW}=1$). Here we assumed $m_s\ll m_h/2$, so that the initial s particle momenta are $p\simeq m_h/2$. As discussed in Refs. [48, 76], it indeed suffices to evaluate Eq. (28) at $a_\mathrm{nrel}$, which is the latest moment when the s particles can reach chemical equilibrium with themselves.

Reminiscent to the standard WIMP case, the final DM abundance only depends on the time of the freeze-out, and therefore the scenario is not sensitive to when the hidden sector thermalization occurs. Thus, the thermalization condition for the s field’s quartic self-interaction strength can be solved from Eq. (28) to be

$$\begin{aligned} \lambda ^{\mathrm{FI}}_s \simeq 6.6\,\lambda _{hs}^{-3/2}\left( \frac{m_s}{\mathrm{GeV}}\right) ^{1/8}\frac{H_{\mathrm{EW}}}{m_h}\,. \end{aligned}$$

(30)

If $\lambda _s<\lambda ^\mathrm{FI}_s$, the final yield is given by Eq. (26); if not, cannibalism has to be taken into account in solving Eq. (27). Therefore, if $\lambda _s>\lambda ^\mathrm{FI}_s$, the s particles thermalize with themselves and the sector exhibits a cannibal phase before the final freeze-out of DM density from the hidden sector heat bath. The time of the dark freeze-out of s particles can be solved in the standard way from Eq. (27) as the time when the 4-to-2 interaction rate equals the Hubble expansion rate

$$\begin{aligned} \left. \frac{\langle \sigma _{4\rightarrow 2}v^3\rangle n_s^3}{H}\right| _{T_s^\text {FO}} = 1, \end{aligned}$$

(31)

where H is given by Eq. (13) and

$$\begin{aligned} n_s(T_s) = \left( \frac{m_sT_s}{2\pi }\right) ^{\frac{3}{2}}e^{-\frac{m_s}{T_s}} = \frac{m_s^3}{(2\pi )^{3/2}}x_s^{-3/2}e^{-x_s} , \end{aligned}$$

(32)

where $T_s$ is the temperature of the hidden sector heat bath which in general is not the same as the SM sector temperature, $T_s\ne T$. Here we also introduced the conventional units $x_s\equiv m_s/T_s$.

The relation between $T_s$ and T can be inferred from entropy conservation, as after the thermalization within the hidden sector the two entropy densities are separately conserved. First, consider the times when the s particles are still relativistic, whence

$$\begin{aligned} \zeta&\equiv \left. \frac{\mathfrak {s}_\mathrm{rad}}{\mathfrak {s}_\mathrm{hid}}\right| _\mathrm{rel} = \frac{g_{*\mathfrak {s}}\,T^3}{T_s^3}= g_{*\mathfrak {s}}\left( \frac{\rho _\mathrm{SM}}{g_*\,\rho _s}\right) ^{3/4}\nonumber \\&= g_{*\mathfrak {s}}\left( \frac{\rho _\mathrm{SM}}{g_*(m_h/2)\,n_s^\mathrm{init}}\right) ^{3/4}, \end{aligned}$$

(33)

where $\mathfrak {s}_\mathrm{rad}$ and $\mathfrak {s}_\mathrm{hid}$ are the SM and hidden sector entropy densities, respectively, and $g_{*\mathfrak {s}}$ corresponds to the relativistic degrees of freedom that contribute to the SM entropy density. On the other hand, between the moment when the s particles became non-relativistic and their final freeze-out, the ratio $\zeta$ is

$$\begin{aligned} \zeta = \left. \frac{\mathfrak {s}_\mathrm{rad}}{\mathfrak {s}_\mathrm{hid}}\right| _\mathrm{nrel} = \frac{2\pi ^2(2\pi )^{3/2}\,g_*(T)}{45}\,\frac{T^3}{m_s^3}\,x_s^{1/2}\,e^{x_s}, \end{aligned}$$

(34)

where we used $\mathfrak {s}_\mathrm{hid} = m_s\,n_s(T_s)/T_s$. By equating Eqs. (33) and (34), one can express the SM sector temperature T as a function of the hidden sector temperature

$$\begin{aligned} T\simeq 1.7\,\lambda _{hs}^{-1/2}\left( \frac{H_\text {EW}}{m_h}\right) ^{1/4}x_s^{-1/6}\,e^{-x_s/3}\,m_s. \end{aligned}$$

(35)

The moment of the dark freeze-out can then be calculated be using Eqs. (31), (32), (13) and (35), which give

$$\begin{aligned} x_s^\mathrm{FO} {=} \frac{17}{10} W\left[ 0.1\,\lambda _s^{16/17}\lambda _{hs}^{3/17}\left( \frac{m_h}{m_s}\right) ^{2/17}\left( \frac{H_\text {EW}}{m_h}\right) ^{-11/34} \right] ,\nonumber \\ \end{aligned}$$

(36)

where $W=W[\lambda _s,\,\lambda _{hs},\,m_s,\,H_\text {EW}]$ is again the 0-branch of the Lambert W function. The final DM abundance after the freeze-out then is

$$\begin{aligned} n_s^\mathrm{final} = \frac{m_s^3}{(2\pi )^{3/2}}(x_s^\mathrm{FO})^{-3/2}e^{-x_s^\mathrm{FO}} , \end{aligned}$$

(37)

from which the DM abundance today can be calculated to be

$$\begin{aligned} \frac{\varOmega _{s}\,h^2}{0.12}&\simeq \,3\times 10^{8}\,g_{*}(T_\mathrm{FO})^{-1/4}\left( \frac{n_{s}^{\mathrm {final}}}{T_\mathrm{FO}^{3}}\right) \nonumber \\&\quad \times \,\left( \frac{H_\mathrm{EW}/m_h}{10^{-16}}\right) ^{-3/2}\,\left( \frac{T_{\mathrm {end}}}{m_h}\right) ^{3/4}\,\left( \frac{m_{s}}{\mathrm {GeV}}\right) , \end{aligned}$$

(38)

as shown in the Appendix A. Using then Eqs. (35) and (37), we get a relation for $x_{s}^{\mathrm {FO}}$ that takes into account the present DM abundance

$$\begin{aligned} x_{s}^{\mathrm {FO}}&\simeq \, 4\times 10^{18}\,g_*^{-1/4}\lambda _{hs}^{3/2}\left( \frac{\varOmega _s h^2}{0.12}\right) ^{-1}\nonumber \\&\quad \times \left( \frac{H_\mathrm{EW}/m_h}{10^{-16}}\right) ^{-9/4}\left( \frac{T_\mathrm{end}}{m_h}\right) ^{3/4}\left( \frac{m_{s}}{\mathrm {GeV}}\right) . \end{aligned}$$

(39)

Equating this result with Eq. (36) then gives the connection between the model parameters $m_{s}$, $\lambda _{s}$, $\lambda _{hs}$, $T_{\mathrm {end}}$, $H_{\mathrm {EW}}$ that yields the correct DM abundance.

Figure 6 shows again slices of the parameter space that give rise to the observed DM relic abundance. The cosmological parameters are fixed and we scan over the particle physics parameters, fixing $\lambda _{hs}=10^{-9}$ in the upper panel and $m_s=1$ GeV in the lower panel. The red bands, corresponding to $\lambda _s>10$ (perturbativity bound) and $\lambda _{hs}\gtrsim 3\times 10^{-5}$ ($\lambda _{hs} < 2m_s^2/v^2$ in order to avoid a spontaneous symmetry breaking in the s direction) are excluded. Again, an increase in the dilution factor due to either an enhancement of the Hubble expansion rate $H_\text {EW}$ or a decrease in the temperature $T_\text {end}$ when the MD era ends has to be compensated with a higher DM abundance at the dark freeze-out. This requires a smaller 4-to-2 annihilation cross-section and hence a small $\lambda _s$.

Fig. 6

DM freeze-in with cannibalism. Parameter space giving rise to the observed DM relic abundance, for $\lambda _{hs}=10^{-9}$ (upper panel) and $m_s=1$ GeV (lower panel). The red regions correspond to the constraints discussed in Sect. 2.2, and the shaded region in the upper panel to the reannihilation regime. Other observational constraints are shown in Fig. 8

Figure 7 depicts the results of an extensive scan over the parameter space for the DM freeze-in without (left column) and with (right column) cannibalism. The blue regions produce the observed DM relic abundance, the red regions correspond to the constraints discussed in Sect. 2.2. Other observational constraints on the scenario will be discussed in Sect. 4.

The plots generalize the results of Figs. 5 and 6. First, the usual RD scenario without cannibalism can be approximately recovered by taking $H_\text {EW}/m_h=H_\text {EW}^\text {rad}/m_h\simeq 1.76\times 10^{-16}$ and $T_\text {end}/m_h=1$, as discussed in Sect. 3.2.1. This corresponds to the black dotted line with $\lambda _{hs}\simeq \mathcal {O}(10^{-11})$. Second, in the MD scenario the Higgs portal can reach much higher values up to $\mathcal {O}(10^{-4})$. Such big values for freeze-in naturally need large dilution factors, characterized by large expansion rates $H_\text {EW}/m_h$ up to $\mathcal {O}(10^{-11})$ and low temperatures for the end of the MD era, $T_\text {end}/m_h$ down to $\mathcal {O}(10^{-8})$. Higher values of $\lambda _{hs}$ cannot be reached, because in the present case thermalization with the SM must be avoided.

Fig. 7

DM freeze-in without (left column) and with (right column) cannibalism. Parameter space giving rise to the observed DM relic abundance. The black dotted line shows the parameters yielding the correct DM abundance in the usual RD scenario. The red regions correspond to the constraints discussed in Sect. 2.2: the SM temperature after the matter-like component has decayed into SM particles must be larger than the BBN temperature and small enough not to not re-trigger DM production, Eq. (6); in a MD Universe $H_{\mathrm {EW}}/m_h >1.76\times 10^{-16}$; the portal coupling has to satisfy $\lambda _{hs}<2\,m_s^2/v^2$ and $\lambda _{hs}<\lambda _{hs}^\text {eq}$ with $\lambda _{hs}^\text {eq}$ given by Eq. (7). The shaded region in panels on the right hand side corresponds to the reannihilation regime. Other observational constraints are shown in Fig. 8

4 Observational properties

Finally, we turn into observational prospects, discussing collider signatures, direct and indirect detection, as well as the observational consequences of DM self-interactions.

4.1 Collider signatures

For small singlet masses, $m_s < m_h/2$, the Higgs can decay efficiently into a pair of DM particles. Thus, the current limits on the invisible Higgs branching ratio (BR$_\text {inv}\lesssim 20\%$ [77]) and the total Higgs decay width ($\varGamma ^\text {tot}\lesssim 22$ MeV [78]) constrain the Higgs portal coupling, $\lambda _{hs}$, by Eq. (23). This constraint applies to both freeze-out and freeze-in scenarios, although typically it can be expected to constrain only the freeze-out case, as usually in freeze-in scenarios the value of $\lambda _{hs}$ required to reproduce the observed DM abundance is orders of magnitudes below these values. Indeed, the collider signatures of frozen-in DM were recently deemed unobservable in Ref. [79]. However, the paper considered only the usual RD case, and in a scenario containing an early phase of rapid expansion, such as in the present paper, the portal coupling can take a much larger value than what is usually encountered in the context of freeze-in. It is therefore not a priori clear whether constraints of the above kind can be neglected or not. We will present them in Sect. 4.4.

Fig. 8

Detection prospects for frozen-out and frozen-in DM with and without cannibalism, as indicated in the figures. The green regions are excluded by different measurements: DM direct detection, invisible Higgs decay, or DM self-interactions. The blue regions give rise to the observed DM relic abundance, the light blue being already in tension with observations. The black thick dashed line corresponds to the bounds that might be reached by next generation direct detection DM experiments. The red regions correspond to the constraints discussed in Sect. 2.2: the SM temperature after the matter-like component has decayed into SM particles must be larger than the BBN temperature and small enough not to not re-trigger DM production, Eq. (6); the DM freeze-out occurs while the s particles are non-relativistic, $x_\mathrm{FO}>3$; in a MD Universe $H_{\mathrm {EW}}/m_h >1.76\times 10^{-16}$; the portal coupling has to satisfy $\lambda _{hs}<2\,m_s^2/v^2$ and $\lambda _{hs}\ge \lambda _{hs}^\text {eq}$ for the freeze-out case and $\lambda _{hs}< \lambda _{hs}^\text {eq}$ for the freeze-in case, with $\lambda _{hs}^\text {eq}$ given by Eq. (7). The shaded region in the lower panels corresponds to the reannihilation regime. The black dotted line shows the parameters yielding the correct DM abundance in the usual RD scenario

In MD cosmologies, the interaction rates required to produce the observed DM abundance via freeze-in could lead to displaced signals at the LHC and future colliders [25]. However, as in our scenario DM is produced via the decay of the Higgs, we will have no exotic signals displaced from the primary vertex.

4.2 Direct and indirect detection signatures

The direct detection constraint is obtained by comparing the spin-independent cross section for the scattering of the DM off of a nucleon,

$$\begin{aligned} \sigma _\text {SI}=\frac{\lambda _{hs}^2\,m_N^4\,f^2}{4\pi \,m_s^2\,m_h^4}, \end{aligned}$$

(40)

to the latest limits on $\sigma _\text {SI}$ provided by PandaX-II [80], LUX [81] and Xenon1T [82]. Here $m_N$ is the nucleon mass and $f\simeq 1/3$ corresponds to the form factor [83,84,85,86]. We also take into account the projected sensitivities of the next generation DM direct detection experiments like LZ [87] and DARWIN [88]. Moreover, multiple experimental setups have recently been suggested for the detection of elastic scatterings of DM in the mass range from keV to MeV [89,90,91,92,93,94,95,96,97,98,99,100,101,102,103]. In particular, the typical DM-electron cross sections for MeV-scale FIMP DM could be tested by some next generation experiments [104,105,106,107].

The current limits from the analysis of gamma-rays coming from dwarf spheroidal galaxies with Fermi-LAT and DES [108,109,110] do not probe relevant parts of our parameter space. In the case of freeze-in, indirect detection signals can be expected in scenarios where the singlet scalar is a mediator and the hidden sector exhibits a richer structure, as recently studied in Ref. [111].

4.3 Dark matter self-interactions

Finally, we consider the observational ramifications of DM self-interactions. Two long-standing puzzles of the collisionless cold DM paradigm are the ‘cusp vs. core’ [112,113,114,115,116,117] and the ‘too-big-to-fail’ [118, 119] problems. These issues are collectively referred to as small scale structure problems of the $\Lambda$CDM model; for a recent review, see Ref. [120]. These tensions can be alleviated if at the scale of dwarf galaxies DM exhibits a large self-scattering cross section, $\sigma$, over DM particle mass, $m_s$, in the range $0.1\lesssim \sigma /m_s\lesssim 10$ cm$^2$/g [121,122,123,124,125,126,127,128,129,130]. Nevertheless, the non observation of an offset between the mass distribution of DM and galaxies in the Bullet Cluster constrains such self-interacting cross section, concretely $\sigma /m_s < 1.25$ cm$^2$/g at $68\%$ CL [131,132,133]. In the limit $m_s\ll m_h$ we have

$$\begin{aligned} \frac{\sigma }{m_{s}}\simeq \frac{9}{32\pi }\,\frac{\lambda _{s}^{2}}{m_{s}^{3}}\lesssim 1.25\,\mathrm {\frac{cm^{2}}{g}}, \end{aligned}$$

(41)

which imposes an important constraint

$$\begin{aligned} \lambda _{s}\lesssim 2\times 10^{2}\,\left( \frac{m_{s}}{\mathrm {GeV}}\right) ^{3/2}, \end{aligned}$$

(42)

which we will show in our results in the next Subsection.

In the present case, no cosmological signatures can be expected. Even though in the case where the singlet scalar never thermalizes with the SM sector the DM generically comprises an isocurvature mode in the CMB fluctuations [72, 73], the relative amount of such perturbations gets strongly diluted due to $\rho _s \ll \rho _\mathrm{tot}$, leaving no observable imprints on the CMB.

4.4 Results

Figure 8 depicts the detection prospects for frozen-out and frozen-in DM, with and without cannibalism. The green regions are excluded by different observations discussed in the above subsections: DM direct detection, invisible Higgs decay or DM self-interactions. The blue regions give rise to the observed DM relic abundance, the light blue region being already in tension with observations. The black thick dashed line corresponds to the bounds that might be reached by next generation direct detection DM experiments. The constraints discussed in Sect. 2.2 are shown in red. Finally, the black dotted line shows the parameters yielding the correct DM abundance in the usual RD scenario.

In the MD scenario, DM direct detection already excludes an important region of the parameter space for the freeze-out case both with and without cannibalism. More interestingly, the next generation of DM direct detection experiments will be able to probe almost the whole region of parameter space compatible with the DM relic abundance, for the freeze-out scenario with $m_s<m_h/2$.

On the other hand, the regions favored by freeze-in could be tangentially probed by next generation of direct detection experiments. A particularly interesting thing in this case is that the effect of non-vanishing self-interactions seems to be crucial in determining whether the scenario can be tested by the next-generation direct detection experiments or not, as shown in the two lower panels of Fig. 8. In the standard RD case the freeze-in scenario obviously does not have any such observational consequences, as the required values of $\lambda _{hs}$ are in that case much smaller regardless of the value of $\lambda _s$. Note that in Ref. [39] we obtained the opposite result, showing that FIMP DM cannot be tested by the next-generation direct detection experiments. This conclusion, however, is due to different assumptions for the decay of the matter-like component, as discussed in Sect. 1. However, observational constraints on DM self-interactions already rule out a corner of the parameter space corresponding to MeV-scale masses regardless of the prospects for direct detection. Finally, the region between the two dashed lines in the case of freeze-in with cannibalism corresponds to 0.1 cm$^2$/g$~<\sigma /m_s<10$ cm$^2$/g, the zone where the small-scale structure tensions can be alleviated.

5 Conclusions

In cosmology, one typically assumes that at early times the Universe was radiation-dominated from the end of inflation. However, there are no indispensable reasons to assume that, and alternative cosmologies not only can lead to interesting observational ramifications but are also well-motivated. For example, an early period of matter domination is still a perfectly viable option.

In this context, we studied different dark matter production mechanisms during an early MD era. We focused first on the usual case where DM is produced by the freeze-out mechanism, corresponding to the WIMP paradigm. Then, the assumption of thermal equilibrium with the SM was relaxed allowing the DM to be produced via the freeze-in mechanism, corresponding to FIMP DM. For these two cases, we took for the first time into account the effects of sizable self-interactions within the hidden sector. Indeed, as we showed in the present context, DM self-interactions can be crucial for the determination of the final DM relic abundance and observational consequences.

When the expansion rate of the Universe differs from the usual radiation-dominated case, it tends to effectively dilute the DM abundance when the era of non-standard expansion ends and the visible sector gets reheated. This means that in case the expansion was faster than in the RD case and the DM particles were initially in thermal equilibrium with the visible sector, they generically have to undergo freeze-out earlier than in the usual RD case, thus resulting in larger DM abundance to match the observed one. In case the DM particles interacted so feebly that they never became part of the SM equilibrium heat bath, the coupling between DM and the visible sector typically has to be orders of magnitude larger than in the usual freeze-in case to compensate the larger expansion rate. As we showed, sizable self-interactions can further complicate this picture. Production of self-interacting DM during a non-standard expansion phase may thus result in important experimental and observational ramifications, as shown in Fig. 8.

In this paper we studied a benchmark scenario where the SM is extended with a real singlet scalar DM, odd under a $\mathbb {Z}_2$ symmetry. It would be interesting to see what are the consequences in other models where, for example, the hidden sector has a richer structure (e.g. sterile neutrinos, gauge structure, etc.) or where the DM is not coupled to the SM via the Higgs portal but via some other portal, for example the $Z'$ or a lepton portal [134,135,136].

Appendix A: Dark matter abundance in the present universe

The DM abundance at present is

$$\begin{aligned} \varOmega _{s}h^2= \frac{\rho _s}{\rho _\mathrm{c}/h^2} = \frac{\xi \, \mathfrak {s}_0}{\rho _\mathrm{c}/h^2} , \end{aligned}$$

(43)

where $\mathfrak {s}_0=2891$ cm$^{-3}$ and $\rho _\mathrm{c}/h^2 = 1.054\times 10^{-5}$ GeV/cm$^3$ are, respectively, the entropy density and critical energy density today [25], and

$$\begin{aligned} \xi \equiv \frac{\rho _s(T_\mathrm{end}')}{\mathfrak {s}(T_\mathrm{end}')} = m_s\frac{n_s(T_\mathrm{end}')}{\mathfrak {s}(T_\mathrm{end}')} = m_s\frac{\chi _s^\infty }{\mathcal {S}(T_\mathrm{end}')}, \end{aligned}$$

(44)

where $\chi _s^\infty \equiv a^3 n_s$ is the comoving DM number density after freeze-in/-out and the SM entropy at the temperature the SM sector gained when the MD ended, $T_\mathrm{end}'$, is given by

$$\begin{aligned} \mathcal {S}(T_\mathrm{end}') = \frac{2\pi ^2}{45}g_{*\mathfrak {s}}(T_\mathrm{end}')\,T_\mathrm{end}'^3\, a_\mathrm{end}^3. \end{aligned}$$

(45)

Only after this point the comoving entropy density in the SM sector is conserved. Note that from this point on, the expansion history of the Universe does not affect the result. In Eq. (45), $a_\mathrm{end}$ can be replaced by $H_\mathrm{EW}$ by using the Friedmann equation, $H_\mathrm{end} \propto H_\mathrm{EW}\,a_\mathrm{end}^{-3/2} \propto T_\mathrm{end}'^2/M_\mathrm{P}$, so that

$$\begin{aligned} a_\mathrm{end}^3 = \left( \frac{90}{\pi ^2 g_*(T_\mathrm{end}')}\right) \left( \frac{M_\mathrm{P}\,H_\mathrm{EW}}{T_\mathrm{end}'^2}\right) ^2. \end{aligned}$$

(46)

We reiterate that we have normalized the scale factor so that $a_\mathrm{EW}=1$.

One can then either substitute the comoving number density $\chi _s^\infty$ into Eq. (44) (as in the case of Eq. (25), which gives the result (26)) or calculate the actual DM number density $n_s(T_\mathrm{end}')$ in Eq. (44) by relating it to the number density at the time the DM production ended

$$\begin{aligned} n_{s}\left( T_\mathrm{end}'\right) =n_{s}^{\mathrm {final}}\left( T_{\mathrm {F}}\right) \frac{g_{*\mathfrak {s}}(T_\mathrm{end})}{g_{*\mathfrak {s}}(T_\mathrm{F})}\,\left( \frac{T_{\mathrm {end}}}{T_\mathrm{F}}\right) ^{3}, \end{aligned}$$

(47)

as in the case of Eqs. (15), (21) and (38). Relating $n_{s}\left( T_\mathrm{end}'\right)$ to $T_\mathrm{end}$ but using $T_\mathrm{end}'$ for the entropy density $\mathfrak {s}$ in Eq. (44) leads to an artificial discontinuity in DM number density. This reflects the fact that we assume that the dominant matter-like component decays instantaneously to the SM sector, heating the SM particles instantaneously from temperature $T_\mathrm{end}$ to a higher temperature $T_\mathrm{end}'$ and simultaneously effectively diluting the DM number density.

The relation between $T_{\mathrm {end}}$ and $T_{\mathrm {end}}'$ can be found as follows. Following Eq. (13), the matter-like component’s energy density can be written as

$$\begin{aligned} \rho _{\mathrm {M}}\left( T\right) =3\,M_{\mathrm {P}}^{2}\,H_{\mathrm {EW}}^{2}\,\left( \frac{T}{m_h}\right) ^{3}\,\left( \frac{g_{*}\left( T\right) }{g_{*}\left( m_h\right) }\right) , \end{aligned}$$

(48)

and the SM energy density in the usual way as

$$\begin{aligned} \rho _{\mathrm {SM}}\left( T\right) =\frac{\pi ^{2}}{30}\,g_{*}\left( T\right) \,T^{4}. \end{aligned}$$

(49)

At $T=T_{\mathrm {end}}$, the matter-like component transfers all of its energy into the SM sector, $\rho _{\mathrm {M}}\left( T_{\mathrm {end}}\right) =\rho _{\mathrm {SM}}\left( T_{\mathrm {end}}'\right)$, so that one finds

$$\begin{aligned} \frac{T_{\mathrm {end}}}{T_{\mathrm {end}}'}&\simeq 0.4\left( \frac{H_\mathrm{EW}/m_h}{10^{-16}}\right) ^{-1/2}\,\left( \frac{T_{\mathrm {end}}}{m_h}\right) ^{1/4}\nonumber \\&\quad \times \left( \frac{g_{*}\left( T_{\mathrm {end}}'\right) }{g_{*}\left( T_{\mathrm {end}}\right) }\,g_{*}\left( m_h\right) \right) ^{1/4}. \end{aligned}$$

(50)

Substituting this results into Eq. (47) and the resulting expression into Eq. (44) then gives the present DM abundance as a function of the model parameters. This procedure gives us the results (21) and (38).

The relation (50) also makes it possible to constraint the duration of the early MD phase. As discussed in Sect. 2.2, we require that the SM temperature after the matter-like component has decayed into SM particles, $T_\mathrm{end}'$, must be larger than the BBN temperature $T_\mathrm{BBN}=4$ MeV, and also that the temperature has to be smaller than either the final freeze-out temperature or smaller than $m_h$ in the freeze-in case in order not to re-trigger the DM yield after the decay of the matter-like component. This is what gives the conditions in Eq. (6). In order to determine the numerical values, we use $g_*(T_\mathrm{end}')=106.75$ for the upper limit and $g_*(T_\mathrm{end}')=10.75$ for the lower limit.