# Simulations of Gravitational Heating Due to Early Matter Domination

Himanish Ganjoo11footnotetext: Corresponding author. and M. Sten Delos

###### Abstract

In cosmologies with an early matter-dominated era (EMDE) prior to Big Bang
nucleosynthesis, the boosted growth of small-scale matter perturbations during
the EMDE leads to microhalo formation long before halos would otherwise begin
to form. For a range of models, halos can even form during the EMDE itself.
These halos would dissipate at the end of the EMDE, releasing their
gravitationally heated dark matter and thereby imprinting a free-streaming
cut-off on the matter power spectrum. We conduct the first cosmological
$N$-body simulations of the formation and evaporation of halos during and
after an EMDE. We show that in these scenarios, the free-streaming cut-off
after the EMDE can be predicted accurately from the linear matter power
spectrum. Although the free streaming can erase much of the EMDE-driven boost
to density perturbations, we use our findings to show that the (re-)formation
of halos after the EMDE nevertheless proceeds before redshift $\sim 1000$.
Early-forming microhalos are a key observational signature of an EMDE, and our
prescription for the impact of gravitational heating will allow studies of the
observational status and prospects of EMDE scenarios to cover a much wider
range of parameters.

##  1 Introduction

The history of the universe between a purported period of inflationary
expansion [1, 2, 3] and Big Bang Nucleosynthesis (BBN) is poorly understood.
Standard cosmology assumes that radiation dominated the Universe between the
end of inflation and the onset of the last matter-dominated epoch. However, a
wide variety of models predict alternative histories of the early universe
[4], including a possible early matter-dominated era (EMDE). EMDEs can arise
due to oscillating string moduli fields whose energy density dominates the
early universe [5] or in theories in which long-lived hidden-sector particles
become the dominant component of the Universe [6, 7, 8, 9, 10, 11]. Such
scenarios also include cannibal domination [12]. Current observations only
require that the heavy species that drives an EMDE decay into Standard Model
(SM) radiation by a reheat temperature ${T}_{\mathrm{RH}}$ larger than a few
MeV [13, 14].

While a potentially weak SM coupling or high mass could leave this early
matter species inaccessible to terrestrial experiments, the EMDE can leave a
detectable gravitational imprint on the late-time Universe. Subhorizon matter
perturbations grow linearly with the scale factor during an EMDE, leading to a
small-scale enhancement to the matter power spectrum that can persist in the
dark matter after the end of the EMDE [15, 16, 17, 18]. This enhancement
results in the formation of nonlinear structures much earlier than in standard
cosmological scenarios [15, 16, 18]. While these small-scale dark matter
microhalos do not appreciably change the large-scale structure of the
universe, their dense cores reside as substructure in current-day halos and
boost the dark matter annihilation rate [18, 19, 20, 21]. In addition, these
structures can be detected via their effects on pulsar timing [22, 23, 24] or
through their microlensing of stars [25, 26, 27, 28].

If the EMDE is sufficiently long and the dark matter sufficiently cold, bound
microhalos can even form and accrete dark matter during the EMDE itself. Such
a scenario has been previously considered in the context of hidden-sector dark
matter models [20] but is also a likely outcome for more minimal dark matter
models (such as the supersymmetric models discussed by Ref. [19]; see Ref.
[21]). As the dominant particles decay, these early-forming structures
evaporate, and the dark matter particles residing in them shoot out at boosted
speeds in random directions. The free streaming that results from this
gravitational heating potentially erases some of the boosted small-scale
power, thereby suppressing subsequent microhalo formation [20, 29]. In order
to predict the observational consequences of these scenarios, it is necessary
to understand the nature of this suppression.

In this paper, we carry out the first cosmological $N$-body simulations of the
formation of structure during an EMDE and its subsequent evaporation as the
EMDE ends. This process was previously studied by Refs. [20, 29] using
analytic approximations. We find that the impact of gravitational heating can
be accurately modeled as a free-streaming cut-off to the matter power
spectrum, as suggested by Ref. [20]. Moreover, even though gravitational
heating is a nonlinear process, the cut-off scale is related in a simple way
to the matter power spectrum calculated in linear theory, which has been
studied in detail for EMDEs (e.g., Refs. [15, 30]). Consequently, the matter
power spectrum long after reheating can be predicted straightforwardly from
such analytic results. As a demonstration, we show that the suppression of
structure due to gravitational heating is sufficiently weak that dark matter
microhalos are expected to begin forming either before or not long after
matter-radiation equality, even in scenarios with an arbitrarily large degree
of structure formation during the EMDE.

This paper is organized as follows. In Section 2, we review the evolution of
the homogeneous background and density perturbations during an EMDE. Section 3
explains the formation of nonlinear structures during an EMDE and the idea of
gravitational heating due to their evaporation. In Section 4, we describe our
suite of $N$-body simulations and show qualitatively how halos form and
subsequently dissipate. Section 5 analyzes the free-streaming cut-off to the
power spectrum that arises in the simulations after reheating and develops a
connection between the scale of this cut-off and the linear power spectrum
during the EMDE. In Section 6, we use this connection to explore how
gravitational heating impacts the (re-)formation of structure long after
reheating, particularly focusing on how the time of microhalo formation is
affected. Finally, Section 7 summarizes our results and looks at future
prospects for extending and applying them. This paper uses natural units
throughout, in which $c=\hbar=1$.

##  2 Background and Perturbation Evolution

###  2.1 The Homogeneous Background

Although we expect the results of this work to be applicable to EMDEs in
general, for concreteness we adopt the class of hidden-sector models
considered in Ref. [30] (hereafter G23). These models consist of a universe
with three fluids: the dark matter ($X$), the $Y$ particles that drive the
EMDE, and the SM radiation (denoted by the subscript $R$). The $X$ and $Y$
particles reside in a hidden sector kinetically decoupled from the SM, with
temperature $T_{\rm hs}$, while the temperature of the SM is $T$. The $X$
particles are nonrelativistic with their energy density evolving with scale
factor as $\rho_{X}\propto a^{-3}$. The $Y$ particles are initially
relativistic but become nonrelativistic as $T_{\rm hs}$ decreases. The
nonrelativistic $Y$ particles come to dominate the energy density of the
Universe and cause an epoch of early matter domination, but they decay into SM
radiation with a rate $\Gamma$. Although we do not appeal to any dark matter-
SM interaction in this work, we note that in the models, the $Y$ particles are
the lightest in the hidden sector and act as mediators between the heavier $X$
and the Standard Model; dark matter annihilates via $XX\rightarrow YY$.

The equations for the homogeneous energy densities of the three components
are:

$\displaystyle\dot{\rho}_{Y}+3H(1+w_{Y})\rho_{Y}$ $\displaystyle=-\Gamma
mn_{Y};$ (2.1a) $\displaystyle\dot{\rho}_{R}+4H\rho_{R}$ $\displaystyle=\Gamma
mn_{Y};$ (2.1b) $\displaystyle\dot{\rho}_{X}+3H\rho_{X}$ $\displaystyle=0,$
(2.1c)

where overdots denote $d/dt$ and $H\equiv\dot{a}/a$. In Eq. (2.1a), $n_{Y}$ is
the number density of $Y$ particles while $w_{Y}$ is their time-varying
equation of state parameter, the ratio of their pressure to their density. As
$w_{Y}$ falls from 1/3 to 0, the $Y$ particles transition from being
relativistic to nonrelativistic. The terms on the RHS of Eqs. (2.1a) and
(2.1b) depend on $mn_{Y}$ instead of $\rho_{Y}$ because the longer lab-frame
lifetimes of faster particles compensate for the higher energies released by
their decays [12]. For details of the methods used to numerically solve for
the evolution of the hidden sector temperature and the density and pressure of
the $Y$ particles, we refer the reader to G23.

Figure 1: The background evolution of the energy densities of the $Y$
particles, SM radiation ($R$) and dark matter ($X$) as a function of scale
factor. The yellow shaded region shows the EMDE. The particular scenario
exhibited here is one in which the SM radiation initially dominates
($\eta>1$); scenarios are also possible in which relativistic $Y$ particles
initially dominate ($\eta<1$).

Figure 1 shows the time evolution of the background densities of the three
fluids, obtained by solving Eqs. (2.1). Equation (2.1b) shows that the SM
radiation density drops as $\rho_{R}\propto a^{-4}$ when when $\Gamma
mn_{Y}\ll H\rho_{R}$. When $\Gamma mn_{Y}$ becomes comparable to $H\rho_{R}$,
$\rho_{R}\propto a^{-3/2}$ due to the entropy injection from the decay of the
$Y$ particles into the visible sector. After the $Y$ particles decay away,
$\rho_{R}\propto a^{-4}$ again.

At the start of our calculations, $\rho_{Y}\propto a^{-4}$ because the $Y$
particles are relativistic. As $T_{\rm hs}$ drops below the $Y$ particle mass
($m$), the $Y$ particles become nonrelativistic. This transition can be
modelled as a broken power law with a sharp pivot point at the scale factor
$a_{p}=ba_{i}T_{\rm hs,i}/m$, where $b=2.70$ for bosonic $Y$ particles and
$3.15$ for fermionic $Y$ particles [30], and the subscript $i$ indicates any
time when the $Y$ particles are relativistic. When the $Y$ particles become
nonrelativistic, $\rho_{Y}\propto a^{-3}$ until $\Gamma/H\sim 1$. When
$\rho_{Y}$ exceeds $\rho_{R}$, the universe becomes matter-dominated. This
phase is shown by the yellow shaded region in Figure 1. When $\Gamma/H$
becomes comparable to unity, the decay of the $Y$ particles begins to
significantly affect their abundance, and their comoving number density falls
rapidly. Radiation domination starts shortly thereafter. This process is
called reheating, and, although it does not happen instantaneously, it is
conventional to define the reheat temperature via the relation

$\Gamma\equiv\sqrt{\frac{8\pi
G}{3}\frac{\pi^{2}}{30}g_{*}({T}_{\mathrm{RH}}){T}_{\mathrm{RH}}^{4}}\,,$
(2.2)

which equates $\Gamma$ to what the Hubble rate would be in a purely radiation-
dominated universe at temperature ${T}_{\mathrm{RH}}$. Another quantity that
is useful to define is the scale factor ${a}_{\mathrm{RH}}$ of reheating,
given by

$\frac{{a}_{\mathrm{RH}}}{a_{0}}=\frac{1}{1.02}\left[\frac{g_{*S}(T_{0})}{g_{*S}(0.204{T}_{\mathrm{RH}})}\right]^{\frac{1}{3}}\left[\frac{T_{0}}{{T}_{\mathrm{RH}}}\right],$
(2.3)

where $a_{0}$ and $T_{0}$ are respectively the scale factor and the radiation
temperature today [30]. Note that both ${T}_{\mathrm{RH}}$ and
${a}_{\mathrm{RH}}$ are defined quantities and
$T({a}_{\mathrm{RH}})\neq{T}_{\mathrm{RH}}$.

The behavior of the homogeneous background universe is controlled by three
parameters. The ratio $\eta$ of the SM and $Y$ energy densities when the $Y$
particles are relativistic controls which component dominates the universe
before the EMDE. The $Y$ particle mass $m$ decides when the $Y$ particles
become nonrelativistic. Together, $\eta$ and $m$ determine when $\rho_{Y}$
exceeds $\rho_{R}$. For a fixed $\eta$, a higher value of $m$ means the $Y$
particles become nonrelativistic earlier and the EMDE begins earlier. For the
same value of $m$, a higher value of $\eta$ means the EMDE begins later. In
addition to these, ${T}_{\mathrm{RH}}$ determines when the EMDE ends by
setting the decay rate of the $Y$ particles — a smaller value of
${T}_{\mathrm{RH}}$ for the same $\eta$ and $m$ means that reheating happens
later and the EMDE is longer.

The case shown in Fig. 1 has $\eta=100$, which implies that SM radiation
dominates the energy budget of the universe before the EMDE. There can be
cases with $\eta<1$, in which $\rho_{Y}>\rho_{R}$ initially and the EMDE
begins when the $Y$ particles become nonrelativistic.

###  2.2 Perturbation Evolution and Power Spectra

Subhorizon matter perturbations evolve linearly with the scale factor in a
matter-dominated era, as opposed to the logarithmic evolution seen during
radiation domination [18]. As a result, in cosmologies with an EMDE, small-
scale modes which enter the horizon during or before the EMDE are enhanced in
amplitude compared to cases without an EMDE. This effect corresponds to an
enhancement of the matter power spectrum on small scales. Additionally, the
relativistic pressure of the $Y$ particles acts against gravity and suppresses
the growth of $Y$ density perturbations $\delta_{Y}$. Modes which enter the
horizon when the $Y$ particles have significant pressure experience inhibited
growth until the $Y$ particle becomes nonrelativistic and the pressure
subsides. This effect leads to a small-scale cut-off in the power spectrum of
$\delta_{Y}$.

