Annihilation Signatures of Hidden Sector Dark Matter Within Early-Forming Microhalos

TL;DR

This paper investigates the formation and annihilation signatures of microhalos from hidden sector dark matter, revealing that their high density leads to detectable gamma-ray signals, which constrain many such models.

Contribution

It provides a detailed calculation of density perturbation evolution and nonlinear collapse of microhalos in hidden sector models, using N-body simulations to predict their properties and observational signatures.

Findings

Microhalos are extremely dense, enhancing dark matter annihilation signals.

Fermi gamma-ray data constrains a wide range of hidden sector dark matter models.

Longer early matter-dominated eras can suppress microhalo formation, affecting detectability.

Abstract

If the dark matter is part of a hidden sector with only very feeble couplings to the Standard Model, the lightest particle in the hidden sector will generically be long-lived and could come to dominate the energy density of the universe prior to the onset of nucleosynthesis. During this early matter-dominated era, density perturbations will grow more quickly than otherwise predicted, leading to a large abundance of sub-earth-mass dark matter microhalos. Since the dark matter does not couple directly to the Standard Model, the minimum halo mass is much smaller than expected for weakly interacting dark matter, and the smallest halos could form during the radiation-dominated era. In this paper, we calculate the evolution of density perturbations within the context of such hidden sector models and use a series of $N$-body simulations to determine the outcome of nonlinear collapse during radiation domination. The resulting microhalos are extremely dense, which leads to very high rates of dark matter annihilation and to large indirect detection signals that resemble those ordinarily predicted for decaying dark matter. We find that the Fermi Collaboration's measurement of the high-latitude gamma-ray background rules out a wide range of parameter space within this class of models. The scenarios that are most difficult to constrain are those that feature a very long early matter-dominated era; if microhalos form prior to the decay of the unstable hidden sector matter, the destruction of these microhalos effectively heats the dark matter, suppressing the later formation of microhalos.

Equations

$$\mathcal{L}\supset -\frac{ϵ}{2}{B}^{\mu \nu }{Z}_{\mu \nu }^{\prime }+g_X{Z}_{\mu }^{\prime }\bar{X}{\gamma }^{\mu }X.$$

$$g_{f_{R,L}}=\epsilon\bigg{(}\frac{m^{2}_{Z^{\prime}}g_{Y}Y_{f_{R,L}}-m^{2}_{Z}g\sin\theta_{W}\cos\theta_{W}Q_{f}}{m^{2}_{Z}-m^{2}_{Z^{\prime}}}\bigg{)},$$

$$\Gamma_{Z^{\prime}}=\sum_{f}\frac{n_{c}m_{Z^{\prime}}\beta_{f}}{12\pi}\bigg{[}g^{2}_{fv}\bigg{(}1+\frac{2m_{f}^{2}}{m^{2}_{Z^{\prime}}}\bigg{)}+g^{2}_{fa}\bigg{(}1-\frac{4m_{f}^{2}}{m^{2}_{Z^{\prime}}}\bigg{)}\bigg{]},$$

$$\tau_{Z^{\prime}}\approx 4\times 10^{-6}\,{\rm s}\,\times\bigg{(}\frac{10^{-10}}{\epsilon}\bigg{)}^{2}\bigg{(}\frac{{\rm TeV}}{m_{Z^{\prime}}}\bigg{)}.$$

$$\frac{1}{2}\sigma v(XX\rightarrow Z^{\prime}Z^{\prime})\approx\frac{\pi\alpha^{2}_{X}}{m^{2}_{X}}\bigg{(}1-\frac{7v^{2}}{12}\bigg{)},$$

$$\xi=\xi_{\rm inf}\bigg{(}\frac{g^{h}_{\star,{\rm inf}}}{g^{h}_{\star}}\bigg{)}^{1/3}\bigg{(}\frac{g_{\star}}{g_{\star,{\rm inf}}}\bigg{)}^{1/3},$$

$${\Omega }_X{h}^2$$

$${\rho }_{Z^{\prime }}=\frac{13\zeta (3)g_{\star }{m}_{Z^{\prime }}{\xi }_{\inf }^3{T}^3}{2{\pi }^2{g}_{\star ,\inf }},$$

$$T_{\mathrm{dom}}$$

$${\Gamma }_{Z^{\prime }}\equiv \sqrt{\frac{8{\pi }^3{g}_{\star }(T_{\mathrm{RH}})}{90}}\frac{T_{\mathrm{RH}}^2}{m_{\mathrm{Pl}}},$$

$$T_{\rm RH}\approx 620\,{\rm MeV}\bigg{(}\frac{\epsilon}{10^{-10}}\bigg{)}\bigg{(}\frac{m_{Z^{\prime}}}{{\rm TeV}}\bigg{)}^{1/2}\bigg{(}\frac{10}{g_{\star}}\bigg{)}^{1/4},$$

$$\frac{T_{\rm dom}}{T_{\rm RH}}\simeq 36\,\times\xi_{\rm inf}^{3}\,\bigg{(}\frac{10^{-10}}{\epsilon}\bigg{)}\bigg{(}\frac{m_{Z^{\prime}}}{{\rm TeV}}\bigg{)}^{1/2}\bigg{(}\frac{g_{\star}}{10}\bigg{)}^{1/4}.$$

$${\Omega }_X{h}^2$$

$$\frac{1}{R}\frac{d^{2}R}{d{t}^2}=-\frac{4\pi G}{3}\left({\rho }_X+{\rho }_{Z^{\prime }}+2{\bar{\rho }}_R-2{\bar{\rho }}_{\Lambda }\right),$$

$${\rho }_X+{\rho }_{Z^{\prime }}=\left[{\bar{\rho }}_{X,i}+\left(\frac{{\bar{\rho }}_{Z^{\prime }}{a}^3}{{\bar{\rho }}_{Z^{\prime },i}{a}_i^3}\right){\bar{\rho }}_{Z^{\prime },i}\right](1+{\delta }_i){\left(\frac{R_i}{R}\right)}^3,$$

$$\frac{dR}{dt}=H_i{R}_i\left(1-\frac{{\delta }_i}{3}\right).$$

$$1+{\delta }_{\mathrm{SC}}=(1+{\delta }_i){\left(\frac{a}{a_i}\right)}^3{\left(\frac{R_i}{R}\right)}^3.$$

$$m=\left(\frac{{\bar{\rho }}_{Z^{\prime }}+{\bar{\rho }}_X}{{\bar{\rho }}_X}\right)m_0,$$

$${\mathcal{P}}_f(k)=\mathcal{P}(k)\exp (-R_f^2{k}^2),$$

$$\rho (r)=\left[A\ln \left(\frac{1+z_c}{1+z}\right)+B\right]{\bar{\rho }}_M\frac{R_f}{r},$$

$$\frac{df}{d\ln M}=\sqrt{\frac{2}{\pi }}\left|\frac{d\ln \sigma }{d\ln M}\right|\frac{{\delta }_c}{\sigma (M,z)}\exp \left[-\frac{{\delta }_c^2}{2{\sigma }^2(M,z)}\right],$$

$$\frac{t_{\Delta M}}{t_{\mathrm{dyn}}}\simeq$$

$$\simeq$$

$$\times \sqrt{\frac{T_{\mathrm{eq}}}{T_{\mathrm{RH}}}+e^{-{\Gamma }_{Z^{\prime }}\Delta t}},$$

$$\frac{a_f}{a_{\mathrm{RH}}}\text{raise}1.29167pt>\text{lower}4.30554pt\sim {\left(\frac{81}{32\pi }\frac{T_{\mathrm{eq}}}{T_{\mathrm{RH}}}\right)}^{1/3}.$$

$${\sigma }_v$$

$$\propto \sqrt{\frac{M_X^{2/3}}{a_f}},$$

$${\lambda }_{\mathrm{fs}}={\int }_{t_{\mathrm{kd}}}^{t_0}\frac{v}{a}\mathrm{d}t={\sigma }_v{a}_{\mathrm{RH}}{\int }_{a_{\mathrm{RH}}}^1\frac{\mathrm{d}a}{a^{3}H(a)}.$$

