Axion minicluster streams in the Solar neighbourhood

In this section, we give further details on the simulations that inspire our initial population of miniclusters and detail the structural parameters that are necessary to model their subsequent disruption. As described in Refs. Eggemeier:2022hqa ; Pierobon:2023ozb we evolve realistic configurations of the axion DM field within the redshifts $10^{16}\lesssim z\lesssim 10^{2}$.

The distribution of miniclusters can either be predicted using the spectrum of density fluctuations $\Delta^{2}$ left behind following the collapse of the axion string-wall network, or they can obtained directly by performing an N-body simulation. Recently, Ref. Pierobon:2023ozb demonstrated how both the mass distribution and internal structures of miniclusters vary when adopting different treatments of the early Universe dynamics—in particular the modelling of axion strings. Due to these outstanding uncertainties, current simulations Eggemeier:2019jsu ; Xiao:2021nkb ; Pierobon:2023ozb are unfortunately unable to conclusively identify a universal density profile for all forming miniclusters. For this study, we have chosen to adopt the early Universe modelling of direct simulations at large string tension (dubbed high-$\kappa$) from Ref. Pierobon:2023ozb , which best describes the string-wall network with respect to previous works.

Uncertainties aside, in more general terms however, recent N-body simulations are able to show that the NFW profile fits well the high-mass end of the mass function corresponding to hierarchical mergers of smaller miniclusters. Whereas on the low-mass end, consisting of halos that remain isolated after their initial collapse, power-law profiles fit much better. These two types of minicluster can be seen in Fig. S1 where we display the halo mass function taken from the N-body simulations presented in Ref. Pierobon:2023ozb . The non-linear regime refers to those miniclusters that promptly collapse and then do not undergo significant mergers, whereas the high mass tail, which grows out of the hierarchical merger of many initially isolated miniclusters, has a spectrum well-described by Gaussian white noise, hence why it is referred to the white-noise regime. We describe how we model each of these now in turn. The high-mass extrapolation is inspired by the N-body simulations, which follow the miniclusters to even higher redshifts of $z=19$, when they fall into the first galactic halos.

The number of stars encountered within an impact parameter of $b_{\rm max}$ is given by the integral along the minicluster trajectory:

To get the stellar number density, $n_{\star}(t)$ we write down a model for the galactic stellar mass density and choose a characteristic star mass $M_{\star}=1~{}M_{\odot}$. 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 Dokuchaev –and in any case, $1~{}M_{\odot}$ is a very typical star mass for standard stellar initial mass functions.

Our stellar density model, as in Refs. Kavanagh:2020gcy ; Dokuchaev , 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 Binney:1996sv ; Bissantz:2001wx ; Kavanagh:2020gcy

where $r=\sqrt{x^{2}+y^{2}}$ is the cylindrical radius in the disk-plane, and $z$ the height above the disk plane. The central density is $\rho_{{\rm bulge},0}\simeq 99.5~{}M_{\odot}/{\rm pc}^{3}$, and the radial scales are $r_{0}=0.075$ kpc and $r_{\rm cut}=2.1$ kpc.

Both the thick and the thin stellar disk profiles are modelled by a double exponential Dehnen:1996fa

where $r^{t,T},z^{t,T}$ are the scale radius and height, and $\Sigma^{t,T}_{\star}$ 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.

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. Kavanagh:2020gcy .

Once an encounter time ($t_{\rm enc}$), star mass ($M_{\star}$), 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 1958ApJ…127…17S :

where $\langle R^{2}\rangle$ is the minicluster mean-squared radius. This can be expressed in terms of a structural parameter, $\alpha$,

For an NFW minicluster with concentration $c=100$, Refs. Kavanagh:2020gcy ; Shen:2022ltx find $\alpha^{2}=0.13$. For the power-law profile, we use $\alpha^{2}=0.27$ Kavanagh:2020gcy .