During the EMDE, the $Y$ particles are the dominant component in the universe
and they cluster to form gravitational wells. The $X$ particles fall into
these gravitational centers. As a result, the $X$ density perturbations
$\delta_{X}$ track $\delta_{Y}$ and the dark matter ($X$) power spectrum
inherits the shape and amplitude of the $\delta_{Y}$ power spectrum, including
the enhancement due to the EMDE and the small-scale cut-off resulting from the
relativistic pressure of the $Y$ particles.

Figure 2: The evolution of perturbations for two modes, normalized to the
primordial potential and plotted as a function of $a/a_{k}$, where $a_{k}$ is
the scale factor of horizon entry for the mode. The mode in the left panel
enters the horizon after the $Y$ particles have become pressureless ($T_{\rm
hs}/m=0.00018$ at horizon entry, corresponding to equation of state parameter
$w_{Y}=0.00017$); there is no suppression of $\delta_{Y}$ for this mode. The
mode in the right panel enters the horizon when the $Y$ particles have
significant pressure ($T_{\rm hs}/m=0.21$ at horizon entry, corresponding to
$w_{Y}=0.145$) because of which $\delta_{Y}$ is suppressed compared to
$\delta_{Y,c}$, the $Y$ density perturbation if the $Y$ particles are
nonrelativistic. Dark matter perturbations $\delta_{X}$ eventually track
$\delta_{Y}$ in both cases.

Figure 2 (taken from G23) illustrates the time evolution of two perturbation
modes in Newtonian gauge, showing $\delta_{X}$ and $\delta_{Y}$ scaled by
$\Phi_{0}$, the primordial metric perturbation in a radiation-dominated
universe. It also shows the evolution of $\delta_{Y,c}/\Phi_{0}$, the $Y$
density perturbation if the $Y$ particles were pressureless. The left panel
shows the evolution of the density perturbations for a mode which enters the
horizon when the $Y$ particles are nearly pressureless, as evidenced by the
equation of state $w_{Y}$ at horizon entry being 0.00017. In this case,
$\delta_{X}$, $\delta_{Y}$ and $\delta_{Y,c}$ evolve like a cold dark matter
perturbation. The right panel shows the evolution of the same three for a mode
which enters the horizon when $w_{Y}=0.18$ and the $Y$ particles have
significant relativistic pressure. Here, the growth of $\delta_{Y}$ is
initially suppressed compared to $\delta_{Y,c}$, before the pressure subsides
and $\delta_{Y}$ starts growing linearly with scale factor during the EMDE. In
addition, $\delta_{X}$ starts tracking $\delta_{Y}$ after some time has
elapsed in the EMDE, on account of falling into the gravitational wells
created by the $Y$ particles.

Figure 3: The dimensionless power spectrum
$\mathcal{P}(k)\equiv[k^{3}/(2\pi^{2})]P(k)$ of dark matter density
perturbations, evaluated in linear theory at ${a}_{\mathrm{RH}}=1.64\times
10^{-11}$ for a cosmology with ${T}_{\mathrm{RH}}=10$ MeV, $\eta=100$
(implying that SM radiation initially dominates), and a bosonic $Y$ particle
of mass $m=1$ TeV that dominates the energy density of the universe during the
EMDE. The vertical dashed line shows ${k}_{\mathrm{RH}}$, roughly the
wavenumber corresponding to the horizon at the end of the EMDE. The dotted
line marks the wavenumber that enters the horizon at the onset of the EMDE.

The exact nature of the effect of the $Y$ particle pressure on the matter
power spectrum was worked out in G23 using a custom Boltzmann solver that
incorporated the effects of the varying equation of state and sound speed of
the $Y$ particles into the perturbation equations for the three fluid
densities and velocities. Figure 3 shows one example solution, showing the
dimensionless $X$ power spectrum $\mathcal{P}(k)\equiv[k^{3}/(2\pi^{2})]P(k)$
evaluated in linear theory:

$\mathcal{P}(k)\equiv
A_{s}\left(\frac{k}{k_{0}}\right)^{n_{s}-1}\left(\frac{4}{9}\right)\left(\frac{\delta_{X}(k)}{\Phi_{0}}\right)^{2},$
(2.4)

where $A_{s}=2.2\times 10^{-9}$, $k_{0}=0.05$ Mpc-1, and $n_{s}=0.965$ are
taken from the Planck 2015 results [31]. The first two factors in Eq. (2.4)
are the power spectrum of primordial curvature perturbations $\zeta$, and the
remaining factors translate this into the matter power spectrum (since the
primordial gravitational potential $\Phi_{0}=-2\zeta/3$). The power spectrum
is shown at ${a}_{\mathrm{RH}}=1.64\times 10^{-11}$ for a case with a bosonic
$Y$ particle of $m=1$ TeV, $\eta=100$ (so that SM radiation initially
dominates), and ${T}_{\mathrm{RH}}=10$ MeV. The vertical dashed line shows
${k}_{\mathrm{RH}}\equiv{a}_{\mathrm{RH}}\Gamma$. For modes with
$k\gtrsim{k}_{\mathrm{RH}}$, the power spectrum is enhanced due to the EMDE.
Modes that enter the horizon during the EMDE (${k}_{\mathrm{RH}}<k<k_{\rm
dom}$) grow linearly with $a$ from horizon entry until ${a}_{\mathrm{RH}}$.
For these modes,
$\delta_{X}(k,{a}_{\mathrm{RH}})\propto(k/{k}_{\mathrm{RH}})^{2}$ and
$\mathcal{P}(k)\propto(k/{k}_{\mathrm{RH}})^{4}$. Modes with $k>k_{\rm dom}$
enter the horizon in the radiation-dominated era before the EMDE and
experience logarithmic growth with $a$ before linear growth during the EMDE.
If the $Y$ particles were pressureless, this growth would imply that
$\mathcal{P}(k)\propto[\ln(k/k_{\rm dom})]^{2}$ for $k>k_{\rm dom}$. This
transition to logarithmic shape at $k_{\rm dom}$ is apparent in Fig. 3. In
contrast, for a case with initial $Y$ domination ($\eta<1$),
$\mathcal{P}(k)\propto(k/{k}_{\mathrm{RH}})^{4}$ for all $k>{k}_{\mathrm{RH}}$
down to the cut-off scale arising due to the $Y$ particle pressure.

Due to the pressure of the $Y$ particles, the growth of modes with $k\gtrsim
10^{8}$ Mpc-1 in Fig. 3 is suppressed. This suppression leads to a
characteristic peak in the power spectrum. For scales smaller than the peak
scale, the amplitude of $\mathcal{P}(k)$ displays an oscillating pattern with
a decaying envelope. This decaying amplitude is because smaller scales enter
the horizon earlier and experience a longer period of suppressed growth.

##  3 Structure Formation During The EMDE and Gravitational Heating

In the scenario shown in Fig. 3, the linearly extrapolated matter power
spectrum significantly exceeds 1 (dashed horizontal line) at reheating.
Consequently, nonlinear structures would have already formed during the EMDE.
According to Press-Schechter theory [32], the fraction $f_{\mathrm{bound}}$ of
matter in halos is

$f_{\mathrm{bound}}(a)=2\int_{\delta_{\rm
c}/\sigma}^{\infty}\mathrm{d}\nu\frac{\mathrm{e}^{-\nu^{2}/2}}{\sqrt{2\pi}}=\mathrm{erfc}\left(\frac{\delta_{\rm
c}}{\sqrt{2}\sigma(a)}\right)$ (3.1)

at the scale factor $a$, where $\sigma(a)=[\int\mathcal{P}(k,a)\mathrm{d}\ln
k]^{1/2}$ is the rms of the unfiltered linear density contrast field,
$\delta_{\rm c}\simeq 1.686$ is the linear threshold for spherical collapse,
and $\mathrm{erfc}(x)\equiv 1-\mathrm{erf}(x)$ is the complementary error
function.222In Press-Schechter theory, particles lying at a point in the
initial conditions become part of a halo at some later time if the linear
density contrast about that point exceeds the threshold, $\delta_{\rm c}$,
when smoothed on some scale. In Eq. (3.1), $f_{\mathrm{bound}}/2$ is the
fraction of points at which the unfiltered density exceeds the collapse
threshold. The factor of 2 accounts for the points for which the density only
exceeds the threshold when smoothed on a nonzero scale [33]. To quantify the
dominance of nonlinear structure formation during the EMDE, we consider
$f_{\mathrm{bound}}({a}_{\mathrm{RH}})$, the fraction of $Y$ particles bound
in halos by the time of reheating. For a case with
$\\{m,\eta,{T}_{\mathrm{RH}}\\}=\\{2\text{ TeV},10,{20}\text{ MeV}\\}$, the
bound fraction is 0.805. That is, more than 80% of the $Y$ particles are in
bound structures at reheating in this case. Since the $X$ particles fall into
the gravitational wells created by the $Y$ particles, we can assume a similar
fraction of $X$ particles are in bound halos at this time.

Bound microhalos that form before reheating are mostly composed of the $Y$
particles, since $\rho_{Y}\gg\rho_{X}$. Around the time of reheating, as the
decay of the $Y$ particles suppresses their abundance, the gravitational
support of these structures vanishes. The $X$ particles that have virialized
in a microhalo cannot adjust their orbits adiabatically to this evaporation of
the halo’s mass, so they leave the halo at boosted speeds in random
directions. This gravitational heating leads to free streaming of dark matter
after reheating, which suppresses density perturbations on comoving length
scales smaller than the free-streaming length $\lambda_{\rm fs}$. Reference
[20] estimated that structures formed after about $10^{-6}{a}_{\mathrm{RH}}$
will dissipate at reheating, releasing their dark matter particles to stream
freely.

If the population of dark matter particles has some (peculiar) velocity
dispersion $\sigma_{v}$, the free-streaming length is given by

$\lambda_{\rm
fs}(t)\sim\int^{t}_{{t}_{\mathrm{RH}}}\frac{\sigma_{v}(t^{\prime})}{a}\mathrm{d}t^{\prime}.$
(3.2)

Assuming that nonlinear growth stops at ${a}_{\mathrm{RH}}$ and the velocity
dispersion redshifts after that as
$\sigma_{v}(a)={\sigma_{v}}_{\mathrm{RH}}{a}_{\mathrm{RH}}/a$, the above
integral can be recast as

$\lambda_{\rm
fs}(a)\sim{\sigma_{v}}_{\mathrm{RH}}{a}_{\mathrm{RH}}\int^{a}_{{a}_{\mathrm{RH}}}\frac{\mathrm{d}a^{\prime}}{a^{\prime
3}H(a^{\prime})}.$ (3.3)

During the period of radiation domination that follows the EMDE, we have
$H(a)\propto a^{-2}$, where we neglect the minor influence of changes in the
relativistic content. This implies
$a^{2}H(a)={a}_{\mathrm{RH}}^{2}H({a}_{\mathrm{RH}})$, using which the
integral is

$\lambda_{\rm
fs}(a)\sim\frac{{\sigma_{v}}_{\mathrm{RH}}}{{a}_{\mathrm{RH}}H({a}_{\mathrm{RH}})}\int^{a}_{{a}_{\mathrm{RH}}}\frac{\mathrm{d}a^{\prime}}{a^{\prime}}=\frac{{\sigma_{v}}_{\mathrm{RH}}}{{a}_{\mathrm{RH}}H({a}_{\mathrm{RH}})}\ln\left[\frac{a}{{a}_{\mathrm{RH}}}\right].$
(3.4)

Finally, defining $k_{\rm fs}=\lambda_{\rm fs}^{-1}$ and using
${k}_{\mathrm{RH}}={a}_{\mathrm{RH}}\Gamma=0.96{a}_{\mathrm{RH}}H({a}_{\mathrm{RH}})$,333The
numerical factor here corresponds to
$\Gamma=H({T}_{\mathrm{RH}})=0.96H({a}_{\mathrm{RH}})$, where $H$ is the
radiation-dominated extrapolation of the Hubble rate, and follows from Eq.
(2.3). where $H({a}_{\mathrm{RH}})$ is the Hubble rate in a radiation-
dominated universe, we can write

$\frac{k_{\rm
fs}}{{k}_{\mathrm{RH}}}\sim\frac{1}{{\sigma_{v}}_{\mathrm{RH}}}\left(\ln\left[\frac{a}{{a}_{\mathrm{RH}}}\right]\right)^{-1}.$
(3.5)

Note that this discussion is conceptual, as ${\sigma_{v}}_{\mathrm{RH}}$ is
set by nonlinear gravitational evolution and cannot be evaluated directly (see
Sec. 5 for a numerically precise analysis).

This free streaming will result in a suppression of power on scales with
$k\gtrsim k_{\rm fs}$, leading to a fall-off in the matter power spectrum
after the EMDE ends. As Eq. (3.5) shows, the cut-off scale grows larger
($k_{\rm fs}$ decreases) with time, as the dark matter particles cover
increasing distances. This growth is logarithmic during the radiation epoch
but comes to a halt after matter-radiation equality, since during matter
domination, the integral $\int da/[a^{3}H(a)]$ converges in the infinite
future.

