Dark Matter Microhalos From Simplified Models
Nikita Blinov, Matthew J. Dolan, Patrick Draper, Jessie Shelton
FFERMILAB-PUB-21-015-AE-T
Dark Matter Microhalos From Simplified Models
Nikita Blinov,
1, 2
Matthew J. Dolan, Patrick Draper, and Jessie Shelton Fermi National Accelerator Laboratory, Batavia, IL 60510, USA Kavli Institute for Cosmological Physics,University of Chicago, Chicago, IL 60637, USA ARC Centre of Excellence for Dark Matter Particle Physics,School of Physics, University of Melbourne, 3010, Australia Illinois Center for Advanced Studies of the Universe,Department of Physics, University of Illinois, Urbana, IL 61801, USA (Dated: February 11, 2021)
Abstract
We introduce simplified models for enhancements in the matter power spectrum at small scalesand study their implications for dark matter substructure and gravitational observables. Thesemodels capture the salient aspects of a variety of early universe scenarios that predict enhancedsmall-scale structure, such as axion-like particle dark matter, light vector dark matter, and epochsof early matter domination. We use a model-independent, semi-analytic treatment to map bumpsin the matter power spectrum to early-forming sub-solar mass dark matter halos and estimatetheir evolution, disruption, and contribution to substructure of clusters and galaxies at late times.We discuss the sensitivity of gravitational observables, including pulsar timing arrays and causticmicrolensing, to both the presence of bumps in the power spectrum and variations in their basicproperties. a r X i v : . [ a s t r o - ph . C O ] F e b ONTENTS
I. Introduction 2II. From Power Spectrum to Microhalos 5A. Linear Evolution 5B. Collapse and the Formation of Microhalos 10C. Scaling Behaviors 12D. Results for Power Law Peaks 14E. Microhalo Substructure, Evolution, and Disruption 19F. Subhalo Distribution 21III. Probes of Small Scale Structure 23A. Caustic Microlensing 23B. Pulsar Timing 25C. Encounter Rate 28IV. Examples 29A. EMD-induced Bumps in the MPS 31V. Discussion 36Acknowledgments 38A. Window Function Dependence 38References 39
I. INTRODUCTION
Many different microphysical models can produce features in the matter power spectrum(MPS) in the early universe. These features correspond to enhancement or suppression ofthe MPS on given scales relative to the predictions of the minimal Cold Dark Matter (CDM)scenario. Sufficiently pronounced deviations on large enough scales can lead to observablesignatures in the late-time matter distribution.In this paper, we study simplified models for modifications of the MPS on scales well2elow the baryon Jeans length. The resulting dark matter (DM) substructure (which werefer to interchangeably as microhalos, minihalos or simply “clumps”) form at high redshift;as a result, they are significantly less massive than small-scale structure often considered inthe literature (e.g., dwarf and satellite galaxies up to 10 M (cid:12) in mass). Instead, our focusis on very small-scale, non-compact structures orders of magnitude lighter than our Sun,which do not host baryonic objects such as stars.We consider a range of localized features in the small-scale MPS and study the resultingobservational signatures in the late-time universe; we focus on signatures that depend only onthe gravitational properties of the DM distribution. Our aim is to develop a simple procedureto approximately predict the properties of small-scale DM substructure from the broadfeatures of an underlying MPS, and to illustrate the sensitivity of gravitational observationtechniques in a model-independent way.Microphysical models that yield bump-like features can be classified according to whetherfluctuations in the density contrast begin with small initial amplitudes, and features arecreated or magnified over time in the process of cosmic expansion, or whether fluctuationsbegin with large initial amplitudes. Example models in the first class include standard darkmatter models with a modified expansion history such that primordial matter perturbationsexperience a period of enhanced growth. Studies of models of this type have addressed theimpact on small scale structure for WIMP DM candidates [1–3], ALPs/axions [4–6], hiddensector DM [7–11], and others. Models in the second class do not rely on large alterations tothe expansion history. Examples include vector DM from inflationary fluctuations [12–15],ALPs/axions with post-inflationary U (1) P Q breaking [16–19], ALPs/axions with parametricresonance enhancement [20], large semi-diffuse solitons, and others.Given this profusion of models that lead to similar phenomenology in small-scale structure,we argue for a simplified model-style approach. This approach allows us to abstract awaythe microphysics and detailed dependence on any specific dark matter model, and studyinstead a series of simple, few-parameter modifications to the matter power spectrum whichefficiently and accurately capture the dominant physical effects of a broad range of models.One can contemplate many different modifications to the MPS. From the microphysicalpoint of view, one of the simplest features to generate is a bump . A bump is characterized bya rise in the power spectrum over a range of wavenumbers k , followed by a peak (which maybe narrow or broad), and finally a short-distance cutoff at higher k . The properties of boththe rising (smaller k ) and falling (larger k ) sides of the bump are ultimately determined by3icrophysics. For example, in models where a bump is generated by a period of modifiedcosmic expansion prior to Big Bang Nucleosynthesis (BBN), the rising slope is determined bythe equation of state of the universe during the epoch of modified expansion. Both a periodof early matter domination (EMD) and kination lead to linear growth in density perturba-tions [1–3, 6, 7, 10, 11, 21–23]. However, EMD leads to a rise in the (dimensionless) powerspectrum that goes like k , while for kination the power spectrum grows like k . The largestscale that inherits this enhanced growth—i.e., the location of the bump—is determined byother aspects of the microphysics, for example, the scale of reheating. On the small-scaleside, the short-distance cutoff can be associated with a free-streaming scale, a Jeans length,a kinetic decoupling scale, cannibal oscillations, and so on.Observational signatures can be divided according to whether they are sensitive primarilyto the gravitational properties of DM or its other, model-dependent interactions. Examplesof non-gravitational search techniques sensitive to very small-scale DM substructure includedirect and indirect detection experiments. In particular, a number of works have explored thepossibility that enhanced DM substructure might boost the sensitivity of indirect searches forDM annihilations [10, 23]. Examples of purely gravitational signatures include stellar [24, 25]and caustic [26] microlensing, and pulsar timing measurements [27–33]. The latter twotechniques are especially sensitive to the diffuse clumps produced by MPS enhancementsconsidered in this work. These can be used to search for DM substructure in a way that islargely independent of the underlying particle properties.This work is organized as follows. In Sec. II we describe the semi-analytic methods weuse to predict clump formation times, masses, sizes, and densities, and their distributions,taking as input the MPS. We study simple power-law models for MPS enhancements, captur-ing the gross features of many of the particle physics models described above, and considerboth isocurvature and adiabatic primordial perturbations, which are associated with differ-ent transfer functions. We use the Press-Schechter ansatz to determine clump propertiesas a function of collapse redshift. In Sec. III we discuss the observability of microhalosproduced from MPS bumps via gravitational probes or in Earth-microhalo encounters (ifa non-gravitational coupling exists); there we briefly describe pulsar timing (recently stud-ied in detail in Ref. [33] in the context of MPS enhancements similar to the ones consideredhere), and mainly focus on estimating caustic microlensing sensitivities and the clump-Earthencounter rates. In Sec. IV we sketch the analysis for three example models motivated by dif-ferent microphysical scenarios and compare with the caustic microlensing sensitivity curves.4n this section we also use the extended Press-Schechter method to estimate the subhalodistribution of the Milky Way in each of the example scenarios. In Sec. V we summarizeand conclude, discussing in particular the limitations of our simplified analysis. A centraltheme here is that we must necessarily make predictions for the DM substructure on verynon-linear scales. It is crucial to check these results using N -body simulations. Given theprevalence of small-scale structure enhancements in microphysical models of DM, togetherwith the prospect of new detection capabilities, such studies are of very high interest. II. FROM POWER SPECTRUM TO MICROHALOS
The late-time distribution of DM substructure can be estimated by evolving the primordialdensity field to collapse. In this paper we employ a combination of semi-analytic techniquesthat enable rapid exploration of clump parameter space. Given an initial power spectrumof DM density fluctuations, we evolve it to the non-linear regime using linear perturbationtheory. This allows us to estimate average properties of the collapsed objects, in particularsize and density, and to use the Press-Schechter formalism to determine the halo distribution.Halo size and density are the key properties to which the observational techniques discussedin Sec. III are sensitive.The semi-analytic prescription outlined above glosses over the non-linear evolution of DMclumps in the presence of baryons and their potential disruption by gravitational interac-tions. The small-scale structures considered in this work are, at best, at the edge of theresolution of state-of-the art cosmological N -body simulations. We content ourselves withsimple estimates of clump disruption in our galaxy, but emphasize that continuing numeri-cal studies are necessary for making more realistic projections for some astrophysical probes.(Other probes, as discussed in Sec. III, are primarily sensitive to clumps moving freely ingalaxy clusters, and are less subject to disruptive processes.) Below we discuss the predictionof late-time DM clump distribution in more detail. A. Linear Evolution
If the primordial density fluctuations of DM are smaller than unity, linear perturbationtheory can be used to evolve them forward in time. The evolution equations that follow from5ovariant stress-energy conservation are [34]˙ δ = − θ + 3 ˙ φ (1a)˙ θ = − ˙ aa θ + k ψ , (1b)where dots indicate derivatives with respect to conformal time, δ is the Fourier-space den-sity contrast with comoving wavenumber k , θ = i(cid:126)k · (cid:126)v is the velocity divergence, a is theFriedmann-Robertson-Walker (FRW) scale factor, and φ and ψ are the metric perturbationsin conformal Newtonian gauge.Different microphysical realizations of DM usually have different evolution equations atsmall scales (i.e., additional terms in the system above that become important for sufficientlylarge k ), leading to non-trivial k dependence for k larger than some cutoff k c . For example,in the standard WIMP case, small scales are affected by the DM coupling to SM radiation,leading to a non-vanishing pressure. This effect leads to washout of structure on scalessmaller than the corresponding Jeans scale. In models of wave-like DM (such as axion-likeparticles or vector DM), k c corresponds to the comoving Compton wavelength of the fieldbelow which the field cannot be localized, also leading to an effective pressure term. In bothexamples, power on scales k > k c is suppressed.Away from these small scales, the evolution of a DM candidate is fully specified by itsequation of state (assumed to be constant during the epochs of interest), and thereforeEq. (1) is universal. The small-scale effects can often be modeled using the generalized DMframework, which modifies Eq. (1) by additional model-dependent pressure and sound speedterms [35]. The modes of interest for structure formation first evolve through a period ofradiation domination (RD, possibly after a period of non-standard cosmology), followed bymatter domination (MD). We will consider evolution during these two phases separately.The formal solution of Eq. (1) during radiation domination is [36] δ = C ln a + C + (cid:90) η [ln a (cid:48) − ln a ] a (cid:48) ˙ a (cid:48) ( k ψ − φ − φ ) dη (cid:48) (2)where η is conformal time. In Newtonian gauge, perturbations are constant prior to hori-zon entry, so C = 0. The value of C depends on whether the modes are adiabatic orisocurvature.For adiabatic perturbations, superhorizon initial conditions give C = 3 φ/
2, so the secondand third terms are comparable; when the mode enters the horizon during radiation domi-nation, the rapidly-decaying gravitational potentials in the third term generate a logarith-6ically growing contribution. Well after horizon entry, the solution in standard cosmologyis well-described by (see Appendix B2 of Ref. [36]) δ = C + I φ ln (cid:18) I aa hor (cid:19) , (3)where I ≈ . I ≈ . a hor /a eq ≈ √ k eq / (2 k ), and φ is the initial (superhorizon) valueof the gravitational potential. This parametrization is valid for a hor (cid:28) a (cid:28) a eq , where a hor is the scale factor at horizon entry of mode k and a eq is the scale factor at matter-radiationequality (MRE).In models with isocurvature initial conditions, C is not correlated with the size of gravi-tational potentials. This means that if C (cid:29) φ , the perturbations do not grow significantlyduring RD. Note that Eq. (3) is correct in both adiabatic and isocurvature cases.The part of the solution in Eq. (2) that is sourced by the gravitational potentials is onlyvalid for modes that enter during standard RD. However, the form of the late-time solution,Eq. (3), holds even for modes that entered during a period of non-standard cosmology, albeitwith different values of C , I and I . These parameters are usually k -dependent and canbe determined by matching the evolution of δ through a period of modified cosmology (andinto RD) onto the ln a and constant terms. For detailed examples of this matching, seeRefs. [1, 6, 22]. For example, if a mode has entered the horizon during a period of earlymatter domination, then C = I φ = 2 φ ( k/k ∗ ) / I = ( k ∗ /k ) [1, 6]. If a mode entersduring a period of kination, then C = I φ = 3 . φ (cid:112) k/k ∗ , I = (cid:113) k ∗ / ( √ k ) [22]. In theseexpressions k ∗ is the comoving horizon size at the end of the period of modified cosmology.Next we consider the evolution of perturbations through matter-radiation equality (MRE)and beyond. Well after the radiation contributions to gravitational potentials have decayed,the system in Eq. (1) can be cast as the Meszaros equation [36, 37], δ (cid:48)(cid:48) + 2 + 3 y y (1 + y ) δ (cid:48) − y (1 + y ) δ = 0 , (4)where y ≡ a/a eq , primes denote derivatives with respect to y , and we set the baryon densityto Ω b = 0 for simplicity. Matching the solution of this equation in the y (cid:28) δ = 32 (cid:20) C + I φ ln (cid:18) I e − a eq a hor (cid:19)(cid:21) U ( y ) − I φU ( y ) , (5) While we mostly follow the notation of Ref. [36] here, note that our parameter I in Eq. (3) slightly differsfrom that of Ref. [36], which incorporates C into I . It is somewhat convenient to keep them separatewhen discussing isocurvature and adiabatic perturbations in a unified way; when we consider the specificcase of adiabatic perturbations alone, we will absorb C into I , in which case I → . U ( y ) = 23 + y, U ( y ) = 458 (cid:20)(cid:18)
23 + y (cid:19) ln √ y + 1 √ y − − (cid:112) y (cid:21) (6)are the growing and decaying solutions of Eq. (4). We again note that the constants C , I ,and I depend on the expansion rate of the universe when the mode entered the horizon,and can be determined by matching as described above.In what follows we will need the autocorrelation function (cid:104) δ (cid:105) . We write it as (cid:104) δ (cid:105) = (cid:18) (cid:19) (cid:2) D iso (cid:104) C (cid:105) + D adi I L (cid:104) φ (cid:105) + 2 D iso D adi I L (cid:104) C φ (cid:105) (cid:3) (7)where L = ln (cid:18) I e − a eq a hor (cid:19) (8)is the logarithm from Eq. (5), and we defined the growth functions: D iso = U , D adi = U − L U . (9)At late times, D iso ≈ D adi ≈ a/a eq ( a/a eq (cid:29) . (10)Let us consider two limits of Eq. (7). First we take the standard case of adiabatic initialconditions and C = 3 φ/
2. In this case, for a/a eq (cid:29)
1, it is easy to see that (cid:104) δ (cid:105) = D adi I L (cid:48) P R , (11)where L (cid:48) = ln (cid:18) I e − a eq a hor (cid:19) , I → .
594 (12)and we related (cid:104) φ (cid:105) to the curvature power spectrum P R ( k ) = (cid:18) (cid:19) (cid:104) φ (cid:105) = 2 π k × A s (cid:18) kk (cid:19) n s − . (13)The 2018 Planck best-fit values of the spectral index and initial amplitude of curvatureperturbations are n s = 0 . ± .
004 and ln(10 A s ) = 3 . ± .
014 at the pivot scale k = 0 .
05 Mpc − [38]. The Meszaros equation (4) can be solved for a non-zero baryon fraction f b as in Ref. [36]; we present the f b = 0 solutions for simplicity but use f b (cid:54) = 0 in our numerical results. I changes because in this particular case it is convenient to absorb C into I . See Footnote 1.
8n the large isocurvature limit C (cid:29) φ , we instead find (cid:104) δ (cid:105) = D iso (cid:18) (cid:19) (cid:104) C (cid:105) . (14)Comparing Eqs. (11) and (14) we see that large isocurvature perturbations do not benefitfrom the radiation driving effect or the logarithmic growth during RD. This difference canbe significant – for example, adiabatic modes that enter the horizon at T ∼ MeV have I L ≈ C and I in concrete models) until they are of O (1), at which point perturbation theory breaks down, signaling the onset of gravitationalcollapse. This, however, does not mean that microhalos begin to self-gravitate at this point.If collapse occurs during RD, the DM provides a negligible contribution to the gravitationalpotential, unless that overdensity is so large that the region around it is locally matter-dominated. This can occur in two distinct ways, corresponding to isocurvature and adiabaticperturbations. In the isocurvature case, the initial density fluctuation can be so large thatlocal MRE is attained well before the global MRE. This was studied in Ref. [16] in a sphericalcollapse model. Adiabatic perturbations, on the other hand, by definition must reach collapsedynamically. Ref. [10] studied the case where this occurs during the logarithmic growthof a perturbation (originally enhanced through a period of EMD). In this example, theperturbation becomes non-linear well before it begins to self-gravitate and virialize (whenlocal MRE is attained).In what follows, we focus on fluctuations that collapse only after global MRE, enabling theuse of standard results on spherical collapse and avoiding the difficulties of collapse withoutvirialization. While not fully general, this restriction still captures many of the scenariosintroduced in Sec. I. Moreover, we expect that the gravitational probes discussed later aremore sensitive to objects that collapse at or after MRE, given the range of mass scales thatthey can probe. Other probes, such as indirect detection, may be significantly enhanced in objects that form earlier [10]. . Collapse and the Formation of Microhalos The Press-Schechter (PS) formalism enables a semi-analytic understanding of structureformation from the evolution of the linear density contrast discussed above [39]. PS pos-tulates that the fraction of matter in collapsed objects of size R at a given time is relatedto the cumulative probability distribution to find δ R > δ c , where δ R is the density contrastsmoothed over scales R and δ c is the critical linear density contrast, i.e., the fractional over-density that would collapse in a non-linear treatment. For spherical collapse during MD, δ c ≈ .
686 [40]. The second key assumption in PS is that the probability distribution of δ R is Gaussian and therefore fully specified by its variance, σ ( z, R ) = (cid:90) d k (2 π ) (cid:104) δ ( a, k ) (cid:105) W ( kR ) , (15)where W is a smoothing kernel and we use redshift z = 1 /a − a interchangeably as independent variables. The comoving smoothing scale R can be translatedto a mass scale M ∝ ¯ ρ R ( ¯ ρ is the average DM density today), with the precise relationshipbetween M and R depending on the choice of W (see Appendix A).The variance σ encodes the evolution of structure as a function of time. A useful quantitythat can be derived from it is the critical mass M ∗ ( z ) such that σ ( z, M ∗ ) = δ c . (16)At a given redshift, M ∗ gives the typical mass of the largest structures that collapse at thistime. In all microphysical models, M ∗ is a monotonically increasing function of 1 / (1 + z ),indicating the hierarchical assembly of smaller clumps into larger and larger structures.The density fluctuation variance σ and the characteristic mass M ∗ can be used to modelthe distribution of microhalos as a function of redshift within the PS framework and itsextensions. We will do this in Sec. II F, but for now we briefly describe the basic propertiesof the earliest-forming microhalos.In the spherical collapse model, the average density of matter in a newly-formed objectdepends only on the collapse redshift z c [40–42]: ρ ( z c ) ≈
178 ¯ ρ ( z c ) ≈ / cm (cid:18) z c (cid:19) , (17) Here “typical” means 1 σ fluctuations at redshift z . Objects of a fixed mass M ∗ ( z ) can form earlier than z , but they must arise from rarer upward fluctuations in the density contrast. ρ ( z c ) is the background density at the time of collapse. Then a characteristic radiusof clumps forming at redshift z c can be assigned as R ∗ ( z c ) = (cid:18) M ∗ ( z c )4 πρ ( z c ) (cid:19) / . (18)We see that objects that form earlier are denser and more compact. In Sec. II E we willargue that earlier-forming microhalos are more likely to survive tidal disruption by stars,other clumps, and the galactic halo.After a clump forms, its early-time evolution depends on whether the bump in the powerspectrum is a localized spike or a broad enhancement [43, 44]. For spiked features, a narrowrange of scales goes nonlinear at approximately the same time, while all other scales mustevolve for a much longer period before crossing the collapse threshold. In this case, theclumps that form first evolve for a long time in relative isolation, and they have featurelessprofiles with peaked inner densities scaling as ρ ∼ r − / . In the case of broad enhancements,larger and larger clumps are continually forming, coalescing out of their predecessors. Thusthe enhanced population of microhalos develops with a substructure of smaller, older clumps,which is expected to survive at least an order of magnitude in redshift after formation [43,44]. The average density profiles in this case are modified by accretion, trending towardthe Navarro-Frenk-White (NFW) [45] profile with ρ ∼ r − . In the examples of Sec. II Dand IV we will consider both narrow and broad enhancements to the MPS. For simplicity,however, in both cases we will model the density profile as NFW. For a narrow spike, thischoice underestimates how compact the resulting objects are, and therefore overestimates thedisruption probability. This choice is thus conservative in estimating the late time abundanceof these objects.The NFW profile is given by ρ ( r ) = 4 ρ s ( r/r s )(1 + r/r s ) , (19)where ρ s and r s are scale mass and scale radius. The mass within a given radius r is then M ( r ) = 16 πρ s r s f ( r/r s ) , (20)where f ( c ) ≡ ln( c + 1) − cc + 1 . (21)The scale mass is the mass within the scale radius: M s ≡ M ( r s ) = 16 πρ s r s f (1) . (22)11e can approximately relate the collapse mass M ∗ and radius R ∗ , given by Eqs. (16) and (18),to r s , M s , and ρ s . Recall that R ∗ is defined as the radius within which the density is 178 ¯ ρ ( z c ),where ¯ ρ ( z c ) is the background density at collapse. This means that M ∗ and R ∗ are closeto the virial quantities M and R , defined for halos with an average density 200 timesthe background density, evaluated at collapse. Simulations of ΛCDM halos suggest that theconcentration parameter c = R /r s shortly after formation is about 2 for Earth-masshalos in ΛCDM [46] (one needs to take the Earth mass c ≈
60 and redshift it back to z ≈
30 where these halos were measured in the first place), implying that M ∗ and R ∗ aresimilar to the NFW scale quantities. Let c ∗ be the concentration parameter at formation,i.e., c ∗ = R ∗ /r s . (23)Let us also assume that R ( z c ) = R ∗ , neglecting the difference between 200 and 178. Thisallows us to find M s and ρ s from c ∗ , R ∗ and M ∗ : M s = f (1) f ( c ∗ ) M ∗ . (24)The scale density is then obtained from Eq. (22). As mentioned above, halos of differentmasses can form at around the same redshift. Thus a more realistic description of themicrohalo population would account for scatter in the concentration parameter c ∗ [47]. Thisscatter reflects both the “rarity” of a collapsing density peak (or, equivalently, the spread information redshifts of halos with similar masses) [48] and the subsequent assembly history.We will focus on a single representative value of c ∗ for simplicity. C. Scaling Behaviors
Here we derive some simple scaling laws that allow estimation of the properties of thelargest clumps forming from an enhancement in the power spectrum. We modify the pri-mordial power spectrum in Eq. (13) by introducing a “bump” function B that enhances thepower spectrum over some range of scales, such that P R ( k ) → P R ( k ) B ( k ). The resultingvariance in the density contrast on a given comoving scale ˆ k at scale factor a is: σ ( a, R = 1 / ˆ k ) = 12 π (cid:90) dk k W ( k/ ˆ k ) | D ( a, k ) | P R ( k ) B ( k ) ∼ A s D ( a, ˆ k ) (cid:32) ˆ kk (cid:33) n s − (cid:32) ˆ k ¯ k (cid:33) n . (25)12
10 50 100 500 100010 - - Figure 1. Comparison of the critical mass M ∗ ( z ) derived by numerically solving Eq. (16) (solid line)with the analytic estimate in Eq. (29) (dashed line). This example corresponds to an adiabaticpower-law enhancement of the MPS defined in Eq. (30) with n = n = 3, k p /k eq = 10 and E = 10 . In the shaded grey regions the estimate of Eq. (29) is not valid, since the power spectrumchanges its logarithmic slope; the low- z boundary is found by setting the factor in parenthesis inEq. (29) to 1, while the high- z boundary corresponds to collapse of the peak scale k p . Here W is a window function that isolates scales of order ˆ k in some way. In the second line weassume that B ∼ ( k/ ¯ k ) n characterizes the nontrivial behavior of B at the largest enhancedlength scales, where ¯ k is the smallest value of k with any enhancement, and we focus onthose scales since they will determine the largest substructures. The relevant scale-dependentgrowth function D ( a, k ) depends on whether the bump is isocurvature or adiabatic: D ( a, ˆ k ) ∼ (cid:18) aa eq (cid:19) L α , L ∼
10 ln( a eq /a hor ) ∼
10 ln(ˆ k/k eq ) , (26)where the scale factor and k dependence follows from Eqs. (11) and (14), and the exponent α is 1 (0) for adiabatic (isocurvature) perturbations. Then the collapse redshift z c (ˆ k ) is determined by setting the critical linear density contrast The factor (cid:104) δ (cid:105) ∝ | D ( a, k ) | P R is sometimes written as | T ( k ) D ( a ) | P ΛCDM , where T ( k ) is a transferfunction, D ( a ) = a , and P ΛCDM ∼ ( k/k eq ) P R (see, e.g., the textbooks [40, 49]). The conventionalnormalization of P ΛCDM is chosen such that T ( k ) ≈ L α ( k eq /k ) ) for modes that enter the horizonafter (or before) matter-radiation equality. These conventions and normalizations are convenient forstudying wavenumbers in the vicinity of k eq , but are somewhat encumbering in our case of interest with k (cid:29) k eq , where the factors of k/k eq always cancel between T ( k ) and P ΛCDM . δ c ∼ (cid:112) A s (cid:32) ˆ k ¯ k (cid:33) n / (cid:32) ˆ kk eq (cid:33) ( n s − / (cid:32) z eq z c (ˆ k ) (cid:33) L α . (27)Turning the perspective around, we can find the mode k ∗ ( z ) that characterizes scales col-lapsing at z . If we associate a mass with k ∗ as M ∗ ∼ ¯ ρ ( a )( k ∗ /a ) − ∼ ¯ ρ /k ∗ and set n s ≈ δ c ∼ (cid:112) A s (cid:18) MM ∗ (cid:19) n / (cid:16) z eq z (cid:17) L α (28)where M ∼ ¯ ρ / ¯ k and M ∗ ( z ) ∼ M (cid:18) z eq L α √ A s zδ c (cid:19) /n . (29)In Fig. 1 we compare the analytic estimate in Eq. (29) (dashed line) with the critical mass M ∗ derived by solving Eq. (16) numerically (solid line) for a particular choice of B definedbelow in Eq. (30). We see that Eq. (29) provides an excellent approximation to the criticalmass over nearly two orders of magnitude in z . In the shaded grey regions the estimateof Eq. (29) is not valid, since the power spectrum changes its logarithmic slope; the low- z boundary is found by setting the factor in parenthesis in Eq. (29) to 1, while the high- z boundary corresponds to collapse of the peak scale k p where the power law growth with n ends (see Eq. (30)). D. Results for Power Law Peaks
We use simplified models to parametrize bumps in the power spectrum. We add power-law enhancements over the primordial ΛCDM baseline in Eq. (13), specifying a wavenumber k p which determines the location of the peak of the bump in the power spectrum and anumerical factor E which gives the size of the enhancement over the value for ΛCDM at thesame value of k p . We assume that the rise and fall in the power spectrum on either side ofpeak are power laws in k , which we allow to have different exponents n and n . The finalform of the peak enhancement can be written as B ( k ) = E × ( k/k p ) n k < k p ( k p /k ) n k ≥ k p , (30)14uch that the full power spectrum is given by P R ( k ) for k < k p / E /n and B ( k ) P R ( k )otherwise. This simple expression captures the main features of a broad range of well-motivated models. For instance, an early epoch of kination predicts n = 1 [22], while earlymatter domination gives n = 4 [1]. The inflationary production of massive scalar spectatorfields can generate 0 ≤ n ≤ n = 3 [12].Meanwhile, altering the falling slope on the right hand side of the peak approximates anumber of physical cutoff mechanisms. A very rapidly falling peak, which we can modelas e.g. n = ∞ , can correspond to a Jeans cutoff, related to the wave nature of a DMcandidate [6] or its coupling to the radiation bath [1]. The other extreme n = 0, where therise is followed by a plateau, can be found in models where a modified expansion period isrelatively short and preceded by radiation domination [51]. (In such cases the plateau shouldreflect logarithmic growth; this is a minor effect that is not captured by the simplified model.)An intermediate regime, where the fall is less steep than the rise, can be found in models ofdark photon DM produced by inflationary fluctuations, where n = 1 [12] or EMD modelswith cannibal interactions, where n = 2 [11]. We discuss these three scenarios in furtherdetail in Sec. IV.This simplified model lets us study the statistical properties of dark matter clumps as afunction of the peak location, peak height, and the power laws on the rising and falling sidesof the peak, for both adiabatic and isocurvature perturbations. These are the most impor-tant data affecting final distributions and observables. We take Eq. (30) to parametrize an“effective” primordial power spectrum, and then use the standard growth functions derivedin Sec. II A to evolve them forward in time. The effect of any non-standard cosmologicalevolution is absorbed into B ( k ), even though it can occur well after inflation. More explicitly,in Eq. (15) we use (cid:104) δ ( a, k ) (cid:105) = D adi ( a ) I L (cid:48) ( k ) B ( k ) P R ( k ) (adiabatic) D iso ( a ) [ B ( k ) − P R ( k ) + D adi ( a ) I L (cid:48) ( k ) P R ( k ) (isocurvature) , (31)where the curvature power spectrum P R is given in Eq. (13) and the various growth factorsare for the standard cosmology. For reference, growth functions D are defined in Eq. (9), I is defined below Eq. (3), and L (cid:48) is given in Eq. (12). Note that when B ( k ) represents anisocurvature enhancement, we assume that the MPS is still the standard adiabatic one atlarge scales; this is enforced in the second line by the presence of B ( k ) −
1, which vanishesfor k < k p / E /n . We now consider varying each of the MPS parameters in Eq. (30) in turn,15 - - - - - - - - - - - Figure 2. Power spectra and the resulting microhalo properties for different peak amplitudes (toprow) and positions (bottom row). Upper left: The dimensionless power spectrum for differentpeak heights for E = 10 i/ with i ∈ {− , − , , , } (corresponding to the gradation of lighterto darker lines) in Eq. 30; peak position k p /k eq = 10 and slopes n = n = 3 are held fixed.Upper right: Minihalo scale density as a function of scale mass, ρ s ( M s ), for these peak heightvariations, showing both adiabatic (solid) and isocurvature (dashed) perturbations. Dashed greylines indicate the collapse epoch z c . Lower left: The power spectrum for different peak positions k p /k eq = 10 i/ with i ∈ {− , − , , , } (corresponding to the gradation of lighter to darkerlines); the enhancement is held fixed at E = 10 . Lower right: As for upper right, showing varyingpeak position. Note the scales of the axes in the right-hand plots differ. measuring the sensitivity of the properties of the resulting microhalos.16n Fig. 2 we vary the height of the peak in the top line as E = 10 i/ for i ∈{− , − , , , } . The peak location is kept fixed at k p /k eq = 10 and the slopes to n = n = 3. These enhancements to the power spectrum are shown in the upper-leftpanel. The upper right panel shows the properties of the collapsed clumps as in the M s - ρ s plane. The scale mass M s and scale density ρ s can be related to the critical mass M ∗ andradius R ∗ assuming a value for the concentration parameter at formation, as described inSec. II B.Each point on the lines in the right-hand panels of Fig. 2 corresponds to a specific value ofthe collapse epoch z c for the clumps of that mass, since z c relates the scale mass to the scaleradius (and hence scale density). The horizontal grey-dashed lines correspond to collapseepochs of z c = 10 , ,
250 and 1000 from bottom to top, respectively. As we discuss inSec. II E clumps which form at redshifts prior to z = 250 are unlikely to be severely disruptedin our galaxy.We present results for both adiabatic (solid) and isocurvature (dashed) perturbations.Increasing the enhancement causes modes of a given k to go nonlinear earlier. Fluctuations onphysical scales to which those modes contribute significantly thus collapse earlier, increasingthe density of the resulting clump. For adiabatic perturbations, enhancements greater than E = 10 lead to the formation of clumps which are likely to survive in the environment of theMilky Way to the present day. Isocurvature perturbations require an enhancement factor of E ∼ . for the same result. Since isocurvature fluctuations do not benefit from logarithmicgrowth or radiation driving during RD (see the discussion in Sec. II A) they require a largerenhancement in order to collapse at the same time as an equivalent adiabatic mode.The lower panels of Fig. 2 vary the peak position as k p /k eq = 10 i/ for i ∈ {− , − , , , } ,keeping the enhancement factor fixed to E = 10 . Shifting the peak toward lower k meansthat larger scales become non-linear earlier on, leading to an increase in the density for agiven mass. This also leads to an increase in the maximum clump mass that receives adensity enhancement over ΛCDM.In Fig. 3 we present results for varying the slope on the left-hand side of the peak (atlower values of k , upper plots), and the right-hand side of the peak (at higher values of k , lower plots). We keep the peak location and enhancement fixed at k/k eq = 10 and E = 10 in these plots. In the upper panels we take n ∈ { , , , } and in the lower panels n ∈ { , , , ∞} , with the other exponent fixed to 3.Fattening the low- k tail causes longer wavelength modes to collapse earlier if we hold k p - - - - - - - - - - - - - - Figure 3. Power spectra and the resulting microhalo properties for different rise (top row) and fallslopes (bottom row). Upper left: Rise slope variations for n ∈ { , , , } (corresponding to thegradation of darker to lighter lines) in Eq. 30, with fixed n = 3, k p /k eq = 10 and E = 10 . Upperright: Minihalo scale density as a function of scale mass, ρ s ( M s ) for the different rise slopes. Dashedgrey lines indicate the collapse epoch z c . Lower left: Power spectra for different peak fall slopes n ∈ { , , , ∞} (corresponding to the gradation of darker to lighter lines), with fixed n = 3, k p /k eq = 10 and E = 10 . Lower right: As for upper right, showing varying peak fall slopes. Notethe scales of the axes in the right-hand plots differ. and E fixed. It increases the largest length scales associated with the MPS enhancement,and thus the mass of the largest objects with density enhancement over the ordinary ΛCDMscenario. On the high- k side of the peak, broadening the cutoff slope with k p and E fixed18ushes the collapse epoch z c of the smallest structures to higher redshifts, resulting in denserobjects.We do not show results for isocurvature enhancements in Fig. 3. These would followthe same general trends of Fig. 2: compared to adiabatic enhancements, isocurvature bumpswould result in smaller characteristic densities and masses due to the absence of the radiationdriving effect at early times. E. Microhalo Substructure, Evolution, and Disruption
The survival of early-forming microhalos until late times is paramount to their possibleobservation using the methods described in Sec. III. Gravitational interactions with otherclumps, the galactic halo, and baryonic objects can transfer energy to a clump, changingits internal structure, stripping away some of its matter or disrupting it completely. Theearly-time hierarchical assembly of smaller structures into larger ones tends to flatten thecore density profiles [43, 44]; this effect is more pronounced in broad MPS enhancements andmotivates the NFW profile we use to model these objects, as discussed in Sec. II B. At latertimes, microhalos accreted onto the host galactic halo can be stripped of a large fraction oftheir mass. This mass loss occurs beyond a tidal radius R t from the microhalo center; R t is estimated as the radius beyond which the tidal force of the host halo is larger than themicrohalo’s self-gravity [52]: R t ( R ) ≈ R (cid:20) M ( R t ) M h ( R ) (cid:21) / , (32)where R is the radius of the microhalo orbit around the host, and M h ( R ) is the host massinterior to R . Thus R t /r s ∼ ( ρ s /ρ h ) / , where ρ h is the average density of the host, so R t /r s (cid:29) r (cid:46) r s . As an explicit example we can estimate the tidal radiusfor microhalos orbiting the MW at R ≈ R ⊕ ≈ . M h ( R ) ≈ M (cid:12) [53]. This isrelevant for the pulsar timing probe of microhalos [33]. We find that R t /r s > z ∼
10. In all examples we consider in this work, the MPSenhancements lead to formation of microhalos at much earlier times, and therefore we donot expect tidal disruption by the host halo to significantly alter these objects at late times. This conclusion might not hold for all microhalos, such as those on highly radial orbits that take themclose to the center of the host. z (cid:29)
10. In the impulse approximation,the specific impulse transferred to the clump constituents from an encounter with a star withmass M (cid:63) and impact parameter b is of order∆ v ∼ (cid:18) GM (cid:63) r vir b (cid:19) × (cid:18) bv rel (cid:19) (33)where the first factor is the characteristic tidal force and the second is the encounter timescale.Therefore the energy transfer relative to the binding energy of the clump is of order [54, 55]∆ EE b ∼ v vir ∆ vv vir ∼ GM ∗ r vir b v rel v vir . (34)The ratio r vir /v vir ∝ √ ρ , where ρ is the clump mass density, so the density is the onlymicrohalo property on which the relative energy transfer depends in this approximation.Setting ∆ E/E b = 1 defines a critical impact parameter [54, 55] b c = √ GρM (cid:63) v rel . (35)Encounters with larger b transfer less energy, but are also more probable. The probabilityto transfer a total energy of order E b over many subcritical encounters is comparable tothe probability of a single critical encounter [54, 55], so the total probability of transferring∆ E ∼ E b is roughly double the probability of a single critical encounter.On the other hand, critical encounters with b < b c may not completely destroy theclump, because the transferred energy may not be efficiently redistributed before the outerlayers are stripped away [52]. In this case, a dense remnant core remains. Nonetheless, theprobability of a critical encounter is a reasonable estimate for the survival rate of largerclumps. In particular, if the clumps contain substructure, the smaller, denser constituentsmay be spilled out by disruptive encounters, with their own structure preserved.We define the total disruption probability P to be twice the critical encounter probabil-ity [54] P = 2 nπb c S (36)20here n is the number of galactic disc crossings and S is the orbit-averaged stellar columnmass density along the clump orbit. For clumps in the Milky Way, n ∼
100 and S ∼ M (cid:12) / pc . Since b c depends only on the clump density, which in turn depends only on theredshift of collapse, we can estimate [57] P ≈ n (cid:18) z c (cid:19) / . (37)Therefore, we arrive at a qualitative estimate that clumps forming before z ≈
250 typicallysurvive, while clumps forming later are typically strongly disrupted. A similar calculation ofthe survival probability has been validated in N -body simulations of a microhalo traversinga dense field of stars [55].While the pulsar timing techniques described in Sec. III probe the MW microhalo distri-bution, caustic microlensing is mainly sensitive to the diffuse microhalo population in galaxycluster lenses. These clumps can be disrupted by the diffuse stars in the cluster. We can useEq. (36) to estimate the disruption probability. Taking some representative cluster param-eters v rel ∼ / s, S ∼ M (cid:12) / kpc , and n ≈
10 we find that disruption is significantonly for objects that formed after z c (cid:46) We therefore expect clumps originating fromMPS enhancements to survive in the diffuse cluster environment.We overlay contours of z c in the plots, both to give a sense for the pace of structureformation once the bump in the power spectrum starts to go nonlinear, and to estimatethe sensitivity of the clump population and the observables to the variations in the survivalcutoff on z c . F. Subhalo Distribution
In this subsection we estimate the microhalo distribution inside galaxies. The densityfluctuation variance in Eq. (15) can be used to construct a global (Universe-averaged) distri-bution function within the Press-Schechter formalism [39]. In this framework, the fraction df of matter in objects of mass in the range [ M, M + dM ] is dfdM ( M, z ) = (cid:114) π δ c M σ (cid:12)(cid:12)(cid:12)(cid:12) d ln σd ln M (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:18) − δ c σ (cid:19) . (38) A more sophisticated treatment of the stellar distribution yields a similar disruption probability [56]. See, e.g., Ref. [58] for numerical studies of the diffuse stellar population in clusters and Ref. [59] for relevantparameters of a cluster where caustic microlensing was observed. dndM ( M, z ) = ¯ ρ M dfdM , (39)where ¯ ρ is the DM density today. Generalizations of the PS ansatz exist (such as Ref. [47,60]) that provide a better fit to N -body simulations. However, these simulations have beenmostly performed for standard power spectra. (Exceptions include, e.g., Refs. [43, 44, 61],which study power spectrum enhancements from a modified cosmology. Post-inflationaryaxion DM models are another example widely studied in simulation; for a recent analysis,see [62].) It is therefore important to validate the predictions of the Press-Schechter ansatzand its extensions for other primordial power spectra.Note that M in Eq. (38) is the mass of recently-formed objects, i.e., those overdensitiesthat are just collapsing at z . Evaluating the mass function at z = 0 as a proxy for thedistribution of microhalos in the MW (or other galaxies or clusters relevant for causticmicrolensing) misses two important effects: disruption of clumps, as described in the previoussection, and the fact that the objects forming now are assembled from smaller objects;disruption may “spill” these constituents such that these “sub-microhalos” are the relevantclumps at late times. A coarse model that side-steps these issues is to simply evaluate thedistribution function at an earlier time, e.g., z ∼ conditional probability distribution of microhalos of mass M identified at z to end upin large-scale overdensities of mass M and (later) redshift z is given by dfdM ( M , z | M , z ) = ∆ δ √ π (∆ S ) / (cid:12)(cid:12)(cid:12)(cid:12) dS dM (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:20) − (∆ δ ) S (cid:21) . (40)Here ∆ δ = δ − δ , with δ i = δ c / D ( z i ), and ∆ S = S − S , with S i = σ ( z i = 0 , M i ) and thedensity variance is evaluated using the sharp- k filter. The quantity ( df /dM ) dM gives the Note that since we want to take z i deep in the MD era, the growth functions for adiabatic and isocurvaturefluctuations are nearly equal, D ≡ D iso = D adi as discussed in Sec. II. M in microhalos in the mass range [ M , M + dM ]. Other distributionscan be obtained from Eq. (40) by adding appropriate factors (e.g., multiplying by M /M gives the number distribution of microhalos as a function of their mass at z ).In Sec. IV we will use Eq. (40) to show that microhalos make up an O (1) fraction of theMW DM mass in several representative models, by taking z = 0, M = M MW ≈ M (cid:12) ,and z ≈ z (cid:46)
250 are disrupted.
III. PROBES OF SMALL SCALE STRUCTUREA. Caustic Microlensing
An interesting gravitational observable was considered in Refs. [66, 67] and applied toaxion miniclusters by Dai and Miralda-Escud´e [26], and in Refs. [6, 20] to other kinds ofsubstructure. The goal is to observe a background star near a galaxy cluster lens causticthat undergoes microlensing by a compact object in the cluster, giving rise to a time-varyinglight curve with a total magnification that can reach a factor of 10 − . During the eventduration τ , the star appears to move a distance d ∼ τ v rel µ , where µ is the magnification and v rel is the star-cluster relative velocity. As the line of sight sweeps over a clumpy dark matterdistribution, the magnification varies, leading to jitter in the light curve. These variationsmay be observable if they are the order of the brightness of the background star when it isnot being enhanced by the caustic.The surface density field Σ of a cluster DM halo is related to the line-of-sight integral of thethree-dimensional density distribution. It therefore inherits fluctuations from the intrinsicgranularity of the DM distribution in the presence of a high abundance of microhalos. Onnon-linear scales the density power spectrum can be computed in the halo model [47] fromthe halo mass function df /d ln M and the Fourier transform of the clump density profile˜ ρ h ( q ; M ): P ρ ( q ) = ¯ ρ (cid:90) dMM dfd ln M | ˜ ρ h ( q ; M ) | , (41)where ¯ ρ is a mean density in some region of the cluster. From this the surface density powerspectrum can be estimated by partitioning the cluster lens into a series of slices and assuming23hey have the same mass function, with the result [26]: P Σ ( q ⊥ ) = ¯Σ (cid:90) dMM dfd ln M | ˜ ρ h ( q ⊥ ; M ) | (42)where q ⊥ is the surface Fourier mode. Here ¯Σ = (cid:82) dL ¯ ρ ( L ) is the mean surface density. Weare interested in a long line of sight through galaxy clusters, of order a Mpc. In the casesanalyzed in [26], the mean surface density along the line of sight is dominated by the cluster,¯Σ = ¯Σ cl (as opposed to DM along the line of sight but outside of the cluster).We will define the sensitivity criterion in terms of the amplitude of fluctuations of thelensing convergence κ = Σ / Σ crit , (43)where Σ crit = c / (4 πGD eff ) D eff = D L D LS /D S . (44)Here D L,S,LS are distances to the lens, the source, and from the source to the lens. Σ crit sets the density scale at which an isolated lens produces multiple images. The convergencepower spectrum is then given by P κ ( q ) = P Σ ( q )Σ crit , (45)or, in dimensionless form, ∆ κ = q P κ π . (46)To obtain rough observational sensitivities we require that ∆ κ is larger than some value10 − − on some scale q (such that the observed brightness fluctuations are O (1) if the mag-nification is 10 − ). Ref. [20] employed a monochromatic simplified assumption, where onefocuses only on clumps of mass M s composing a fraction f of the dark matter ( df /d ln M = f M δ ( M − M s ) for some constant 0 < f < κ can be expressed analytically,∆ κ = 1log(2 / √ e ) √ Σ cl f M s qr s g ( qr s )Σ crit r s , (47)24here g is a function appearing in the Fourier transform of the NFW profile. The quantity qr s g ( qr s ) reaches a maximum of 0.35 at qr s = 0 .
77. In evaluating sensitivities we use D eff ∼ Gpc to compute Σ crit and set Σ cl = 0 . crit . These values are consistent with the observationof the magnified star LS1 reported in [59].To estimate the sensitivity of caustic microlensing, we require that max q ∆ κ ( q ) > − and impose three extra criteria [20]: the lensed light should sweep over many clumps as thestar traverses d to enable a statistical treatment; the clumps should be smaller than d tocause significant density fluctuations during the microlensing event; and the clumps shouldbe large enough that the fluctuations they induce in the light curve are not washed out bythe finite size of the lensed source star.The first criterion, that there are many clumps along the line of sight, can be imposed byrequiring f π ( d/ Σ cl /M s >
10 or so. This results in a sharp cutoff of M s ∼ − M (cid:12) . Thesecond criterion, r s < d , can be imposed by picking representative values for τ , v rel , and µ ,leading to d ∼ AU [20]. For the third criterion, we take the minimum sensitivity length ∼
10 AU used in Ref. [26] and require r s >
10 AU.These criteria and the requirement on ∆ κ lead to the sensitivity curves on the M s , ρ s planeshown in Fig. 4. The thickness of the bands comes from varying f from 0 . B. Pulsar Timing
Dark matter clumps can also be probed with searches using pulsar timing arrays (PTAs).This possibility has recently been studied in detail by [31–33]. Two types of signal arepossible. The first is a Doppler-shift in the frequency of the pulsar as a clump passes near25 igure 4. Estimate of the sensitivity of caustic microlensing for a monochromatic microhalo distri-bution in the plane of microhalo scale mass M s and scale density ρ s (assuming NFW profiles forthese clumps). The bands correspond to varying the fraction of DM in microhalos between 0 . the pulsar or the Earth. The second type is a Shapiro time-delay if the clump traverses theline-of-sight between the Earth and the pulsar. Depending on how many clumps traverse thevicinity of a pulsar in the array, and how long it takes these signals can be further subdividedinto static, dynamic, and stochastic searches.While initial work indicated excellent prospects for constraining the properties of clumps,it was realised in Ref. [32] that some of the DM signal is absorbed by the fit of the pulsartiming model. Incorporating this effect leads to nearly an order of magnitude decrease in thesize of the Doppler signal, and a factor of around 4 in the Shapiro signal. The prospects fordetection at SKA become more challenging, and so [32, 33] also discuss the capabilities somemore futuristic PTAs. The most important parameters in the analysis of [31–33] are the PTAobservation period T , timing residuals t RMS , cadence ∆ T , distance d , and number of pulsars N P . The optimistic parameter sets considered in [31–33] primarily involve decreasing thetiming residual to order 10 ns and increasing the number of pulsars in the array to order 1000.There are also increases in the observation time to 30 years and decreasing the observationcadence to 1 week.Current PTAs include the Parkes Pulsar Timing Array [68] (20 pulsars), the European26ulsar Timing Array [69] (42 pulsars), and the NANOGrav 11-year data set [70] (45 pulsars).There is some overlap between the arrays, so only 73 of these are independent. Some of thepulsars were combined in the International Pulsar Timing Array data release 1 (IPTA-dr1) [71], which used 49 MSPs rather than 73. The more recent IPTA data release 2 [72](IPTA-dr2) includes 65 pulsars.The Square Kilometre Array (SKA) is a next-generation radio-telescope that will com-mence construction later in the current decade. Ref. [73] estimates that an SKA PTA couldinclude 200 pulsars with 50 ns timing residual and a cadence of 2 weeks. These parametersare consistent with other work [74, 75] - specifically, that SKA2 could discover up to 3000millisecond pulsars (MSPs), and that 5-10% of these MSPs are appropriate for use in tim-ing arrays. Results assuming a 20 year observation period at SKA would therefore appeararound 2050.The limitation on the number of pulsars is set by the total number of MSPs in the MW(estimated to be around 30,000 [76, 77]) and the fraction suitable for timing arrays. If thelatter is of order 5-10%, achieving a timing array with 1000 pulsars could require discoveringan order one fraction of the MSPs in the galaxy.The current best timing residuals are around 50 ns. Indeed, there is one pulsar in theNANOGrav-11 dataset with t RMS = 30 ns, although the average in that array is around 10times higher. The Parkes telescope can achieve residuals of 100 ns on only a few pulsars. Onecould scale the Parkes parameters to the recently commissioned FAST telescope, which wouldlead to the conclusion that FAST could achieve residuals of 1-10 ns on some pulsars [78].However, [78] indicates this target will be difficult to reach; a “realistic” timing array onFAST is considered to be 50 pulsars with residuals of around 100 ns.The main challenge to achieving low timing residuals is noise from a variety of sources [79].In particular, jitter noise on short time scales due to the intrinsic variability of the shape ofindividual pulsars is already becoming an issue, although there are proposals that may beable to deal with it. On longer time scales there are timing noise and irregularities (“rednoise”), and effects from the interstellar medium (ISM). Scattering off the ISM causes adispersive delay in the signal. For a static source this can be manageable, but millisecondpulsars tend to have relatively high velocities. Consequently over a 20-30 year observationtime this dispersion can have an impact. Whether all of these noise sources can be overcometo allow timing residuals of O (10) ns is as of yet unclear.Recently [33] has considered the prospects for the observation of dark matter substructure27rom a number of different dark matter scenarios, including early matter domination and darkphotons which we also consider here. For early matter domination and assuming optimisticPTA parameters reheating temperatures less than 1 GeV can be probed. For dark photondark matter and optimistic parameters PTAs may be sensitive to masses less than 10 − eV(note that for masses (cid:46) − eV inflationary production of dark photons cannot saturatethe observed DM relic density due to bounds on the scale of inflation [80]). C. Encounter Rate
While the previous two subsections focused on gravitational signatures of microhalos,DM substructure can also impact terrestrial direct detection searches if a non-gravitationalcoupling to Standard Model particles exists. If most of the dark matter in the galaxy isbound into clumps that formed at high redshift, then the encounter rate of DM constituentswith Earth can be drastically different from ΛCDM. We first estimate the encounter rateassuming (1) the distribution of the largest free clumps is dominated by a single characteristiccollapse redshift z c , and (2) the clump properties are well-approximated by the scaling lawsof Sec. II C.For a monochromatic (i.e., delta-function-like) microhalo mass distribution the encounterrate is τ ≡ Γ − = 1 nσv rel (48)where the Earth-microhalo relative velocity is v rel ∼ − , the local microhalo number densityis n = ρ DM /M ∗ ( z c ), and the cross section is σ = πR ∗ ( z c ) . Putting together the pieces, usingEq. (18) for R ∗ , and dropping O (1) numbers, we can express the rate as τ ≈ M ∗ ( z c ) ρ ( z c ) ρ DM v rel . (49)Here ρ ( z c ) and M ∗ ( z c ) are given by Eqs. (17) and (29) respectively, so in this approximationthe rate is a function of z c , the exponent n characterizing the power law enhancement, andthe mass ¯ M corresponding to the maximum enhanced scale in the matter power spectrum.A nontrivial distribution of microhalos masses is easily incorporated to give a differentialencounter rate d Γ enc d ln M = ρ DM M dfd ln M σv rel , (50)where df /d ln M can be estimated from EPS as in Eq. (40), and the cross-section σ implicitlydepends on M as above. 28 V. EXAMPLES
In this section we apply the results discussed in the previous sections to three specificmicrophysical scenarios: two EMD-like models, and dark photon DM produced through in-flationary fluctuations. In Figs. 5 and 6, we show power spectra and quantify the resultingmicrohalo properties (mass and density), comparing them to the sensitivity of caustic mi-crolensing. In Fig. 7 we estimate the microhalo distribution functions, which are then used tocompute lensing convergence power spectra (Figs. 8 and 9) and Earth-microhalo encounterrates (Fig. 10). Below, we give a brief summary of the models and these results, and then,for illustration, we explain in greater detail the features of early matter domination (arisingin the first two models) in subsection IV A.The first model, denoted “EMD+cutoff”, exhibits a k rise in the power spectrum begin-ning around k/k eq ∼ and ending in a small-scale cutoff. This power spectrum can begenerated by a period of early matter domination with a reheating temperature of around 10MeV and a Jeans length providing the short distance cutoff. The particular spectrum shownin Fig. 5 corresponds, for example, to an ALP dark matter particle with mass around 10 − eV [6]. In this case, the modified expansion history allows the relic density to be saturatedwithout fine-tuning the ALP misalignment angle [81]. The ρ s ( M s ) curve is well within theestimated sensitivity of caustic microlensing for clumps collapsing before z ∼ k rise in the power spectrum,followed by a plateau. In this case, the transition from ΛCDM to k again corresponds toEMD with a reheating temperature of around 10 MeV. The transition from k to the plateaucorresponds to the onset of EMD, prior to which the universe is assumed to be radiationdominated. (Logarithmic growth during RD is neglected in the toy model.) A short distancecutoff is typically present at higher k , but we extend the plateau to high scales to maximizethe difference with the EMD+cutoff model. An example power spectrum of this kind arisesin the cannibal models studied in [51]. The large scale behavior of the ρ s ( M s ) curve is thesame as that of the cutoff model, while the greater power at small scales leads to higherdensities at low masses, further strengthening the lensing sensitivity.The third model, denoted “Dark Photon” and shown in Fig. 6, exhibits a k rise inpower starting around k/k eq ∼ and reaching O (1) before descending in a shallow k − While the presence of a small scale cutoff is expected due to the Jeans scale, the steepness used here isonly a toy model, and it is probably less sharp in a complete treatment. Exponentially falling cutoffs canbe achieved, however, in particle DM models with a non-negligible free-streaming scale. ∼ − eV produced near the end ofinflation [12]. Although the boost in power is much larger than in the previous scenarios,it is isocurvature and occurs at smaller scales. Resultingly, the ρ s ( M s ) curve skirts thecaustic microlensing projection, only exhibiting significant enhancement below 10 − M (cid:12) .This suggests more refined analysis of the microlensing sensitivity is warranted.Fig. 7 shows the Milky Way subhalo mass distribution derived from extended Press-Schechter for each of the three models. These distributions are evaluated at high redshifts z = 250 and z = 100, in the left and right panels, respectively. These are particularlyrelevant because the rough estimate of clump survival from stellar encounters depends onlyon the clump density, which in turn depends only on the redshift of collapse. Therefore,by choosing z ∼ O (100) (and z = 0, M = M MW ) we obtain an estimate for the DMclump distribution in the galaxy today. This amounts to the coarse approximation that allmicrohalos forming later than z are disrupted, while those that form earlier all survive.Next, we evaluate lensing convergence power spectra. These are shown in Fig. 8 forthe EMD-inspired models and in Fig. 9 for the dark photon model. In each panel weshow the result of evaluating Eq. (46) using the EPS distribution with M = 10 M (cid:12) and z ∼ z . Since caustic microlensing is a statistical probe, it is potentiallysensitive to density inhomogeneities on a wide range of scales, including sub-substructurewhose existence is not captured by the EPS mass function and the use of a smooth NFWprofile in Eq. (45). We therefore show the convergence spectrum for z = 1000 and z = 100to illustrate its range. We note again that the high- z stellar disruption cutoff relevantfor the MW does not apply to the clusters relevant for caustic microlensing. Ref. [26]argued that convergence fluctuations with a magnitude ∆ κ (cid:38) − − − in the range q ∈ (10 − , .
1) AU − are potentially observable. We see that the denser, earlier-formingstructures result in larger fluctuation amplitude, but on smaller scales (larger q ). In eachpanel of Figs. 8 and 9 we compare the power spectra from realistic microhalo distributions tothose from monochromatic distributions centered at M = M ∗ ( z ), as evaluated in Eq. (47).We see that the latter provides a reasonable, order-of-magnitude approximation to the morerealistic distributions. Comparing the two EMD-like examples in Eq. (8), we note minutedifferences in the large- q behavior of the power spectra for the cutoff and plateau cases.The EMD+Plateau example features more power at large k which translates in a larger This is approximately the lowest dark photon mass for which inflationary production can saturate theobserved relic abundance of DM given current bounds on the scale of inflation. q in theconvergence PS. Note that the power spectra shown in Figs. 8 and 9 depend on microphysicalparameters which determine the position of the MPS peak and high- k cutoff; for example,increasing T RH in the EMD models, or the dark photon mass, shifts the MPS peak to larger k , which would also translate the convergence power spectra to higher q .We conclude this section with an estimate of Earth-microhalo encounter rates in ourbenchmark scenarios. The result of evaluating Eq. (50) with the EPS microhalo distributionfunctions in Fig. 7 (with M = M MW , z = 0, and z = 250) is shown in Fig. 10. We seethat enhanced substructure leads to Earth-microhalo encounter rates of ∼ once in 10 to 10 years. The higher rates correspond to power spectra with more power at large k , leadingto the formation of lighter microhalos (which have a large number density for fixed totalDM energy density). Since we expect that our treatment of disruption is conservative, truerates are likely to be even smaller, since more DM mass is locked up in somewhat heaviermicrohalos. A. EMD-induced Bumps in the MPS
The simplified power spectra in Fig. 5 are motivated by scenarios where the cosmologicalexpansion history is modified by a period of EMD, ending shortly before Big Bang Nucle-osynthesis (BBN). Here, for illustration, we sketch the enhancement generated by EMD ingreater detail.Various models can give rise to an EMD era with significant impact on the DM powerspectrum. Relic hidden radiation baths, for example, can come to temporarily dominate theenergy density of the universe once the lightest dark particle becomes nonrelativistic [1, 7].Alternatively, EMD can arise due to the nonthermal production of unstable states, as inmisalignment production of moduli [82, 83]. Obtaining enhanced small-scale structure froman EMD era requires DM to be kinetically decoupled from the SM radiation bath followingreheating [1]. While this condition is often not satisfied for WIMPs [3], it can easily besatisfied in models with nonthermal DM (e.g. [81]), or theories where DM arises from ahidden radiation bath that is thermally decoupled from the SM (e.g. [7, 11]).During EMD, perturbations grow linearly with scale factor between horizon entry andreheating, resulting in a total enhancement of a density perturbation by a factor a RH /a hor .Since the scale factor at horizon entry scales as a hor ∼ /k , EMD gives rise to a power31 � � �� � �� � �� � �� �� �� - �� �� - � �� - � �� - � �� � �� � �� � �� � �� �� �� - �� �� - � �� - � �� - � Figure 5. Upper left: Example power spectrum resulting from a period of EMD with 10 MeVreheating temperature and a sharp cutoff on small-scale power. Such a cutoff could arise, forexample, from the Jeans length in an ALP dark matter model. In this case the corresponding ALPmass is about a nano-eV. Upper right: The ρ s − M s relation for this model and the projectionsfor caustic microlensing. The dashed blue curve shows the relation for clumps that form from a2 σ fluctuation in the density contrast. Dashed grey lines indicate the collapse epoch z c . Lowerleft: Example power spectrum resulting from a period of EMD with 10 MeV reheating temperatureand a plateau at small scales. Approximate plateaus arise when any cutoff occurs on much smallerscales than the horizon size at the start of the modified expansion history. Lower right: The ρ s − M s relation for this model and the projections for caustic microlensing. spectrum that grows like k . The largest scales that can receive EMD enhancement depend32 - - Figure 6. Left: Example power spectrum for 10 − eV dark photon dark matter produced duringinflation. Right: The ρ s − M s relation for this model and the projections for caustic microlensing.The dashed blue curve shows the relation for clumps that form from a 2 σ fluctuation in the densitycontrast. Dashed grey lines indicate the collapse epoch z c . on the reheating temperature where EMD terminates. Both BBN and the Cosmic MicrowaveBackground (through N eff ) place lower bounds on this temperature, requiring T RH (cid:38) M ∗ ∼ ¯ k − ∼ T − .If T RH ∼
10 MeV – few hundred MeV, the resulting clumps are in the sensitivity region forcaustic microlensing.The smallest scales that can receive EMD enhancement are also indirectly affected by T RH . Modes that go nonlinear and collapse during EMD form halos that are dominantlycomposed of the metastable particle, and they are obliterated by the reheating process [10].Thus, for adiabatic perturbations that begin at the level of one part in 10 , there is aneffective upper limit on the amplitude of the EMD-induced enhancement of 10 for modesthat enter the horizon during EMD, or 10 for modes that enter during an earlier period ofRD (since modes grow by a factor of ∼
10 upon horizon entry during RD).The location of the peak is primarily determined by the small-scale cutoff, k cut . Herethere are a range of possibilities, depending on the microphysics of both DM and the physicsgenerating the EMD era. For example, free-streaming, collisional damping, and diffusiondamping are all typically modeled with a Gaussian cutoff, T ( k ) ∝ exp( − k /k cut ) [86, 87]. Forfree-streaming, k cut is given by the particle free-streaming horizon, k − cut = λ F S ; for collisional33 - - - - - - - - - - - - - - - - Figure 7. Sample subhalo distribution functions of microhalos in identified at z = 250 (leftpanel) or z = 100 (right panel) that end up in a MW-mass galaxy by z = 0. The three curvescorrespond the EMD and dark photon scenarios described in Sec. IV. The vertical dashed lines givethe characteristic mass M ∗ ( z ) in each model (defined in Eq. 16). The two EMD examples have T RH ≈
10 MeV, while the dark photon curves correspond to m A (cid:48) ≈ − eV. - - × - × - × - × - - - × - × - × - × - Figure 8. Convergence power spectra for the two EMD models (with a high- k cutoff, left, andwithout, right) with T RH ≈
10 MeV, evaluated using the subhalo distributions at z = 100 and z = 1000. Larger reheating temperatures shift the spectra to larger wavenumbers q . The solid anddotted curves assume an Extended Press-Schechter (EPS) and Delta Function subhalo mass distri-butions respectively. Convergence fluctuations with a magnitude ∆ κ (cid:38) − − − in the range q ∈ (10 − , .
1) AU − are potentially observable using microlensing of highly magnified stars [26]. - × - × - × - × - Figure 9. Convergence power spectra for the dark photon model with m A (cid:48) ≈ − eV, evaluatedusing the subhalo distributions at z = 100 and z = 1000. The solid and dotted curves assume anExtended Press-Schechter (EPS) and Delta Function subhalo mass distributions. Larger dark pho-ton masses shift the spectra to larger wavenumbers q . Convergence fluctuations with a magnitude∆ κ (cid:38) − − − in the range q ∈ (10 − , .
1) AU − are potentially observable using microlensingof highly magnified stars [26]. - - - - - - - - Figure 10. Differential encounter rate as a function of the subhalo mass for EMD-inspired PSenhancements with and without a small-scale cutoff, and the dark photon model described in thetext. m φ if DM is predominantly coupled to φ , where φ is themetastable species responsible for generating EMD [7, 10] (i.e., k cut = a ( T = m φ ) H ( T = m φ )in [10]); if DM is instead coupled to the SM radiation bath the temperature of kineticdecoupling, T kd , dictates collisional damping scale [3] (i.e., k cut = a ( T = T kd ) H ( T = T kd )).For another example, in the case of cannibal self-interactions, the cutoff is a more gradualfunction of k , with an envelope falling off as T ( k ) ∼ ( k cut /k ) [11].The nature of the peak – whether it is a relatively sharp feature or a plateau – depends onwhether k i is smaller or larger than k cut , where k i is the wavenumber of the mode that entersthe horizon at the onset of EMD. If k i > k cut , the feature is typically sharp, and k i playslittle role. If, however, k i < k cut , then modes with k cut > k > k i enter the horizon duringa prior period of radiation domination. They experience logarithmic growth, followed bylinear growth during EMD. The resulting feature in the power spectrum is a broad plateau,rather than a sharp bump, eventually cut off near k cut [51]. V. DISCUSSION
Features in the matter power spectrum arise in many different models of dark matterand early universe cosmology, and, as a result, they yield unique probes of both new particlephysics and the universe prior to BBN. We have endeavored to give a simplified, phenomeno-logical presentation of the map from enhancements in the power spectrum to the properties ofdark matter clumps at late times and their impact on gravitational observables, particularlycaustic microlensing [26]. In general, detailed features of the primordial power spectrum aresmoothed out in late-time observables, so that simple broken power law bumps in the MPSare sufficient to model more physical enhancements. The small- k power law, and whetherthe perturbations are isocurvature or adiabatic, determines many features of the largest,most observable clumps, if they survive until late times. Increasing the slope of the powerlaw increases the collapse redshift for a given scale, which in turn increases the density anddecreases the radius of the resulting collapsed objects. Extending the power law to smaller k results in larger, more massive, less dense clumps. If the power law extends to large enough k , neither the peak nor the cutoff strongly affects observable prospects.For more localized enhancements, the details of the peak and cutoff affect the propertiesof the densest clumps. These details are important in some cases and not in others. Forexample, for caustic microlensing, differences in the MPS cutoff slope lead only to minute36ifferences in the shape of the convergence power spectrum (its amplitude is also modified,but this is degenerate with the fraction of DM in microhalos). In fact, in our treatmentthis observable appears to be well approximated using a monochromatic microhalo massdistribution which only depends on the location of the MPS peak and its falloff at small k .In contrast, Earth-microhalo encounter rates can be sensitive to the very small-scale behaviorof the MPS, strongly influencing direct detection prospects. This is because shallower slopesin the MPS beyond the peak lead to a large abundance of light microhalos, thereby enhancingthe probability of an encounter.Our analysis relies on a number of assumptions and simplifications which provide numer-ous directions for future development. We conclude by enumerating some of these points.We have only considered clumps that collapse after matter-radiation equality. However,there are other possibilities which lead to structures forming at earlier epochs. One way thiscan happen is collapse during the radiation-dominated period that precedes matter-radiationequality. The properties of these clumps depend on the free-streaming length of the darkmatter, and hence precisely how cold it is. Alternatively, there may be an epoch of non-standard cosmological evolution such as early matter domination (EMD), when structuresmay also collapse. Clumps that collapse during EMD may disappear at reheating due to thedecay of the field responsible for the EMD. This introduces an aspect of model-dependencewhich is not captured by our phenomenological approach. A detailed study of structureformation where both possibilities occur was performed in Ref. [10].We have focused on clumps whose masses lie within a few orders of magnitude of the massof the Earth. Partly this is because that is the regime of maximum sensitivity for causticmicrolensing. Pulsar timing searches also have sensitivity in this region and at higher masses.It is unclear what searches would be sensitive for clumps with masses substantially less thanthis. The existence of lower mass clumps is also more sensitive to UV physics and themechanism which determines the relic density. For example, ALP clumps from a period ofEMD have a lower bound on their masses associated with the size of the mass inside thehorizon when the ALP field starts oscillating [6].To predict the properties of the clumps we rely on the Press-Schechter formalism. This isbased on assumptions of Gaussian fluctuations, which may not hold for the large enhance-ments of the power spectra we consider. The Press-Schechter picture tells us that larger halosare formed by the hierarchical merges of smaller clump halos. The detailed treatment of thisprocess, and the ultimate fate of the clumps involved requires N-body simulations which37re particularly challenging due to the large hierarchies of scales involved. Consequently inour calculations of microlensing sensitivity we assume either a “monochromatic” microhalomass function or one derived from a high-redshift snapshot of an extended Press-Schechtermass function; the true late-time distribution is certainly much more complex. We assumean NFW-type profile for the clumps. Studies of axion miniclusters [88] and ultra-compactmicrohalos [43, 44] both support NFW (or slightly-cuspier-than-NFW) profiles for the firstforming halos.Our projection for the caustic microlensing sensitivity in Fig. 4 is quite simplified, andbased on [6, 20]. We have assumed the existence of a single “typical” cluster for our projec-tions. This could be improved by considering an ensemble of clusters whose properties arealready known. We also have not considered the impact of astrophysical uncertainties andpotential backgrounds inherent in the cluster, such as clumps of gas, although these effectshave been argued elsewhere to be unimportant [26].Finally, we do not address the inverse problem. The extent to which any details of theclump mass distribution and profiles (and the early universe physics responsible) can bereconstructed from extended observations using future telescopes is an interesting questionwhich we leave for future work. ACKNOWLEDGMENTS
We thank Adrienne Erickcek and Antonella Palmese for useful conversations. Thismanuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of HighEnergy Physics. MJD is supported by the Australian Research Council. The work of PD issupported by the US Department of Energy under Grant No. DE-SC0015655. The work ofJS is supported by the US Department of Energy under Grant No. DE-SC0017840.
Appendix A: Window Function Dependence
The semi-analytical results in the body of the paper depend on the choice of the windowfunction W used to filter the density contrast – see the definition of the density variance σ in Eq. (15). All subsequent results, such as typical collapsing masses M ∗ (Eq. (16))and the subhalo distribution function depend on this choice. In certain cases, the window38unction choice is physically or theoretically motivated: the real space top hat filter scaleis straightforward to interpret in terms of a physical length/mass scale [64], while a Fouriertop hat filter simplifies calculations in the extended Press-Schechter formalism [63]. In thecontext of sharply-rising dimensionless power spectra considered here, the position-space tophat can lead to divergences when computing σ , while the Fourier top hat always yields finiteresults. It is therefore important to check whether different choices of window functions canchange our conclusions. In this Appendix we show that this is not the case – qualitativelysimilar results are obtained for many different choices of W .We consider three different window functions: a spherical top hat (TH) with a Gaussiancut off at large k , a Fourier space top hat (FTH) and a Gaussian (G): W TH ( kR ) = 3(sin kR − kR cos kR )( kR ) e − α ( kR ) / (A1a) W FTH ( kR ) = Θ(1 − kR ) (A1b) W G ( kR ) = e − ( kR ) / , (A1c)where we take α = 10 − following Ref. [1]. The mass scale corresponding to R is then givenby M = (4 π/ , π , (2 π ) / ) × ¯ ρ m, R for TH, FTH and Gaussian filters, respectively. In Fig. 11 we show the result of varying W for two representative MPS enhancementsinspired by early matter domination: one has a sharp cutoff at k = k peak , and another aplateau for k > k peak . These are described in Sec. IV. Different choices of W give estimates ofthe fluctuation variance and typical collapse mass that differ at most by a factor of ∼ − [1] A. L. Erickcek and K. Sigurdson, Reheating Effects in the Matter Power Spectrum andImplications for Substructure , Phys. Rev.
D84 (2011) 083503, [ ].[2] G. Barenboim and J. Rasero,
Structure Formation during an early period of matterdomination , JHEP (2014) 138, [ ]. The equality is only approximate for TH with a Gaussian cutoff. Moreover, the normalizations of theFTH and Gaussian filters are ambiguous – see, e.g., Refs. [64, 89, 90] for discussions of these issues. - - - - - - - - - - Figure 11. Impact of the window function choice on density fluctuation variance (left column,defined in Eq. 15) and typical mass collapsing at z (right column, defined in Eq. 16). The two rowscorrespond to Early Matter Domination-inspired models described in Sec. IV, with and withouta sharp cutoff at small scales. The solid, dashed and dotted lines correspond to the Top Hat,Fourier Top Hat and Gaussian window functions defined in Eq. A1. In the left column we havefixed z = 1000.[3] J. Fan, O. ¨Ozsoy and S. Watson, Nonthermal histories and implications for structureformation , Phys. Rev.
D90 (2014) 043536, [ ].[4] L. Visinelli and J. Redondo,
Axion Miniclusters in Modified Cosmological Histories , Phys.Rev. D (2020) 023008, [ ].[5] A. E. Nelson and H. Xiao,
Axion Cosmology with Early Matter Domination , Phys. Rev.
D98 (2018) 063516, [ ].
6] N. Blinov, M. J. Dolan and P. Draper,
Imprints of the Early Universe on Axion Dark MatterSubstructure , Phys. Rev. D (2020) 035002, [ ].[7] Y. Zhang,
Long-lived Light Mediator to Dark Matter and Primordial Small Scale Spectrum , JCAP (2015) 008, [ ].[8] J. A. Dror, E. Kuflik and W. H. Ng,
Codecaying Dark Matter , Phys. Rev. Lett. (2016)211801, [ ].[9] J. A. Dror, E. Kuflik, B. Melcher and S. Watson,
Concentrated dark matter: Enhancedsmall-scale structure from codecaying dark matter , Phys. Rev. D (2018) 063524,[ ].[10] C. Blanco, M. S. Delos, A. L. Erickcek and D. Hooper, Annihilation Signatures of HiddenSector Dark Matter Within Early-Forming Microhalos , Phys. Rev. D (2019) 103010,[ ].[11] A. L. Erickcek, P. Ralegankar and J. Shelton,
Cannibal domination and the matter powerspectrum , .[12] P. W. Graham, J. Mardon and S. Rajendran, Vector Dark Matter from InflationaryFluctuations , Phys. Rev. D (2016) 103520, [ ].[13] Y. Ema, K. Nakayama and Y. Tang, Production of Purely Gravitational Dark Matter: TheCase of Fermion and Vector Boson , JHEP (2019) 060, [ ].[14] A. Ahmed, B. Grzadkowski and A. Socha, Gravitational production of vector dark matter , JHEP (2020) 059, [ ].[15] E. W. Kolb and A. J. Long, Completely Dark Photons from Gravitational ParticleProduction During Inflation , .[16] E. W. Kolb and I. I. Tkachev, Large amplitude isothermal fluctuations and high density darkmatter clumps , Phys. Rev.
D50 (1994) 769–773, [ astro-ph/9403011 ].[17] M. Fairbairn, D. J. E. Marsh, J. Quevillon and S. Rozier,
Structure formation andmicrolensing with axion miniclusters , Phys. Rev. D (2018) 083502, [ ].[18] A. Vaquero, J. Redondo and J. Stadler, Early seeds of axion miniclusters , JCAP (2019)012, [ ].[19] M. Buschmann, J. W. Foster and B. R. Safdi, Early-Universe Simulations of theCosmological Axion , Phys. Rev. Lett. (2020) 161103, [ ].[20] A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson and K. Van Tilburg,
The Large-Misalignment Mechanism for the Formation of Compact Axion Structures: ignatures from the QCD Axion to Fuzzy Dark Matter , .[21] A. L. Erickcek, The Dark Matter Annihilation Boost from Low-Temperature Reheating , Phys.Rev.
D92 (2015) 103505, [ ].[22] K. Redmond, A. Trezza and A. L. Erickcek,
Growth of Dark Matter Perturbations duringKination , Phys. Rev. D (2018) 063504, [ ].[23] K. Kadota and J. Silk, Boosting small-scale structure via primordial black holes andimplications for sub-GeV dark matter annihilation , .[24] D. Croon, D. McKeen and N. Raj, Gravitational microlensing by dark matter in extendedstructures , Phys. Rev. D (2020) 083013, [ ].[25] D. Croon, D. McKeen, N. Raj and Z. Wang,
Subaru-HSC through a different lens:Microlensing by extended dark matter structures , Phys. Rev. D (2020) 083021,[ ].[26] L. Dai and J. Miralda-Escud´e,
Gravitational Lensing Signatures of Axion Dark MatterMinihalos in Highly Magnified Stars , Astron. J. (2020) 49, [ ].[27] E. R. Siegel, M. Hertzberg and J. Fry,
Probing Dark Matter Substructure with PulsarTiming , Mon. Not. Roy. Astron. Soc. (2007) 879, [ astro-ph/0702546 ].[28] S. Baghram, N. Afshordi and K. M. Zurek,
Prospects for Detecting Dark Matter HaloSubstructure with Pulsar Timing , Phys. Rev. D (2011) 043511, [ ].[29] K. Kashiyama and M. Oguri, Detectability of Small-Scale Dark Matter Clumps with PulsarTiming Arrays , .[30] H. A. Clark, G. F. Lewis and P. Scott, Investigating dark matter substructure with pulsartiming – I. Constraints on ultracompact minihaloes , Mon. Not. Roy. Astron. Soc. (2016)1394–1401, [ ].[31] J. A. Dror, H. Ramani, T. Trickle and K. M. Zurek,
Pulsar Timing Probes of PrimordialBlack Holes and Subhalos , Phys. Rev. D (2019) 023003, [ ].[32] H. Ramani, T. Trickle and K. M. Zurek,
Observability of Dark Matter Substructure withPulsar Timing Correlations , .[33] V. S. Lee, A. Mitridate, T. Trickle and K. M. Zurek, Probing Small-Scale Power Spectra withPulsar Timing Arrays , .[34] C.-P. Ma and E. Bertschinger, Cosmological perturbation theory in the synchronous andconformal Newtonian gauges , Astrophys. J. (1995) 7–25, [ astro-ph/9506072 ].[35] W. Hu,
Structure formation with generalized dark matter , Astrophys. J. (1998) 485–494, astro-ph/9801234 ].[36] W. Hu and N. Sugiyama, Small scale cosmological perturbations: An Analytic approach , Astrophys. J. (1996) 542–570, [ astro-ph/9510117 ].[37] P. Meszaros,
The behaviour of point masses in an expanding cosmological substratum , Astron.Astrophys. (1974) 225–228.[38] Planck collaboration, N. Aghanim et al.,
Planck 2018 results. VI. Cosmological parameters , .[39] W. H. Press and P. Schechter, Formation of galaxies and clusters of galaxies by selfsimilargravitational condensation , Astrophys. J. (1974) 425–438.[40] H. Mo, F. C. van den Bosch and S. White,
Galaxy Formation and Evolution . CambridgeUniversity Press, 2010.[41] P. J. E. Peebles,
The large-scale structure of the universe . Princeton University Press, 1980.[42] S. Cole and C. G. Lacey,
The Structure of dark matter halos in hierarchical clusteringmodels , Mon. Not. Roy. Astron. Soc. (1996) 716, [ astro-ph/9510147 ].[43] M. S. Delos, A. L. Erickcek, A. P. Bailey and M. A. Alvarez,
Are ultracompact minihalosreally ultracompact? , Phys. Rev. D (2018) 041303, [ ].[44] M. S. Delos, A. L. Erickcek, A. P. Bailey and M. A. Alvarez, Density profiles of ultracompactminihalos: Implications for constraining the primordial power spectrum , Phys. Rev. D (2018) 063527, [ ].[45] J. F. Navarro, C. S. Frenk and S. D. White, The Structure of cold dark matter halos , Astrophys. J. (1996) 563–575, [ astro-ph/9508025 ].[46] M. A. S´anchez-Conde and F. Prada,
The flattening of the concentration–mass relationtowards low halo masses and its implications for the annihilation signal boost , Mon. Not.Roy. Astron. Soc. (2014) 2271–2277, [ ].[47] A. Cooray and R. K. Sheth,
Halo Models of Large Scale Structure , Phys. Rept. (2002)1–129, [ astro-ph/0206508 ].[48] R. H. Wechsler, J. S. Bullock, J. R. Primack, A. V. Kravtsov and A. Dekel,
Concentrationsof dark halos from their assembly histories , Astrophys. J. (2002) 52–70,[ astro-ph/0108151 ].[49] D. H. Lyth and A. R. Liddle,
The primordial density perturbation: Cosmology, inflation andthe origin of structure . Cambridge, UK: Cambridge Univ. Pr. (2009) 497 p, 2009.[50] D. J. H. Chung, E. W. Kolb, A. Riotto and L. Senatore,
Isocurvature constraints on ravitationally produced superheavy dark matter , Phys. Rev. D (2005) 023511,[ astro-ph/0411468 ].[51] A. Erickcek, P. Ralegankar and J. Shelton, “Cannibal Cutoffs in Early Matter-DominatedEras.” to appear.[52] F. C. van den Bosch, G. Ogiya, O. Hahn and A. Burkert, Disruption of Dark MatterSubstructure: Fact or Fiction? , Mon. Not. Roy. Astron. Soc. (2018) 3043–3066,[ ].[53] M. Cautun, A. Ben´ıtez-Llambay, A. J. Deason, C. S. Frenk, A. Fattahi, F. A. G´omez et al.,
The milky way total mass profile as inferred from Gaia DR2 , MNRAS (Apr., 2020)4291–4313, [ ].[54] T. Goerdt, O. Y. Gnedin, B. Moore, J. Diemand and J. Stadel,
The survival and disruptionof CDM micro-haloes: Implications for direct and indirect detection experiments , Mon. Not.Roy. Astron. Soc. (2007) 191–198, [ astro-ph/0608495 ].[55] A. Schneider, L. Krauss and B. Moore,
Impact of Dark Matter Microhalos on Signatures forDirect and Indirect Detection , Phys. Rev. D (2010) 063525, [ ].[56] V. Dokuchaev, Y. Eroshenko and I. Tkachev, Destruction of axion miniclusters in theGalaxy , J. Exp. Theor. Phys. (2017) 434–442, [ ].[57] P. Tinyakov, I. Tkachev and K. Zioutas,
Tidal streams from axion miniclusters and directaxion searches , JCAP (2016) 035, [ ].[58] A. Pillepich et al., First results from the IllustrisTNG simulations: the stellar mass contentof groups and clusters of galaxies , Mon. Not. Roy. Astron. Soc. (2018) 648–675,[ ].[59] P. L. Kelly et al.,
Extreme magnification of an individual star at redshift 1.5 by agalaxy-cluster lens , Nature Astron. (2018) 334–342, [ ].[60] R. K. Sheth and G. Tormen, Large scale bias and the peak background split , Mon. Not. Roy.Astron. Soc. (1999) 119, [ astro-ph/9901122 ].[61] M. S. Delos, M. Bruff and A. L. Erickcek,
Predicting the density profiles of the first halos , Phys. Rev. D (2019) 023523, [ ].[62] H. Xiao, I. Williams and M. McQuinn,
Simulations of Axion Minihalos , .[63] J. R. Bond, S. Cole, G. Efstathiou and N. Kaiser, Excursion set mass functions forhierarchical Gaussian fluctuations , Astrophys. J. (1991) 440.[64] C. G. Lacey and S. Cole,
Merger rates in hierarchical models of galaxy formation , Mon. Not. oy. Astron. Soc. (1993) 627–649.[65] A. R. Zentner, The Excursion Set Theory of Halo Mass Functions, Halo Clustering, andHalo Growth , Int. J. Mod. Phys.
D16 (2007) 763–816, [ astro-ph/0611454 ].[66] J. M. Diego et al.,
Dark Matter under the Microscope: Constraining Compact Dark Matterwith Caustic Crossing Events , Astrophys. J. (2018) 25, [ ].[67] M. Oguri, J. M. Diego, N. Kaiser, P. L. Kelly and T. Broadhurst,
Understanding causticcrossings in giant arcs: characteristic scales, event rates, and constraints on compact darkmatter , Phys. Rev. D (2018) 023518, [ ].[68] D. Reardon et al., Timing analysis for 20 millisecond pulsars in the Parkes Pulsar TimingArray , Mon. Not. Roy. Astron. Soc. (2016) 1751–1769, [ ].[69] G. Desvignes et al.,
High-precision timing of 42 millisecond pulsars with the European PulsarTiming Array , Mon. Not. Roy. Astron. Soc. (2016) 3341–3380, [ ].[70]
NANOGrav collaboration, Z. Arzoumanian et al.,
The NANOGrav 11-year Data Set:High-precision timing of 45 Millisecond Pulsars , Astrophys. J. Suppl. (2018) 37,[ ].[71] J. Verbiest et al.,
The International Pulsar Timing Array: First Data Release , Mon. Not.Roy. Astron. Soc. (2016) 1267–1288, [ ].[72] B. Perera et al.,
The International Pulsar Timing Array: Second data release , Mon. Not.Roy. Astron. Soc. (2019) 4666–4687, [ ].[73] P. A. Rosado, A. Sesana and J. Gair,
Expected properties of the first gravitational wave signaldetected with pulsar timing arrays , Mon. Not. Roy. Astron. Soc. (2015) 2417–2433,[ ].[74] E. Keane et al.,
A Cosmic Census of Radio Pulsars with the SKA , PoS
AASKA14 (2015)040, [ ].[75] B. Stappers, E. Keane, M. Kramer, A. Possenti and I. Stairs,
The prospects of pulsar timingwith new-generation radio telescopes and the Square Kilometre Array , Phil. Trans. Roy. Soc.Lond. A (2018) 20170293.[76] D. Lorimer,
Binary and Millisecond Pulsars , Living Rev. Rel. (2008) 8, [ ].[77] D. R. Lorimer, The Galactic Millisecond Pulsar Population , IAU Symp. (2013) 237–242,[ ].[78] G. Hobbs, S. Dai, R. Manchester, R. Shannon, M. Kerr, K. Lee et al.,
The Role of FAST inPulsar Timing Arrays , Res. Astron. Astrophys. (2019) 020, [ ].
79] G. Janssen et al.,
Gravitational wave astronomy with the SKA , PoS
AASKA14 (2015) 037,[ ].[80]
Planck collaboration, Y. Akrami et al.,
Planck 2018 results. X. Constraints on inflation , Astron. Astrophys. (2020) A10, [ ].[81] N. Blinov, M. J. Dolan, P. Draper and J. Kozaczuk,
Dark matter targets for axionlikeparticle searches , Phys. Rev. D (2019) 015049, [ ].[82] T. Banks, D. B. Kaplan and A. E. Nelson,
Cosmological implications of dynamicalsupersymmetry breaking , Phys. Rev. D (1994) 779–787, [ hep-ph/9308292 ].[83] B. de Carlos, J. Casas, F. Quevedo and E. Roulet, Model independent properties andcosmological implications of the dilaton and moduli sectors of 4-d strings , Phys. Lett. B (1993) 447–456, [ hep-ph/9308325 ].[84] P. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti,
Bounds on verylow reheating scenarios after Planck , Phys. Rev. D (2015) 123534, [ ].[85] T. Hasegawa, N. Hiroshima, K. Kohri, R. S. Hansen, T. Tram and S. Hannestad, MeV-scalereheating temperature and thermalization of oscillating neutrinos by radiative and hadronicdecays of massive particles , JCAP (2019) 012, [ ].[86] A. M. Green, S. Hofmann and D. J. Schwarz, The First wimpy halos , JCAP (2005)003, [ astro-ph/0503387 ].[87] E. Bertschinger,
The Effects of Cold Dark Matter Decoupling and Pair Annihilation onCosmological Perturbations , Phys. Rev. D (2006) 063509, [ astro-ph/0607319 ].[88] B. Eggemeier, J. Redondo, K. Dolag, J. C. Niemeyer and A. Vaquero, First Simulations ofAxion Minicluster Halos , Phys. Rev. Lett. (2020) 041301, [ ].[89] A. J. Benson, A. Farahi, S. Cole, L. A. Moustakas, A. Jenkins, M. Lovell et al.,
Dark matterhalo merger histories beyond cold dark matter - I. Methods and application to warm darkmatter , MNRAS (Jan., 2013) 1774–1789, [ ].[90] A. Schneider, R. E. Smith and D. Reed,
Halo Mass Function and the Free Streaming Scale , Mon. Not. Roy. Astron. Soc. (2013) 1573, [ ].].