$$\Gamma =\frac{⟨\sigma v⟩}{2m_X^2}\sum _i\int {\rho }_{X,i}^2(\vec{r}){\mathrm{d}}^3\vec{r},$$

$$\Gamma =\frac{⟨\sigma v⟩}{2m_{\chi }^2}{N}_h{J}_h,$$

$$\Gamma =\frac{⟨\sigma v⟩}{2m_X^2}\int {\mathrm{d}}^3\vec{r}{\int }_0^{M_{\mathrm{RH}}}\mathrm{d}\ln M_f\frac{{\rho }_m(\vec{r})}{M_f}{\frac{df}{d\ln M}|}_{z_f}{J}_h(M_f),$$

$$J_h$$

Full text

††preprint: FERMILAB-PUB-19-249-A

I Introduction

The null results of direct detection experiments Aprile et al. (2018); Akerib et al. (2017); Cui et al. (2017), as well as searches for new physics at the Large Hadron Collider Sirunyan et al. (2019, 2018a, 2018b, 2018c, 2018d, 2018e); Aaboud et al. (2018a, b, c, d, 2017), have provided motivation to consider scenarios in which the dark matter does not couple directly to the particle content of the Standard Model (SM). In particular, it is plausible that the dark matter particles do not couple to the SM, but instead annihilate into one or more unstable species that reside within a hidden sector Pospelov et al. (2008); Arkani-Hamed et al. (2009); Abdullah et al. (2014); Berlin et al. (2014); Martin et al. (2014); Hooper et al. (2012). If these annihilation products are short-lived, the evolution of the early universe will closely resemble that found in the standard case of thermal freeze-out. If they are long-lived, however, such hidden sector annihilation products could become non-relativistic and come to dominate the energy density of the universe before they decay. In such a scenario, the decays of these particles would reheat the universe, producing a new bath of SM radiation and diluting the abundances of all previously formed relics, including that of the dark matter Berlin et al. (2016a, b) (see also Refs. Dror et al. (2016, 2018)).

Long lifetimes are easily realized in hidden sector models, in particular if the decaying particle is the lightest hidden sector state. Furthermore, if the coupling between the hidden sector and the SM is very small, these two sectors will remain decoupled from one another throughout thermal history. Such a picture is a relatively generic feature of models in which the dark matter freezes out through annihilations within a highly decoupled hidden sector. While a variety of simple and well-motivated hidden sector dark matter scenarios have been proposed, it has often proven difficult to constrain or otherwise test these models. In particular, the feeble couplings between the hidden and SM sectors makes the prospects for dark matter searches at colliders and direct detection experiments very bleak.

The prospects for indirect detection in hidden sector dark matter models are less uniformly negative, but they depend strongly on the thermal history of the hidden sector. In scenarios in which the unstable hidden sector particles decay while they are still relativistic, the dark matter freezes out during radiation domination, requiring an annihilation cross section near $\sigma v\simeq 2\times 10^{-26}$ cm${}^{3}/$s in order to yield the measured dark matter abundance Steigman et al. (2012). In this case, measurements by gamma-ray telescopes and cosmic-ray detectors can constrain a wide range of such models, for dark matter masses up to $\sim$ $100$ GeV Albert et al. (2017); Cholis et al. (2019); Bergstrom et al. (2013). If the hidden and SM sectors are more weakly coupled, however, the eventual out-of-equilibrium decay of the hidden sector would dilute the dark matter, requiring a smaller dark matter annihilation cross section to obtain the measured abundance and suppressing the predicted indirect detection signal.

In this article, we study the growth of small-scale structure in hidden sector dark matter scenarios and discuss the impact of such structure on the prospects for indirect detection. During an early matter-dominated era (EMDE), matter density perturbations grow linearly, in contrast to the logarithmic growth that occurs during standard radiation domination Erickcek and Sigurdson (2011); Barenboim and Rasero (2014); Fan et al. (2014); Erickcek (2015). This rapid growth enhances small-scale inhomogeneity and can trigger the formation of numerous sub-earth-mass microhalos long before structure would form in the absence of an EMDE. The microhalo abundance at early times is primarily limited by the temperature of the dark matter particles: the velocity dispersion of dark matter erases perturbations below the free streaming scale. For weakly interacting dark matter, the first microhalos generally form at redshifts less than 1000 Erickcek et al. (2016), but microhalos can form much earlier in hidden sector scenarios because the dark matter particles are not in kinetic equilibrium with the SM and are therefore colder.

In many scenarios we consider, perturbation growth during the EMDE leads to structure formation long before the universe becomes matter dominated, and we develop a theory of gravitational collapse during radiation domination to predict the properties of these microhalos. Although the presence of these microhalos would not have any discernible impact on the observable structure of our universe, the extraordinarily high dark matter annihilation rates within their dense cores leads to detectable fluxes of gamma rays and cosmic rays. We also point out that the angular distribution of this gamma-ray signal is distinct from that predicted from ordinary annihilating dark matter, instead mimicking that normally associated with decaying dark matter particles. We show that the annihilation signal can be described by an effective dark matter lifetime, and we use the Fermi Collaboration’s measurement of the high-latitude gamma-ray background Ackermann et al. (2015) and resulting constraints on dark matter decay Blanco and Hooper (2019) to constrain these scenarios.

The lack of dark matter free streaming also opens up the possibility that even smaller microhalos could form during the EMDE if it lasts long enough for perturbations to collapse. We show that these halos would be destroyed when most of their content decays into relativistic particles. The freed dark matter particles would be released with much larger velocities than they had before their infall into a halo, and their subsequent free streaming suppresses the formation of microhalos after the EMDE. This gravitational heating effectively erases the signatures of the EMDE, making these scenarios difficult to constrain.

In the following section, we describe the thermal history of the early universe within the context of a vector portal hidden sector model, including the freeze-out of the dark matter, the subsequent EMDE, and the reheating of the SM bath through the decay of the unstable hidden sector matter. We then discuss in Secs. III and  IV the evolution of density perturbations during and after the EMDE, leading to the formation of a large population of sub-earth-mass microhalos. In Sec. V we discuss how structure formation during the EMDE suppresses microhalo formation after the EMDE. Our main results are presented in Sec. VI, in which we calculate the dark matter’s annihilation rate and place constraints on this class of models from measurements of the high-latitude gamma-ray background. Finally, we summarize our results and conclusions in Sec. VII.

II Vector Portal Dark Matter

A variety of hidden sector models have been described in the literature, including those in which the hidden sector couples to the SM through the renormalizable interactions known as the vector portal ($B^{\mu\nu}$) Pospelov et al. (2008); Krolikowski (2008), the Higgs portal ($H^{\dagger}H$) Pospelov et al. (2008); Burgess et al. (2001); Davoudiasl et al. (2005); Bird et al. (2006); Kim and Lee (2007); Finkbeiner and Weiner (2007); D’Eramo (2007); March-Russell et al. (2008), and the lepton portal ($H^{\dagger}L$) Pospelov et al. (2008); Bai and Berger (2014). For concreteness, we will focus in this study on a vector portal model, although most of our results apply to each of these cases. We take our dark matter candidate to be a stable Dirac fermion, $X$, that couples to a gauge boson, $Z^{\prime}$:

[TABLE]

The first term in this expression describes the kinetic mixing between the $Z^{\prime}$ and SM $Z$/photon Holdom (1986). Through this kinetic mixing, the $Z^{\prime}$ acquires the following couplings to SM fermions:

[TABLE]

where $\theta_{W}$ is the weak mixing angle, $m_{Z}$ is the mass of the SM $Z$ and $g_{Y}$ and $g$ are the $U(1)_{Y}$ and $SU(2)_{L}$ gauge couplings, respectively. As a result of these couplings, the $Z^{\prime}$ decays to SM fermions with the following width:

[TABLE]

where $n_{c}$ is the number of colors of the final state fermions and $\beta_{f}\equiv\sqrt{1-4m^{2}_{f}/m^{2}_{Z^{\prime}}}$. For $m_{Z^{\prime}}\gg m_{Z},m_{f}$, this reduces to a lifetime given by

[TABLE]

Dark matter annihilation proceeds through $t$-channel $X$ exchange in this model to produce a $Z^{\prime}Z^{\prime}$ pair, with a cross section that is given by