As the free-streaming scale grows larger, the enhancement to the power
spectrum resulting from the EMDE is erased up to larger scales. Yet as long as
$k_{\rm fs}/{k}_{\mathrm{RH}}>1$ around the time of matter-radiation equality,
some of the enhancement to the power spectrum due to the EMDE remains even
after the impact of gravitational heating, and microhalos can start forming
much earlier than they do in standard $\Lambda$CDM scenarios. Since the
density of a halo forming at $a_{f}$ scales as $a_{f}^{-3}$, these early-
forming structures will have cores much denser than halos which form in
standard cosmologies, greatly boosting their annihilation signals while also
increasing their susceptibility to detection by gravitational means. Accurate
modeling of the free-streaming scale and the manner in which it suppresses the
matter power spectrum is therefore important for developing observational
constraints on these EMDE scenarios.

Reference [20] previously estimated the cut-off to the power spectrum imposed
by gravitational heating by considering the virial velocities of the halos
that form during the EMDE. We improve upon this treatment by using
cosmological $N$-body simulations, which we describe next.

##  4 Simulations

To analyze the formation of nonlinear structures during the EMDE and the
impact of their subsequent evaporation on the matter power spectrum, we used a
modified version of the $N$-body simulation code GADGET-2 [34]. The code was
changed to include SM radiation as in Refs. [35, 36] by adding homogeneous
density terms to the expression for the Hubble rate in the time integrals that
are used to calculate the drifts and kicks in the simulations. The decaying
$Y$ particles are modeled in the same way at the homogeneous level, and their
inhomogeneity is accounted for by scaling the efficiency of gravitational
kicks proportionally to the time-varying quantity
$(\bar{\rho}_{Y}+\bar{\rho}_{X})/\bar{\rho}_{X}$;444More specifically, we
scale the time integrand for gravitational kicks. This approach avoids
introducing artifacts associated with the finite time-step resolution. this is
conceptually equivalent to scaling the simulation particle mass $M$ as

$M=\left(\frac{\bar{\rho}_{Y}+\bar{\rho}_{X}}{\bar{\rho}_{X}}\right)M_{0},$
(4.1)

where $M_{0}$ is the simulation particle mass at late times. This scheme was
developed for and tested in Ref. [20], and similar methods have been employed
in simulations of decaying dark matter [37, 38, 39]. Since the $X$ and $Y$
particles are coupled gravitationally during the EMDE, we represent both
species together in each simulation particle.

We used our modified code to run a suite of nine cosmological simulations. The
initial conditions in each case were set by sampling a random initial density
field from the power spectrum of $Y$ particles obtained using the Boltzmann
solver described in G23. Simulation particles were initially on a uniform
grid, and their positions and velocities were perturbed using the Zel’dovich
approximation.

We began all of our simulations during the EMDE, at $a=0.02{a}_{\mathrm{RH}}$,
with the exception of two cases (YD1 and YD2), which began at the earlier time
$a=0.01{a}_{\mathrm{RH}}$ and $0.005{a}_{\mathrm{RH}}$, respectively. The
starting times were chosen to ensure that the matter power spectrum is
initially in the linear regime, as much as possible. Periodic boundary
conditions were assumed and the comoving box sizes were set to half the size
of the comoving horizon at the starting time. The RMS density fluctation at
the scale of the mean interparticle distance,
$\sigma=[\int\mathcal{P}(k,a_{\rm start})\mathrm{d}\log k]^{1/2}$ in each box
was less than 0.15 at the starting time, with the exception of YD2, for which
it was 0.35. The comoving softening length for the simulation particles was
set to 0.03 times the initial grid spacing.

Case |  $m$ (TeV) | $\eta$ |  ${T}_{\mathrm{RH}}$ (GeV) | ${a}_{\mathrm{RH}}$ |  ${k}_{\mathrm{RH}}$ (1/Mpc) | Bound % at ${a}_{\mathrm{RH}}$
---|---|---|---|---|---|---
RD1 | 2 | 70 | 0.02 | $8.21\times 10^{-12}$ | $2.36\times 10^{5}$ | 59.66
RD2 | 2 | 10 | 0.02 | $8.21\times 10^{-12}$ | $2.36\times 10^{5}$ | 80.54
RD3 | 2 | 10 | 0.05 | $3.29\times 10^{-12}$ | $6.68\times 10^{5}$ | 36.02
YD1 | 2 | 0.01 | 0.02 | $8.21\times 10^{-12}$ | $2.36\times 10^{5}$ | 86.72
YD2 | 2 | 0.01 | 0.005 | $3.31\times 10^{-11}$ | $5.77\times 10^{4}$ | 97.23
YD3 | 2 | 0.01 | 0.069 | $2.37\times 10^{-12}$ | $9.51\times 10^{5}$ | 32.84
YD4 | 20 | 0.01 | 0.5 | $2.85\times 10^{-13}$ | $1.26\times 10^{7}$ | 27.46
RD4 | 2 | 400 | 0.02 | $8.21\times 10^{-12}$ | $2.36\times 10^{5}$ | 14.61
RD5 | 2 | 100 | 0.05 | $3.29\times 10^{-12}$ | $6.68\times 10^{5}$ | 1.81
Table 1: List of simulations with the model parameters $m$, $\eta$ and
${T}_{\mathrm{RH}}$, the values of ${a}_{\mathrm{RH}}$ and
${k}_{\mathrm{RH}}$, and the percentage of matter bound in halos at reheating,
evaluated using Eq. (3.1).

Table 1 lists the hidden-sector model parameters
$\\{m,\eta,{T}_{\mathrm{RH}}\\}$ and the details of the various simulation
runs. The three parameters were varied to cover cases with different amounts
of predicted bound structure at ${a}_{\mathrm{RH}}$ to check the impact of
differing levels of structure formation on the post-reheating cut-off that
arises due to gravitational heating. Our simulations go from a case with 1%
bound structure at reheating to a case with nearly all of the matter bound in
halos at reheating. The range of models is also associated with a variety of
different shapes of the linear matter power spectrum during the EMDE.

|
---|---
|
Figure 4: Density plots of slices of thickness $4.5\times 10^{-5}$ kpc/$h$
from the RD2 simulation box, for which ${a}_{\mathrm{RH}}=8.21\times
10^{-12}$. The box size is $9.01\times 10^{-4}$ kpc/$h$; the x-y plane is
shown from 0 to $5\times 10^{-4}$ kpc/$h$ . The four slices show the formation
and erasure of structure as the box evolves towards and past reheating.

Figure 4 shows slices of the RD2 simulation box at various times. The slices
have a thickness of 0.05 times the box size and extend across 0.45 times the
box size in the other directions; the particles are plotted as dots. The top
left panel shows the box at $0.10{a}_{\mathrm{RH}}$, when the perturbations in
the density field are small. The top right panel shows the box just after
reheating, at $1.13{a}_{\mathrm{RH}}$, displaying the formation of a web of
structure with clumps and filaments. The bottom left panel displays the box at
$11.27{a}_{\mathrm{RH}}$. By this time, the structure formed at reheating is
already being erased, with the clearly defined web visible at
$1.13{a}_{\mathrm{RH}}$ blurred and smoothed out. By $389.24{a}_{\mathrm{RH}}$
(bottom right panel), the structure formed during the EMDE has been almost
completely wiped out. The progression of the box shown via these four
snapshots illustrates the formation of structures during the EMDE and their
subsequent dissipation after reheating.

In Figure 5, we show the power spectra of the simulation particles for the
four snapshots corresponding to the four panels in Fig. 4. The panels show
$\mathcal{P}(k)$ from the simulation as dashed lines and the linear theory
$\mathcal{P}(k)$ as solid lines. For $a=0.10{a}_{\mathrm{RH}}$, the smallest
scales in the box are already nonlinear, as the power spectrum on these scales
is higher than unity (shown by the grey dotted line). At
$a=1.13{a}_{\mathrm{RH}}$, the scales with $k>10^{8}$ $h$/Mpc in the box have
become nonlinear, which corresponds to the formation of the web of structures
seen in the top right panel of 4.

At $a=11.27{a}_{\mathrm{RH}}$, as the bottom left panel of Fig. 4 shows,
erstwhile bound structures have mostly dissipated. This is shown in the fall-
off of the simulation power spectrum in the bottom left panel of Fig. 5. The
simulation power spectrum is much lower than the linear version for almost all
scales. This suppression results from the free streaming of the $X$ particles
after reheating. The bottom right panel of Fig. 5 shows that the peak of the
post-reheating power spectrum has moved to larger scales by
$389.24{a}_{\mathrm{RH}}$, confirming that the free-streaming cut-off has
moved to larger scales with time. The slight rise in the simulation power
spectrum at the highest $k$ in the bottom two panels of Fig. 5 arises due to
the Poisson noise of the discrete simulation particles. It is also apparent
from Fig. 5 that the free-streaming cut-off function is much shallower than
the Gaussian cut-off used in Ref. [20], allowing significant power to persist
more than a decade in $k$ beyond the scale at which suppression begins.

Figure 5: The dimensionless power spectrum $\mathcal{P}(k)$ taken from the RD2
simulation box (dashed lines) compared to the linear theory $\mathcal{P}(k)$
(solid lines). The times of the simulation snapshots correspond to the four
panels of Fig. 4. The grey dotted line marks $\mathcal{P}=1$.

##  5 Analyzing The Free-Streaming Cut-off

We quantify the free-streaming cut-off induced by gravitational heating as the
ratio between the simulation and linear power spectrum as a function of time
for each of our simulation runs, given by555Since $\mathcal{P}_{\rm sim}$ is
the nonlinear power spectrum, the $R_{\rm fs}$ inferred from Eq. (5.1) could
be influenced by nonlinear effects instead of only free streaming. However, we
find that this is not a significant concern. For example, in Fig. 7,
$\mathcal{P}_{\rm sim}$ peaks at 0.4 in the first panel and 0.2 in the last,
and the form of the cut-off is nevertheless essentially the same between them.

$R_{\rm fs}(k,a)\equiv\sqrt{\frac{\mathcal{P}_{\rm sim}(k,a)}{\mathcal{P}_{\rm
lin}(k,a)}},$ (5.1)

for $a>10{a}_{\mathrm{RH}}$. Here we evaluate the linear-theory power spectrum
$\mathcal{P}_{\rm lin}(k,a)$ by using the Boltzmann solver described in G23 to
obtain $\mathcal{P}_{\rm lin}(k,10{a}_{\mathrm{RH}})$ and extrapolating it to
later times using $\mathcal{P}_{\rm
lin}(k,a)\propto[\ln(a/{a}_{\mathrm{RH}})]^{2}$.

We analyzed snapshots from $11{a}_{\mathrm{RH}}$ to $790{a}_{\mathrm{RH}}$ for
each case. Figure 6 shows $R_{\rm fs}$ at different times taken from the RD5
simulation box. As $a/{a}_{\mathrm{RH}}$ increases, $R_{\rm fs}(k,a)$
decreases and the wavenumber at which $R_{\rm fs}=0.5$ (dashed line) moves to
smaller $k$, indicating that the free-streaming cut-off to the power spectrum
moves to larger scales. The noise at low $k$ (also visible in Fig. 5) is due
to the low number of modes sampled at scales close to the simulation box size
and is not related to free streaming.

Figure 6: Free-streaming cut-off $R_{\rm fs}=\sqrt{\mathcal{P}_{\rm
sim}/\mathcal{P}_{\rm lin}}$ from the RD5 simulation box, shown at different
times after reheating. The grey dashed horizontal line indicates $R_{\rm
fs}=0.5$.

For each snapshot, $R_{\rm fs}(k,a)$ was computed and fit to a function of the
form

$\frac{1}{1+(k/k_{\rm fs})^{n}},$ (5.2)

where $k_{\rm fs}$ and $n$ were treated as free parameters. We found that
$n\approx 2.5$ fits the fall-off in $R_{\rm fs}$ for the range of our cases,
so we subsequently fixed $n=2.5$ and fit this function to the simulation
$R_{\rm fs}$ while treating only $k_{\rm fs}$ as a fitting parameter. Figure 7
shows the $R_{\rm fs}$ obtained from the RD1 simulation (orange dots) at three
different times. The blue curve in each panel shows the fit to Eq. (5.2) with
$n=2.5$. The bottom panels show the relative error between the simulation-
based $R_{\rm fs}$ values and the fit function. For the range of wavenumbers
in which $R_{\rm fs}$ drops, the error stays roughly within a 10% bound,
exceeding 10% only when $R_{\rm fs}\simeq 0.25$. As in Fig. 6, the noise at
low $k$ is an artifact of the finite box size and is unrelated to free
streaming.

Figure 7: $R_{\rm fs}=\sqrt{\mathcal{P}_{\rm sim}/\mathcal{P}_{\rm lin}}$ from
the RD1 simulation box, shown at different times after reheating. The blue
curves show the function given by Eq. (5.2) with $n=2.5$, while the bottom
panels show the relative error between the $R_{\rm fs}$ from simulations and
the blue curves.

