Stellar Disruption of Axion Miniclusters in the Milky Way

Axion miniclusters are dense bound structures of dark matter axions that are predicted to form in the post-inflationary Peccei-Quinn symmetry breaking scenario. Although dense, miniclusters can easily be perturbed or even become unbound by interactions with baryonic objects such as stars. Here, we characterize the spatial distribution and properties of miniclusters in the Milky Way (MW) today after undergoing these stellar interactions throughout their lifetime. We do this by performing a suite of Monte Carlo simulations which track the miniclusters' structure and, in particular, accounts for partial disruption and mass loss through successive interactions. We consider two density profiles - Navarro-Frenk-White (NFW) and Power-law (PL) - for the individual miniclusters in order to bracket the uncertainties on the minicluster population today due to their uncertain formation history. For our fiducial analysis at the Solar position, we find a survival probability of 99% for miniclusters with PL profiles and 46% for those with NFW profiles. Our work extends previous estimates of this local survival probability to the entire MW. We find that towards the Galactic center, the survival probabilities drop drastically. Although we present results for a particular initial halo mass function, our simulations can be easily recast to different models using the provided data and code. Finally, we comment on the impact of our results on lensing, direct, and indirect detection.


I. INTRODUCTION
Both the dark matter (DM) and Strong-CP problems can be solved by introducing a new global symmetry into the Standard Model (SM) of particle physics [1,2].This new Peccei-Quinn (PQ) symmetry U PQ (1) predicts a hypothetical particle known as the QCD axion [3,4].Unlike weakly interacting massive particles, DM axions are typically much lighter than the rest of the SM [5][6][7].Their production mechanism must therefore rely upon non-thermal processes to ensure they are non-relativistic at the time of recombination.These non-thermal processes generically produce gravitationally bound clumps of axions known as axion miniclusters.In this paper, we characterize the degree to which tidal interactions can change the properties of these miniclusters over the lifetime of the Milky Way (MW).
The production of QCD axion DM is tightly connected to the thermal history of the Universe.After the U PQ (1) symmetry is spontaneously broken, the axion field relaxes towards the bottom of its potential.When the Universe has cooled down to the QCD phase transition, non-perturbative QCD instantons lead to the explicit breaking of the PQ symmetry [44,45], giving rise to a new CP-conserving minimum in the potential.After the QCD phase transition, the axion field undergoes coherent oscillations about this minimum, damped by the Hubble expansion rate H.This process is known as the vacuum realignment mechanism [5][6][7].The DM energy density in this scenario is stored in the coherent oscillations of the axion condensate and depends on the initial value of the axion field when the PQ symmetry breaks, which is parametrized by the initial misalignment angle θ i .We generally expect the axion energy density to be proportional to θ2 i except around θ i ∼ π where the non-harmonic terms in the axion potential become important [46][47][48][49][50].The properties of the axion condensate crucially depend on whether the spontaneous breaking of the U PQ (1) symmetry occurs before or after the end of inflation.In the pre-inflationary scenario, the value of θ i is uniquely selected over the whole observable Universe.
Here, we consider the opposite scenario in which the PQ symmetry is broken after the end of inflationthe post-inflationary scenario.In this case, the initial misalignment angle θ i takes different values in different patches of the observable universe, since no patch has been selected by the inflationary process.In this post-inflationary scenario, self-gravitating substructures called axion miniclusters (AMCs) [51][52][53] are expected to form.Moderate O (1) overdensities initially lead to the formation of minicluster 'seeds' [54] which later collapse into gravitationally bound AMCs at around matterradiation equality [51][52][53]55]. 2 Instead, AMCs cannot form in the pre-inflationary scenario even when the initial conditions of the QCD axion field are extremely finetuned [56].
Significant progress has been made towards solving the early universe dynamics of the axion; numerical simulations loosely constrain the fraction of cold DM axions in these bound structures f AMC to be O(1%−100%) [54,57].This fraction f AMC plays a fundamental role in the prospects for axion DM detection.For example, if most of the DM is bound in AMCs, the probability of having a substantial DM density near Earth may drop drastically, due to the rarity of Earth-AMC encounters [58], making direct detection methods ineffective.Similarly, the encounter rate of AMCs with a single NS is likely to be low, rendering radio observations of individual NSs ineffective in the search for axion-photon conversion.Fortunately, encounters between miniclusters and the MW population of NSs can give rise to interesting transient radio signals, as we show in our companion paper Ref. [59].
It is possible to assess the fraction of cold axions bound in AMCs through femto-lensing induced by individual miniclusters [60,61], micro-lensing from minicluster halos formed after matter-radiation equality from hierarchical merging [62][63][64], and minicluster lensing of highly magnified stars [65].These studies typically treat f AMC as a constant, but miniclusters interacting with their environment in fact cause f AMC to become dependent on both time and spatial position.Tidal interactions of miniclusters with larger host halos, with each other, and with condensed baryonic objects all play a pivotal role in the survival of miniclusters [66].In this paper, we quantify the degree to which interactions between miniclusters and stars can change the characteristics of AMCs today.We focus on the MW, where stars are abundant and constitute the dominant disruption mechanism [67,68].In particular, we present a formalism that describes the reaction of an AMC's internal state to an interaction with a star.We then use this formalism to run Monte Carlo simulations of AMCs as they orbit the MW and interact with the stellar population.For computational simplicity, we make a number of simplifying assumptions.Firstly, we do not concurrently evolve AMCs through structure formation and stellar disruption.Secondly, we assume that the MW has been in a steady state since its formation.Despite these assumptions our results represent a fundamental step towards quantifying the importance of tidal stellar interactions for the distribution of AMCs in the MW.As we argue in § IV A, relaxing our assumptions will not change our overall conclusions but future work should quantitatively address these issues.
In addition to AMCs, axion stars (ASs) -another class of axionic astrophysical object -are expected to form and remain stable over cosmological times (in particular the dilute branch of axion stars [69]).Importantly, they may readily form within miniclusters, producing a central core [70][71][72].For simplicity, we do not simultaneously consider both AMCs and ASs.Instead, we make a cut on the minicluster parameter space in order to focus on those AMCs for which the density profile is known most reliably (as described in § II E).
This paper is structured as follows: in Sec.II we describe the initial distributions of AMCs in the MW that represent the starting point of our work.Section III discusses the dynamics of successive stellar encounters.In Sec.IV we discuss our Monte Carlo simulations and how to interpret the results.We then discuss the minicluster population today in Sec.V and how this applies to observational results in Sec.VI.Finally, we discuss future work and conclude in Sec.VII.Throughout this paper, we limit our discussion to axions which constitute 100% of the DM and only consider the QCD axion.Specifically, we assume a KSVZ-like axion 3 of mass m a = 20 µeV, 4 and we comment on results for different masses in Sec.VII.Our results can be extended to other axion models, provided that the distributions of AMC masses and overdensities are modified accordingly.All code associated with this work (and the companion paper, Ref. [59]) is available online at github.com/bradkav/axion-miniclusters[76].

II. MINICLUSTERS IN THE MILKY WAY
Cold axions are produced in the early Universe through non-thermal processes at the time t osc at which the axion field begins oscillating, when the Hubble rate is of the order of the axion mass.The present number density of axions is given by where a is the scale factor (set to unity today), and ρ is the axion energy density.The subscript 'osc' indicates the value of the parameter at t osc .Equating the energy density of axions today to the observed DM abundance fixes the value of the axion mass.The computation of ρ osc requires one to simulate the dynamics of the topological defects associated with the spontaneous breaking of the PQ symmetry, which is a significant challenge.Here, we follow recent literature on the subject and we fix m a = 20 µeV [57,75], although a wide range of masses is still possible [77].
An individual AMC can be characterized by an initial overdensity parameter δ, discussed in § II A, and an initial mass described in § II B. Due to the randomness of the initial conditions of the axion field over causallyconnected patches, AMCs are formed with a range of overdensity parameters and masses.We note that predictions for the AMC properties are still under debate in the literature and we attempt to highlight these uncertainties throughout.We emphasize, however, that our framework can be straightforwardly re-cast under different assumptions for the initial distribution of AMCs, as we discuss in Sec.V.

A. Distribution of Overdensities
After the DM axion field has started to oscillate, its mean background density scales with temperature as ρa (T ) = 3T eq s(T )/4, where s(T ) is the entropy density at temperature T and T eq is the temperature at matterradiation equality.A density fluctuation ρ a > ρa decouples from the cosmological expansion at the temperature T δ = (1+δ) T eq , where δ ≡ (ρ a − ρa )/ρ a is the overdensity parameter.After this, the overdense region undergoes non-linear gravitational collapse, becoming gravitationally bound into an AMC [52,60].Assuming spherical collapse, the density of a virialized minicluster is [78] where ρ eq is the average matter density at matterradiation equality.The AMC density ρ AMC (δ) does not depend on the interaction of the axion with matter, nor on the axion mass.For this reason, we expect that our results will be unchanged for AMCs formed from an axionlike field, as long as the axion-like field makes up the entirety of the DM.The distribution of overdensities df AMC /dδ can be assessed through numerical simulations [57,78].In this paper, we adopt the expression for df AMC /dδ used in Ref. [57] and we span the range of values δ ∈ [0.1, 20] (corresponding to characteristic densities ρ AMC ∈ 10 2 , 2 × 10 10 M pc −3 ).For completeness, we report the formula used in Appendix C. Mapping out the correlation between δ and the AMC mass M AMC is currently challenging, as the simulations used to derive df AMC /dδ stop at matter-radiation equality [57], while we are interested in the mass function in the late Universe.In the following, we assume that there is no correlation between df AMC /dδ and the AMC mass distribution introduced below, though our formalism can be straightforwardly extended to incorporate such correlations.

