Dark Matter Microhalos in the Solar Neighborhood: Pulsar Timing Signatures of Early Matter Domination

Pulsar timing provides a sensitive probe of small-scale structure. Gravitational perturbations arising from an inhomogeneous environment could manifest as detectable perturbations in the pulsation phase. Consequently, pulsar timing arrays have been proposed as a probe of dark matter substructure on mass scales as small as $10^{-11} M_\odot$. Since the small-scale mass distribution is connected to early-Universe physics, pulsar timing can therefore constrain the thermal history prior to Big Bang nucleosynthesis (BBN), a period that remains largely unprobed. We explore here the prospects for pulsar timing arrays to detect the dark substructure imprinted by a period of early matter domination (EMD) prior to BBN. EMD amplifies density variations, leading to a population of highly dense sub-Earth-mass dark matter microhalos. We use recently developed semianalytic models to characterize the distribution of EMD-induced microhalos, and we evaluate the extent to which the pulsar timing distortions caused by these microhalos can be detected. Broadly, we find that sub-0.1-$\mu$s timing noise residuals are necessary to probe EMD. However, with 10-ns residuals, a pulsar timing array with just 70 pulsars could detect the evidence of an EMD epoch with 20 years of observation time if the reheat temperature is of order 10 MeV. With 40 years of observation time, pulsar timing arrays could probe EMD reheat temperatures as high as 150 MeV.


I. INTRODUCTION
Pulsar timing presents an exquisitely sensitive probe of gravitational perturbations. Due to the highly regular periodicity of many pulsars' emissions, tiny disturbances can manifest as detectable shifts in their pulsation phases. For this reason, pulsar timing is indispensable in the search for lowfrequency gravitational waves, for which pulsar timing arrays (PTAs)-which search for correlated phase shifts-represent the most sensitive detection methodology [1]. However, PTAs can also probe perturbations arising from an inhomogeneous matter environment. In particular, pulsar timing holds the potential to detect primordial black holes [2][3][4] and dark matter subhalos [5][6][7][8][9][10][11][12] as small as 10 −11 M , which lie far beyond the reach of other astrophysical probes of dark matter.
In this article, we explore the potential for pulsar timing to probe an early matter-dominated epoch through its imprint on small-scale dark structure. The boost to density variations causes a large fraction of the dark matter to become bound into highly dense sub-earth-mass microhalos [51]. If a microhalo passes near a pulsar (or the Earth), it perturbs the body's motion and potentially leads to a detectable pulsar timing distortion. Sufficiently compact microhalos that cross the line of sight to a pulsar can also perturb the incoming light directly. Recently, Ref. [11] studied microhalo detection prospects using pulsar timing arrays; they considered EMD-induced microhalos (among others) and concluded that an EMD could be probed if it ends at a reheat temperature T RH below about 1 GeV. Our research builds on the analysis of Ref. [11]. In particular, microhalos detectable by pulsar timing are near the Galactic disk, making them susceptible to disruption by tidal forces and high-speed encounters. Within this context, we employ recently developed semianalytic descriptions of microhalo formation and evolution to precisely characterize the EMD-induced microhalo distribution that manifests within this environment. We also explore the parameter space for early matter-dominated epochs more exhaustively.
To describe microhalo formation, we employ the approach developed in Ref. [52]. This prescription maps peaks in the initial density field onto collapsed halos at later times; the properties of the peak predict the internal structure of the resulting halo. The next step is to model how microhalos respond to the Galactic disk's disruptive environment, and for this purpose we employ the model presented in Ref. [53] to describe subhalo tidal evolution and that presented in Ref. [54] to describe the impact of high-speed encounters. We previously used these models in Ref. [55] to characterize another potential signature of EMD-boosted dark matter annihilation-and we also employ here the model refinements developed for that work. PYTHON codes that implement the microhalo models employed in this work are publicly available. 1 In this paper, we parametrize the impact of an arbitrary EMD scenario (following Ref. [56]) in terms of its reheat temperature T RH and the ratio k cut /k RH between the dark matter free-streaming wavenumber and the wavenumber that enters the horizon at reheating. Both of these parameters are necessary to characterize the imprint of early matter domination: the former sets the scales at which density variations are boosted, while the latter sets the amplitude of the boost. On the observational side we consider the impact of the pulsar count N P and observation duration t obs as well as the rms timing noise residual t res . We find that, roughly, timing residuals t res 0.1 µs are needed to detect EMD epochs. If t res = 10 ns, EMD-induced microhalos could be detected with about N P = 70 pulsars and t obs 20 years if T RH 20 MeV and k cut /k RH 30. Higher T RH and lower k cut /k RH can be reached with larger t obs and N P , although microhalo detection prospects are relatively insensitive to N P . With hundreds to thousands of pulsars, reheat temperatures as high as T RH 150 MeV or cutoff ratios as low as k cut /k RH 8 could be probed in t obs = 40 years.
This article is organized as follows. Section II describes the matter power spectrum that arises from an EMD scenario. In Sec. III, we characterize the microhalo distribution resulting from such a power spectrum. In Sec. IV, we treat the evolution of microhalos due to tidal forces and encounters within the Galactic potential. Section V outlines our method for assessing the detection prospects of a given microhalo population using pulsar timing, which is largely the same as the method presented in Ref. [11]. These detection prospects, which constitute the main results of this article, are presented in Sec. VI. Section VII presents our conclusions. Finally, Appendix A further details Sec. IV's computation of the microhalo-stellar encounter distribution, while Appendix B refines the model in Ref. [54] to better describe a halo's long-term response to a high-speed encounter.