[TABLE]

where $v$ is the relative velocity between the annihilating particles and $\alpha_{X}\equiv g^{2}_{X}/4\pi$.

For \epsilon\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}3\times 10^{-8}\times(T/{\rm GeV})^{1/2}\,(g_{\star}/75)^{1/4}, where $g_{\star}$ is the effective number of degrees-of-freedom at temperature $T$, interactions of the $Z^{\prime}$ will maintain equilibrium in the early universe between the hidden sector and the SM (see Appendix 7 of Ref. Berlin et al. (2016a)). For smaller values of $\epsilon$, however, the hidden sector will be entirely decoupled from the SM bath, allowing the two sectors to evolve independently and maintain their own temperatures. We define $\xi\equiv T_{h}/T$ as the ratio of the temperatures of the hidden and SM sectors. Following from entropy conservation, the value of $\xi$ evolves as

[TABLE]

where $g_{\star}$ and $g^{h}_{\star}$ are the numbers of effective relativistic degrees-of-freedom in the SM and hidden sectors, respectively, and the “inf” subscript denotes the initial value after inflation.

If the $Z^{\prime}$ population remains close to its equilibrium value during dark matter freeze-out (for details, see Ref. Berlin et al. (2016b)) and decays before becoming non-relativistic, the resulting dark matter thermal relic abundance is given by

[TABLE]

where $a$ and $b$ are terms in the expansion of the dark matter annihilation cross section, $\sigma v/2\approx a+bv^{2}+\mathcal{O}(v^{4})$, and $g^{\rm eff}_{\star}\equiv g_{\star}+g^{h}_{\star}\,\xi^{4}$ at freeze-out. $x_{f}$ is defined as the mass of $X$ divided by the SM temperature at freeze-out, and is found to have a value of $\sim 20\times\xi$ over a wide range of parameters.

The dark matter relic abundance can be very different from that described in the above expression if the two sectors are not maintained in equilibrium in the early universe. In this case, the energy density of the $Z^{\prime}$ population after they have become non-relativistic and prior to their decays is given by

[TABLE]

where $T$ is the temperature of the SM bath and $g_{\star,\mathrm{inf}}\simeq 106$. Included in this expression is factor of $13/6$ which is equal to one plus the ratio of the internal degrees of freedom associated with the $X$ and $Z^{\prime}$ (for details, see Ref. Berlin et al. (2016b)). The EMDE begins when this density exceeds that of the SM radiation bath ($\rho_{Z^{\prime}}>\rho_{\rm SM}=\pi^{2}g_{\star}T^{4}/30$), which occurs at the following temperature:

[TABLE]

The EMDE ends when $Z^{\prime}$ particles decay into SM particles, and this decay rate defines the reheat temperature:111Note that this definition of $T_{\rm RH}$ is equal to $1.22\times T_{f}$ as used in Ref. Berlin et al. (2016a).

[TABLE]

which roughly corresponds to the temperature at the onset of radiation domination following the EMDE. Setting this equal to the decay rate given in Eq. (3), we find that the SM bath is heated to the following temperature:

[TABLE]

and

[TABLE]

The decays of the $Z^{\prime}$ population dilute the relic abundance of dark matter by a factor of $T_{\rm dom}/T_{\rm RH}$, leading to the following final dark matter relic abundance Berlin et al. (2016a, b):

[TABLE]

where $\langle g_{\star}\rangle$ denotes the time-averaged value over the period of decay Kolb and Turner (1990).

III Growth of Perturbations during an Early Matter-Dominated Era

In this section, we study the growth of density perturbations during an era in the early universe in which the total energy density is dominated by non-relativistic hidden sector particles. Within this context, we consider a universe that contains three components: the dark matter, $X$, the unstable but relatively long-lived $Z^{\prime}$, and the relativistic SM bath. We assume that the universe was initially dominated by SM particles and that, during this epoch, the dark matter thermally decoupled from the $Z^{\prime}$ population while the $Z^{\prime}$ particles were still relativistic. The $Z^{\prime}$ population then became nonrelativistic and came to dominate the universe, leading to an EMDE. The duration of the EMDE and the subsequent reheating temperature are determined by the properties of the $Z^{\prime}$ particles, as described in the previous section.

The evolution of density perturbations during an EMDE has been studied extensively Erickcek and Sigurdson (2011); Barenboim and Rasero (2014); Fan et al. (2014); Erickcek (2015). In this study, we use the suite of linear perturbation equations in Newtonian gauge derived in Ref. Erickcek and Sigurdson (2011) for a pressureless scalar field that decays into relativistic particles and dark matter. In our case, however, the pressureless scalar field is replaced by the $Z^{\prime}$ particles, and these particles do not decay into dark matter particles. We also start our perturbation evolution prior to the EMDE, and we do not assume that all of the radiation present is sourced by the decay of the $Z^{\prime}$ population. Our initial conditions, therefore, differ from those used in Ref. Erickcek and Sigurdson (2011). In particular, assuming that the $Z^{\prime}$ population is nonrelativistic, superhorizon adiabatic initial conditions during radiation domination demand that the fractional density perturbation in the SM particles, $\delta_{R}\equiv(\rho_{R}-\bar{\rho}_{R})/\bar{\rho}_{R}$, is related to the curvature perturbation ($\delta_{R}=2\Phi$, where $\Phi$ is defined as in Ref. Erickcek and Sigurdson (2011)) and to the fractional perturbations in the densities of $X$ and $Z^{\prime}$ particles: $\delta_{X}=\delta_{Z^{\prime}}=(3/4)\delta_{R}$. While the perturbation mode is outside the horizon, the comoving divergence of each fluid’s physical velocity, $\theta\equiv a\,\partial_{i}(dx^{i}/dt)$, is equal to $-(\Phi/2)k^{2}/(aH)^{2}$, where $k$ is the wave number of the mode and $H$ is the Hubble parameter.

After the $Z^{\prime}$ population has become nonrelativistic and prior to the EMDE, $\delta_{X}$ and $\delta_{Z^{\prime}}$ each grow logarithmically with the scale factor after the perturbation mode enters the horizon (when $k=aH$). This growth is triggered by the gravitational pull toward overdense regions at horizon entry, but the decay of subhorizon gravitational potential perturbations during radiation domination implies that these forces quickly fade. The subsequent logarithmic growth in $\delta$ shown in Fig. 1 is a consequence of the continued drift of the dark matter particles toward regions that were initially overdense.222In the absence of gravitational forces, comoving velocities decay as $a^{-2}$, so the comoving displacement of a particle, $\vec{s}$, is proportional to ${\int a^{-2}dt=\int a^{-2}da/(aH)\propto\int da/a}$, since $H\propto a^{-2}$ during radiation domination. Thus, $\vec{s}\propto\ln a$, and since $\delta=-\nabla\cdot\vec{s}$ at linear order, $\delta\propto\ln a$. When the $Z^{\prime}$ population comes to dominate the universe at the start of the EMDE, the gravitational forces reappear and pull both the $X$ and $Z^{\prime}$ particles toward the regions that have more $Z^{\prime}$ particles than average, prompting both $\delta_{X}$ and $\delta_{Z^{\prime}}$ to transition to linear growth. At the end of the EMDE, $\rho_{Z^{\prime}}\simeq 0$, and so the perturbations in $\rho_{Z^{\prime}}$ are no longer relevant. As shown in Fig. 1, $\delta_{X}$ transitions to logarithmic growth after the EMDE and continues to grow logarithmically with the scale factor until it either approaches unity, which signifies the breakdown of linear theory, or the universe becomes matter dominated, at which point $\delta_{X}$ resumes linear growth. Fig. 1 also illustrates that the perturbation mode continues to grow rapidly for several e-folds of expansion after the EMDE ends, which implies that even a short EMDE can significantly enhance the amplitude of dark matter density fluctuations. This post-EMDE growth reflects the fact that the particles were accelerated toward the overdense regions during the EMDE, leading to a faster drift toward those regions after the EMDE ends.