This fitting procedure yields a value of $k_{\rm fs}$ for each value of $a$.
Figure 8 shows the values of $k_{\rm fs}/{k}_{\mathrm{RH}}$ obtained by this
fitting for the RD1, RD3 and YD1 cases as red crosses, plotted against
$[\ln(a/{a}_{\mathrm{RH}})]^{-1}$. Guided by Eq. (3.5), we fit the array of
$k_{\rm fs}$ values obtained to the function

$\frac{k_{\rm fs}}{{k}_{\mathrm{RH}}}=\alpha[\ln(a/{a}_{\mathrm{RH}})]^{-1},$
(5.3)

treating $\alpha$ as a fitting parameter. The dark blue lines in Fig. 8 show
the best fit lines using the above function. A total of 18 snapshots from
$a\simeq 106{a}_{\mathrm{RH}}$ to $a\simeq 790{a}_{\mathrm{RH}}$ were used for
this linear fit, in order to avoid the transients resulting from the decay of
the $Y$ particle population and fit the free-streaming scale values only
during radiation domination. For each of our nine simulation cases, we obtain
the best fit value of $\alpha$, showing in Table 2.

Figure 8: The best fit values of $k_{\rm fs}$ for the RD1, RD3 and YD1 simulation runs, with $[\ln(a/{a}_{\mathrm{RH}})]^{-1}$ on the x-axis. The red crosses show the $k_{\rm fs}$ values obtained by fitting the $R_{\rm fs}$ from the simulations to Eq. (5.2) with $n=2.5$. The dark blue lines show the best fit of these $k_{\rm fs}$ values to a function proportional to $[\ln(a/{a}_{\mathrm{RH}})]^{-1}$.  Case |  $\sigma_{\rm lin,RH}$ (km/s) | Bound % at ${a}_{\mathrm{RH}}$ | $\alpha$ |  $\sigma_{\rm lin,RH}\alpha/10^{5}$ (km/s)
---|---|---|---|---
RD1 | 851.26 | 59.66 | 786.17 | 6.69
RD2 | 1738.77 | 80.54 | 593.51 | 10.32
RD3 | 898.54 | 36.02 | 741.59 | 6.66
YD1 | 2841.66 | 86.72 | 539.78 | 15.34
YD2 | 7276.62 | 97.23 | 554.65 | 40.36
YD3 | 1168.94 | 32.84 | 611.03 | 7.14
YD4 | 1081.34 | 27.46 | 630.61 | 6.82
RD4 | 385.80 | 14.61 | 1650.05 | 6.37
RD5 | 377.40 | 1.81 | 1636.21 | 6.18
Table 2: The proportionality constant $\alpha$ relating $k_{\rm
fs}/{k}_{\mathrm{RH}}$ and $[\ln(a/{a}_{\mathrm{RH}})]^{-1}$ for the various
simulation runs. Also listed are the linear velocity dispersion $\sigma_{\rm
lin,RH}$ and the product of linear velocity dispersion and $\alpha$. Note that
$\sigma_{\rm lin,RH}$ is calculated using Eq. (5.4) and converted to km/s from
natural units.

###  5.1 Relating the Cut-Off To The Linear Power Spectrum

According to Eq. (3.5), the value of $\alpha$ is inversely proportional to the
average velocity dispersion of the dark matter particles at reheating, denoted
by ${\sigma_{v}}_{\mathrm{RH}}$. To test this relation, we use the linear
power spectrum evaluated at ${a}_{\mathrm{RH}}$ to calculate the linear
velocity dispersion for each of our nine cases. The power spectrum at
${a}_{\mathrm{RH}}$ determines the linear velocity dispersion, which we take
to be

$\sigma_{\rm
lin,RH}\equiv{k}_{\mathrm{RH}}\left(\int_{0}^{\infty}\frac{\mathrm{d}k}{k}\frac{\mathcal{P}(k,{a}_{\mathrm{RH}})}{k^{2}}\right)^{1/2}$
(5.4)

based on the derivation given in Appendix A (we leave out the numerical factor
for simplicity). Table 2 shows the linear velocity dispersion for each case.
Also shown is the product $\sigma_{\rm lin,RH}\alpha$.

Two distinct classes can be identified from Table 2. For the six cases RD1,
RD3, YD3, YD4, RD4 and RD5, the product $\sigma_{\rm lin,RH}\alpha$ has an
average of $6.63\times 10^{5}$ km/s, with the maximum deviation of the product
from the average being 7.7% (YD3). Within this error, for these cases,
$\alpha\propto 1/\sigma_{\rm lin,RH}$. For the other three cases, the product
values are higher than those of the other six cases, and the product increases
as the predicted bound fraction at reheating increases.

The three cases with the high $\sigma_{\rm lin,RH}\alpha$ values are those
with significant structure having formed at ${a}_{\mathrm{RH}}$. Since Eq.
(5.4) is valid in the linear regime, it is not expected to hold in these
scenarios. As particles become bound into virialized structures, their
velocities stop growing, so in general the correct nonlinear velocity
dispersion is smaller than the linear extrapolation. The most massive halos
have the largest virial velocities, so we expect them to dominate the
nonlinear contributions to the velocity dispersion.

This motivates a truncated integral for the velocity dispersion, in which the
integrand of Eq. (5.4) stops contributing for modes with $k>k_{*}$, where
$\mathcal{P}(k_{*},{a}_{\mathrm{RH}})$ crosses some threshold value
$\mathcal{C}$. We define the truncated linear velocity dispersion to be

$\sigma_{\rm
t,RH}\equiv{k}_{\mathrm{RH}}\left(\int_{0}^{\infty}\frac{\mathrm{d}k}{k}\frac{\mathcal{P}_{\rm
t}(k,{a}_{\mathrm{RH}})}{k^{2}}\right)^{1/2},$ (5.5)

where

$\mathcal{P}_{\rm t}\equiv\begin{cases}\mathcal{P},&\
\mathcal{P}<\mathcal{C}\\\ \mathcal{C}^{2}/\mathcal{P},&\
\mathcal{P}>\mathcal{C}\end{cases}$ (5.6)

is the truncated linear power spectrum. The only motivation of this form is
that it suppresses the contribution of modes with $\mathcal{P}\gg\mathcal{C}$
in a continuous way and fits the simulation results well, as we will see.

Case | Bound % at ${a}_{\mathrm{RH}}$ |  $\sigma_{\rm t,RH}$ (km/s) | $\alpha$ |  $\sigma_{\rm t,RH}\alpha/10^{5}$ (km/s) |  % diff. from avg. $\sigma_{\rm t,RH}\alpha$
---|---|---|---|---|---
RD1 | 59.66 | 851.26 | 786.17 | 6.69 | 0.04
RD2 | 80.54 | 1130.70 | 593.51 | 6.71 | 0.31
RD3 | 36.02 | 898.54 | 741.59 | 6.66 | 0.40
YD1 | 86.72 | 1245.33 | 539.78 | 6.72 | 0.48
YD2 | 97.23 | 1252.44 | 554.65 | 6.95 | 3.84
YD3 | 32.84 | 1168.94 | 611.03 | 7.14 | 6.77
YD4 | 27.46 | 1081.34 | 630.61 | 6.82 | 1.93
RD4 | 14.61 | 385.80 | 1650.05 | 6.37 | 4.84
RD5 | 1.81 | 377.40 | 1636.21 | 6.18 | 7.70
Table 3: The truncated linear velocity dispersion given by Eq. (5.5) with a
threshold of $\mathcal{C}=4.3$ (c.f. Table 2). The second column from the
right lists the product of $\sigma_{\rm t,RH}$ and $\alpha$. The rightmost
column shows the percentage by which the product differs from the average
product in the nine cases. The percentage difference is within 7.7%,
demonstrating that $\sigma_{\rm t,RH}$ is inversely proportional to $\alpha$.

Various values of $\mathcal{C}$ were tried to obtain values of $\sigma_{\rm
t,RH}$ that would bring the product $\sigma_{\rm t,RH}\alpha$ for the three
high-bound-fraction cases closer to the products in the six other cases. We
found $\mathcal{C}=4.3$ to be an appropriate threshold. Table 3 shows the
truncated linear velocity dispersion and the products $\sigma_{\rm
t,RH}\alpha$ for all nine cases with $\mathcal{C}=4.3$. This truncation of the
linear velocity dispersion does not affect the six cases in which the bound
fraction at reheating is less than $80\%$, but the values of $\sigma_{\rm
t,RH}\alpha$ are close to each other for all the nine cases, with a maximum
deviation of 7.7% from the mean value.

This prescription of calculating $\sigma_{\rm t,RH}$ yields a velocity
dispersion that is inversely proportional to $\alpha$ within 7.7% error for
our nine cases, covering the wide range of scenarios, with the bound
percentages at ${a}_{\mathrm{RH}}$ ranging from 1.81% to 97.23%. For the range
of cases studied in this work, the free-streaming cut-off is given by

$\frac{k_{\rm fs}(a)}{{k}_{\mathrm{RH}}}=\frac{2.23}{\sigma_{\rm
t,RH}}\left[\ln\left(\frac{a}{{a}_{\mathrm{RH}}}\right)\right]^{-1},$ (5.7)

where the numerical factor is the average of the $\sigma_{\rm t,RH}\alpha$
values from Table 3 (converted to units of $c$).

For the RD4 and RD5 cases, since the amount of structure formed at
${a}_{\mathrm{RH}}$ is not very high, the free-streaming fall-off in $R_{\rm
fs}$ is a little shallower than a $(k/k_{\rm fs})^{-2.5}$ power law. Figure 9
shows $R_{\rm fs}$ from the RD4 simulations compared to the fitting form given
by $(1+(k_{\rm fs}/{k}_{\mathrm{RH}})^{2.5})^{-1}$. At
$66.22{a}_{\mathrm{RH}}$, the slopes of the fall-off from the simulation and
the fitting function are noticeably different. At $889.59{a}_{\mathrm{RH}}$,
the slopes are much closer to each other. The relative error between the
$R_{\rm fs}$ from the simulation and the best fit curve remains within 10% for
$R_{\rm fs}\gtrsim 0.25$. Cases with such low amounts of bound structure at
${a}_{\mathrm{RH}}$ present the boundaries down to which the $(1+(k_{\rm
fs}/{k}_{\mathrm{RH}})^{2.5})^{-1}$ fitting form is valid. For cases with
lower bound percentages at ${a}_{\mathrm{RH}}$, the fall-off in $R_{\rm fs}$
is much shallower than $(k/k_{\rm fs})^{-2.5}$; we leave the detailed study of
these cases to future work.

Figure 9: $R_{\rm fs}=\sqrt{\mathcal{P}_{\rm sim}/\mathcal{P}_{\rm lin}}$ from
the RD4 simulation box, shown at different times after reheating. The blue
curves show the function given by Eq. (5.2) with $n=2.5$, while the bottom
panels show the relative error between the $R_{\rm fs}$ from simulations and
the blue curves.

Finally, we remark that due to discreteness noise, it is difficult to resolve
the precise form of the free-streaming cut-off for $k\gg k_{\rm fs}$. This
regime can nevertheless be important, since the linear power spectrum arising
from an EMDE can grow so rapidly as a function of $k$. Another possibility is
to derive the form of the cut-off from the velocity distribution in the
simulations. We use this approach in Appendix B to obtain an alternative
expression for the cut-off function $R_{\rm fs}$ and the free-streaming scale
$k_{\rm fs}$. For the remainder of this article, however, we will continue to
use the simple cut-off expression given by Eq. (5.2) (with free-streaming
scale given by Eq. 5.7).

##  6 Structure Formation After Reheating

With the shape and scale of the free-streaming cut-off ascertained, we can
obtain the matter power spectrum in the post-EMDE era in scenarios with
gravitational heating. The power spectrum in turn determines the microhalo
distribution that arises in connection with the scenario.

We first make a technical note. The expression $\lambda_{\rm
fs}\propto\ln(a/{a}_{\mathrm{RH}})$ is derived under the assumption of
radiation domination. During the transition to late matter domination,

$\displaystyle\lambda_{\rm fs}$
$\displaystyle\propto\int_{{a}_{\mathrm{RH}}}^{a}\frac{\mathrm{d}a^{\prime}}{a^{\prime}\sqrt{a_{\rm
eq}+a^{\prime}}}$
$\displaystyle\propto\ln\left[\left(\frac{1+\sqrt{1+{a}_{\mathrm{RH}}/a_{\rm
eq}}}{1+\sqrt{1+a/a_{\rm
eq}}}\right)^{2}\frac{a}{{a}_{\mathrm{RH}}}\right]\simeq\ln\left[\left(\frac{2}{1+\sqrt{1+a/a_{\rm
eq}}}\right)^{2}\frac{a}{{a}_{\mathrm{RH}}}\right]$ (6.1)

