Disruption of Dark Matter Minihaloes in the Milky Way environment: Implications for Axion Miniclusters and Early Matter Domination

Many theories of dark matter beyond the Weakly Interacting Massive Particles (WIMP) paradigm feature an enhanced matter power spectrum on sub-parsec scales, leading to the formation of dense dark matter minihaloes. Future local observations are promising to search for and constrain such substructures. The survival probability of these dense minihaloes in the Milky Way environment is crucial for interpreting local observations. In this work, we investigate two environmental effects: stellar disruption and (smooth) tidal disruption. These two mechanisms are studied using semi-analytic models and idealized N-body simulations. For stellar disruption, we perform a series of N-body simulations of isolated minihalo-star encounters to test and calibrate analytic models of stellar encounters before applying the model to the realistic Milky Way disk environment. For tidal disruption, we perform N-body simulations to confirm the effectiveness of the analytic treatment. Finally, we propose a framework to combine the hierarchical assembly and infall of minihaloes to the Milky Way with the late-time disruption mechanisms. We make predictions for the mass functions of minihaloes in the Milky Way. The survival fraction of dense dark matter minihaloes, e.g. for axion miniclusters and minihaloes from Early Matter Domination, is $\sim 60\%$ with the relatively low-mass, compact population surviving. The survival fraction is insensitive to the detailed model parameters. We discuss various implications of the framework and future direct detection prospects.


INTRODUCTION
The gravitational clustering of dark matter has been well measured on galactic scales and super-galactic scales and is consistent with a nearly scale-invariant spectrum of primordial fluctuations (e.g.Aghanim et al. 2020).However, the matter power spectrum on extremely small scales ( pc −1 ), which corresponds to sub-planetary-mass structures, is still weakly constrained and is sensitive to both the nature of dark matter and the thermal history of the early Universe.There have been proposals to detect small-scale structures in the mass range ∼ 10 −13 -10 2 M in the future with Pulsar Timing Arrays (PTAs; e.g., Siegel et al. 2007;Baghram et al. 2011;Dror et al. 2019;Ramani et al. 2020;Lee et al. 2021a,b) and lensing effects (e.g., Kolb & Tkachev 1996;Metcalf & Madau 2001;Diaz Rivero et al. 2018;Fairbairn et al. 2018;Katz et al. 2018;Van Tilburg et al. 2018;Dai & Miralda-Escudé 2020).
Many well-motivated dark matter theories can leave unique fingerprints on the primordial perturbations at small scales ( pc −1 ), such as the quantum chromodynamics (QCD) axion/axion-like particles (ALPs) with the Peccei-Quinn (PQ) symmetry (Peccei & Quinn 1977) broken after inflation (e.g., Hogan & Rees 1988;Kolb & Tkachev 1993;Kolb & Tkachev 1994;Zurek et al. 2007), Early Matter Domination (EMD; e.g., Erickcek & Sigurdson 2011;Fan et al. 2014) and vector dark matter produced during inflation (e.g., Nelson & Scholtz 2011;Graham et al. 2016).Therefore, small dark matter substructures provide unique insights into the microphysics of dark matter.The model space of interest here differs from that in more common WIMP-like collisionless cold dark matter (CDM) models.In those models, adiabatic fluctuations produced at the end of inflation can also seed small CDM We note that the "small scale" here is fundamentally different from the smallscale problem of CDM (at kpc scale) discussed in astrophysical studies (see review by Bullock & Boylan-Kolchin 2017).subhalos down to the kinetic decoupling and free-streaming limit ( ∼ O (pc −1 )), roughly corresponding to the Earth mass (∼ 10 −6 M ; e.g., Hofmann et al. 2001;Berezinsky et al. 2003;Green et al. 2005;Loeb & Zaldarriaga 2005).The evolution of these subhalos in the Milky Way environment has been studied in the past (e.g., Angus & Zhao 2007;Goerdt et al. 2007;Green & Goodwin 2007;Zhao et al. 2007;Schneider et al. 2010;Berezinsky et al. 2014;Delos 2019;Facchinetti et al. 2022).However, WIMP-like CDM formed its first non-linear structures at relatively late times with correspondingly low density (e.g. ∼ 60  eq as shown in Green et al. 2005), and the minihalos are thus subject to significant disruption due to tidal stripping and disk shocking after falling onto their host halos (e.g., Ostriker et al. 1972;Gnedin et al. 1999;Goerdt et al. 2007;Zhao et al. 2007;Schneider et al. 2010).The typical minihalos of WIMP-like CDM are out of reach for PTAs and other observations we discuss here (e.g., Lee et al. 2021a).
On the other hand, dark matter minihalos in the theories we consider here formed in the early Universe.Therefore, they are much denser and less likely to be disrupted by tidal forces than WIMP-like CDM subhalos.For instance, a (pseudo)scalar field (e.g. the QCD axion) with the PQ symmetry broken after inflation (e.g., Hogan & Rees 1988;Zurek et al. 2007, hereafter called the "post-inflationary axion") can induce order-unity isocurvature fluctuations on the horizon scale during the symmetry breaking.Regions with orderunity overdensities tend to collapse gravitationally very early, even before matter-radiation equality ( eq ∼ 3000), into small axion miniclusters (AMCs).The miniclusters underwent subsequent hierarchical clustering until the large-scale adiabatic perturbations intervened.As another example, EMD models can introduce a non-standard thermal history that is not constrained by any current data, but adiabatic fluctuations within the horizon can grow during this early period of matter domination.The characteristic mass of dark matter minihalos formed in early matter domination models is determined by the reheating temperature of this period.In vector dark matter models, the longitudinal modes of the vector DM produced at the end of inflation give rise to a peak in the matter power spectrum on small scales, with the scale directly determined by the dark matter particle mass (Graham et al. 2016;Lee et al. 2021a).Those models have interesting dynamics in the early Universe, which are difficult to probe directly.However, the remnants of those early Universe dynamics, dark matter minihalos, may be detectable in local observations (e.g., Dror et al. 2019;Ramani et al. 2020;Lee et al. 2021a).
While the evolution of (some versions of) these dark matter minihalos in the non-linear regime has been studied both semi-analytically with the Press-Schechter model (Zurek et al. 2007;Fairbairn et al. 2018;Enander et al. 2017;Blinov et al. 2020;Lee et al. 2021a;Blinov et al. 2021) and numerically with N-body simulations (Zurek et al. 2007;Buschmann et al. 2020;Eggemeier et al. 2020;Xiao et al. 2021), the gravitational interactions between dark matter minihalos and the large-scale dark matter structures or baryonic structures are not well studied, due to the large dynamic range and nonlinear behaviors involved.Dark matter minihalos formed from non-standard early universe dynamics can be as light as ∼ 10 −12 M , while the Milky Way has a halo mass of ∼ 10 12 M .Therefore it is challenging to resolve these small structures while simultaneously simulating the dynamics of the largest structures.However, it is critical to study the survival probability of minihalos in the Milky Way and to determine the prospects of detecting such structures in the local environment.This aspect is actively studied in the literature.For example, Kavanagh et al. (2020) made an analytic estimate of the effects of stellar disruption on AMCs as well as the survival fraction of miniclusters.Dandoy et al. (2022) proposed a quantum mechanical description of AMCs and studied the impact of stellar encounters using standard perturbation theories.However, these studies neglect some important non-linear effects, such as the realistic minihalo concentration distribution, multiple mechanisms of disruption, and their interplay.In this paper, we improve these estimates with detailed numerical simulations and semi-analytic treatment to combine all the expected dominant disruption terms.To deal with the enormous dynamic range involved, our strategy is to numerically simulate the dynamics of the dark matter minihalos in individual encounters (varying e.g.halo parameters and impact parameters), using these to build detailed semi-analytic models which can be used to treat the largescale behavior.This allows us to capture the key non-linear physics on small scales while making realistic predictions for the overall behavior of dark matter minihalos in the Milky Way galaxy (in ways that can be generalized, in principle, to a broad class of minihalo-like models).Broadly speaking, the disruption of dark matter minihalos in the Milky Way can be divided into two parts: stellar disruption and tidal disruption.The stellar disruption term is sensitive to close encounters with individual stars, while the tidal term depends on the gradient of the collective gravitational potential on large scales.We study these separately using N-body simulations and then combine them using semi-analytic models.
This paper is organized as follows.In Section 2, we discuss the minihalo mass function, mass-concentration relation, and the analytic models for stellar disruptions due to individual stellar encounters and tidal disruptions.In Section 3, we present our numerical results from a series of idealized Nbody simulations for isolated minihalos encountering a star and simulations of tidal stripping.In Section 5, we apply the semi-analytic model calibrated using idealized simulations to a realistic Milky Way environment and combine both stellar disruption and tidal disruption.In Section 6, we study the survival fraction of dark matter minihalos by applying the framework to different physical models.
We assume a ΛCDM cosmology with parameters given as ℎ = 0.697, Ω m = 0.2814, and Ω Λ = 0.7186, and adopt scalar spectral index   = 0.9667.These are consistent with the recent Planck results (Aghanim et al. 2020), and our conclusions are relatively insensitive to variations in these parameters (compared to the much larger uncertainties in e.g.minihalo properties from different physical models).

