Simulations of Gravitational Heating Due to Early Matter Domination

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.


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 matterdominated 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 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 freestreaming 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 = ℏ = 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 hs , while the temperature of the SM is T .The X particles are nonrelativistic with their energy density evolving with scale factor as ρ X ∝ a −3 .The Y particles are initially relativistic but become nonrelativistic as T 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 Γ.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 → Y Y .
The equations for the homogeneous energy densities of the three components are: where overdots denote d/dt and H ≡ ȧ/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 ρ 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 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 ρ R ∝ a −4 when when Γmn Y ≪ Hρ R .When Γmn Y becomes comparable to Hρ R , ρ R ∝ 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, ρ R ∝ a −4 again.
At the start of our calculations, ρ Y ∝ a −4 because the Y particles are relativistic.As T 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 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, ρ Y ∝ a −3 until Γ/H ∼ 1.When ρ Y exceeds ρ R , the universe becomes matter-dominated.This phase is shown by the yellow shaded region in Figure 1.When Γ/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 which equates Γ to what the Hubble rate would be in a purely radiation-dominated universe at temperature T RH .Another quantity that is useful to define is the scale factor a RH of reheating, given by where a 0 and T 0 are respectively the scale factor and the radiation temperature today [30].
Note that both T RH and a RH are defined quantities and T (a RH ) ̸ = T RH .
The behavior of the homogeneous background universe is controlled by three parameters.The ratio η 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, η and m determine when ρ Y exceeds ρ R .For a fixed η, 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 η means the EMDE begins later.In addition to these, T RH determines when the EMDE ends by setting the decay rate of the Y particles -a smaller value of T RH for the same η and m means that reheating happens later and the EMDE is longer.
The case shown in Fig. 1 has η = 100, which implies that SM radiation dominates the energy budget of the universe before the EMDE.There can be cases with η < 1, in which ρ Y > ρ R initially and the EMDE begins when the Y particles become nonrelativistic.

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 δ 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 δ 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 δ X track δ Y and the dark matter (X) power spectrum inherits the shape and amplitude of the δ 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 (taken from G23) illustrates the time evolution of two perturbation modes in Newtonian gauge, showing δ X and δ Y scaled by Φ 0 , the primordial metric perturbation in a radiation-dominated universe.It also shows the evolution of δ Y,c /Φ 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, δ X , δ Y and δ 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 δ Y is initially suppressed compared to δ Y,c , before the pressure subsides and δ Y starts growing linearly with scale factor during the EMDE.In addition, δ X starts tracking δ Y after some time has elapsed in the EMDE, on account of falling into the gravitational wells created by the Y particles.
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 P(k) ≡ [k 3 /(2π 2 )]P (k) evaluated in linear theory: where A s = 2.2 × 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 ζ, and the remaining factors translate this into the matter power spectrum (since the primordial gravitational potential Φ 0 = −2ζ/3).The power spectrum is shown at a RH = 1.64 × 10 −11 for a case with a bosonic Y particle of m = 1 TeV, η = 100 (so that SM radiation initially dominates), and T RH = 10 MeV.The vertical dashed line shows k RH ≡ a RH Γ.For modes with k ≳ k RH , the power spectrum is enhanced due to the EMDE.Modes that enter the horizon during the EMDE (k RH < k < k dom ) grow linearly with a from horizon entry until a RH .For these modes, δ X (k, a RH ) ∝ (k/k RH ) 2 and P(k) ∝ (k/k RH ) 4 .Modes with k > k 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 This transition to logarithmic shape at k dom is apparent in Fig. 3.In contrast, for a case with initial Y domination (η < 1), P(k) ∝ (k/k RH ) 4 for all k > k 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 ≳ 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 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.

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 bound of matter in halos is at the scale factor a, where σ(a) = [ P(k, a)d ln k]1/2 is the rms of the unfiltered linear density contrast field, δ c ≃ 1.686 is the linear threshold for spherical collapse, and erfc(x) ≡ 1 − erf(x) is the complementary error function. 1To quantify the dominance of nonlinear structure formation during the EMDE, we consider f bound (a RH ), the fraction of Y particles bound in halos by the time of reheating.For a case with {m, η, T RH } = {2 TeV, 10, 20 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 ρ Y ≫ ρ 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 λ fs .Reference [20] estimated that structures formed after about 10 −6 a 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 σ v , the free-streaming length is given by Assuming that nonlinear growth stops at a RH and the velocity dispersion redshifts after that as σ v (a) = σ vRH a RH /a, the above integral can be recast as During the period of radiation domination that follows the EMDE, we have H(a) ∝ a −2 , where we neglect the minor influence of changes in the relativistic content.This implies a 2 H(a) = a 2 RH H(a RH ), using which the integral is Finally, defining k fs = λ −1 fs and using k RH = a RH Γ = 0.96a RH H(a RH ), 2 where H(a RH ) is the Hubble rate in a radiation-dominated universe, we can write Note that this discussion is conceptual, as σ vRH 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 ≳ k 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 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 da/[a 3 H(a)] converges in the infinite future.
the EMDE remains even after the impact of gravitational heating, and microhalos can start forming much earlier than they do in standard ΛCDM scenarios.Since the density of a halo forming at a f scales as a −3 f , 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.

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 (ρ Y + ρX )/ρ X ; 3 this is conceptually equivalent to scaling the simulation particle mass M as 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.02a RH , with the exception of two cases (YD1 and YD2), which began at the earlier time a = 0.01a RH and 0.005a 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, σ = [ P(k, a start )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) , and the percentage of matter bound in halos at reheating, evaluated using Eq.(3.1).
Table 1 lists the hidden-sector model parameters {m, η, T 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 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.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 P(k) from the simulation as dashed lines and the linear theory P(k) as solid lines.For a = 0.10a 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.13aRH , 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.27aRH , 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.24aRH , 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.

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 for a > 10a RH .Here we evaluate the linear-theory power spectrum P lin (k, a) by using the Boltzmann solver described in G23 to obtain P lin (k, 10a RH ) and extrapolating it to later times using P lin (k, a) ∝ [ln(a/a RH )] 2 .
We analyzed snapshots from 11a RH to 790a RH for each case.Figure 6 shows R fs at different times taken from the RD5 simulation box.As a/a RH increases, R fs (k, a) decreases and the wavenumber at which R 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.
For each snapshot, R fs (k, a) was computed and fit to a function of the form where k fs and n were treated as free parameters.We found that n ≈ 2.5 fits the fall-off in R fs for the range of our cases, so we subsequently fixed n = 2.5 and fit this function to the simulation R fs while treating only k fs as a fitting parameter.Figure 7 shows the R fs obtained  As in Fig. 6, the noise at low k is an artifact of the finite box size and is unrelated to free streaming.
This fitting procedure yields a value of k fs for each value of a. Figure 8 shows the values of k fs /k RH obtained by this fitting for the RD1, RD3 and YD1 cases as red crosses, plotted against [ln(a/a RH )] −1 .Guided by Eq. (3.5), we fit the array of k fs values obtained to the function treating α 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 ≃ 106a RH to a ≃ 790a 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 α, showing in Table 2.

Relating the Cut-Off To The Linear Power Spectrum
According to Eq. (3.5), the value of α is inversely proportional to the average velocity dispersion of the dark matter particles at reheating, denoted by σ vRH .To test this relation, we use the linear power spectrum evaluated at a RH to calculate the linear velocity dispersion for each of our nine cases.The power spectrum at a RH determines the linear velocity dispersion, which we take to be (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 σ lin,RH α.
Two distinct classes can be identified from Table 2.For the six cases RD1, RD3, YD3, YD4, RD4 and RD5, the product σ lin,RH α has an average of 6.63 × 10 5 km/s, with the maximum deviation of the product from the average being 7.7% (YD3 for these cases, α ∝ 1/σ 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 σ lin,RH α values are those with significant structure having formed at a 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 P(k * , a RH ) crosses some threshold value C. We define the truncated linear velocity dispersion to be where is the truncated linear power spectrum.The only motivation of this form is that it suppresses the contribution of modes with P ≫ C in a continuous way and fits the simulation results well, as we will see.
Various values of C were tried to obtain values of σ t,RH that would bring the product σ t,RH α for the three high-bound-fraction cases closer to the products in the six other cases.We found C = 4.3 to be an appropriate threshold.Table 3 shows the truncated linear velocity dispersion and the products σ t,RH α for all nine cases with 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 σ t,RH α are close to each other for all the nine cases, with a maximum deviation of 7.7% from the mean value.
Case Bound % at a RH σ t,RH (km/s) α σ t,RH α/10 The truncated linear velocity dispersion given by Eq. (5.5) with a threshold of C = 4.3 (c.f.Table 2).The second column from the right lists the product of σ t,RH and α.
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 σ t,RH is inversely proportional to α.
This prescription of calculating σ t,RH yields a velocity dispersion that is inversely proportional to α within 7.7% error for our nine cases, covering the wide range of scenarios, with the bound percentages at a 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 where the numerical factor is the average of the σ t,RH α values from Table 3 (converted to units of c).
For the RD4 and RD5 cases, since the amount of structure formed at a RH is not very high, the free-streaming fall-off in R fs is a little shallower than a (k/k fs ) −2.5 power law. Figure 9 shows R fs from the RD4 simulations compared to the fitting form given by (1 + (k fs /k RH ) 2.5 ) −1 .At 66.22a RH , the slopes of the fall-off from the simulation and the fitting function are noticeably different.At 889.59a RH , the slopes are much closer to each other.The relative error between the R fs from the simulation and the best fit curve remains within 10% for R fs ≳ 0.25.Cases with such low amounts of bound structure at a RH present the boundaries down to which the (1 + (k fs /k RH ) 2.5 ) −1 fitting form is valid.For cases with lower bound percentages at a RH , the fall-off in R fs is much shallower than (k/k fs ) −2.5 ; we leave the detailed study of these cases to future work.
Finally, we remark that due to discreteness noise, it is difficult to resolve the precise form of the free-streaming cut-off for k ≫ k 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 fs and the free-streaming scale k 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).

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 λ fs ∝ ln(a/a RH ) is derived under the assumption of radiation domination.During the transition to late matter domination, (compare Eq. 3.4), where we assume a eq ≫ a RH and again neglect changes to the radiation content.This expression is approximately ln(a/a RH ) when a ≪ a eq but asymptotes to ln(4a eq /a RH ) in the a ≫ a eq limit.We therefore evaluate the free-streaming cutoff arising from gravitational heating by replacing ln(a/a RH ) with the logarithm in Eq. (6.1) in the expression for the cut-off scale, Eq. (5.7).

The linear matter power spectrum
Figure 10 shows the ΛCDM power spectrum P at a eq (black solid curve).We first evaluated the linear power spectrum using CAMB at z = 500 and translated it back in time to a 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 Λ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 η 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 RH as the green dashed curve, but its value of η 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 2 fs from Eq. (5.2) with cut-off scale k fs given by Eq. (5.7) (modified as described after Eq. 6.1) at a = a 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 freestreaming effects is roughly four orders of magnitude larger than the ΛCDM counterpart at the same scale.The Λ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 eq .
It is noteworthy that with gravitational heating accounted for, all three EMDE power spectra in Fig. 10 peak around P ∼ 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 P ∼ 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 σ t,RH ∼ k RH /k * , where k * is the nonlinear scale (where P ∼ 1) at reheating.Equation (5.7) then implies that k fs ∼ k * / ln(a/a RH ).Shortly after reheating, the power spectrum thus peaks at P ∼ 1 (see also Fig. 5), since k fs is close to the nonlinear scale k * .Over time, k 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 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 P ∼ [ln(a/a RH )] 2 (k fs /k * ) n ∼ [ln(a/a 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 P ∼ 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 P(k pk , a eq )/P(k RH , a eq ), the factor by which density variations at the scale of the peak in the power spectrum exceed density variations at k 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 P(k pk , a eq )/P(k RH , a eq ) rises from ∼ 10 3.5 , its gravitational heating-suppressed value drops from ∼ 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, η and T RH is varied while the others are held fixed.The red and green curves are obtained by varying η; as indicated by the comparison between the blue and green dashed curves in Fig. 10, smaller values of η raise the peak value of 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 P(k) (see Eq. 5.4), which explains why the gravitational heating-suppressed value of P(k pk , a eq )/P(k RH , a 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 η is large, increasing m raises the peak value of 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 P raises it while moving the peak to smaller scales lowers itwhich is why the purple curve in Fig. 11 drops more shallowly.
Figure 11 suggests that the free-streaming-suppressed value of P(k pk , a eq )/P(k RH , a eq ) does not drop below a factor of ∼ 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 P(k) ≃ 0.15A s (k/k 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 (P ≲ 1), it converges to a finite value that evaluates to σ t,RH ≃ 0.73A 1/4 s ≃ 0.005 (1500 km/s), assuming A s = 2.2×10 −9 .The free-streaming scale at matter-radiation equality is then k fs ∼ σ −1 t,RH [ln(a eq /a RH )] −1 k RH , and thus P(k pk , a eq )/P(k RH , a eq ) ∼ (k fs /k RH ) 2 ∼ σ −2 t,RH [ln(a eq /a RH )] −2 .Since σ t,RH ≃ 0.005 The logarithm of P(k pk , a eq )/P(k RH , a 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 η for different values of m, while the purple curve corresponds to varying m with fixed η.
and ln(a eq /a RH ) ∼ O (10), small-scale density variations remain strongly enhanced even in this limiting scenario.

Dark matter halo formation
To illustrate how the gravitational heating cut-off affects halo formation, we consider the halo mass function dn/d ln M , the number density of bound in logarithmic bins of mass, at different times around the epoch of matter-radiation equality.According to Press-Schechter theory [32], where σ(M ) is the rms density contrast smoothed on the mass scale M and ρX is the comoving density of dark matter.Following common practice for cosmologies with a small-scale cutoff to the power spectrum, we use a sharp k-space smoothing filter and evaluate σ 2 (M ) = k M 0 P(k)dk/k with the smoothing wavenumber k M defined by M ≡ (4π/3)ρ 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, δ c , which is valid for mixed matter-radiation domination.
Figure 12 shows the mass function for the case m = 2 TeV, η = 50 and T 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 eq .For a < a eq , the mass function shown by the solid lines grows as more small-scale halos form.For a > a eq , the same mass function decreases with time for M ≲ 10 −8 M ⊙ 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 Mpc −3 at a eq if a Λ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 fs ≈ 40k 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.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 eq for the three cases in Fig. 10

+ z
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 σ 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 Λ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.which we now label using σ eq , the value that the rms linear density fluctuation would reach at a = a eq if gravitational heating were neglected.That is, where P emde is the power spectrum without the free-streaming cut-off (corresponding to the dashed curves in Fig. 10).σ 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 σ 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 σ 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 σ 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 Λ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 ≈ 10a eq in cases with an EMDE.
Figure 14 shows, as a function of σ 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 σ = 0.608δ c (σ = 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, η and T 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 σ eq ≳ 40).In this regime, increasing σ 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 η causes a 10 to change more mildly as a function of σ eq , compared to varying m, due to how the scale of the peak in 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 eq (redshift 1000) in all cases.
The dotted curve in Fig. 14 shows how a 10 varies with σ eq when gravitational heating is neglected.Under this assumption, σ 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 σ eq leads to earlier halo formation, and when σ eq ≳ 1, a significant halo population forms during the radiation era,5 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 ≲ σ eq ≲ 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 σ eq results in earlier halo formation until σ eq ∼ O (10), beyond which point the trend should reverse.The scale factor at which 10% of the dark matter is predicted to be bound in halos according to Press-Schechter theory, plotted against σ 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 η (red and green) and m (purple) to vary σ 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 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 σ eq regardless of the power spectrum shape.

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,24,28].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 P(k) during the subsequent radiation epoch by a factor that is fit well by the function [1 + (k/k fs ) 2.5 ] −2 (or an alternative form given in Appendix B).The free-streaming wavenumber k 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 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 eq .In cases where ∼ 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 eq .Without an EMDE, the formation of a similar amount of bound structure would not occur until a ∼ 100a eq (z ∼ 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 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 ∼ 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 ∼ 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-boundfraction 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 radiationdominated epoch, and further work is separately needed to understand the halos that arise this rescaled version of the velocity distribution.This means that -up to a rescaling in kthe cut-off function R 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, where f (x, 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 ∂f /∂t + iv for some initial time t = 0. Now approximate that the distribution function is separable, , the velocity distribution f v is independent of position.Integrating Eq. (B.3) over velocities then yields where fv is the Fourier transform of f v .The k-space density ρ(k, t) is evidently suppressed, in relation to ρ(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 v/σ, where σ = ⟨v 2 ⟩ 1/2 is the velocity dispersion, converges quickly after reheating to a form that is well approximated by f (x) = 3.98 × 10 −11 e (x/14.9) 1.94 − 0.953 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 f (x) = 2.43e (x/3.69)(B.8) Figure 16 shows the R fs obtained from the YD1 simulation, with the function given by Eq. (B.7) and k 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 fs ≳ 0.2.
For the RD4 and RD5 cases, the fall-off in R fs is much shallower than in the other seven cases because of lower structure formation before a 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 (≳ 20%), the exponential fit works well.
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 lowvelocity 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.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 with + and − for i = 1 and 2, respectively, where f b ≃ 0.157 is the fraction of the late matter density that is in baryons.Note that when a ≫ a eq , the hypergeometric function approaches 1, so D 1 ∝ a µ 1 becomes the growing mode and D 2 ∝ a µ 2 the decaying mode.In the radiation-dominated a ≪ a eq limit, Ref. [46] showed that Meanwhile, if we normalize the growth function so that D(a) → A 1 (a/a eq ) µ 1 ≡ a µ 1 (C.9) in the matter-dominated a ≫ a eq limit, then A 1 = a µ 1 eq . (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.8).

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.

Figure 1 :
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 (η > 1); scenarios are also possible in which relativistic Y particles initially dominate (η < 1).

Figure 2 :
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 hs /m = 0.00018 at horizon entry, corresponding to equation of state parameter w Y = 0.00017); there is no suppression of δ Y for this mode.The mode in the right panel enters the horizon when the Y particles have significant pressure (T hs /m = 0.21 at horizon entry, corresponding to w Y = 0.145) because of which δ Y is suppressed compared to δ Y,c , the Y density perturbation if the Y particles are nonrelativistic.Dark matter perturbations δ X eventually track δ Y in both cases.

Figure 3 :
Figure 3: The dimensionless power spectrum P(k) ≡ [k 3 /(2π 2 )]P (k) of dark matter density perturbations, evaluated in linear theory at a RH = 1.64×10 −11 for a cosmology with T RH = 10 MeV, η = 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 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.

Figure 4
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.10a RH , when the perturbations in the density field are small.The top right panel shows the box just after reheating, at 1.13a RH , displaying the formation of a web of structure with clumps and filaments.The bottom left panel displays the box at 11.27a RH .By this time, the structure formed at reheating is already being erased, with the clearly defined web visible at 1.13a RH blurred and smoothed out.By 389.24aRH (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.

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

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

Figure 6 :
Figure 6: Free-streaming cut-off R fs = P sim /P lin from the RD5 simulation box, shown at different times after reheating.The grey dashed horizontal line indicates R fs = 0.5.

Figure 7 :Figure 8 :
Figure7: R fs = P sim /P 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 fs from simulations and the blue curves.

Figure 9 :
Figure 9: R fs = P sim /P 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 fs from simulations and the blue curves.

Figure 10 :
Figure10: The Λ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 eq .
Figure 11: The logarithm of P(k pk , a eq )/P(k RH , a 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 η for different values of m, while the purple curve corresponds to varying m with fixed η.

3 a 1 a/a eq = 10 a/a eq = 20 aFigure 12 :
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, η = 50, T 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.

m = 2 Figure 14 :
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 σ 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 η (red and green) and m (purple) to vary σ 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 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 σ eq regardless of the power spectrum shape.
), where x = v/σ.This distribution has a significantly more pronounced lowvelocity tail than the Maxwell-Boltzmann distribution (faint dashed curve).

Figure 15 :
Figure15: 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 σ = ⟨v 2 ⟩ 1/2 .The distribution evidently does not vary after about 3a 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.

Figure 16 :
Figure16: R fs = P sim /P lin from the simulation box, shown at different times after reheating.The blue curves show the function given by Eq. (B.7) with k fs given by Eq. (B.8), while the bottom panels show the relative error between the R fs from simulations and the blue curves.

4 ) for i = 1 and 2 , where 2 F 1
is the hypergeometric function.Here µ

Table 1 :
Bound % at a RH List of simulations with the model parameters m, η and T RH , the values of a RH and k RH

Table 2 :
).Within this error, Case σ lin,RH (km/s) Bound % at a RH The proportionality constant α relating k fs /k RH and [ln(a/a RH )] −1 for the various simulation runs.Also listed are the linear velocity dispersion σ lin,RH and the product of linear velocity dispersion and α.Note that σ lin,RH is calculated using Eq.(5.4) and converted to km/s from natural units.

Table 3 :
, 4399 − 1.43The strong low-velocity support in Eq. (B.5) evidently leads to a long tail in the Fourier transform; f is significantly larger for x ≫ 1 than the Fourier transformed Maxwell-Boltzmann distribution (faint dashed curve), which is just a Gaussian function. Pe 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.