B. Initial Halo Mass Function
The characteristic comoving number density of AMCs per logarithmic mass interval is described by the halo mass function (HMF).The HMF at matter-radiation equality (at redshift z eq ) can be assessed by evolving the PQ field from the moment at which the PQ symmetry breaks until z eq [54,57].The high-end tail of the HMF shows an exponential cutoff which, at recombination, is placed at around the largest mass of the AMCs, M max (z eq ) ≈ M 0 [54,79,80].The characteristic mass M 0 is associated with the axion energy density contained within a Hubble horizon at t osc [81] where M is the Solar mass.
Perturbations in the axion density continue to grow after matter-radiation equality, so that the HMF evolves under hierarchical structure formation.N -body simulations following AMC structures from recombination to z ≈ 99 lead to an HMF dP/d ln M AMC ∝ M AMC γ , with a characteristic slope γ ∼ −0.7 [79].This result corroborates the semi-analytic solution obtained by using the Press-Schechter formalism [82], which finds that the mass function at late times scales as M −0.68 for small masses, and as M −0.35 for large masses, over the mass interval 10 −15 M/M 10 −9 [64].Earlier work found the slope ∼ −0.5 [62,63].
As structure formation proceeds, the high-end cutoff of the HMF evolves according to the Press-Schechter analysis, since AMCs of mass M AMC > M 0 form through hierarchical structure formation from the early seeds of mass M AMC ≤ M 0 .On the other hand, the assessment of the low-end mass cutoff is challenging [83,Sec. 5.2].Both the semi-analytic formalism and the numerical assessment of the low-end tail of the mass distribution M AMC M 0 show limitations due to a number of factors.For example, fluctuations in this regime are not Gaussian so the Press-Schechter formalism cannot be used; in addition, resolving the power spectrum at such small scales is a numerical challenge [84].Note that neither the numerical simulations of the early Universe, nor the Press-Schechter formalism account for the possible presence of ASs, which we discuss in § II E. (5) Note that the lower end of this mass range will be suppressed by our AS cut, as described in § II E. The characteristic radius R AMC for an AMC of mass M AMC is of the order of The HMF in Eq. ( 4) has been obtained without considering the effects of the tidal stripping of AMCs due to nearby stars in the MW or due to the mean Galactic field.The distribution in Eq. ( 4) gives us the initial HMF.In Sec.III, we assess the effects of tidal stripping on the population of AMCs, which effectively modifies the initial HMF to yield the true HMF today.Because the HMF is a falling Power-law, the lightest AMCs will dominate the MW population, meaning that our results could in principle be sensitive to the low-mass cut-off M min .While M min lies below the minimum mass that passes our AS cut criteria, it may still affect the AMC population through its influence on the normalization of the HMF.9), and NFW, Eq. ( 8).Vertical dashed lines show the truncation radii RAMC.We fix the characteristic mass and density to MAMC = 10 −10 M and ρAMC = 106 M pc −3 (δ ≈ 1.55) respectively.In the NFW case, the mean density is much lower and the AMC is much larger.

C. Distribution of AMCs in the Galaxy
Given the DM density profile in the MW ρ DM (r), we model the spatial distribution of the number density of AMCs as where M AMC is the mean AMC mass before disruption is accounted for.Using the HMF in Eq. ( 4), we obtain the value M AMC = 1.4 × 10 −14 M .Here, we set f AMC = 100% (though all our results can be trivially rescaled).
The distribution of the DM density in the Galaxy ρ DM (r) is modelled according to a Navarro-Frenk-White (NFW) density profile [85], where we set the parameters ρ s = 0.014 M pc −3 and r s = 16.1 kpc [86].

D. Density Profiles of AMCs
For the internal density profile of the AMCs ρ int (R), we consider two different models, namely i) a self-similar Power-law (PL) profile, and ii) an NFW profile as in Eq. (8).These density profiles are illustrated in Fig. 1. 6  The density profile of an AMC for case i) is described by [63,87] where Θ (x) is the Heaviside step function.We truncate the PL profile at a radius We fix ρ s r 9/4 s = ρ AMC (δ)(R PL AMC ) 9/4 /4 [63], to give mean density ρ AMC (δ) and the correct total mass for the AMC.Such PL profiles are expected from models of secondary infall [88] and have also been observed in numerical simulations of the gravitational collapse of AMCs [55].The PL profile likely describes the density of AMCs at formation and therefore would be most suitable for the lightest AMCs which have not undergone growth through mergers.
The most massive AMCs are formed hierarchically through mergers of smaller AMCs, as is generically expected for cold DM substructures [62,63,79].This motivates model ii), in which the density is described by an NFW profile, defined in Eq. ( 8).The exact correspondence between (M AMC , δ) and the NFW parameters is somewhat arbitrary.We follow Ref. [63] and make the identifications ρ s = ρ AMC (δ) and where the function f NFW (c) = ln(1 + c) − c/(1 + c) is defined in terms of a 'truncation' parameter c = R AMC /r s , which we fix to c = 100. 7The AMC profile is truncated at the radius where the numerical result holds for δ = 1.With these choices, the mean density enclosed within r s is 3f NFW (1)ρ AMC (δ) ≈ 0.58ρ AMC (δ) and the total mass enclosed within R AMC is M AMC .
Numerical simulations suggest AMC concentrations of c ∼ O(100) at z = 99 [79], growing roughly as (1+z) −1 to c ∼ O(10 4 ) today [64].However, AMCs with such large concentrations would have a density in their outskirts which is many orders of magnitude lower than that of the host halo of the MW.We therefore fix c = 100 to account for the fact that such diffuse AMCs would have been rapidly tidally truncated.This tidal truncation will lead to a reduction in the AMC mass compared to the initial mass described in § II B. As described in more detail in § IV A and Appendix A, this amounts to a mass loss of 5-40% (depending on the initial AMC mass).For NFW AMCs, we therefore correct the initial HMF to account for this mass loss, as described in § V A. Note also that as we show in Appendix E, fixing the truncation parameter to c = 10 4 instead should have a minimal impact on our formalism and results.
Even fixing c = 100, the above choices lead to NFW miniclusters which are much more dilute than the PL case, allowing us to also explore a more conservative scenario.For a given M AMC and ρ AMC , AMCs described by this NFW profile (with c = 100) will have a mean internal density which is O(10 5 ) times lower than the corresponding PL profile, as illustrated in Fig. 1.
Recent N-body simulation [79] suggest that the transition from direct collapse (PL-like profiles) to hierarchical structure formation (NFW-like profiles) should occur at around M ∼ 10 −13 M for an axion of mass m a = 20 µeV.However, dedicated cosmological simulations are required to confirm the detailed behaviour of AMC density profiles as a function of M AMC .We therefore perform our analysis assuming that either all AMCs have PL profiles or that all have NFW profiles, in the spirit of bracketing the uncertainties on the final AMC properties.

E. Axion Stars
Non-relativisitic ASs are described by solitonic solutions to the Schrödinger-Poisson equation and are expected to form in the central regions of AMCs in the right conditions.In particular, the central density of an AMC must be high enough to allow two-to-two processes to cool their inner core and lead to the formation of a Bose-Einstein condensate [53,89].This process has been observed in recent numerical simulations [70][71][72] and has shown a characteristic core-halo mass relation [90] where we have evaluated the expression today and ignored O(1) factors.The corresponding radius is given by [90] The inverse relationship between the AS's mass and radius leads to a problematic scenario for low-mass AMCs in which the central AS's radius would be larger than that of the corresponding AMCs.To avoid this unphysical description of an AMC, we perform a cut on the overall population in which we only consider miniclusters with a radius larger than the radius of the corresponding AS at its center, i.e.
These results will be referred to as the 'AS cut'.
For the smallest overdensity parameter that we consider δ = 0.1, we find that no AMCs with masses below 5.0 × 10 −16 M pass the AS cut for PL profiles, while no AMCs below 1.6 × 10 −18 M pass in the case of NFW profiles.In both cases, these minimum masses exceed the value of M min in Eq. (5).Starting from the initial population of AMCs described in this section, we then find the fraction of AMCs which pass the AS cut is f PL cut = 2.7 × 10 −4 for PL density profiles and f NFW cut = 1.5 × 10 −2 for NFW profiles.The difference between the two density profiles arises because for a fixed mass M AMC and characteristic density ρ AMC , PL profiles are more compact and the AS radius in Eq. ( 14) is more likely to exceed the AMC radius.We emphasize that current numerical simulations (for example, Refs.[57,71,79]) do not have sufficient resolution to observe the formation of ASs in the lightest AMCs and therefore their existence and evolution has not yet been confirmed.However, our aim is to cut out AMC-AS systems which are likely to be most problematic.Even with this cut, it is also possible that the central density core produced by the presence of an AS may affect the stability of AMCs to tidal perturbations.The treatment of light AMC-AS systems requires dedicated study and is left to future work.