(compare Eq. 3.4), where we assume $a_{\rm eq}\gg{a}_{\mathrm{RH}}$ and again
neglect changes to the radiation content. This expression is approximately
$\ln(a/{a}_{\mathrm{RH}})$ when $a\ll a_{\rm eq}$ but asymptotes to
$\ln(4a_{\rm eq}/{a}_{\mathrm{RH}})$ in the $a\gg a_{\rm eq}$ limit. We
therefore evaluate the free-streaming cutoff arising from gravitational
heating by replacing $\ln(a/{a}_{\mathrm{RH}})$ with the logarithm in Eq. (6)
in the expression for the cut-off scale, Eq. (5.7).

###  6.1 The linear matter power spectrum

Figure 10 shows the $\Lambda$CDM power spectrum $\mathcal{P}$ at $a_{\rm eq}$
(black solid curve). We first evaluated the linear power spectrum using CAMB
at $z=500$ and translated it back in time to $a_{\rm eq}$ via the growth
factor derived in Appendix C. This growth factor ignores the scale-dependent
evolution of modes that enter the horizon after reheating, but the amplitudes
of these modes are not affected by the EMDE, so they are largely irrelevant
when evaluating how the EMDE influences structure formation. The dashed
colored curves show the power spectra without the effects of free-streaming
for three EMDE cosmologies, computed by multiplying the $\Lambda$CDM power
spectrum by the transfer functions from G23. The orange and green dashed
curves are different only in the value of the $Y$ particle mass $m$. The shape
of the power spectrum peak is similar for these two cases, since the value of
$\eta$ is the same. The case with the lower $m$ (orange) has a larger cut-off
scale from the pressure of the $Y$ particles compared to the green power
spectrum. In contrast, the blue dashed curve is associated with the same $m$
and ${T}_{\mathrm{RH}}$ as the green dashed curve, but its value of $\eta$ is
much higher. This difference leads to a much shallower peak for the blue
dashed curve.

The solid curves show the power spectra for the same EMDE cases with the free-
streaming cut-off included by multiplying the dashed curves by the cut-off
function $R_{\rm fs}^{2}$ from Eq. (5.2) with cut-off scale $k_{\rm fs}$ given
by Eq. (5.7) (modified as described after Eq. 6) at $a=a_{\rm eq}$. In all
three cases, the power spectrum peaks move to much larger scales and smaller
amplitudes due to the effects of gravitational heating. Nonetheless, a
significant portion of the EMDE enhancement is retained in all cases. For
example, the case with the green curves has a 65% predicted bound fraction at
reheating, and yet the power spectrum peak after including free-streaming
effects is roughly four orders of magnitude larger than the $\Lambda$CDM
counterpart at the same scale.

Figure 10: The $\Lambda$CDM power spectrum (black) along with three EMDE-
enhanced power spectra (dashed) and the same power spectra with the free-
streaming cut-off (solid), all evaluated at $a=a_{\rm eq}$.

It is noteworthy that with gravitational heating accounted for, all three EMDE
power spectra in Fig. 10 peak around $\mathcal{P}\sim 1$ at matter-radiation
equality. In fact, this is generally expected as long as the dimensionless
linear power spectrum during the EMDE scales more steeply than
$\mathcal{P}\sim k^{2}$, because then the velocity dispersion integral in Eq.
(5.4) is dominated by high-$k$ contributions. That is, for such power spectra,
modes that are just becoming nonlinear at reheating (as opposed to modes that
are still linear) contribute the most to the velocity dispersion. This leads
to $\sigma_{\rm t,RH}\sim{k}_{\mathrm{RH}}/k_{*}$, where $k_{*}$ is the
nonlinear scale (where $\mathcal{P}\sim 1$) at reheating. Equation (5.7) then
implies that $k_{\rm fs}\sim k_{*}/\ln(a/{a}_{\mathrm{RH}})$. Shortly after
reheating, the power spectrum thus peaks at $\mathcal{P}\sim 1$ (see also Fig.
5), since $k_{\rm fs}$ is close to the nonlinear scale $k_{*}$. Over time,
$k_{\rm fs}$ shifts to logarithmically larger scales, and the power spectrum
also grows logarithmically (per Eq. C.3). The former effect tends to win,
causing the peak in the power spectrum to slightly drop over time (as is
evident in Fig. 5). For example, if $\mathcal{P}(k)=(k/k_{*})^{n}$ at
reheating, then the combination of linear growth and continued free streaming
results in the power spectrum peaking at a value around
$\mathcal{P}\sim[\ln(a/{a}_{\mathrm{RH}})]^{2}(k_{\rm
fs}/k_{*})^{n}\sim[\ln(a/{a}_{\mathrm{RH}})]^{2-n}$ at later times during
radiation domination, which drops over time as long as $n>2$. But since this
drop is only logarithmic, we expect that the power spectrum still peaks at
some value $\mathcal{P}\sim 1$ at the onset of late matter domination.

Nevertheless, Fig. 10 shows a tendency where if the power spectrum without the
gravitational heating suppression (dashed curves) peaks at a higher value, the
peak with the suppression accounted for (solid curves) is lower. Figure 11
demonstrates this trend more broadly. Here we plot the logarithm of
$\sqrt{\mathcal{P}(k_{\rm pk},a_{\rm eq})/\mathcal{P}({k}_{\mathrm{RH}},a_{\rm
eq})}$, the factor by which density variations at the scale of the peak in the
power spectrum exceed density variations at ${k}_{\mathrm{RH}}$ (which are not
affected by the EMDE). This ratio quantifies the degree to which the EMDE
boosts the power spectrum, and we plot its value with gravitational heating
accounted for as a function of its value with gravitational heating neglected.
Generally, as the unsuppressed value of $\sqrt{\mathcal{P}(k_{\rm pk},a_{\rm
eq})/\mathcal{P}({k}_{\mathrm{RH}},a_{\rm eq})}$ rises from $\sim 10^{3.5}$,
its gravitational heating-suppressed value drops from $\sim 10^{2.4}$.

More precisely, each of the three curves in Fig. 11 represents a family of
EMDE models where one parameter out of $m$, $\eta$ and ${T}_{\mathrm{RH}}$ is
varied while the others are held fixed. The red and green curves are obtained
by varying $\eta$; as indicated by the comparison between the blue and green
dashed curves in Fig. 10, smaller values of $\eta$ raise the peak value of
$\mathcal{P}$ while shifting it to larger scales. Both of these effects raise
the velocity dispersion during the EMDE, since it involves an integral over
$k^{-2}\mathcal{P}(k)$ (see Eq. 5.4), which explains why the gravitational
heating-suppressed value of $\sqrt{\mathcal{P}(k_{\rm pk},a_{\rm
eq})/\mathcal{P}({k}_{\mathrm{RH}},a_{\rm eq})}$ in Fig. 11 drops sharply as a
function of its unsuppressed value. In contrast, the purple curve is obtained
by varying the $Y$ particle mass $m$; as illustrated by the orange and green
dashed curves in Fig. 10, when $\eta$ is large, increasing $m$ raises the peak
value of $\mathcal{P}$ while shifting it to smaller (not larger) scales. These
two changes to the power spectrum influence the velocity dispersion during the
EMDE oppositely – increasing $\mathcal{P}$ raises it while moving the peak to
smaller scales lowers it – which is why the purple curve in Fig. 11 drops more
shallowly.

Figure 11: The logarithm of $\sqrt{\mathcal{P}(k_{\rm pk},a_{\rm
eq})/\mathcal{P}({k}_{\mathrm{RH}},a_{\rm eq})}$, which shows the peak
enhancement to the power spectrum. The x-axis shows the peak with the effects
of gravitational heating ignored, while the y-axis shows the peak with the
free-streaming cut-off included for the same EMDE cosmology. The red and green
curves are obtained by varying $\eta$ for different values of $m$, while the
purple curve corresponds to varying $m$ with fixed $\eta$.

Figure 11 suggests that the free-streaming-suppressed value of
$\sqrt{\mathcal{P}(k_{\rm pk},a_{\rm eq})/\mathcal{P}({k}_{\mathrm{RH}},a_{\rm
eq})}$ does not drop below a factor of $\sim 10^{1.6}$. The presence of such a
limit is explained because it corresponds to the scenario where the
unsuppressed linear power spectrum at reheating is $\mathcal{P}(k)\simeq
0.15A_{s}(k/{k}_{\mathrm{RH}})^{4}$ down to arbitrarily small scales (high
$k$; compare Fig. 3). Here $A_{s}$ is the primordial spectral amplitude; we
neglect the spectral tilt for simplicity. Since the truncated linear velocity
dispersion relevant to gravitational heating, given by Eq. (5.5), involves an
integral only over scales that are linear or mildly nonlinear
($\mathcal{P}\lesssim 1$), it converges to a finite value that evaluates to
$\sigma_{\mathrm{t,RH}}\simeq 0.73A_{s}^{1/4}\simeq 0.005$ (1500 km/s),
assuming $A_{s}=2.2\times 10^{-9}$. The free-streaming scale at matter-
radiation equality is then $k_{\rm
fs}\sim\sigma_{\mathrm{t,RH}}^{-1}[\ln(a_{\rm
eq}/{a}_{\mathrm{RH}})]^{-1}{k}_{\mathrm{RH}}$, and thus
$\sqrt{\mathcal{P}(k_{\rm pk},a_{\rm eq})/\mathcal{P}({k}_{\mathrm{RH}},a_{\rm
eq})}\sim(k_{\rm
fs}/{k}_{\mathrm{RH}})^{2}\sim\sigma_{\mathrm{t,RH}}^{-2}[\ln(a_{\rm
eq}/{a}_{\mathrm{RH}})]^{-2}$. Since $\sigma_{\mathrm{t,RH}}\simeq 0.005$ and
$\ln(a_{\rm eq}/{a}_{\mathrm{RH}})\sim\mathcal{O}(10)$, small-scale density
variations remain strongly enhanced even in this limiting scenario.

###  6.2 Dark matter halo formation

To illustrate how the gravitational heating cut-off affects halo formation, we
consider the halo mass function $\mathrm{d}n/\mathrm{d}\ln M$, the number
density of bound halos in logarithmic bins of mass, at different times around
the epoch of matter-radiation equality. According to Press-Schechter theory
[32],

$\frac{\mathrm{d}n}{\mathrm{d}\ln
M}=\sqrt{\frac{2}{\pi}}\frac{\bar{\rho}_{X}}{M}\left|\frac{\mathrm{d}\ln\sigma}{\mathrm{d}\ln
M}\right|\frac{\delta_{c}}{\sigma(M)}\exp\left[-\frac{\delta_{c}^{2}}{2\sigma^{2}(M)}\right],$
(6.2)

where $\sigma(M)$ is the rms density contrast smoothed on the mass scale $M$
and $\bar{\rho}_{X}$ is the comoving density of dark matter. Following common
practice for cosmologies with a small-scale cut-off to the power spectrum, we
use a sharp $k$-space smoothing filter and evaluate
$\sigma^{2}(M)=\int_{0}^{k_{M}}\mathcal{P}(k)\mathrm{d}k/k$ with the smoothing
wavenumber $k_{M}$ defined by $M\equiv(4\pi/3)\bar{\rho}_{X}(2.5/k_{M})^{3}$;
the factor of 2.5 is motivated in Ref. [40] by a match to simulations, since
there is no natural relationship between length and mass scales for this
smoothing filter. Since we will focus on times close to matter-radiation
equality, we use the time-dependent fitting function given in Ref. [20] for
the spherical collapse threshold, $\delta_{c}$, which is valid for mixed
matter-radiation domination.

Figure 12 shows the mass function for the case $m=2$ TeV, $\eta=50$ and
${T}_{\mathrm{RH}}=20$ MeV, the same cosmology as the green curves in Figure
10. The solid lines show the mass function with the gravitational heating cut-
off imposed on the power spectrum using the results of Sec 5. The dashed lines
show the version without the cut-off, with the EMDE enhancement intact. The
dashed lines are significantly higher than their solid counterparts for low
masses, showing the formation of a significantly higher number of halos at
small scales. This phenomenon is strongly inhibited by the effects of heating.
The solid lines show a high density of bound halos of sub-Earth mass, even as
early as $0.1a_{\rm eq}$. For $a<a_{\rm eq}$, the mass function shown by the
solid lines grows as more small-scale halos form. For $a>a_{\rm eq}$, the same
mass function decreases with time for $M\lesssim 10^{-8}M_{\odot}$ and
increases for larger masses as smaller halos merge to make larger ones. For
comparison, the mass function at the scales shown in Figure 12 is lower than
$10^{-150}\rm{Mpc}^{-3}$ at $a_{\rm eq}$ if a $\Lambda$CDM power spectrum with
a Gaussian cut-off at $10^{6}$ Mpc-1 is used. The cut-off wavenumber imposed
by gravitational heating is $k_{\rm fs}\approx 40{k}_{\mathrm{RH}}$ for this
EMDE scenario, which shows that even with a small portion of the EMDE-induced
bump persisting, the number of small-scale halos is significantly higher than
in cases without an EMDE.