Our suite of perturbation equations does not include baryons, which means that it does not accurately model the evolution of $\delta_{X}$ after the baryons become nonrelativistic. Baryons do not participate in structure growth on the small scales that enter the horizon prior to reheating (e.g. Bertschinger, 2006), and the growth of $\delta_{X}$ is suppressed as a result. To account for this suppression, we solve our perturbation equations until logarithmic growth is established after reheating ($a\simeq 1000\,a_{\mathrm{RH}}$, where $a_{\rm RH}$ is the scale factor at the time of reheating). We then match this evolution to the solutions to the Meszaros equation provided in Ref. Hu and Sugiyama (1996) that apply while the baryons are coupled to the photons and therefore do not participate in structure growth. We use these solutions to evolve $\delta_{X}$ past $a=1000\,a_{\mathrm{RH}}$.

Perturbation modes that enter the horizon while the $Z^{\prime}$ population is still relativistic are not accurately described by our suite of perturbation equations. While the $Z^{\prime}$ particles are relativistic, their pressure will prevent perturbation growth, and $\delta_{Z^{\prime}}$ will oscillate. Since the $X$ particles remain kinetically coupled to the $Z^{\prime}$ particles, they will inherit these oscillations, but the oscillations in $\delta_{X}$ will be damped Bertschinger (2006). We account for the suppression of perturbations due to these damped oscillations by introducing an exponential cut-off in the matter power spectrum: $P(k)\propto\exp[-k^{2}/k_{Z^{\prime}}^{2}]$, where $k_{Z^{\prime}}$ is the wavenumber of the mode that enters the horizon when the temperature of the hidden sector equals $m_{Z^{\prime}}$. In Sec. VI, we explore how altering this cut-off scale impacts the microhalo population and the resulting gamma-ray constraints.

IV Structure Formation Following an Early Matter-Dominated Era

The perturbation mode whose evolution is shown in Fig. 1 enters the nonlinear regime while the universe is radiation dominated. To track the evolution of this perturbation, we generalize the spherical collapse model Peebles (1980); Padmanabhan and Subramanian (1993) so that we can apply it to cosmologies that include an EMDE. The evolution of a spherical shell of matter that surrounds an overdense region that has an average matter density that is $(1+\delta)$ times the background density is governed by the second Friedmann equation: its physical radius, $R$, obeys

[TABLE]

where $\bar{\rho}_{\Lambda}$ and $\bar{\rho}_{R}$ are the energy densities associated with the cosmological constant and the relativistic SM particles respectively.333For simplicity, we do not include changes to the radiation’s composition in our spherical collapse formulation. To account for such variations, the $2\bar{\rho}_{R}$ term should be replaced with $\bar{\rho}_{R}+3\bar{P}_{R}$, where $\bar{P}_{R}$ is the radiation pressure. These densities are unperturbed; we neglect the perturbations in the relativistic particles because these perturbations are suppressed after reheating Erickcek and Sigurdson (2011); Barenboim and Rasero (2014); Fan et al. (2014). The average density of $X$ and $Z^{\prime}$ particles inside the shell is given by

[TABLE]

where a subscript $i$ indicates that the quantity is evaluated at some initial time, and $\bar{\rho}$ refers to the homogenous background density. This equation accounts for the decay of the $Z^{\prime}$ particles by multiplying the mass of $Z^{\prime}$ particles enclosed within the shell by the change in their comoving number density. The initial time is chosen to be during the EMDE such that $\delta_{i}\ll 1$ according to linear theory. With this selection, the Zel’dovich approximation Zel’dovich (1970) provides an initial condition for $dR/dt$:

[TABLE]

Since the shells do not cross, any value of $R_{i}$ gives the same evolution for the density within the shell:

[TABLE]

Fig. 1 shows the evolution of $\delta_{\mathrm{SC}}$ that is obtained by solving Eq. (14) for $R(a)$ with these initial conditions. When the value of $\delta_{\mathrm{SC}}$ is much less than one, the spherical collapse model matches the linear theory, but $\delta_{\mathrm{SC}}$ increases faster than the linear solution as $\delta$ approaches unity and eventually diverges at a finite time. This divergence indicates that the spherical shell has collapsed to $R=0$, and we denote value of the scale factor at this time as $a_{c}$. We also define $\delta_{c}\equiv\delta(a_{c})$ as predicted by linear theory. For the perturbation mode shown in Fig. 1, $\delta_{c}=2.8$, which is significantly greater than its value during matter domination $(\delta_{c}=1.686)$. Thus Fig. 1 illustrates that gravitational collapse does occur while the universe is radiation dominated, but it occurs later than it would in a matter-dominated universe. By considering different values for $\delta_{i}$, corresponding to different perturbation wavenumbers, we can determine how $\delta_{c}$ varies with $a_{c}$: it equals 1.68 during the EMDE and then increases up to $\delta_{c}\simeq 2.5$ before falling back to $\sim\!\!1.68$ at a_{c}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}a_{\mathrm{eq}}. While the precise shape of $\delta_{c}(a_{c})$ depends on $T_{\mathrm{RH}}$, $\delta_{c}$ for $a>10^{-6}$ is nearly the same for all $T_{\mathrm{RH}}>2$ MeV. In Appendix B, we plot and provide a fitting function for $\delta_{c}(a_{c},T_{\mathrm{RH}})$.

Since dark matter particles do not experience any gravitational forces while the universe is radiation dominated, the collapse of a spherical shell does not necessarily indicate that a dark matter halo has formed; it is possible that the particles in that shell pass through the origin and then continue drifting out of the overdense region. To determine the outcome of gravitational collapse during radiation domination, we used the cosmological simulation code GADGET-2 Springel (2005) to simulate the evolution of overdense regions. We modified the code to include both radiation Delos et al. (2018a, b) and decaying $Z^{\prime}$ particles. These components are modeled at the homogeneous level by adding corresponding terms to the Hubble rate appearing in the time integrals used for particle drifts and gravitational kicks. The clustering capacity of the $Z^{\prime}$ particles is modeled by scaling the mass, $m$, of each simulation particle as

[TABLE]

where $m_{0}$ is the particle mass at late times Suto (1987); Enqvist et al. (2015); Dakin et al. (2019). Effectively, each simulation particle represents $Z^{\prime}$ and $X$ particles together.444As a technical matter, this modification is also implemented by modifying a time integral. In particular, the time integrand used for gravitational kicks is scaled by $m/m_{0}$. We verified that this simulation matches the predictions of linear theory while $\delta\ll 1$.

The volume of our simulation box is ($10^{-9}h^{-1}$ Mpc)3 and we start the simulation at $a=3\times 10^{-13}$ in a scenario with $T_{\mathrm{RH}}=22$ MeV, which implies that the EMDE ends and radiation domination begins at $a=10^{-11}$. We assume that the dimensionless power spectrum scales as $\mathcal{P}(k)\propto\ln(k/k_{0})^{2}$ with $k_{0}\simeq 6\times 10^{6}h\,\mathrm{Mpc^{-1}}$ for the relevant perturbation wavenumbers; that is, we consider scales that entered the horizon during a period of radiation domination prior to the EMDE. We then draw individual proto-halo density peaks from the Gaussian-filtered power spectrum,

[TABLE]

using the density profiles of Ref. Bardeen et al. (1986) with smoothing scale $R_{f}=10^{-8}h^{-1}$ kpc. To add some realism to this otherwise smooth density peak, we also superpose small-scale Gaussian noise drawn from the oppositely filtered power spectrum.

We simulated the evolution of overdense regions with several different initial density profiles, three-dimensional shapes, and collapse redshifts, $z_{c}$, where $1+z_{c}=a_{c}^{-1}$ as predicted by the spherical collapse model. For overdense regions that collapse well after reheating, these simulations revealed a universal evolution that varied very little between the simulations. First, the density profile of the overdense region developes a shallow cusp with ${\mathrm{d}\ln\rho/\mathrm{d}\ln r\simeq-1}$ at $z\simeq z_{c}$. The density profile of this proto-halo is well-fit by the function

[TABLE]