Initial mass functions and concentrations of small-scale structures
The two models discussed in the introduction that lead to the formation of dense minihalos in the early Universe are: i. (pseudo)scalars (e.g. the QCD axion) with symmetry broken after inflation (Hogan & Rees 1988;Zurek et al. 2007), and ii. the EMD model (Erickcek & Sigurdson 2011).These two models are physically well-motivated and are also representatives of models enhancing the matter power spectrum at small scales.Although the cosmological perturbations in those models have different origins, they share common features.They are generated at extremely small scales before matter-radiation equality, forming dark matter substructures as light as 10 −12 M ubiquitously, decoupled from the usual adiabatic fluctuations.The model parameters are the axion mass  a for the AMC model and the reheating temperature  rh for the EMD model.We vary the axion mass from 1.25 eV to 500 eV, which is roughly the mass window that can produce the correct dark matter relic abundance, given uncertainties introduced by the axion emission from strings as shown in Buschmann et al. (2022); Gorghetto et al. (2021).On the contrary, the reheating temperature is loosely constrained although  rh cannot be below a few MeV.Otherwise, Big Bang Nucleosynthesis will be spoiled (see e.g.Kawasaki et al. 1999).We study the value from 15 MeV to 120 MeV.We choose the fiducial values  a = 25eV and  rh = 30 MeV, such that the minihalo mass function in both models peaks at ∼ 10 −8 -10 −7 M .
The formation of dense substructures from these smallscale perturbations occurs at early times (see Lee et al. 2021a for an analytic study and Xiao et al. 2021 for simulations of axion miniclusters).The formation and hierarchical mergers of minihalos are described by the redshift-dependent mass function d 0 /d (), calibrated by simulations in Xiao et al. (2021).Meanwhile, the minihalos can fall onto CDM substructures formed from adiabatic perturbations.The redshift of infall,  i , which occurs after the minihalo formation, defines the time when they stop merging or accretion, owing to the high virial velocities in the normal CDM halos, and their properties (mass, concentration) remain unchanged afterward.In general, the halo concentration hits a floor at the massive end as the adiabatic CDM power spectrum takes over.At the low-mass end, since one would not expect minihalos to form before matter-radiation equality, a cap to the concentration appears at around (1 Therefore, the final mass function, including minihalos falling onto the CDM structures, can be expressed as where  CDM col () is the collapse fraction of normal CDM halos formed from adiabatic fluctuations.The probability of infall , where we have assumed there is no dynamic decoupling between collapsed minihalos and the overall dark matter content.
We can use the Press-Schechter model (Press & Schechter 1974) to compute the collapse fraction where  c = 1.686 is the critical overdensity for spherical collapse,  () is the growth function,  2 CDM () is the variance of the adiabatic fluctuations calculated using the Code for Anisotropies in the Microwave Background (C ; Lewis et al. 2000;Howlett et al. 2012;Lewis 2021), and  min is the smallest CDM structure we consider formed from adiabatic fluctuations.We take  min to be 10 −2 M , corresponding to a scale where the CDM power spectrum dominates while the enhanced small-scale perturbations become subdominant.The result is insensitive to  min as () depends logarithmically on  at small scales.There could be even smaller CDM halos that can form, which will increase the collapse fraction  CDM col (), especially at high redshifts.However, below a certain mass, the normal CDM halos may have comparable masses to the enhanced substructures, and our assumptions will no longer hold.A detailed study of this infall process will ultimately be required for more detailed predictions and we leave it for future study.
The method described above allows us to compute the collapse fraction at various redshifts analytically and obtain the final mass function for the enhanced substructures.This will also give us the distribution of the infall redshift  i , P ( i ) , which determines the time when the structural changes of minihalos halt and thus the average density of minihalos.Including minihalos falling onto CDM halos, we obtain the initial (pre-disruption) minihalo mass function at  = 0 from different models, as shown in Figure 1.
Regarding the concentration of minihalos, we adopt the model in Lee et al. (2021a) to calculate the concentration from the power spectra of different models.For each mass, it evaluates the redshift when the corresponding primordial fluctuation mode collapse ( c ) and assumes () ∝ (1 +  c ()), which is similar to the CDM results (e.g.Bullock et al. 2001).The only difference here is that we assume the merger or smooth accretion of minihalos stops at  =  i instead of  = 0. Therefore, effectively, we have () ∝ (1 +  c ())/(1 +  i ).
In Figure 1, we show the initial (pre-disruption) mass func-C documentation Integrated over the entire mass range, P ( i ) ∝ d  CDM col ()/d is purely determined by CDM cosmology.However, if looking at minihalos at  = 0 in a certain mass range, the conditional probability distribution of  i will be mass dependent.In the hierarchical formation of minihalos, more massive ones form at later times.Therefore, more massive minihalos found at  = 0 tend to have lower  i .
tions and the mass-concentration relations of minihalos in different physics models.At the massive end, the adiabatic CDM power spectrum will dominate and one should reproduce the mass-concentration relation of normal CDM halos.There is a cap for concentration at the low-mass end since one will not expect minihalo formation before matter-radiation equality.The axion mass for the AMC models or the reheating temperature of the EMD models only creates constant mass shifts in the mass function and the mass-concentration relation.For comparison, the typical concentration of WIMP-like CDM halos studied in e.g. Green & Goodwin (2007) and Delos (2019) are  ∼ 1 -20 at the redshift when minihalos are initialized (usually  ∼ 60−100 in these studies).In terms of the central density, they are comparable to the relatively massive minihalos studied here.

Disruption in late-time evolution
Owing to the ultra-compact structure of the minihalos, they are largely immune to external perturbations through their evolutionary history after decoupling from the Hubble flow.However, after falling into a massive host system like the Milky Way, non-linear gravitational interactions with the host halo and the dense baryonic structures in the host halo could lead to significant disruption of minihalos.The two leading disruption mechanisms are tidal disruption from the host halo (and the baryonic disk) and close encounters with stars (referred to as stellar disruption).The relevant spatial and time scales on which these two mechanisms operate are drastically different.In this section, we will review the analytic model developed in the literature.We also note that the term "disruption" used in this paper describes the mass loss of minihalo from external perturbations at various levels and is not restricted to the case where a minihalo is completely "destroyed" as in some literature.

Stellar disruption
First, we consider the consequence of the encounter between a minihalo and a star.The virial radius of minihalos with a mass of interest (∼ 10 −10 M ) is of the order of 0.01 pc ∼ 2000 AU, which is still much larger than the radius of main sequence stars even though we are considering minihalos.Therefore, for simplicity, stars can be treated as point-like objects during encounters.In addition, after a stellar encounter, the structure of the minihalo cannot immediately react to the energy imparted during the encounter.It takes roughly a dynamical time for the minihalos to relax to the final state after disruption, which is given by where ρmh is the average density of the minihalo (assumed to be Δ c times the critical density of the Universe at  i , see the discussion after Equation 11).The dynamical time  dyn is comparable to the Hubble time scale for  i = 0, but is much shorter than that if the minihalos fall into CDM structures at a high redshift.The duration of the star-minihalo encounter, however, is several orders of magnitude shorter than  dyn .Therefore, the impulse approximation holds and the encounter can be treated as an instantaneous interaction (Spitzer 1958).
In the distant-tide approximation (when the impact parameter is much larger than minihalo size), the imparted energy from a single star encounter can be expressed as (Spitzer 1958) where Δ is the increase of internal energy of the minihalo,  mh ( mh ) is the virial mass (radius) of the minihalo,  2 represents the mean-squared radius of particles with respect to the center of the minihalo. * , , and  * are the mass, impact parameter, and relative velocity of the stellar object, respectively.The mean-squared radius can be parameterized as  2 =  2  2 mh , where  is a dimensionless parameter determined by the density profile () of the dark matter minihalo Assuming that the minihalo has the NFW profile (Navarro et al. 1996(Navarro et al. , 1997)), one obtains where  is the concentration number of the minihalo.However, when the impact parameter becomes comparable to the size of the minihalo, the distant-tide approximation made by Equation 4 breaks down.The strong  −4 dependence of Δ will disappear once the star passes through the minihalo and the disruption effect is suppressed (e.g., Gerhard & Fall 1983;Moore 1993;Carr & Sakellariadou 1999;Green & Goodwin 2007).For a single encounter, Green & Goodwin (2007) proposed a more general treatment of the imparted energy calibrated using simulations where  s =  b (2/3) 1/2  mh is the transition radius, which is close to the physical size of the minihalo up to a factor determined by structure parameters, and  b is an order-unity correction factor we introduce to be determined by our simulations, which will be discussed in Section 3.  is another structural parameter defined as where the NFW profile is assumed and  c is the smallest radius that the profile extends to, which we assume to be 0.01  s .In principle, an axion star formed in the center of a minihalo can provide a natural physical scale for  c .However, the size of an axion star sensitively depends on particle physics parameters as well as the growth rate of axion stars (e.g., Visinelli et al. 2018;Helfer et al. 2017;Chen et al. 2021).
We note that the choice of  c has weak effects since (1) it only appears in the logarithm; (2) the uncertainties of  c and  have been effectively captured by the free parameter  b ; (3)  in fact only matters for rare close encounters.Compared to full analytic calculations under the impulse and distanttide approximations, Equation 7gives better agreement with simulations in the transitional regime.Roughly speaking, disruption of a minihalo is expected to occur when the increase in internal energy of the minihalo given by Equation 4 exceeds the binding energy of the minihalo Here  is a dimensionless parameter again determined by the mass profile of the dark matter halo.For the NFW profile, it takes the form (Mo et al. 1998) Utilizing Equation 7and 9, we obtain the normalized energy input to the minihalo as ) where ρmh ≡  mh /(4  3 mh /3) is the average density of the minihalo.Assuming the minihalos are in virial equilibrium with respect to the background density at redshift  i (the infall The original formula proposed in Green & Goodwin (2007), who were studying low concentration halos, does not have the correction term (equivalently  b = 1).However, minihalos with higher concentrations are studied in this work and we find the correction term is necessary empirically to fit the simulation results.redshift), we obtain ρmh = Δ c  crit ( i ), where  crit ( i ) is the critical density of the Universe at  i and Δ c = 200 is the critical overdensity of collapsed objects (neglecting corrections from e.g.Ω Λ ,  crit ( i ) ∼ (1 +  i ) 3 ).Since ρmh is independent of minihalo mass, the energy input in the large- case in Equation 11 will be independent of minihalo mass while having a strong dependence on the impact parameter.The only free parameter that is left to be determined by simulations in Equation 11 is the correction factor  b .
From Equation 11(assuming the  >  s case), we can compute a characteristic impact parameter when Δ/ b = 1 which gives a crude estimate of the condition for the destruction of minihalos.It is worth noting that the average density of a virialized halo is much larger when formed at higher redshifts, given as ρmh ∝ (1 +  i ) 3 .Minihalos that collapsed and fell into the host halo earlier in cosmic time should be less vulnerable to stellar disruptions due to higher central densities.Although  min serves as an indicator for "significant" disruption of the minihalo, the actual mass loss of the minihalo (as a function of Δ/ b ) after a single encounter should be calibrated by numerical simulations presented in the following section.

Tidal disruption
Another important disruption mechanism for minihalos is tidal disruption.After the minihalos fall into the Milky Way halo, they experience tidal forces from the Milky Way dark matter halo and the Galactic disk.In contrast to stellar disruption, which sensitively depends on close encounters with stars, the tidal forces are determined by gravity at galactic scales, dominated by the collective effects of the smooth gravitational potential rather than the fluctuating component from individual objects.Moreover, rather than an impulsive event, tidal disruption is secular and usually treated differently from stellar disruptions.
In a sufficiently strong and smooth external tidal field, the outskirts of the minihalo will be stripped away where the tidal force exceeds the self-gravity of the minihalo.This tidal radius is given by (e.g., King 1962;Taylor & Babul 2001;Zentner & Bullock 2003a) where  is the galactocentric distance of the minihalo,  MW ( < ) is the mass of the Milky Way (including both dark and baryonic matter) enclosed within radius  and  mh ( <  t ) is the minihalo mass enclosed within   . t (R) is the instantaneous tangential speed of the minihalo, which equals the circular velocity  c () when the minihalo is on a circular orbit.The minihalo mass outside the tidal radius will be stripped away roughly over a dynamical time scale of the host, and the instantaneous mass-loss rate is often expressed as (e.g., Taffoni et al. 2003;Zentner & Bullock 2003b;Oguri & Lee 2004;Pullen et al. 2014) where  ts is the characteristic time scale of tidal stripping, which is proportional to the dynamical time scale of the host halo (not the minihalo) where ρhost () is the averaged density of the host halo within radius  and  is a constant fudge factor of order unity that is found to be ∼ 0.5 -3 in several previous studies (e.g., Zentner & Bullock 2003b;Zentner et al. 2005;Pullen et al. 2014;van den Bosch et al. 2018;Green et al. 2021;Errani et al. 2023).If the tidal stripping occurs in the Solar neighborhood, the time scale will be much shorter than the lifetime of minihalo in the Milky Way environment.As a simple estimation, the characteristic formation time (defined as when d log /d log  falls below a threshold as proposed in Wechsler et al. 2002) of a Milky Way-mass halo is about  = 0.3 (corresponding to  ∼ 2.5 and lookback time 10 Gyr).If the minihalos are dynamically coupled to the smooth dark matter content accreted by the Milky Way halo, the typical lifetime of minihalos in the Milky Way environment will be of the same order.Meanwhile, since  ts is approximately the orbital period of the minihalo at the pericenter, a reasonable assumption is that the minihalo outskirts will be tidally disrupted "immediately" in the first pericenter passage and before any form of stellar disruptions takes place.

