The Effects Of Relativistic Hidden Sector Particles on the Matter Power Spectrum

If dark matter resides in a hidden sector minimally coupled to the Standard Model, another particle within the hidden sector might dominate the energy density of the early universe temporarily, causing an early matter-dominated era (EMDE). During an EMDE, matter perturbations grow more rapidly than they would in a period of radiation domination, which leads to the formation of microhalos much earlier than they would form in standard cosmological scenarios. These microhalos boost the dark matter annihilation signal, but this boost is highly sensitive to the small-scale cut-off in the matter power spectrum. If the dark matter is sufficiently cold, this cut-off is set by the relativistic pressure of the particle that dominates the hidden sector. We determine the evolution of dark matter density perturbations in this scenario, obtaining the power spectrum at the end of the EMDE. We analyze the suppression of perturbations due to the relativistic pressure of the dominant hidden sector particle and express the cut-off scale and peak scale for which the matter power spectrum is maximized in terms of the properties of this particle. We also supply transfer functions to relate the matter power spectrum with a small-scale cut-off resulting from the pressure of the dominant hidden sector particle to the matter power spectrum that results from a cold hidden sector. These transfer functions facilitate the quick computation of accurate matter power spectra in EMDE scenarios with initially hot hidden sectors and allow us to identify which models significantly enhance the microhalo abundance.

Although these microhalos do not affect the large-scale structure of the universe, they boost dark matter (DM) annihilation rates, potentially producing detectable gamma-ray signals [32][33][34][35]. The DM annihilation signal is highly sensitive to the small-scale cut-off in the matter power spectrum because the cut-off scale sets the formation times and central densities of the microhalos that form due to an EMDE [32,35,36]. For instance, changing the cut-off scale by a factor of two causes the DM annihilation boost to increase by two orders of magnitude [35]. Therefore, an accurate calculation of this small-scale cut-off is key to observationally constraining scenarios with an EMDE. In this work, we determine the smallscale cut-off scale that results from the relativistic pressure of the particle that dominates the hidden sector.
If the particle that dominates the energy density of the universe during the EMDE is initially relativistic, the growth of density perturbations is inhibited for modes that enter the horizon while the particle has significant pressure. We obtain exact solutions of the evolution of perturbations during an EMDE caused by a massive particle (which we call Y ) in the hidden sector. We include the process by which this particle transitions from relativistic to nonrelativistic behavior before dominating the energy content of the universe. While perturbation equations for a relativistic hidden sector particle have been solved previously for a single set of parameters [21], we provide analytical expressions for the power spectrum peak and cut-off scales in terms of the statistics of the Y particles and the initial ratio of densities of Y and Standard Model (SM) particles. It is also possible that the Y particle experiences cannibalistic number-changing interactions that alter the evolution of its pressure; the resulting cut-off to the matter power spectrum was computed in Refs. [26,28]. Our analysis of how the Y particle generates a cut-off in the matter power spectrum in the absence of such interactions completes our understanding of how the pressure of the particle that dominates the energy density during the EMDE inhibits the growth of dark matter perturbations during the EMDE.
We provide fitting forms for transfer functions between the cases with a hot and cold hidden sector. These transfer functions facilitate the easy computation of the power spectrum cut-off caused by the pressure suppression of density perturbations. We also consider how our transfer functions change the boost factor calculations presented in Ref. [34] (hereafter B19), in which the cut-off in the power spectrum was taken to be a Gaussian function of wavenumber with the cut-off scale set as the horizon scale when the mass of the dominant hidden sector particle is equal to the hidden sector temperature. Finally, we use our transfer functions to determine which EMDE scenarios generate observable enhancements to the microhalo population.
This paper is organized as follows. In section 2, we study the evolution of the different components of the universe in our model, including the density, sound speed and equation of state of the Y particles as they transition from being relativistic to nonrelativistic. In section 3, the evolution of the density perturbations in the Y particles and dark matter before, during, and after the EMDE is determined, and the suppression of perturbation growth due to the pressure of the Y particles is analyzed. In section 4, we present expressions for the wavenumber of the peak scale, for which the matter power spectrum is maximized. In section 5, we provide fitting forms for transfer functions for the computation of the matter power spectrum in scenarios with an initially relativistic particle dominating the hidden sector. Section 6 presents calculations of the dark matter annihilation boost and the power spectrum peak height using our transfer functions; we also discuss prospects for detecting the microhalos generated in EMDE cosmologies. Our results are summarized in section 7. The full calculation of the density, pressure, and sound speed of the Y particles is presented in Appendix A. Appendix B contains the derivations of several relations between the parameters that describe the EMDE and the properties of the Y particle. The equations that govern the evolution of perturbations and their initial conditions are detailed in Appendix C. Finally, we provide an online application for computing EMDE-enhanced power spectra with the accurate smallscale cut-off that is described in Appendix D. This paper uses natural units throughout, in which c = = k B = 1.

Evolution Of The Homogeneous Background
Our model considers a universe with three components: dark matter X; the thermal bath of relativistic SM particles, which we call radiation (denoted by the subscript R); and a particle Y with mass m that decays into SM particles. X and Y live in a hidden sector that is thermally decoupled from the Standard Model and has its own temperature T hs . The Y particles are initially relativistic but transition to nonrelativistic behavior as the temperature of the hidden sector decreases. We assume that the X particles have frozen out before our calculations begin and are nonrelativistic with m X T hs and ρ X (a) ∝ a −3 .
We first establish the evolution of the homogeneous energy densities of the various components of our model. We begin our calculations at scale factor a i , which is chosen such that T hs,i ≡ T hs (a i ) = 300m, so that the Y particles are initially relativistic. The initial SM density is set by the parameter η ≡ ρ R (a i )/ρ Y (a i ). The Y particles are weakly coupled to the SM particles with a decay rate Γ. Such couplings of the hidden sector to the Standard Model can arise via various renormalizable interactions, including the lepton portal [14,37], the Higgs portal [14,38,39], and the vector portal [14,40]. To obtain the evolution of the energy densities of these three components, the coupled equations for ρ X , ρ Y and ρ R are solved numerically:ρ (2.1a) ρ R + 4Hρ R = Γmn Y ; (2.1b) ρ X + 3Hρ X = 0, (2.1c) where overdots denote d/dt and H ≡ȧ/a. In Eq. (2.1a), n Y is the number density of Y particles, and w Y is their time-varying equation of state parameter, defined as the ratio between pressure and density, w Y ≡ P Y /ρ Y . The time evolution of w Y encodes the transition from relativistic to nonrelativistic behavior for the Y particles, which we solve for exactly; the process is detailed in Appendix A. 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 [28]. We assume that X and Y particles are coupled only gravitationally, with no momentum exchange between the two species. However, the effects of additional couplings are discussed in Sec. 4.1. Figure 1 shows the solutions to Eqs. (2.1) for a chosen set of parameters. The transition from ρ Y (a) ∝ a −4 to ρ Y ∝ a −3 can be modeled by a broken power law with a pivot scale factor given by a p /a i = bT hs,i /m where b depends only on the statistics of the Y particles. We find that b is 2.70 for bosons and 3.15 for fermions; these values of b are derived in Appendix A. It follows from Eq. (2.1b) that ρ R ∝ g * (T )T 4 ∝ a −4 when Γmn Y Hρ R , where g * (T ) is the relativistic degrees of freedom contributing to the energy density of relativistic SM particles. However, all our analytical results assume that entropy is conserved in the visible sector when Γmn Y Hρ R , so that g * S (T )a 3 T 3 is constant, where g * S is the relativistic Figure 1. The background evolution of the energy densities of the Y particles, SM radiation and dark matter (X) as a function of scale factor, for parameters m = 2 TeV and η = 1000. The pivot scale factor a p marks the transition from ρ Y ∝ a −4 to ρ Y ∝ a −3 . The yellow shaded region shows the EMDE, which begins at the scale factor a dom . At the end of the EMDE, ρ Y rapidly decreases, and the universe becomes radiation dominated.
degrees of freedom contributing to the entropy density of the SM bath. When Γmn Y exceeds 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.
The EMDE, indicated by the yellow shaded region in Figure 1, starts when ρ Y exceeds ρ R at the scale factor a dom . We parameterize this point by the temperature of the SM radiation T dom , so that ρ R (a dom ) = (π 2 /30)g * (T dom )T 4 dom . We show in Appendix B that T dom can be expressed in terms of our model parameters as where g equals the number of degrees of freedom of the Y particles, T i is the temperature of the SM radiation at a i , and f is 1 if the Y particles are bosons and 7/8 if they are fermions.
The EMDE lasts until Γ/H becomes comparable to unity. After this point, the comoving number density of the Y particles starts decreasing rapidly. Shortly thereafter, ρ Y becomes negligible and the universe transitions to radiation domination. This transition, called reheating, is not an instantaneous process, but we find it useful to define a reheating temperature T RH in terms of the decay rate as which sets Γ equal to the Hubble rate in a purely radiation-dominated universe at temperature T RH . It is also useful to define a RH as the scale factor at which Note that T RH is the quantity defined in Eq. (2.3) and does not equal T (a RH ).
In our broken-power-law model, (2.5) To relate a RH to the scale factor today (a 0 ), we note that there is negligible transfer of entropy from the decay of the Y particles to the SM radiation for a > 5a RH . We find numerically that T (5a RH ) = 0.204T RH and use entropy conservation from 5a RH to a 0 to express where T 0 is the temperature of radiation in the Universe today.