The energy injection must then be compared to the binding energy of the minicluster, $E_{b}$, which is given by the familiar $\sim GM^{2}/R$ up to some factor $\beta$, that encodes details of the internal structure. We define it via,

where the structural parameter is (see, e.g., Refs. Kavanagh:2020gcy ; Dandoy:2022prp .)

where $M_{\rm 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 Shen:2022ltx . However, since most encounters inject $\Delta E\ll E_{b}$, we deal with perturbations giving rise to small mass losses $\Delta M$.

In the distant tide approximation, the normalised injected energy increases proportionally to $R^{2}$:

where $\mathcal{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 $\Delta\mathcal{E}$ Kavanagh:2020gcy

In our computation, we use a numerical interpolation of this mass loss function, $\Delta M(\Delta E)$, given in Ref. Kavanagh:2020gcy .²²2Data available at github.com/bradkav/axion-miniclusters. The functional dependence of the mass loss is in rough agreement with Ref. Shen:2022ltx , who instead performed idealised N-body simulations of an NFW halo interacting with a $1~{}M_{\odot}$ star. We notice mild disagreement only for $\Delta E/E_{b}\lesssim 10^{-2}$ and both models approach $\Delta M/M\sim 1$ as expected for very large energy injections $\Delta E/E_{b}\gg 1$ which can totally disrupt the minicluster, though these occur rarely. Reference Shen:2022ltx limits the treatment to NFW initial profiles only, while the modelling of Ref. Kavanagh:2020gcy 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\sigma_{\rm mc}^{2}$) and binding energy ($E_{b}$):

where $\kappa=1.73,\,3.54$—for isolated and merged miniclusters respectively—parameterises the velocity dispersion as follows,

Following energy conservation, the final minicluster energy after $\Delta E$ has been injected can be written as,

$f_{\rm ej}(\Delta E)$ is the fraction of the imparted energy that becomes unbound and $E^{\rm unb}_{i}(\Delta 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 $\Delta M(\Delta E)$ we adopt interpolations of $f_{\rm ej}(\Delta E)$ and $E^{\rm unb}_{i}(\Delta E)$ provided by Ref. Kavanagh:2020gcy .

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_{\rm 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 $\Delta t=t_{\rm MW}/10^{4}$, which is chosen to be shorter than the relaxation timescale $t_{\rm rel}\sim(G\rho_{\rm mc})^{-1/2}\sim\mathcal{O}(10\,{\rm 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 $\mathcal{A}=1.34$ and $\zeta=0.07$ Jiang:2014nsa , $t_{\rm dyn}=2.4$ Gyr is the dynamical time of the MW, and $M_{\rm MW}=10^{12}~{}M_{\odot}$. The prescription for this is discussed in Appendix A of Ref. Kavanagh:2020gcy (see also Ref. vandenBosch:2004zs ). 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_{\rm MW})^{\zeta}$ 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, $\rho_{\star}(t)$ for a random orbit ending at the Solar position today $t=t_{\rm MW}$ (we zoom in on the initial 5 Gyr for visual clarity, although in our calculation we evolve each minicluster until today). The orbit has a mean galactocentric radius of $\sqrt{X^{2}+Y^{2}}=22$ kpc and a mean inclination of $\sin^{-1}(Z/\sqrt{X^{2}+Y^{2}})=21^{\circ}$. In the upper panel, we show the gradual loss in mass as a function of time normalised to the minicluster’s initial mass. The miniclusters are randomly sampled in mass from their mass function (see S.M. I), but all follow the same orbit, an $X$-$Y$ projection of which is shown in the inset panel to the right. The mass loss exhibits a step-like behaviour for this orbit due to the exponential enhancement in the stellar density whenever the orbit crosses the disk plane $Z=0$. We see here clearly that the merged NFW-profile miniclusters lose a much larger fraction of their mass and over a much shorter timescale—particularly in the first disk-crossing. In contrast, many isolated power-law miniclusters survive with an $\mathcal{O}(1)$ fraction of their mass today.

To provide further intuition behind the results presented in the main text, in this section, we make some back-of-the-envelope estimates that explain why we should expect the DM density to be mostly replenished by the disruption of the most massive miniclusters. Doing this will also allow us to evaluate the possible effects (if any) that changing our modelling assumptions might have on our main result.

Let us start with a simplified description of the situation wherein the DM is composed entirely of merged miniclusters, which we assume are all disrupted in their entirety. This is already very close to our fiducial analysis since most of the mass is in the high-mass merged miniclusters anyway and, as we have seen (e.g. in Fig. S2), the stellar disruption process causes them to lose the majority of their mass quite quickly. If we now further simplify the situation and assume all the resulting streams have the same volume $V_{\rm str}$ and the same mass $M$ then the expected number of streams overlapping at a single point is,

where in the second step we have rewritten the local galactic DM density we are required to saturate as $\rho_{\rm DM}=N_{\rm 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: $\langle\rho_{\rm obs}\rangle=\langle N_{\rm str}\rangle\rho_{\rm str}=N_{\rm str}M/V_{\rm str}$. Notice that if we substitute in $\langle N_{\rm str}\rangle$ from above we get back $\langle\rho_{\rm obs}\rangle=\rho_{\rm DM}$. We enforce this because when averaging the DM densities observed at many different points in a volume $V\gg V_{\rm str}$ should give us back DM density defined on that scale.

So the possible values of $\rho_{\rm obs}$ follow a binomial distribution with a mean of $\rho_{\rm DM}$. When $\langle N_{\rm str}\rangle$ is large, the central limit theorem dictates that the distribution of $\rho_{\rm obs}$ is a Gaussian centred on $\rho_{\rm DM}$. However in the opposite extreme, when $\langle N_{\rm str}\rangle\lesssim 1$ (for example if the miniclusters were not disrupted for whatever reason) the mean is still $\rho_{\rm DM}$, but only because the high probability to observe $\rho_{\rm obs}=0$ is balanced against the very rare probability to observe a hugely enhanced density $\rho_{\rm obs}=\rho_{\rm str}\gg\rho_{\rm DM}$ because the stream volumes would have to be very small to cause $N_{\rm str}\lesssim 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 $\rho_{\rm DM}$ by construction³³3that is, in this simple example where our streams are objects of the same size and have constant internal densities.. Instead, we must consider the distribution of $\rho_{\rm obs}$, which is dictated by $N_{\rm 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\sigma_{\rm mc}(t_{\rm MW}-t_{\rm 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 $\gtrsim$ sign because the minicluster velocity dispersion gives us a lower limit on the stream’s velocity dispersion. The former is expressed as $\sigma_{\rm mc}^{2}=\kappa 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_{\odot})^{1/3}$ (as in Eq.(S7)), we see that the dependence on $M$ in the stream number $N_{\rm str}=\rho_{\rm DM}V_{\rm str}/M$ drops out entirely, giving us,

where we assign $\rho_{\rm mc}=\rho_{\rm DM}(1-f_{\rm void})\approx 0.01~{}M_{\odot}\,{\rm 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 taken the average mass-loss-weighted encounter time from our simulations of $t_{\rm enc}\sim 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 $\rho_{\rm obs}$, the key point to highlight is that we are quite firmly in the regime $N_{\rm str}\gg 1$, especially because we have fixed the stream length $\sigma_{\rm mc}(t_{\rm MW}-t_{\rm 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 $\sigma_{\rm 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_{\rm str}\gg 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 $\rho_{\rm obs}$ slightly smaller than $\rho_{\rm 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, $\Sigma\rho_{\rm 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 $\sigma_{\rm str}=\sigma_{\rm 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 $\langle N_{\rm str}\rangle$).

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 Springel:2020plp code. From our recent N-body simulations (see Sec. I) 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 $\sim 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.