III. TIDAL STRIPPING OF AXION MINICLUSTERS
AMCs can be disrupted by their encounters with stars [91] as well as by tidal interactions with the gravitational field of the disk [67].In this section, we aim to model the interactions of stars with AMCs.Importantly, we model and track all interactions, including those that do not lead to the total disruption of an AMC.Through many successive weak interactions, these AMC's can lose mass and potentially have greatly enlarged radii.This population of perturbed AMCs may result in quantitatively distinct observational signatures when compared to an unperturbed population (see Sec. VI and Ref. [59]). 8 First, we describe how to treat an individual AMC going through a series of interactions.We then discuss in Sec.IV our Monte Carlo procedure to model a population of AMCs being perturbed. 8These encounters may leave a stream of axions behind them which can also lead to features in direct detection experiments [68,92], but here we focus on the properties of surviving AMCs.

RAMC
< l a t e x i t s h a 1 _ b a s e 6 4 = " 9 B L S b h A E u S t R 3 0 9 S 6 p I f l T R n e I U 3 5 9 F 5 c d 6 d j 0 V r w c l n j u E P n M 8 f t U m M 4 g = = < / l a t e x i t > FIG.
2. Illustration of an AMC-star encounter taking place at a distance r from the Galactic center.An AMC with radius RAMC is passed by a star of mass M with impact parameter b and relative velocity V .The energy ∆E injected into the AMC is given in Eq. ( 16).

A. Encounter Dynamics
Stars are dense objects with relatively small radii.Similarly, AMCs are small, meaning that the large majority of encounters will occur when the separation between these objects is significantly larger than their physical size.We therefore work in the 'distant-tide' approximation [93].In this approximation, a minicluster of mass M AMC going through an encounter with a stellar object would increase its internal energy by a quantity [94] (see also Refs.[95][96][97][98]): where M is the mass of the stellar object, b is the impact parameter of the interaction, V is the relative velocity between the objects, and the mean squared radius R 2 accounts for the mass distribution inside the AMCs [95].We parametrize the mean squared radius as R 2 = α 2 R AMC 2 , with α 2 = 3/11 ≈ 0.27 for the PL profile and α 2 ≈ 0.13 for the NFW profile.Such an encounter is illustrated in Fig. 2.
The size of the energy injection, as described by Eq. ( 16), should be compared with the binding energy of the AMC, which we write as E bind = βGM AMC 2 /R AMC .The O(1) prefactor β depends on the internal density profile for which we find β = 1.5 for the PL profile and β = 3.46 for the NFW profile.There are then two distinct regimes for the energy injection:9 • An encounter with a sufficiently small impact parameter will inject more energy than the binding energy E bind of the AMC leading to complete disruption.
• An encounter with a large impact parameter that simply injects energy into the AMC but does not completely unbind it.
The first regime (∆E E bind ) can be re-expressed as b b min , where we have defined the minimal impact parameter that does not entirely disrupt the minicluster as Here, ρ = 3M AMC /(4πR AMC (δ) 3 ) is the mean density of the AMC.An encounter between a PL AMC with δ = 1 and a perturbing object of mass M = 1 M with a relative velocity V = 10 −3 c gives b min ≈ 0.01 pc, which is much larger than the typical size of an AMC, as is required by the distant-tide approximation.Note that this expression depends only on the density of the minicluster, and not on its size or mass separately.Indeed, the fractional energy injected ∆E/E bind depends on the AMC properties only through the mean density, and we therefore expect that the behaviour of the AMCs under perturbations should be independent of M AMC .As pointed out in § II D, for a given mass M AMC and overdensity δ, the mean density of an AMC is significantly lower for NFW profiles than for PL profiles.As we will see, AMCs with NFW profiles are much more easily disrupted than their PL counterparts.
The second regime occurs for larger values of the impact parameter b > b min .In this regime, a single encounter does not completely unbind the AMC, but energy injected through multiple encounters can lead to massloss or a change in radius and may eventually disrupt the AMC.We study this second regime in more detail below.

B. Perturbing the Miniclusters
To estimate the mass loss from a minicluster when it is perturbed, we study the evolution of the phase space distribution function of axions in the minicluster: For isotropic, spherically symmetric distributions of particles, the distribution function depends only on their specific relative energy where Ψ(R) = −Φ(R) is the gravitational potential relative to the boundary at infinity [93,Ch. 4.3].
For spherically symmetric systems in equilibrium, the potential is a monotonic function of the radius R, meaning that the density profile can be expressed as a function of Ψ, ρ(R) = ρ(Ψ(R)).The distribution function can then be determined from the density profile using the Eddington inversion method [93, p. 290]: As discussed in § II D, we will consider two possible density profiles for the AMCs: Power-law (PL) and NFW.
For the PL profile, the distribution function can be computed analytically (see e.g.Ref. [99]) while the NFW distribution must be computed numerically.See Appendix D for more details.Consider then a perturbation to the minicluster of total size ∆E.Given that the critical impact parameter for disrupting the minicluster is already much larger than the minicluster size, we will assume that b R AMC .Under this condition, the average energy injected per unit mass increases with distance from the AMC center as R 2 [94,95] and can be written: Particles with E < 0 immediately after the perturbation can be considered unbound; numerical simulations of the disruption of stellar clusters suggest that the subsequent relaxation of the system should not substantially change the fraction of particles which are unbound [100,101].So in order to compute the mass loss we need only calculate the minicluster mass in particles with energy E < ∆E.This is given by: where the escape speed is v max (R) = 2Ψ(R) and we have also used v = 2(Ψ(R) − E) and dE = −v dv.All calculations in this section are done assuming a fixed potential Ψ(R) as given just before the interaction, which is justified in the limit where the mass loss is small and mostly happening at the outskirts of the object.Note that we have defined the distribution function f (E) for an isolated AMC which extends to infinity.In practice, in Eq. ( 22), we implement a hard truncation of the AMC at a radius R = R AMC .We assume therefore that particles at R > R AMC are no longer bound to the AMC but are instead bound to the diffuse halo of the MW.Note also that our definition of the distribution function is not strictly consistent with a truncated density profile, as we have calculated the potential assuming that the AMC extends to infinity (for convenience).Physically, this means that particles near R = R max are moving more quickly (and are therefore more easily unbound) than they would be in a self-consistent model for the miniclusters.We therefore consider this calculation as conservative from the perspective of AMC disruption.Nonetheless, the error induced by this approximation should be small because, as we will see, these density profiles are close to being in virial equilibrium.
In Fig. 3, we plot the mass loss as a function of the size of the perturbation ∆E, expressed in terms of the binding energy E bind = βGM AMC 2 /R AMC .The fraction of mass lost in an encounter grows with the size of the perturbation, tending slowly to ∆M ∼ M for ∆E E bind .The 'flattening' in ∆M/M occurs because energy is predominantly injected into particles in the outskirts of the AMC which are only loosely bound.Very large amounts of energy are required to strip away the tightly bound particles close to the center.Our calculations are in line with results from N -body simulations of stellar clusters, which show a mass loss of 20-30% for energy injections ∆E ∼ E bind (see e.g.Fig. 6 of Ref. [101]).
Once some mass has been stripped away from the minicluster, we must compute the properties of the surviving object.The total energy of the minicluster is: The velocity dispersion squared σ 2 can be computed from the distribution function and may be parametrized as σ 2 = κGM AMC /R AMC .The total energy can then be written allowing the energy of the AMC to be related to its mass and radius.Objects in virial equilibrium should have κ = β.As we have noted above, the artificially truncated profiles we consider are not strictly in equilibrium, leading to values of κ = 1.15 β for PL profiles and κ = 1.02 β for NFW profiles.In Table I we collect the numerical values of the prefactors obtained for each expression.Energy conservation implies that  where the subscripts i and f denote quantities defined just before and after the interaction.Superscripts 'bound' or 'unbound' refer to the respective subsets of particles, as defined through Eq. ( 22), and we use

