Axion minicluster streams in the Solar neighbourhood

A consequence of QCD axion dark matter being born after inflation is the emergence of small-scale substructures known as miniclusters. Although miniclusters merge to form minihalos, this intrinsic granularity is expected to remain imprinted on small scales in our galaxy, leading to potentially damaging consequences for the campaign to detect axions directly on Earth. This picture, however, is modified when one takes into account the fact that encounters with stars will tidally strip mass from the miniclusters, creating pc-long tidal streams that act to refill the dark matter distribution. Here we ask whether or not this stripping rescues experimental prospects from the worst-case scenario in which the majority of axions remain bound up in unobservably small miniclusters. We find that the density sampled by terrestrial experiment on mpc-scales will be, on average, around 70–90% of the average local DM density, and at a typical point in the solar neighbourhood, we expect most of the dark matter to be comprised of debris from O (10 2 –10 3 ) overlapping streams. If haloscopes can measure the axion signal with high-enough frequency resolution, then these streams are revealed in the form of an intrinsically spiky lineshape, in stark contrast with the standard assumption of a smooth, featureless Maxwellian distribution—a unique prediction that constitutes a way for experiments to distinguish between pre and post-inflationary axion cosmologies.

It has been known for many years [24][25][26] that the complicated multi-scale dynamics of the axion field that necessitate the constriction of these simulations, also implies that the DM distribution in this scenario will inherit inhomogeneities on scales set by the horizon at the QCD phase transition.This results in the majority of DM becoming bound inside planetary-mass structures called miniclusters [24][25][26][27][28][29][30][31][32][33][34][35][36][37].To put it plainly: miniclusters would be too sparsely distributed for us to have a realistic chance of encountering one, so prospects for detecting DM axions in the lab rest upon whether or not miniclusters survive in our galaxy today.
In this letter, we address the detectability of DM axions by quantifying their distribution in the solar neighbourhood.Broadly speaking, DM in this scenario can be thought of in terms of three distinct populations.Firstly, there are the axions that never end up inside miniclusters to begin with, existing instead in the "minivoids" between them [38].Secondly, there are the miniclusters themselves, which lock up more than 80% of the mass of DM prior to galaxy formation [39].Lastly, there is the minicluster debris-axions tidally stripped from their hosts as they orbit the Milky Way (MW) [40][41][42][43][44].The worst-case scenario for direct detection experiments (referred to as "haloscopes") is if we consider the axions only from the minivoids [38].However, tidal disruption will act to spread the DM across a wider volume [40], so observational prospects could be rescued if we account for this.
To that end, we have performed Monte-Carlo simulations that begin with a realistic population of miniclusters taken from early-Universe simulations, and are then evolved as they orbit the MW.For each minicluster, we calculate how much mass is stripped from it, and over what length scale this mass is spread so that we can build a model for the resulting stream network.Figure 1 shows an illustration of the stages of this calculation.We then use these results to create example haloscope signals, as in Fig. 2. Initial miniclusters.-Tostart our calculation, we need to know the initial distribution of minicluster masses and density profiles.We source this information from the most recent N-body simulations of minicluster formation, which themselves begin from initial conditions left over from lattice simulations of axions around the QCD phase transition [45].The most pertinent result of these simulations is the two types of miniclusters that form-what we refer to as merged miniclusters, and isolated miniclusters.The former result from hierarchical merging and develop Navarro-Frenk-White (NFW) density profiles [45][46][47], while the latter form from the prompt collapse of isolated overdensities and then do not undergo many substantial mergers-retaining the powerlaw profiles usually associated with self-similar collapse, ρ ∝ r − α [27].In our baseline set of results, we adopt α = 2.71 taken from fitting their density profiles at the latest times in our simulation [45].
Monte Carlo simulation.-Afterdrawing a sample of miniclusters, we propagate their orbits around the galaxy.We are interested here in the set of all possible galactic orbits that end in the Solar neighbourhood today.Velocities are sampled according to the Standard Halo Model (SHM) by drawing from an isotropic 3D Gaussian with width σ halo = v circ / √ 2 where v circ = 233 km/s is the circular speed of the Solar orbit [66].We truncate the velocity distribution at v esc = 528 km/s [67] to discount orbits unbound to the MW.We integrate orbits over a duration t MW = 13.5 Gyr using galpy [68], adopting the commonly-used potential "MWPotential2014".
For each orbit, we evaluate the variation in the local stellar number density n ⋆ (t) felt by the minicluster.Our galactic model includes a central bulge and the thin/thick disks-described in S.M. II.The total number of stars encountered by the minicluster is the integral, where v(t) is the minicluster velocity, and b max some maximum impact parameter between the minicluster and a star.We want to set b max to be as large as possible to capture all possible disruption but without being so large as to add unnecessary computational burden by having large numbers of negligible encounters.As in Ref. [42], we find that b max = 0.1 pc strikes a good balance.Increasing this number by a factor of 10 does not change our results.The N enc encounters are labelled by a set of encounter times {t i enc }, drawn from a probability distribution proportional to the integrand in Eq.(1).The impact parameter, b of each encounter, are drawn randomly from inside a circle of radius b max .
Stars appear to miniclusters as point-like objects, which means we may work in the distant tide approximation, where the impact parameter between the minicluster and a star is much larger than the minicluster radius, b ≫ R. The energy injected by an encounter with a star of mass M ⋆ under this approximation is [69,70], where v rel is the minicluster-star relative velocity, and ⟨R 2 ⟩ is the minicluster mean-squared radius (see S.M. III for how the latter is calculated).Most encounters inject ∆E ≪ E b so we must deal with perturbations, which do not totally disrupt the minicluster in one go, but lead to a series of mass losses ∆M i .The procedure then is to execute i = 1, ..., N enc successive perturbations over the minicluster's orbit, repeatedly updating the mass and radius that it relaxes to, until it is either fully unbound or we reach the end of the orbit.The recipe for this procedure is given in S.M. III.For com-putational efficiency, we group together the very large number of encounters occurring during a disk-crossing event, and only allow the minicluster to relax to a new profile and radius over its relaxation time between diskcrossings.This is justified here because the relaxation time [t rel ∼ (Gρ mc ) −1/2 ∼ O(10 Myr)] is shorter than the timescale between disk-crossings where major encounters occur: O(10 − 100 Myr).We find that isolated miniclusters are relatively stable against disruption, and the majority survive with some mass intact.The high-mass merged miniclusters are much more easily disrupted, with around half of them losing > 99% of their mass by today.Stream formation.-Whena certain amount of the minicluster's mass is stripped, where does it go?Although this unbound mass shell has effectively evaporated off of the minicluster, it will still retain its host's centre-of-mass orbital velocity, and so continue along the same orbit in the vicinity.The tidal field of the Milky Way will continue to act on these unbound particles, causing the debris to elongate into a stream in the direction of the orbit.This is seen generically for tidally disrupted structures at all scales in astrophysics, but has also been simulated for DM microhalos on the scales we are interested in here in Ref. [70].We use this study as inspiration for our stream model.We assume that an unbound mass shell turns into a tidal stream, with a leading tail that advances beyond the original host minicluster and a trailing tail that lags behind.By today, a given mass shell evaporating in some encounter at time t i enc will have turned into a stream of minimum length given by the minicluster velocity dispersion [70], It is important to note that the stream will likely be longer than this today due to the energy injection and further tidal heating, potentially by a factor of 10 compared to ∼ σ mc t [70].We explore levels of additional elongation in S.M IV, but for reasons that will become clear, taking the smallest possible length that the streams could grow to is the most conservative option for the question we are trying to answer.
To model the stream formation, we will make the assumption that the mass lost in each tidal disruption event goes into a cylinder the same radius as the original minicluster, R, and with a Gaussian density profile running along the stream.The total stream is then comprised of the sum of all of these individual evaporated mass shells occurring at different times.We can parameterise this in terms of a function ρ str (ℓ), where ℓ is a coordinate that runs along the stream.The length scales for the miniclusters (weighted by the ∆M in each segment) are in the range 7.3 +3.2 −4.0 pc for merged miniclusters and 0.28 +0.34  −0.2 pc for isolated miniclusters.
The local density in minicluster streams.-Wewill now present an argument for the degree to which the local DM density at our position in the galaxy is replenished by the disruption of the miniclusters into long streams.First, imagine a sphere around the Sun of radius 100 pc.This is the order-of-magnitude scale within which we have strong evidence for a galactic DM density of ρ DM ≈ 0.4 GeV/cm 3 ≈ 0.01 M ⊙ pc −3 [71].Within this sphere, the total mass of DM is M DM = 4.2 × 10 4 M ⊙ .From Ref. [38], we know around f void = 8% of this should be comprised of an ambient density of unbound axionsthe minivoids.Therefore M DM (1 − f void ) is the total mass of axions that were initially bound in miniclusters inside this volume.Using our baseline mass function, this implies a total of around N mc ∼ 10 14 isolated+merged miniclusters.
Let us assume the final volume occupied by each stream after the disruption process is V str ≈ πR 2 mc l 95 .We define l 95 to be the length within which 95% of the stream mass is contained, calculated using Eq.( 3).If we have N mc miniclusters inside a volume V = 4/3πr 3 local and those miniclusters each occupy a volume V str after disruption, the expected number of streams overlapping a random point inside the sphere is then, where the ratio of the two volumes gives us the probability that a particular stream overlaps our chosen point.
For our baseline set of assumptions, we find that this number is 246 ± 15 (statistical error).Varying our assumptions-for example the density profiles, the concentration of the NFW merged miniclusters, and/or accounting for additional stream heating-we obtain numbers in the range ∼ 200 to 6000, see S.M. IV.
With this number, we can now re-sample from our distribution of streams O(10 2 − 10 3 ) times, with probabilities weighted by V i str .For each stream, we also randomly choose our position within it ℓ ∈ [−l 95 /2, l 95 /2] to get the value of the DM density that it contributes from Eq. ( 3).Repeating this process many times, we find that the sum of all the individual densities adds up to a total where ρ DM ≈ 0.4 GeV/cm 3 is the usual DM density inferred on much larger scales.So when added to the 8% of the DM density originally filling the minivoids, this is a replenishment of almost 90% of ρ DM -a substantial improvement in the prospects for direct detection.
We emphasise here that this final number is generally insensitive to many of our cruder assumptions.The first point to state is that the debris from disrupted miniclusters is dominated almost entirely by the merged ones.
These streams constitute the most DM by mass, and have the highest probability of intersecting our position because they cover more volume.
The second key point is that the expected value of the DM density (Eq.5) is robust against re-parameterisations as long as we are in the regime when N str ≫ 1.That said, the number of expected streams is somewhat arbitrary, but given that it is very safely ≫ 1 for any reasonable choice of parameters, we show explicitly in S.M. IV that the typical (i.e.median, or mean) DM density they add up to does not depend on this number because we have conserved the total DM mass.Instead, the number of streams affects the variance in ρ str .It can also be shown (see S.M. IV) that as long as the minicluster mass and radius are related like R ∝ M 1/3 then the mass-dependence drops out entirely, meaning assumptions about the mass function also do not affect N str , leaving only a dependence on the NFW concentration parameter.
Quantities that affect N str include the typical radii that the miniclusters are truncated at (related to the concentration parameter for NFW profiles), and the extent to which streams are additionally elongated beyond the minimal expectation σ mc t.If we allow the miniclusters to be larger in radius, or the streams to be longer, then the typical number of streams overlapping at our position increases to O( 103 )-however this means that the variance in the mean density our position ends up smaller.Because of this, our baseline parameterisation is the most conservative-allowing for processes that can further strip the miniclusters or elongate streams only acts to suppress the variation in the value stated in Eq. (5). Figure S3 in the S.M. shows the impact of these uncertainties quantitatively.
Put together, the only uncertainty remaining that could change our results substantially is if the merged miniclusters do not continue to evolve towards smooth NFW profiles: If instead some of them remain as "clusters of miniclusters", then it is possible that they are more resilient.We discuss this issue in S.M. V.The takeaway from our N-body simulation halo finder is that the majority of the axions inside merged miniclusters are indeed attached to their host halo as opposed to internal subhalos.Impact on axion haloscopes.-Wenow know how many streams we expect to have around us at any one time, and so we can draw a sample of this size and predict how this would look in a haloscope.First, we build the velocity distribution of axions.In the voids, we can safely assume the halo is described by the smooth, fully-virialised SHM [66] modelled as an isotropic Gaussian with width σ void = v circ / √ 2 ≈ 165 km/s.An experiment observes this after a Galilean boost into our frame of reference, moving at a velocity v lab (t) with respect to the Galactic centre.
We take v lab = (11.1,235.2, 7.3) km/s + v ⊕ (t) [66] in the same galactocentric coordinate system as Fig. 1, where v ⊕ (t) is the Earth velocity, see e.g.[72,73].Following past literature [74][75][76][77][78][79], we take the stream velocity distribution to have the same Gaussian form, except we boost by v lab − v str to account for the stream's velocity, and set the width to σ str .Because the values of σ str are small, the stream signals are extremely narrowband in frequency.
We now relate this velocity distribution to the signal measured by a typical axion haloscope known as the lineshape.The axion field's oscillations are highly coherent, in keeping with its description as cold DM.Within a "coherence time", the axion field will appear to oscillate at a single frequency ω ≈ m a (1 + v 2 /2) where v ∼ 10 −3 c is some speed drawn from f (v).This frequency will then evolve on timescales longer than coherence time, depending on the spread in velocities: 01 ms (100µeV/m a ).If oscillations are measured over timescales much longer than this, then the discrete Fourier transform of that measurement will have a spectrum related to the distribution of component frequencies.For a measurement time, T int , the power spectrum S(ω) has a frequency resolution of ∆ω = 2π/T int .Therefore, if T int > 10 6 × 2π/m a (i.e.longer than a million oscillation periods) then we expect the axion's lineshape to be resolved [76,[80][81][82].
To illustrate this, we construct a simplified version of a signal power spectrum by discretising the distribution of frequencies into bins of width ∆ω [76,82,83].Figure 2 compares the smooth Maxwellian case (as in, for example, pre-inflationary axions) against the case where there are streams in the signal.The two will become strikingly different once the signal is integrated for timescales longer than 10 8 oscillation periods.
As well as showing up as sharp features in shortduration measurements, the lineshape will also evolve in amplitude over human timescales as we traverse the streams.Sampling over all of our streams and randomising our trajectory through them, we find that they would typically persist for around ∆t = 30 +23 −11 years (median and 68% containment).However, since we expect 10 2 -10 3 present in the lineshape at one time, the timescale over which the signal is expected to vary is on the order of weeks.Ultra-narrowband axion signals like these are already being searched for by haloscope collaborations [84][85][86], so in light of our results, we advocate that these efforts continue.
Conclusions.-We have evaluated the extent to which axion haloscopes are doomed to never discover the axion in the post-inflationary scenario because axions find themselves bound up inside of small substructures.By modelling the tidal disruption of these miniclusters by stars, we have found that the density of DM around us in the solar neighbourhood is refilled to around 80% of the 100-pc-scale average density value ρ DM ≈ 0.4 GeV cm −3 usually adopted by experiments.Combined with an estimate of the leftover density of DM in minivoids [38], this boosts the signal up to an impressive 90% of the commonly assumed value.In other words: axion haloscopes may not be doomed.The isolated minicluster mass spectrum is derived from simulations that had initial conditions calculated by a high-string-tension lattice simulation (high-κ, as in Ref. [45]).On the other hand, the high-mass extrapolated mass function for hierarchically merged miniclusters is inspired by N-body simulations evolving miniclusters up to z = 19.This latter regime, which has an initial power spectrum well-modelled as Gaussian white noise, extends up to ∼ 5 × 10 −7 M⊙.The total mass function for all miniclusters and minicluster halos (MCHs) is shown as a dashed line, while the coloured lines show the individual mass functions for minicluster halos comprised of different numbers of individual miniclusters.Our final results are heavily dominated by the upper end of this mass function comprised of the mergers of many miniclusters, however the dependence on details like the maximum mass and the exact power law turn out to not have a significant impact on our results.

Merged miniclusters
Much more important are the merged minicluster halos, forming in the white-noise regime, corresponding to fluctuations on scales k ≪ L −1 1 .Defining M 1 as the characteristic mass defined as 1 where m a is the axion mass, merged miniclusters will have masses M ≳ M 1 , for which the spectrum is well-described by Gaussian white-noise [30].Merged minicluster halos around and above the characteristic mass evolve towards a Navarro-Frank-White (NFW) [63] profile, which is given by with r s , ρ s scale radius and scale density, respectively.To model this population, we use numerical results of the white-noise N-body simulations from Ref. [47], which found a scale radius and concentration as a function of the minicluster mass and redshift, z.
From this, the radius at which the profile is truncated to contain the mass M can be found to be [47] which we adopt to model our minicluster, taking z = 19.Note that the overall scale for this radius is almost a factor of three larger than that assumed by Kavanagh et al. [42] in their minicluster disruption study.Their value is equivalent to a concentration parameter of only c = 100, which they argue would be a result of rapid initial truncation by tidal stripping as the miniclusters originally fell into the MW halo-this procedure we will also be implementing in our mass-loss pipeline.Nonetheless, this overall truncation scale of the profile ends up impacting some of our results.Instead of parameterising this uncertainty in terms of c, which does not directly enter our results, we instead parameterise the truncation radius by varying the prefactor in the expression above.We define this to be R 10 and will show results as a function of this quantity in Sec. .For our NFW halos, the individual mean densities are sampled from a normal distribution with mean and variance according to the N-body data [45].We sample masses in the range M ∈ [10 −12 , 5 × 10 −7 ] M ⊙ according to halo-mass function dn/d log M ∝ M γ with γ = −0.5.The upper mass cut-off has been estimated as ∼ 0.03⟨ρ⟩L 3  1 z 2 eq from the Press-Schechter estimate of the variance obtained from numerical simulations in [45].
Before moving on, it is worth remarking here briefly that our results are largely insensitive to the question of whether or not miniclusters host small solitonic cores called axion stars [62].Numerous recent studies have suggested that all gravitationally bound structures of a bosonic field ought to develop these cores [46,60,61], which may survive and even grow during mergers [59].Assuming the core-halo mass relation seen in simulations of the Schrodinger-Poisson system [58] is universal (i.e. down to minicluster masses [46]), we can evaluate it at the present day to yield the following relationship with the minicluster mass, Some intrinsic diversity away from this relation has been seen in some simulations, but we are primarily interested in the ratio between the core and halo mass at an order-of-magnitude level [55][56][57].
Axion stars have presented a potentially problematic theoretical uncertainty in previous studies that focused on the survival of miniclusters.Since the axion star could potentially prevent miniclusters from being totally stripped.However, we notice here that the total mass contained in the axion star is vastly subdominant, and all of our results depend dominantly on the loosely bound outer layers of the high-mass merged miniclusters which are stripped early on in the evolution.Hence we do not need to worry about whether or not miniclusters host axion stars, and accounting for presence in the cores of our miniclusters, or implementing cuts on stripped miniclusters which would host unphysically large axion stars (as in Ref. [42]), would not change any of our results.

GALACTIC STELLAR DENSITY MODEL
The number of stars encountered within an impact parameter of b max is given by the integral along the minicluster trajectory: To get the stellar number density, n ⋆ (t) we write down a model for the galactic stellar mass density and choose a characteristic star mass M ⋆ = 1 M ⊙ .Because the disruption probability is proportional to the number density of stars times their mass, including a detailed stellar mass function is not expected to change our conclusions if we keep the volume and total mass fixed [41]-and in any case, 1 M ⊙ is a very typical star mass for standard stellar initial mass functions.Our stellar density model, as in Refs.[41,42], consists of a bulge and disk, of which the latter has two components, labelled as thick and thin: The bulge density profile is modelled as a truncated power-law [42,53,54] where r = x 2 + y 2 is the cylindrical radius in the disk-plane, and z the height above the disk plane.The central density is ρ bulge,0 ≃ 99.5 M ⊙ /pc 3 , and the radial scales are r 0 = 0.075 kpc and r cut = 2.1 kpc.Both the thick and the thin stellar disk profiles are modelled by a double exponential [52] where r t,T , z t,T are the scale radius and height, and Σ t,T ⋆ are the stellar surface densities: Our results are not sensitive to the detailed parameterisation of the stellar density, and varying any of these numbers within up to a few tens of percent leaves our results generally unaffected.

MINICLUSTER DISRUPTION COMPUTATION
Here we provide a brief account of the disruption computation, though for further discussion of the motivation and justification behind this treatment, we refer to Ref. [42].
Once an encounter time (t enc ), star mass (M ⋆ ), and impact parameter (b) have been chosen, the first step is to compute the energy injected into the minicluster by that encounter.Under the distant tide and the impulse approximations assumptions, the energy injection can be written as [48]: where ⟨R 2 ⟩ is the minicluster mean-squared radius.This can be expressed in terms of a structural parameter, α, For an NFW minicluster with concentration c = 100, Refs.[42,44] find α 2 = 0.13.For the power-law profile, we use α 2 = 0.27 [42].The energy injection must then be compared to the binding energy of the minicluster, E b , which is given by the familiar ∼ GM 2 /R up to some factor β, that encodes details of the internal structure.We define it via, where the structural parameter is (see, e.g., Refs.[42,65].) where M enc (r) is the mass enclosed inside radius r.If the energy injected into the minicluster by the star's tidal field exceeds the minicluster's binding energy, we expect the minicluster to lose and order-1 fraction of its mass [44].However, since most encounters inject ∆E ≪ E b , we deal with perturbations giving rise to small mass losses ∆M .
In the distant tide approximation, the normalised injected energy increases proportionally to R 2 : where E = −E/m a is the energy per unit mass.The mass that is then left unbound to the cluster following this injection is then given by the total mass in axions with energies smaller than ∆E [42] ∆M In our computation, we use a numerical interpolation of this mass loss function, ∆M (∆E), given in Ref. [42]. 2 The functional dependence of the mass loss is in rough agreement with Ref. [44], who instead performed idealised N-body simulations of an NFW halo interacting with a 1 M ⊙ star.We notice mild disagreement only for ∆E/E b ≲ 10 −2 and both models approach ∆M/M ∼ 1 as expected for very large energy injections ∆E/E b ≫ 1 which can totally disrupt the minicluster, though these occur rarely.Reference [44] limits the treatment to NFW initial profiles only, while the modelling of Ref. [42] can be extended to different profiles.The initial energy of the minicluster before the energy injection is given by the sum of the kinetic energy (1/2M σ 2 mc ) and binding energy (E b ): where κ = 1.73, 3.54-for isolated and merged miniclusters respectively-parameterises the velocity dispersion as follows, Following energy conservation, the final minicluster energy after ∆E has been injected can be written as, f ej (∆E) is the fraction of the imparted energy that becomes unbound and E unb i (∆E) the energy originally possessed by those axions that became unbound.Bound axions are those satisfying the integral cutoff in Eq.(S20), and unbound axions the reverse.Like with the function ∆M (∆E) we adopt interpolations of f ej (∆E) and E unb i (∆E) provided by Ref. [42].
After each encounter, and before moving to the next one, the minicluster mass and radius are updated to, We continue to perturb each minicluster through its N enc encounters until either it is fully unbound (E f > 0) or we have reached the end of the orbit.Because many of the encounters occur very rapidly, we group together encounters in time bins of ∆t = t MW /10 4 , which is chosen to be shorter than the relaxation timescale t rel ∼ (Gρ mc ) −1/2 ∼ O(10 Myr).This makes the procedure more computationally efficient whilst also ensuring that the relaxation to a new profile is only allowed to occur whenever there is sufficient time for it to occur (usually this is between disk crossings).In addition to these successive mass-loss events, we also implement a second mass-loss procedure running in parallel to the stellar disruption in which the minicluster is continually stripped by the tidal field of the MW halo.The relevant formula for the rate of mass loss through this process is, where A = 1.34 and ζ = 0.07 [50], t dyn = 2.4 Gyr is the dynamical time of the MW, and M MW = 10 12 M ⊙ .The prescription for this is discussed in Appendix A of Ref. [42] (see also Ref. [51]).Stripping by the MW halo contributes towards some of the early mass-loss experienced by the most massive miniclusters in our simulation, but is suppressed heavily by (M/M MW ) ζ for smaller miniclusters.This also means it becomes increasingly negligible towards the end of the simulation.In Fig. S2 we show the stellar density as a function of time, ρ ⋆ (t) for a random orbit ending at the Solar position today t = t MW (we zoom in on the initial 5 Gyr for visual clarity, although in our calculation we evolve each ρ DM = N mc M/V .The expectation value for the density that we measure on Earth is then just the sum of all of the individual stream densities: ⟨ρ obs ⟩ = ⟨N str ⟩ρ str = N str M/V str .Notice that if we substitute in ⟨N str ⟩ from above we get back ⟨ρ obs ⟩ = ρ DM .We enforce this because when averaging the DM densities observed at many different points in a volume V ≫ V str should give us back DM density defined on that scale. So the possible values of ρ obs follow a binomial distribution with a mean of ρ DM .When ⟨N str ⟩ is large, the central limit theorem dictates that the distribution of ρ obs is a Gaussian centred on ρ DM .However in the opposite extreme, when ⟨N str ⟩ ≲ 1 (for example if the miniclusters were not disrupted for whatever reason) the mean is still ρ DM , but only because the high probability to observe ρ obs = 0 is balanced against the very rare probability to observe a hugely enhanced density ρ obs = ρ str ≫ ρ DM because the stream volumes would have to be very small to cause N str ≲ 1.
So for this reason, the pertinent result is not the expectation value of the density, because as a result of conserving the total mass of DM on large scales, this is always ρ DM by construction 3 .Instead, we must consider the distribution of ρ obs , which is dictated by N str .
So to see which regime we are in, whilst also evaluating the effects of changing some of our assumed input parameters, let us develop a simple analytic estimate for the expected number of observed streams.Starting with the formula for the stream volume, parametrically this can be modelled as, where 4σ mc (t MW − t enc ) is roughly the stream length assuming a disruption happens (as we observe here) very quickly, i.e. after only a few disk crossings.We use a ≳ sign because the minicluster velocity dispersion gives us a lower limit on the stream's velocity dispersion.The former is expressed as σ 2 mc = κGM/R, as in Eq.(S22) above.If we now plug in our formula for the radii of merged miniclusters R(M ) = R 10 (M/10 −10 M ⊙ ) 1/3 (as in Eq.(S7)), we see that the dependence on M in the stream number N str = ρ DM V str /M drops out entirely, giving us, where we assign ρ mc = ρ DM (1 − f void ) ≈ 0.01 M ⊙ pc −3 as the portion of the DM density that was originally bound inside miniclusters (i.e.not in the minivoids).To evaluate this expression, we have the average mass-lossweighted encounter time from our simulations of t enc ∼ 3 Gyr, which gives us a number strikingly close to our fiducial numerical result.So referring back to our statements above about the distribution of ρ obs , the key point to highlight is that we are quite firmly in the regime N str ≫ 1, especially because we have fixed the stream length σ mc (t MW − t enc ) to be the smallest that we expect it to grow to.We may account for additional heating of the axions to dispersions larger than σ mc , but will only act to increase this number, further reinforcing the fact that we expect there to be a large number of streams overlapping at our position.
This estimate also reveals a general insensitivity of the stream number to many model parameters.So we expect N str ≫ 1 to remain the case if we allow the mass function, density profiles, mass-radius relation, and abundance of merged miniclusters to vary over any reasonable range.We can see this from our simple analytic estimate, but have confirmed this remains the case in our full simulation as well, which accounts for variable disruption probabilities as a function of M , varying disruption times as a function of the minicluster orbit, as well as the varying densities along the stream.In fact, the reason our numerical analysis leads to a median ρ obs slightly smaller than ρ DM is because of these additional details.Nevertheless, our full numerical results are also generally insensitive to most of our more crude assumptions.To show a few explicit examples of factors that cause the result to change the most, Fig. S3 shows multiple versions of our final result-the distribution in the sum of the stream densities, Σρ str -for multiple different choices for the initial sizes of the miniclusters R 10 (or equivalently, their concentration parameters, c), as well as for different stream heating factors beyond the most conservative case where σ str = σ mc .We see that even with these alternative choices, the median value of the density is mostly unchanged, and these quantities instead only affect the variance in the distribution (as expected because these quantities affect ⟨N str ⟩).
FIG. S3.Independence of our final result on the minicluster mass-radius relation, parameterised by the radius scale R10, and how much the stream's velocity dispersion is additionally heated from its initial value given by the minicluster's velocity dispersion.In all cases we see that the median value of observed density is roughly the same, with the primary role of these parameters being to change the expected number of overlapping streams, and hence the variance in the density, as opposed to its expectation value.The number of streams observed in the three examples ranges from Nstr = 63, 246, and 547, from the top to bottom, and a factor of 10 larger for the σstr = 10σmc case.

IMPACT OF MINICLUSTER HALO SUBSTRUCTURE
As argued in the previous section, the conclusion that we expect the local DM density to be almost saturated after minicluster disruption can be considered reasonably robust, assuming our result that the merged miniclusters lose the majority of their mass holds up.However, it is this final issue that is really the only factor that could change our result in any substantial way.If the merged miniclusters are in fact more stable than the NFW halos that we have modelled them as, then it is possible they hold on to more of their mass and not spread it over such large volumes.One way in which this might not be captured in our treatment is if the merged miniclusters are not smooth NFW profiles but rather retain a degree of substructure from the isolated miniclusters that created them.We identify this as the primary remaining theoretical uncertainty that should be resolved with high-resolution (i.e.longer timescale) simulations that can evolve the minicluster halos to lower redshifts.
So are our simulated minicluster halos actually just clusters of compact miniclusters?Minicluster halos are partially made by merging and accretion of isolated MCs.It is crucial for what follows to know the fraction of the mass of the MCH that is smoothly distributed and the fraction that it is still retained in compact MCs.This can be extracted from the numerical simulations to some extent using the subfind algorithm implemented in the gadget-4 [49] code.From our recent N-body simulations (see Sec. ) we calculate the fraction of axions contained in subhalos for each host (i.e., minicluster halo) identified by the friend-of-friends (FOF) algorithm.In this calculation we however do not consider the main virialised halo (the first subhalo within the FOF group), since the latter will not contain itself additional substructure.As seen in Fig. S4, we find that only ∼ 10% of the MCHs are composed of smaller virialised orbiting miniclusters.As expected, we confirm that high-mass MCHs have comparatively more substructure than low-mass MCHs.This lends further confidence that our modelling of these larger miniclusters as smooth NFW profiles can be considered reasonably accurate.

-FIG. 1 .
FIG.1.Schematic of our study.Left: We begin by modelling the orbits of many miniclusters all ending at the Solar position today, (X, Y, Z) ≈ (8, 0, 0) kpc.For each orbit, we draw random encounters with stars from a MW stellar density model that includes the thin and thick disks, and the bulge.We colour the minicluster orbits by the integral of this density along the orbit.For clarity, we display only the last 50 Myr of the orbit.Centre: Zooming in on a 100-pc-radius sphere around the Sun, we illustrate how the axions stripped from their host minicluster are elongated into tidal streams of O(pc) length.Right: Zoom in on ∼mpc scales relevant for experiments where the network of tidal streams sums to give the local density in axions.
FIG. S1.Halo mass functions resulting from N-body simulations.The isolated minicluster mass spectrum is derived from simulations that had initial conditions calculated by a high-string-tension lattice simulation (high-κ, as in Ref.[45]).On the other hand, the high-mass extrapolated mass function for hierarchically merged miniclusters is inspired by N-body simulations evolving miniclusters up to z = 19.This latter regime, which has an initial power spectrum well-modelled as Gaussian white noise, extends up to ∼ 5 × 10 −7 M⊙.The total mass function for all miniclusters and minicluster halos (MCHs) is shown as a dashed line, while the coloured lines show the individual mass functions for minicluster halos comprised of different numbers of individual miniclusters.Our final results are heavily dominated by the upper end of this mass function comprised of the mergers of many miniclusters, however the dependence on details like the maximum mass and the exact power law turn out to not have a significant impact on our results.
FIG.2.An illustration of how minicluster streams would manifest in the experimental signature of axions called the "lineshape".The frequency resolution is inversely proportional to the duration of the observation, Tint, from which this spectrum is obtained via a discrete Fourier transform.If Tint is short, then the lineshape is indistinguishable from a smooth halo.However, samples longer than ∼ 10 8 oscillation periods have sufficient resolution to identify the streams.