Evolution Of Perturbations
The Einstein equations are perturbed to obtain the equations for the evolution of the density contrast δ ≡ (ρ −ρ)/ρ (whereρ denotes homogeneous background density) and the velocity dispersion θ ≡ a∂ i v i for each fluid, where v i = dx i /dt. We work in the Newtonian gauge, in which the metric is given by We neglect anisotropic stress and set ψ = −φ. The perturbation equations and initial conditions are provided in Appendix C. Figure 2 shows the time evolution of |δ i |/Φ 0 , where Φ 0 is the primordial metric perturbation in a radiation-dominated universe and i denotes the three fluids in our model. Also shown is the evolution of the Y density perturbation if the Y particles were pressureless (δ Y,c ). The left panel shows a mode that enters the horizon after the Y particles have become nonrelativistic, with T hs /m = 0.00018 at horizon entry. In the absence of pressure, subhorizon density perturbations in Y grow logarithmically with scale factor during radiation domination and linearly during the EMDE. After the EMDE, radiation domination resumes and δ X and δ Y start growing logarithmically. For this mode, δ Y coincides with δ Y,c because the Y particles are already pressureless when the mode enters the horizon. In contrast, the right panel of Figure 2 shows a mode that enters the horizon when the Y particles have significant pressure, with w Y = 0.14 and T hs /m = 0.21 at horizon entry. For this mode, the growth of δ Y is suppressed compared to that of δ Y,c until the Y particles become pressureless. As a result, δ Y starts linear growth later than δ Y,c and δ Y < δ Y,c at the end of the EMDE.
The right panel of Figure 2 also shows how the evolution of δ X is affected by the pressure of the Y particles. When the mode enters the horizon during radiation domination, δ X starts to grow logarithmically with the scale factor. The pressure of the Y particles delays the onset of linear growth during the EMDE because the Y particles are not as clustered as they would have been if δ Y had also grown logarithmically prior to the EMDE. Instead of growing linearly with scale factor throughout the EMDE, δ X converges to δ Y because the X particles fall into Figure 2. The evolution of perturbations for two modes, 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); 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) because of which δ Y is suppressed compared to δ Y,c , the Y density perturbation if the Y particles are nonrelativistic. the gravitational wells generated by the Y particles. Due to this convergence, we will focus hereafter on analyzing the behavior of δ Y .
To quantify which scales undergo growth suppression, we consider the continuity and Euler equations for the evolution of density and velocity perturbations in the Y particles along with the Poisson equation. Since the comoving number density of Y particles remains constant until Γ H at the end of the EMDE, we can neglect the decay terms when the pressure of the Y particles is significant. We then have the following equations (taken from Appendix C): where c 2 sY = δP Y /δρ Y is the sound speed of the Y particles (see Appendix A). In Eq. (3.2c), the contribution of the dark matter term (δ X ρ X ) on the RHS is neglected because ρ X ρ Y . Working in the subhorizon limit where k aH and using To obtain the evolution of δ Y , we neglect the derivative of φ in Eq. (3.2a) as it is small compared to θ Y /(a 2 H) and neglect the ( (3.4) As Figure 2 shows, δ R begins oscillating shortly after the mode enters the horizon. The gravitational contribution of the δ R ρ R term on the RHS of Eq. (3.4) thus averages to zero and the term can be ignored. We can then express Eq. (3.4) as where we define the time-varying Jeans wavenumber When the Y particles are relativistic, k J is roughly proportional to a −1 because c 2 sY is constant and ρ Y ∝ a −4 . As the Y particles become colder, k J increases proportional to a 1/2 because c 2 sY ∝ a −2 and ρ Y ∝ a −3 . This behavior is shown in the top panel of Figure 3, where the black line shows the Jeans length λ J ≡ k −1 J . The sign of the coefficient of δ Y in Eq. (3.5) determines whether δ Y grows or oscillates. Figure 3 illustrates the contrast between the growing and oscillating solutions. The top panel shows the comoving length scales (k −1 ) corresponding to two different modes, plotted relative to the Jeans length. The bottom panel shows the time evolution of δ Y for the two modes. The thin lines show the evolution of each mode if the Y particles are treated as nonrelativistic (δ Y,c ). When k < k J (so that k −1 > λ J ), the coefficient of δ Y in Eq. (3.5) is negative, which leads to a growing solution for δ Y . The mode represented by the red dashed line in Figure 3 is such an example; its wavelength is always larger than the Jeans length. The bottom panel shows how the amplitude for this mode grows logarithmically with a during radiation domination and then grows linearly with a during the EMDE. In contrast, for the mode indicated by the purple dot-dashed line, the perturbation amplitude oscillates when k −1 < λ J and starts growing when k −1 > λ J . Since δ Y starts growing only when the Jeans length becomes smaller than the mode wavelength, δ Y is reduced compared to δ Y,c .
The suppression of perturbation modes that enter the Jeans horizon is readily apparent in Figure 4, which shows δ Y and δ Y,c evaluated at a RH as a function of wavenumber scaled by k RH ≡ a RH Γ. If the Y particles are always pressureless, modes that enter the horizon before the EMDE grow logarithmically with a during radiation domination and then linearly with a during the EMDE, so that δ Y,c (k > k dom , a RH ) ∝ ln(k/k dom ), where k dom ≡ a dom H(a dom ) is the horizon wavenumber at a dom . Modes that enter the horizon during the EMDE grow linearly with a from horizon entry until a RH , so that δ Y,c (k RH < k < k dom , a RH ) ∝ (k/k RH ) 2 . The shape of δ Y,c (k) after the EMDE only depends on the ratio k dom /k RH , which (as shown in Appendix B) can be expressed in terms of our model parameters as If the Y particles are initially relativistic, the growth of perturbations is suppressed for scales close to or smaller than the maximum value of the Jeans length (shown by the wavenumber k J,min = λ −1 J,max in Figure 4). For these modes, δ Y does not begin to grow until the Jeans length becomes smaller than the mode's wavelength. As a result, δ Y at a RH is increasingly suppressed compared to δ Y,c as k increases, as the blue curve in Figure 4 shows. The suppression leads to a peak in δ Y (k) at the wavenumber k pk . For k > k pk , modes start growing not only later, but also at different points in the oscillation cycles of their amplitudes. This leads to an oscillation pattern in δ Y (k) with a decaying envelope.
B19 modeled the suppression of modes that enter the horizon when the Y particle is relativistic by multiplying δ Y,c (k) by exp[−k 2 /(2k 2 y )], where k y is the wavenumber of the mode that enters the horizon when m = T hs . In Appendix B, we derive expressions for k y /k dom for a universe with η > 1. Using the expression for k y /k dom from Eq. (B.14) with Eq. (3.7) yields where g * y = g * (T (a y )), with a y /a i = T hs,i /m. The cut-off used by B19 does not describe δ Y (k) accurately: Figure 4 shows that δ Y falls off at smaller wavenumbers than k y . In section 4, we derive the model dependence of the actual peak and cut-off scales of δ Y (k).

The Peak Scale
In order to determine the observational signatures of an EMDE, it is necessary to evaluate the location and amplitude of the peak in the matter power spectrum, since this peak sets the masses, formation times, and central densities of the first microhalos [41,42]. In this section, we provide expressions for the peak wavenumber k pk for which δ Y (k)/Φ 0 is maximized.
Due to the gravitational coupling between X and Y particles during the EMDE, the peak wavenumber of δ X (k) is generally very close to that of δ Y (k). However, the peaks are not exactly equal in all cases. The relative closeness of the peaks of δ Y (k, a RH ) and δ X (k, a RH ) depends on the duration of the EMDE, quantified by k dom /k RH . Figure 5 shows δ X (k, a RH ) and δ Y (k, a RH ) for three different EMDE durations. The leftmost panel shows a short EMDE with k dom /k RH = 4.8, in which case the peak wavenumbers of δ Y and δ X differ by 10% at the end of the EMDE. This difference arises because the EMDE is too short for δ X and δ Y to become equal for modes close to the peak wavenumbers. For scales smaller than the second peak in the left panel of Figure 5, δ Y oscillates throughout the EMDE because the Jeans length does not fall below the comoving wavelengths of these modes before a RH . As a result, the Y particles do not cluster and never exert a coherent gravitational pull on the X particles. The X particles drift during the EMDE and δ X does not approach δ Y . For longer EMDEs, the peaks of δ X and δ Y are nearly identical. The middle panel of Figure 5 shows the case with k dom /k RH = 17.7, for which the peak wavenumbers of δ X and δ Y differ by 1.3%. For k dom /k RH = 32.6, this discrepancy between the peak scales falls to 0.4%. Therefore, the peak in δ Y generally matches the corresponding peak in the matter power spectrum after the EMDE.

The Effect Of Kinetic Coupling In The Hidden Sector
Thus far, we have assumed that the Y particles and the dark matter X are coupled only gravitationally. In this section, we explore how scatterings between X and Y particles affect the peak amplitude and scale of δ X and δ Y .
If the X and Y particles are initially kept in kinetic equilibrium through a scattering process, the momentum transfer rate (dp/dt)/p to the X particles from this scattering is given by n Y (m/m X ) σv , where σv is the velocity-averaged scattering cross section. This interaction modifies the Euler equation for the velocity perturbations in the dark matter [43]: where the prime denotes d/d ln a. The corresponding coupling term in the Euler equation for θ Y is suppressed by a factor of ρ X /ρ Y and can be neglected. The coupling strength is parameterized by the scale factor of kinetic decoupling a kd , which is defined by the relation n Y (a kd ) σv = H(a kd ).
To study the effect of this kinetic coupling, we consider three examples with η = 300, k dom /k RH = 36, and 1. no kinetic coupling between X and Y particles, 2. kinetic coupling with a kd = 0.5a dom = 150a p , such that the X and Y particles decouple before the EMDE starts but after the Y particles have become cold, and 3. kinetic coupling with a kd = 1.4a RH , such that the X and Y particles remain coupled until after the Y particles have decayed into SM radiation.
While n Y σv H, θ X θ Y . As the hidden sector temperature decreases, the Y particles become nonrelativistic. In this regime, we can use the results of Ref. [43] for the momentum transfer rate for the collision of two nonrelativistic particles and take σv ∝ √ T hs . For all these cases, we choose m X /m = 50 so we can assume the DM particle is much heavier than the Y particle and w X ≈ 0. This ensures that the evolution of ρ X is given by Eq. (2.1c). Figure 6. Density perturbations as a function of wavenumber, evaluated at a RH , with η = 300, k dom /k RH = 36 and m X /m = 50. δ X (k, a RH ) is plotted for cases with three different kinetic coupling strengths. The labels mark the scale factor of kinetic decoupling a kd for the different coupling scenarios. For the case with no coupling (leftmost panel), the peak wavenumbers of δ X and δ Y differ by 0.4%. Figure 6 shows δ Y (k, a RH ) and δ X (k, a RH ) for the three cases mentioned above. The amplitude and the location of the first peak remain the same between the cases. Gravitational coupling during the EMDE is strong enough to make δ X and δ Y converge to within a 0.5% difference for wavenumbers close to k pk , even without kinetic coupling. The left panel of Figure 6 also shows that δ Y and δ X differ for scales smaller than the peak scale in the case with no kinetic coupling. For these modes, δ Y oscillates for a portion the EMDE, and the Y particles do not exert a coherent gravitational force on the DM until δ Y stops oscillating. As a result, δ X has not fully converged to δ Y at a RH . Comparing the left and middle panels of Figure 6, it is apparent that scatterings tie δ X to δ Y for these modes. Therefore non-gravitational interactions between the X and Y particles only serve to tighten the correspondence between δ Y and δ X and do not significantly affect the matter power spectrum.

Scenarios With Initially Subdominant Y particles
To derive an analytical expression for k pk when the Y particle is initially subdominant (η > 1), we adopt the approach used to find the peak scale for cannibalistic hidden sector particles [28]. The peak scale enters the horizon after the Y particles have become nonrelativistic but before the onset of the EMDE. Since δ R oscillates rapidly after the peak scale enters the horizon, the term proportional to δ R ρ R in Eq. (3.4) does not affect the evolution of δ Y . In addition, the term proportional to ρ Y /(ρ Y + ρ R ) in the coefficient of δ Y on the LHS is negligible since ρ Y ρ R prior to the EMDE. Using entropy conservation in the visible sector, we can write a 2 H = a 2 i H(a i )[g * (T (a i ))/g * (T (a))] 1/6 . The first term in the coefficient of dδ Y /da in Eq. (3.4) is then proportional to d ln g * (T (a))/d ln a, which is negligible. We also set 1 + w Y to unity, since w Y is small compared to 1 and decreases as a −2 when the modes close to the peak scale enter the horizon. With these approximations, Eq. (3.4) can be written as In Appendix A, we present a piecewise model for c 2 sY (a): c 2 sY (a) = 0.33a 2 pc /a 2 for a 1.4a p . Here, a pc = 1.43a p for bosons and 1.41a p for fermions, where a p = bT hs,i a i /m is the pivot scale factor for broken power law that models ρ Y (a) and b is 2.70 for bosonic Y particles and 3.15 for fermionic Y particles. The different factors 1.41 and 1.43 reflect that the Y particles have slightly different pressure for the same value of a/a p if their statistics are different, which leads to the pivot points for their sound speed being at slightly different values of a/a p . Since the peak scale enters the horizon when a 1.4a p , we can use c 2 sY (a) = 0.33a 2 pc /a 2 in Eq. (4.2), which then describes a simple harmonic oscillator in ln a with the k-dependent frequency ω k = βa p kg i ), where g a = g * (T (a)), g i = g * (T (a i )) and β = 0.82 for bosons and β = 0.80 for fermions. The factor of g 1/6 a introduces a slight timedependence into ω k ; we neglect this variation and set g a = g k ≡ g * (T (a k )) when solving Eq. (4.3). The solution is Since ω k encodes the effect of the relativistic pressure of the Y particles, the expression for δ Y for small ω k should match the evolution of cold dark matter in radiation domination [44]: δ X (a) = AΦ 0 ln(Ba/a k ) with A = 9.11 and B = 0.594. The coefficients A 1 and A 2 are determined by evaluating Eq. (4.3) when ω k 1 and matching it to this function. Prior to the EMDE, ln(a/a k ) 10 for modes near the peak scale, and the argument within the sine in Eq. (4.3) is small compared to unity if ω k 1. Using the approximation that sin x x for x 1, it follows that A 1 = AΦ 0 /ω k and A 2 = B.
The peak wavenumber k pk can be found by maximizing the amplitude δ Y (k, a dom ). Using the expressions for A 1 and A 2 , we have where the second equality results from using the expression for a k /a dom in radiation domination from Eq. (B.11). Neglecting the weak k-dependence of (g k /g dom ) 1/6 and Φ 0 while setting the derivative of Eq. (4.4) with respect to k equal to zero implies tan ω pk ln √ 2Bk pk k dom g pk g dom 1 6 = ω pk 1 + ln √ 2Bk pk k dom g pk g dom 1 6 .
Since k pk is an extremum of δ Y (k), the tangent function on the LHS is well-described by a Taylor expansion to second order around k pk . Using this expansion and solving the resulting equation for k pk yields √ 2B k pk k dom g pk g dom where W is the Lambert W-function and We express a i H i /k dom = (a dom /( √ 2a i ))(g dom /g i ) 1/6 (using Eq. (B.11)) and use the definition of a dom /a p from Eq. (B.5) to simplify the dependence of r on the model parameters. Substituting B = 0.594, we have Finally, we use the expression for k y /k dom from Eq. (B.14) to eliminate k dom from Eq. (4.6) and obtain In the above expression, g pk = g * (T (a pk )) and g dom = g * (T dom ). The RHS of this expression includes an additional factor of 1.08 that brings the k pk values into better agreement with those obtained from the numerical solutions of the perturbation equations.
To establish the relation between the peak scale and our model parameters, the peak scale can be rewritten in physical units. Using Eqs. (2.3) and (2.6), and substituting T 0 = 2.726 K and g * S (T 0 ) = 3.91 , k RH ≡ a RH Γ can be expressed as Using the definition of k y /k RH from Eq. (3.8), the peak wavenumber is where we have ignored the variation of g * before the EMDE for simplicity. In Eq. (4.11), k pk /a 0 depends on T RH because the reheat temperature determines when the EMDE ends and thus affects the expansion history of the Universe after the peak scale enters the horizon.
The points in Figure 7 show k pk /k y for different η values as determined from the numerical solutions for the evolution of δ Y , while the solid lines for η > 1 show k pk /k y from Eq. (4.9) with g * (T ) = 100. The analytical expression explains the variation of k pk /k y with η and predicts the peak scale of δ Y (k, a RH ) to within 3% of the numerically determined peak scale for η ≥ 100. As η decreases, the peak scale enters the horizon closer to the pivot point of c 2 sY . Since the asymptotic late-time expression for c 2 sY was used in the derivation of Eq. (4.9), its prediction for k pk diverges from the numerically determined peak wavenumber for η < 100.