II. DENSITY VARIATIONS ARISING FROM EMD
We first review the spectrum of density variations that results from an early matter-dominated phase. During EMD, the dominant species, which we call φ, gravitationally clusters and creates gravitational potential wells. If the dark matter is both nonrelativistic and kinetically decoupled from any relativistic species (such as the Standard Model) by a temperature of about 2T RH [55], then it falls into these potential wells and inherits the density variations in φ. After reheating, these density variations persist in the dark matter.
Reference [51] developed a prescription to quantify the power spectrum P(k) of density fluctuations imprinted onto the dark matter by a general EMD epoch. The form of P(k) is influenced by two parameters: the reheat temperature T RH 1 https://github.com/delos/microhalo-models EMD-induced (dimensionless) dark matter power spectra at z = 300, computed using linear theory. We show different values of the reheat temperature TRH and the dark matter free-streaming cutoff scale kcut, expressed as the ratio kcut/kRH. TRH sets the scales at which EMD boosts density variations, while kcut/kRH sets the maximum amplitude of the boost. For comparison, the solid black curve indicates the cold dark matter power spectrum without early matter domination. We also show the mass scale M associated with each wavenumber, which roughly corresponds to the masses of halos forming from density variations of that scale. For kcut/kRH 20, density variations are already nonlinear (P(k) 1), implying microhalos are already forming by z = 300. associated with the end of EMD; and the dark matter freestreaming scale, which sets a cutoff wavenumber k cut . 2 Figure 1 shows the dark matter power spectrum predicted by this prescription for several different EMD scenarios. These power spectra are computed in linear theory, i.e., assuming P(k) 1.
A key observation is that if the ratio k cut /k RH is larger than about 20 between k cut and the wavenumber, k RH , which enters the horizon at reheating, then density variations are already becoming nonlinear (P 1) by the redshift z = 300 at which the power spectra are plotted. Consequently, collapsed halos are already beginning to form at this early time. We also emphasize the fundamental importance of the cutoff scale k cut to the EMD epoch's imprint on density variations: if k cut < k RH , then free streaming erases any trace of EMD from the dark matter distribution.
We focus in this study on the reheat-temperature range 3 MeV < T RH < 180 MeV. According to Ref. [14], the effective number of neutrino species reflected in the CMB constrains T RH > 4.7 MeV, which motivates our lower bound. Meanwhile, high reheat temperatures T RH are associated with smaller mass scales (see Fig. 1), which are more difficult to detect through pulsar timing, and we will see later that re-heat temperatures T RH 180 MeV are broadly inaccessible. With respect to the cutoff parameter k cut /k RH , 3 we limit our consideration to 5 < k cut /k RH < 40. We will see that microhalos arising in scenarios with k cut /k RH < 5 are not dense enough to be detected. On the other hand, larger values k cut /k RH > 40 lead to the collapse of overdense regions (and likely halo formation [58]) deep in the radiation-dominated epoch, a scenario for which our modeling of halo populations is not calibrated. Nevertheless, as we discuss in Sec. VI, our conclusions are likely applicable to many scenarios with k cut /k RH > 40.