IDEALIZED SIMULATIONS FOR STELLAR ENCOUNTERS AND TIDAL STRIPPING
In this section, we use N-body simulations to systematically study the stellar disruption and tidal stripping effects for isolated minihalos initialized with the NFW profile.The goal is to test the analytic models described in Section 2 and calibrate them against various minihalo parameters.We show the dark matter surface density distribution of the minihalo.The left panel shows the initial dark matter distribution.Right after the encounter, a large fraction of dark matter has been heated up by the gravitational interaction with the passing star and become unbounded.Although this energy is transferred impulsively, it requires roughly the minihalo dynamical time ( dyn 1 Gyr) for the disruption to be reflected in the minihalo density distribution, as shown in the middle and right panel.
We perform a suite of N-body simulations for minihalos initialized with the NFW profile by varying minihalo concentrations, masses, and impact parameters of the encountering stars.The simulations adopt the code G (Hopkins 2015), which has been widely used in cosmological N-body or hydrodynamical simulations (e.g., Hopkins et al. 2014Hopkins et al. , 2018;;Feldmann et al. 2022).The simulations aim to test and calibrate analytic predictions of the minihalo mass fraction disrupted in stellar encounters.The relaxation process of minihalos and the evolution of density profiles are also studied, which provides important insights into modeling the disruption effects under multiple encounters.
In these simulations, an isolated minihalo (composed of collisionless N-body dark matter particles) is initialized at  = 0 with a star (represented by a single point-mass particle) having mass  * = 1 M at a large distance (10 pc) moving towards the minihalo with a relative velocity of 200 km s −1 .Since the code is Lagrangian, it makes no difference whether the star or the minihalo is moving and in which frame we solve the dynamics equations.In Figure 2, the evolution of a minihalo during disruption is visualized.The encounter with a star will impart a certain amount of energy into the minihalo, disrupting the halo outskirts at roughly a minihalo dynamical time.Our default simulations resolve the minihalo with 10 6 dark matter particles initialized in an equilibrium NFW G documentation The simulations are not cosmological, but the redshift value is required for initializing the NFW halo.The response of minihalos to encounters we test in this section is not sensitive to the density normalization or the redshift we set up the minihalo.halo following the method in Springel (2005).We generate the initial condition using the package IC (Herpich et al. 2017).To systematically study stellar disruption effects, we vary the impact parameter of the star, the minihalo mass, and the minihalo concentration (or equivalently scale radius).The gravitational softening length of the dark matter particles is taken to be 10 −9 kpc, which is small enough to resolve the dense core of a minihalo with mass 10 −10 M , concentration  = 100, and scale radius  s = 9.6 × 10 −8 kpc.The timestepping in the simulation requires more careful consideration since most disruption occurs when the distance between the star and the minihalo is around its minimum.We choose the maximum size of the timestep to be 10 −8 Gyr even when the star is still at a large distance (10 pc) from the minihalo and study the subsequent disruption as the star moves closer towards the minihalo.It is worth noting that this timestep is about an order of magnitude smaller than the crossing time of the star, ∼  mh / * so that the trajectory of the star around the minihalo can be well resolved.After the close encounter, we can relax the upper limit on the timestep to integrate the relaxation of the minihalo after the stellar disruption to arbitrarily long times (at a fairly low computational cost).A detailed discussion of numerical convergence is given in Appendix B. Experimenting with the maximum timestep and/or the gravitational softening length shows that smaller values produce essentially identical results (with larger computational costs).However, order-of-magnitude larger values of softening length can risk allowing the simulation to "over smooth" gravity or take excessively large timesteps for some particles, which "overshoot" the very brief duration of the encounter (for detailed numerical tests, see Hopkins et al. 2018 andGrudić &Hopkins 2020).

Disruption under different impact parameters and calibration of the response function
We first run three sets of simulations for encounters with fixed minihalo mass and concentration but different impact parameters.In all these simulations, we set a fixed minihalo mass  mh = 10 −10 M , star mass  * = 1 M and initial star-minihalo relative velocity  * = 200 km/s.The halo concentration has been fixed to be  = 10, 30, 100 for each set of simulations, respectively.The goal is to characterize the relation between the mass loss of the minihalo after a stellar encounter and the normalized energy input.The imparted energy can be related to the impact parameter with Equation 11.After the minihalo becomes fully relaxed after the stellar encounter (at  = 1 Gyr >  dyn ( = 0)), dark matter particles, with kinetic energy (in the center-of-momentum frame) larger than the absolute value of their gravitational potential energy, are identified as unbound and disrupted.The remaining mass of the minihalo is measured as the fraction of bound dark matter particles after minihalo relaxation, and we have verified that the remaining minihalo mass has converged at the time of measurement.
In Figure 3, we show the minihalo mass loss as a function of the normalized imparted energy for different choices of halo concentrations.The imparted energy is calculated using Equation 11.The mass loss is negligible when Δ/ b 1, and quickly increases in a power-law fashion with respect to Δ/ b .In this regime, minihalos with low concentrations are more vulnerable to stellar encounters with steeper powerlaw slopes.In the following, the mass loss as a function of the imparted energy and halo concentration will be referred to as the response curve, F (Δ/ b , ).We propose the following functional form to fit the response curve where Δ/ b should be evaluated using Equation 11(the free order-unity parameter  b is yet to be determined, but the calibration here is done in the  >  s regime, so  b does not have any real impact).After exploring the response curve of each choice of minihalo concentration, we propose the following fitting formula for the parameters () and  () Then we perform the least-square fits jointly on the results of all sets of simulations.The best-fit parameters are  = 0.987,  1 = −0.8(fixed), 2 = −0.586, 3 = −0.034, 0 = −0.583, 1 = −0.559.The best-fit model is also shown in Figure 3 and is in good agreement with the simulation results of minihalos with various concentrations.This best-fit model of the response curve will be the foundation we use to understand the disruption of minihalos with different masses or concentrations.For comparison, we also show the analytical model in Kavanagh et al. (2020) (where  = 100 is fixed) and the semi-analytical model developed in Delos (2019) (which is calibrated for large Δ/ b values).Our model agrees reasonably well with their results for shared dynamical ranges and minihalo parameters.

Disruption under different halo concentrations
Given the calibrated response function, we next run an additional set of simulations with a fixed impact parameter  = 0.05 pc but for minihalos with different concentrations.The goal is to further validate this semi-analytic model (the analytic calculation of imparted energy plus simulation calibrated response function) for minihalos with various compactness.The minihalo mass is still fixed to  mh = 10 −10 M and the properties of the encountering star are the same as in Section 3.1.1.The halo concentration will affect the mass loss from stellar encounters in two ways.First, the structural parameters , , and  all have an explicit dependence on minihalo concentration, which will propagate to the calculation of energy imparted (i.e.Δ/ b ∝  2 ()/() ∼ ln ()/ when   s and  1).Secondly, the response function also has a strong dependence on concentration (see Equation 16), especially when Δ/ b 1. Less concentrated minihalos will become increasingly vulnerable to disruptive stellar encounters.Therefore, it is non-trivial for the semi-analytic model to correctly capture the concentration dependence of minihalo mass loss.In Figure 4, we show the mass loss of minihalo versus minihalo concentration.The best-fit semianalytic model agrees with the simulation results over a wide range of concentrations, even extrapolating into regimes not covered in our original calibration step above.

Disruption under different minihalo masses
One remaining ingredient of the semi-analytic model that needs to be calibrated is the disruption behavior in the   s regime.According to Equation 11, the imparted energy will stop rising as  decreases when a characteristic scale  s is reached.The free, order-unity correction factor should be calibrated by simulations.To fulfill that, we run an additional set of simulations for minihalos of different masses ranging from 10 −10 to 10 −3 M but fixing the impact parameter  = 0.05pc and minihalo concentration  = 100.In low-mass minihalos,  s is small enough that the encounter is in the  >  s regime, where mass loss is independent of minihalo mass.As minihalo mass increases (and  s increases), test cases with massive minihalos enter the  <  s regime, where stellar disruptions are suppressed.In Figure 5, we show the mass loss as a function of minihalo mass.A transition of the response of the minihalo occurs, and we use the mass of this transition to calibrate the free parameter to be  b = 6.The predicted mass loss using the best-fit response curve for  = 100 halos is shown with the solid line in the figure .In general, the semi-analytic model gives the correct location and shape of the transition at  ∼  s .The mass loss remains a constant at low minihalo masses as indicated by Equation 11in the   s regime.A sharp transition occurs at   s such that the mass-loss rate reduces to almost zero at high minihalo masses.

Density profiles
The density profile will not immediately change after the close encounter with stars, but will gradually relax to the final minihalo profile within a few minihalo dynamical times (e.g.Delos 2019).This is mainly because the close encounter with a star occurs on a timescale much shorter than the minihalo dynamical time.
In Figure 6, we show the evolution of minihalo density profiles for minihalos with mass 10 −10 M and concentration  = 100.In the top (bottom) panel, we show the case with impact parameter  = 2 × 10 −5 pc (5 × 10 −5 pc).The density profile is shown as   3 , which is proportional to Δ/Δ log , the contribution to total mass per unit logarithmic interval of radius.The vertical shaded region shows the convergence radius for collisionless particles based on the Power et al. (2003) criterion (this is roughly the radius interior to which the numerical two-body relaxation time drops below the Hubble time, a conservative indication of where N-body integration error could be significant).In the  = 5 × 10 −5 pc case, the outskirts of the halo are predominantly disrupted while the core of the halo remains relatively unperturbed, leaving approximately 60% of the minihalo mass disrupted.Although the central density exhibits a small decrease, its contribution to the total mass loss is negligible and the scale of this decrease is close to the convergence radius of dark matter properties, which makes it hard to distinguish the decrease from a numerical artifact.At the outskirts of the minihalo, the density profile turns up, corresponding to a shell of unbound particles (heated by the encounter) propagating outward.For the  = 2 × 10 −5 pc case, due to higher imparted energy, the disruption is more significant with over 80% of the minihalo mass disrupted.However, the behavior of the final density profile is rather similar to the previous case.The density profile of minihalo after disruption can be well-fitted by a broken power-law profile with best-fit asymptotic slopes around −4 analogous to that of the Hernquist profile.The "model fit" curve provides an asymptotic limit of the density profile after infinite time.
The upturn at the outskirt of the minihalo represents a propagating shell of unbound particles.As this upturn propagates towards the outer part of the minihalo, the density profile gradually converges to the asymptotic limit.
It is worth noting the remaining minihalos never relax to a new NFW profile.Instead, the final density profile (after the shell of unbound particles escapes) can be well described by a broken power law of the form with the best-fit asymptotic slope  = 3.2 for  = 2×10 −5 kpc and  = 3.3 for  = 5 × 10 −5 kpc.This slope implies that the density profile is close to a Hernquist (1990) profile ( = 3), or in a more general sense, the -profile family (Dehnen 1993;Tremaine et al. 1994) with an asymptotic slope of −4.Unsurprisingly, given the scale-free physics involved, this is similar to well-studied simulations of impulsive high-mass ratio galaxy-galaxy encounters (e.g., Hernquist & Quinn 1988;Barnes & Hernquist 1992;Hopkins et al. 2008Hopkins et al. , 2009;;Boylan-Kolchin et al. 2005), particularly the structure of shell galaxies (Hernquist & Quinn 1987;Hernquist & Spergel 1992) and the closely-related slopes of the marginally-bound layer of particles at apocenter in cosmological simulations defined as the halo "splashback" radius (e.g., Diemer & Kravtsov 2014;More et al. 2015).Generically, this slope arises from relaxation after dynamical mass ejection events, which (by definition) excite some material to a broad distribution of energies crossing the specific binding energy E = 0, as the outer slope  +1 = 4 corresponds to the only finite-mass asymptotic power-law distribution function which is continuous through E = 0 (Hernquist et al. 1993).