PL NFW Expression
We assume that unbound particles are removed instantaneously to infinity.Energies include both the kinetic energy of the particles as well as their potential energy.
We can estimate E unbound f by taking the initial (preinteraction) energy of the subset of particles that is going to be unbound E unbound i and adding the energy that is transferred to these particles during the interaction.This yields where f ej is the fraction of the total injected energy that goes to unbound particles.Putting these components together, we can write the final energy of the AMC after collision as The fraction f ej can be calculated by performing a similar calculation as for the mass loss in Eq. ( 22), but weighting the integral by the amount of energy injected at each radius ∆E(R), as given in Eq. (21).We find that f ej depends on the size of the perturbation ∆E and is typically a factor of a few times larger than the fraction of mass ejected ∆M/M AMC , as illustrated in Fig. 3.The initial energy of the particles which will eventually be unbound E unbound i can be taken as the sum of the kinetic energy of these particles plus the change in the binding energy from removing these particles.The final binding energy is calculated self-consistently from the density profile immediately after the interaction (see Appendix D for more details).
We assume that after the perturbation, the AMC will relax on a short timescale to have the same density profile (and for NFW profiles, the same truncation parameter c = 100), but described by a new mass and radius.This assumption is made for computational simplicity -however, there is some indication from N-body simulations that perturbed DM substructures do retain a universal density profile [98,101]. 10The final mass is M f AMC = M i AMC − ∆M , while the final radius can be calculated from the total energy using Eq. ( 24): We note that the internal relaxation time scale t rel ∼ R AMC /σ v ∼ 10 4 yr of the AMCs is several orders of magnitude shorter than the average time between substantial encounters, ∆t enc ∼ t MW /N enc ∼ O(10 6 ) yr, where t MW ≈ 13.5×10 9 yr is the age of the MW and N enc ∼ 10 4 is the typical number of encounters (see Fig. 7).We therefore assume that there should be sufficient time for AMCs to relax between successive encounters.
The assumption that the AMCs will have the same density profile after the perturbation allows us to follow the evolution of a large number of AMCs under many perturbations in a computationally feasible way.However, we note that this assumption is conservative.From Eq. ( 21), energy is injected into and mass is lost predominantly from the outer parts of the AMC.The remnant should therefore be more dense after the perturbation, making it more resistant to further perturbations.Our assumption therefore leads to a smaller number of surviving AMCs than a more detailed (but computationally expensive) calculation.
In Fig. 4, we illustrate the evolution of the mass, radius, and binding energy of a minicluster under repeated perturbations.We fix the size of each perturbation to be 1000 times smaller than the initial binding energy, ∆E = 10 −3 E bind,i .In the PL case (left panel), we see that with each perturbation, mass is lost but energy is also injected and as a result the AMC becomes larger in size and is eventually disrupted entirely.This typically takes place with fewer encounters than expected without accounting for the evolution of the minicluster properties (which would be O(10 3 ) in this case, due to the fixed size of the perturbations).Miniclusters with NFW profiles (right panel) show a similar behaviour.Crucially, when ∆E becomes comparable to the remaining binding energy (after O(450) encounters in this example) the AMC radius begins to decrease.This is because for large ∆E/E bind , almost all of the injected energy is carried away by the ejected mass (f ej tends to one in Fig. 3), leaving the remnant smaller and more dense than before the interaction.This behaviour -seen also in studies of tidal shocks in globular clusters [102] -emphasizes why we must carefully account for the redistribution of energy through Eq. ( 27).

IV. MONTE CARLO SIMULATIONS
Having described how individual miniclusters are perturbed, we now want to understand the overall population of AMCs and their interactions over the history of the MW.We therefore run Monte Carlo simulations for this population of AMCs.For simplicity, we make a distinct split between structure formation and the stellar disruption of AMCs -these two processes would typically happen concurrently (see § IV A).Our simulations are therefore initialized with the expected properties of an AMC population today, at z = 0. We will discuss how removing this assumption could affect our results in § IV A. N -body simulations of a population of AMCs may be required to fully understand the details of the population today.Nevertheless, we attempt to capture the relevant physics in a simple and interpretable process.
We perform two different simulations using either the PL or NFW AMC density profiles.In reality, we expect the majority of AMCs to show an NFW profile for M M 0 and a PL profile for M M 0 , with the intermediate mass region being populated by both.Our simulations therefore attempt to bound the range of possible AMC populations today.We run the simulations for t MW = 13.5 billion years, therefore allowing for the maximum amount of disruption within the lifetime of the MW.The stellar distribution in the MW is assumed to be static in time and is modeled as a bulge plus a disk, as described in Appendix B. Again, we will discuss how changing these assumptions could change our results in § IV A.
Our simulated miniclusters are distributed in a spherically symmetric halo and follow elliptic orbits with a focus at the Galactic center, where eccentricities follow the distribution shown in the top panel of Fig. 5 [103].We show later in Fig. 10 that this reproduces the Galactic NFW density profile.For a given orbit, the time variation of the galactocentric radius r of the orbit is described by the expression where a and e are the semi-major axis and the eccentricity of the elliptic orbit respectively.For simplicity, we take M encl to be the mass of the MW enclosed within a sphere of radius a given by the NFW profile [85].The orbital radius is bounded by the values a(1 − e) ≤ r ≤ a(1 + e) and the period is given by T orb = 2π a 3 / (GM encl ).The probability of finding an AMC at a specific radius at a particular instant in time is given by In Fig. 6 we show the binned distribution of P (r) for the values of the eccentricity e = 0.1 (blue), e = 0.5 (orange), and e = 0.9 (green).For small values of e, P Importantly, it shows that the majority of AMCs found at a particular radius will have similar semi-major axes (r ∼ a) and relatively low eccentricities e 0.4.At the same time, Fig. 5 shows a non-negligible fraction of the population of AMCs found at a particular galactocentric radius will have highly eccentric orbits and a variety of semi-major axes.
Our Monte Carlo simulations range over a logarithmic grid of values for the semi-major axis a ∈ [0.1, 50] kpc.
For each value of a, the simulation proceeds as follows: Top panel shows the eccentricity probability distribution for the orbits of AMCs taken from Ref. [103].Bottom panel shows the joint probability that an AMC will have a particular semi-major axis and eccentricity given a particular galactocentric radius r = 8 kpc.
1. We generate a set of N AMC = 10 5 AMCs where the mass is sampled from a log-flat distribution between M min and M max and in overdensity according to the distribution df AMC /dδ in Eq. (C1).The log-flat sampling was used to ensure that we had a sufficiently large number of high mass AMCs in our simulations.We will discuss how to re-weight the results of the simulations (and how we apply the AS cut) to recover the true distribution in Sec.V. We draw the eccentricity of each AMC orbit from P (e), which we treat as independent of the galactocentric radius.Finally, each AMC is also given a random inclination angle ψ with respect to the Galactic plane, uniformly sampled within 2. For each AMC orbit, we compute the number density of the stellar field encountered by the AMC as a function of time n (t) over a full orbit (see Appendix B).We then evaluate the total number of encounters with stellar objects over an orbit as where b max is the maximum impact parameter that we consider (see below) and V AMC (t) is the velocity of the AMC as a function of time.The total number of interactions is then given as N int = N × (t MW /T orb ).We truncate the total number of interactions at N cut = 10 6 which we justify below.
3. For each AMC, we sample N int relative velocities and impact parameters.From these we compute a list of energies E list to be injected with each encounter, using Eq. ( 16).The relative velocities are calculated by sampling from the integrand of Eq. ( 31) to compute a distribution of the most likely interaction times.These interaction times can then be converted into a list of radii and AMC velocities using the Vis-Visa equation [104].To obtain the encounter velocities, we then add a random 3-D velocity drawn from the local stellar velocity distribution, which we take as a Maxwell-Boltzmann distribution with dispersion σ v = GM encl (r)/r.The impact parameter is randomly drawn from the probability distribution de-fined as We fix the maximum impact parameter b max such that the energy injected is 1/N cut times smaller than the initial binding energy of the minicluster, ∆E(b max ) = E bind /N cut = 10 −6 E bind .The truncation at N cut is therefore not physically relevant for completely disrupted AMCs but is sufficiently large to capture the effects of partial disruption.We find typical values of b max ∼ 10 −2 pc and b max ∼ 10 −1 pc for PL and NFW profiles respectively.
4. We then iteratively perturb each AMC, through E list , using the prescription described in § III A.
We recompute the new radius and mass after each iteration.As we note in § III B, in some cases the density of an NFW minicluster can increase after an encounter, making it more resistant to further disruption.When this happens, we recompute b max using the procedure described in the previous step.
We then recompute the number of interactions in the remaining simulation time and truncate this at N cut .If the total AMC energy, given in Eq. ( 23), climbs above zero, we consider it to be completely disrupted and remove it from the simulation.Note that we do not keep track of which AMCs pass the AS cut during the simulations, instead applying the AS cut to the distribution of perturbed AMCs, as described in § V A.
Histograms of the number of interactions as a function of the galactocentric radius can be seen in the upper two panels of Fig. 7 -to simplify the discussion, we show results for AMCs on circular orbits only.Note that here we count the number of simulated interactions for each AMC, stopping either at the end of the simulation time or when the AMC is disrupted.There is a clear difference between the PL (top) and NFW (middle) simulations which can be seen as a distinct shift to larger numbers of interactions for NFW profiles for a fixed galactocentric radius.The shift originates from the smaller average density of the NFW miniclusters, which is reflected as a lower average binding energy.The reduced binding energy leads to more interactions above the threshold of ∆E/E bind > 10 −6 .In the NFW case, we also see that the number of interactions rarely extends beyond 10 4 .This is because the AMCs are either completely disrupted or are stripped to leave a high density remnant (for which further interactions above the threshold of ∆E/E bind > 10 −6 are rare).
The lower panel of Fig. 7 shows the value of ∆E/E bind for a sample of 10 5 interactions in the Monte Carlo simulation.The majority of interactions inject only a small amount of energy, while about 1 in 1000 interactions are strong enough to give rise to substantial mass loss ∆E ∼ E bind .This further justifies our cut on the maximum number of interactions N cut = 10 6 , as this is much larger than the ∼ 1000 interactions which would typically be required to unbind an AMC.

Number of interactions
Circular orbits, PL profile FIG. 7. The top and central panels show the number of stellar interactions with AMCs on circular orbits at various galactocentric radii.Compared to AMCs with PL density profiles, miniclusters with NFW density profiles undergo significantly more interactions with ∆E/E bind at a fixed galactocentric radius.For illustration, in the bottom panel we show the distribution of injected energies for PL profiles (which is approximately independent of galactocentric radii).