Figure 12: The halo mass function, or differential number density of halos per
logarithmic mass interval, as a function of mass at different times around
matter-radiation equality for an EMDE scenario with $m=2$ TeV, $\eta=50$,
${T}_{\mathrm{RH}}=20$ MeV. The solid curves show the mass function evaluated
with the free-streaming cut-off, while the dashed curves show the version
calculated without the cut-off, with the EMDE enhancement intact.  Figure 13:
The bound fraction calculated using Press-Schechter theory at different times
after matter-radiation equality for three EMDE scenarios, which are the same
ones considered in Fig. 10. The amount of small-scale enhancement to the power
spectrum due to the EMDE is quantified by $\sigma_{\rm eq}$. The dashed lines
show the bound fractions if no free-streaming cut-off is imposed on the power
spectrum after the EMDE, while the solid lines show the bound fraction with
the cut-off. For comparison, we also show a standard $\Lambda$CDM scenario
(black curve), where we impose a Gaussian cut-off at $10^{6}$ Mpc-1
(corresponding to the smallest microhalos being roughly earth-mass; e.g. Ref.
[41]), although the result is only very weakly sensitive to this choice.

In cases with higher enhancement on small scales due to the EMDE, more
structure forms in the EMDE and the free-streaming scale at matter-radiation
equality is larger. Therefore, the higher the EMDE enhancement, the greater
the suppression of structure formation after the EMDE due to gravitational
heating. Figure 13 illustrates this point. The fraction of dark matter in
bound halos (Eq. 3.1) is plotted versus $a/a_{\rm eq}$ for the three cases in
Fig. 10, which we now label using $\sigma_{\rm eq}$, the value that the rms
linear density fluctuation would reach at $a=a_{\rm eq}$ if gravitational
heating were neglected. That is,

$\sigma_{\rm eq}\equiv\left[\int\mathcal{P}_{\rm emde}(k,a_{\rm
eq})\mathrm{d}\ln k\right]^{1/2},$ (6.3)

where $\mathcal{P}_{\rm emde}$ is the power spectrum without the free-
streaming cut-off (corresponding to the dashed curves in Fig. 10).
$\sigma_{\rm eq}$ thus quantifies the extent of the enhanced perturbation
growth during an EMDE.

The dashed curves show the predicted bound fraction for each case with
gravitational heating neglected. The bound fraction predicted under this
assumption at any given time is always larger with higher $\sigma_{\rm eq}$,
and in all three cases nearly all of the matter is predicted to be bound in
halos by the time of matter-radiation equality. The solid curves show the same
cases with the free-streaming cut-off from gravitational heating now imposed.
For higher $\sigma_{\rm eq}$, the bound fraction shown by the solid curves at
a given time is lower. This inversion shows the effect of gravitational
heating: more structure formation during the EMDE (higher $\sigma_{\rm eq}$)
leads to a larger free-streaming scale and hence more suppression of structure
formation after the EMDE. For comparison, the solid black line shows the
predicted bound fraction without the EMDE; here a $\Lambda$CDM power spectrum
is assumed with a Gaussian cut-off at $10^{6}$ Mpc-1, which corresponds to the
smallest microhalos being earth-mass (see [41] for an example). Even with the
gravitational heating-induced suppression, the bound fraction is much higher
at $a\approx 10a_{\rm eq}$ in cases with an EMDE.

Figure 14 shows, as a function of $\sigma_{\rm eq}$, the time $a_{10}$ at
which 10% of the dark matter is predicted to be bound in halos after matter-
radiation equality. By Eq. (3.1), this happens when the rms linear density
contrast is $\sigma=0.608\delta_{\rm c}$ ($\sigma=1.025$ during matter
domination) when gravitational heating is accounted for. The three colors here
correspond to the same three families of EMDE model considered in Fig. 11,
each obtained by changing one of $m$, $\eta$ and ${T}_{\mathrm{RH}}$ while
keeping the others fixed. The solid curves account for the free-streaming cut-
off imposed by gravitational heating and mark the scenarios where at least 1
percent of the matter is bound into halos by reheating (corresponding to
$\sigma_{\rm eq}\gtrsim 40$). In this regime, increasing $\sigma_{\rm eq}$
raises the bound fraction at reheating, leading to significant gravitational
heating and the consequent erasure of some of the EMDE-induced boost to the
power spectrum. This erasure causes a reduction in the peak amplitude of the
power spectrum (as we showed previously in Fig. 11), and this results in later
formation of nonlinear structure after reheating. The same trends are apparent
in Fig. 14 as in Fig. 11; varying $\eta$ causes $a_{10}$ to change more mildly
as a function of $\sigma_{\rm eq}$, compared to varying $m$, due to how the
scale of the peak in $\mathcal{P}(k)$ shifts as a result of these changes. As
discussed above, the EMDE-induced enhancement to small-scale density
variations remains large even in scenarios with an arbitrarily large degree of
structure formation during the EMDE. That outcome is reflected in Fig. 14 in
how $a_{10}$ remains below about $3a_{\rm eq}$ (redshift 1000) in all cases.

Figure 14:  The scale factor at which 10% of the dark matter is predicted to
be bound in halos according to Press-Schechter theory, plotted against
$\sigma_{\rm eq}$ (Eq. 6.3), which quantifies the extent to which the small-
scale power spectrum is enhanced during an EMDE. The curves are obtained by
varying $\eta$ (red and green) and $m$ (purple) to vary $\sigma_{\rm eq}$ by
varying the shape and peak of the EMDE-enhanced power spectrum. The dashed
curves show cases in which the bound fraction at ${a}_{\mathrm{RH}}$ $<0.01$,
the regime in which our free-streaming cut-off prescription is untested. The
dotted line shows $a_{10}$ for the same EMDE cases without the free-streaming
cut-off, which only depends on $\sigma_{\rm eq}$ regardless of the power
spectrum shape.

The dotted curve in Fig. 14 shows how $a_{10}$ varies with $\sigma_{\rm eq}$
when gravitational heating is neglected. Under this assumption, $\sigma_{\rm
eq}$ completely determines the unsmoothed rms density contrast and hence the
bound fraction during the late matter era (see Eq. 3.1), so there is only one
curve. Higher $\sigma_{\rm eq}$ leads to earlier halo formation, and when
$\sigma_{\rm eq}\gtrsim 1$, a significant halo population forms during the
radiation era,666But deviations from spherical symmetry are particularly
important for halo formation during radiation domination [42, 43, 44], so our
use of the spherical collapse threshold may yield inaccurate results in this
regime. but only as long as no significant nonlinear structures arise during
the EMDE. Further work is needed to understand how the dotted curve (with no
gravitational heating) transitions into the solid curves (with gravitational
heating) in the $1\lesssim\sigma_{\rm eq}\lesssim 40$ regime. As an
illustration, the thin dashed curves in Fig. 14 use our gravitational heating
prescription to extend the solid curves to include cases with less than 1% of
the matter bound in halos at reheating. However, the prescription is untested
in this regime, and indeed it is unclear whether the underlying principles
should continue to hold. If density variations remain linear during the EMDE,
then even though particle velocities are boosted, they remain associated with
coherent fluid motions rather than the velocity dispersion that is required to
produce a free-streaming cut-off. Nevertheless, the qualitative outcome should
hold that higher $\sigma_{\rm eq}$ results in earlier halo formation until
$\sigma_{\rm eq}\sim\mathcal{O}(10)$, beyond which point the trend should
reverse.

##  7 Summary And Discussion

In cosmologies with an early matter-dominated era (EMDE) prior to Big Bang
Nucleosynthesis, the small-scale power spectrum is enhanced due to the linear
growth of subhorizon matter perturbations with scale factor [18]. This
enhancement causes the formation of dense microhalos much earlier than bound
structures form in non-EMDE scenarios [15, 16, 18]. These new sub-Earth-mass
structures boost the dark matter annihilation signal [18, 20, 21] and can be
detected gravitationally via caustic microlensing [25, 26, 27, 28] and pulsar
timing arrays [23, 28, 24]. However, if perturbations are sufficiently
boosted, microhalos can form during the EMDE itself. They would dissipate when
the EMDE ends, releasing their dark matter particles at boosted speeds in
random directions [20]. This gravitational heating imprints a free-streaming
cut-off on the matter power spectrum after the EMDE, which suppresses some of
the EMDE-induced enhancement to the power spectrum. Ascertaining the nature of
this free-streaming cut-off is thus essential for understanding the halo
populations that arise from these EMDE scenarios and hence their observational
constraints and prospects.

To analyze the shape and scale of the free-streaming cut-off, we employed a
modified version of GADGET-2 to run the first suite of $N$-body simulations of
the process of microhalo formation during the EMDE and their evaporation at
the end of the EMDE. We found the formation of a rich web of small-scale
structure during the EMDE, which gets wiped out as structures evaporate when
the particle driving the EMDE decays and radiation domination starts. By
analyzing the power spectra of our simulations during the subsequent radiation
era, we found a free-streaming cut-off that evolved to larger scales as time
elapsed. The various simulation runs adopted different EMDE model parameters
in order to cover a range of scenarios, including variations in the level of
structure formation during the EMDE.

We found that free streaming due to gravitational heating suppresses the
matter power spectrum $\mathcal{P}(k)$ during the subsequent radiation epoch
by a factor that is fit well by the function $[1+(k/k_{\rm fs})^{2.5}]^{-2}$
(or an alternative form given in Appendix B). The free-streaming wavenumber
$k_{\rm fs}$ can be straightforwardly calculated by means of Eq. (5.7); it is
tightly connected to an expression given by Eq. (5.5) for the velocity
dispersion during the EMDE. Given the linear power spectrum due to an EMDE
(expressions for which are provided in Ref. [30]), these results enable easy
analytic computation of the impact of gravitational heating.

Finally, we used our results to explore how gravitational heating impacts
structure formation during the late matter epoch. Free streaming erases a
large portion of the EMDE-induced boost to the linear matter power spectrum;
its peak value at the time of matter-radiation equality, $a_{\rm eq}$, can be
4-5 orders of magnitude smaller than its peak value at the end of the EMDE.
Despite this, sufficient power typically persists to form a significant
microhalo population around $a_{\rm eq}$. In cases where $\sim 1$ percent of
the matter is bound in halos by the end of the EMDE, about 10 percent of the
dark matter is expected to be in microhalos already as early as $a=0.1a_{\rm
eq}$. Without an EMDE, the formation of a similar amount of bound structure
would not occur until $a\sim 100a_{\rm eq}$ ($z\sim 30$). As the level of
structure formation during the EMDE is raised, the free-streaming scale due to
gravitational heating also grows, which results in greater erasure of the
EMDE-induced boost to the power spectrum and hence later microhalo formation
times after the EMDE. However, even in cases where most of the matter is bound
in halos by the end of the EMDE, microhalo formation afterward occurs at most
marginally later than $a_{\rm eq}$.

Our prescription for the impact of gravitational heating will allow
explorations of the observational accessibility of EMDE scenarios to be
extended to a far wider range of parameters. For example, a significant
portion of the parameter space that Ref. [30] considered in relation to
observational probes is impacted by gravitational heating. For scenarios in
which the maximum enhancement of density variations would be $\sim 10^{4.4}$,
which Ref. [30] took to be the regime where gravitational heating begins to be
important, we found that the enhancement is already suppressed by free
streaming down to the $\sim 10^{2}$ level (see Fig. 11). That outcome makes
these scenarios compatible with the Fermi collaboration’s measurement of the
isotropic gamma-ray background [45], whereas without gravitational heating
they would have been ruled out. Gravitational heating is of great relevance to
the observational prospects of EMDE scenarios, both via dark matter
annihilation and through gravitational detection of microhalos, but we leave a
fuller analysis for future work.

A limitation to our results is that we only explored scenarios in which the
fraction of matter bound into halos by the end of the EMDE exceeds about 1%.
Also, for bound fractions around the 10% level, the form of the cut-off in the
power spectrum due to gravitational heating can already begin to differ from
the parametrization that is valid for high-bound-fraction cases (see Fig. 9).
Further work is needed to establish the impact of gravitational heating in
cases with only marginal structure formation during the EMDE. We remark that
this is also the regime in which microhalo formation is expected to begin deep
in the radiation-dominated epoch, and further work is separately needed to
understand the halos that arise in this way [44].

##  Appendix A Calculating the linear velocity dispersion

This Appendix shows how the velocity dispersion is connected to the power
spectrum at linear order in perturbations. Let $\bm{x}$ be the comoving
Lagrangian coordinate. If the linear density contrast field is
$\delta(\bm{x})$, the comoving displacement field $\bm{s}(\bm{x})$ is related
to $\delta(\bm{x})$ by $\nabla\cdot\bm{s}=-\delta$ at linear order. In Fourier
space, this implies

$\bm{s}(\bm{k})=\frac{\mathrm{i}\bm{k}}{\bm{k}^{2}}\delta(\bm{k})$ (A.1)

if $\bm{s}$ is irrotational. During matter domination, the comoving velocity
field is