Numerical test of tidal stripping
To validate the semi-analytic description of tidal disruptions of minihalos in Section 2.2.2, we perform an idealized simulation of a minihalo traveling through the analytic gravitational potential of a Milky Way-mass host halo ( vir = 10 12 M ).The host halo profile is modeled as an NFW profile with concentration  = 12 (e.g., Klypin et al. 2002;McMillan 2011;Deason et al. 2012;Bland-Hawthorn & Gerhard 2016).The minihalo is assumed on a circular orbit with  =  8 kpc.In Figure 7, we show the evolution of the enclosed mass profile and density profile of the minihalo.The tidal radius calculated using Equation 13 is shown as the red vertical dashed line.After evolving for about 0.5 Gyr (for reference, the dynamical time scale of the host halo at  = 8 kpc is  host dyn ∼ 0.1 Gyr, see Equation 15), the enclosed mass profile starts flattening and eventually plateaus outside the analytically evaluated tidal radius.In Figure 8, we show the mass evolutionary history of this minihalo, and specifically the mass inside and outside the tidal radius.The mass within the tidal radius is almost immune to tidal stripping, with 75% of the mass remaining after 0.5 Gyr.On the other hand, the mass outside the tidal radius exhibits an exponential decay after about 1  host dyn .The mass loss of the minihalo can be well represented by Equation 14, implying  ( >  t ) ∼  −/  host dyn with a fudge factor  close to 1.8 for the minihalo tested, in broad agreement with previous numerical studies of more massive CDM subhalos (e.g., Zentner & Bullock 2003b;Pullen et al. 2014;van den Bosch et al. 2018) The simulation presented in this section demonstrates that the semi-analytic treatment of tidal disruption works reason- ably well in predicting the mass loss and the tidal stripping time scale.Given the order unity fudge factor  found here, the tidal stripping time scale  ts in Equation 15 is at least an order of magnitude smaller than the Hubble time.Therefore, tidal stripping can be treated as an "instantaneous" process in modeling the cosmological evolution of minihalos.

MODELING MULTIPLE ENCOUNTERS
In the preceding section, we introduced a model that describes the mass loss of minihalo resulting from individual stellar interactions.In this section, we will present a model that deals with the cumulative effects of multiple encounters.(the same one as in Figure 7).The total mass here is the sum of dark matter mass within the initial virial radius of the minihalo.In addition, we show the mass inside and outside the tidal radius.After approximately a dynamical time of the host halo, the mass outside the tidal radius starts an exponential decay as described by Equation 14with the fudge factor  1.8.The mass inside the tidal radius is marginally affected.
To begin, we examine the amalgamation of encounters occurring within the dynamical time of the minihalo.A relevant illustration of such an event would be a minihalo traversing the Milky Way disk.The typical crossing time (given by Equation 22) is significantly shorter than the dynamical time of the minihalo (given by Equation 3).Minihalo does not have enough time to fully relax and realize the mass loss before the subsequent encounter occurs.
The change in the velocity of a particle within a minihalo of size  at position  (relative to the center of the minihalo) due to an encounter with a perturber of mass  * moving with relative velocity v * at an impact parameter b (when  ) can be expressed as (Spitzer 1958;Green & Goodwin 2007) where e b is the unit vector in the direction of b, e z is the unit vector perpendicular to b and u (e z ≡ e b × e u ).For a particle in the minihalo moving with velocity , the change in energy is given by The ensemble average of  over all particles in the minihalo can be obtained in two steps.First, at a fixed position r, an average over velocities v is taken, where the first term in Equation 20 vanishes if the particle velocities are assumed to be isotropic.Second, an average over r is taken, where the assumption of spherical symmetry leads to the final form of Δ in Equation 4. When multiple encounters occur in a short enough time compared to the dynamical time of the minihalo, particles in the minihalo can be considered frozen.The aggregation of encounters is simply v = Σv i .According to Equation 20, the aggregated  will have two types of terms: v i • v i and v i • v j , where the former is simply energy injection from each individual encounters while the latter involves two different encounters When the total number of encounters aggregated is large, we can take the ensemble average over encounter parameters of the i-th and j-th events.The result depends on the local velocity distribution of stars with respect to the minihalo.Two typical scenarios are considered in the following.
(1) An isotropic velocity field: Both e b and e z are isotropic.
Assuming the i-th and j-th encounters are independent, the ensemble averages (over the orientation of incident stars) of the four terms in Equation 21 are the same and exactly cancel each other.(2) A stream of stars: All the stars are moving in the same direction.However, the impact parameter b is still isotropic in the cross-section of the stream.Therefore, e b and e z are isotropic in this two-dimensional plane.The ensemble averages of the four terms in Equation 21 are the same and cancel each other.For a realistic application, if we consider a minihalo moving through a cloud of stars with an isotropic velocity distribution, the relative velocity field of stars with respect to the minihalo is a linear combination of the two scenarios above.Due to the linear nature of the ensemble average, we can conclude that the v i • v j terms can be ignored.The v i • v i terms are left and  tot is simply the summation of  from individual encounters.

Multiple encounters after full relaxation
A more complex scenario arises when the minihalo has enough time to fully relax before the next encounter, such as multiple disk crossings with an orbital period comparable to or longer than the relaxation time.In such cases, the disruption model calibrated on single encounters needs to be revised because the internal structure of the minihalo can differ significantly from the initial condition.In Section 3, we demonstrated that the density profile post-disruption is no longer NFW-like and features steeper slopes at the outskirts.
One approach often used in the literature is to integrate the effective "optical" depth.This relies on the approximation that minihalos are completely destroyed when the energy imparted exceeds a certain threshold.The number fraction of minihalos "survived" scales exponentially with the "optical" depth of encounters (see the definition in Equation 25).This approach has been used in e.g.Schneider et al. (2010) for CDM minihalos.We will discuss it more and use it as a back-of-the-envelope estimate in the following section.The second approach is to "resolve" individual encounters and process them sequentially, relying on a model to describe the internal states of minihalos after the encounter(s).Delos (2019) proposed a self-similar profile for post-encounter minihalos and calibrated disruption models accordingly.This method works reasonably well in the large d/ b regime, as shown in the comparison in Figure 3.However, it is unknown whether it is accurate in the limit we are interested in.The typical d/ b for a single stellar disk crossing is of order 10 −2 (see Equation 31), which is outside the calibration range of this model .
Here we propose an alternative and intuitive way to treat multiple encounters which directly builds upon the knowledge about single encounters.This method requires no assumptions on parameterizing the density profiles of postdisruption minihalos and works reasonably well in a quasistatic limit.Similar to what we derived for multiple encounters within a dynamical time, we accumulate the injected energy Δ for successive encounters and use the response function F (Δ tot / b ) to obtain the total mass disruption fraction.The validity of this treatment is supported by the following observations.
• The state of a minihalo after disruption could be well described by a single parameter.This parameter can be the total energy or the gravitational binding energy It should be noted that the definition of d/ b in Delos ( 2019) differs from the one used in this paper, as they normalize d by the binding energy of the central dark matter core.
or some characteristic density, etc.An example of such a description is given by Delos (2019).They found a self-similar density profile for minihalos after disruption with two parameters  s and  s .A scaling relation between the two parameters was also discovered.Therefore, the density profile of the minihalo is effectively determined by one single parameter.This parameter can be the total energy of the minihalo because the total energy is uniquely related to the density profile assuming virial equilibrium.
• For the typical d/ b ∼ 10 −2 of a single stellar disk crossing, the disruption happens in a quasi-static fashion.In this limit, each encounter only causes a marginal disruption to the least bound particles.The change in the total energy of the minihalo can be expressed as where  ej is the fraction of imparted energy ends up in unbound particles.In the quasi-static limit, the shell of unbound particles that eventually escape the minihalo will have negligible energy at infinity,  ej Δ +  init unbound → 0. This is supported by both our simulations and the analytical calculations shown in Kavanagh et al. (2020).For example, when Δ/ b ∼ 10 −2 , the corresponding  ej is ∼ 10 −2 and  init unbound / i is ∼ 10 −4 .Neglecting second order terms, we have  final bound −  init bound Δ.
Combining the two points listed above, the final state of the minihalo is determined by the final energy of the minihalo, which is only related to the total imparted energy Δ and has no dependence on the exact history of encounters.Using this approach, the additional disruption from successive encounters is sensitive to the slope of the response function at large Δ tot / b .A unique feature of the calibrated response curve is that subsequent encounters with the same energy injection will rapidly become less important (Δ mh / mh from each encounter rapidly decreases as the total number of encounters increases), especially for minihalos with large concentrations.The physical interpretation is that the particles which would be unbound by such an energetic encounter have already been unbound in previous encounters (the loosely bound particles at the outskirt have already been stripped and the remaining particles are more tightly bound and less vulnerable to new encounters).
To validate this treatment, We perform a series of numerical tests.We take a minihalo with  mh = 10 −10 M ,  = 100 and set the mass and velocity of the star as 1 M and 200 km/s.In the first test, we pick an impact parameter such that the injected energy from a single event is Δ/ b ∼ 0.1.We simulate eight repetitive encounters with Δ/ b = 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2.The minihalo is allowed to fully relax between these encounters.In Figure 9, we compare the mass loss curve of successive encounters with the response curve we calibrated for single encounters.The error is at 0.05 dex and shows no sign of error accumulation in multiple encounters.In our real application, the maximum Δ/ b from individual disk crossing will be at the level 10 −2 (an estimate is given in Equation 31 and we have explicitly checked in our numerical sampling of minihalos).Therefore, the approach we propose can reasonably accurately approximate the disruption fraction from multiple disk crossings.In the second test, we pick an impact parameter such that the injected energy from a single event is Δ/ b ∼ 10 and repeat the encounter five times.This is a limit where all our assumptions above should fail.In this case, the response curve calibrated underpredicts the amount of mass loss, and the prediction error clearly inflates as the number of encounters increases.