A. Assumptions and Caveats
Milky Way Properties -We have so far neglected the disruption of AMCs due to tidal stripping by the host halo of the Galaxy.The impact of tidal stripping can be quantified by considering the tidal radius of the AMCs, the distance from the center of the AMC at which tidal forces from the MW halo become comparable to the self-gravity of the bound AMC [105].We find that for PL profiles, the tidal radius is several orders of magnitude larger than the physical radius of the AMC, making them robust to tidal stripping.Instead, NFW profiles with c = 100 may be comparable in size to their tidal radius, especially at small Galactocentric radii.It is therefore likely that NFW AMCs have undergone some tidal stripping by this mechanism, from their initial concentration of c ∼ 10 4 , to reach our assumed concentration of c = 100.We therefore apply a correction of 5 − 40% to the initial mass of NFW AMCs, to account for this massloss.However, given the small size of these corrections, we conclude that complete disruption by this mechanism is unlikely and subdominant to stellar disruption (see also Ref. [105]).Full details concerning tidal stripping due to the MW halo are given in Appendix A.
The orbital motion of an individual AMC is influenced by dynamical friction exerted by the MW [93, p. 644].The orbit of an AMC with virial velocity v(r) at the galactocentric radius r decays with a timescale where ξ(r) ∼ O (10) is a dimensionless function and ρ(r) is the background density of the MW at the orbital radius r.Setting ρ(r) to the Galactic NFW distribution, with the corresponding expression for the virial velocity, Eq. ( 33) gives t frc t MW for M 10 −5 M and r 0.1 pc.Conversely, the orbits of the heaviest AMCs would be destabilized at very small Galactocentric radii.However, as we will see, tidal disruption by objects in the stellar bulge would disrupt these AMCs well before t MW .We can therefore ignore the effect of orbital decay.
Throughout this work, we only include tidal interactions with stars.In particular, we account for ∼ 1011 stars and fix their mass to be 1 M .Since the stellar mass function is relatively steep, 11 the vast majority of stars are around 1 M .We therefore expect that considering different stellar masses will produce only a small correction to our results.In addition, we have not considered tidal interactions with a separate population of NSs or white dwarfs which, despite having a mass of the same order as that of a typical star, are at least an order of magnitude less numerous [107].Again we expect the corrections to our simulation results to be small when accounting for these additional astrophysical objects.
Finally, we neglect variations in the stellar density over the lifetime of the MW.Since the lifetime of a Solar mass star is O(10 10 ) years, the majority of stars born early in the MW's history will have finished their life cycle by today.It is therefore important to understand whether changes to the stellar abundance could affect our results.Luckily, the star formation rate is much larger than the death rate [108], meaning that the stellar density has been increasing throughout the lifetime of the Galaxy.By using a stellar density as measured today we are therefore overestimating the amount of tidal stripping that could happen over the lifetime of the Galaxy.Nevertheless, the total stellar mass of galaxies like the MW is thought to be relatively constant (within a factor of two) since z ∼ 1 [109].Future work should adopt a time varying model of the MW which follows the cosmological star formation history [110].
Structure Formation -As mentioned above, for computational simplicity we made a distinct split between the hierarchical structure formation that these AMCs will undergo and their tidal stripping through interactions with stars.The interplay of these two physical effects requires a detailed study which we leave to future work.Nevertheless, our split represents a conservative approach since we allow for the maximum amount of tidal stripping to occur for all AMC masses.Lighter AMC's are more abundant and are likely to have experienced fewer merger events than their heavier counterparts.These lighter AMC's have therefore been present in our Galaxy for the longest period of time -our procedure should therefore be a good reflection of the tidal stripping for lighter AMCs.Heavier miniclusters, on the other hand, are less abundant and have been gradually merging throughout the history of our Galaxy.Mergers gradually increase the maximum mass that an AMC can achieve.For the heaviest AMCs, our simulations overestimate the amount of tidal stripping that may have occurred by today.On the other hand, Figs. 5 and 7 of Ref. [63] show that AMCs with M 10 5 M 0 ≈ 10 −6 M collapsed before z ≈ 10 and therefore substantially before the formation of the MW.For these lower masses our simulations should capture the effects of stellar tidal interactions very well throughout the MW halo.
For the most massive AMCs, from 10 −6 M up to M max ≈ 5 × 10 −5 M , there is still some uncertainty.At smaller galactocentric radii the survival probability is low (r 10 kpc for NFW profiles and r 3 kpc for PL profiles).For high mass AMCs (which formed through mergers) in these regions, we neglect the effects of concurrent structure formation and stellar interactions, leaving this to future work.In the outskirts of the MW halo r 30 kpc, stellar encounters are quite rare and will therefore not dramatically affect the growth of more massive AMCs.Fortunately, the high mass AMCs in the inner regions of the MW are a tiny proportion of the total number of miniclusters (O 10 −9 ) so we therefore con-clude that our simulations are accurate for the majority of the AMC population.
Mutual AMC collisions -Throughout our simulations, we do not consider the mutual interactions between miniclusters.To see that this is a good approximation, we estimate the encounter rate of two AMCs at the Galactocentric radius r as Γ ∼ n(r)v(r) R 2 where n(r) ∼ ρ MW (r)/ M AMC and v(r) is the virial velocity associated with the NFW profile.These encounters occur rather frequently, for instance at r ≈ 4 kpc the encounter rate is Γ ≈ 10 5 years −1 .For comparison, the same computation for the encounter of a minicluster with a star yields Γ ≈ 10 19 years −1 .However, mutual AMC encounters only deposit a small amount of energy during an interaction -this can be seen from Eq. ( 16) which scales as the square of the mass of the perturbing object.The mean mass of an AMC is ∼ 10 −14 M , meaning that a typical AMC-AMC interaction will deposit 10 −28 times less energy than a typical AMC interaction with a star.Therefore, despite their large interaction rate, mutual AMC encounters will not significantly contribute to the tidal disruption of miniclusters.
In addition to AMC collisions that lead to tidal disruption, mutual miniclusters interactions can also lead to mergers.Importantly, the increased background density within the MW with respect to the critical density will cause these merger interactions to happen more regularly than in typical simulations.Fortunately, this effect will be most prominent in the Galactic center where the density is largest, but also where stellar disruption will be dominant.We leave a complete study of this effect to future work.
Minicluster Substructure -As discussed in Sec.II D, we expect the most massive AMCs to form from hierarchical mergers of lighter ones.Numerical simulations of AMC clustering show that an internal granular structure is expected on top of an overall NFW distribution [79].We have not considered such a granular substructure, since lighter AMCs are expected to be dissolved within the larger minicluster over the lifetime of the MW due to the effects of dynamical friction [65].

V. AXION MINICLUSTERS TODAY
In this section we first discuss how to construct the true AMC population distributions today from our Monte Carlo simulations followed by a discussion of the results.Note that although we present results for a specific choice of HMF, the reconstruction procedure allows us to use the same Monte Carlo results for arbitrary HMFs.In particular, this could include changes in the mass cutoffs M min and M max , changes to the AS cut, or the introduction of correlations between the minicluster mass M AMC and overdensities δ.The full suite of results and the corresponding code to reinterpret them can be found at github.com/bradkav/axion-miniclusters[76].We show the probability distributions of the mass MAMC, radius RAMC, and mean density ρ.We assume these to have NFW internal density profiles and to be on circular orbits with a galactocentric radius of r = 6.96 kpc, leading to a survival probability psurv = 0.91.These probability distributions are obtained through the reconstruction procedure described in § V A. The grey shading indicates the regions of the parameter space that are removed by the AS cut.The white lines shows the smallest value of the corresponding parameter that passes the AS cut.