where $A\simeq 8.4$, $B\simeq 12.4$, and $\bar{\rho}_{M}$ is the cosmological background matter density. When the average matter density within $R_{f}$ exceeds the homogeneous radiation density, the proto-halo re-collapses to form a much denser halo. This evolution is depicted in Fig. 2. In contrast, overdense regions that collapse too soon after reheating do not re-collapse to form a halo because the dark matter particles are moving fast enough to pass through the overdense region and escape the gravitational forces induced by the growing matter density. For T_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}10 MeV, our simulations indicate that regions with z_{c}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}10^{8} re-collapse.

Since our simulations only include a single overdense region with vacuum boundary conditions, our simulated halos do not develop the double power-law density profiles seen in cosmological simulations that start from a Gaussian density field. However, since the re-collapse of the proto-halo occurs in a locally matter-dominated environment, is seems reasonable to assume that the resulting halo should be similar to a halo that forms shortly after matter-radiation equality. Refs. Delos et al. (2018a, b) showed that halos forming at $z\simeq 1000$ in cosmological simulations with Gaussian initial conditions realized from an enhanced scale-invariant power spectrum have Navarro-Frenk-White (NFW) profiles Navarro et al. (1997). We assume that the halos that form during radiation domination also develop NFW profiles and that the concentration parameter of the resulting halo is determined by the redshift of re-collapse, not $z_{c}$. We use the density profile given in Eq. (20) to determine the re-collapse redshift, $z_{rc}$, from $z_{c}$ ($z_{rc}$ is the redshift at which the average matter density within $R_{f}$ exceeds the radiation density, unless that is already true at $z_{c}$, in which case $z_{rc}=z_{c}$). To be conservative, we only consider z_{c}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}10^{6}, thereby ensuring that the collapsing region forms a gravitationally bound halo (as demonstrated in Fig. 2 for $z_{c}=3.9\times 10^{7}$). Restricting ourselves to halo formation with $z_{c}<10^{6}$ also gives the halos extra time to establish NFW profiles after they re-collapse, and we explore how possible further delays in halo formation affect our gamma-ray constraints in Sec. VI.

V Structure Formation during the EMDE: The Impact of Gravitational Heating

If the duration of the EMDE is sufficiently long, then the density perturbations among the $X$ and $Z^{\prime}$ particles will become large enough to form gravitationally bound structures prior to reheating. The characteristic initial density perturbation is $\delta\simeq 10^{-5}$ at horizon entry. Both $\delta_{X}$ and $\delta_{Z^{\prime}}$ grow by roughly a factor of 100 at horizon entry and then linearly with the scale factor during the EMDE (as seen in Fig. 1), so we expect significant structure formation if a_{\mathrm{RH}}/a_{\mathrm{dom}}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}10^{3}. Since $\rho_{Z^{\prime}}\propto a^{-3}$ during the EMDE and $\rho_{Z^{\prime}}(a_{\mathrm{dom}})=\rho_{\mathrm{SM}}(a_{\mathrm{dom}})$ and $\rho_{Z^{\prime}}(a_{\mathrm{RH}})\simeq\rho_{\mathrm{SM}}(a_{\mathrm{RH}})$, $a_{\mathrm{RH}}/a_{\mathrm{dom}}\simeq(g_{*,\mathrm{dom}}/g_{*,\mathrm{RH}})^{1/3}(T_{\mathrm{dom}}/T_{\mathrm{RH}})^{4/3}$. Therefore, we expect that most of the hidden sector particles will be in halos at reheating if T_{\mathrm{dom}}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}200\,T_{\mathrm{RH}}.

A more precise determination of the halo abundance can be obtained from the Press-Schechter formalism Press and Schechter (1974), which provides an expression for the the fraction of dark matter that is bound into halos as a function of halo mass:

[TABLE]

where $\sigma^{2}(M,z)$ is the variance of $\delta\rho/\bar{\rho}$ in spheres that contain mass $M$ in $X$ and SM particles evaluated using a sharp-$k$ filter at $k_{s}=2.5/[3M/(4\pi\rho_{m,0})]^{1/3}$, as advocated by Ref. Benson et al. (2013). Since we are interested in structure formation during a matter-dominated era, $\delta_{c}=1.686$. We integrate Eq. (21) over all halo masses to obtain the total fraction of matter contained in halos at reheating. This computation confirms that halos are common for T_{\mathrm{dom}}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}200\,T_{\mathrm{RH}}. For example, if $T_{\mathrm{dom}}=470$ GeV and $T_{\mathrm{RH}}=1$ GeV, 20% of the matter is bound into halos at the time of reheating.

The halos present at reheating are predominantly composed of $Z^{\prime}$ particles; while $H\gg\Gamma_{Z^{\prime}}$, ${\rho_{X}}/{\rho_{Z^{\prime}}}\simeq T_{\mathrm{eq}}/T_{\mathrm{RH}}$. When the $Z^{\prime}$ population decays into relativistic SM particles, their energy will leave the halo. The remaining $X$ particles will escape the halo if their orbits cannot adiabatically adjust to the reduction in the halo’s mass. The characteristic dynamical time of a halo is $t_{\mathrm{dyn}}\simeq(G\rho)^{-1/2}$, while the characteristic mass-loss timescale is $t_{\Delta M}=M/\dot{M}=(\rho_{X}+\rho_{Z^{\prime}})/(\Gamma_{Z^{\prime}}\rho_{Z^{\prime}})$. For a halo that forms at $a_{f}$, the ratio of these timescales is as follows:

[TABLE]

where $\Delta t$ is the time since the halo’s formation. This ratio falls below unity for all halos with

[TABLE]

Therefore, we expect that all halos that form at scale factors greater than $10^{-6}a_{\mathrm{RH}}$ will be destroyed at reheating and that the $X$ particles within these halos will be released to free stream with the velocity they had within the halo. In this way, the formation of structure during the EMDE effectively heats the dark matter, imparting the $X$ particles with much larger velocities than they otherwise would have had.

The typical particle velocity within a halo is given by the virial velocity dispersion, $\sigma_{v}\simeq\sqrt{GM/R}$, where $R$ is radius within which the average density is 200 times the background density and $M$ is the mass enclosed in this radius. From this, it follows that

[TABLE]

where $M_{X}$ is the mass in $X$ particles. The formation time of a halo is determined by the amplitude of density fluctuations with $k\simeq[3M_{X}/(4\pi\rho_{X,0})]^{-1/3}$; halos with this mass will typically form when $\delta\simeq\delta_{c}$ for this wavenumber. For modes that enter the horizon during an EMDE, $\delta_{X}\propto k^{2}$ Erickcek and Sigurdson (2011). Since these modes grow linearly with the scale factor during the EMDE, $a_{f}\propto k^{-2}\propto M_{X}^{2/3}$ and $\sigma_{v}\simeq 2000$ km/s, independent of halo mass. The amplitudes of modes that enter the horizon before the EMDE starts do not increase as rapidly with increasing $k$ because these modes only grow logarithmically between horizon entry and the start of the EMDE. As a result, $a_{f}$ does not decrease as fast with increasing $M_{X}$, and $\sigma_{v}$ increases as $M_{X}$ increases. Therefore, the largest halos that form during the EMDE will always have the largest virial velocities.

After the dark matter particles are released from the halos that evaporate at reheating, they will free stream in random directions, thereby suppressing any density perturbations with wavelengths less than the free streaming horizon, $\lambda_{\mathrm{fs}}$. If all the dark matter particles have the same velocity dispersion, the matter power spectrum is suppressed by a factor of $\exp[-k^{2}\lambda_{\mathrm{fs}}^{2}]$ Green et al. (2005); Bertschinger (2006); Boyarsky et al. (2009), where

[TABLE]

To probe how this gravitational heating of dark matter affects the formation of microhalos after the EMDE, we make two estimates. For the “optimistic” estimate, we only consider the coldest 80% of the dark matter. If less than 20% of the $X$ particles are bound into halos at reheating, we impose no free streaming cutoff to the matter power spectrum. Otherwise, we find the halo mass, $M_{20}$, such that 20% of the $X$ particles are bound into halos with $M>M_{20}$ at reheating. At $a=0.8\,a_{\mathrm{RH}}$ about 10% of the matter is bound into such halos, so we take this to be the formation time of these halos. We then use Eq. (24) to determine the virial velocity within a halo with a mass equal to $M_{20}$ that formed at $a_{f}=0.8\,a_{\mathrm{RH}}$ and apply a $\exp[-k^{2}\lambda_{\mathrm{fs}}^{2}]$ cutoff to the matter power spectrum with $\lambda_{\mathrm{fs}}$ given by Eq. (25). We also assume that only 80% of the dark matter participates in structure formation after the EMDE and reduce the growth rate of perturbations accordingly. For the “pessimistic” estimate, we follow this same procedure for the coldest 90% of the dark matter. Since $M_{10}>M_{20}$, this procedure yields a larger value of $\sigma_{v}$ and suppresses perturbations on larger scales.