Scenarios With Initial Y -Domination
If η < 1, ρ R remains subdominant until reheating. The EMDE begins when ρ Y starts decreasing proportional to a −3 at a = a p , and this pivot also determines which modes are suppressed Figure 7. The peak wavenumber k pk that maximizes δ Y (k), scaled by k y and plotted against η, the initial value of ρ Y /ρ R . The dots represent the k pk determined from numerical data for m = 2 TeV and the lines represent the analytical predictions for k pk /k y , given by Eq. (4.9) for η > 1 and by Eq. (4.13) for η < 1. The plotted quantity k pk /k y depends only on η and the statistics of the Y particles. Different colors show cases with different Y particle statistics. This plot assumes g * (T ) = g * S (T ) = 100. The black crosses show k pk /k y for a boson Y particle for m = 200 GeV; they overlap with the black dots, confirming that k pk /k y is independent of m.
by the relativistic pressure of the Y particles. The numerical solutions to the perturbation equations for η < 0.1 indicate that k pk enters the horizon while the Y particles are still relativistic (a pk < a p ) and that a pk = a p √ 1 + η/γ, where γ = 2.055 and 2.065 for bosonic and fermionic Y particles, respectively. The factor γ accounts for a slight difference between k pk for fermionic and bosonic Y particles, which arises because the Y particles have slightly lower pressure at a given value of a/a p if they are fermions compared to if they are bosons.
For a < a p , H(a) ∝ a −2 and thus k ∝ a −1 k . Therefore, k pk /k p = γ/ To express k p in terms of k y , we again use the scaling k ∝ a −1 k , which applies since a y < a p . Using this scaling yields k y /k p = a p /a y = b, so that k pk k y = γ b This prediction for k pk /k y is shown by the solid lines in Figure 7 for η < 1. The value given by Eq. (4.13) agrees with the peak scale to within 1% for η < 0.1. For η > 0.1, k pk enters the horizon after c 2 sY begins to decrease. This makes k pk /k y diverge from the prediction of Eq. (4.13), which is valid for cases in which c 2 sY (a pk ) 1/3. We can use Eq. (4.12) and k y /k p = b in conjunction with the definition of k RH from Eq. (4.10) to express k pk in physical units as (4.14)