DISRUPTIONS IN THE REALISTIC MILKY WAY ENVIRONMENT
In the sections above, we described the physical process and consequence of stellar encounters and the effects of tidal fields.In the following, we will apply the model to minihalos in the realistic Milky Way environment.

Multiple encounters in the Milky Way disk
Disk stars are the dominant component of the stellar populations in the Milky Way, and thus will be the main contributor to the disruption of dark matter substructures.When a minihalo passes through the stellar disk, it will encounter a slab of stars within a short timescale where  d is the scale height of the Milky Way thin disk and  ⊥ mh is the relative velocity of the minihalo perpendicular to the disk plane.As discussed in Section 4, since this disk passage time is much shorter than the dynamical time of dark matter minihalos, the internal structure of the minihalos will not have enough time to relax between successive stellar encounters during a single passage.Therefore, a series of encounters with disk stars can be effectively considered as one encounter with the injected energy accumulated over the passage.

Back-of-the-envelope model
First, we will consider a simplified model assuming that a minihalo will be completely destroyed after a single encounter with an impact parameter smaller than  min given by Equation 12, and it will remain unperturbed if the impact parameter is larger than  min .Note that we aim for a back-of-the-envelope estimate with this simplified model and this model will not be used for our final analysis.Consider a minihalo with mass  mh and an average density ρmh , moving through a field of stars with the differential number per unit mass   * = d * /(d * d 3 x) characterized by the stellar presentday mass function (PDMF), one can choose a normalization such that the total stellar mass density is ∫  *   * d * =  * .We assume that an encounter with impact parameter  has a probability  = (,  mh , ,  * , v, ...) to destroy the minihalo, where v = v mh − v * , and that the stars have a locally Maxwell-Boltzmann velocity distribution function with velocity dispersion  * which is independent of stellar mass.The destruction rate of the cluster is Taking the destruction probability  to be a step function between  = 0 for   min and  =  0 ∼ 1 for   min as assumed by this simplified model, after the integration, we obtain Integrating this over time for a given minihalo gives a survival probability  survive = exp (−) where  ≡ ∫  d = ∫ dℓ/ mh along the minihalo trajectory.The time-integrated destruction "optical depth"  is dominated by the minihalo time in the disk.If a minihalo stayed in the disk over the entire Hubble time, then where we use the asymptotic approximation  2 ()/() ∼ 1/ at large c.This implies the complete destruction of the minihalo if it is always in the disk.On the other hand, if a minihalo only has a single passage through the disk, then , where Σ * is the stellar surface density of the Milky Way disk, which is reasonably well-fit by an exponential profile In real cases, the survival mass of the minihalo is far from negligible when Δ/ b ∼ 1.The back-of-the-envelope model certainly breaks.This fact is shown in our calibrated response curve in Figure 5  Σ * ≈ Σ 0 exp (−/ d ) with Σ 0 = 816.6M pc −2 and  d = 2.9 kpc (Σ 0 = 209.5 M pc −2 and  d = 3.31 kpc) for the thin (thick) disk of the Milky Way (e.g., McMillan 2011McMillan , 2017)).If we are only interested in the Solar neighborhood ( ≈ 8 kpc, see Bland-Hawthorn & Gerhard 2016 for a review of measurements on  ), the total surface density of both the thin and thick disks is Σ * ∼ 70 M pc −2 .Therefore, the impact of a single passage through the disk is insignificant for reasonably high halo concentrations.
Besides the two extreme cases above, the more realistic scenario is successive passages through the disk, where the number of passages of a minihalo through the disk is a critical factor.Assuming circular orbits, the number of disk passages for minihalos in the Solar neighborhood is approximately  Hubble / circ ∼ O (100), leading to an O (1) effect after accumulation.Therefore, those minihalos in the Solar neighborhood should experience significant disruption.To make a quantitative prediction for the remaining mass and structure of minihalos, we need a more complete model, which will be described in the following.

General response model
In general, the mass loss of minihalos in single (or multiple successive) encounters is characterized by the response curve 1 − Δ mh / mh = F (Δ tot / b , ) was found in Section 3. Assuming the surface density of stars in the disk is large enough so that we can neglect stellar shot noise (discussed further below), the total energy injection Δ tot / b when crossing a localized slab of stars can be calculated by integrating over all possible impact parameters, masses, and velocities of stars where Δ/ b is evaluated using Equation 11,  * is the onedimensional velocity dispersion of stars and  (v * ) is the Maxwell-Boltzmann velocity distribution function .
The integral in the velocity space is evaluated in the  mh  * or  mh  * regime first.The leading order term in the asymptotic limits can be connected using the expression, 1/(  2 * +  2 mh ).The relative error between this expression and the true result is suppressed by the factor (  * / mh ) 2 when  mh  * and vice versa.
is a characteristic mass that depends on the PDMF of stars in the Milky Way, A simple approximation that gives a good fit to the Milky Way data is to take the Kroupa IMF (Kroupa 2001) as the PDMF at  * ≤ 1 M (since the evolutionary effects are small at low masses) with a power-law cutoff   * ∼  −4.5 * at  * > 1 M (e.g., Scalo 1986;Kroupa et al. 1993;Sollima 2019).After the integration, we obtain   0.6 M , which is insensitive to the minimum or maximum star mass assumed.  can be viewed as the characteristic mass of the most effective disruptor, which can vary in different environments.For example, a clumpy medium may have significantly higher   and thus stronger disruption effects.The integration in Equation 27is carried out to infinite distances, where in principle the localized quantities we define in a disk patch no longer apply.This will not affect the results significantly since the contribution from distant stars is suppressed by the 1/ 4 dependence of energy imparted.It is, however, still important to note that Equation 27 will no longer be valid if the stellar surface density is small enough that shot noise becomes important, i.e. when  2 c Σ * /  ∼ 1.If we take   = 0.6 M as estimated above, we obtain the cut-off impact parameter as  c = 0.044 pc If shot noise becomes dominant,  c  s , the total energy injection becomes Connecting the behavior at  c  s , the general solution of the accumulated energy injection during one disk passage can be written as More accurately, we can calculate  c at each value of the stellar mass  * , using the same PDMF and   c ( * ) d * (>  * )/dArea = 1, insert this into Equation 27and then integrate over all masses numerically to define an appropriately-weighted   .Doing so, we find that this gives an "effective"  c which is only ∼ 8% larger than what we obtain using   c Σ * /  = 1.
where we use the same asymptotic approximation of  2 ()/() as in Equation 25.Typically, one has  s >  c for massive minihalos with low concentrations and  s <  c for low-mass minihalos with high concentrations.It is also worth noting that   will cancel out when     , implying that clumpiness of the medium only matters for minihalos with physical sizes comparable to or larger than the typical spacing of the disruptors.
When the Galactic disk is dominated by stars, our treatment by integrating the cumulative perturbations from disk stars to large distances (rather than considering only close encounters like in the back-of-the-envelope model) should be physically equivalent to the disk shocking effect (e.g., Ostriker et al. 1972;Binney & Tremaine 1987;Gnedin et al. 1999;Stref & Lavalle 2017) studied in the depletion of substructures in the Milky Way (e.g., D'Onghia et al. 2010;Stref & Lavalle 2017;Facchinetti et al. 2022).The disk shocking calculation implies (e.g., Binney & Tremaine 1987;Stref & Lavalle 2017) where  z,disk ∝ Σ disk Σ * is the gravitational acceleration at disk vicinity,  mh is the internal velocity dispersion of the minihalo which should scale as √︁   mh / mh .However, taking our Equation 31 to a "smooth" disk limit (  being infinitesimal), we do not get the disk shocking limit naturally.
The key difference is that we assume that all individual stellar encounter events are independent.We evaluate the energy injection in each independent event first before summing up, while the disk shocking evaluates the tidal field of the entire baryonic disk simultaneously before estimating the momentum and energy injection.This independence assumption we made is supported by the fact that the disk scale height,  ∼ O (100) pc, is at least two orders of magnitude larger than the characteristic impact parameter  s ∼  mh of the minihalo, so individual encounters should operate locally in an independent fashion.If we abandon the independence assumption, the break in energy injection at  s in Equation 7(which comes from resolving individual encounters) will disappear.We can integrate the 1/ 4 law until reaching  c and would obtain which is totally consistent with the disk shocking calculation.Depending on the orbit of the minihalo, it can pass through the stellar disk multiple times after falling into the host.When the orbital time is much smaller than the relaxation time of the minihalo, the combined effect can again be modeled as one single passage with the accumulated injected energy.The total number of passages and stellar surface densities at the encounter point should take an ensemble average of all possible stellar orbits passing the location of the detector.In Appendix C, we calculate the relevant correction factors from the ensemble average of all possible orbits with a simplified model.The accumulated injected energy will ultimately be fed to the response function F (ΣΔ tot / b ) to calculate the mass loss.Even if the time between passages is comparable to or larger than the minihalo relaxation time, as discussed in Section 4, it is still reasonably accurate to use the accumulated injected energy and the response curve to approximate the disruption fraction.

Semi-analytic model to combine stellar and tidal disruptions
As the minihalo moves closer to the Galactic center, both  t (Equation 13) and  ts () (Equation 15) will decrease sharply.Therefore, the total mass loss of an infalling minihalo is dominated by its pericenter passages.For minihalos of interest for detection (e.g. in the Solar neighborhood  8kpc), the tidal-stripping time scale during a pericenter passage at  is  (100) Myr, which is of the same order as the minihalo orbital time and much shorter than the lifetime of the minihalo in the host (∼  Hubble ).We assume that the mass of an infalling minihalo outside the tidal radius will be quickly stripped away during the first few pericenter passages before the impact of stellar disruptions start to accumulate.
For simplicity, we evaluate the tidal radius at the target radius  obs of observation, assuming a circular orbit where the Milky Way halo mass distribution is modeled as an NFW profile for dark matter plus a Hernquist profile (Hernquist 1990) for the stellar content where  nfw () ≡ ln (1 + ) − /(1 + ), the host halo parameters are  = 12,  dm = 10 12 M (e.g., Klypin et al. 2002;McMillan 2011;Deason et al. 2012;Bland-Hawthorn & Gerhard 2016).For baryon properties, abundance-matching studies have shown that a Milky-Way mass system typically has a stellar-to-total-mass ratio of  b 0.01  dm (e.g., Moster et al. 2013) and the half mass radius  * 1/2 0.02  vir (e.g., Somerville et al. 2018).The scale radius of the Hernquist profile is related to the stellar-half-mass-radius as  = 0.414 * 1/2 (Hernquist 1990).
The post-stripping density profile of the minihalo is assumed to be a truncated NFW profile at  t , which is equivalent to a normal NFW profile with effective virial radius, concentration, and overdensity as The mass loss due to the tidal disruption is The stripped minihalo forms the initial condition for the following stellar disruptions.Therefore, for a minihalo observed at  obs , the energy accumulated from stellar encounters can be written as (following Equation 31 but considering multiple passages through the disk during the lifetime of the minihalo) where  p is the number of passages through the stellar disk, o denotes averaging over an ensemble of minihalos observed at  obs with all possible orbits.Therefore,  p represents the averaged number of passages over all possible orbits. circ p is the number of passages assuming the minihalo is on a circular orbit with radius  calculated based on the Hubble time  Hubble and the circular orbit period  circ ( obs ). N p characterizes the deviation of  p from this circular orbit estimation.In Equation 38, Σ * is the stellar surface density where the minihalo crossed the disk (the surface density profile Σ * () is given below Equation 26).x denotes averaging over all past disk crossings given the orbit of the minihalo.The correction factor  Σ * characterizes the deviation of the averaged Σ * at all past encounter locations for all possible orbits from Σ * ( obs ).
accounts for the increased effective stellar surface density when the minihalo is not passing perpendicular to the disk, see Appendix C for details).In Appendix C,  N p and  Σ * are estimated based on the orbital model of an isothermal halo.
The combined effect of   ,  N p and  Σ * on Δ tot / b is O (10) at the Solar neighborhood.The velocity term  2 * +  2 mh has a weak dependence on  obs and minihalo orbits, so it is assumed to be the constant value (250 km/s) 2 for simplicity.
We note that  c has an implicit dependence on the surface density at the encounter.When  c  s (when Σ * is large or  mh is small), Δ tot / b will be proportional to Σ 2 * and the correction factor  Σ * should be replaced with  Σ 2 * (see the calculation in Appendix C).To properly account for this, we model the transition from  Σ * to  Σ 2 * empirically as where  = 3 is assumed.We note that the value of  or the detailed form of the transition does not affect the postdisruption mass function in any significant way.