Both of these estimates are conservative in that most of the dark matter does not free stream to this extent, but it is also possible that the actual cut-off in the matter power spectrum differs from $\lambda_{\mathrm{fs}}$. Unfortunately, the boost factor from microhalos that form after an EMDE is extremely sensitive to the cut-off scale in the matter power spectrum; changing this scale by a factor of two can change the boost factor by more than an order of magnitude Erickcek (2015). The contrast between the optimistic and pessimistic estimates gives some idea of the uncertainty associated with gravitational heating, but it is certainly possible that both approaches overestimate how much gravitational heating suppresses the microhalo population. More comprehensive numerical simulations of halo evolution during and after the EMDE could allow us to extend our constraints to larger values of $T_{\mathrm{dom}}/T_{\mathrm{RH}}$.

VI Microhalo-dominated Dark Matter Annihilation

If most of the dark matter is concentrated within a population of microhalos, then the dark matter annihilation rate in a given volume of space will be given by

[TABLE]

where the sum is over the microhalos within the volume and $\rho_{X,i}(\vec{r})$ is the density profile of microhalo, $i$. If, for illustration, we make the simplifying assumption that each microhalo is identical, we find that the dark matter annihilation rate is proportional to the number of microhalos, $N_{h}$, within that volume:

[TABLE]

where $J_{h}\equiv\int\rho_{X,h}^{2}(\vec{r})\,\mathrm{d}^{3}\vec{r}$.

To determine the fraction of dark matter that is bound into halos as a function of halo mass, we use the Press-Schechter formalism Press and Schechter (1974), as described in the previous section. If the halos present at a redshift $z_{f}$ survive to the present day, then the annihilation rate in a volume will be given by

[TABLE]

where $\rho_{m}$ is the matter density and $M_{\rm RH}$ is the mass enclosed in the cosmological horizon at reheating. We only consider the microhalos with $M<M_{\mathrm{RH}}$ because these microhalos form much earlier than larger microhalos and are therefore more likely to survive accretion into larger halos Erickcek (2015).

If the microhalos present at redshift $z_{f}$ have NFW profiles with concentration $c=r_{\mathrm{200}}/r_{s}$, where $r_{\mathrm{200}}$ is the radius within which the average matter density is 200$\bar{\rho}_{M}$, then

[TABLE]

where $f_{X}$ is the fraction of mass that is dark matter. We assume that all microhalos have $c=2$ at a redshift $z_{f}$ because this is the lowest concentration seen in microhalo simulations Anderhalden and Diemand (2013); Ishiyama (2014). We will revisit this assumption when we discuss uncertainties in our computation of gamma-ray emission later in this section.

Inserting Eq. (29) into Eq. (28) gives

[TABLE]

If we define $f_{\mathrm{tot}}$ to be the fraction of mass contained in halos with $M<M_{\mathrm{RH}}$, the annihilation rate in a volume divided by the dark matter mass within that volume is

[TABLE]

where $B_{0}=13.76\times{200}(1+z_{f})^{3}f_{\mathrm{tot}}(z_{f})/3$ is the boost due to the microhalos. As in Ref. Erickcek (2015), we choose $z_{f}$ such that $B_{0}$ is maximized. While this estimate of the boost factor assumes that all the microhalos present at a redshift of $z_{f}$ survive until the present day, it also neglects the annihilation boost due to any subhalos present within these microhalos and any microhalos that form at $z<z_{f}$. Since $B_{0}$ only depends on the fraction of dark matter bound into microhalos, and not on the size of the microhalos, it does not depend strongly on the reheat temperature, but it is very sensitive to the ratio $k_{\mathrm{cut}}/k_{\mathrm{RH}}$, where $k_{\mathrm{cut}}$ is the scale that sets the cut-off in the matter power spectrum, $P(k)\propto\exp[-k^{2}/k_{\mathrm{cut}}^{2}]$  Erickcek (2015). In the absence of gravitational heating, $k_{\mathrm{cut}}$ is the scale that enters the horizon when $T_{h}=m_{Z^{\prime}}$. Otherwise it is the free-streaming scale determined by the velocities of the $X$ particles.

In Fig. 3, we plot the value of the boost factor, $B_{0}$, as a function of $T_{\rm dom}/T_{\rm RH}$, for the case of $\alpha_{X}=0.01$ (dashed) and $\alpha_{X}=0.1$ (solid). At each point, we adopt $\xi_{\rm inf}=1$, $m_{X}/m_{Z^{\prime}}=20$, and select the values of $\epsilon$ and $m_{X}$ such that we obtain the desired dark matter abundance, $\Omega_{X}h^{2}\simeq 0.11$. In each case the grey and blue curves represent the results found following our default and conservative procedures, respectively (as described later in this section). Applying our computation of $B_{0}$ to a power spectrum derived assuming no deviations from radiation domination in the early universe gives $B_{0}\simeq 10^{6}$ for a very broad range of minimum halo masses. In scenarios with a very long EMDE, gravitational heating significantly suppresses the microhalo population, reducing the boost factors to this value. We also see in Fig. 3 that very short EMDEs still increase the boost factor: the presence of the $Z^{\prime}$ particles enhances the growth of dark matter perturbations even if they never fully dominate the energy content of the universe.

The difference between the solid ($\alpha_{X}=0.01$) and dashed ($\alpha_{X}=0.1$) curves in Fig. 3 reflects the fact that $B_{0}$ depends on both $T_{\mathrm{RH}}$ and $T_{\mathrm{dom}}$. The enhancement to small-scale perturbations due to the EMDE only depends on the ratio $k_{\mathrm{dom}}/k_{\mathrm{RH}}$, where $k_{\mathrm{dom}}$ is the wave number of the mode that enters the horizon when $\rho_{Z^{\prime}}=\rho_{R}$ and $k_{\mathrm{RH}}$ is the wavenumber of the mode that enters the horizon at reheating, as defined in Ref. Erickcek (2015). This ratio is determined by $T_{\mathrm{dom}}$ and $T_{\mathrm{RH}}$:

[TABLE]

for ${T_{\mathrm{dom}}}/{T_{\mathrm{RH}}}\gtrsim 5$. For fixed $k_{\mathrm{dom}}/k_{\mathrm{RH}}$, $B_{0}$ still has a weak dependence on $T_{\mathrm{RH}}$ due to the red tilt of the primordial power spectrum and the duration of the radiation-dominated era, but this slight variation is less than a factor of 5 for the range of $T_{\mathrm{RH}}$ we consider.

Since the microhalos are expected to track the dark matter density (as opposed to the square of the density), the emission signatures in these scenarios are more akin to decaying dark matter than traditional annihilating dark matter. With this in mind, we will use measurements of the high-latitude gamma-ray background to constrain these models. The annihilation rate per mass given by Eq. (31) can be converted to an effective lifetime by considering the rate of jet production in a fixed volume containing dark matter mass, $M_{X}$.555We use the phrase “jet production” here to denote the production of any primary annihilation product, and can include leptons which do not result in the formation of a QCD jet. Dark matter annihilations produce jets of energy $m_{X}$ at a rate of $2(\Gamma/M_{X})M_{X}$. To get the same jets from decaying dark matter, the particles must have mass $2m_{X}$, which implies that the jet production rate is $(2/\tau)[M_{X}/(2m_{X})]$. Equating these rates yields the following effective lifetime:

[TABLE]

This lifetime should be compared to the lower bound on the dark matter lifetime for particles with mass $2m_{X}$.