A. Reconstructing Physical Properties
Each of our Monte Carlo samples corresponds to N = 10 5 AMCs with orbits of a given semi-major axis a.In order to calculate the survival probability and AMC properties as a function of galactocentric radius r, we must assign each sample a weight w, proportional to the time that AMC spends at a given value of r, from Eq. (30).With this, we essentially 'smear' each AMC sample over a range of radii r, allowing us to estimate the properties as a function of r (instead of a).We assume that the initial number density of miniclusters as a function of semi-major axis n AMC (a) follows an NFW profile.For a single AMC with semi-major axis a i and eccentricity e i , the weight assigned in some radial bin ∆r is then: The first term accounts for the fact that the Monte Carlo samples are not uniformly distributed in a, but concentrated on logarithmically-spaced grid points, with spacing ∆a i .The second term is the assumed initial probability distribution for the semi-major axis P (a), with n AMC defined in Eq. ( 7).The final term is the fraction of time spent at a given radius, defined in Eq. ( 30), averaged over the radial bin of interest.The AMC number density at a given Galactocentric radius can be obtained by summing over the weights of all AMCs (with potential contributions from all values of a).This smearing procedure gives rise to an approximately NFW profile as a function of galactocentric radius n AMC (r), as shown by the black dashed line in Fig. 10.The number density of AMCs at a galactocentric radius r can be written as: where the survival probability is defined as the ratio of the number of surviving AMCs N surv (r) to the number N initial (r) at a given radius: p surv (r) = N surv (r)/N initial (r).The probability distribution for the minicluster mass and radius at a given galactocentric radius P (M, R|r) can be extracted from the Monte Carlo simulations. 12ur Monte Carlo simulations were performed with a log-flat distribution of AMC masses.The results must then be adapted to reflect the true mass function of the AMCs.From our simulations, we obtain the final density ρ and the mass-loss fraction ν = M f /M i of each AMC in a sample. 13The simulations confirm that the distribution of ρ and ν do not depend on the initial AMC mass but only on the initial density, as discussed in § III A.
For a given initial mass function P i (M i ), the final distribution of masses can be obtained by integrating over all possible initial masses M i , selecting only those which produce the correct final mass M f .The joint distribution of AMC mass and density after disruption can then be written as: Here, P (ρ, ν|r) is the final probability distribution for the density and mass-loss fraction at a given galactocentric radius.We can now write the final mass function as a Monte Carlo integral: where we have replaced the integral by a sum over the N surviving AMCs, with properties (ρ k , ν k ), distributed according to P (ρ, ν).The index k runs over the AMCs in the sample, with the weights w k calculated at a given galactocentric radius, as described above.Here, AMCs which have been completely disrupted are excluded from the integral (or the corresponding sum).The distribution of final radii R f can also be written as: and similarly for any other distribution of interest.We fix the initial mass function according to Eq. ( 4), with slope γ = −0.7.For NFW AMCs, we also apply a correction of 5 − 40% to the initial AMC mass to account for tidal stripping due to the MW host halo.This is implemented directly in the definition of P i (M ) for NFW AMCs.Full details are given in Appendix A.
Since the disruption process is mainly sensitive to the densities of the AMCs (and not the AMC masses), we do not expect a significant difference to our results when changing the HMF.For consistency, we therefore re-run the entire pipeline using a different slope γ = −0.5 and find only a ∼ 10% increase in the survival probability (compared to γ = −0.7)for NFW AMCs at the Solar position.This is not due to a change in the disruption properties of the AMCs but rather a change in the fraction of AMCs passing the final AS cut.There is no appreciable change for PL miniclusters.We therefore do not consider γ = −0.5 further.
In our numerical results, we assume that all axions are bound in AMCs (f AMC = 1) with masses between M min and M max , given in Eq. ( 5).For our 'AS cut' results we use the total perturbed sample of AMCs and require that the AS radius be smaller than the AMC radius.This reduced sample is then compared to the unperturbed sample with the same cut applied (as described in § II E).As mentioned in § II E, this cut effectively reduces the fraction of axions bound in AMCs.
The AS radius in Eq. ( 14) can be re-written as where we define the constants R = 1.7 × 10 −6 pc and M = 10 −16 M .The AS cut therefore requires that where we calculate the AS radius using the initial AMC mass, assuming that the properties of the central AS are unaffected by perturbations.In Eq. ( 37) and Eq. ( 38), this cut is equivalent to summing only over those samples which satisfy

B. Results
In Fig. 8, we show an example of the reconstructed probability distributions of AMC properties at the end of our Monte Carlo simulations.Specifically, we show results for AMCs with NFW density profiles and, for simplicity, on circular orbits with Galactocentric radius of r = 6.96 kpc.In this case, only around 9% of miniclusters are destroyed -this can be seen in Fig. 9 from the olive dashed line -but the properties of those which survive are drastically altered.The results in Fig. 8 do not include the AS cut.However, the shaded areas indicate regions of the parameter space which are progressively removed by this cut.In particular, no AMCs pass the AS cut to the left of the vertical white lines.
Focusing on the left panel of Fig. 8, we see that the lowmass tail of the distribution of AMC masses extends to lower values after disruption is taken into account.This is due to mass loss from miniclusters which are perturbed but not completely disrupted (see § III B).In some cases, the final mass of an AMC is reduced by several orders of magnitude compared to its initial mass.We have verified also that the fractional mass loss does not depend on the initial mass of the minicluster but only on its density.
The right panel in Fig. 8 demonstrates the strong dependence of the disruption on the initial AMC density.From our discussion in Sec.III, we expect the amount of energy injected per encounter to scale as ∆E/E bind ∼ R 3 AMC /M AMC ∼ 1/ρ AMC .Indeed, we see that very dense miniclusters undergo little disruption.Instead, less dense miniclusters (e.g.around 1 M pc −3 ) have a low survival probability and those which survive lose mass.If the AMCs undergo only small perturbations, they may increase in radius through energy injection, leading to a tail of diffuse miniclusters (down to ∼ 10 −3 M pc −3 ).Instead, large perturbations can cause substantial mass loss from the AMCs, but very little energy is injected into the remnant.This leads to an overall increase in the typical AMC density.The distribution of AMC radii (central panel in Fig. 8) follows from this same argument.The typical AMC is more dense after accounting for stellar interactions, leading to a reduction in the AMC radius.
In Fig. 9, we show the survival probability of our simulated miniclusters as a function of both their semi-major axes (left) and galactocentric radii (right) for both PL and NFW internal density profiles.The right panel shows results for both eccentric and circular orbits where the former is constructed using the prescription described in In the panel we show the survival probability for both circular (dotted) and eccentric (dashed) orbits.The eccentric orbits are smeared out according to the proportion of time they spend at a given galactocentric radius (as described in § V A).Finally, we show the survival probability for miniclusters with eccentric orbits that pass our AS cut (solid line) where we have normalized this to be one at large galactocentric radii.§ V A. In both panels we see that there is a characteristic transition from the high stellar density inner region r O(1) kpc (where the number of interactions between stars and AMCs is so high that almost all miniclusters are completely disrupted) to a low stellar density outer region r O( 10) kpc (where interactions are rare).This change in the number of interactions can be seen clearly in Fig. 7. Focusing on the left panel of Fig. 9 we see a distinct shift of the transition region to larger radii from the PL to NFW density profiles.This shift comes from the enhanced number of interactions at a fixed galactocentric radius for NFW miniclusters -this can also be seen in Fig. 7.
In the right panel of Fig. 9, we see this same distinct shift from PL to NFW profiles for both eccentric and circular orbits.Moving from circular to eccentric orbits produces a smearing of the transition region caused by the distributions of semi-major axes and eccentricities that contribute to the AMCs at a particular galactocentric radius, as shown in Fig. 5.For example, at low galactocentric radii, the AMC density will have contributions from both quasi-circular orbits which spend a large amount of time in dense stellar regions and highly eccentric orbits that spend the majority of their time in low density stellar environments at larger radii.
We also show results for the survival probability of AMCs which pass the AS cut (solid lines).In this case, we normalize the survival probability to one at large galactocentric radii, which is equivalent to factoring out the initial fraction of AMCs which pass the cut, f PL cut = 2.7 × 10 −4 for PL density profiles and f NFW cut = 1.5 × 10 −2 for NFW profiles.The survival probability with the AS cut is always smaller than without the cut.This is because stellar perturbations can strip the AMCs until their radius drops below the corresponding AS radius.This demonstrates that the AMC properties can be substantially altered by stellar interactions.At the Solar radius, we find a survival probability (including the AS cut) of 99% for AMCs with PL profiles and 46% for AMCs with NFW profiles.
Finally, Fig. 10 shows the density of AMCs as a function of the galactocentric radius.Firstly, we show that we are able to correctly reconstruct the Galactic NFW profile (gray dotted line) using the initial sample of unperturbed miniclusters with NFW-distributed semi-major axes (gray solid line).The fraction of AMCs that are removed through the AS cut is indicated by the reduced normalization of the dot-dashed lines.The density of perturbed miniclusters passing the AS cut is shown by the solid lines.The PL profile AMCs (blue solid line) maintain a population down to small radii -these AMCs are still likely to have undergone many interactions and therefore may have significantly different properties to those at the start of the simulation.In contrast, the NFW miniclusters (olive solid line) show a sharp reduction in density at around r ∼ 10 kpc.

VI. APPLICATIONS
Many experimental probes of axion DM require one to make assumptions about both the large scale halo distribution of DM and its structure on small scales.
Here we briefly investigate the phenomenological consequences of our results for a few primary observational channels.Density of AMCs ρ(r) [M pc −3 ]   Galactic density profile Unperturbed AMCs Power-law NFW Perturbed AMCs Unperturbed AMCs + AS cut Perturbed AMCs + AS cut FIG.10.Density of simulated AMCs before taking into account stellar interactions (solid gray), compared with the expected NFW profile for the DM halo (dotted gray).The spatial distribution of AMCs with PL and NFW internal density profiles are shown as blue and olive lines respectively.The density profiles are extracted from the Monte Carlo simulations according to the weighting procedure in Sec.V. We also show the distribution of unperturbed and perturbed AMCs passing the AS cut (dot-dashed and solid lines respectively).