$\dot{\bm{s}}(\bm{k})=H\bm{s}(\bm{k})=H\frac{\mathrm{i}\bm{k}}{\bm{k}^{2}}\delta(\bm{k})$
(A.2)

by the Zel’dovich approximation, where $H$ is the Hubble rate. Using the
inverse Fourier transform
$\bm{s}(\bm{x})\equiv(2\pi)^{-3}\int\mathrm{d}^{3}\bm{k}\mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{x}}\bm{s}(\bm{k})$,
we write the squared comoving velocity dispersion as

$\displaystyle\langle[\dot{\bm{s}}(\bm{x})]^{2}\rangle$
$\displaystyle=H^{2}\int\frac{\mathrm{d}^{3}\bm{k}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}\bm{k}^{\prime}}{(2\pi)^{3}}\mathrm{e}^{-\mathrm{i}(\bm{k}-\bm{k}^{\prime})\cdot\bm{x}}\langle\bm{s}^{*}(\bm{k})\cdot\bm{s}(\bm{k}^{\prime})\rangle$
$\displaystyle=H^{2}\int\frac{\mathrm{d}^{3}\bm{k}}{(2\pi)^{3}}\int\frac{\mathrm{d}^{3}\bm{k}^{\prime}}{(2\pi)^{3}}\mathrm{e}^{-\mathrm{i}(\bm{k}-\bm{k}^{\prime})\cdot\bm{x}}\frac{\bm{k}\cdot\bm{k}^{\prime}}{\bm{k}^{2}\bm{k}^{\prime
2}}\langle\delta^{*}(\bm{k})\delta(\bm{k}^{\prime})\rangle.$ (A.3)

But
$\langle\delta^{*}(\bm{k})\delta(\bm{k}^{\prime})\rangle=(2\pi)^{3}\delta_{\mathrm{D}}^{3}(\bm{k}-\bm{k}^{\prime})P(k)$,
where $P(k)$ is the power spectrum of $\delta$ and $\delta_{\mathrm{D}}^{3}$
is the three-dimensional Dirac delta function. Thus,

$\displaystyle\langle[\dot{\bm{s}}(\bm{x})]^{2}\rangle$
$\displaystyle=H^{2}\int\frac{\mathrm{d}^{3}\bm{k}}{(2\pi)^{3}}\frac{P(k)}{k^{2}}=H^{2}\int_{0}^{\infty}\frac{\mathrm{d}k}{k}\frac{\mathcal{P}(k)}{k^{2}},$
(A.4)

where $\mathcal{P}(k)\equiv k^{3}P(k)/(2\pi^{2})$ is the dimensionless power
spectrum. The peculiar velocity dispersion is then

$\sigma_{\rm
lin}=a\langle[\dot{\bm{s}}(\bm{x})]^{2}\rangle^{1/2}=aH\left(\int_{0}^{\infty}\frac{\mathrm{d}k}{k}\frac{\mathcal{P}(k)}{k^{2}}\right)^{1/2},$
(A.5)

and at reheating it is

$\sigma_{\rm
lin}|_{a={a}_{\mathrm{RH}}}=1.04{k}_{\mathrm{RH}}\left(\int_{0}^{\infty}\frac{\mathrm{d}k}{k}\frac{\mathcal{P}(k,{a}_{\mathrm{RH}})}{k^{2}}\right)^{1/2},$
(A.6)

since
$(aH)|_{a={a}_{\mathrm{RH}}}=1.04{a}_{\mathrm{RH}}\Gamma=1.04{k}_{\mathrm{RH}}$.
The numerical factor here follows from the separate definitions of
${T}_{\mathrm{RH}}$ and ${a}_{\mathrm{RH}}$; see Eq. (2.3) and footnote 3.

##  Appendix B Predicting the Free-Streaming Cut-off Using The Velocity
Distribution

In Section 5, we presented a prescription for the free-streaming cut-off that
arises from gravitational heating. That prescription, given by Eqs. (5.2) and
(5.7), was motivated solely by the ratio
$\mathcal{P}_{\mathrm{sim}}(k)/\mathcal{P}_{\mathrm{lin}}(k)$ of the matter
power spectrum in the simulations to that evaluated using linear theory alone.
Here we present an alternative prescription, in which we additionally employ
information about the velocity distribution in the simulations.

As long as particles cover distances proportional to their streaming
velocities, the distribution of particle displacements due to free streaming
is just a rescaled version of the velocity distribution. Consequently, the
impact of free streaming is to convolve the density field with this rescaled
version of the velocity distribution. This means that – up to a rescaling in
$k$ – the cut-off function $R_{\rm fs}$ is the Fourier transform of the
velocity distribution.

A more explicit way to see this is to start with the collisionless Boltzmann
equation without gravity,

$\frac{\mathrm{d}f}{\mathrm{d}t}=\frac{\partial f}{\partial
t}+v_{i}\frac{\partial f}{\partial x_{i}}=0,$ (B.1)

where $f(\bm{x},\bm{v},t)$ is the distribution function and the repeated index
is summed over. We neglect cosmic expansion for simplicity, but it is easy to
check that the conclusion is unchanged if expansion is accounted for. In
Fourier space

$\partial f/\partial t+\mathrm{i}\bm{v}\cdot\bm{k}f=0,$ (B.2)

the solution to which is

$f(\bm{k},\bm{v},t)=f(\bm{k},\bm{v},0)\mathrm{e}^{-\mathrm{i}\bm{k}\cdot\bm{v}t}$
(B.3)

for some initial time $t=0$. Now approximate that the distribution function is
separable, $f(\bm{k},\bm{v},t)=\rho(\bm{k},a)f_{v}(\bm{v})$, i.e., the
velocity distribution $f_{v}$ is independent of position. Integrating Eq.
(B.3) over velocities then yields

$\displaystyle\rho(\bm{k},t)$
$\displaystyle=\rho(\bm{k},0)\int\mathrm{d}^{3}\bm{v}f_{v}(\bm{v})\mathrm{e}^{-\mathrm{i}\bm{k}\cdot\bm{v}t}=\rho(\bm{k},0)\tilde{f}_{v}(\bm{k}t),$
(B.4)

where $\tilde{f}_{v}$ is the Fourier transform of $f_{v}$. The $k$-space
density $\rho(\bm{k},t)$ is evidently suppressed, in relation to
$\rho(\bm{k},0)$, by the Fourier transform of the velocity distribution,
scaled in wavenumber by a time-dependent factor. Accounting for cosmic
expansion changes the time dependence but otherwise leaves the result
unchanged.

The velocity distribution in the RD2 simulation is shown in the left panel of
Fig. 15. The distribution of normalized velocities $\bm{v}/\sigma$, where
$\sigma=\langle\bm{v}^{2}\rangle^{1/2}$ is the velocity dispersion, converges
quickly after reheating to a form that is well approximated by

$f(\bm{x})=3.98\times
10^{-11}\left[\mathrm{e}^{(x/14.9)^{1.94}}-0.953\right]^{-7.06}$ (B.5)

(dotted curve), where $\bm{x}=\bm{v}/\sigma$. This distribution has a
significantly more pronounced low-velocity tail than the Maxwell-Boltzmann
distribution (faint dashed curve).

Figure 15: The velocity distribution (left) and its Fourier transform (right)
in the RD2 simulation at several times (different colors). We show the
spherical average and normalize to the velocity dispersion
$\sigma=\langle\bm{v}^{2}\rangle^{1/2}$. The distribution evidently does not
vary after about $3{a}_{\mathrm{RH}}$. We show as black dotted curves the fit
to the late-time velocity distribution, given by Eq. (B.5), in the left panel,
and a fit to its Fourier transform, Eq. (B.6), in the right panel. For
comparison, the gray dashed curves show the Maxwell-Boltzmann distribution.

The right panels of Fig. 15 show the three-dimensional Fourier transforms of
the velocity distributions in the left panels. At late times, the Fourier
transformed distribution is well approximated by

$\tilde{f}(\bm{x})=\left[2.43\mathrm{e}^{(x/3.69)^{1.99}}-1.43\right]^{-2.81}$
(B.6)

(dotted curve). The strong low-velocity support in Eq. (B.5) evidently leads
to a long tail in the Fourier transform; $\tilde{f}$ is significantly larger
for $x\gg 1$ than the Fourier transformed Maxwell-Boltzmann distribution
(faint dashed curve), which is just a Gaussian function. Per the discussion
above, this implies that gravitational heating leads to a longer tail in the
cut-off to the power spectrum, in qualitative agreement with the findings in
Sec. 5.

Motivated by Eq. (B.6), we fit the simulation cut-off $R_{\rm fs}(k,a)$ (Eq.
5.1) to a function of the form

$\left[2.43\mathrm{e}^{(k/k_{\rm fs})^{1.99}}-1.43\right]^{-2.81}.$ (B.7)

We find that the above fitting form works well for the RD1, RD2, RD3, YD1,
YD2, YD3, and YD4 simulations. The cut-off scale obtained by this fitting
follows the relation

$\frac{k_{\rm fs}(a)}{{k}_{\mathrm{RH}}}=\frac{7.41}{\sigma_{\rm
t,RH}}\left[\ln\left(\frac{a}{{a}_{\mathrm{RH}}}\right)\right]^{-1}.$ (B.8)

Figure 16 shows the $R_{\rm fs}$ obtained from the YD1 simulation, with the
function given by Eq. (B.7) and $k_{\rm fs}$ given by Eq. (B.8). The error
between the exponential fitting function and the simulation-based cut-off
(lower panel) remains within 10% for modes for which $R_{\rm fs}\gtrsim 0.2$.
For the RD4 and RD5 cases, the fall-off in $R_{\rm fs}$ is much shallower than
in the other seven cases because of lower structure formation before
${a}_{\mathrm{RH}}$. The prescription in this section does not adequately fit
the free-streaming cut-off in these cases. In cases with significant structure
formation during the EMDE ($\gtrsim 20\%$), the exponential fit works well.

Figure 16: $R_{\rm fs}=\sqrt{\mathcal{P}_{\rm sim}/\mathcal{P}_{\rm lin}}$
from the YD1 simulation box, shown at different times after reheating. The
blue curves show the function given by Eq. (B.7) with $k_{\rm fs}$ given by
Eq. (B.8), while the bottom panels show the relative error between the $R_{\rm
fs}$ from simulations and the blue curves.

Compared to the simpler form in Eq. (5.2), the cut-off function in Eq. (B.7)
is not an obviously better fit to the simulation results, but it has a
stronger theoretical motivation. However, several simplifying assumptions were
made in associating the cut-off with the velocity distribution. First, the
velocity distribution was assumed to be the same everywhere. In practice, the
velocities arise from structure formation and are hence inhomogeneously
distributed. Although particles in the high-velocity tail should rapidly
homogenize, the low-velocity portion of Fig. 15 (left panel) may include a
significant contribution from coherent large-scale motions, which would not
suppress power. Additionally, free streaming is assumed to begin
instantaneously, whereas in practice the $Y$ particles decay over an extended
time interval. These considerations limit the degree to which the velocity
distribution can be assumed to predict the correct power spectrum cut-off. For
this reason, we continue to adopt in the main text the simpler form obtained
in Sec. 5 from the simulation power spectra alone.

##  Appendix C Linear growth function in mixed matter-radiation domination

Here we derive an expression that is valid long after reheating for the growth
of perturbations that entered the horizon before reheating. During the
radiation-dominated epoch that follows an EMDE, we can write the general
growth function as

$D(a)\propto\ln(Ba/{a}_{\mathrm{RH}})$ (C.1)

for subhorizon modes, where $B$ is a constant and ${a}_{\mathrm{RH}}$ is
defined as in Eq. (2.3). Note that

$B=\mathrm{e}^{(\mathrm{d}\ln D/\mathrm{d}\ln a)^{-1}}/(a/{a}_{\mathrm{RH}}).$
(C.2)

We numerically evaluate $D(a)$ starting during an arbitrary EMDE at
$a=10^{-6}{a}_{\mathrm{RH}}$ and find that the value of $B$ given by Eq. (C.2)
rapidly settles for $a\gtrsim 4{a}_{\mathrm{RH}}$ to a value of $B\simeq
1.733$. That is, the growth function during the radiation-dominated epoch
following an EMDE is

$D(a)\propto\ln(1.733a/{a}_{\mathrm{RH}})$ (C.3)

for modes that are already subhorizon.

At the small scales relevant to EMDE-enhanced structure, baryons are a
homogeneous background even during the late matter epoch. Reference [46]
derived the linear growth solutions during the transition to late matter
domination under this assumption; they are

$D_{i}(a)=(1+a/a_{\rm
eq})^{\mu_{i}}\,_{2}F_{1}\\!\left(-\mu_{i},\frac{1}{2}-\mu_{i};\frac{1}{2}-2\mu_{i};\frac{1}{1+a/a_{\rm
eq}}\right)$ (C.4)

for $i=1$ and 2, where ${}_{2}F_{1}$ is the hypergeometric function. Here

$\mu_{i}=\pm\frac{5}{4}\sqrt{1-\frac{24}{25}f_{\rm b}}-\frac{1}{4}$ (C.5)