III. THE INITIAL MICROHALO DISTRIBUTION
To understand the present-day imprint of EMD, the next step is to advance from the linear-theory power spectrum to the population of (nonlinear) collapsed dark matter halos. To model halo populations in this way, it is common to use Press-Schechter theory [59,60]. In this model, a halo of mass M is associated with a region of Lagrangian space (the space of initial comoving particle positions) of mass M inside which the average linear-theory density contrast δ(x) ≡ δρ/ρ exceeds some critical value of order 1. Mathematically, the linear-theory density field δ(x) is convolved with a filter that smooths it over the mass scale M .
However, the smallest dark matter microhalos form directly from the collapse of peaks in the linear-theory density field δ(x). To study them, it is more natural to consider the unfiltered density field and simply map each peak therein to a collapsed halo. This procedure carries the added benefit that each halo's density profile can be predicted from the properties of the associated peak [52], so we do not need to rely on empirical concentration-mass relations (e.g., Ref. [61]) that can be difficult to generalize to arbitrary power spectra.
Therefore we choose to characterize the microhalo population using the peak-based prescription that Ref. [52] introduced. 4 To begin, we exploit the statistics of Gaussian fields [62] to generate a random sample of peaks along with their density profiles δ(q), where q is the comoving radius. With these peak density profiles, we estimate the density profiles of the resulting collapsed halos using a secondary infall model [63,64]. Specifically, we use the simulation-tuned 3 Note that the two scales k RH and kcut are naturally connected despite being associated with entirely different physics (k RH with the φ decay rate and kcut with dark matter's interaction properties). To achieve the observed dark matter abundance, a thermal relic must thermally decouple from the radiation bath not too long before reheating; see Ref. [55]. The free-streaming scale kcut is set by the the time of kinetic decoupling, which is in turn related to the time of thermal decoupling through the particular dark matter microphysics. Further details of the kcut-k RH connection are laid out in Refs. [56,57], the latter of which suggests that kcut/k RH as high as 200 can be plausible for supersymmetric dark matter. However, the presence of a radiation-dominated epoch prior to EMD can break the kcut-k RH connection and allow kcut k RH (see, e.g., Ref. [58]). 4 This prescription was also employed in the context of EMD by Ref. [55] to study a different potential observable signature-enhanced annihilation radiation. "turnaround" model of Ref. [52] (refined as described in Appendix A of Ref. [55]) to predict the radius r max at which the circular velocity is maximized and the corresponding enclosed mass M max . Some peaks do not have a predicted collapse time due to their high ellipticity, and we neglect these peaks. 5 Additionally, some peaks produce halos whose predicted r max values are so large as to be populated by material initially lying well beyond the scales boosted by early matter domination. These predictions represent conventional CDM halos unboosted by EMD, so they are irrelevant to pulsar timing observations (e.g., Ref. [11]), and we neglect them for computational convenience. 6 While microhalos form with density profiles that asymptote to ρ ∝ r −3/2 at small radii [52,[66][67][68][69][70][71][72][73][74], successive mergers drive their inner density cusps toward a shallower ρ ∝ r −1 scaling [52,[73][74][75][76]. Therefore, we assume microhalos develop the Navarro-Frenk-White (NFW) profile [77,78], where the scale parameters r s and ρ s are set to reproduce the predicted r max and M max . No precise description of the impact of microhalo mergers has yet been presented; the above procedure yields a conservative estimate, as mergers between microhalos are understood to raise r max and M max [52]. The above prescription sets the distribution of microhalo density profiles. We set their total number density in the following way. The cosmological microhalo number densityn is simply the number density of initial density peaks (which can be computed from P(k) using, e.g., the methods of Ref. [62]) scaled to account for the aforementioned removal of some peaks from consideration. The number density of microhalos in the Sun's vicinity is then n = (ρ local /ρ)n, where (ρ local /ρ) 3 × 10 5 is the factor by which the local dark matter density exceeds the cosmological mean.
These procedures do not fully account for the impact of microhalo-microhalo mergers, a topic that we leave for future study. These mergers tend to raise microhalo masses, as noted above, while also reducing their number count. While these two effects have opposite implications on observational prospects, we argue that neglecting both of them represents a conservative choice. Mergers have two relevant effects: (1) Mergers boost individual microhalo masses while preserving their overall mass content (i.e., the fraction of total dark matter mass residing in microhalos).
Reference [11], exploring constraints from pulsar timing on monochromatic halo mass functions, found that when condition (2) holds-concentrations are held fixed-the parameter change (1)-increasing the halo mass while preserving the fraction of total dark matter mass residing in microhaloscan only strengthen pulsar timing constraints. Consequently, mergers between microhalos can only boost their detectability in pulsar timing. We also note that the preservation of microhalo internal densities implies that their susceptibility to disruptive effects within the Galaxy, which we treat in the next section, is unaltered by mergers. Figures 2 and 3 show (as dashed lines) the initial microhalo distributions, computed as described above, for two different EMD scenarios. We show the distributions in NFW scale density ρ NFW s and halo mass M , where we take M to be the virial mass at z = 2, roughly the time at which microhalos might be expected to accrete onto the Galaxy. 7 The bulk of the microhalos form from density variations that are just above the freestreaming scale, and these halos comprise the main cluster in ρ − M space. Beyond this, there is a tail of larger and lessdense halos that are associated with increasingly large-scale density variations.

IV. DISRUPTION OF MICROHALOS
In isolation, microhalos arising from EMD maintain static internal density profiles, only growing outward as they accrete material. In that case, the above predictions would remain accurate today. However, microhalos relevant to pulsar timing are in the vicinity of the Galactic disk at the present time. These microhalos accreted onto the Galactic halo and were subjected to evolution due to tidal forces and encounters with other objects. We now discuss how we treat these effects.

A. Tidal evolution
A common treatment for tidal evolution (e.g., Ref. [11]) is to simply truncate a subhalo's density profile at its "tidal radius", the radius beyond which tidal forces from the host are stronger than the subhalo's self-gravity. While this procedure supplies a simple approximation, in reality tidal stripping is a continuous process as the subhalo revirializes in response to mass loss (e.g., Ref. [80]). Instead, we use the model from Appendix E of Ref. [53], which predicts the time evolution of r max and M max for an NFW subhalo orbiting an NFW host and is tuned to match idealized N -body simulations. The Galactic potential is not NFW, owing to the contribution of baryons, so as an approximation, we pick the NFW potential that best reproduces the observationally modeled Galactic rotation curves in Ref. [81] for Galactocentric radii smaller than 25 kpc. This NFW potential has scale radius r s = 3.4 kpc and ρ s = 4.0 × 10 8 M /kpc 3 . 8 This choice of host potential represents an approximation, and we note the following limitations. Galactic rotation curves are valid within the disk plane, but the dark matter halo-and hence microhalo orbits-typically extend far above and below this plane, where the potential is weaker. In this respect our choice is likely to overestimate the impact of tidal stripping. Meanwhile, abrupt features within the host potential heat subhalos, increasing the rate at which they lose mass. While we account for encounters with individual disk stars below, the encounter with the Galactic disk's bulk potential can itself be a significant factor in microhalo disruption. This concern suggests we may risk underestimating the impact of tidal stripping. However, Ref. [82] estimates that the impact of individual stellar encounters dominates over that of coherent disk shocking for halos smaller than earth mass. As Fig. 1 indicates, EMD-induced microhalos are sub-Earth mass if T RH 10 MeV. In any event, since our tidal evolution model is only calibrated for NFW host potentials, we leave to future work more precise modeling of microhalo survival within the Galactic potential.
We now sample the orbital parameters of microhalos at the solar radius r 0 = 8 kpc under the approximation that the microhalo velocity distribution is isotropic. Specifically, we use rejection methods to sample the orbital energy E from the un- is the phase-space distribution function and Φ is the Galactic potential. A fitting form for f (E) is drawn from Ref. [83]. The microhalo velocity's magnitude |V | immediately follows. To describe the direction of the microhalo's velocity V , we uniformly sample µ ≡ cos θ, where θ is the angle between V and the Galactic radius vector, and the angle φ between the microhalo's tangential velocity projection and the velocity V LSR of the local standard of rest. Since faster microhalos are encountered more frequently, we additionally use rejection methods to weight our orbit sample by V rel ≡ |V − V LSR |. 9 In this way, we associate a randomly sampled orbit with each initial microhalo from Sec. III. We now apply the tidal evolution model in Ref. [53] assuming a duration of 10 Gyr, which corresponds to accretion onto the Galaxy at a redshift of about z = 2. This approximation does not have a major impact; large changes in the microhalo accretion redshift correspond to only modest changes in the tidal evolution duration. and M separately. In all plots, the dashed curves show the initial microhalo population before disruptive effects are accounted for (equivalently, before accretion onto a larger halo). The bulk of these halos are tightly clustered in ρ NFW s -M space because they form from the smallest-scale density variations above the free-streaming cutoff, but a tail of halos forming from larger-scale fluctuations is visible. Tidal forces from a host halo strip off the microhalos' weakly bound outskirts, resulting in halo populations with similar ρ NFW s but significantly reduced M (dotted lines). Also including the impact of stellar encounters greatly spreads out the distribution (solid curves), producing a long tail of "unlucky" low-mass halos. Note that stellar encounters raise ρ NFW s only in the sense that stripping off a halo's low-density outskirts raises its average density. Stellar encounters still reduce a halo's density at any given radius.

B. Stellar encounters
The next step is to consider encounters with other objects, such as stars or other microhalos. Reference [55] found that the impact of encounters with other microhalos was subdominant to tidal evolution and stellar encounters even for the dark matter-dominated Draco dwarf spheroidal galaxy, so we can reasonably assume it is even more subdominant for microhalos in the solar neighborhood. However, stellar encounters are likely to have a significant, and perhaps even dominant, impact on microhalo evolution.
To handle stellar encounters, a common approximation (e.g., Ref. [85]) is to compare the energy injected into a microhalo by these encounters to the halo's total binding energy. However, as noted by Ref. [80], this computation is not clearly connected to the question of halo survival because these energy injections are inefficient: the most weakly bound particles receive the most energy in an encounter. Instead, we use the model in Ref. [54] to treat the impact of stellar encounters; this model precisely predicts the evolution of a halo's density profile due to impulsive point-particle encounters. To apply this model it is necessary that for each microhalo we sample a series of stellar encounters, and we do so as follows. Let n * (R) be the number density of stars at the galactocentric The initial halo population (dashed lines) exhibits coherently lower halo masses and higher internal densities but is otherwise the same as Fig. 2. However, the higher internal halo densities make this population less susceptible to disruptive tidal effects and stellar encounters (dotted and solid lines).
position R, and suppose a microhalo's orbital position and velocity are described by R(t) and V (t), respectively. The differential number of stellar encounters, per impact parameter b and time t, is whereV and f (R; V * ) is the distribution of stellar velocities V * at position R. Appendix A shows how Eq. (3) can be rapidly evaluated.
We use the Galactic model in Ref. [86] to describe the distributions, n * (R) and f (R; V * ), of stellar positions and velocities, but for simplicity we cylindrically symmetrize the Galactic bar and assume that the stellar velocity dispersion is isotropic. 10 With these distributions established, we use inverse transforms to sample stellar encounter times t and impact parameters b from Eq. (2). We then sample the relative velocities |V * − V | of these encounters using the microhalo's orbital information and the stellar velocity dispersions. Finally, we sample the masses of the stars that our microhalos encounter from a Kroupa initial mass function with minimum mass 0.01 M (we include brown dwarfs) and a high-mass index of 2.7 [87].
The above considerations determine the stellar encounter distribution for each microhalo, and the next step is to apply the stellar encounter model of Ref. [54]. We additionally fol-low the suggestions of Ref. [55], which carried out combined tidal evolution and stellar encounter simulations. To wit, we do not model any relaxation of microhalos between encounters, instead simply summing the (relative) energy injections from all encounters, and we apply the full sequence of stellar encounters after the tidal evolution. The tidal evolution model predicts for each halo its maximum circular velocity v max and associated radius r max , and we relate these quantities back to a density profile by assuming that the profile after tidal evolution resembles the profile after stellar encounters, which, as discussed in Ref. [54], is (almost) universal. We then apply the stellar encounter model with this density profile as a starting point. The post-encounter density profile given in Ref. [54] is ρ = ρ s (r s /r) exp [(−1/α) (r/r s ) α ] with α = 0.78, but as we discuss in Appendix B, this density profile becomes unphysically steep at large radii. We instead propose the form β , α = 0.78, and β = 5. Here, 2 F 1 is the hypergeometric function. This profile transitions to ρ ∝ r −4 at large radii instead of becoming arbitrarily steep, and it agrees with the simulation we carry out in Appendix B. Figures 2 and 3 show the impact of tidal effects (dotted lines) and stellar encounters (solid lines) on the microhalo population for T RH = 32 MeV and two different freestreaming scales. The post-encounter density profile, Eq. (4), has finite total mass, so we do not assume any radial truncation in determining each halo's mass. Additionally, to fairly compare characteristic density values between the disrupted and initial halo populations, we convert the post-encounter density profile's scale density ρ s into the equivalent NFW scale density ρ NFW s = ρ s /1.17 (see Ref. [54]). Microhalos evidently suffer a large drop in mass, which can be attributed to the loss of the halo's weakly bound outskirts. Much of this mass loss actually reflects the shift in density profile from NFW to Eq. (4), but highly dense microhalos might not fully transition their density profiles, so we likely overestimate microhalo disruption when k cut /k RH is large. 11 Meanwhile, the characteristic density ρ NFW s of microhalos increases because when the low-density outskirts are stripped, the average density of the halo rises. 12 This change is relatively modest and arises predominantly due to stellar encounters and not tidal stripping. 11 In particular, in the kcut/k RH = 30 case the tidal-stripping model predicts essentially no evolution in the scale parameters of the microhalo density profiles, so almost all of the change between the dashed and dotted lines in Fig. 3 arises from the change in the density profile's form. Thus, it is likely that for kcut/k RH 30 the assumption that density profiles fully transition into the new form represents an overestimate. However, it should be noted that the total microhalo mass M does not enter into the pulsartiming analysis in Sec. V except for halos that remain very distant from the pulsar. The pulsar-timing signal is dominated by the closest encounters [11], for which only the halo mass interior to the pulsar's separation from the halo is relevant. 12 As can be seen in Ref. [54], the density at any given radius within a halo always drops due to stellar encounters even as the halo's characteristic density rises.
The stochastic nature of stellar encounters also spreads out the microhalo distribution, producing a long tail of low-mass halos that suffered more disruption by chance.

V. PULSAR TIMING SIGNALS
With microhalo populations established, the remaining step is to determine the extent to which they can be detected by PTAs. In particular, we consider here the timing distortions that arise from the Doppler effect when microhalos perturb pulsar motions. While there is also a Doppler effect associated with perturbations to the Earth's motion, Ref. [11] found that this effect generally produces weaker constraints (although it has the potential to probe smaller masses). Timing distortions could also arise from the Shapiro effect when microhalos cross the line of sight to a pulsar, but this effect is most sensitive to mass scales M 10 −5 M [10], which are larger than the masses of EMD-induced microhalos (see Fig. 1).
We follow the procedure of Ref. [11] and use a modified version of the associated Monte Carlo simulation code. 13 This code prepares a random array of pulsar positions, and for each pulsar, it randomly samples encounters with microhalos. Pulsar velocities (relative to the local standard of rest) are neglected because they are much smaller than microhalo velocities (see footnote 9). We modified the code to sample microhalos from the distribution computed in Secs. III and IV. The remainder of this section constitutes a review of the computational procedure described in Ref. [11] through which we evaluate the the pulsar timing signal that arises from such a microhalo distribution.
A pulsar's Doppler shift due to microhalo encounters manifests into an accumulated shift δφ in the pulse phase. For a pulsar with frequency ν and earth-pulsar direction vectord, the phase shift due to an encounter with a point object (like a primordial black hole) with mass M , velocity v, and impact parameter b, is (where c is the speed of light). Here, x ≡ (b/v)(t −t), wherẽ t is the time of closest approach when the vector from the pulsar to the object is b. The phase shift due to an extended object must in general be computed numerically, but to avoid the computational expense, we follow Ref. [11] and conservatively approximate the phase shift due to a microhalo encounter as δφ δφ pt F(r min ).
Here, r min ≥ b is the microhalo's closest approach to the pulsar during the observing time, and the "form factor" 13 The original code resides at the URL https://github.com/ szehiml/dm-pta-mc, and our modified code is located at https: //github.com/delos/dm-pta-mc. where M (r) is the enclosed mass profile, is the factor by which the microhalo's finite extent scales its gravitational field at the distance r, relative to a point object.
We sum the microhalo-induced phase shift over all of the randomly sampled microhalos in the pulsar's neighborhood. Next, we subtract a quadratic fit, φ fit (t) = φ 0 + φ 1 t + φ 2 t 2 , reflecting that these low-order terms are degenerate with the pulsar's intrinsic behavior. 14 We plot in Fig. 4 an example phase-shift signal from a single pulsar due to a microhalo distribution arising from the EMD scenario with T RH = 32 MeV and k cut /k RH = 30. A range of observation durations t obs are shown (for the same microhalo encounters); the signal is different in each case due to the subtraction of the quadratic fit. While these phase shifts δφ are evaluated using the approximate expression in Eq. (6), we also compare the exact phase shift computed using numerical integration (dotted black curve, nearly overlapping the solid black curve). The difference is marginal, which suggests that the approximation does not have a large impact.
We repeat the above calculation for all N P pulsars under observation. Additionally, to account for stochasticity, we repeat this process for 1000 "universes" each with independently sampled microhalo distributions. From each pulsar's 14 In particular, we follow Refs. [9][10][11] in assuming thatν and higher derivatives of the pulsar frequency arise solely due to perturbations by substructure. In practice contributions toν can arise from the pulsar's intrinsic spindown, its unperturbed line-of-sight velocity, or a combination of the two; these contributions can be estimated and are of orderν/ν ∼ 10 −31 s −2 [88].
phase-shift signal, we compute the signal-to-noise ratio (SNR) as described in Ref. [11] assuming optimal filtering, and the SNR associated with each "universe" is taken to be the largest SNR associated with any pulsar therein. The overall SNR is then the tenth percentile of the SNRs of the 1000 "universes". Finally, we evaluate the significance level σ associated with that SNR. That is, σ is the significance level with which we can exclude the possibility that the given SNR arose from noise alone. These choices are identical to those made in Ref. [11].

VI. RESULTS AND DISCUSSION
We initially fix the observational cadence to 1 week and the rms timing residual to 10 ns, and we assume that the timing noise is white. These are the choices made in Ref. [11], but we discuss their impact and validity later. However, we vary the observing time t obs and the number of pulsars N P under observation, and we also vary the cosmology parameters T RH and k cut /k RH . We thus have a four-dimensional parameter space, and each point in this space is associated with a particular pulsar timing SNR and hence statistical significance, as described in Sec. V. From a more practical point of view, for any given pulsar count N P and observation duration t obs , we may target a minimum significance level for a pulsar timing signal. There is a two-dimensional subspace of cosmological parameters T RH and k cut /k RH that satisfy this constraint and hence can be probed by the given observational parameters.

A. Impact of pulsar count and observation time
We first explore the impact of the observation duration t obs and pulsar count N P on the range of EMD scenarios that can be probed. We begin by varying N P with fixed t obs = 20 years. The upper panel of Fig. 5 shows, as colors, the number of pulsars N P needed to detect microhalos arising from the given T RH and k cut /k RH at 2σ significance. With this t obs , the detection of microhalos arising from the EMD scenario with T RH = 3 MeV and k cut /k RH = 40 requires at least N P 70 pulsars, a number comparable to the count in existing PTAs (e.g., Ref. [89]). Meanwhile, as N P is raised, higher reheat temperatures T RH and lower cutoff ratios k cut /k RH can be probed. This trend is easy to understand: higher T RH leads to microhalos of lower mass, while lower k cut /k RH leads to microhalos of lower internal density; both of these trends make the microhalos more difficult to detect. Scenarios with T RH 70 MeV or k cut /k RH 13 require more than 2000 pulsars if t obs = 20 years.
We next explore what observation duration t obs is necessary if we are able to observe N P = 200 pulsars. The lower panel of Fig. 5 shows, again as colors, the t obs needed to detect microhalos arising from the EMD scenario with the given T RH and k cut /k RH at 2σ significance. With this pulsar count, at least t obs 15 years (black) are required to detect microhalos arising from the EMD scenario with T RH = 3 MeV and k cut /k RH = 40. Similarly to N P , raising t obs allows higher T RH and lower k cut /k RH to be probed. EMD scenarios with k cut /k RH 10 or T RH 100 MeV require more than 40 years of observation time if N P = 200. Naïvely one might expect N P and t obs to have a similar quantitative impact, because the effective spatial volume that pulsar timing probes is proportional to both. That is, the number of microhalos passing near any given pulsar is proportional to t obs , so the total number of microhalos passing near pulsars is proportional N P t obs . However, by comparing the two panels of Fig. 5 we notice that raising the observation duration t obs has a much greater impact than raising the pulsar count N P on the range of EMD scenarios that can be probed. For instance, raising t obs from 20 years to 30 years has about as much impact as boosting N P by a factor of 10 from 200 to 2000.
The reason for observation time's power is that beyond boosting the number of microhalo encounters, t obs has another advantage that is illustrated in Fig. 4. This figure suggests that raising t obs increases not only the duration of the signal but also its amplitude. The mathematical reason for this effect is the subtraction of the quadratic fit. Physically, more observation time allows us to better distinguish the impact of microhalos-which induce perturbations to the phase shift δφ at higher-than-quadratic order-from the pulsar's in-trinsic quadratic-order phase shift. 15 Put another way, the time scale associated with microhalo encounters is of order years (see Fig. 4), so tens of years of observation time are needed to fully resolve them.
B. Detection prospects for EMD reheat temperatures and cutoff ratios Figure 5 indicates that when k cut /k RH 30, the detection prospects of an EMD scenario are nearly independent of k cut /k RH , and when when T RH 20 MeV, the detection prospects are nearly independent of T RH . Both of these trends can be understood in light of Fig. 2 of Ref. [11]. Low T RH corresponds to large microhalo mass scales, and this figure indicates that when the mass scale of microhalos exceeds roughly 10 −7 M , their detection prospects become mass independent. Meanwhile, high k cut /k RH leads to microhalos of high internal density, which means that these halos are highly centrally concentrated. In this case we approach the pointmass limit, in which microhalo detection prospects become independent of their level of central concentration.
Thus, in these regimes we can express EMD detection prospects as a function of T RH alone and k cut /k RH alone, respectively. Accordingly, the upper panel of Fig. 6 shows the number N P of pulsars and observation duration t obs needed to detect microhalos (at 2σ significance) arising from an EMD scenario with a given T RH in the case where k cut /k RH 30. 16 Evidently, microhalos arising from reheat temperatures as high as T RH 170 MeV can be detected if N P = 2000 and t obs = 40 years. We also remark that to the extent that a direct comparison can be made, our results are similar to the EMD detection prospects presented in Ref. [11]. 17 Similarly, the lower panel of Fig. 6 shows the N P and t obs required to detect microhalos (at 2σ significance) arising from an EMD scenario with a given k cut /k RH if T RH 20 MeV. In this case, microhalos associated with cutoff ratios as low as k cut /k RH 8 can be detected if N P = 2000 and t obs = 40 years. 15 In practice, long observation durations can also increase the impact of timing noise if the timing residuals are correlated over long time periods. For instance, the contribution from any intrinsic frequency second derivative, ν, could lead to such correlations. The analysis here does not reflect this possibility because we follow Ref. [11] in assuming, for simplicity, that the pulsars' intrinsic timing noise is white. 16 For kcut/k RH 40, microhalos can form significantly before the onset of the last matter-dominated epoch. Halos can form during radiation domination and would be even denser than those that form later [58], so our conclusions still apply in this case. However, if kcut/k RH is sufficiently large that halos form during the EMD epoch, their evaporation at reheating can suppress small-scale structure [58,90]. Further analysis is needed to determine microhalo detection prospects in these cosmologies. 17 For instance, in Fig. 8 of Ref. [11] the observation duration is fixed at t obs = 30 years, and when N P = 1000, 2σ significance is achieved when T RH is slightly over 100 MeV, a result that is also evident in our

C. Impact of timing noise and cadence
Up to this point, we have fixed the observational cadence at t cad = 1 week and the rms timing noise residual at t res = 10 ns. This is an optimistic timing-noise assumption; timing noise residuals in present PTAs are on the order of µs or greater (e.g., Ref. [89]), although precision modeling can mitigate this noise (e.g., Ref. [91]). In Fig. 7 we test how EMD detection prospects change if t res is increased. This figure shows, for t res ranging from 10 to 50 ns, the maximum T RH and minimum k cut /k RH detectable for given t obs and N P . Evidently, the rms timing residual t res is a highly important variable, and t res 0.1 µs is necessary to detect microhalos arising from EMD using pulsar timing if t obs 40 years and N P 800.
Under the assumptions made in Sec. V (and Ref. [11]), decreasing the observational cadence t cad can counteract the impact of timing noise because the signal-to-noise ratio scales as SNR ∝ t −1 res t −1/2 cad (see Eq. 14 of Ref. [11]). That is, since the noise is assumed to be white (uncorrelated), its impact can be reduced to an arbitrary extent by simply taking more measurements. In practice, however, a significant contribution to the timing noise is red noise associated with the pulsar's spin (e.g., Refs. [92,93]), which exhibits temporal correlations. For a red noise spectrum there is less advantage to raising the observational cadence. We also remark that the timescale associated with a typical microhalo encounter is about a year if T RH ∼ 200 MeV and is longer for lower T RH . Any cadence t cad much shorter than a year is sufficient to temporally resolve these encounters. Consequently, the cadence plays no other role in our calculation than to control the noise level per the above discussion.

VII. CONCLUSION
In this article, we explored the prospects for pulsar timing arrays to probe early matter-dominated eras through their impact on small-scale dark structure. We found that in order to detect EMD-induced dark matter microhalos, pulsar timing noise residuals must be brought below roughly the 0.1 µs level. Additionally, due to the long timescale associated with microhalo-pulsar encounters (of order years), tens of years of observing time are necessary. However, with 10 ns timing residuals, EMD scenarios can begin to be probed with about 20 years of observation time and as few as 70 pulsars. This is roughly the pulsar count in existing arrays [89]. With 40 years of observing time and hundreds to thousands of pulsars, EMD reheat temperatures T RH up to approximately 150 MeV can be reached. Note that detection prospects are only weakly sensitive to the pulsar count. Figure 6 shows our main results. In addition to the reheat temperature T RH , which sets the mass scale for EMDinduced microhalos, we also explore the detection prospects with respect to the dark matter's free streaming scale k cut . Specifically we consider the ratio k cut /k RH , where k RH is the wavenumber that enters the horizon at reheating. The combination k cut /k RH is fundamentally important to the microhalo population induced by early matter domination because it sets microhalos' internal density scale (see Figs. 1-3). For k cut /k RH 30 microhalos are close to the point-mass limit, for the purpose of pulsar timing distortions caused by the Doppler effect. In this case the detectable range of T RH is insensitive to k cut /k RH . For smaller k cut /k RH microhalos' spatial extent matters, and k cut /k RH as small as about 8 can be probed given again hundreds to thousands of pulsars and 40 years of observing time (see Fig. 6).
This work refines the EMD detection prospects presented in Ref. [11]. The major new feature of our analysis is the use of the recently developed semianalytic microhalo models put forth in Refs. [52][53][54]. These models describe the formation and evolution of the first and smallest dark matter halos, and they were specifically built for the purpose of understanding the microhalo populations that arise from EMD scenarios and other cosmologies that feature boosted small-scale inhomogeneity. The application of these models is discussed in Secs. III and IV, and PYTHON codes that implement them are publicly available. 18 Compared to Ref. [11], we also explore more extensively the space of cosmological and observational parameters. Consequently, we are able to clarify EMD detection prospects more precisely and to do so in terms of both T RH and k cut /k RH .
We close with more general remarks about the capacity for PTAs to probe the Universe's early history. The results of a simplified analysis in Ref. [11] suggest that pulsar timing is sensitive to halo masses above approximately 2 × 10 −8 M (see Fig. 2 of that article). Through the size of the cosmological horizon, this mass scale corresponds to a time when the temperature of the Universe was roughly 150 MeV. Consequently, one can make a general statement that PTAs are able to probe the early Universe up to a temperature of about 150 MeV. EMD is not the only early-Universe scenario that results in boosted small-scale power. For instance, pulsar timing could also probe the possibility of domination by a cannibal species with number-changing interactions, which leaves a different spectral imprint [94,95], as long as it occurred below 150 MeV. If the dark matter is an axion and its Peccei-Quinn phase transition occurred below 150 MeV, then the resulting axion miniclusters [96][97][98][99][100][101][102][103] could also be detected using pulsar timing. We anticipate that the procedure through which we characterized the microhalo population arising from EMD-and hence the resulting pulsar timing signals-could be straightforwardly applied to these other scenarios as well. 18 https://github.com/delos/microhalo-models.
In this appendix we evaluate the integral in Eq. (3) that expresses the average relative velocityV rel (V , R) between a microhalo at position R and velocity V and the stars it encounters. Since we assume a spherically symmetric potential, a microhalo's orbit is planar with a static angle φ to the Galactic disk plane. The microhalo's position and velocity within this plane are described by the radius R, angle θ, and associated velocity components V r and V θ , all of which are functions of time. We can decompose the halo's velocity instead into components V = cos φ cos θ 1 + tan 2 θ cos 2 φ −1/2 V θ (A1) parallel to the mean stellar motion, which is assumed to be circular within the Galactic plane, and V ⊥ = V 2 r + sin 2 φ 1 + tan 2 θ cos 2 φ V 2 Here, F (a) ≡ d 3 x (2π) 3/2 (x − a) 2 + y 2 + z 2 e −|x| 2 /2 (A4) with x ≡ (x, y, z), which can be rapidly evaluated using an interpolation table. We also note that F (a) a + 1/a to within one part in 10 4 when a > 3.
with α = 0.78, as the universal density profile of a halo initially possessing an NFW profile after arbitrarily many encounters. However, this profile becomes arbitrarily steep at ) as it responds to an impulsive stellar encounter at time t = 0. Colored curves show the profile in time intervals of 3t dyn , with t dyn defined as in Ref. [54]. The solid black line shows Eq. (B1), the expression proposed in Ref. [54] as a universal fit to the post-encounter density profile. This expression is fitted (with ρs and rs allowed to vary) to the late-time simulated density profile up to r = 4r NFW s . The dashed line shows the modified profile, Eq. (B2), with the same ρs and rs (not an independent fit); it is identical to Eq. (B1) at small radii but approaches ρ ∝ r −4 at large radii, as it should [104]. Moreover, tuned such that β = 5, it also accurately matches the temporally converged part of the simulated density profile at all radii. large radii, whereas an analytic argument by Ref. [104] suggests that the post-encounter profile does not steepen beyond ρ ∝ r −4 . The simulations in Ref. [54] each covered a time duration equal to only about 10 times the dynamical timescale within the halo's scale radius, and within that duration the large-radius density profile did not stabilize sufficiently to test this argument. Consequently, in this appendix we carry out a new impulsive encounter simulation using the methodology of Ref. [54]. The simulated microhalo initially has an NFW density profile with scale parameters r s,0 and ρ s,0 , and it undergoes an encounter with a star of mass M * at impact parameter b and relative velocity V such that M * /(V b 2 ) = 0.52 ρ s,0 /G and b = 16r s,0 . This description completely characterizes an encounter in the limit that the encounter is impulsive. We follow the halo's response for almost 60 dynamical time intervals. Figure 8 shows how the microhalo's density profile evolves in this simulation. To fit the time-stabilized part of this density profile at late times, we propose the form β , α = 0.78, and β = 5. 19 Here, 2 F 1 is the hypergeometric function. This profile is identical to Eq. (B1) at small radii but transitions to ρ ∝ r −4 at large radii. The suddenness of the transition is controlled by the parameter β, and we fix it through comparison with our simulation. As shown in Fig. 8, Eq. (B1) (solid black curve) fails to match the post-encounter density profile at large radii because the latter does not steepen beyond ρ ∝ r −4 at large radii. Instead, Eq. (B2) (dashed curve) accurately describes the post-encounter density profile.