Transfer Functions
Solving the Boltzmann equations with an initially relativistic Y particle is computationally expensive. To facilitate the computation of the matter power spectrum in such hidden-sector cosmologies, we present analytical transfer functions that relate δ Y (k) to δ Y,c (k). The transfer function is defined as 3 , thus obtaining the evolution of density perturbations in cold particles for the same value of ρ Y (a) during the EMDE. If the relativistic Y particles initially dominate the universe, making the Y particles cold radically alters the evolution of the Hubble rate. To avoid conflating the effects of changing the Hubble rate with the effects of the Y particles' pressure, we use an analytical expression for δ Y,c when computing T (k) when η < 1, as described in Sec. 5.2. We focus on transfer functions for the Y density perturbations because they are less sensitive to the duration of the EMDE than transfer functions for DM perturbations would be. In most cases though, these transfer functions can be applied directly to the DM power spectrum. Figure 8 shows the correspondence of T (k) and T X (k) ≡ δ X (k, a RH )/δ X,c (k, a RH ) for k dom /k RH = 17.7 (where δ X,c is the DM perturbation if the Y particles were pressureless). The bottom panel shows the relative error between T and T X , which remains within 5% for T (k) ≥ 0.25. Longer EMDEs lead to even closer agreement between T X (k) and T (k).
We wish to fit a functional form to T (k) that accurately models the transfer function. As can be seen in Figure 9, the oscillatory pattern in δ Y (k) for k > k pk has a much lower amplitude than δ Y (k pk ): δ Y at the second peak scale is 0.25δ Y (k pk ) for both η < 1 and η > 1. Since perturbations at the scale of the first peak will collapse long before modes on smaller scales, we do not expect perturbations with k > k pk to significantly affect the microhalo population. We will therefore prioritize accurately describing the first peak in δ Y (k), while neglecting the smaller peaks at k > k pk . The function that best fits T (k) and accurately describes the first peak in δ Y (k) is where both n and k cut are fitting parameters. Figure 9. Comparing the amplitude of the first and second peaks in δ Y (k). For both η < 1 and η > 1, the second peak amplitude 0.25δ Y (k pk ). Figure 10 shows T (k) for m = 2 TeV, T RH = 20 MeV and η = 500. The transfer function equals unity for k k J,min . As k increases beyond k J,min , T (k) falls off in amplitude as δ Y is increasingly suppressed relative to δ Y,c . After the fall-off, T (k) shows oscillations in k that reflect the small-scale decaying oscillations of δ Y (k) due to the pressure of the Y particles. Figure 10 also shows the transfer function used in B19, given by exp[−k 2 /(2k 2 y )], and we see that our transfer function falls off at comparatively smaller k values. The dashed orange curve in Figure 10 shows the fit to T (k) using the function given by Eq. (5.2), with fit parameters k cut /k RH = 6539 and n = 2.7. The bottom panel shows the percentage error between T (k) and the fitting function. At k = k pk , the value of the fitting function is within 1% of the numerical value of T (k).

Scenarios With Initially Subdominant Y particles
For η between 3 and 1500, the best-fit values for n are between 2.60 and 2.78. Since this variation is small for a range of η that spans nearly two orders of magnitude, we fix n = 2.7 for η > 1. With n fixed, we derive an expression for the cut-off scale k cut by relating it to the peak scale evaluated in Sec. 4.2. For k ≥ 10k dom , the numerical solutions for δ Y,c (k) follow a logarithmic function of a dom /a k , where the mode k enters the horizon at a k : where Eq. (B.11) was used to express a dom /a k in terms of k/k dom . Differentiating δ Y (k) = δ Y,c (k)T (k) with respect to k, while ignoring the weak k-dependence of g * (T (a k )) and using the expression for T (k) given by Eq. (5.2) with n = 2.7 generates a relation between k cut and the peak wavenumber k pk that maximizes δ Y (k): k cut k pk = 2.7 log 0.21 g * (T (a pk )) g * (T dom ) Using Eqs. (4.6) and (4.7) and the expression for k pk /k y from Eq. (4.9), we have The numerical solutions for three values of η are shown in Figure 11 along with the curves given by Eq. (5.2) with n = 2.7 and k cut calculated using Eq. (5.5). The percentage errors between the functional forms and T (k) are plotted in the bottom panel. For η = 1500, which was the maximum value of η for which T (k) was computed, the functional form with n = 2.7 is within 1% of the numerical value of T (k) at k pk . The percentage error remains less than 8% for T (k) ≥ 0.25. The numerically determined k cut is within 4% of the analytical expression given by Eq. (5.5) for η > 50. For smaller values of η, k cut 10k dom , and the expression given by Eq. (5.3) becomes increasingly inaccurate. In addition, our analytical prediction of k pk diverges from the peak wavenumber in the range η < 50. Thus, the prediction of k cut given by Eq. (5.5) becomes inaccurate. For 1.1 ≤ η ≤ 50, we empirically find that a power law describes the variation of k cut /k y with η: This expression predicts the cut-off scale to within 2.5% error for 1.1 ≤ η ≤ 50.
In the η > 1 regime in Figure 12 we show the numerically determined k cut divided by k y as points plotted for different values of η for cases when the Y particles are bosons (red) and fermions (black). For η > 50, the expressions given by Eq. (5.5) are plotted as the solid curves. For 1 < η < 50, the power law fits given by Eq. (5.6) are plotted as the dashed curves. The numerically determined k cut /k y values are shown for m = 2 TeV (dots) and m = 200 GeV (crosses). The overlap of the dots and crosses demonstrates that the validity of the power law fit of Eq. (5.6) is independent of m.

Scenarios With Initial Y -Domination
If η < 1, the universe is initially dominated by the energy density of the Y particles, and δ Y,c (k, a RH ) ∝ a RH /a k ∝ (k/k RH ) 2 for all k k RH . From our numerical solutions, we find (5.7) Here, δ Y,c is evaluated by setting the initial conditions outlined in Appendix C, and Φ 0m = 9Φ 0 /10 is the primordial metric perturbation in a matter-dominated universe. This expression for δ Y,c (k, a RH ) is used to evaluate T (k) when η < 1.
We fit the functional form exp[−(k/k cut ) n ] to the transfer function, treating n and k cut as free parameters. From fitting T (k) for 10 −3 < η < 1, the η-dependence of n can be summarized as The value of n falls with increasing η in the range 0.1 ≤ η < 1 as the contribution of the SM radiation density to the Hubble rate becomes increasingly significant. To find an analytical expression for k cut , we use Eq. (5.7) with T (k) = exp[−(k/k cut ) n ] and maximize δ Y (k) = T (k)δ Y,c (k) with respect to k to obtain the peak wavenumber k pk in terms of k cut : k cut k pk = n 2 1/n .