Monte Carlo sampling of the minihalos
We are ready to implement all the physics of disruption discussed above to a sample of minihalos and track their mass loss.We model the evolution of the minihalo population in the Milky Way halo following the steps below: • We initialized a Monte Carlo sample of minihalos.
First, we construct a grid of infall redshifts from , and compute the infall probability at each redshift point as Δ  CDM col ( i ) (using the matter power spectrum from adiabatic CDM fluctuations on small scales, see the discussion in Section 2.1).Then, at each redshift point, the minihalo masses are sampled uniformly over the dynamical range 10 −14 to 10 −3 M .The number densities of these minihalos are calculated following the redshift-dependent pre-infall mass function given in Section 2.1 and Appendix A. The weight of each individually sampled minihalo is proportional to the product of the number density and the infall probability at  i .The minihalo concentrations are calculated following the mass-concentration relation given in Section 2.1.These sampled physical properties represent the initial status of the minihalos upon falling into the Milky Way host.
• Tidal stripping and structural corrections are applied to the sampled minihalos as described in Equation 36and 37. Since the  t solved from Equation 34 after normalizing over  s is independent of minihalo mass, we calculate  t / s on a grid of  and  i for several different choices of  obs and prepare them as lookup tables for efficient interpolation of the tidal radius.
• After the tidal stripping and the implementation of relevant structural corrections, we apply the stellar disruption with the mass loss given by F (Δ tot / b , ).
The cumulative energy injection from stellar encounters, Δ tot / b , is evaluated with Equation 38 using the structural parameters corrected after tidal disruptions.We include the orbit corrections derived in Appendix C. In general, the disruptions taken together induce approximately a 30% decrease in the peak value of the mass function and shift the mass of the peak by roughly half an order of magnitude.The massive end is more strongly affected by disruption than the low-mass end.Bottom: We show the integrated number (left) and mass (right) of minihalos before and after the disruption.The typical survival fraction of minihalos with  mh ≥ 10 −12 M is 83% in terms of number and about 58% in terms of mass.Stellar disruption is the dominant disruption mechanism through the entire mass range of interest.
All the results demonstrated below have passed the convergence tests over hyperparameters in the sampling approach above, including the maximum sample redshift  max i , the mass range of sampling, the number of redshift grid points, the number of samples at each redshift and the resolution of tidal correction grid.

Post-disruption mass functions
In Figure 10, we show the mass function of minihalos (from AMC with the fiducial  a = 25 eV as an example) at the Solar neighborhood ( obs = 8 kpc).The mass functions are presented as the matter mass fraction (with respect to the total dark matter mass) in minihalos per unit logarithm interval (dex) of minihalo mass.We have assumed the axion gives the correct dark matter relic abundance.We present the initial (pre-disruption) mass function, the post-tidal/stellar disruption mass function, the mass function considering the combined effect of tidal and stellar disruptions, and the mass function aggregating tidal and stellar disruptions linearly.With both tidal and stellar disruption, the peak of the mass function is reduced by about 30% accompanied by roughly half an order of magnitude mass shift of the peak.The low-mass end of minihalos is less disrupted than the massive end because lighter minihalos generically form earlier in these models and are more concentrated.Comparing the two disruption mechanisms, the stellar disruption dominates over the entire mass range, which agrees with the conclusion of previous studies of AMC (e.g., Kavanagh et al. 2020;Lee et al. 2021a).It is worth noting that tidal disruption can alter the structural parameters and enhance the average density of minihalos, which could lead to non-linear effects on the following stellar disruption.However, the effect is relatively weak on mass functions in our experiments.We note that a spike shows up in the post-disruption mass function at  mh ∼ 10 −7 M when considering only the stellar disruption.The spike is related to a turnover in the distribution of the minihalo initial concentrations and will be discussed in the following section.In the bottom panels of Figure 10, we show the integrated fraction of disrupted minihalos in number and total mass, respectively.The "survival fraction" quoted hereafter is defined as the fraction of the integrated mass/number of minihalos above a certain mass threshold retained after disruptions.Both the mass and number survival fractions show a plateau as the mass threshold becomes sufficiently low compared to the peak of the initial mass function.Thus we pick  lim mh = 10 −12 M as the limit to measure the "overall" survival fraction, which is insensitive to  lim mh whatsoever.The overall survival fraction of minihalos with  mh ≥ 10 −12 M is about 83% in terms of number and about 58% in terms of mass.The survival fraction will quickly diminish to 30% at  mh 10 −7 M .Again the dominance of the stellar disruption is manifest.
In Figure 11, we show the post-disruption mass function of minihalos at different galactocentric distances,  obs = 4, 8, 16 kpc, for the  a = 25 eV AMC model.We find similar behavior of the post-disruption mass functions at different target radii.In all cases, the massive end is more severely disrupted, and stellar disruption dominates.The disruption at smaller galactocentric distances is stronger primarily due to enhanced stellar surface densities, resulting in the peak of the mass function shifting towards lower masses.The geometric assumptions of our model will break at the central bulge ( obs 1 kpc) of the galaxy.The effective stellar density can be much higher than implied by the disk model, and minihalos will spend most of their lifetime in dense stellar environments.According to the back-of-the-envelope estimation given in Equation 25, the survival probability in this scenario should diminish to zero with moderately high stellar densities.

Impacts of different disruption mechanisms
To better illustrate the impact of tidal disruption, in the top panels of Figure 12, we show the mass distribution of minihalos on the plane of  versus  mh before and after the tidal disruption.The initial mass concentration relation has a big scatter driven by the distribution of  i .At the massive end, the scatter approaches zero since all of the minihalos there have  i ∼  c ∼ 0. Tidal disruptions (along with corrections on minihalo structural parameters, see Equation 36) effectively lower the concentration of minihalos, especially at the massive end.We need to note that the average densities of minihalos are enhanced in the meantime, so the net effects over following stellar disruptions will depend on the trade-off between lowered concentrations and enhanced densities (see Equation 38).
In Figure 13, we show the median mass loss versus the initial mass of minihalo decomposed by the disruption mechanism.The stellar disruption is the dominant mechanism over the entire mass range of interest.However, the tidal disruptions do have non-linear effects in modifying the structural parameters of minihalos prior to stellar encounters.The net effect over stellar disruptions is a trade-off between lowered effective concentrations and enhanced average densities of minihalos.It suppresses (promotes) the stellar disruption at  mh < 10 −6 M ( mh > 10 −6 M ).Nevertheless, the nonlinear effects do not affect much the mass functions or the overall survival probability of minihalos shown in Figure 10 since most of the differences show up at the decreasing edge of the mass function.Another interesting feature to note is the peak of the mass loss curve due to stellar disruptions at around  mh ∼ 10 −5 M (10 −6 if tidal disruption is not introduced).The peak is related to the turnover feature (at the same mass scale) of the mass-concentration distribution

Disruption for different physics models
Here we explore the disruption of minihalos in different physics models summarized in Section 2.1.These models feature different initial mass functions and mass-concentration relations.
In the top panel of Figure 14, we present the pre-and post-disruption mass functions of minihalos in the AMC model with axion mass  a = 1.25, 25, 500 eV.For the EMD case, we show the model with reheating temperature  rh = 15, 30, 60 MeV.The numerical sampling experiments are all conducted at  obs = 8 kpc.Since both the initial mass functions and mass-concentration relations are almost selfsimilar a horizontal mass shift (and the relative mass change only depends on  and  mh ), the post-disruption mass functions are also similar albeit with a horizontal shift.The minihalo abundance reduction from disruptions has similar patterns between the AMC and EMD models.In the bottom panel of Figure 14, we show the mass survival fraction of minihalos.Regardless of the model choice, the overall survival fraction is about 60% for minihalos with  mh > 10 −12 M .It is visible that the AMC (EMD) model with higher  a ( rh ) suffers from stronger disruption at the massive end.This is due to the lower minihalo concentration at the same mass in these models (as shown in Figure 1).The EMD models exhibit sharper decreases in survival fraction at the massive end.The mass function of EMD models features a steeper decline at the massive end, and therefore, at a fixed mass, the same level of a horizontal shift of mass function will give rise to a larger decrease in minihalo abundance.
From the comparisons shown in this section, we can conclude that the model variations have little impact on the mass function and survival fraction of minihalos up to the mass shift noted previously.The estimated mass survival fraction of minihalos in the Solar neighborhood is about 60%.