In Fig. 4, we plot the effective lifetime of the dark matter as a function of $m_{X}$, for the same parameter values (and the same line types) as adopted in Fig. 3. These curves do not extend to arbitrarily low $m_{X}$ values because we only show results for $k_{\rm dom}/k_{\rm RH}>1.5$, which corresponds to requiring that the density of $Z^{\prime}$ particles exceeds the density of SM particles at some point prior to their decay. We limit our analysis to scenarios with an EMDE because the annihilation rate scales with the number density of microhalos (and thus with the density of dark matter) only if the boost factors are very large. In a scenario without an EMDE, constraints based on gamma-ray observations of dwarf spheroidal galaxies Albert et al. (2017) should instead be applied to the model.

The uncertainties in our computation of the boost factor, $B_{0}$, arise from two main sources: the cut-off in the matter power spectrum ($P(k)\propto\exp[-k^{2}/k_{\mathrm{cut}}^{2}]$) and the characterization of the halo density profiles. The first cut-off scale we consider is set by the mass of the $Z^{\prime}$ particle: $k_{\mathrm{cut}}=k_{Z^{\prime}}$ where $k_{Z^{\prime}}$ is the wavenumber of the mode that enters the horizon when the temperature of the hidden sector is equal to $m_{Z^{\prime}}$. This cut-off only approximates the suppression of perturbations on scales that enter the horizon while the $Z^{\prime}$ population is relativistic and it is possible that the actual power spectrum has a slightly different shape. This change will not affect scenarios with 30\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}70 because most of the dark matter is already bound into microhalos at $a=10^{-6}$ in these scenarios, so changing the power spectrum has a minimal effect on the value of $B_{0}$ obtained from the halo population at this time.

If modes that enter the horizon while $T_{h}\simeq m_{Z^{\prime}}$ are not as damped as this model predicts, $B_{0}$ will increase for k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}30 because microhalos will form earlier. However, more structure will also form prior to reheating, so gravitational heating will suppress $B_{0}$ at smaller values of $k_{\mathrm{dom}}/k_{\mathrm{RH}}$. For example, if $k_{\mathrm{cut}}=2k_{Z^{\prime}}$, gravitational heating will suppress $B_{0}$ for k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}70 instead of k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}80. Conversely, if $k_{\mathrm{cut}}=k_{Z^{\prime}}/2$, the impact of gravitational heating will be pushed to larger values of $k_{\mathrm{dom}}/k_{\mathrm{RH}}$, but $B_{0}$ will be suppressed by about an order of magnitude for k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}30 due to the later formation of microhalos. For the conservative $B_{0}$ calculation, we adopt $k_{\mathrm{cut}}=k_{Z^{\prime}}/2$ for k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}30 to account for this suppression.

The second cut-off scale we consider arises from the suppression of perturbations following gravitational heating. We computed this suppression based on the evolution of the coldest 80% and 90% of the dark matter to illustrate the uncertainty associated with structure formation durning the EMDE, but it is possible that both of these calculations overestimate the suppression of perturbations due to gravitational heating. Simulations of halo evaporation at reheating could reveal that the free-streaming of particles released from halos does not result in $k_{\mathrm{cut}}=\lambda_{\mathrm{fs}}^{-1}$ as we have assumed. For scenarios with gravitational heating, $k_{\mathrm{RH}}\lambda_{\mathrm{fs}}^{-1}\simeq 40$; changing $k_{\mathrm{cut}}$ by a factor of two can change $B_{0}$ by up to two orders of magnitude but would still not alter our constraints. Determining the onset of gravitational heating is a more significant source of uncertainty. In our optimistic estimate, gravitational heating dramatically suppresses $B_{0}$ for k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}90 because $20\%$ of the dark matter is gravitationally heated in these scenarios. It is worth noting, however, that less than half of the dark matter is gravitational heated for k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}130. Therefore, it is possible that the formation of microhalos in these scenarios is not completely suppressed by the free streaming of released dark matter particles. For our default $B_{0}$ computation, we base the free streaming horizon on the coldest $80\%$ of the dark matter, while the conservative computation uses the coldest 90% of the dark matter. The difference in $B_{0}$ between these two approaches is roughly equivalent to increasing $k_{\mathrm{cut}}$ to $2k_{Z^{\prime}}$ while basing the free-streaming horizon on the coldest 80% of the dark matter.

Turning to the characterization of the halo density profiles, our procedure for calculating $B_{0}$ relies on the assumption that the boost factor is dominated by halos present at a certain redshift and that all halos have NFW profiles with concentration $c=2$ at this time. We choose the redshift that maximizes $B_{0}$ by striking the optimal compromise between increased microhalo abundance and increased microhalo density. We then assume that these halos survive to the present day and that any halos that form subsequently make a minimal contribution to the boost factor. Remarkably, the boost factor computed by this simple approach applied to a standard matter power spectrum (i.e. no deviations from radiation domination at early times) nearly matches the results of more complicated halo-model approaches that use concentration-mass relations and (sub)halo mass functions calibrated using $N$-body simulations. The $\zeta(z)$ boost factor used in Ref. The Fermi LAT Collaboration (2015) is related to $B_{0}$ by $B_{0}=\zeta(z)(1+z)^{3}$. For a minimum halo mass of $10^{-6}h^{-1}\,M_{\odot}$, $\zeta(z)(1+z)^{3}\simeq 4\times 10^{5}$, whereas our method applied to a standard power spectrum yields $B_{0}\simeq 10^{6}$. The discrepancy of a factor of $\sim 2-3$ could easily be attributed to the destruction of subhalos.

Nevertheless, we should question the density profiles assumed by our calculation. The adoption of an NFW profile is probably overly conservative, as numerous simulations have indicated that the first generation of halos have steeper inner density profiles Anderhalden and Diemand (2013); Ishiyama (2014); Delos et al. (2018a, b, 2019). However, simulations also show that it takes considerable time after the halo’s formation for the NFW profile to stabilize: a halo that forms at $z=1000$ from initial conditions drawn from an enhanced plateau in the power spectrum has an NFW profile with $c=2$ at $z\simeq 400$ Delos et al. (2018a, b). The time interval between these two redshifts corresponds to $10\,t_{\mathrm{dyn}}$, where $t_{\mathrm{dyn}}=\sqrt{3\pi^{2}/(16G\rho_{\mathrm{vir}})}$ with $\rho_{\mathrm{vir}}=200\bar{\rho}_{M}(z_{f})$ is the dynamical time of the halo at formation. Waiting ten dynamical times after $z_{f}$ to evaluate the density profile reduces $B_{0}$ by a factor of 10 for halos that form during matter domination, which would apply to k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}10. For halos that form during radiation domination, the reduction factor depends on $z_{f}=z_{\mathrm{rc}}$, increasing to 100 for $a_{c}=10^{-6}$. However, if most of the dark matter is bound into halos at $a=10^{-6}$, as the case for 30\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}k_{\mathrm{dom}}/k_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{<\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}70, then it is safe to assume that halos were present at much earlier times and that our estimate of $B_{0}$ based on the halos present at $a_{f}=10^{-6}$ underestimates the true annihilation rate. Nevertheless, for our conservative $B_{0}$ computation, we reduce all $B_{0}$ values by the reduction factor implied by delaying the establishment of the NFW profile with $c=2$ ten dynamical times after the halos form.

The most stringent and broadly applicable constraints on this class of models are derived from the measurements of the Fermi Gamma-Ray Space Telescope. To derive these constraints, we follow the procedure described in Ref. Blanco and Hooper (2019), modified as described in Appendix A of this paper. The black curves shown in Fig. 4 represent the lower limits on $\tau_{\rm eff}$ in this model, as derived from Fermi’s measurement of the high-latitude gamma-ray background Ackermann et al. (2015). The dashed (solid) curves derive from the default (more conservative) procedure as described in Ref. Blanco and Hooper (2019). When these constraints are compared to the predicted values of $\tau_{\rm eff}$, we find that a wide range of dark matter masses are already ruled out. In particular, for the case of $\alpha_{X}=0.01$ (0.1), this excludes any scenario with a significant EMDE for masses up to $m_{X}\sim 6-7$ TeV (60-70 TeV). For higher masses, gravitational heating suppresses the abundance of microhalos, reducing the boost factors to a level consistent with Fermi’s measurements.