(5.9)
Substituting k pk /k y from Eq. (4.13), we have For η < 0.1, n = 2.2, which can be substituted in Eq. (5.10) to obtain k cut /k y = 1.04γ/(b √ 1 + η). As η increases from 0.1 to 1, n decreases linearly and k pk /k y increases relative to the prediction of Eq. (4.13). Since the cut-off scale follows the relation k cut ∝ k pk (n/2) 1/n , the rise of k pk /k y nearly cancels out the effect of n decreasing. As a result, predicts k cut to within 2% error even for 0.1 < η < 1. Figure 12 shows k cut /k y in the range η < 1 for the cases when the Y particles are bosons (red) and fermions (black). The solid lines for η < 1 show the predictions of Eq. (5.11). The discontinuity at η = 1 between the analytical predictions given by Eqs. (5.6) and (5.11) arises because Eq. (5.7) for δ Y,c (k) was used to evaluate T (k) for cases with η < 1. Unlike the δ Y,c used in T (k) for η > 1, Eq. (5.7) neglects the contribution of the SM radiation density to the Hubble rate.

The Peak Amplitude and Observational Prospects
Having determined the transfer functions and the cut-off scale, we can estimate how the EMDE impacts the dark matter annihilation rate today. Following the same procedure as Ref. [34], we use the Press-Schechter formalism [45] to obtain the abundance of microhalos, and then we calculate the annihilation rate per volume assuming that the microhalos have an NFW profile with concentration c = 2 at their formation time. The increase of the annihilation rate due to microhalo formation is quantified by the boost factor B(z) ≡ ρ 2 X / ρ X 2 − 1. The resulting boost factor initially increases with time as more halos form, but then it starts to decrease as the earliest-forming microhalos are absorbed into larger halos and their Press-Schechter abundance decreases. The first microhalos are very dense and are expected to survive their absorption into larger halos [46][47][48], so we take the maximum value of B(z) to be the boost factor today (B 0 ). Figure 14 shows B 0 as a function of k dom /k RH for η = 12.69, which corresponds to a hidden sector containing vector Y particles (g = 3) and Dirac fermion X particles (g X = 4) that kinetically decoupled from the SM particles while the Y and X particles were relativistic. The solid lines show B 0 calculated using the transfer functions given by this work for three T RH values, while the dashed lines show B 0 for the same reheat temperatures calculated using the exp[−k 2 /(2k 2 y )] transfer function assumed in Ref. [34]. As the duration of the EMDE increases, earlier structure formation leads to higher boost factors. The growth of density perturbations during the EMDE can even lead to halo formation prior to matter-radiation equality. Following Ref. [34], we only consider halos that form at redshifts less than 10 6 , and  [34] are used. Using our transfer functions increases the value of k dom /k RH that produces a given value of B 0 because our transfer functions suppress perturbations on larger scales, so a longer EMDE is required to generate the same enhancement to the matter power spectrum. this restriction is responsible for the plateau at B 0 2 × 10 18 . The boost increases with increasing T RH at fixed k dom /k RH due to a longer period of logarithmic growth for δ X after the EMDE.
If the EMDE is long enough, microhalos form before reheating. These microhalos are dominated by Y particles, so they dissipate when the Y particles decay [34]. The released dark matter particles have randomly oriented velocities with magnitudes boosted by nonlinear structure formation. This gravitational heating imposes a free-streaming cut-off on the power spectrum after reheating that reduces the boost factor to the standard ΛCDM prediction, B 0 ∼ 10 6 . We implement gravitational heating following the "optimistic" approach from Ref. [34]: a free-streaming cut-off is applied to the matter power spectrum based on the minimum virial velocity of the halos that contain 20% of the dark matter at the end of the EMDE, and no free-streaming cut-off is imposed if less than 20% of the dark matter is bound into halos during the EMDE. The sharp decrease in B 0 due to gravitational heating can be seen at k dom /k RH 100 in Figure 14. The slight red tilt of the primordial power spectrum causes the reduction of the boost due to gravitational heating to move to higher values of k dom /k RH for higher T RH . Figure 14 shows that B 0 for k dom /k RH 80 is smaller for our transfer functions compared to those from Ref. [34]; our transfer functions suppress longer-wavelength perturbations, which reduces the amplitude of the peak in the power spectrum and delays the formation of bound structures. For k dom /k RH 80, the comparatively fewer microhalos predicted by our transfer functions at reheating implies that the reduction of the boost due to gravitational heating happens for larger k dom /k RH (longer EMDEs) compared to when the transfer functions from Ref. [34] are used. Figure 14 demonstrates that the annihilation boost does not strongly depend on T RH , but it is highly sensitive the duration of the EMDE, which sets the peak amplitude of the matter power spectrum for a fixed value of η. While the peak scale controls the size of the first microhalos that form during or after an EMDE, the peak amplitude determines their formation times because gravitational collapse occurs when δ 1.68. The central density of a halo forming at a f scales as a −3 f . Consequently, structures form earlier and have denser cores if the power spectrum has a higher peak [42], which yields larger annihilation boosts up to the point that the peak becomes high enough that halos form during the EMDE.
If the Y particles initially dominate the universe, the amplitude of the peak in the matter power spectrum depends only on the duration of EMDE: a longer EMDE implies a longer period of linear perturbation growth, translating to a higher peak in the power spectrum. If the Y particles are initially subdominant, then the peak amplitude also depends on how long the universe remains radiation dominated after the Y particles become nonrelativistic. Figure  15 demonstrates that a dom /a p depends exclusively on η: a dom /a p remains the same if m is varied while η is held fixed. We use the transfer functions derived in the previous section to calculate δ pk ≡ δ Y (k pk ) and evaluate observational prospects in terms of η and the duration of the EMDE.
If the Y particle is initially subdominant (η > 1), the following fitting function describes δ Y,c (k, a RH ) well for k > 10k RH : where q = k/k dom and p(q) = [1 + 1.11q + (0.94q) 2 + (0.63q) 3 + (0.45q) 4 ] −1/4 . The peak amplitude is δ Y (k pk ) = δ Y,c (k pk )T (k pk ), where T (k) = exp[−(k/k cut ) 2.7 ]. Equation (6.1) shows that δ Y,c (k pk ) is separable into (k dom /k RH ) 2 , which sets the duration of the EMDE, times a function of k pk /k dom . The ratio k pk /k dom depends only on η and the Y particle statistics, as shown by Eq. (4.6). Furthermore, k pk /k cut also depends only on η and the statistics of the Y particles, as illustrated by Eqs. (5.5) and (5.6). The peak scale is given by Eq. (4.9), and Eq. (5.5) provides k cut for η > 50. For 1 < η < 50, the power-law fit of Eq. (5.6) gives k cut , and the peak scale is well-described by the fit where α = 1.82 and b = 2.70 for bosonic Y particles and α = 1.84 and b = 3.15 for fermionic Y particles. This prescription for calculating δ pk matches the numerically determined maximum of δ Y (k) to within 4%. As expected, δ pk ≡ δ Y (k pk ) depends on k dom /k RH and η. Figure 16 shows δ pk /δ Y (k RH ) at a RH for a bosonic Y particle, where δ Y (k RH , a RH ) 3.05Φ 0 is the value of δ Y at the scale at which the power spectrum begins deviating from the power spectrum in scenarios without an EMDE, and we have continued to neglect the scale-dependence of Φ 0 . Since δ Y (k) δ X (k), and δ X (k) only logarithmically increases with k for modes that enter the horizon during radiation domination, this ratio nearly equals the maximum enhancement to the power spectrum. Figure 16 shows that δ pk increases with increasing k dom /k RH , as expected. Figure 16 also shows that δ pk is rather sensitive to η for η 10, but that sensitivity wanes as η increases. As can be seen in Figure 4, δ Y (k) continues to rise steeply with k for k k dom and plateaus when k 10k dom . Therefore, δ pk sharply depends on k pk if k pk 10k dom but then becomes less sensitive to k pk as increasing η increases k pk /k dom . The variation of the peak height with η thus becomes weaker, as indicated by the contours in Figure 16 becoming increasingly vertical as η increases.
If the Y particles dominate the universe before they become nonrelativistic (η < 1), the expression for δ Y,c from Eq. (5.7) implies that where k pk /k RH was split into (k y /k RH )(k pk /k y ). The expressions for k pk /k y and k cut /k y can be taken from Eqs. (4.13) and (5.11) respectively, while n is given by Eq. (5.8). The peak amplitude can be predicted to within 4% error for η < 1 using these expressions.
Equation (6.3) shows that δ pk is proportional to (k y /k RH ) 2 ; this ratio sets the duration of the EMDE. Apart from k y /k RH , the other factors in Eq. (6.3) depend only on η and the Y particle statistics. For η < 0.1, the η-dependence becomes negligible as the universe becomes increasingly Y -dominated at the time of horizon entry of the peak scale. Figure 17 shows the peak enhancement to δ Y at the end of the EMDE as a function of k y /k RH and η for a bosonic Y particle. The contours show that δ pk increases with k y /k RH as the EMDE becomes longer Figure 16. Contours of log 10 [δ pk /δ(k RH )], the maximum enhancement to δ Y at the end of the EMDE, as a function of η and k dom /k RH for cases in which the universe is initially dominated by SM radiation. This plot takes g * (T ) = 100. The red dotted region shows cases in which microhalos are erased at the end of the EMDE and gravitational heating suppresses structure formation after the EMDE. The thick dashed contours mark the largest k dom /k RH values that are compatible with IGRB observations for three values of the DM annihilation cross-section and m X = 10 6 GeV, with the labels showing the σv value considered. For instance, the white region to the left of the dashed contour with log 10 [δ pk /δ(k RH )] = 1.72 is allowed for σv ≥ 10 −9 GeV −2 . The deep blue region is excluded if σv ≥ 10 −15 GeV −2 . The yellow hatched region can be probed with pulsar timing arrays with 100 pulsars observed weekly for 25 years for cases with T RH 20 MeV. and that the enhancement is independent of η for η < 0.1.
For very long EMDEs, B 0 falls below 10 8 due to gravitational heating (as can be seen in Figure 14). All scenarios with higher power spectrum peaks will be similarly affected by the destruction of bound structures during reheating. In Figures 16 and 17, the red dotted regions indicate points with δ pk /δ(k RH ) 10 4.4 , which is the peak enhancement for {η, k dom /k RH } = {12.69, 150}. These areas mark the parameter combinations for which gravitational heating reduces the abundance of microhalos after the EMDE.
If dark matter is a thermal relic, then it is possible to constrain EMDE cosmologies using limits on the dark matter annihilation rate. Reference [34] calculated an annihilation rate per mass in a given volume if the dark matter resides within a population of microhalos: where σv is the velocity-averaged DM annihilation cross-section and M X is the total dark matter mass in the volume. Since the annihilation rate in these scenarios is proportional to the number density of microhalos and thus the dark matter density, the annihilation signal produced is similar to that from decaying dark matter. Assuming that annihilation and decay events both produce two primary particles, the annihilation rate can be related to an effective DM lifetime by equating the rate of particle production from annihilation, 2M X (Γ DM /M X ), to the particle production rate from decaying DM of mass 2m X and lifetime τ , given by (2/τ )[M X /(2m X )]. This yields τ −1 eff = 2m X (Γ DM /M X ), which should be compared to the bounds on the lifetime of DM particles with mass 2m X .
Using Fermi-LAT observations of the Isotropic Gamma Ray Background (IGRB) [49], Ref. [50] established τ eff 10 28 seconds for a variety of decay channels and DM masses ranging from 10 GeV to 10 9 GeV. By connecting B 0 to peak height, we can estimate the allowed regions of the EMDE parameter space based on this lower limit on τ eff . Using Ω X h 2 = 0.12 and m X = 10 6 GeV, the IGRB constraint on τ eff translates to B 0 10 10 for σv = 10 −9 GeV −2 , which is close to the canonical WIMP DM cross-section of 2.2 × 10 −26 cm 3 s −1 [51]. From Figure 14, we see that {η, k dom /k RH } = {12.69, 7} corresponds to B 0 = 10 10 . The white unshaded region to the bottom left in Figure 16 marks all points with δ pk values that are smaller than δ pk for {η, k dom /k RH } = {12.69, 7}, marking the allowed parameter space based on this constraint. The white unshaded region to the left in Figure 17 marks the allowed area of the η-k y /k RH space for cases with η < 1 based on the same peak height limits.
If the DM particle freezes out during or before the EMDE, then a smaller annihilation cross section is required to generate the observed DM density [22,52,53]. Lowering the crosssection increases the upper limit on B 0 , expanding the allowed parameter space. Considering σv = 10 −12 GeV −2 , the light blue regions to the left of the thick dashed contour line with log 10 [δ pk /δ(k RH )] = 2.39 are allowed in addition to the white regions in Figures 16 and 17. With σv = 10 −15 GeV −2 , the allowed region increases to include the medium blue regions to the left of the contour lines with log 10 [δ pk /δ(k RH )] = 2.73, with only the deep blue regions excluded in both plots.
The annihilation contours in Figures 16 and 17 are presented as estimates because they assume that the relation between B 0 and δ pk for a single value of η can be extended to other η values. Power spectra with the same peak heights have different peak scales for different values of η, but the fact that B 0 is largely insensitive to T RH indicates that B 0 does not depend strongly on k pk . Changing η also changes the shape of the power spectrum around its peak. The impact of peak shape on B 0 has not been extensively studied, but the values for B 0 computed in Ref. [35] for a power spectrum with η 1 and the B 0 values for η = 12.69 in Figure 14 differ by less than an order of magnitude for power spectra with the same δ pk . We conclude that the B 0 − δ pk relation for η = 12.69 provides a strong indication of which δ pk values can be ruled out by limits on the dark matter annihilation rate.
It is also possible to detect the microhalos that form after an EMDE through their gravitational influence. Pulsar timing arrays (PTAs) are promising probes of EMDE cosmologies [54,55]: PTAs are sensitive to both the Shapiro time delays as signals pass through microhalos and the Doppler shifts that result when a microhalo pulls on a pulsar, with the latter being most sensitive to sub-earth-mass microhalos [56]. With weekly observations and an RMS timing residual of 10 ns, Ref. [55] showed that microhalos resulting from EMDE-enhanced power spectra with T RH 20 MeV and k cut /k RH > 20 can be detected at 2σ significance if 100 pulsars are observed for 25 years or if 1000 pulsars are observed for 15 years.
Reference [55] used power spectra from initially Y -dominated EMDEs [29] with a Gaussian cut-off. A cut-off given by k cut /k RH = 20 on their power spectra implies that δ(k) peaks at around 24k RH with a value close to 27δ(k RH ). Consequently, power spectra with δ(24k RH ) 27δ(k RH ) will produce microhalos that have similar detection prospects to those produced by initially Y -dominated EMDE scenarios with k cut /k RH 20. Such cases are marked by the yellow hatched regions in Figures 16 and 17; we expect that these EMDE scenarios with T RH 20 MeV will generate signals that are detectable by the PTAs described above.
Another possible method of observing the microhalos resulting from an EMDE comes from how they impact the magnification of stars that pass behind the lensing caustics of galaxy clusters [57][58][59]. As the star passes through the caustic, fluctuations in the dark matter density generate variations in the star's brightness, which can be used to detect sub-earthmass microhalos. Reference [42] identified the ranges of microhalo masses and central densities that can be detected using this method by imposing lower bounds on the magnitude of the observed brightness fluctuations and on the abundance of microhalos. They demonstrated that microhalos that meet their detection criteria are generated by a power spectrum that rises as k 4 for k RH < k < k pk and decreases sharply for k > k pk , reaching a peak enhancement of 10 4 times the ΛCDM power spectrum at k pk ≈ 10 8.5 k eq , where k eq is the horizon scale at matter-radiation equality. Employing our transfer functions, the parameters T RH = 4 MeV with η = 0.01 and k y /k RH = 100 generate a power spectrum with a k 4 rise before a peak at the scale k pk ≈ 10 8.5 k eq . The peak enhancement is δ pk /δ(k RH ) = 126, corresponding to a power spectrum peak enhancement factor of ≈ 10 4 . Using T RH = 4 MeV, η = 8 and k dom /k RH = 20 generates a power spectrum with a similar peak scale and peak enhancement. Although the scaling is not strictly k 4 for k RH < k < k pk in this case, the power spectrum is logarithmic for only a narrow range of k near the peak, making this power spectrum roughly similar to one that rises as k 4 before the peak. For similar peak enhancements, power spectra with a k 4 rise before the peak and larger k pk values from cases with T RH up to a few hundred MeV also result in microhalos that can be detected using caustic microlensing observations [42].