Galactic survival fraction
In Figure 15, we show the survival (mass) fraction of minihalos after stellar and tidal disruption as a function of the radius of observation.We choose the fiducial AMC and EMD models (with  a = 25 eV and  rh = 30 MeV) for comparison here.The number of surviving minihalos is evaluated by integrating the mass function in three mass bins: log ( mh /M ) ∈ [−12, −10], (−10, −8], (−8, +∞).For both models, the survival fraction of minihalos significantly drops with increasing minihalo mass and decreasing galactocentric distance.Quantitatively, more than 70% of minihalos in the mass range log ( mh /M ) ∈ [−12, −10] survive at any radius of observation, as opposed to 50% ( 30%) survival fraction of minihalos with log ( mh /M ) > −8 at  obs 8 kpc (4 kpc).Low-mass minihalos are less vulnerable to disruption with more concentrated structures, although the survival fraction still decreases sharply at small galactocentric radii due to enhanced stellar surface density.We do not extend this to  < 4 kpc.The bulge component will start to have an impact, lowering the survivability of minihalos even further.
Even though the survival fraction of the most massive minihaloes in the Solar neighborhood can drop significantly because of their diluted structures from the hierarchical assembly, we expect a high overall survival fraction, 60%, of QCD axion miniclusters or EMD minihalos, which is dominated by the concentrated minihalos with low masses.Local measurements like PTAs will be sensitive to minihalos after disruption since they can potentially probe a mass fraction (d  /d log  mh ) well below ∼ 10% (e.g., Dror et al. 2019;Ramani et al. 2020;Lee et al. 2021a).Most of the surviving minihalos are cores with very high densities.We have shown in Figure 6 in Section 3 that the stellar disruption has a limited impact on the central core but primarily destroys the outer shells of the minihalo.The density profile at the outskirts of the minihalo remains steeper than the NFW profile even after the minihalo has fully relaxed from the previous encounter.Direct detection signals can be promoted further for these concentrated remnants (e.g. Lee et al. 2021a).In addition, specifically for the axion scenario, we expect a nonnegligible fraction of free axions in the Solar neighborhood.These free axions could impact direct detection signals in axion haloscope experiments (e.g., Asztalos et al. 2010;Graham et al. 2015;Du et al. 2018).In summary, our findings are encouraging for the prospects of axion direct detection in the post-inflationary scenario as well as minihalos formed from EMD.We leave the calculations of astrophysical signals of dark matter minihalos in the Milky Way for a follow-up study.

CONCLUSIONS
In this paper, we systematically studied the environmental effects on dark matter minihalos (as light as 10 −12 M ) after infall to the Milky Way.Due to the large dynamic range, it is impossible to simultaneously track the evolution of minihalos in a standard cosmological simulation of the Milky Way-mass galaxy.Therefore, we developed a framework to combine small-scale, idealized N-body simulations with an analytic model of the Milky Way galaxy to make these sorts of predictions tractable.
Stellar disruption and tidal disruption are the most critical environmental effects.For stellar disruption, the analytic expressions in the literature are inadequate to accurately capture the minihalo mass loss after stellar encounters.To address this, we developed a semi-analytic model calibrated by a suite of N-body simulations of idealized encounters between a star and a minihalo, varying impact parameters, minihalo masses, and concentrations.On the other hand, for tidal disruption, our N-body simulations show that the disruption effects can be predicted accurately by relatively simple analytic models.To connect to galactic scales, we apply a simplified orbital model of minihalos and derive an orbital-averaged treatment of stellar and tidal disruptions.We perform Monte Carlo simulations to model the mass evolution of minihalos with various masses, concentrations, and infall histories.The framework allows us to make predictions for the survival fraction of minihalos in the Milky Way as a function of galactocentric radius and halo parameters.This paper focuses on the miniclusters of post-inflationary axions and minihalos in the Early Matter Domination.However, the same framework can be easily applied to other types of dark matter substructure models.
Our major findings can be summarized below: • When the energy imparted (Δ) during a stellar encounter is much smaller than the binding energy of the minihalo ( b ), the mass loss is relatively insensitive to minihalo concentration.However, when Δ/ b exceeds unity and continues to increase, minihalos with lower concentration will experience significantly larger mass loss in stellar encounters.This aspect has not been properly considered in many previous studies of dark matter substructures.We propose a simple and intuitive way to model successive encounters by aggregating the imparted energy from individual encounters and applying the same response curve calibrated for single encounters.This method works reasonably well when the imparted energy from individual encounters Δ/ b 0.1 (quasi-static limit).
• For tidal disruption, we confirm the existence of a tidal radius outside which minihalos are largely stripped.The analytical formula works remarkably well even for extremely light minihalos.On the other hand, the stellar disruption results in a shell of marginally-bound particles propagating outwards, leaving a characteristic density profile similar to the Hernquist profile at large radii.The stellar disruption can also slightly reduce the minihalo central densities.
• We show that there could be non-trivial non-linear interactions/combined effects between these disruption processes.The first tidal or stellar disruption will leave dense cores of the original minihalos, which are typically less vulnerable to future disruptions.Therefore, the assumptions on the structures of minihalos need to be revised.Our analytic method includes analytical corrections on minihalo structural parameters after tidal disruptions.
• Applying our methods to well-motivated models like post-inflationary AMC and EMD minihalos, we find a relatively stable mass (number) survival fraction of ∼ 60% (∼ 80%) at the Solar neighborhood.The numbers are insensitive to the limiting minihalo mass we define as survival.The number is also insensitive to model choices.The survival fraction shows a strong dependence on the galactocentric radii, especially for massive minihalos.For example, the mass survival fraction of minihalos above 10 −8 M can drop to 50% ( 30%) at  obs < 8 kpc (4 kpc).Over the entire mass range of interest, stellar disruption is the dominant disruption mechanism.The remaining minihalos are abundant and dense enough to give direct detection signals in e.g.upcoming PTA observations.
The follow-up paper will focus on the direct detection signals of minihalos.In the future, the framework built in this paper can be applied to a spectrum of interesting topics in particle astrophysics, such as the dark matter annihilation rate in minihalos and axion minicluster-neutron star encounters.
We thank Andrea Mitridate for useful discussions on the mass-concentration relation and Gabriel Aguiar for the collaboration in the early stages of this work.HX is supported in part by the United States Department of Energy (DOE) under grant number DE-SC0011637.KZ is supported by a Simons Investigator award and the U.S. DOE, Office of Science, Office of High Energy Physics, under Award No. DE-SC0011632.Support for XS and PFH was provided by the National Science Foundation (NSF) Research Grants 1911233, 20009234, 2108318, NSF CAREER grant 1455342, the National Aeronautics and Space Administration (NASA) grants 80NSSC18K0562, HST-AR-15800.Numerical calculations were run on the Caltech computing cluster "Wheeler", allocations AST21010 and AST20016 supported by the NSF and the Texas Advanced Computing Center (TACC), and NASA HEC SMD-16-7592.The simulation data of this work was generated and stored on the computing system "Wheeler" at California Institute of Technology.The code for the semi-analytic model and the summary of simulation results are available at the project repository.The raw data of the idealized simulations will be shared on reasonable request to the corresponding author.profile are robust.We have also verified that the results are robust to particle number at our fiducial resolution.

B.2. Time-stepping
Since the crossing time of a stellar encounter is orders of magnitude shorter than the internal dynamical time of the minihalo, we need sufficiently small timesteps to resolve the trajectory of the star in the vicinity of the minihalo.The time-stepping parameter in our fiducial simulations is capped at 10 −8 Gyr, which is roughly 1/5  mh / * , where  mh is the minihalo radius and the corresponding minihalo mass is 10 −10 M .In an additional run, the time-stepping parameter is changed to 2 × 10 −8 Gyr while other parameters including the minihalo parameters are exactly the same.We obtain the same disruption fraction, mass profile, and energy change after the halo is fully relaxed.

C. ORBITAL MODEL OF MINIHALOS
The analytic calculation of accumulated energy in Equation 31 comes from only one encounter with the minihalo velocities perpendicular to the plane of the Galactic disk.In reality, minihalos after infall to the Milky Way halo will typically cross the stellar disks multiple times at various locations.To measure the accumulated energy input from a series of disk crossings, we need to evaluate the total number of passages through the disk, the stellar surface density where the encounter occurs, and the angle of incidence, with an ensemble average over all possible orbits for the minihalos eventually found at the observed radius.
For simplicity, we adopt a singular isothermal sphere model following the method in van den Bosch et al. (1999) to estimate the uncertainties related to the orbits of minihalos.The density and potential of the system are given by () =  2 c 4 2 , Φ() =  2 c ln (/ 0 ), (C7) where  c is the constant circular velocity assumed to be 200 km/s and  0 is the zero-potential reference point chosen to be 10 kpc.For a bound test particle in this potential, its halocentric distance  will oscillate between the peri and apocenter with the period where  and  are the energy and angular momentum of the test particle,  1 and  2 are the peri and apocenter of its orbit, which can be determined by solving the equation where  ≤  max =  obs √︁ 2( − Φ( obs ))/ c ().Following van den Bosch et al. (1999), we perform a Monte Carlo sampling of particles in the phase space based on this distribution function.For each sample particle, we compute its orbital period , pericenter and apocenter distances  1 and  2 based on the look-up table created earlier.We note that because of the self-similar nature of the isothermal sphere, the value of  1 and  2 with respect to  obs is independent of  obs , similarly for / circ ( obs ).In Figure 16, we show the probability distribution function of eccentricities, pericenter, and apocenter distances, and orbital periods of sampled minihalo orbits.The distribution should be self-similar for any target radius of observation.In the left panel of Figure 16, we compare the distribution of eccentricity derived here with that from van den Bosch et al. (1999) which match perfectly.
In the right panel of Figure 16, we compute the mean value of  p / circ p =  circ / which gives us the correction factor for the number of passages through the stellar disk (compared to the circular orbit case) as   p ≡  p o / circ p 1.3.o denotes averaging over all possible minihalo orbits.This correction factor is independent of the target radius of observation.
Assuming spherical symmetry, the orbit of a test particle will be confined in a plane and the precession of the orbit will eventually lead to a rosette-like pattern.The phase of the precession when the orbit crosses the disk plane is random.Considering a large ensemble of dark matter particles, the distribution of the radial location where the encounter with the stellar disk takes place will be the same as the probability of the presence of the particle at that distance, and thus equivalent to the time-averaged radial distance of the test particle.Therefore, we have the averaged surface density given the orbit parameter ,  as where x denotes averaging over all past disk crossings and the stellar surface density profile Σ * () is given below Equation 26 following the measurements in McMillan (2011McMillan ( , 2017)).The average of the second-order term Σ 2 * x can be obtained similarly.In the top panel of Figure 17, we show the distribution of Σ * x /Σ * ( obs ) at the target radius 8 kpc and the correction factor  Σ * ( obs ) ≡ Σ * x /Σ * ( obs ) o (averaging over all possible orbits) is about 1.16.We compute the value of  Σ * and  Σ 2 * at several different  obs from 2 to 16 kpc and find that both  Σ * ( obs ) and  Σ 2 * ( obs ) can be fitted by the functional form   +( obs / c )  .The best-fit parameters are  = 0.106,  = 2.03,  c = 12.961,  = 2.048 for  Σ * and  = 0.318,  = 0.781,  c = 5.740,  = 1.628 for  Σ 2 * .An additional correction comes from the enhanced surface density when the minihalo trajectory is not perpendicular to the disk plane.To the leading order, the effective surface density along the trajectory of the incident minihalo (as well as the time duration the minihalo stays in the disk) should scale as 1/cos.Assuming the velocities of dark matter particles are isotropic, the correction factor is where we have imposed a cut-off at cos =  d / d with  d and  d the scale height and length of the disk, assuming to be 400 pc and 3 kpc, respectively.Particles with even smaller incidence angles stay in the disk and will become completely disrupted (see Equation 26) which will have no impact on the averaged energy imparted.Combining the two effects above, we obtain the correction factor for the effective stellar surface density.