In Fig. 5, we plot some of the features of and constraints on this model, again for the case of $\xi_{\rm inf}=1$ and $m_{X}/m_{Z^{\prime}}=20$. The dashed curves represent the regions of this plane that yield a relic abundance equal to the measured dark matter density, $\Omega_{X}h^{2}\simeq 0.11$, for choices of $\alpha_{X}=0.01$ or $\alpha_{X}=0.1$. In the upper regions of this plane, the hidden sector never dominates the energy density of the early universe and there is no EMDE. For smaller values of $\epsilon$, however, the $Z^{\prime}$ population is sufficiently long-lived that it comes to dominate the energy density and subsequently reheat the SM bath through its decays. The blue curves represent contours of constant $T_{\rm dom}/T_{\rm RH}$ (see Eq. 12), while the red curves are constant in $T_{\rm RH}$ (see Eq. 11). The lower left region of the plane is ruled out, as we require that the temperature at the onset of radiation domination exceeds 3 MeV in order to generate the neutrinos that are needed to produce the observed light element abundances Kawasaki et al. (1999, 2000); Hannestad (2004); Ichikawa et al. (2005), as well as the observed anisotropies in the cosmic microwave background and large-scale density perturbations Ichikawa et al. (2007); de Bernardis et al. (2008).

The grey region of Fig. 5 is currently ruled out by Fermi’s measurement of the high-latitude gamma-ray background, based on our conservative computation of $B_{0}$. Across this region of parameter space, the EMDE leads to the formation and survival of a large population of microhalos, resulting in large boost factors that are inconsistent with Fermi’s observations. Using our conservative $B_{0}$ values, short EMDEs ($T_{\mathrm{dom}}/T_{\mathrm{RH}}\lesssim 10$) do not sufficiently enhance the dark matter annihilation rate to contradict Fermi observations for $m_{X}\gtrsim 500$ GeV, but these scenarios are ruled out if we use our default $B_{0}$ computation, as seen in Fig. 4. In the white region in the lower right portion of this figure, however, gravitational heating suppresses the abundance of microhalos, reducing the boost factors and resulting gamma-ray emission to acceptable levels for both our conservative and default computations.

VII Summary and Conclusions

If the dark matter is part of a hidden sector that is only very feebly coupled to the Standard Model, the lightest particle within this sector would generically be expected to be long-lived. In such a model, these long-lived particles could come to dominate the energy density of the universe prior to their decays. During such an era of early matter domination, density perturbations would grow more quickly than predicted in standard radiation-dominated cosmology, leading to the formation of a large population of sub-earth-mass dark matter microhalos.

In this study, we have calculated the evolution of density perturbations in a class of hidden sector dark matter models and estimated how the enhancement of these perturbations during the EMDE affects the microhalo abundance. The key difference between this analysis and earlier considerations of microhalo formation following an EMDE with weakly interacting dark matter Erickcek (2015); Erickcek et al. (2016) derives from the fact that the dark matter within our hidden sector models does not couple directly to the particles in the Standard Model and is therefore much colder. As a result, it is possible to form microhalos while the universe is still radiation dominated, and we have employed both spherical collapse theory and $N$-body simulations to determine the formation time and properties of these microhalos. We have also studied the case in which the decay of the unstable hidden sector particles occurs after gravitationally bound structures have begun to form. The disruption of these halos at reheating effectively heats the dark matter and suppresses subsequent microhalo formation.

The existence of a large microhalo population can cause the dark matter to annihilate at a very high rate, producing detectable fluxes of gamma rays and other annihilation products. Unlike with typical dark matter annihilation, the annihilation rate in a microhalo-dominated scenario scales with the density of dark matter (as opposed to the square of the density), leading to indirect detection signals that resemble those ordinarily predicted from decaying dark matter. Using Fermi’s measurement of the high-latitude gamma-ray background, we are able to rule out a wide range of parameter space within this class of models. In particular, if there was an early matter-dominated era that persisted for less than a factor of $\sim\!\!10^{3}$ in scale factor, we predict that the dark matter will largely be bound within microhalos, leading to very large annihilation boost factors. In contrast, we find that scenarios with even longer early matter-dominated eras cannot be excluded, as the microhalo abundance is efficiently suppressed in this case by the gravitational heating of the dark matter prior to reheating.

Acknowledgements.

This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of High Energy Physics. The simulations for this work were carried out on the Dogwood computing cluster at the University of North Carolina at Chapel Hill. M. S. D. and A. L. E. were partially supported by NSF Grant No. PHY-1752752.

Appendix A Gamma-Ray Constraints

The angular distribution of the gamma rays produced through dark matter annihilating in microhalos is not the same as in the standard (smoothly distributed) scenario, but instead mimics the morphology predicted for the case of decaying dark matter particles. Over the range of dark matter models considered in this study, the strongest constraints are those derived from the Fermi Gamma-Ray Space Telescope’s measurement of the high-latitude gamma-ray background. In evaluating the constraints on the dark matter’s effective lifetime, we follow closely the approach described in Ref. Blanco and Hooper (2019), altering only the branching fractions that are appropriate for the model under consideration and comparing the effective lifetime to the lower bound on the dark matter lifetime for particles with mass of $2m_{X}$.

In Fig. 6, we plot the contribution to the isotropic gamma-ray background from dark matter annihilations in the hidden sector model under consideration for the case of $m_{X}/m_{Z^{\prime}}=20$. In each frame, the overall normalization (and corresponding dark matter annihilation rate) has been set to the maximum value allowed by the data (at the 95% confidence level), as evaluated following Ref. Blanco and Hooper (2019).

In Fig. 7, we plot the limits on the effective lifetime of the dark matter particles in this model (again, for the case of $m_{X}/m_{Z^{\prime}}=20$). Following Ref. Blanco and Hooper (2019), the solid curves treat the systematic errors as entirely independent and uncorrelated. At the other extreme, the dashed curve has been derived assuming that the systematic errors are fully correlated, moving upward or downward together in unison.

Appendix B Threshold for Collapse During Radiation Domination

During radiation domination, the linear collapse threshold is larger than its value during matter domination, $\delta_{c}>\delta_{c,0}\equiv 1.686$ . The behavior of $\delta_{c}$ is a function of both the scale factor at collapse as well as the reheat temperature, as shown in Fig. 8. In this appendix, we present a useful fitting formula for $\delta_{c}$.

To begin, we define $a_{\mathrm{RH}}$ as the scale factor corresponding to temperature $T_{\mathrm{RH}}$ assuming radiation domination. This definition implies that

[TABLE]

where $T_{0}=$2.3\text{\times}{10}^{-4}\text{,}\mathrm{e}\mathrm{V}$$ is the temperature today. Note that this differs from other definitions of $a_{\mathrm{RH}}$ that appear in the literature, and we only use it as a parameter for our fitting formula. The following fitting formula for $\delta_{\mathrm{c}}$ is accurate to within 0.004 for $2\text{ MeV}<T_{\mathrm{RH}}<1\text{ TeV}$ as long as $a_{c}\lesssim 0.1$ (dark energy is neglected) and ${a_{c}\gg T_{0}/T_{\mathrm{dom}}}$ (radiation prior to the EMDE is neglected):

[TABLE]

where $a_{\mathrm{eq}}=\Omega_{\mathrm{r}}/\Omega_{\mathrm{m}}$. The functions $f(x)$, $g(x)$ and $h(x)$ are defined as follows:

[TABLE]

[TABLE]

and

[TABLE]

This fitting form was calibrated assuming no changes in the relativistic content of the universe, and incorporating them reduces $\delta_{c}$ by of order 1%. The difference is largest if reheating occurs while all Standard Model particles are still relativistic (T_{\mathrm{RH}}\mathrel{\raise 1.29167pt\hbox{>\kern-7.5pt\lower 4.30554pt\hbox{\sim}}}1 TeV) and collapse occurs late in radiation domination ($a_{c}\sim 10^{-5}$), in which case $\delta_{c}$ is smaller by about 2%. Reducing $\delta_{c}$ increases the microhalo abundance and the resulting annihilation boost factor, but a $2\%$ reduction in $\delta_{c}$ increases $B_{0}$ by less than 20%.