Summary and Discussion
The linear growth of dark matter perturbations during an early matter-dominated era (EMDE) leads to the formation of microhalos much earlier than in standard cosmologies [29,30,32]. These dense microhalos may be detected gravitationally by upcoming pulsar timing arrays [42,54,55] and through their impact on stellar microlensing events in galaxy clusters [42,[57][58][59]. They can also boost the dark matter annihilation rate by several orders of magnitude [32,34,35]. Perturbation growth is suppressed for modes that enter the horizon while the particle that dominates the universe during the EMDE has significant relativistic pressure. The DM power spectrum manifests this suppression as a small-scale cut-off, which strongly affects the DM annihilation signal [32,35]. The small-scale cut-off also impacts the prospects of detecting the structures formed in EMDE cosmologies via pulsar timing arrays [42,54,55] and caustic microlensing [42]. It is therefore important to accurately calculate this cut-off scale, so that EMDE scenarios with initially hot hidden sectors may be tested against observational data. In this paper, we have investigated the small-scale cut-off in the matter power spectrum that results from the relativistic initial state of the particle responsible for the EMDE.
We employed a custom Boltzmann solver to calculate the evolution of perturbations in a universe with an initially relativistic hidden sector particle (Y ). We found that the evolution of subhorizon perturbations in the Y particle density (δ Y ) depends on the wavelength of the perturbation mode compared to a time-varying Jeans length. This Jeans length is set by the sound speed of the Y particles, and it increases while they are relativistic and then starts decreasing after they transition to nonrelativistic behavior. As long as the Jeans length is greater than a perturbation mode's wavelength, δ Y oscillates, while it grows when the Jeans length drops below the mode wavelength. Therefore, linear growth during the EMDE starts later for smaller-scale modes. This suppression of growth due to relativistic pressure generates a peak in the power spectrum of δ Y : for wavelengths smaller than the peak scale, the power spectrum falls off in amplitude due to the delayed onset of growth during the EMDE, whereas longer wavelength modes have less time to grow during the EMDE because they enter the horizon later. This peak is inherited by the dark matter power spectrum as dark matter particles fall into the gravitational wells created by the clustered Y particles during the EMDE.
To describe how the relativistic pressure of the Y particles affects the matter power spectrum, we provided transfer functions that relate the matter perturbations of initially cold and hot hidden sectors. These transfer functions generate the matter power spectrum following an EMDE arising from an initially hot hidden sector without the cumbersome calculation of the density evolution of the hidden sector particle as it transitions from relativistic to non-relativistic behavior. The transfer functions take the form exp[−(k/k cut ) n ], where n depends on ρ SM /ρ Y when the Y particles were relativistic (η) and k cut is the cut-off scale. We found that k cut /k y is a function of η and the Y particle statistics, where k y is the wavenumber of the mode that enters the horizon when the hidden sector temperature equals the Y particle mass m. The ratio k y /k RH , where k RH is defined as the horizon wavenumber at the end of the EMDE, depends on η and is proportional to (m/T RH ) 2/3 . We found that k cut is smaller than k y , which was used as an estimate of the cut-off scale in Ref. [34]. Our result also disproves the claim in Ref. [21] that the horizon scale at the start of the EMDE sets the cut-off scale.
The cut-off scale determines the power spectrum peak height, which sets the formation times and central densities of the first microhalos. The peak height δ pk depends on the EMDE duration and η. Longer EMDEs translate to larger δ pk since they involve longer periods of linear perturbation growth. For η < 1, δ pk ∝ (k y /k RH ) 2 . For η > 1, δ pk ∝ (k dom /k RH ) 2 , where k dom is the horizon wavenumber at the start of the EMDE. If η < 0.1, the peak height is independent of η because the subdominant SM radiation density does not affect the evolution of perturbations prior to the end of the EMDE. For η > 1, the peak height depends on η because η determines how long it takes the Y particle to dominate the universe after it becomes nonrelativistic. Relating the peak height to η and the EMDE duration enables the discussion of observational prospects and constraints in the parameter space of hidden-sector EMDE histories.
If the peak is high enough for microhalos to form during the EMDE, the evaporation of these microhalos at reheating causes the ejection of DM particles at high speeds in random directions. This gravitational heating leads to a free-streaming cut-off on the power spectrum after the EMDE. The exact evolution of this free-streaming cut-off and its relation to the abundance of microhalos that formed during the EMDE is unknown, with recent studies [60] even suggesting that the remnants of evaporated halos may re-collapse into bound structures around the epoch of matter-radiation equality. We identified the regions of parameter space where 20% or more of the dark matter is gravitationally heated; the affected parameter space has peak enhancement δ pk /δ(k RH ) 10 4.4 . This corresponds roughly to cases with η 1/4 k dom /k RH 250 for 1 < η 1000 and k y /k RH 800 for η < 1.
Since the microhalos that form after an EMDE track the dark matter density, the annihilation rate within microhalos can be compared to the rate of particle production from decaying dark matter to define an effective DM lifetime. We used constraints on the dark matter lifetime [50] based on the Fermi-LAT observations of the Isotropic Gamma Ray Background (IGRB) [49] to derive bounds on the dark matter annihilation boost B 0 . By connecting the bounds on B 0 to the peak height, we identified the allowed regions of the parameter space of hidden-sector EMDE histories. Assuming a DM mass of 10 6 GeV with an annihilation crosssection close to the canonical value of 10 −9 GeV −2 , the IGRB constraint allows cases obeying η 1/3 k dom /k RH 25 for η > 1, or cases with k y /k RH 35 for η < 1. Smaller cross-sections are required to match the currently observed DM relic abundance if the DM freezes out during or before an EMDE; the allowed parameter space expands for these lower cross-section values and for higher values of DM mass. Since k cut < k y , our transfer functions yield less structure formation for the same EMDE duration compared to Ref. [34]. For cases not involving gravitational heating, we therefore obtain smaller annihilation boost factors for the same EMDE duration. In addition, this reduced structure formation also delays the onset of gravitational heating, which happens for longer EMDEs compared to Ref. [34].
We also found that a large portion of the parameter space of hidden-sector EMDEs can be probed with the pulsar timing arrays discussed in Ref. [55]. For example, weekly observations of 100 pulsars for 25 years would detect microhalos generated from EMDEs with T RH 20 MeV, 30 k y /k RH 800, and η < 0.1, where the upper limit on k y /k RH comes from the uncertainty associated with the disruption to the post-EMDE power spectrum due to gravitational heating. If η > 1, the same PTA observations would detect microhalos resulting from EMDEs with 13 η 1/4 k dom /k RH 250 and T RH 20 MeV. Furthermore, EMDE power spectra for reheat temperatures less than O(100 MeV) with peaks that are enhanced by a factor of 10 4 relative to the standard ΛCDM power spectrum lead to microhalos that produce detectable brightness fluctuations when stars pass through the lensing caustics of galaxy clusters [42].
Our calculation of the small-scale power spectrum cut-off that results from the relativistic pressure yields a more accurate mapping between the properties of EMDE cosmologies and the observable signals that can help detect or constrain them. Our work thus improves our ability to probe the microscopic properties of hidden sectors and the expansion history of the early universe.

A The Evolution Of The Homogeneous Hidden Sector Background
The Y particles that dominate the energy density of the universe during the EMDE are initially relativistic and transition to a pressureless state as the hidden sector temperature decreases. This appendix presents calculations for the evolution of the equation of state, pressure, and density of the Y particles.

A.1 Method
We use energy conservation and number density conservation to formulate a system of coupled differential equations for quantities related to the hidden sector temperature T hs and the chemical potential of the Y particles, denoted by µ. We will assume here that the Y particle has g degrees of freedom and write the energy density ρ Y , pressure P Y , and number density n Y as thermodynamic integrals: where the ±1 in the denominator denotes fermions (upper sign) or bosons (lower sign). These integrals can be expressed in terms of dimensionless quantities: z ≡ p/T hs , x ≡ m/T hs , We also introduce the notation Conservation of energy density and number density implẏ where w Y (x, ∆) ≡ P Y (x, ∆)/ρ Y (x, ∆) and overdots denote proper time derivatives. Note that we have ignored the decay of the Y particles because an EMDE only occurs when the Y particles transition to nonrelativistic behavior well before they decay. To transform Eqs. (A.4) into differential equations for x and ∆, we expressρ Y andṅ Y in terms ofẋ and∆. For ρ Y , we haveρ And similarly,ṅ where .
Substituting these definitions in Eqs. (A.4) and isolatingẋ and∆ yieldṡ Equations (A.9) are solved to obtain the hidden sector temperature T hs and the chemical potential µ of the Y particles as a function of time. With T hs and µ obtained, the time evolution of the Y particle density and pressure can be calculated using Eqs. (A.1a) and (A.1b) respectively. Finally, the equation of state w Y and the sound speed c 2 sY can be computed;

A.2 Modeling The Transition From Relativistic to Nonrelativistic Behavior
The evolution of ρ Y can be modeled by a broken power law with a pivot point a p . Since the Y particles become nonrelativistic long before their comoving density is altered by their decays, an expression for a p can be obtained by conserving n Y a 3 through the transition from relativistic to nonrelativistic behavior.
We use the ansatz a p /a i = bT hs,i /m, where T hs,i is the hidden sector temperature at a i . Let us assume that the Y particles have become fully nonrelativistic at scale factor a nr . Conserving particle number implies n Y (a i )a 3 i = n Y (a nr )a 3 nr . Since the Y particles are nonrelativistic at a nr , we can write The evolution of w Y and c 2 sY can also be described by broken power laws. Both these quantities are equal to 1/3 when the Y particles are relativistic and are proportional to a −2 when the Y particles become nonrelativistic. This behavior is illustrated for a case with m = 1 TeV and T hs,i = 200m by the blue solid curves in Fig. 18. We find that w Y and c 2 sY are well-described by the functional form where a b is the bending scale factor where the function transitions from the early-time power law to the late-time power law and D models the width of the transition.
Treating a b and D as fit parameters, the numerical solutions for w Y and c 2 sY were fit to the above function for a range of masses and values of T hs,i /m for both boson and fermion Y particles. The best fit values are presented in Table 1. The orange dashed lines in Figure 18 show the functions of the form given by Eq. (A.13) with the best fit values of a b and D given in Table 1. The bottom panels show the relative error between the numerical solution and the best fit functions. The error stays within 1.5% for w Y and 1.2% for c 2 sY and stays within 0.2% at late times for both quantities.  Table 1. Bottom: The relative error between the numerical solution and the fit function, given by 1 − Numerical/Fit.
We also used the fitting function for w Y with the best fit parameters and integrated Eq. (A.4a) to obtain ρ Y to compare it with the numerical solution of ρ Y . The relative error between the numerical and integrated ρ Y peaks at 0.5% and stays constant at 0.2% at late times. Our fitting forms for w Y can be used to obtain the time evolution of ρ Y to within this error. In the derivation of the peak wavenumber in Sec. 4, we also use a piecewise model for c 2 sY , in which c 2 sY is approximated as a sharply broken power law with a pivot point, so that with a pc = 1.43a p for bosonic Y particles and a pc = 1.41a p for fermionic Y particles, where a p is the pivot scale factor for the evolution of ρ Y (a).

B Relating The Start Of The EMDE To Model Parameters
In a universe that is initially dominated by relativistic SM particles, the EMDE starts when the energy density of the Y particles exceeds the energy density of the SM particles at a scale factor a dom , when the SM temperature is T dom . Here, we derive a few important expressions for quantities related to the start of the EMDE in terms of the parameters of our model: m (the Y particle mass), η (the ratio of the initial energy densities of SM radiation and the Y particles), the reheat temperature T RH defined in Eq. (2.3), and b, which is 2.70 and 3.15 if the Y particles are bosons and fermions, respectively. In the following, g denotes the degrees of freedom of the Y particles, and f is 1 or 7/8 for boson or fermion Y particles, respectively.
We first evaluate a dom /a p , where a p is the pivot scale factor for the broken power law followed by ρ Y (a). Since entropy is conserved for the SM radiation before the Y particle decays become significant, we have g * (T )a 3 T 3 = constant, where we assume g * (T ) = g * S (T ). As a result, ρ R ∝ g * (T (a)) −1/3 a −4 . Then, ρ R (a p ) ρ dom = g where g dom = g * (T dom ), g p = g * (T (a p )) and ρ dom = ρ R (a dom ) = ρ Y (a dom ) = (π 2 /30)g * (T dom )T 4 dom . Furthermore, since ρ Y ∝ a −3 from a p to a dom , we can write Similarly, we can express ρ R (a p ) = ρ R (a i )[g i /g p ] 1/3 [a i /a p ] 4 , where g i = g * (T i ). Since ρ Y ∝ a −4 from a i to a p , ρ Y (a p ) = ρ Y (a i )[a i /a p ] 4 = η −1 ρ R (a i )[a i /a p ] 4 . Combining the previous two expressions yields ρ R (a p ) ρ Y (a p ) = η g i g p a dom a p = g i g dom 1 3 η. (B.5) The above relation can be used to express T dom in terms of our model parameters. We can use ρ Y ∝ a −4 from a i to a p and ρ Y ∝ a −3 from a p to a dom to write (B.6) According to our model for the evolution of ρ Y (a), a p /a i = bT hs,i /m and ρ Y (a i ) = (π 2 /30)gf T 4 hs,i . In addition, we use the expression for a dom /a p from Eq. (B.5) and equate ρ Y (a dom ) to (π 2 /30)g * (T dom )T 4 dom to obtain Next, we derive an expression for k dom /k RH ≡ (aH) a dom /(a RH Γ). We divide Eq. (B.5) by Eq. (2.5) to get a dom /a RH , substitute Γ from Eq. (2.3) and use ρ dom from Eq. (B.8) in H(a dom ) = 2(8πG/3)ρ dom to obtain We also find it useful to derive an expression for k/k dom = (aH) a k /(aH) a dom for a mode k that enters the horizon at a k during the period of radiation domination before the EMDE. Since the energy densities of the Y particles and the radiation are equal at a dom , we have H 2 (a dom ) = 2 × (8πG/3)ρ R (a dom ). It follows that a dom a k 2 , (B.10) where g k = g * (T (a k )), and k k dom = 1 √ 2 g dom g k η √ 2 (B.12) for a universe with η > 1.
Finally, we obtain an expression for k y /k dom = a y H(a y )/(a dom H(a dom )) where a y is the scale factor at which T hs = m. For this derivation, we relax the assumption of radiation domination before the EMDE because ρ Y contributes significantly to H(a y ) for η 10. Since the Y particles are relativistic at a y and T hs ∝ a −1 for a i < a < a y , we can express a y = a i T hs,i /m = a p /b. Using Eq. (B.5), this yields a dom a y = g i g dom 1 3 bη. (B.13) Next, we can express H 2 (a y ) = 8πG[ρ Y (a y ) + ρ R (a y )]/3 = 8πGρ R (a y )[1 + η −1 ]/3. Using the scaling ρ R (a) ∝ g(T (a)) −1/3 a −4 with the definition of ρ dom from Eq. (B.8) in H 2 (a dom ) = 16πGρ dom /3, and the expression for a dom /a y from Eq. (B.13), we obtain k y k dom = g 2 i g y g dom (B.14)

C Perturbation Equations
We work in the Newtonian gauge: ds 2 = −(1 + 2ψ)dt 2 + a 2 (t)(1 + 2φ)(dx 2 + dy 2 + dz 2 ). (C.1) Ignoring anisotropic stress, we have ψ = −φ. In the absence of decays, the general equations for the density contrast δ ≡ (ρ −ρ)/ρ and velocity dispersion θ ≡ a∂ i dv i /dt of a fluid for a For a universe that is initially dominated by nonrelativistic Y particles, the initial conditions are similar to those in matter domination. For superhorizon modes, the primordial curvature perturbation is related to the metric perturbation as ζ 0 = 5Φ 0m /3. The primordial curvature perturbations for all species, given by Eq. (C.4), are set equal to each other, yielding The initial conditions for θ i are given by Eq. (C.6) with Φ 0m replacing Φ 0 .

D EMDE Power Spectrum Application
The EMDE modifies the matter power spectrum for modes that enter the horizon during or before the EMDE (k > k RH ). For an EMDE that results from cold Y particles dominating the universe after inflation, this modification to the power spectrum was described by Ref. [29]. For k < 0.05k RH , the power spectrum remains the same. For k > 0.05k RH , δ(k) → R(k)δ(k), where (D.1) In this equation, a k is the scale factor of horizon entry for mode k and a eq and k eq are the scale factor and horizon wavenumber at matter-radiation equality, respectively. The values of f 1 and f 2 are determined by the baryon fraction f b ≡ ρ bar /(ρ b + ρ matter ): a eq a k = √ 2k k eq 1 + k k If an epoch of SM radiation domination precedes the EMDE, modes with k > k dom grow logarithmically with scale factor after entering the horizon and before the EMDE. This modifies δ(k > k dom ). From our fitting function for δ Y,c (k > 10k RH ) given by Eq. (5.7), we find that this modification is modeled by the scale-dependent factor R RD (q) = ln(1 + 0.22q) 0.22q [1 + 1.11q + (0.94q) 2 + (0.63q) 3 + (0.45q) 4 ] −1/4 , (D. 5) where q = k/k dom and R RD = 1 for η < 1. Finally, the small-scale cut-off can be imposed on δ(k) using our transfer functions T (k) = exp[−(k/k cut ) n ] from Section 5. In summary, the combined effect of the EMDE, an epoch of radiation domination before the EMDE, and the small-scale cut-off due to the relativistic pressure of the Y particles modifies δ(k) by a factor R EMD (k) = R(k)R RD (k)T (k). (D.6) The above expression for R EMD (k) is valid at all times after the EMDE ends. We provide an online application for the easy computation and visualization of R EMD 1 . The calculations of the peak and cut-off scales in the application neglect the variation of g * , the number of relativistic degrees of freedom in the SM radiation, before the EMDE. The parameters T RH , η and k dom /k RH or k y /k RH can be varied by the user, and the output is downloadable as a table.