Figure 1 .
Figure1.Top: Initial (pre-disruption) mass function of minihalos in different physics models (  is the mass fraction with respect to the total dark matter mass in the Universe).Here we show the mass functions for the AMC models with axion mass  a = 1.25, 25, 500 eV and the EMD models with reheating temperature  rh = 15, 30, 60, 120 MeV.The model parameter variations manifest as constant horizontal shifts of the mass function (see Appendix A for details).The EMD models exhibit sharper peaks than the AMC models.Bottom: Mass-concentration relation of minihalos in different models shown in the top panel.The mass-concentration relation in the EMD model is peaked due to its peaked matter power spectrum.In general, the halo concentration hits a floor at the massive end as the adiabatic CDM power spectrum takes over.At the low-mass end, since one would not expect minihalos to form before matter-radiation equality, a cap to the concentration appears at around (1 +  i ) = 10 4 .

Figure 2 .
Figure2.The evolution of a minihalo during stellar disruption is visualized.This is a simulation of an encounter between a minihalo of mass  mh = 10 −10 M and a star of mass  * = 1 M , with the impact parameter  = 10 −4 kpc.The star passed by in the y-direction of the image.We show the dark matter surface density distribution of the minihalo.The left panel shows the initial dark matter distribution.Right after the encounter, a large fraction of dark matter has been heated up by the gravitational interaction with the passing star and become unbounded.Although this energy is transferred impulsively, it requires roughly the minihalo dynamical time ( dyn 1 Gyr) for the disruption to be reflected in the minihalo density distribution, as shown in the middle and right panel.

Figure 3 .
Figure 3. Fractional mass loss as a function of the normalized imparted energy from stellar encounters with minihalos.Results are shown for minihalo concentrations  = 10, 30, 100, 500.The simulation results and the best-fit model are shown in solid points and lines.For comparison, the purple curve is the analytic prediction from Kavanagh et al. (2020) for  = 100 minihalos.The dot-dashed lines show the semi-analytical model developed in Delos (2019).The zoom-in subplot in the bottom left corner shows the transitional regime around Δ/ b ∼ 1.The asymptotic behavior of the response curve at large Δ/ b has a significant concentration dependence.The response curves calibrated from our simulations are in reasonably good agreement with results in Kavanagh et al. (2020) and Delos (2019) in the shared dynamical ranges.

Figure 4 .
Figure 4. Fractional mass loss as a function of concentration with impact parameter  = 0.05 pc.The blue solid curve is the semianalytic model prediction while the data points are obtained from simulations.

Figure 5 .
Figure 5. Fractional mass loss as a function of minihalo mass with concentration  = 100 and impact parameter  = 0.05 pc.The black solid points are simulation results.The blue curve is the semianalytic prediction using the best-fit response curve.The vertical line indicates where the transition impact parameter  s is reached.The semi-analytical model accurately predicts the location and the shape of this transition.

Figure 6 .
Figure6.The evolution of minihalo density profiles after a close encounter with a star.For the tests here, the minihalo mass is set as 10 −10 M initialized with the NFW profile at redshift zero (with virial radius  mh ∼ 0.01 pc and dynamical time scale  dyn ∼ 2.2 Gyr).The halo concentration is  = 100.The impact parameter is  = 2 × 10 −5 kpc and  = 5 × 10 −5 kpc in the top and bottom panels, respectively.The vertical shaded regions indicate the range of convergence radii at different times.After the close encounter with the star, the outskirts of the minihalos are dominantly disrupted in less than a dynamical time.Eventually 81% (60%) of the minihalo mass is disrupted for the  = 2 × 10 −5 kpc ( = 5 × 10 −5 kpc) case.The density profile of minihalo after disruption can be well-fitted by a broken power-law profile with best-fit asymptotic slopes around −4 analogous to that of the Hernquist profile.The "model fit" curve provides an asymptotic limit of the density profile after infinite time.The upturn at the outskirt of the minihalo represents a propagating shell of unbound particles.As this upturn propagates towards the outer part of the minihalo, the density profile gradually converges to the asymptotic limit.

Figure 7 .
Figure7.Top: The enclosed mass profile of a minihalo with  = 10 −10 M at  = 0 − 0.5 Gyr.The minihalo is on a circular orbit at  =  = 8 kpc in a Milky Way-mass dark matter halo.For reference, the dynamical time of the host halo within this radius is about 0.1 Gyr.The gray and red vertical lines indicate the scale radius and the tidal radius of the minihalo.The mass outside the tidal radius is stripped away at roughly the dynamical time scale of the host halo while the mass inside is marginally perturbed.Bottom: The density profile of the minihalo at the same time as the top panel.The plot illustrates the behavior of the matter distribution at the outskirts of the minihalo during tidal stripping.

4. 1 .Figure 8 .
Figure8.Mass evolution history of the minihalo in tidal disruption (the same one as in Figure7).The total mass here is the sum of dark matter mass within the initial virial radius of the minihalo.In addition, we show the mass inside and outside the tidal radius.After approximately a dynamical time of the host halo, the mass outside the tidal radius starts an exponential decay as described by Equation14with the fudge factor  1.8.The mass inside the tidal radius is marginally affected.

Figure 9 .
Figure 9. Mass loss of the minihalo from single/multiple encounters.The blue line and open circles show the response curve of  = 100 minihalos we calibrated in Section 3. The green circles show the minihalo mass loss as a function of accumulated imparted energy from successive encounters (with Δ/ b ∼ 10 −1 ).The red circles show the same test but with typical Δ/ b ∼ 10.Our treatment of accumulating imparted energy from independent encounters works reasonably well when individual Δ/ b 1, but will underpredict the mass loss when Δ/ b 1.
and has been highlighted in studies on CDM substructures (e.g.van den Bosch et al. 2018).

Figure 10 .
Figure 10.Top: Mass function of minihalos (from axion miniclusters) at the Solar neighborhood ( obs 8 kpc).We assume the AMC model with  a = 25 eV.The mass function before disruption is shown as the gray dashed line.The mass function after processing only tidal (or stellar) disruption is shown as the blue (red) solid line.The mass function post-disruption, combining both tidal and stellar disruption, is shown as the solid black line.The mass function post-disruption, combining tidal and stellar disruption linearly, is shown as the dashed black line.In general, the disruptions taken together induce approximately a 30% decrease in the peak value of the mass function and shift the mass of the peak by roughly half an order of magnitude.The massive end is more strongly affected by disruption than the low-mass end.Bottom: We show the integrated number (left) and mass (right) of minihalos before and after the disruption.The typical survival fraction of minihalos with  mh ≥ 10 −12 M is 83% in terms of number and about 58% in terms of mass.Stellar disruption is the dominant disruption mechanism through the entire mass range of interest.
Figure 11.Mass function of minihalos (from AMC with  a = 25 eV) observed at three galactocentric distances ( obs = 4, 8, 16 kpc).A similar reduction pattern in the mass density of minihalos is found at each distance, with stronger disruption at smaller radii and the massive end of the mass function.The peak of the postdisruption mass function is slightly shifted towards lower masses at smaller radii.Stellar disruption is the dominant disruption mechanism at all radii.

Figure 12 .
Figure 12.Top: Mass distribution of sampled minihalos in the phase space ( versus  mh ) before disruption.The mass concentration relation from Lee et al. (2021a) is shown with the white dashed line (assuming  i = 0).The down scatter of concentration on this plane is driven by the spread of the minihalo infall redshift,  i .The scatter shrinks to zero towards the massive end, since most of the massive minihalos have  i ∼ 0. Bottom: Mass distribution of sampled minihalos in the phase space after tidal disruption.Effectively, tidal disruption lowers the concentration of minihalos to  t / s while increasing the averaged densities of minihalos (since the low-density outskirt has been shredded).The effect is significant in massive, low-concentration minihalos.

Figure 13 .
Figure13.Median mass loss induced by disruption versus the initial mass of minihalos.We have decomposed the mass loss to tidal and stellar disruptions.The red dashed line shows the mass loss due to stellar disruption when there are no prior tidal corrections.The black dashed line shows the total mass loss if adding stellar and tidal terms linearly.In all cases, stellar disruption is the dominant mechanism.The non-linear effects of tidal disruption suppress stellar disruption at  mh < 10 −6 M , but promote that at  mh > 10 −6 M .Nevertheless, this does not affect the mass functions shown in Figure10, since the differences show up at the mass where the density of minihalos is dropping rapidly.

Figure 14 .
Figure14.Top: Mass function of minihalos in different models before (dashed) and after (solid) disruption.Two different models are studied: post-inflationary AMCs and EMD.The reheating temperature for the EMD era is  rh = 15, 30, 60 GeV while the axion mass in the post-inflationary AMC scenario is chosen as   = 1.25, 25, 500 eV.These parameters are purely chosen for illustrative purposes and one can use other parameters which will shift the mass range, but the shape of the mass function will remain the same.Bottom: The integrated mass survival fraction of minihalos for the models shown in the top panel.Regardless of the model choices, the overall survival fraction of minihalos with  mh > 10 −12 M is stably around 60%.

Figure 15 .
Figure 15.Mass survival fraction of minihalos as a function of galactocentric distance.In addition to the overall survival fraction after stellar and tidal disruptions, we also show the survival (mass) fraction of minihalos in three mass bins: log ( mh /M ) ∈ [−12, −10], (−10, −8], (−8, ∞).The survival fraction has a strong dependence on minihalo mass.Massive minihalos are more severely disrupted, especially at small galactocentric radii.Less than a half of the minihalo in the (−8, ∞) bin can survive at  obs 8 kpc.On the other hand, low-mass minihalos in the [−12, −10] bin have 70% survival fraction at any radius.
to the scale-free nature of the singular isothermal profile, the equations can be simplified by defining the maximum angular momentum  c () =  c () c , where  c () is the radius of the circular orbit with energy  given by  c () =  0 exp , the circularity parameter is defined as  = / c () and Equation C9 is reduced to 1  ≡ / c () and the two solutions correspond to the peri and apocenter distances, which are independent of the energy of the test particle if normalized by  c ().Similarly, if normalized by  c ()/ c ,  also becomes independent of .The look-up tables of ,  1 , and  2 with respect to  are computed numerically.Here we do not consider the scattering and mergers of minihalos within the parent halo and treat their orbits as unperturbed from various relaxation mechanisms.Assuming spherical symmetry, the ergodic phase-space distribution function  () can be derived through the Ed- ≡ −Φ and  ≡ −.For the singular isothermal profile, the solution is simply the Boltzmann distribution() =  exp (−2/ 2 c ), (C13)where  is a constant normalization factor.Given a target radius for observation  obs , the normalized phase-space probability density function of dark matter particles localized around  obs is (, )|  obs = 4  2 obs ( obs ) ()  √︁ 2( − Φ( obs )) −  2 / obs ,(C14) where  ≥ Φ( obs ) and  ≤ 2 obs √︁  − Φ( obs ) are required.Replacing  as  c (), we obtain  (, )|  obs = 4  obs ( obs ) ()  c ()