A. Lensing
Gravitational lensing has been used to constrain the fraction of DM in the form of faint compact objects [111].Surveys actively search for microlensing events caused by objects passing through the line of sight between the Earth and stars in target structures like the Magellanic System [112][113][114][115][116], the inner Galactic bulge [114][115][116], and M31 [117,118].
For simplicity, we limit our discussion to microlensing from a point-like lens of mass M with a point-like background source -in reality both the lensing AMC and the background galaxy are extended objects.To assess the impact of our Monte Carlo results we compute the difference between the expected number of lensed events, perturbed and unperturbed, for a variety of observational targets.We leave a more detailed analysis for future work.
We denote the observer-source, observer-lens, and lens-source distances as D S , D L ≡ x D S , and D LS = (1 − x) D S respectively with 0 ≤ x ≤ 1.We expect the set of microlensing events over some observation time to be Poisson distributed with the expected number of events given by [113] The differential event rate for a single source star with respect to distance and event time, d 2 Γ/(dx dt), depends on the mass distribution of AMCs, their velocity distribution, and is proportional to the number density of lenses TABLE II.Parameters used for the source location and halo modelling.We report the distances of the LMC [119] and M31.We also specify the NFW parameters used for the MW [120], for the LMC [121], and for M31 [122].
along the line of sight n AMC (x).Here we treat n AMC (x) as a proxy for the expected number of lensing events.We consider microlensing events from sources residing in the following targets: the MW Galactic bulge, M31, and the Large Magellanic Cloud (LMC).Miniclusters in either the halo of the MW or in the target's halo could lead to lensing events.We model the DM halo distribution in each target galaxy and in the MW according to the NFW density profile in Eq. (8).In Table II we report the distances between the Solar System and the source D S , together with the Galactic longitude and the latitude ( g , b g ) and the parameters used to model the NFW profile.For the LMC and M31 we use the AMC survival probability obtained in the MW but rescaled by the scale factor r s of the target galaxy.The distribution of AMCs between the Solar System and the source is given by the sum of the profiles along the line of sight.
Results are shown in Fig. 11 for the LMC (top), the MW bulge (middle), and M31 (bottom).We show the number density of AMCs as a function of the ratio of the lens distance to the source distance x = D L /D S for the perturbed population derived from our simulations (solid line), assuming NFW (olive) and PL (blue) AMC density profiles, including the AS cut, and using the same color code as in Fig. 9.For comparison, we also show the results when perturbations are neglected (dot-dashed line).
In Fig. 11, considering observations of the LMC (top panel), the number density of the NFW AMC population is strongly suppressed both towards the inner region of the MW (x 1) and the LMC (x ≈ 1), with respect to the unperturbed one.This behavior is evident also for AMCs along the line of sight towards the Galactic center (middle panel), for both PL and NFW profiles.For M31 (bottom panel) the difference between the populations is not clearly visible, except near the center of M31 where the number density is significantly affected.This is because we fix the line of sight to end at the center of the target galaxy, where the distribution of perturbed AMCs drops significantly.Fixing the line of sight to point at other regions of the target galaxies, where the survival probability of AMCs is higher, would lead to smaller differences.Note that for PL profile AMCs, the number density closely follows the unperturbed distribution for x 1, i.e. within the MW.Compared to those with PL profiles, NFW miniclusters show a more dramatic reduc-tion in the number density along the line of sight for all three sources.This is due to the difference in survival probabilities shown in Fig. 9.
Using Eq. ( 42), we denote the fractional decrease in the expected number of lensing events before and after tidal interactions as δ N ≡ ∆ Nex / Nex .For sources in M31 we expect δ N ≈ 0.5% (PL) or δ N ≈ 18% (NFW), where the discrepancy is almost entirely due to the AMC disruption in the target galaxy.When considering the LMC as a source, we find δ N ≈ 1% (PL) or δ N ≈ 32% (NFW).The largest difference is clearly seen in searches towards the MW bulge where δ N ≈ 12% (PL) or δ N ≈ 92% (NFW).Even when the fractional decrease in the number of events is small, the properties of these events could change dramatically.For instance, the duration of the microlensing events can be significantly shortened since the surviving AMCs are each less massive than in the unperturbed case.A more careful analysis would be required to determine the predicted properties of AMC microlensing events.

B. Direct Detection
The main direct detection experiments that would be affected by the presence of AMCs are 'haloscopes'.Haloscopes convert cosmic axions into a detectable signal inside a resonant cavity immersed in a strong magnetic field [9,10].Once the settings for the cavity have been fixed, the power output P a from the conversion of the cosmic axions in the magnetic field of the cavity is proportional to the local energy density of axions ρ a times g 2 aγγ [9,10].The value of the local energy density has contributions from a smooth component in the Solar neighborhood and from substructures like streams and AMCs.This means that for a fixed value of g aγγ , P a can still vary significantly.In particular, an AMC passing near the Earth would enhance the power output of the cavity for a short amount of time.This is expected to be a rather rare event, as estimates show that an encounter would occur only every 10 4 − 10 6 years [58,81].
Here, we compute the encounter rate for AMCs of radius R AMC with the geometrical cross section of the Earth σ = πR 2 AMC .In principle, a correction due to gravitational focusing becomes important for miniclusters with radii smaller than R = 2G N M ⊕ /σ 2 u ≈ 10 km, where M ⊕ is the mass of the Earth and σ u ≈ O (200 km/s) is the velocity dispersion at the Earth's location.Since miniclusters are much larger than R, we neglect this focusing contribution.The encounter rate between the Earth and the population of miniclusters is where n AMC (r ) is the local number density of AMCs, r = 8.33 kpc is the distance of the Solar System from the MW Galactic center [123] and σu (r ) is the velocityaveraged cross section weighted with the probability distribution of R AMC at r .11.The number density of AMCs (in kpc −3 ) as a function of the distance of the lens x = DL/DS for the unperturbed population (dot-dashed line) and for the total population of perturbed AMCs (solid lines).We show the results for target sources placed in the LMC (top), in the MW bulge (middle), and in M31 (bottom).For the both perturbed and unperturbed populations, we have applied to AS cut.
The local number density is given by Eq. ( 7) times the fraction f cut that accounts for the AS cut in the HMF given in § II E. We compute σu (r ) from the distribution of AMC radii at r = r obtained from the Monte Carlo simulations for both the unperturbed and the perturbed populations, taking into account the full distribution of masses and densities.The number densities of perturbed AMCs are also weighted by the sur-vival probability in Fig. 9. Two competing effects combine in the computation of the encounter rate for the perturbed population.In general, successive perturbations tend to puff up AMCs, making their mean radius larger and encounters with Earth more probable.On the other hand, stellar encounters destroy some of the AMCs (or strip them to below the AS cut), thereby lowering the chance of encounters.For the PL profile the two effects compensate, leading to the same encounter rate Γ ≈ 4 × 10 6 years −1 for both the unperturbed and the perturbed AMC distributions.For the NFW profile, the first effect dominates, leading to Γ ≈ 10 3 years −1 for unperturbed AMCs and Γ ≈ 4 × 10 3 years −1 for perturbed AMCs.The difference in the magnitude of the results between the NFW and the PL populations is due to the larger fraction of NFW AMCs which pass the AS cut as well as the larger typical radius for NFW profiles.
Recall that the results for Γ are inversely proportional to the AMC fraction, which we set to f AMC = 1.
In Fig. 12, we show the enhancement of the local energy density of axions, ρ a /ρ , during an encounter with an AMC, as a function of time.We fix ρ = 0.45 GeV cm −3 and we show time in days.The energy density ρ a is given by Eq. ( 9) for miniclusters with a PL profile (top panel) and by Eq. ( 8) for miniclusters with an NFW profile (bottom panel).We have considered an overdensity δ = 1 and masses M = 10 −10 M (red and black lines) and M = 10 −12 M (blue line), with an impact parameter b = 0.1 R AMC (red and blue lines) and b = 0.5 R AMC (black line).The maximum enhancement is regulated solely by the ratio between the density of the AMC, which is proportional to ρ AMC (δ) in Eq. ( 2), and ρ , while the mass of the AMC controls the duration of the encounter.Encounters with PL profile miniclusters are shorter and have greater enhancements than those with NFW miniclusters.
The small rate and short duration of these interactions makes AMC encounters with Earth mostly irrelevant for direct axion searches.We find that an encounter has a typical duration of O(10) days, giving a probability of O(10 −5 ) and O(10 −8 ) that the Earth is currently inside such an AMC for NFW and PL AMCs respectively.Nevertheless, the AMC fraction f AMC and its evolution in the MW is an important quantity for direct detection.If all axions are locked up in AMCs, direct searches may be ineffective.In addition, the disruption of miniclusters can lead to streams of axions, which can produce an important contribution to the local density [68,92].We will assess these streams more carefully in future work.

C. Indirect Detection
Axion-photon conversion in astrophysical magnetic fields can lead to potentially detectable radio-wave signals.For example, axions can resonantly convert into photons in the magnetic field of a NS [38] and non- 12. The enhancement in the local axion density due to the encounter of Earth with an AMC, in units of ρ = 0.45 GeV/cm 3 , as a function of the duration of the encounter in days.We separately show the impact of a minicluster that follows a PL profile (top panel) and an NFW profile (bottom panel).We have considered an overdensity δ = 1 and masses M = 10 −10 M (red and black lines) and M = 10 −12 M (blue line), with an impact parameter b = 0.1 RAMC (red and blue lines) and b = 0.5 RAMC (black line).We have considered AMCs which have not been perturbed.As we show in the main text, these AMC encounters with Earth would be rare.
resonantly in the Galactic center of target galaxies [124].The latter process is severely suppressed in favor of photon-photon pair production [125][126][127].The search for these radio signals is complementary to laboratory experiments and might help disentangle the coupling g aγγ from the cosmic axion abundance.
The photon production rate in these astrophysical environments depends on the local number density of axions, which is enhanced in the presence of an AMC.In Ref. [59], we compute the expected signals from encounters between miniclusters and NSs.