with $+$ and $-$ for $i=1$ and $2$, respectively, where $f_{\rm b}\simeq
0.157$ is the fraction of the late matter density that is in baryons. Note
that when $a\gg a_{\rm eq}$, the hypergeometric function approaches $1$, so
$D_{1}\propto a^{\mu_{1}}$ becomes the growing mode and $D_{2}\propto
a^{\mu_{2}}$ the decaying mode. In the radiation-dominated $a\ll a_{\rm eq}$
limit, Ref. [46] showed that

$D_{i}(a)\to-\frac{\Gamma(1/2-2\mu_{i})}{\Gamma(-\mu_{i})\Gamma(1/2-\mu_{i})}\left[\ln(a/a_{\rm
eq})+\psi(-\mu_{i})+\psi(1/2-\mu_{i})-2\psi(1)\right],$ (C.6)

where $\Gamma(x)$ is the gamma function and
$\psi(x)\equiv\mathrm{d}\ln\Gamma(x)/\mathrm{d}x$ is the digamma function.

Let us write the growth function as

$D(a)=A_{1}D_{1}(a)+A_{2}D_{2}(a).$ (C.7)

In the radiation-dominated $a\ll a_{\rm eq}$ limit, we can match Eqs. (C.3)
and (C.7) as well as their first derivatives to obtain

$\displaystyle A_{2}$
$\displaystyle=-A_{1}\frac{\Gamma(-\mu_{2})\Gamma(1/2-\mu_{2})\Gamma(1/2-2\mu_{1})}{\Gamma(-\mu_{1})\Gamma(1/2-\mu_{1})\Gamma(1/2-2\mu_{2})}$
$\displaystyle\hphantom{=}\times\frac{\ln(1.733a_{\rm
eq}/{a}_{\mathrm{RH}})-\psi(-\mu_{1})-\psi(1/2-\mu_{1})+2\psi(1)}{\ln(1.733a_{\rm
eq}/{a}_{\mathrm{RH}})-\psi(-\mu_{2})-\psi(1/2-\mu_{2})+2\psi(1)}.$ (C.8)

Meanwhile, if we normalize the growth function so that

$D(a)\to A_{1}(a/a_{\rm eq})^{\mu_{1}}\equiv a^{\mu_{1}}$ (C.9)

in the matter-dominated $a\gg a_{\rm eq}$ limit, then

$A_{1}=a_{\rm eq}^{\mu_{1}}.$ (C.10)

The linear growth function during mixed matter-radiation domination, for modes
that entered the horizon prior to the last radiation-dominated epoch, is thus
given by the combination Eq. (C.7) of the solutions given by Eq. (C.4) with
coefficients given by Eqs. (C.10) and (C).

##  Appendix D Code

The routines for the time-dependent critical density threshold taken from Ref.
[20], the EMDE transfer functions from Ref. [30], the growth factor used in
this paper, and the free-streaming cut-off are available here.

## Acknowledgments

H.G. is supported by NSF Grant AST-2108931. H.G. thanks Adrienne Erickcek for
guiding him through the early stages of this project, and Katherine Mack for
her constant support and helpful comments. Simulations for this work were run
on the Henry2 / Hazel High Performance Computing Cluster at North Carolina
State University.

## References

  * [1] A.H. Guth, _Inflationary universe: A possible solution to the horizon and flatness problems_ , _Phys. Rev. D_ 23 (1981) 347.
  * [2] A.D. Linde, _A new inflationary universe scenario: A possible solution of the horizon, flatness, homogeneity, isotropy and primordial monopole problems_ , _Physics Letters B_ 108 (1982) 389.
  * [3] A. Albrecht and P.J. Steinhardt, _Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking_ , _Phys. Rev. Lett._ 48 (1982) 1220.
  * [4] R. Allahverdi et al., _The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe_ , _The Open Journal of Astrophysics_ 4 (2021) 1 [2006.16182].
  * [5] G. Kane, K. Sinha and S. Watson, _Cosmological moduli and the post-inflationary universe: A critical review_ , _International Journal of Modern Physics D_ 24 (2015) 1530022 [1502.07746].
  * [6] X. Chen and S.H.H. Tye, _Heating in brane inflation and hidden dark matter_ , _J. Cosmology Astropart. Phys._ 2006 (2006) 011 [hep-th/0602136].
  * [7] Y. Zhang, _Long-lived light mediator to dark matter and primordial small scale spectrum_ , _J. Cosmology Astropart. Phys._ 2015 (2015) 008 [1502.06983].
  * [8] A. Berlin, D. Hooper and G. Krnjaic, _PeV-scale dark matter as a thermal relic of a decoupled sector_ , _Physics Letters B_ 760 (2016) 106 [1602.08490].
  * [9] A. Berlin, D. Hooper and G. Krnjaic, _Thermal dark matter from a highly decoupled sector_ , _Phys. Rev. D_ 94 (2016) 095019 [1609.02555].
  * [10] J.A. Dror, E. Kuflik and W.H. Ng, _Codecaying dark matter_ , _Phys. Rev. Lett._ 117 (2016) 211801.
  * [11] J.A. Dror, E. Kuflik, B. Melcher and S. Watson, _Concentrated dark matter: Enhanced small-scale structure from codecaying dark matter_ , _Phys. Rev. D_ 97 (2018) 063524 [1711.04773].
  * [12] A.L. Erickcek, P. Ralegankar and J. Shelton, _Cannibalism’s lingering imprint on the matter power spectrum_ , _J. Cosmology Astropart. Phys._ 2022 (2022) 017 [2106.09041].
  * [13] P.F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti, _Bounds on very low reheating scenarios after Planck_ , _Phys. Rev. D_ 92 (2015) 123534 [1511.00672].
  * [14] T. Hasegawa, N. Hiroshima, K. Kohri, R.S.L. Hansen, T. Tram and S. Hannestad, _MeV-scale reheating temperature and thermalization of oscillating neutrinos by radiative and hadronic decays of massive particles_ , _J. Cosmol. Astropart. Phys._ 2019 (2019) 012 [1908.10189].
  * [15] A.L. Erickcek and K. Sigurdson, _Reheating effects in the matter power spectrum and implications for substructure_ , _Phys. Rev. D_ 84 (2011) 083503 [1106.0536].
  * [16] G. Barenboim and J. Rasero, _Structure formation during an early period of matter domination_ , _Journal of High Energy Physics_ 2014 (2014) 138 [1311.4034].
  * [17] J. Fan, O. Özsoy and S. Watson, _Nonthermal histories and implications for structure formation_ , _Phys. Rev. D_ 90 (2014) 043536 [1405.7373].
  * [18] A.L. Erickcek, _The dark matter annihilation boost from low-temperature reheating_ , _Phys. Rev. D_ 92 (2015) 103505 [1504.03335].
  * [19] A.L. Erickcek, K. Sinha and S. Watson, _Bringing isolated dark matter out of isolation: Late-time reheating and indirect detection_ , _Phys. Rev. D_ 94 (2016) 063502 [1510.04291].
  * [20] C. Blanco, M.S. Delos, A.L. Erickcek and D. Hooper, _Annihilation signatures of hidden sector dark matter within early-forming microhalos_ , _Phys. Rev. D_ 100 (2019) 103010 [1906.00010].
  * [21] M.S. Delos, T. Linden and A.L. Erickcek, _Breaking a dark degeneracy: The gamma-ray signature of early matter domination_ , _Phys. Rev. D_ 100 (2019) 123546 [1910.08553].
  * [22] H. Ramani, T. Trickle and K.M. Zurek, _Observability of Dark Matter Substructure with Pulsar Timing Correlations_ , _JCAP_ 12 (2020) 033 [2005.03030].
  * [23] V.S.H. Lee, A. Mitridate, T. Trickle and K.M. Zurek, _Probing small-scale power spectra with pulsar timing arrays_ , _Journal of High Energy Physics_ 2021 (2021) 28 [2012.09857].
  * [24] M.S. Delos and T. Linden, _Dark matter microhalos in the solar neighborhood: Pulsar timing signatures of early matter domination_ , _Phys. Rev. D_ 105 (2022) 123514 [2109.03240].
  * [25] M. Oguri, J.M. Diego, N. Kaiser, P.L. Kelly and T. Broadhurst, _Understanding caustic crossings in giant arcs: Characteristic scales, event rates, and constraints on compact dark matter_ , _Phys. Rev. D_ 97 (2018) 023518 [1710.00148].
  * [26] J.M. Diego, N. Kaiser, T. Broadhurst, P.L. Kelly, S. Rodney, T. Morishita et al., _Dark Matter under the Microscope: Constraining Compact Dark Matter with Caustic Crossing Events_ , _ApJ_ 857 (2018) 25 [1706.10281].
  * [27] L. Dai and J. Miralda-Escudé, _Gravitational Lensing Signatures of Axion Dark Matter Minihalos in Highly Magnified Stars_ , _AJ_ 159 (2020) 49 [1908.01773].
  * [28] N. Blinov, M.J. Dolan, P. Draper and J. Shelton, _Dark matter microhalos from simplified models_ , _Phys. Rev. D_ 103 (2021) 103514 [2102.05070].
  * [29] G. Barenboim, N. Blinov and A. Stebbins, _Smallest remnants of early matter domination_ , _J. Cosmology Astropart. Phys._ 2021 (2021) 026 [2107.10293].
  * [30] H. Ganjoo, A.L. Erickcek, W. Lin and K.J. Mack, _The effects of relativistic hidden sector particles on the matter power spectrum_ , _J. Cosmology Astropart. Phys._ 2023 (2023) 004 [2209.02735].
  * [31] Planck Collaboration, _Planck 2015 results. XIII. Cosmological parameters_ , _A &A_ 594 (2016) A13 [1502.01589].
  * [32] W.H. Press and P. Schechter, _Formation of Galaxies and Clusters of Galaxies by Self-Similar Gravitational Condensation_ , _ApJ_ 187 (1974) 425.
  * [33] J.R. Bond, S. Cole, G. Efstathiou and N. Kaiser, _Excursion Set Mass Functions for Hierarchical Gaussian Fluctuations_ , _ApJ_ 379 (1991) 440.
  * [34] V. Springel, _The cosmological simulation code GADGET-2_ , _MNRAS_ 364 (2005) 1105 [astro-ph/0505010].
  * [35] M.S. Delos, A.L. Erickcek, A.P. Bailey and M.A. Alvarez, _Are ultracompact minihalos really ultracompact?_ , _Phys. Rev. D_ 97 (2018) 041303 [1712.05421].
  * [36] M.S. Delos, A.L. Erickcek, A.P. Bailey and M.A. Alvarez, _Density profiles of ultracompact minihalos: Implications for constraining the primordial power spectrum_ , _Phys. Rev. D_ 98 (2018) 063527 [1806.07389].
  * [37] Y. Suto, _Correlation Function in Decaying Particle Cosmology_ , _ApJ_ 321 (1987) 36.
  * [38] K. Enqvist, S. Nadathur, T. Sekiguchi and T. Takahashi, _Decaying dark matter and the tension in $\sigma$8_, _J. Cosmology Astropart. Phys._ 2015 (2015) 067 [1505.05511].
  * [39] J. Dakin, S. Hannestad and T. Tram, _Fully relativistic treatment of decaying cold dark matter in N-body simulations_ , _J. Cosmology Astropart. Phys._ 2019 (2019) 032 [1904.11773].
  * [40] A.J. Benson, A. Farahi, S. Cole, L.A. Moustakas, A. Jenkins, M. Lovell et al., _Dark matter halo merger histories beyond cold dark matter - I. Methods and application to warm dark matter_ , _MNRAS_ 428 (2013) 1774 [1209.3018].
  * [41] A.M. Green, S. Hofmann and D.J. Schwarz, _The power spectrum of SUSY-CDM on subgalactic scales_ , _MNRAS_ 353 (2004) L23 [astro-ph/0309621].
  * [42] V.S. Berezinsky, V.I. Dokuchaev and Y.N. Eroshenko, _Formation and internal structure of superdense dark matter clumps and ultracompact minihaloes_ , _J. Cosmology Astropart. Phys._ 2013 (2013) 059 [1308.6742].
  * [43] M.S. Delos and J. Silk, _Ultradense dark matter haloes accompany primordial black holes_ , _MNRAS_ 520 (2023) 4370 [2210.04904].
  * [44] M.S. Delos and G. Franciolini, _Lensing constraints on ultradense dark matter halos_ , _Phys. Rev. D_ 107 (2023) 083505 [2301.13171].
  * [45] M. Ackermann, M. Ajello, A. Albert, W.B. Atwood, L. Baldini, J. Ballet et al., _The Spectrum of Isotropic Diffuse Gamma-Ray Emission between 100 MeV and 820 GeV_ , _ApJ_ 799 (2015) 86 [1410.3696].
  * [46] W. Hu and N. Sugiyama, _Small-Scale Cosmological Perturbations: an Analytic Approach_ , _ApJ_ 471 (1996) 542 [astro-ph/9510117].