VII. DISCUSSION AND CONCLUSION
In this paper we quantified the degree to which tidal interactions with stars can affect the final distributions of axion miniclusters (AMCs) in the Milky Way (MW).We performed Monte Carlo simulations for AMCs on circular and eccentric orbits, and with Power-law (PL) and NFW density profiles.Importantly, we quantified the survival probability of AMCs as a function of the galactocentric radius (Fig. 9), showing that in the inner regions of the MW, r O(1) kpc, AMCs are substantially depleted.At larger radii, r O(10) kpc, AMCs have a high probability of surviving.
The interactions of AMCs with stars can also alter the properties of the surviving AMCs (Fig. 8) through partial mass loss and energy injection.We have presented the procedure for deriving these properties using the results of our Monte Carlo simulations (with accompanying code and distributions at github.com/bradkav/axionminiclusters[76]).This allows our results to be reinterpreted for an arbitrary initial distribution of AMCs in the MW.
As discussed in § IV A, for computational simplicity we made a number of assumptions.We reiterate two assumptions here as they are of primary importance for future work.Firstly, we decoupled structure formation from our tidal interactions by initially populating the MW with AMCs following a halo mass function evaluated today.Fortunately, the majority of AMCs (in particular the less massive ones) form before the formation of the MW halo, meaning that we expect only small changes to our final distributions of miniclusters if we were to follow both structure formation and tidal interactions.This concurrent evolution should be investigated more precisely in future work.Unfortunately, it may be difficult to simultaneously obtain the necessary scale and resolution needed to simulate both stars and AMCs.Secondly, we assumed a static MW halo with a fixed stellar population as measured today (see Fig. 15).Future work should replace this assumption with a co-evolving stellar population that follows the cosmic star formation rate [110] and an evolving MW halo model such as Ref. [128].
Throughout the disruption calculations, we neglected the effect of axion stars (ASs), which may form in the centers of AMCs.The formation and evolution of these ASs is still uncertain.Since we do not account for these ASs, we instead place a cut on the AMC parameter space which requires the central AS to have a smaller radius than the host AMC (referred to as 'AS cut' and described in § II E).For those that pass the AS cut, we neglect the potential effects of a solitonic core on the stability of the minicluster.Future work should study the response of an AMC-AS system to tidal perturbations and extend our prescription to account for this more complete description of the axion sub-structure population.
Throughout this paper, we assumed that the entire population of AMCs in the MW had either a PL or NFW internal density profile.As discussed in § II D and § IV A, it is still unclear which density profile best describes the overall population of miniclusters -in fact it is likely that the high-mass miniclusters have more NFW-like profiles, while the lowest-mass miniclusters have PL-like profiles.Here, we have tried to bracket the uncertainty on the final distributions of miniclusters by using a very concentrated profile (PL) which is robust to perturbations and a more loosely bound profile (NFW) which is more easily disrupted.More work is needed to assess how the structure of these miniclusters evolves and, in particular, which AMC masses are better described by NFW or PL profiles (or perhaps some intermediate profile).By construction, our results are re-interpretable, allowing us to use this information to directly build more accurate descriptions of AMCs today when it becomes available.
In the post-inflationary scenario considered, the mass of the QCD axion has to be tuned to a specific value ma in order to reproduce the present DM abundance.If we were to set the mass of the QCD axion to be heavier than ma , the axion would not be the dominant component of DM in the Universe.This means that the predictions for the axion mass enclosed at t osc , i.e. the AMC mass, and the distribution of overdensities would have to be recomputed from new simulations taking into account the growth of axion overdensities around the dominant DM component.Such scenarios have never been considered in the literature, nor are they taken into account in the present work.On the other hand, an axion of mass m a < ma is not allowed, as its present energy density would be larger than what is observed for DM.For this reason, we have chosen m a = ma throughout the paper.Another subtlety lies in the possible range of ma .We have set ma = 20 µeV following recent estimates in the literature [57,75].However, uncertainties in the numerical computations hint at a range ma = O (10 − 100) µeV where the KSVZ axion mass could lie, see Ref. [35].Possible scaling violations in the string dynamics could also suggest that a large number of axions are produced from decaying strings in the early Universe, leading to a value of the axion mass ma ≈ 500 µeV that differs from the one we have set here [77].Assuming a different value of ma would lead to a different characteristic mass M 0 , which would cause the related expressions for M min and M max to be modified.However the analysis we have presented could be straightforwardly re-applied -as we have stressed -and the results would not change qualitatively.
Previous work has considered the effects of tidal interactions on AMCs [66,68].In particular, Ref. [66] considered AMCs with PL profiles ∝ r −1.8 and found that around 2-5% are destroyed at the Solar position.At the same position, we find that >99% of AMCs with a PL profile survive, while those with NFW profiles are more easily stripped with a survival probability ∼ 94% (see Fig. 9 without AS cut).This difference can easily be explained by the fact that our PL profile is significantly steeper than the one considered in Ref. [66], making our AMCs more robust to perturbations from stel-lar encounters.We also build upon Ref. [66] in two key ways: firstly, we extend the calculation of the properties of AMCs beyond the Solar position to the entire Galaxy; and secondly, we account for the injection of energy from non-fatal interactions, allowing us to more realistically describe disruption probabilities and the properties of AMCs today.
In our companion paper [59], we use the results of our Monte Carlo simulations to predict indirect signals from AMCs interacting with neutron stars (NSs).More precisely we: i) simulate the expected encounter rate for NSs passing through AMCs; and ii) estimate the expected radio signal from the conversion of axions within the minicluster into photons as they interact with the NS magnetosphere.These results are complementary to the continuous radio emission that is expected to come from axions in the Galactic halo interacting with NSs [38][39][40][41].
As we showed in Sec.VI, correctly calculating the AMC population today can have significant implications for a variety of observational channels.It is therefore of upmost importance to understand the interactions of these structures with their environment.Our results represent a fundamental step towards characterizing these interactions within the MW. on the mean AMC density), so we do not expect a mass-loss of O(10 − 40)% due to the host halo to substantially impact the survival probability of AMCs.

MFIG. 1 .
FIG.1.Models for the internal density profile of AMCs which we consider in this work: Power-law, Eq. (9), and NFW, Eq. (8).Vertical dashed lines show the truncation radii RAMC.We fix the characteristic mass and density to MAMC = 10 −10 M and ρAMC = 10 6 M pc −3 (δ ≈ 1.55) respectively.In the NFW case, the mean density is much lower and the AMC is much larger.

) 10 − 5 3 ∆EFIG. 3 .
FIG.3.Fractional mass loss (solid lines) for AMCs with Power-law and NFW profiles as a function of the size of the perturbation ∆E, in units of the binding energy E bind .Dashed lines show the fraction of the injected energy which is carried away by the ejected mass fej, while dotted lines show the fraction of the initial AMC energy stored in particles which are eventually unbound in the interaction.

3 MFIG. 4 .
FIG.4.Fractional change in the minicluster properties x = {M, R, E bind } under repeated perturbations for AMCs with Powerlaw profiles (left panel) and NFW profiles (right panel).We fix the size of each perturbation to be 1000 times smaller than the initial binding energy, ∆E = 10 −3 E bind,i .
r) approaches a delta function at r = a, i.e. a circular orbit.The values of P (r) are larger at the boundaries r = a(1 ± e) because the radial motion of the minicluster vanishes at either apsis and so dt/dr diverges.Figure5shows the joint probability distribution for a and e given a particular galactocentric radius (assuming the distribution for e shown in the top panel of Fig.5[103]).

1 P
FIG. 5.Top panel shows the eccentricity probability distribution for the orbits of AMCs taken from Ref.[103].Bottom panel shows the joint probability that an AMC will have a particular semi-major axis and eccentricity given a particular galactocentric radius r = 8 kpc.

9 FIG. 6 .
FIG.6.The binned probability of finding an AMC at a particular radius r with a fixed semi-major axis a = 1.We show distributions for eccentricity values: e = 0.1 (blue), e = 0.5 (orange), and e = 0.9 (green).

10 − 22 6 MFIG. 8 .
FIG.8.Example of AMC properties from our Monte Carlo simulations before (black dashed) and after (solid olive) disruption.We show the probability distributions of the mass MAMC, radius RAMC, and mean density ρ.We assume these to have NFW internal density profiles and to be on circular orbits with a galactocentric radius of r = 6.96 kpc, leading to a survival probability psurv = 0.91.These probability distributions are obtained through the reconstruction procedure described in § V A. The grey shading indicates the regions of the parameter space that are removed by the AS cut.The white lines shows the smallest value of the corresponding parameter that passes the AS cut.

FIG. 9 .
FIG.9.Survival probability of AMCs as a function of their galactocentric semi-major axes (left) and galactocentric radii (right).The vertical dashed line marks the position of the Solar System.In the panel we show the survival probability for both circular (dotted) and eccentric (dashed) orbits.The eccentric orbits are smeared out according to the proportion of time they spend at a given galactocentric radius (as described in § V A).Finally, we show the survival probability for miniclusters with eccentric orbits that pass our AS cut (solid line) where we have normalized this to be one at large galactocentric radii.

TABLE I
. Numerical values of the prefactors used for characterizing the AMC properties: the mean-squared radius R 2 ; binding energy E bind ; and velocity dispersion squared σ 2 .