<figure><img src="image_1.png" /></figure>

HAL
open science

# Tidal stripping of dark matter subhalos by baryons from analytical perspectives: disk shocking and encounters with stars

Gaétan Facchinetti, Martin Stref, Julien Lavalle

► To cite this version:

Gaétan Facchinetti, Martin Stref, Julien Lavalle. Tidal stripping of dark matter subhalos by baryons from analytical perspectives: disk shocking and encounters with stars. 2022. hal-03539651

HAL Id: hal-03539651

https://hal.science/hal-03539651

Preprint submitted on 21 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Tidal stripping of dark matter subhalos by baryons from analytical perspectives: disk shocking and encounters with stars

Gaétan Facchinetti,¹,* Martin Stref,²,† and Julien Lavalle¹,‡

¹Laboratoire Univers & Particules de Montpellier (LUPM), CNRS & Université de Montpellier (UMR-5299), Place Eugène Bataillon, F-34095 Montpellier Cedex 05 — France

²LAPTh, Université Savoie Mont Blanc & CNRS, Chemin de Bellevue, 74941 Annecy Cedex — France

The cold dark matter (CDM) scenario predicts that galactic halos should host a huge amount of subhalos possibly as light as or lighter than planets, depending on the nature of dark matter. Predicting their abundance and distribution on such small scales has important implications for dark matter searches and searches for subhalos themselves, which could provide a decisive test of the CDM paradigm. A major difficulty in subhalo population model building is to account for the gravitational stripping induced by baryons, which strongly impact on the overall dynamics within the scale radii of galaxies. In this paper, we focus on these “baryonic” tides from analytical perspectives, summarizing previous work on galactic disk shocking, and thoroughly revisiting the impact of individual encounters with stars. For the latter, we go beyond the reference calculation of Gerhard and Fall (1983) to deal with penetrative encounters, and provide new analytical results. Based upon a full statistical analysis of subhalo energy change during multiple stellar encounters possibly occurring during disk crossing, we show how subhalos lighter than ~1 $M_\odot$ are very efficiently pruned by stellar encounters, and how that modifies their mass function in a stellar environment. If reasonably resilient, surviving subhalos have lost all their mass but the inner cusp, with a tidal mass function strongly departing from the cosmological one; otherwise, their number density can drop by an order of magnitude at the solar position in the Milky Way with respect to disk-shocking effects only. For illustration, we integrate these results into our analytical subhalo population model. They can easily be incorporated to any other analytical or numerical approach. This study complements those using cosmological simulations, which cannot resolve dark matter subhalos on such small scales.

PACS numbers: 12.60.-i,95.35.+d,98.35.Gi

## I. INTRODUCTION

The cold dark matter (CDM) scenario is tied to hierarchical structure formation [1–6], in which small-scale halos much smaller than typical galaxies collapse first in a denser universe and then larger and larger halos assemble through mergers and accretion. Consequently, the distribution of dark matter (DM) in galactic halos such as that of the Milky Way (MW) is expected to exhibit inhomogeneities in the form of smaller structures spanning a wide range of masses. Although subject to tidal stripping, a significant fraction of these subhalos are to survive and populate their host galaxies in number, as explicitly verified in cosmological simulations down to their numerical resolution limits [7–10].

Paradoxically enough, the structuring of CDM on small scales may also lead to a mismatch between theoretical predictions and observations, all this being termed as the “CDM small-scale crisis” (see e.g. Ref. [11] and references therein). If the related core-cusp [12] and diversity [13] problems are certainly serious issues, especially when contrasted with the impressive regularity observed

in some scaling relations with baryons [14], other aspects related to subhalos (counting, etc.) may, as for them, find reasonable explanations as being related to baryonic effects or feedback [15]. It is obviously necessary to inspect DM-only solutions to these potential problems (see e.g. Ref. [16] or [17]), but it is not less important to improve our understanding and description of CDM physics itself on small scales to prepare for additional tests. In this respect, having reliable predictions of the properties of subhalo populations in galaxies, especially of those subhalos light enough not to host baryons, is of particular interest. Indeed, subhalos can imprint gravitational signatures [18, 19], or boost potential DM annihilation signals [20–35], turn into intermittent local DM winds [36, 37], hence providing additional ways to test the CDM scenario.

High-resolution cosmological simulations provide very important clues to understand the formation and evolution of subhalos, but are limited by two aspects: (i) finite spatial or mass resolution, and (ii) the fact they can hardly be matched to specific real galaxies with constrained kinematic properties and specific histories. The former aspect may have strong impact on the way subhalo properties are inferred from simulations (see e.g. Refs. [38–40]), while the latter makes it potentially dangerous to blindly extrapolate simulations’ results (for instance subhalo spatial distributions, mass functions, etc.) to specific objects like the MW [30]. Even with im-

<footnote>* gaetan.facchinetti@umontpellier.fr</footnote>

<footnote>† martin.stref@lapth.cnrs.fr</footnote>

<footnote>‡ lavalle@in2p3.fr</footnote>
proved resolution and empirical implementation of baryonic physics, cosmological simulations will hardly be able to probe most of the substructure mass, which could theoretically reside in subhalos with virial masses as light as ~ $10^{-12}$ M$_{\odot}$ [41–44].

Deepening our physical understanding of the outcomes of simulations is therefore desirable to consistently interpolate their properties down to smaller scales or onto real objects. In the meantime, it is important to develop alternative though complementary analytical or semi-analytical approaches, since these can deal with scales unresolved by simulations, and are also well suited to study other effects like changes in cosmological parameters, in the primordial power spectrum, etc. These alternative approaches are particularly interesting to investigate the effects of subhalos in DM searches and to conceive related tests of the CDM scenario itself [18, 19, 29, 30, 45–52].

In this paper, we will resort to analytical methods to study those gravitational tides experienced by subhalos and generated by the baryonic components of galaxies, which are expected to strongly affect the subhalo properties within the scale radii of galaxies. This notably concerns regions where DM and/or subhalo searches are currently conducted. We will address two different physical phenomena with two different timescales. First, we will briefly review the pruning of subhalos induced by tidal shocks triggered by crossings of galactic disks in spiral galaxies, called disk-shocking effects. During such crossings, which may last for rather long times with respect to the deep inner orbital timescale in subhalos, the stars and gas confined into disks act collectively as a smooth gravitational field. The analytical procedure that we present to account for disk shocking was actually developed in a previous work [30], which we slightly refine here. The second phenomenon is technically more involved, and concerns the tidal stripping induced by individual encounters of subhalos with stars as they pass by each other. The effects of such encounters on subhalos, which occur on much shorter timescales, have already been considered in both analytical [51, 53, 54] and numerical studies [55–58], in which they were shown to be significant. Here, we aim at revisiting this physical problem, notably by improving over an earlier reference calculation meant to describe a singular encounter and presented in Ref. [59] (GF83 henceforth), and widely used in subsequent literature. We further aim at gauging the impact of these baryonic tidal effects on the whole subhalo population of a template galactic host. To do so, we will integrate these new results in the analytical subhalo population model that we designed in a previous work [30] (SL17 hereafter), tuned to describe the subhalo population of the MW consistently with kinematically constrained Galactic mass models comprising both DM and baryons [60].

The paper is organized as follows. In Sect. II we shortly introduce the very bases of the SL17 model, which are more detailed in App. A, and describe the way disk-shocking effects can be analytically described and ac-

<figure><img src="image_2.png" /><figcaption>FIG. 1. Number density of subhalo with respect to the position from the GC for a minimal subhalo cosmological mass of $m_{200}^{\min}=10^{-10}$ M$_{\odot}$ and a mass index $\alpha = 1.9$. The effect of the disk are switched on in the case (sm. + disk) and switched off in the (sm. only) configuration. In comparison is shown the cosmological distribution) in green (superposed to the red dash-dotted curve).</figcaption></figure>

counted for in subhalo models. In Sect. III, we turn to individual stellar encounters, and present the computation of the total energy kick felt by test particles bound to a subhalo that passes by a single star. Then, in Sect. IV, we address the computation of the impact parameter distribution and the probability of encounter before evaluating the total energy kick induced by one crossing of the stellar disk. The consequences in terms of the SL17 population model are discussed and illustrated in Sect. V, before concluding in Sect. VI. Further technical details are given in the appendix sections.

## II. DISK-SHOCKING EFFECTS ON A SUBHALO POPULATION

In this section, we shortly review the effect of disk shocking in the SL17 model. We first summarize the SL17 subhalo population model (more details can be found in App. A), and briefly discuss the gravitational tides generated by the global gravitational potential of the host galaxy (here the MW) before addressing disk shocking. We then propose an easy way to quickly implement these effects, induced by smooth gravitational potentials, into subhalo population models. We illustrate our results with the SL17 model.

The semi-analytical SL17 model [30] (see also Refs. [33,
35, 61]) was designed to consistently incorporate a smooth DM halo, baryonic components (disk/s, bulge/s), together with a subhalo population covering the entire mass range allowed by particle DM models, into a global galactic mass model for a target galaxy. It was calibrated for the MW from the mass model constrained on kinematic data by McMillan in Ref. [60]. The model assumptions are as follows: (i) subhalos are building blocks of galactic halos; (ii) if they were hard spheres, they would simply spatially track the overall DM density profile, and retain their initial properties (cosmological mass function, concentration distribution, inner density profiles, etc. — which can be estimated by investigating the properties of field subhalos); (iii) tidal stripping and mergers are responsible for altering and depleting them. The SL17 model therefore consists in evolving a subhalo population following an overall spherically symmetric DM density profile, starting from initial cosmological properties (virial masses $m = m_{200}$, concentration $c = c_{200}$, position $R$), and then plug in tidal effects to redistribute the DM stripped away from subhalos to the smooth DM component. The model is statistical in essence, because subhalos are described by probability density functions (PDFs) for their position $R$, mass $m$, and concentration $c$. Tidal stripping is then responsible for spatially-dependent mass losses, which make subhalos with initially the same virial mass $m$ end with different tidal (hence physical) masses $m_t(m, c, R)$ further depending on both position and concentration. This procedure makes the overall parametric PDF fully intricate, as tidal masses exhibit strong dependencies on the other parameters. Although the SL17 model can be derived for any assumption as for the inner density profiles of subhalos, we will assume Navarro-Frenk-White (NFW) [62] inner profiles in this paper.

Additionally, the SL17 model allows for tidal disruption of subhalos based on a rather simple criterion inspired from numerical studies [63], in which it was shown that subhalos tidally pruned down to within their scale radii, $r_s$, would actually be destroyed. By defining $x_t$ as the ratio of the subhalo tidal radius $r_t(m, c, R)$ to the scale radius $r_s$, then it is possible to fix a threshold $\epsilon_t$ below which a subhalo is destroyed (i.e. if $x_t < \epsilon_t$). Although early dedicated studies seemed to indicate $\epsilon_t \sim 1$ [63], this result has been strongly questioned in more recent studies [38–40, 64, 65], where it was shown that estimates of the tidal disruption efficiency could be biased by numerical effects, which means that $\epsilon_t$ could very likely take much smaller values. For the sake of completeness, we adopt two reference choices in this paper:

$$ \begin{cases} \epsilon_t = 1 & (\text{for fragile subhalos}) \\ \epsilon_t = 0.01 & (\text{for resilient subhalos}) \end{cases} . \quad (1) $$

As a result of tidal stripping and disruption, the most concentrated subhalos are found to be the most resistant, as expected. Thus, based on tidal stripping effects, the SL17 model is able to quite naturally explain the flattening of the subhalo number density and its further

depletion as one gets closer and closer to the galactic centers, which is observed in cosmological simulations [9, 66, 67], as well as the spatial dependence of the mass-concentration relation [26, 67, 68]. All the SL17 model concentrates in providing the differential number density of subhalos $n_{sub}$, which can be expressed either in terms of virial mass $m$, or in terms of tidal mass $m_t$, such that

$$ \frac{d^2 n_{\text{sub}}(m, c, R)}{dc \, dm_t} = \int dm \frac{d^2 n_{\text{sub}}(R, c, m)}{dc \, dm} \delta(m - m_t(m, c, R)) $$

where all relevant parameters appear explicitly.

Two distinct tidal effects are accounted for in the original SL17 model. The first and simplest effect is the tidal stripping of subhalos by the smooth gravitational potential of the whole Galaxy. This can be modeled by assigning to each subhalo a tidal radius solution to the following equation

$$ r_{t, \text{sm}} = R \left( \frac{m(r_t)}{3 M(R) \left(1 - \frac{1}{3} \frac{d \ln M}{d \ln R}\right)} \right)^{1/3} \quad (3) $$

where $M(R)$ is the total mass of the Galaxy within radius $R$ (including both the DM and baryons) [69]. This relation, based on the extrapolation of the Roche criterion to diffuse objects, has been shown to nicely correlate with simulation results [9, 70].

The second effect, due to baryons this time, and that we want to review here, is the gravitational shock induced at each crossing of the disk, called disk shocking. This gravitational shocking is generated by the smooth potential of the Galactic disk, inside which gas and stars act collectively. Indeed, test-mass particles (e.g. DM particles) bound inside a subhalo that crosses the disk experience a velocity kick due to the “rapidly” changing gravitational potential. In order to evaluate this kick, we use the impulsive approximation: the crossing is considered fast enough so that particles are frozen in the frame of the subhalo. When averaged over a radial shell in the subhalo, it is given by [71]

$$ \Delta \mathbf{v}_d = \frac{2 g_d}{\sqrt{3} v_z} r \hat{\mathbf{e}}_z . \quad (4) $$

This expression depends on $\hat{\mathbf{e}}_z$, the unit vector normal to the galactic plane, $v_z$ the associated subhalo velocity component, and $g_d$ the gravitational acceleration due to the potential of the disk. This velocity kick can be translated into a kick in kinetic energy per unit of particle mass as

$$ \Delta E_d = \frac{1}{2} \left[ (\mathbf{v} + \Delta \mathbf{v}_d)^2 - \mathbf{v}^2 \right] = \frac{1}{2} |\Delta \mathbf{v}_d|^2 + \mathbf{v} \cdot \Delta \mathbf{v}_d . \quad (5) $$

with $\mathbf{v}$ the initial velocity in the frame of the subhalo. Averaging over an initial velocity distribution, the second term vanishes so that one considers only $\Delta E_d \sim (\Delta \mathbf{v}_d)^2/2$. However, the impulse approximation often breaks down in the case of disk shocking, and adiabatic invariance has to be accounted for — if the inner orbital period is short enough with respect to crossing time, then
particles are further protected against stripping by virtue of angular momentum conservation. In the end, according to Ref. [72], the energy gain is

$$ \Delta E_d = \frac{1}{2} (\Delta \mathbf{v}_d)^2 A_1(\eta_d) = \frac{4 g_d^2}{3 v_z^2} r^2 A_1(\eta_d), \quad (6) $$

where $A_1(\eta) = (1 + \eta^2)^{-3/2} \le 1$ is a corrective factor. The adiabatic parameter for disk shocking is given by $\eta_d = t_d \omega \ge 0$ where $t_d = H_d/v_z$ is the crossing time with $H_d = 0.9$ kpc the thickness of the disk. The orbital frequency of DM particles $\omega$ is approximated by $\omega = \sigma_{\text{sub}}/r$, with $\sigma_{\text{sub}}$ the velocity dispersion in the subhalo evaluated using Jeans' equation<fnref n="1" /> [73]. Whenever $\eta_d \gg 1$ i.e. adiabatic shielding is efficient, the energy kick is suppressed by the corrective factor $A_1$.

The tidal radius is then evaluated recursively. The number of disk crossings $N_{\text{cross}}$ is computed with the assumption that subhalo orbits are circular. The algorithm starts with $r_{t,0} = r_{t,\text{sm}}$ given by Eq. (3), and for every crossing it evaluates a new value of $r_t$ by the requirement that if the energy kick in a shell is greater than the gravitational potential of the structure at that position, then the entire shell is removed. More precisely, we make explicit the dependencies on the radius and on the tidal extension $r_t$ of the energy gain function by writing $\Delta E_d = \Delta E_d(r, r_t)$. We denote by $\Phi(r, r_t)$ the gravitational potential

$$ \Phi(r, r_t) = - \int_r^{r_t} dr' \frac{Gm(r')}{r'^2} \quad (7) $$

and the successive tidal radii are evaluated by solving for $r_{t,i+1}$ in the equation

$$ \Delta E_d (r_{t,i+1}, r_{t,i}) = |\Phi (r_{t,i+1}, r_{t,i})| \quad (8) $$

for all $i \in [1, N_{\text{cross}}]$. The tidal radius today is defined by $r_t \equiv r_{t,N_{\text{cross}}}$. From Eq. (A5) we plotted, in Fig. 1, the resulting subhalo number density (integrated over concentration)

$$ n_{\text{sub}}(R) = \int dm_t \frac{dn_{\text{sub}}}{dm_t}, \quad (9) $$

with the distance $R$ to the GC. Two configurations are considered, whether disk shocking effects are switched

<figure><img src="image_3.png" /><figcaption>FIG. 2. Geometry of the problem. The clump in represented as a sphere with center $C$ and radius $r_t$.</figcaption></figure>

on or not. For comparison, the cosmological/unevolved distribution is also represented. When subhalos are fragile (disruption parameter $\epsilon_t = 1$), we find a strong suppression due to disk shocking toward the center of the Galaxy compared with the cosmological distribution and with smooth stripping induced by the overall halo only. The impact is less important with smaller and smaller values of $\epsilon_t$.

The next sections are dedicated to evaluate the effects of stellar encounters on the tidal radius and subsequently on the subhalo population. The relative impact of disk shocking and stellar encounters will be compared in Sec. V.

### III. SINGLE STAR ENCOUNTER

We focus on the case of an encounter between one subhalo and a single star. The main objective is to compute the total energy received by every particle in the subhalo during the crossing. We first describe a complete parameterization of the problem before moving on, in a second part, to the computation and end up with our results and see how they compare to previous results of the literature.

Our calculation starts by closely following the original work by Spitzer [74] and its extension by Gerhard and Fall [59] (hereafter GF83), but then it extends it even further to penetrative encounters. The geometry of the encounter between the subhalo and the star is summarized in Fig. 2. The star $S$ is a point-like object with mass $m_*$ while the subhalo has a radial extension (tidal radius) $r_t$, a mass $m_t \equiv m(r_t)$ and his center of mass is called $C$. We assume the DM is spherically distributed around $C$ and this spherical symmetry is maintained for the duration of the encounter. The center of mass of the entire system defines an inertial frame and we introduce a fixed point $O$ in that frame. By Newton's second law, we have for a DM particle $M_i$

$$ \frac{d^2}{dt^2} \mathbf{OM}_i = - \frac{G m_*}{S M_i^3} \mathbf{SM}_i + \sum_{j \neq i} \frac{G m_p}{M_i M_j^3} \mathbf{M}_i \mathbf{M}_j, \quad (10) $$

where $m_p$ is the mass of a DM particle. The second term

<footnote n="1">¹While this expression was used in SL17, $v_z$ and $\omega$ are computed slightly differently in this new analysis. Indeed $v_z$ was evaluated according to the circular velocity of the subhalo assuming an isothermal profile for the Galactic DM distribution. Here, in order to be consistent with the rest of the study, we approximate this quantity as the average of a Maxwell-Boltzmann, isotropic, distribution with an NFW profile for the total Galactic density. This yields $v_z \sim \sqrt{2/\pi\sigma_c(R)}$, where $\sigma_c(R)$ is the velocity dispersion computed from Jeans' equation. Moreover, $\omega$ was evaluated for an isothermal sphere, while here the true profile of the subhalos is used.</footnote>
on the right-hand side of this equation accounts for the self-gravity of the subhalo. In the following, we are going to assume that the typical period of DM particles inside the clump is much larger than the duration of the encounter so the internal dynamics of the subhalo is effectively frozen and the self-gravity can be neglected. This is called the *impulse approximation* [74] and we discuss its validity further in Sec. IV C. The dynamics of the clump is governed by

$$ \frac{d^2}{dt^2} \mathbf{OC} = \frac{m_p}{m_t} \left\{ -\sum_i \frac{G m_\star}{S M_i^3} \mathbf{SM}_i \right\}. \quad (11) $$

We introduce the positions with respect to the center of the clump $\mathbf{r} \equiv \mathbf{CM}$ and $\mathbf{r}_\star \equiv \mathbf{CS}$, and the velocity of a DM particle with respect to $C$

$$ \mathbf{v}_i \equiv \frac{d}{dt} \mathbf{r}_i, \quad (12) $$

which obeys

$$ \frac{d\mathbf{v}_i}{dt} = -\frac{G m_\star}{|\mathbf{r}_i - \mathbf{r}_\star|^3} (\mathbf{r}_i - \mathbf{r}_\star) + \frac{M_\star}{m_t} \sum_j \frac{G m_p}{|\mathbf{r}_j - \mathbf{r}_\star|^3} (\mathbf{r}_j - \mathbf{r}_\star) \mathbf{\delta}\mathbf{\beta} $$

Taking the continuous limit, this becomes

$$ \frac{d\mathbf{v}}{dt} = -\frac{G m_\star}{|\mathbf{r} - \mathbf{r}_\star|^3} (\mathbf{r} - \mathbf{r}_\star) + \frac{M_\star}{m_t} \int \frac{\rho(r')}{|\mathbf{r}' - \mathbf{r}_\star|^3} (\mathbf{r}' - \mathbf{r}_\star) d^3\mathbf{r}'. \quad (14) $$

The second term on the right-hand side can be further simplified by using the spherical symmetry to get

$$ \frac{d\mathbf{v}}{dt} = -\frac{G m_\star}{|\mathbf{r} - \mathbf{r}_\star|^3} (\mathbf{r} - \mathbf{r}_\star) - \frac{G M_\star m(r_\star)}{m_t} \frac{\mathbf{r}_\star}{r_\star^3}, \quad (15) $$

where $m(r_\star)$ is the subhalo mass within $r_\star$. We are interested in the net change in velocity

$$ \delta \mathbf{v} = \int_{-\infty}^{+\infty} \frac{d\mathbf{v}}{dt} dt. \quad (16) $$

To go further, we assume the encounter happens at high enough speed so the trajectory of the subhalo can be approximated by a straight line. In that case, we have $\mathbf{r}_\star(t) = -\mathbf{b} - \mathbf{v}_r t$ where $\mathbf{b}$ is the impact vector directed from $S$ to $C$ at the time of closest approach and $\mathbf{v}_r$ is the constant relative velocity. Integration over time leads to

$$ \delta \mathbf{v} = -\frac{G M_\star}{v_r} \frac{\mathbf{r} + \mathbf{b} - \hat{\mathbf{e}}_v (\mathbf{r} \cdot \hat{\mathbf{e}}_v)}{r^2 - (\mathbf{r} \cdot \hat{\mathbf{e}}_v)^2 + b^2 + 2 \mathbf{r} \cdot \mathbf{b}} + \frac{2 G M_\star}{v_r b^2} I(b, r_t) \mathbf{b}. \quad (17) $$

In this expression, we have introduced a unitary vector $\hat{\mathbf{e}}_v \equiv \mathbf{v}_r / v_r$ and the following integral

$$ I(b, r_t) \equiv \int_0^\infty (1 + x^2)^{-3/2} \frac{m(b\sqrt{1+x^2})}{m(r_t)} dx. \quad (18) $$

This integral verifies $0 \le I(b, r_t) \le 1$ and is equal to 1 when $b > r_t$. Our result in Eq. (17) can be compared to the work of GF83. In their study, the authors of GF83 are considering a galaxy perturbed by another galaxy, both objects having an extended Plummer density profile [75], and they derive an expression in two limiting cases: $b \ll r_t$ and $b \gg r_t$ where $r_t$ is the extension of the subject galaxy. They propose an interpolation between these two asymptotic cases to get an expression for any $r$. Taking the point-like limit for the perturbing galaxy, their expression becomes identical to ours for both $b \ll r$ and $b \gg r_t$. However our expression in Eq. (17) is valid for any values of $b, r$ and $r_t$ hence no interpolation is required. Having computed the velocity gain, we can translate it into an energy gain per unit mass

$$ \delta E = \frac{1}{2} \left[ (\mathbf{v} + \delta \mathbf{v})^2 - \mathbf{v}^2 \right] = \frac{1}{2} (\delta \mathbf{v})^2 + \mathbf{v} \cdot \delta \mathbf{v} \quad (19) $$

Let us discuss the first term on the right-hand side of Eq. (19). We have

$$ (\delta \mathbf{v})^2 = \left( \frac{2 G m_\star}{v_r b} \right)^2 \left[ I^2 + \frac{b^2(1-2I) - 2I \mathbf{r} \cdot \mathbf{b}}{(\mathbf{r} + \mathbf{b})^2 - (\mathbf{r} \cdot \hat{\mathbf{e}}_v)^2} \right] \quad (20) $$

This shows that the energy gain is not spherically distributed around the center of the clump. In the relative velocity direction, $(\delta \mathbf{v})^2$ does not depend on radius and is proportional to $(I - 1)^2$. In every other direction, we have $(\delta \mathbf{v})^2 \propto (I - 1)^2$ where $r \ll b$ and $(\delta \mathbf{v})^2 \propto I^2$ where $r \gg b$. This is illustrated in Fig. 3 where we show $(\delta \mathbf{v})^2$ along several directions in a case where $b \ll r_t$ for a subhalo with an NFW profile. The direction-dependence of $(\delta \mathbf{v})^2$ means the subhalo loses its spherical symmetry after encountering the star. Instead of describing this deformation, we derive a spherically-symmetric approximation for $(\delta \mathbf{v})^2$. The DM particles that encounter the star during at some point in time gain an infinite amount of energy as seen in the direction opposite to $\mathbf{b}$ (dotted-dashed red line in Fig. 3). To regularize this divergence, we replace all cosines by their average on the sphere i.e. $\mathbf{r} \cdot \mathbf{b} \approx 0$ and $(\mathbf{r} \cdot \hat{\mathbf{e}}_v)^2 \approx r^2/3$. The same simplification is performed by GF83 in their calculation. This leads to

$$ (\delta \mathbf{v})^2 \approx \left( \frac{2 G m_\star}{v_r b} \right)^2 \left[ I^2 + \frac{3(1-2I)}{3+2(r/b)^2} \right]. \quad (21) $$

Under this approximation, $(\delta \mathbf{v})^2$ is always finite and only depends on $r$. This solution is shown as the black curve in Fig. 3 along with the GF83 solution in magenta. We see that our solution reproduces the asymptotic behavior at both large and small radii in all directions (at the exception of the large radii behavior along the $\mathbf{v}_r$ direction). On the other hand, the GF83 solution fails at reproducing the correct asymptotic limits in this case. In the opposite case where $b \gg r_t$, our solution agrees with GF83.

It is common practice in the literature to work with the ratio of the total integrated kinetic energy kick of
the subhalo over its binding energy. For comparison we therefore introduce

$$ \delta E_{\text{int}} = 2\pi \int_{0}^{r_t} dr r^2 \rho(r) (\delta \mathbf{v})^2 \quad (22) $$

and the binding energy

$$ U = 4\pi G_N \int_{0}^{r_t} m(r) \rho(r) r dr. \quad (23) $$

The ratio of these two quantities is represented in Fig. 4 with respect to the impact parameter. It scales as $b^{-4}$ when $b \gg r_s$ and as a constant (up to a small logarithmic correction) when $b \ll r_s$. Therefore one recovers the behavior introduced in [76] and used in the context of dark matter subhalos in Refs. [54, 58] (as well as in Ref. [55] in the large $b$ limit):

$$ \frac{\delta E_{\text{int}}}{U} \sim \frac{G_N m_*^2}{v_r^2 \rho_s (\mu b + r_s)^4} \quad (24) $$

with $\mu$ a parameter. Albeit Ref. [54] provides an estimate for $\mu$ it is ill-defined for a NFW profile. Here we find $\mu = 213$ for $r_t/r_s = 10^{-2}$, $\mu = 3.57$ for $r_t/r_s = 1$ and $\mu = 0.228$ for $r_t/r_s = 10^2$. The dash-dotted curve corresponds to the characteristic binding energy introduced in Ref. [58] (referred to as D19) where the author assumes a slightly different shape $\delta E_{\text{int}} \propto 1/(b^4 + r_s^4)$. If the tidal radius is not smaller than the scale radius our solution provides better agreement with Eq. (24) than the GF83 solution.

Let us now discuss the second term on the right-hand side of Eq. (19). Since $\delta\mathbf{v}$ is independent of $\mathbf{v}$, this term averages to zero over the velocity distribution of DM particles inside the subhalo, however it does contribute to a scatter around $\langle\delta E\rangle$. If the velocity distribution at any point in the subhalo follows an isotropic distribution with variance $\langle v^2 \rangle = 3\sigma_{\text{sub}}^2$ then we have

$$ \langle \delta E^2 \rangle - \langle \delta E \rangle^2 = \sigma_{\text{sub}}^2 |\delta \mathbf{v}|^2. \quad (25) $$

The velocity dispersion $\sigma_{\text{sub}}$ can be computed by solving the Jeans equation, see App. B. Eq. (25) shows that the scatter is important compared to the average when $\sigma_{\text{sub}}$ is large compared to $|\delta\mathbf{v}|$. The energy gain is compared to the potential energy in Fig. 5 for a subhalo with parameters $\rho_s = 5.8 \times 10^8 \text{ M}_\odot/\text{kpc}^3$, $r_s = 3.5 \times 10^{-3} \text{ pc}$, $r_t = 10 r_s$ encountering a star of mass $M_* = 1 \text{ M}_\odot$ with an impact parameter $b = 2 \times 10^{-2} \text{ pc}$ and relative velocity $v_r = 200 \text{ km/s}$. This figure shows that the scatter is important at intermediate radii in the clump, while it is small in the central regions and near the edge where the velocity dispersion is small.

## IV. EFFECT OF THE STELLAR POPULATION ON A SINGLE SUBHALO

This section is devoted to the description of the cumulative effect of multiple stellar encounters as a DM subhalo crosses the Galactic disk.

<figure><img src="image_4.png" /><figcaption>FIG. 3. ($\delta\mathbf{v}$)$^2$ along $\mathbf{b}$ (solid red), in the direction opposite to $\mathbf{b}$ (dotted-dashed red), along $\mathbf{v}_r$ (green) and along $\hat{\mathbf{e}}_\perp$ (blue), compared to the prediction by Ref. [59] (magenta) and to the approximation in Eq. (21) (black). The tidal radius of the subhalo is set to $r_t = 100 \times b$. Gaetan: change $M_*$ for $m_*$ in the y-axis label</figcaption></figure>

### A. The stellar population

Given a subhalo crossing the stellar disk, we want to know what is the distribution of impact parameters for stellar encounters. Assuming the disk is an infinite, homogeneous slab with surface density $\Sigma_*$, and the subhalo moves along a straight line which makes an angle $\theta$ with respect to the perpendicular to the disk, then the number of encounters with parameters between $b$ and $b+db$ is

$$ dN = \frac{\Sigma_*}{\bar{m}_*} \frac{2\pi b db}{\cos(\theta)} \quad (26) $$

where $\bar{m}_*$ is the average mass of a disk star. This distribution evidently diverges for $\cos(\theta) = 0$ i.e. orbits within the disk. The infinite and homogeneous assumptions can be dropped to get a finite distribution everywhere however the final expression is not analytical and the computation is much more involved, as shown in App. D. Comparisons shows that Eq. (26) is a good approximation as long as $\cos(\theta) \gg b/R_d$ where $R_d$ is the disk scale radius. This condition is satisfied for the vast majority of subhalos so Eq. (26) is used in the following. Our calculation of the energy gain in Sec. III is not valid for arbitrary impact parameters. It is assumed that the encounter is isolated so the impact parameter must be smaller than the typical distance between stars. To compute this distance, we need a model for the Galactic disk. We use the Milky-Way mass model established by McMillan [60]. In
<figure><img src="image_5.png" /><figcaption>FIG. 4. Ratio of the total kinetic energy integrated on the entire profile over the binding energy for various values of the tidal radius. The solid curves are obtained with our ansatz of Eq. (21), dashed curve are obtained using GF83 result and the dotted lines correspond to a comparison with the usually adopted shape $\propto (\mu b + r_s)^{-4}$. Comparison is made with the characteristic binding energy introduced in [58] (here denoted D19).</figcaption></figure>

this model, two exponential stellar disks (thick and thin) are fitted against a number of observational constraints, along with a DM halo, a stellar bulge and two gaseous disks. The best-fit parameters for the stellar disks are given in Tab. I. Given the axisymmetric mass density of stars $\rho_*(R, z)$, we define the maximal impact parameter

$$b_{\text{max}}(R) = \frac{\int_{-\infty}^{+\infty} \rho_{*} \frac{1}{2} \left( \frac{\rho_{*}}{m_{*}} \right)^{-1/3} dz}{\int_{-\infty}^{+\infty} \rho_{*} dz} \\ = \frac{\overline{m_{*}}^{1/3}}{\Sigma_{*}(R)} \int_{0}^{+\infty} \rho_{*}^{2/3}(R, z) dz. \qquad (27)$$

Using an average mass $\overline{m}_* \simeq 0.17 M_{\odot}$ [77], we find $b_{\text{max}}(8\,\text{kpc}) \simeq 1.1\,\text{pc}$. We also make a straight-line-trajectory assumption when computing $\delta E$. This is reasonable only if the kinetic energy in the center-of-mass frame is much larger than the potential energy $T \gg |W|$, with

$$T = \frac{1}{2} \frac{m_{*} m_{t}}{m_{*} + m_{t}} v_{r}^{2} \qquad (28)$$

and

$$W = -G m_{*} \left[ \frac{m(r_{*})}{r_{*}} + \Theta(r_{t} - r_{*}) \int_{r_{*}}^{r_{t}} \frac{\rho(r)}{r} d^{3} \mathbf{r} \right] (29)$$

<figure><img src="image_6.png" /><figcaption>FIG. 5. Energy gain (red curve) for a subhalo with parameters $\rho_s = 5.8 \times 10^8 M_{\odot}/\text{kpc}^3$, $r_s = 3.5 \times 10^{-3} \text{pc}$, $r_t = 10 r_s$ encountering a star of mass $M_* = 1 M_{\odot}$ with an impact parameter $b = 2 \times 10^{-2} \text{pc}$ and relative velocity $v_r = 200 \text{km/s}$. The red-shaded area shows the $1\sigma$ scatter around $\delta E/m$ due to the inner velocity dispersion. The blue curve shows the gravitational potential defined in Eq. (B1).</figcaption></figure>

$W$ is minimal when $r_* = b$ thus the condition $T \gg |W|$ defines a minimal impact parameter $b_{\text{min}}$. This parameters is shown for several subhalo masses in Fig. 6. We see that $b_{\text{min}}$ is much smaller than $b_{\text{max}}$ unless the relative velocity is smaller than 0.1 km/s. Since the typical velocity of subhalos in the Galaxy is of order 100 km/s, we have $b_{\text{min}} = 0$ is most cases and the assumption of straight-line trajectory is verified. The total number of encounter is

$$\mathcal{N} = \frac{\Sigma_{*}}{\overline{m_{*}}} \frac{\pi}{\cos(\theta)} (b_{\text{max}}^{2} - b_{\text{min}}^{2}). \qquad (30)$$

At 8 kpc, we find $\mathcal{N} \simeq 2346 \times (0.5/\cos(\theta))$.

Also needed is the relative velocity between stars and subhalos. The velocity of subhalos is assumed to follow an isotropic Maxwell Boltzmann distribution with dispersion $\sigma(R)$ that can be computed using Jean's equation – where the baryonic contribution to the potential has been added according to App. B. Moreover, assuming that stars follows circular trajectories at a velocity $v_*(R)$, we get the relative speed distribution

$$f(v_r) = \sqrt{\frac{2}{\pi}} \frac{v_r}{\sigma v_*} \sinh\left(\frac{v_r v_*}{\sigma^2}\right) e^{-(v_*^2 + v_r^2)/(2\sigma^2)} \quad (31)$$

and the mean relative speed

$$\bar{v}_r = \sigma \sqrt{\frac{2}{\pi}} \left\{ e^{-X^2} + \frac{\sqrt{\pi}}{2}(1 + 2X^2) \frac{\text{erf}(X)}{X} \right\} \quad (32)$$
<figure><img src="image_7.png" /><figcaption>FIG. 6. Minimal impact parameter to avoid capture by the star.</figcaption></figure>

where $X = v_*/(\sqrt{2}\sigma)$. At 8 kpc, the relative speed is $\bar{v}_r \simeq 334$ km/s. The last ingredient is the probability distribution of stellar mass $p_{m_*}$ used to compute $\bar{m}_*$ and taken from Ref. [78]. We are now equipped to define the total energy kick received by a subhalo when it crosses the entire stellar disk.

## B. Total energy kick and scatter

When crossing the stellar disk a subhalo encounters $\mathcal{N}$ stars, each with a different impact parameter. A particle inside the subhalo receives a series of velocity kicks $\{\delta\mathbf{v}_i\}_{1\le i\le \mathcal{N}}$. We assume that the subhalo does not have time to relax between encounters, hence the total velocity kick is given as the sum

$$ \Delta\mathbf{v} = \sum_{i=1}^{\mathcal{N}} \delta\mathbf{v}_i. \qquad (33) $$

Similarly to Eq. (19), the associated total energy kick per units of mass is given by

$$ \Delta E = \frac{1}{2} |\Delta \mathbf{v}|^2 + \mathbf{v} \cdot \Delta \mathbf{v}. \qquad (34) $$

Because encounters are characterized by the statistical distribution of impact parameters and stellar masses, all vectors $\delta\mathbf{v}_i$ and $\Delta\mathbf{v}$ are random variables. In the following, we first show how it is possible to evaluate a PDF for $\Delta\mathbf{v}$. The sequence $\{\delta\mathbf{v}_i\}_{1\le i\le \mathcal{N}}$ behaves as a random walk in "velocity space". From Eq. (17) every $\delta\mathbf{v}_i$ is confined

in the same plane, perpendicular to the relative velocity vector. The random walk is thus two-dimensional$^2$. When $\mathcal{N}$ is large enough it can be described as a Brownian motion through the Central Limit (CL) theorem. As the random walk is isotropic, the PDF for $\Delta\mathbf{v}$ is

$$ p_{\Delta\mathbf{v}}(\Delta\mathbf{v}) \approx \frac{1}{\pi \mathcal{N} \delta v^2} e^{-\frac{|\Delta\mathbf{v}|^2}{\mathcal{N} \delta v^2}} \qquad (35) $$

and only depends on the second moment of $\delta\mathbf{v}$,

$$ \overline{\delta v^2} = \int dm_* p_{m_*}(m_*) \int_{b_{\min}}^{b_{\max}} db p_b(b) (\delta\mathbf{v})^2. \qquad (36) $$

with $(\delta\mathbf{v})^2$ given in Eq. (21), $p_b(b) = (d\mathcal{N}/db)/\mathcal{N}$ and $p_{m_*}(m_*)$ taken from Ref. [77]. Note that to gain time in the numerical evaluation we often approximate $p_{m_*}(m_*) \sim \delta_D(m_* - \bar{m}_*)$. A straightforward first evaluation for $\Delta E$ is obtained by considering the average value of $\Delta\mathbf{v}$. It yields

$$ \overline{\Delta E} \sim \frac{1}{2} |\overline{\Delta \mathbf{v}}|^2 = \frac{1}{2} \mathcal{N} \overline{\delta v^2}. \qquad (37) $$

However, one may wonder about the impact of the second term in Eq. (34). In order to properly take it into account the solution is to derive the full PDF for $\Delta E$. The initial velocity distribution being approximated as a Maxwell-Boltzmann with velocity dispersion $\sigma_{\text{sub}}$ – c.f. App. E 3,

$$ p_{\Delta E}(\Delta E) = \frac{\exp\left(\frac{\Delta E}{2\sigma_{\text{sub}}^2} - \frac{|\Delta E|}{2\sigma_{\text{sub}}^2} \frac{\sqrt{1+s^2}}{s}\right)}{4\sigma_{\text{sub}}^2 s\sqrt{1+s^2}} \qquad (38) $$

where $s^2 \equiv \mathcal{N}\overline{\delta v^2}/(8\sigma_{\text{sub}}^2) = \overline{\Delta E}/(4\sigma_{\text{sub}}^2)$ is a normalised ratio of the variance of $|\Delta\mathbf{v}|^2$ (or also the average in that case) and the variance of the initial velocity $|\mathbf{v}|$. The associated scatter is

$$ \sigma_{\Delta E} = \overline{\Delta E} \sqrt{1 + \frac{1}{2s^2}}. \qquad (39) $$

This distribution is plotted in the left panel of Fig. 7 in terms of the associated centred reduced variable. When $\sigma_{\text{sub}}$ is large, $s$ is small and the distribution is symmetric with respect to the average. It shifts toward lower values for small $\sigma_{\text{sub}}$. In order to properly take into account the dispersion in energy kick and the shift in the distribution, one should evaluate a new density profile for the subhalo by properly removing the particle with a final kinetic energy (resp. final velocity) greater than the gravitational potential (resp. escape velocity) – this

<footnote>$^2$As the stars have their own velocity, the relative velocity vector may vary from one encounter to another if the velocity of the subhalo is not high enough. Therefore the perpendicular plane is not fixed and the random walk is not strictly two-dimensional. However for simplicity and because it should not impact the results by order of magnitudes as seen in Fig. 7 we stick to the two-dimensional hypothesis for the results shown in the main text.</footnote>
<figure><img src="image_8.png" /><figcaption>FIG. 7. In both panels the solid curves are obtained by assuming $\Delta v$ constructed out of a 2D isotropic random walk in velocity space and dashed curve are for a 3D isotropic random walk. **Left panel:** Probability density function of the centred reduced total energy kick for different values of the parameter $s$. We denoted by $\sigma_{\Delta E}$ the standard deviation of $\Delta E$. At small values of $s$ the PDF is symmetric around the average. This can be understood as the dominant dispersion comes from the uncertainty on the initial velocity that, because of the Maxwell-Boltzmann distribution, is symmetric. When $s$ grows however, the PDF is peaked on negative values, so that energies lower than the average are more probable. **Right panel:** The maximal energy that at least a fraction $q$ of the particles gain over the average energy kick with respect to $q$ for different value of $s$. The inset is a zoom on the region $0.2 < q < 0.8$ for $s > 0.5$. The fact that at large values of $s$ the PDF of the right panel is no longer symmetric with respect to the average is the reason why $\text{Med}(\Delta E) < \bar{\Delta E}$ in all configurations.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Complementary cumulative distribution function of the energy kick received after one crossing of the disk due to the encounter with stars $\Delta E/m$. In blue is shown the true value (solid) and the approximation using the CL theorem (dashed-dotted) for $b_{\min} \sim 0$. In red are similar curve imposing a lower cut-off on the distribution of impact parameters $b_{\min} \rightarrow b_0(Q = 0.2)$. The vertical green dashed line is the value of the gravitational potential $\psi$. **Left panel:** In the inner part of the subhalo $r/r_s = 0.01$. **Right panel:** In the outskirts of the subhalo $r/r_s = 2.34$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 9. Left panel. Median energy kick received by at least half of the particles in terms of the radius. Upper panel: In blue is shown the true value (solid) and the approximation using the CL theorem (long-sort dashed) for $b_{min} \sim 0$. In red are similar curve imposing a lower cut-off on the distribution of impact parameters $b_{min} \rightarrow b_0(Q = 0.2)$. The green dashed curve is the value of the gravitational potential of the subhalo. Lower panel: Value of the parameter $s$ in the two configuration with (blue) and without (red) the cut-off on the impact parameters. The magenta dotted curve is the limit $s = 0.5$. In addition as the energy kick becomes sizeable when it is of the same order than the gravitational potential it is indeed when $s \gtrsim 1$ that the effect of stars become important. This can be understood easily as the inequality $\sigma < \psi$ obtained from the expression of $\sigma_{sub}$ given by Jean's equation implies directly $\Delta E \sim |\Phi| \Rightarrow s > 0.5$. Right panel. Same figure for the energy kick standard deviation. Note that the dash-dotted and solid red curve are not distinguishable (one has good convergence of the CL theorem when the impact parameter range is truncated for the dispersion).</figcaption></figure>

the impact of one disk crossing on a single subhalo and cannot be applicable in the framework of SL17. In the following we focus on defining one energy kick value as an estimate of the energy kick felt by all particles in each shell.

Let us introduce $\Delta E(q)$ the maximal total energy kick that is received by at least a fraction $q$ of the particles, in a given shell – the full expression is given in App. E3 for a Maxwellian distribution of initial particles velocity. In particular, the median energy gain is $\text{Med}(\Delta E) \equiv \Delta E(q = 0.5)$. One can show that for any value of $s$,

$$ \frac{1}{2} < \frac{\text{Med}(\Delta E)}{\bar{\Delta E}} < \ln(2) \approx 0.69, \quad (40) $$

so mean and median are always close. The ratio $\Delta E(q)/\bar{\Delta E}$ is plotted in the right panel of Fig. 7 for different values of $s$. If $s$ is large enough ($s \gtrsim 0.5$), i.e. if the effect of the encounters is relevant, then $\Delta E(0.25 < q < 0.75)$ and $\bar{\Delta E}$ are always close to each other $\Delta E(q)/\bar{\Delta E} \in [0.1, 2]$. Hence, it is both physically meaningful and convenient to define the energy kick for all particles at a given distance from the center of the subhalo as being $\text{Med}(\Delta E)$ and it can be well-approximated as $\Delta E \sim \kappa \bar{\Delta E}$ with the coefficient $\kappa \sim 0.7$.

So far, all the results are built on the fact that the number of encountered stars per crossing $N$ is large enough to apply the CL theorem as if it were infinite. As a matter of fact, according to the position inside the subhalo, this is not necessarily the case. Indeed, when $b_{min}$ is close to 0, as the velocity kick satisfies $(\delta\mathbf{v})^2 \propto 1/b^4$ for the innermost particles, it gets tremendously large. However, because $p_b(b) \propto b$, in the majority of the disk crossings, on a total of $N \sim 10^2 - 10^5$ encounters, none of them has an impact parameter $b \ll b_{max}$. Therefore when taken blindly the CL theorem overestimates the energy kick felt by the particles in the innermost part of the subhalo. In order to illustrate and quantify this effect, let us look at a striking example and therefore focus on a small subhalo. More precisely consider a typical subhalo before its first crossing of the stellar disk, with a typical mass $m_t = 1.6 \times 10^{-9} \, \text{M}_\odot$, scale radius $r_s = 7.1 \times 10^{-7} \, \text{kpc}$ and tidal radius $r_t = 2.34192 \times r_s$ at a distance $R = 8 \, \text{kpc}$ from the center of the MW. All these values are consistent with a subhalo that has only been smoothly stripped by the potential of the Galaxy. Its relative velocity with the stars is given by the average value $\bar{v}_t(R = 8 \, \text{kpc}) = 334 \, \text{km.s}^{-1}$ and its inclination is given by $\cos\theta = 0.5$ so that it encounters $N \sim 2346$ stars. The goal here is to determine the true PDF for $\Delta E$ and compare it to Eq. (38).
However, even though the number of encounters is not high enough to have proper convergence to the CL distribution, it is still too high to allow a fully analytical computation. Indeed, that would require to evaluate $N$ convolutions of the PDF of $\delta v$ that are not even possible numerically even with the use of Fourier transforms. A Monte-Carlo (MC) algorithm is better suited for the task and with a total of $5 \times 10^5$ draws we achieve convergence to the true PDF of $\Delta E$. In Fig. 8 we represent the complementary cumulative distribution function (CCDF) of $\Delta E$ (as it is more convenient to display than the PDF and contains the same information)

$$ \bar{F}(> \Delta E) \equiv \int_{\Delta E}^{+\infty} p_{\Delta E}(\Delta E') d(\Delta E'), \quad (41) $$

for two radii, one in the inner part of the subhalo $r = 10^{-2} \times r_s$ and one in the outskirts $r = 2.34 \times r_s$. The solid curves show the MC results while the dashed-dotted lines are the CL theorem expectations. As predicted the discrepancy is less pronounced in the outskirts of the subhalo and the approximated distribution is shifted toward much higher values of $\Delta E$ than the true distribution for the innermost particles. In Fig. 9 we show the evolution of the median energy kick and its dispersion in terms of the radius inside the substructure. For the same subhalo, the left panel clearly shows that the expectation of the CL theorem for $\text{Med}(\Delta E)$, in dash blue, overshoots the gravitational potential on the entire range of radii, while the true value, in solid blue, only crosses it on the outskirts, making the CL result unusable.

Unfortunately, even though a MC algorithm is convenient to treat the aforementioned example, it is too greedy in terms of computation time to be used for the study of a full subhalo population. Because only encounters with small impact parameters are responsible for the convergence issue while they have very small chances to occur, a solution consists in truncating the impact parameter range from below. We detail the method in the following. For one crossing of the disk we denote by $b_0$ the minimal impact parameter. From the PDF of impact parameter, the PDF of $b_0$ is given by

$$ p_{b_0}(b_0) = \frac{2N}{(1 - \beta^2)^N} \frac{b_0}{b_{\max}^2} \left[ 1 - \left( \frac{b_0}{b_{\max}} \right)^2 \right]^{N-1} \quad (42) $$

with $\beta \equiv b_{\min}/b_{\max} \ll 1$. We introduce $b_0(Q)$ defined so that only a fraction $Q$ of the disk crossings has an impact parameter lower than $b_0(Q)$. It is given by

$$ b_0(Q) = b_{\max} \left[ 1 - (1 - \beta^2)(1 - Q)^{\frac{1}{N}} \right]^{\frac{1}{2}}. \quad (43) $$

Therefore, the best way to recover a pseudo-convergence to the CL distribution without losing too much information on the real distribution is to enforce that the minimal impact parameter is no longer $b_{\min}$ but $b_0(Q \sim 0.2)$. The value of $Q = 0.2$ is tuned by hand

<figure><img src="image_11.png" /><figcaption>FIG. 10. Evolution with time of the mass fraction (blue) and tidal radius (red) of a subhalo of initial virial mass $8 \times 10^{-6} M_{\odot}$. This mass correspond to the the rescaled virial mass today studied in AZ07 (Ref. [56]) assuming fixed scale density and scale radius. The initial radius of the structure is fixed to $r_{200}$ and the energy kick induced by star encounters at each crossing is fixed by $\Delta E = \kappa \sqrt{\delta v^2}/2$ with $\kappa \sim 0.7$ and $\delta v^2$ given by the integral in Eq. (36). The concentration $c_{200} = 84.7$ correspond to that of AZ07 rescaled to today assuming fixed scale radius and scale density. The concentration $c_{200} = 56.8$ is the median concentration from [79].</figcaption></figure>

in order to find a correct agreement between the CL and the MC results. When $b_{\min} = 0$ one can show that $b_0(Q) \sim 1.6 \times b_{\max}/N$ and scales as the inverse of the number of encountered stars. The outcome of these procedure are shown by the red curves in Fig. 8. Convergence is not exact but the new CL result is now much closer to the true CCDF in solid blue. In the left panel of Fig. 9, the CL result upgraded with the truncation of impact parameter at $b_0(0.2)$ (in dashed red) provides a good estimation of the true median (in solid blue). Nonetheless, the dispersion is not well recovered in the right panel: the true energy kick has a much higher dispersion. This was expected as the non-convergence of the CL comes from the large discrepancy between the median and the average which follows from a large dispersion. Therefore, one should keep in mind that, defining a precise energy kick for the entire population of particle in one shell of a subhalo, is not trivial and the chosen definition could severely impact the results. Hereafter we stick to the definition using the median as it is a well physically grounded prescription. Now, before moving on to results, let us discuss the validity of the impulse approximation.
<figure><img src="image_12.png" /></figure>

FIG. 11. Evolution with mass of the disruption time for a subhalo trapped at 8 kpc of the GC in the stellar disk. Two cases are considered: stars only where the initial radius is the virial radius today (blue) and stars+smouth were the initial radius is the Jacobi radius. Comparison is made with GG07 (Ref. [54]) in green. They work with the fraction of total energy kick over the binding energy showed in Eq. (24) and consider an approximate expression made of two asymptotic behaviours according to $b$ with a transition at $b = b_s$. In order to parameterize some levels of uncertainties they consider two mass behaviour for this parameter, which are represented here.

### C. Validity of the impulse approximation

So far, our calculations have rested on the impulse approximation which states that the motion of the DM particles within the subhalo can be neglected for the duration of the encounter with a star. The validity of this assumption can be checked by comparing the typical orbital timescale of DM particles to the duration of the encounter $t_{\text{col}} \sim b/v_r$. At a given position $r$ within the clump, the orbital frequency is given by $\omega(r) = \sigma_{\text{sub}}(r)/r$ thus the encounter is impulsive if $t_{\text{col}}\omega(r) \ll 1$ everywhere in the subhalo. For a cuspy halo profile like NFW, the orbital frequency diverges at $r=0$ so the impulse approximation necessarily breaks down at some $r > 0$ but this radius might be very small compared to the scale radius $r_s$ of the structure. Indeed, using the maximal impact parameter defined in Eq. (27) and the mean relative speed in Eq. (32), we find that $t_{\text{col}}\omega(10^{-3}r_s) < 1$ regardless of the mass of the subhalo (we fix the concentration using the mass-concentration relation in [79]). At 8 kpc, for a subhalo with infall mass $10^{-10} M_{\odot}$, we find that $t_{\text{col}}\omega(10^{-3}r_s) \simeq 1$ for $v_r \simeq 20 \text{ km/s}$. For a single encounter, using the distribution in Eq. (31), we

<figure><img src="image_13.png" /></figure>

FIG. 12. Ratio of the physical mass $m_t$ over the cosmological virial mass $m_{200}$ with respect to $m_{200}$ for several configurations of tidal effects and three different concentrations. Are considered a scenario when only smooth stripping is included, with no effects from baryons (dotted), a scenario where only the star encounters (resp. disk shocking) are taken into account on top of the smooth stripping (dashed, resp. solid) and a scenario where all effects are included (dash-dotted).

find that the probability for the relative speed to be less than 20 km/s is 0.02%. We conclude that the impulse approximation is valid for the overwhelming majority of encounters down to radii as small as $10^{-3} r_s$.

The situation is quite different when considering the complimentary effect of disk shocking *i.e.* the gravitational shocking induced by the smooth potential of the disk rather than each individual star. In that latter case, the duration of the encounter is the disk-crossing time $t_{\text{cross}} \sim z_d/v_z$ where $z_d$ is the typical scale-height of the disk and $v_z$ the subhalo crossing speed along the perpendicular direction. Because $z_d \gg b_{\text{max}}$, the impulse approximation breaks down more often in this case and adiabatic invariance must be accounted for to get accurate results [80, 81].

### D. Results and comparisons with previous work

In conclusion, the energy kick induced by stars encounter on a total disk crossing is defined as $\Delta E \sim \text{Med}(\Delta E) = \kappa \sqrt{\delta v^2}/2$ with $\kappa \sim 0.7$ and $\delta v^2$ given by the integral in Eq. (36) truncated from below by imposing the cut-off $b_{\text{min}} \rightarrow b_0(Q=0.2)$. The total impact of star encounters during several disk crossings on one subhalo can be evaluated by replacing the value of $\Delta E$ in Eq. (8) by the value of the median discussed above.
Using this procedure we compared, in Fig. 10, the evolution, with time, of the bound mass fraction of a subhalo with initial radius given by the virial radius of the simulation results in Ref. [56] (hereafter AZ07). As they considered a subhalo at $z \sim 26$ with concentration ($c_{vir} = 2$)) and virial mass $m_{vir} = 10^{-6}$ M$_{\odot}$, to use our formalism we roughly rescaled its size to $z \sim 0$ by assuming that the scale radius and scale density remain constants (which should hold approximately [82]). We found a concentration $c_{200} = 84.7$, a mass $m_{200} = 8 \times 10^{-6}$ M$_{\odot}$ and a tidal radius $r_{200} = 0.42$ pc. In the same figure we also show the evolution of the bound mass fraction for a subhalo of same virial mass but median concentration picked in [79]. We consider that the subhalo enters the galactic disk with an inclination $\cos\theta = 0.5$ and has a relative velocity with the stars $\bar{v}_r(8\,\text{kpc}) = 334\,\text{km.s}^{-1}$. Even if stripping becomes less and less efficient with time in our work, all three results are still in good agreement. Moreover, while we have derived an analytical estimate for the number of star encounters per crossing and for the number of crossings, the authors of AZ07 used a refined model of the Galaxy to evaluate these two quantities with better precision. In addition, the discrepancies can be further understood as we assume a sharp and fixed truncation in radius (parameterized by $r_t$) without change of the density profile at every crossings such as it realistically should be [55, 58].

In Fig. 11 we show the time necessary to completely destroy a subhalo with a mass $m_{200}$ (and median concentration) trapped and immobile in the Galactic disk at 8 kpc of the GC. For comparison with Ref. [54] (hereafter GG07) we also considered both NFW and Plummer density profiles for the subhalo and we assumed stars with a relative velocity in the range $v_r \sim (270 \pm 3 \times 25) \, \text{km.s}^{-1}$. We focus on two cases, one where the initial radius of the subhalo is given by its virial radius at $z \sim 0$ and one where the initial radius is set from the smooth stripping. Comparing the first scenario to GG07 we observe that the orders of magnitudes and the general behavior of $t_{dis}$ with the mass are similar. Even if it does not impact the conclusion, let us point out, nevertheless, that the comparison is somewhat biased as they considered subhalos at $z \sim 26$ similarly to AZ07 and therefore there can be a mismatch between the definition of mass and virial radius for the same substructures. Another caveat is that, here, we make the crude approximation that the gravitational potential of the subhalo does not change during the entire time it stays within the disk, which may artificially lower the value of $t_{dis}$. Eventually, we note here that it is quicker to destroy subhalos with a NFW profile and a virial mass $m \lesssim 10^{-4}$ M$_{\odot}$ than it is for a Plummer profile of the same mass.

In the above analysis, we considered resilient subhalos by choosing $\epsilon_t = 10^{-2}$ as disruption parameter. However this choice appears not to have a strong impact as when the number of encountered stars becomes large ($N > 10^4$) there is a sudden transition between two distinct behaviors of the energy kick with the radius owing

to the dependence of $b_{min}$ in $N$: from $\Delta E \propto (r/r_s)^2$, as plotted in Fig. 9, to $\Delta E \propto \text{cst.}$ larger than the gravitational potential. Therefore subhalos experience a rapid transition with $N$ between having a new tidal radius $r_t$ close to the initial boundaries and a complete dispersal of all their particles (even in the inner shells).

## V. IMPACT OF STARS ON THE SUBHALO POPULATION

In this section, we incorporate the effect of stars into the SL17 model in addition to smooth stripping and disk shocking. In a first step we describe how we compute the combined effect of individual encounters and disk shocking. Then, in a second step, we show our results for the impact on the subhalo mass function and the total number density.

### A. Combination of the different stripping effects

Let us now discuss the inclusion of the gravitational shocking by individual stars into that framework. In the impulsive approximation limit, the total energy kick is the combination $\Delta\mathbf{v} + \Delta\mathbf{v}_d$ for which a PDF could be formally derived. Because of adiabatic corrections in disk shocking, however, this is not possible. Nevertheless one can always write the total energy gain as

$$ \begin{aligned} \Delta E_{\text{tot}} &= \frac{1}{2}(\Delta\mathbf{v} + \Delta\mathbf{v}_d)^2 + \mathbf{v} \cdot (\Delta\mathbf{v} + \Delta\mathbf{v}_d) \\ &= \Delta E + \Delta E_d + \Delta\mathbf{v} \cdot \Delta\mathbf{v}_d. \end{aligned} \quad (44) $$

In order to circumvent our ignorance of the true distribution for $\Delta\mathbf{v}_d$ with adiabatic corrections we make the assumption that $\Delta\mathbf{v} \cdot \Delta\mathbf{v}_d \sim 0$. This approximation is well-justified in the case of a subhalos having a normal incidence with respect to the disk. Indeed, in the 2D random walk scenario, $\Delta\mathbf{v}$ is then parallel to the disk and $\Delta\mathbf{v}_d$ is perpendicular to it. Therefore one approximates the total energy kick as

$$ \Delta E_{\text{tot}} \approx \frac{1}{2} \left[ 0.7 N \overline{\delta v^2} + |\Delta \mathbf{v}_d|^2 A_1(\eta_d) \right] \quad (45) $$

More details on the total distribution of $\Delta E_{tot}$ are given in App. F and support this definition. Moreover, for efficiency, $(\overline{\delta v})^2$ is evaluated assuming a typical subhalo entering the galactic disk with an average inclination $\cos\theta = 1/2$ and an average relative velocity with the stars $\bar{v}_r(R)$ given in Eq. (32). In the following, we recompute the tidal radius of subhalos by replacing the value of the kinetic energy kick in Eq. (8), that only takes into account the smooth and disk shocking effect, by the above definition. We then show how tidal stripping is impacted by the single encounters with stars.
<figure><img src="image_14.png" /><figcaption>FIG. 13. Mass functions taking into account four different stripping configurations, at different distances from the GC, $R=1$, 8, 15 kpc and for two different values of the mass index $\alpha = 2.0$ on the left and $\alpha = 1.9$ on the right. The minimal cosmological mass is set to $m_{200}^{\min} = 10^{-10} M_{\odot}$. The red (resp. blue) curves correspond to the resilient (resp. fragile) subhalos configuration. The green curve is the unevolved mass function where $r_t = r_{200}$. Notice that at $R=1$ kpc there are no solid and dashed blue curves as disk shocking effects destroy all the population of subhalos if they are fragile.</figcaption></figure>

## B. Results

We consider four different configurations for the evaluation of the tidal effects. (i) *smooth only*: the tidal radius is entirely defined by Eq. (3). (ii) *smooth+stars*: on top of the smooth effect only the individual encounters with stars is included. (iii) *smooth+disk*: on top of the smooth effect only the disk shocking effect is included. (iv) *smooth+stars+disk*: all effects are taken into

account. In Fig. 12 we show the evolution of the final tidal mass in terms of the original cosmological mass for the different stripping configurations and for three concentrations at a distance $R = 8$ kpc from the GC. The dominant effect of baryons on small subhalos with initial mass $m_{200} \lesssim 1 M_{\odot}$ are the individual encounters. On the contrary, for larger subhalos baryonic stripping is mostly due to the disk shocking.

The total mass function for the different stripping cases and for different input parameters of the subhalo model
<figure><img src="image_15.png" /><figcaption>FIG. 14. Number density of subhalos taking into account four different stripping configurations, at different distances from the GC, $R=1, 8, 15$ kpc and for two different values of the mass index $\alpha = 2.0$ on the left and $\alpha = 1.9$ on the right. The minimal cosmological mass is set to $m_{200}^{\min} = 10^{-10}$ M$_{\odot}$. The red (resp. blue) corresponds to the resilient (resp. fragile) subhalos configuration. The green curve is the cosmological number density without tidal effects. Notice that at $R=1$ kpc there is no blue circle marker as disk shocking effect destroy all the population of subhalo if they are fragile.</figcaption></figure>

are plotted in Fig. 13 for a minimal cosmological mass $m_{200}^{\min} = 10^{-10}$ M$_{\odot}$. While baryonic effects are mild and even negligible in the outer region of the disk, e.g. at a distance $R = 15$ kpc, they have more and more impact toward the GC. At $R = 1$ kpc the mass functions are strongly suppressed and offset toward small masses, especially because of stellar encounters. Indeed, in the resilient subhalo scenario (with a disruption parameter $\epsilon_r = 10^{-2}$), they reduce by 6 order or magnitude the mass function for $10^{-10} < m_t \lesssim 10^{-6}$ M$_{\odot}$ (in comparison to the smooth only case) and populate the mass range much below the minimal cosmological mass. The disk shocking effects only produce an equivalent reduction of 4 orders of magnitude. In the fragile case ($\epsilon_r = 1$), disk shocking effects disrupt all subhalos and so do stellar encounters at low masses. At $R = 8$ kpc one notices a similar effect with an almost 2 (resp. 4) order magnitude suppression due to stellar encounters and a 1 (resp. 2) order of magnitude suppression due to disk shocking at low masses for resilient subhalos (resp. fragile). The causes for the strength of stellar encounters in the center are two-fold: close to the GC subhalos cross the disk more often and the stellar density is higher, reducing the interstellar distances and the impact parameters, therefore enhancing the kinetic energy kicks.

In Fig. 14 we show the number density of subhalo at every position from the GC. As seen above, stellar encounters strongly impact the low mass range that is also

the most populated. Consequently, the effects can be important. We recover the conclusions that further than 12 kpc from the GC the baryonic effects are negligible. At 8 kpc effects are already sizeable. Although the number of subhalos in a resilient population is not impacted, fragile subhalos have a number density reduced by one order of magnitude in comparison to the scenario with disk shocking and smooth stripping only (as seen with the dashed and solid curves). The difference grows toward the GC and at $R = 1$ kpc fragile subhalos are all destroyed and resilient subhalos have a population divided by 100 because of star encounters.

## VI. DISCUSSION AND CONCLUSION

In this work we have refined the theoretical analysis of dynamical effects induced on DM particles inside a subhalo during its encounter with a star. In particular we have derived a new solution for the estimated velocity kick received by particles in every radial shell of the subhalo. We have then studied the impact of successive stellar encounters for a subhalo crossing the stellar disk. We have given a specific attention to the cumulative energy kick received by particles in every shells of the subhalo, and we have performed Monte-Carlo simulations that we have compared to our analytical estimates. Simulations
show that a careful treatment of the impact parameter distribution is necessary in order not to overestimate the kinetic energy, as very close but very unlikely encounters tend to dominate average quantities. They also evidence that the energy kick distribution is broad. This means to the median energy kick might not be a very reliable estimate of the energy received by a single subhalo, however it can still be used to gauge the effect on the overall MW subhalo population. We have computed the MW subhalo mass function using the SL17 analytical subhalo model and showed that stellar encounters have a sizeable effect in the inner 10 kpc on subhalos with mass $\lesssim 1 M_\odot$. This mass selection differs from other tidal interactions (smooth tides and disk shocking) which strip subhalos based on their concentration. We have accounted for theoretical uncertainties in the tidal disruption of subhalos. If subhalos are fragile, as found in cosmological simulations, their number density is strongly depleted by stars. On the other hand, if subhalos are resilient, as suggested by theoretical arguments and dedicated numerical studies [40, 64, 65], the mass function is shifted to lower masses.

A caveat of our calculation is the assumption that the shape of the subhalo density profile is preserved and stars simply induce a sharp truncation at the tidal radius. In fact we expect relaxation to modify the internal structure after each encounter [58, 63, 83]. Nevertheless qualitative comparisons with previous numerical studies on stellar encounters seem to be in agreement.

We also describe stellar encounters and disk shocking as if they were independent effects while in fact they are two sides of the same coin. Disk shocking accounts for the average disk potential while star shocking accounts for the disk granularity. A more detailed picture can be achieved through the stochastic formalism [84–88].

However, these approaches are more involved and less straightforward to incorporate in our semi-analytical subhalo population model.

The effect of stars can have a substantial impact on the prospects for (local) DM searches. For instance, the probability that a subhalo passes through the Earth and enhances the local density by a non-negligible amount changes, which may be important for DM direct detection experiments [37] or even play a role for a possible detection of the smallest remnants with space-based laser interferometers [89]. Assuming that DM self-annihilates, the presence of the subhalos boosts the local annihilation rate [21, 26, 27, 29, 32]. Their intense stripping due to stars should also modify the amplitude of this boost. In addition, the stripped DM should form dark streams [57]. A back of the envelope computation shows that, for $m_{\min} = 10^{-10} M_\odot$, $10^3$ to $10^5$ subhalos have crossed the Solar System and may have formed streams in the last 10 Gyr. It is therefore highly likely that the solar system be surrounded by a large amount of them. If detectable, they would give an interesting probe of the DM fine-grained structuring. Eventually, one could also consider the heating of stars as a potential signature of the presence of subhalos [90–92].

## ACKNOWLEDGMENTS

This work has been supported by funding from the ANR project ANR-18-CE31-0006 (GaDaMa), from the national CNRS-INSU programs PNHE and PNCG, and from European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement N° 860881-HIDDeN — in addition to recurrent funding by CNRS-IN2P3 and the University of Montpellier.

[1] P. J. E. Peebles, ApJL 263, L1 (1982).

[2] G. R. Blumenthal, S. M. Faber, J. R. Primack, and M. J. Rees, Nature (London) 311, 517 (1984).

[3] W. H. Press and P. Schechter, Astrophys. J. 187, 425 (1974).

[4] J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser, Astrophys. J. 379, 440 (1991).

[5] C. Lacey and S. Cole, MNRAS 262, 627 (1993).

[6] H. Mo, F. C. van den Bosch, and S. White, Galaxy Formation and Evolution (Cambridge University Press, 2010).

[7] L. Gao, S. D. M. White, A. Jenkins, F. Stoehr, and V. Springel, MNRAS 355, 819 (2004), astro-ph/0404589.

[8] J. Diemand, M. Kuhlen, and P. Madau, Astrophys. J. 657, 262 (2007), astro-ph/0611370.

[9] V. Springel, J. Wang, M. Vogelsberger, A. Ludlow, A. Jenkins, A. Helmi, J. F. Navarro, C. S. Frenk, and S. D. M. White, MNRAS 391, 1685 (2008), arXiv:0809.0898.

[10] R. E. Angulo, C. G. Lacey, C. M. Baugh, and C. S. Frenk, MNRAS 399, 983 (2009), arXiv:0810.2177 [astro-ph].

[11] J. S. Bullock and M. Boylan-Kolchin, ARA&A 55, 10.1146/annurev-astro-091916-055313 (2017), arXiv:1707.04256.

[12] W. J. G. de Blok, Advances in Astronomy 2010, 789293 (2010), arXiv:0910.3538.

[13] K. A. Oman, J. F. Navarro, A. Fattahi, C. S. Frenk, T. Sawala, S. D. M. White, R. Bower, R. A. Crain, M. Furlong, M. Schaller, J. Schaye, and T. Theuns, MNRAS 452, 3650 (2015), arXiv:1504.01437.

[14] F. Lelli, S. S. McGaugh, J. M. Schombert, and M. S. Pawlowski, Astrophys. J. 836, 152 (2017), arXiv:1610.08981.

[15] J. Zavala and C. S. Frenk, Galaxies 7, 81 (2019), arXiv:1907.11775 [astro-ph.CO].

[16] W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000), astro-ph/0003365.

[17] D. N. Spergel and P. J. Steinhardt, Phys. Rev. Lett. 84,
3760 (2000), astro-ph/9909386.

[18] K. Van Tilburg, A.-M. Taki, and N. Weiner, J. Cosmology Astropart. Phys. 7, 041 (2018), arXiv:1804.01991.

[19] J. A. Dror, H. Ramani, T. Trickle, and K. M. Zurek, Phys. Rev. D 100, 023003 (2019), arXiv:1901.04490.

[20] J. Silk and A. Stebbins, Astrophys. J. 411, 439 (1993).

[21] L. Bergström, J. Edsjö, P. Gondolo, and P. Ullio, Phys. Rev. D 59, 043506 (1999), astro-ph/9806072.

[22] J. Lavalle, Q. Yuan, D. Maurin, and X.-J. Bi, A&A 479, 427 (2008), arXiv:0709.3634.

[23] S. Ando, Phys. Rev. D 80, 023520 (2009), arXiv:0903.4685 [astro-ph.CO].

[24] M. R. Buckley and D. Hooper, Phys. Rev. D 82, 063501 (2010), arXiv:1004.1644 [hep-ph].

[25] T. Ishiyama, J. Makino, and T. Ebisuzaki, ApJL 723, L195 (2010), arXiv:1006.3392 [astro-ph.CO].

[26] L. Pieri, J. Lavalle, G. Bertone, and E. Branchini, Phys. Rev. D 83, 023518 (2011), arXiv:0908.0195 [astro-ph.HE].

[27] J. Lavalle and P. Salati, Comptes Rendus Physique 13, 740 (2012), arXiv:1205.1004 [astro-ph.HE].

[28] A. Berlin and D. Hooper, Phys. Rev. D 89, 016014 (2014), arXiv:1309.0525 [hep-ph].

[29] R. Bartels and S. Ando, Phys. Rev. D 92, 123508 (2015), arXiv:1507.08656.

[30] M. Stref and J. Lavalle, Phys. Rev. D 95, 063003 (2017), arXiv:1610.02233.

[31] F. Calore, V. De Romeri, M. Di Mauro, F. Donato, and F. Marinacci, Phys. Rev. D 96, 063009 (2017), arXiv:1611.03503 [astro-ph.HE].

[32] S. Ando, T. Ishiyama, and N. Hiroshima, Galaxies 7, 68 (2019), arXiv:1903.11427.

[33] M. Hütten, M. Stref, C. Combet, J. Lavalle, and D. Maurin, Galaxies 7, 60 (2019), arXiv:1904.10935 [astro-ph.HE].

[34] J. Coronado-Blázquez, M. A. Sánchez-Conde, A. Domínguez, A. Aguirre-Santaella, M. Di Mauro, N. Mirabal, D. Nieto, and E. Charles, J. Cosmology Astropart. Phys. 2019, 020 (2019), arXiv:1906.11896 [astro-ph.HE].

[35] G. Facchinetti, J. Lavalle, and M. Stref, arXiv e-prints, arXiv:2007.10392 (2020), arXiv:2007.10392 [astro-ph.HE].

[36] A. M. Green, Phys. Rev. D 66, 083003 (2002), arXiv:astro-ph/0207366 [astro-ph].

[37] A. Ibarra, B. J. Kavanagh, and A. Rappelt, J. Cosmology Astropart. Phys. 2019, 013 (2019), arXiv:1908.00747 [hep-ph].

[38] F. C. van den Bosch, MNRAS 468, 885 (2017), arXiv:1611.02657.

[39] F. C. van den Bosch and G. Ogiya, MNRAS 475, 4066 (2018), arXiv:1801.05427.

[40] F. C. van den Bosch, G. Ogiya, O. Hahn, and A. Burkert, MNRAS 474, 3043 (2018), arXiv:1711.05276 [astro-ph.GA].

[41] S. Hofmann, D. J. Schwarz, and H. Stöcker, Phys. Rev. D 64, 083507 (2001), arXiv:astro-ph/0104173 [astro-ph].

[42] A. M. Green, S. Hofmann, and D. J. Schwarz, J. Cosmology Astropart. Phys. 8, 003 (2005), astro-ph/0503387.

[43] T. Bringmann and S. Hofmann, J. Cosmology Astropart. Phys. 4, 016 (2007), hep-ph/0612238.

[44] T. Bringmann, New Journal of Physics 11, 105027 (2009), arXiv:0903.0189 [astro-ph.CO].

[45] V. Berezinsky, V. Dokuchaev, and Y. Eroshenko, Phys. Rev. D 68, 103003 (2003), astro-ph/0301551.

[46] F. C. van den Bosch, G. Tormen, and C. Giocoli, MNRAS 359, 1029 (2005), astro-ph/0409201.

[47] J. Peñarrubia and A. J. Benson, MNRAS 364, 977 (2005), arXiv:astro-ph/0412370 [astro-ph].

[48] A. R. Zentner, A. A. Berlind, J. S. Bullock, A. V. Kravtsov, and R. H. Wechsler, Astrophys. J. 624, 505 (2005), astro-ph/0411586.

[49] A. J. Benson, New A 17, 175 (2012), arXiv:1008.1786.

[50] J. Zavala and N. Afshordi, MNRAS 441, 1317 (2014), arXiv:1308.1098.

[51] V. S. Berezinsky, V. I. Dokuchaev, and Y. N. Eroshenko, Physics Uspekhi 57, 1 (2014), arXiv:1405.2204 [astro-ph.HE].

[52] N. Hiroshima, S. Ando, and T. Ishiyama, Phys. Rev. D 97, 123002 (2018), arXiv:1803.07691.

[53] A. Schneider, L. Krauss, and B. Moore, Phys. Rev. D 82, 063525 (2010), arXiv:1004.5432 [astro-ph.GA].

[54] A. M. Green and S. P. Goodwin, MNRAS 375, 1111 (2007), astro-ph/0604142.

[55] T. Goerdt, O. Y. Gnedin, B. Moore, J. Diemand, and J. Stadel, MNRAS 375, 191 (2007), astro-ph/0608495.

[56] G. W. Angus and H. Zhao, MNRAS 375, 1146 (2007), astro-ph/0608580.

[57] H. Zhao, D. Hooper, G. W. Angus, J. E. Taylor, and J. Silk, Astrophys. J. 654, 697 (2007), astro-ph/0508215.

[58] M. Sten Delos, Phys. Rev. D 100, 083529 (2019), arXiv:1907.13133.

[59] O. E. Gerhard and S. M. Fall, MNRAS 203, 1253 (1983).

[60] P. J. McMillan, MNRAS 465, 76 (2017), arXiv:1608.00971.

[61] M. Stref, T. Lacroix, and J. Lavalle, Galaxies 7, 65 (2019), arXiv:1905.02008.

[62] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 462, 563 (1996), astro-ph/9508025.

[63] E. Hayashi, J. F. Navarro, J. E. Taylor, J. Stadel, and T. Quinn, Astrophys. J. 584, 541 (2003), astro-ph/0203004.

[64] R. Errani and J. Peñarrubia, MNRAS 491, 4591 (2020), arXiv:1906.01642.

[65] R. Errani and J. F. Navarro, MNRAS 505, 18 (2021), arXiv:2011.07077 [astro-ph.GA].

[66] J. Diemand, B. Moore, and J. Stadel, MNRAS 352, 535 (2004), astro-ph/0402160.

[67] J. Diemand, M. Kuhlen, P. Madau, M. Zemp, B. Moore, D. Potter, and J. Stadel, Nature (London) 454, 735 (2008), arXiv:0805.1244.

[68] Á. Moliné, M. A. Sánchez-Conde, S. Palomares-Ruiz, and F. Prada, MNRAS 466, 4974 (2017), arXiv:1603.04057.

[69] J. Binney and S. Tremaine, Galactic Dynamics, 2nd ed., Princeton series in astrophysics (Princeton University Press, Princeton, NJ USA, 2008., 2008).

[70] G. Tormen, A. Diaferio, and D. Syer, MNRAS 299, 728 (1998), astro-ph/9712222.

[71] J. P. Ostriker, L. Spitzer, Jr., and R. A. Chevalier, ApJL 176, L51 (1972).

[72] O. Y. Gnedin and J. P. Ostriker, Astrophys. J. 513, 626 (1999), astro-ph/9902326.

[73] J. H. Jeans, Philosophical Transactions of the Royal Society
ciety of London Series A **199**, 1 (1902)

[74] L. Spitzer, Jr., Astrophys. J. **127**, 17 (1958).

[75] H. C. Plummer, MNRAS **71**, 460 (1911).

[76] B. Moore, ApJL **413**, L93 (1993), astro-ph/9306004.

[77] G. Chabrier, ApJL **586**, L133 (2003), arXiv:astro-ph/0302511 [astro-ph].

[78] G. Chabrier, PASP **115**, 763 (2003), arXiv:astro-ph/0304382 [astro-ph].

[79] M. A. Sánchez-Conde and F. Prada, MNRAS **442**, 2271 (2014), arXiv:1312.1729.

[80] M. D. Weinberg, AJ **108**, 1398 (1994), astro-ph/9404015.

[81] M. D. Weinberg, AJ **108**, 1403 (1994), astro-ph/9404016.

[82] B. Diemer and M. Joyce, Astrophys. J. **871**, 168 (2019), arXiv:1809.07326.

[83] J. Peñarrubia, J. F. Navarro, and A. W. McConnachie, Astrophys. J. **673**, 226 (2008), arXiv:0708.3087.

[84] S. Chandrasekhar, Astrophys. J. **94**, 511 (1941).

[85] S. Chandrasekhar, *Principles of stellar dynamics* (Chicago, Ill., The University of Chicago press [1942], 1942).

[86] H. E. Kandrup, Phys. Rep. **63**, 1 (1980).

[87] J. Peñarrubia, MNRAS **474**, 1482 (2018), arXiv:1710.06443.

[88] J. Peñarrubia, MNRAS **484**, 5409 (2019), arXiv:1901.11536 [astro-ph.GA].

[89] A. W. Adams and J. S. Bloom, arXiv e-prints, astro-ph/0405266 (2004), arXiv:astro-ph/0405266 [astro-ph].

[90] R. G. Carlberg, Astrophys. J. **748**, 20 (2012), arXiv:1109.6022.

[91] R. Feldmann and D. Spolyar, MNRAS **446**, 1000 (2015), 1310.2243.

[92] M. Petač, arXiv e-prints, arXiv:1910.02492 (2019), arXiv:1910.02492 [astro-ph.GA].

[93] The Planck Collaboration et al., A&A **641**, A6 (2020), arXiv:1807.06209 [astro-ph.CO].

[94] G. L. Bryan and M. L. Norman, Astrophys. J. **495**, 80 (1998), astro-ph/9710107.

[95] Q. Zhu, L. Hernquist, F. Marinacci, V. Springel, and Y. Li, MNRAS **466**, 3876 (2017), arXiv:1701.05933.

[96] J. Diemand, M. Kuhlen, and P. Madau, Astrophys. J. **649**, 1 (2006), astro-ph/0603250.

[97] C. Lacey and S. Cole, MNRAS **271**, 676 (1994), arXiv:astro-ph/9402069 [astro-ph].

[98] R. K. Sheth, H. J. Mo, and G. Tormen, MNRAS **323**, 1 (2001), astro-ph/9907024.

[99] E. Bertschinger, Phys. Rev. D **74**, 063509 (2006), astro-ph/0607319.

[100] R. H. Wechsler, J. S. Bullock, J. R. Primack, A. V. Kravtsov, and A. Dekel, Astrophys. J. **568**, 52 (2002), astro-ph/0108151.

[101] A. V. Macciò, A. A. Dutton, F. C. van den Bosch, B. Moore, D. Potter, and J. Stadel, MNRAS **378**, 55 (2007), astro-ph/0608157.

[102] A. V. Macciò, A. A. Dutton, and F. C. van den Bosch, MNRAS **391**, 1940 (2008), arXiv:0805.1926.

[103] Y. P. Jing, Astrophys. J. **535**, 30 (2000), astro-ph/9901340.

[104] J. S. Bullock, A. Dekel, T. S. Kolatt, A. V. Kravtsov, A. A. Klypin, C. Porciani, and J. R. Primack, Astrophys. J. **555**, 240 (2001), astro-ph/0011001.

[105] T. Lacroix, M. Stref, and J. Lavalle, J. Cosmology Astropart. Phys. **9**, 040 (2018), arXiv:1805.02403 [astro-ph.GA].
Appendix A: The SL17 subhalo population model: a statistical semi-analytical model

The SL17 model [30] was motivated by the need to build a global galactic mass model including both a smooth DM component and a subhalo population aside from baryonic components, easy to make consistent with potential observational constraints, and in which tidal effects would be calculated from the very components of the model itself. The main constraint that was imposed from the beginning was that the sum of the smooth DM density profile and of the smoothed overall density profile of the subhalo component should give the global DM halo profile, i.e. the one that can be constrained from observational data (see e.g. Ref. [60]).

In the SL17 model the mass density of each subhalo is described by its inner profile. In this study we use a standard NFW profile if not said otherwise. We also consider a Plummer [75] profile for some applications. The mass density at a distance r from the center of the subhalo can be parameterized, in both cases, under the form

$$ \rho(x \equiv r/r_s) = \rho_s x^{-\gamma_{\rho}} [1 + x^{\alpha_{\rho}}]^{\frac{\gamma_{\rho}-\beta_{\rho}}{\alpha_{\rho}}}, \quad (A1) $$

with $r_s$ and $\rho_s$ the scale radius and the scale density respectively, and $x$ the dimensionless radius. In the NFW case, $(\alpha_{\rho}, \beta_{\rho}, \gamma_{\rho}) = (1, 3, 1)$ while in the Plummer case $(\alpha_{\rho}, \beta_{\rho}, \gamma_{\rho}) = (2, 5, 0)$. Henceforth, a subhalo is characterised by three quantities, $r_s$, $\rho_s$ as well as its distance $R$ from the Galaxy Center (GC) – circular orbits are assumed. Conveniently, it is also possible to describe the profile from cosmological parameters: the virial mass and concentration. The virial mass, denoted $m_{\Delta}$, corresponds to the mass contained inside a radius $r_{\Delta}$ over which the subhalo has an average density equal to $\Delta$ times the critical density $\rho_{\text{crit}} = 3H_0^2/(8\pi G)$, with $H_0 = 67.4$ km/s/Mpc, the Hubble parameter [93]. This yields $m_{\Delta} = (4/3)\pi r_{\Delta}^2 \Delta \rho_{\text{crit}}$. The concentration is defined by $c_{\Delta} = r_{\Delta}/r_s$, and there is a one to one relationship between the couples $(m_{\Delta}, c_{\Delta})$ and $(\rho_s, r_s)$. In practice the value $\Delta = 200$ is used [94] as it is a good approximation for the critical over-density of subhalos when they virialize, in the matter-dominated Universe.

Owing to several dynamical effects, subhalos are tidally pruned. Their physical tidal extension is not defined by the cosmological size $r_{200}$ they would have in a flat background, but by their tidal radius $r_t$. According to results of cosmological simulations [9, 38, 63, 66, 67, 70] we expect subhalos that are stripped too much (i.e. that have too small a tidal radius) to be destroyed. In the model, this is implemented by the criterion

$$ \begin{cases} x_t = r_t/r_s \ge \epsilon_t \quad \Rightarrow & \text{the subhalo survives} \\ x_t = r_t/r_s < \epsilon_t \quad \Rightarrow & \text{the subhalo is disrupted} \end{cases} \quad (A2) $$

that relies on the value of $\epsilon_t$, treated as a fixed constant input. The lower this coefficient is, the more resilient subhalos are to tidal stripping. Two values are considered in the following: $\epsilon_t = 1$ in agreement with cosmological simulation where subhalos are rather fragile, and $\epsilon_t = 10^{-2}$ following the semi-analytical studies of [39, 40, 64, 65] where cuspy subhalos are shown to be instead very resilient to tides.

The SL17 model does not only describe individual sub-halos but their entire population using a joint probability distribution function (PDF) on the virial mass, concentration and position of all subhalos. Assuming that all clumps are independent from each other this global PDF can be factorized into $N_{\text{tot}}$ one-point PDFs, with $N_{\text{tot}}$ being the total number of surviving subhalos. The value of $N_{\text{tot}}$ is normalised consistently against DM only, numerical simulations (more precisely the Via Lactea DM only results [67]). The one point PDF is given as

$$ p_t(m_{200}, c_{200}, R) \equiv \frac{1}{K_t} p_R(R) p_m(m_{200}) \times p_c(c_{200}|m_{200}) \Theta \left[ \frac{r_t}{r_s}(m_{200}, c_{200}, R) - \epsilon_t \right] \quad (A3) $$

where $K_t$ is a normalisation parameter to have a probability of one if integrated on the entire parameter space. The PDF for the position $p_R$ is obtained considering that the subhalo spatial distribution follows the global profile of the total DM halo. The PDF for the mass is obtained trough cosmological mass function. Cosmological simulation exhibit power-law mass functions $\propto m^{-\alpha}$ with a mass index $\alpha \lesssim 2$ [8, 9, 67, 95, 96]. This is theoretically backed-up by the Press-Schechter formalism and its extension [3–5, 97, 98], even if the small mass range is still weakly constrained today. Therefore we set $p_m(m_{200}) \propto m_{200}^{-\alpha}$ with $\alpha \in [1.9, 2.0]$. Besides, the virial mass must be bounded from below by $m_{200}^{\min}$, here set as a free parameter of the model. Within a thermal DM particle model the minimal mass is fixed by kinetic decoupling in the Early Universe [41–44, 99] and can go down to $10^{-12}$ M$_{\odot}$. Eventually, the PDF for the concentration is a log-normal [100–104], whose median, given in Ref. [79], is fitted against numerical simulations. The Heaviside function $\Theta$, which encodes the subhalo disruption, leads to the entanglement of the latter three PDFs
through the dependency of $r_t$ on $m_{200}$, $c_{200}$ and R. The full PDF $p_t$ is therefore a complicated, non-separable function.

With this formalism it is possible to describe more precisely the decomposition of the density of the total DM halo $\rho_{\text{tot}}$, as the sum

$$ \rho_{\text{tot}}(R) = \rho_{\text{sm}}(R) + \int dm_t m_t \frac{dn_{\text{sub}}}{dm_t}(R) \qquad (A4) $$

where $\rho_{\text{sm}}$ the density of DM in the smooth component of the halo and the integral corresponds to the contribution of subhalos. The function $dn_{\text{sub}}/dm_t$ is the local evolved (i.e. after stripping) subhalo mass function. It is related to the subhalo PDF through

$$ \frac{dn_{\text{sub}}}{dm_t}(R) = N_{\text{tot}} \int dc_{200} \, dm_{200} \, p_t(m_{200}, c_{200}, R) \times \delta_D[m_t - m[r_t(m_{200}, c_{200}, R)]], \qquad (A5) $$

where $m(r_t)$ is the mass of the subhalo within the tidal radius $r_t$. Note that $m_t$ is the tidal mass of the subhalo i.e. its physical mass as opposed to its cosmological virial mass $m_{200}$. Tidal effects come into play through $p_t$ (which encodes whether subhalos are destroyed or not) and directly through the tidal radius $r_t$. The next section is dedicated to describe the evaluation of $r_t$ at different radii, masses and concentrations as implemented in SL17.

## Appendix B: Velocity dispersion

To compute the velocity dispersion $\sigma$, we start from the Jeans equation for a spherical system

$$ \frac{1}{\rho} \frac{\partial(\rho \langle v_r^2 \rangle)}{\partial r} + 2 \frac{\beta}{r} \langle v_r^2 \rangle = -\frac{d\Phi}{dr}, \quad \text{where} \quad \Phi(r) = -\int_r^{r_t} \frac{G m(r')}{r'^2} dr' \qquad (B1) $$

is the gravitational potential with $r_t$ the tidal radius of the subhalo, and

$$ \beta(r) \equiv 1 - \frac{\langle v_{\theta}^2 \rangle + \langle v_{\phi}^2 \rangle}{2 \langle v_r^2 \rangle} \qquad (B2) $$

is the anisotropy parameter. Here isotropy is assumed, therefore we have $\beta = 0$ and $\sigma_{\text{sub}}^2 = \langle v_r^2 \rangle$ thus

$$ \sigma_{\text{sub}}^2(r) = \frac{G}{\rho} \int_r^{r_t} \frac{\rho(r') m(r')}{r'^2} dr'. \qquad (B3) $$

The same approach is used to compute the velocity distribution of subhalos in the dark halo of the Galaxy. One major difference is that baryons now contribute to the potential $\Phi$ and the mass of the system is

$$ m_{\text{tot}}(R) = m_{\text{DM}}(R) + \int_{|\mathbf{r}'|<R} \rho_b(\mathbf{r}') d^3\mathbf{r}'. \qquad (B4) $$

where $\rho_b$ is the baryonic mass density, which is axisymmetric rather than spherical. The DM velocity variance is then

$$ \sigma^2(R) = \frac{G}{\rho} \int_R^{R_{\text{max}}} \frac{\rho(R') m_{\text{tot}}(R')}{R'^2} dR', \qquad (B5) $$

where the radial extension of the dark halo is fixed to $R_{\text{max}} = 500$ kpc.

## Appendix C: The stellar disks

We use the mass model of McMillan [60] where two stellar disks with density

$$ \rho_d(R, z) = \frac{\Sigma_0}{2 z_d} \exp \left( -\frac{R}{R_d} - \frac{|z|}{z_d} \right) \qquad (C1) $$

are fitted against a number of observational constraints. The best-fit parameters are shown in Tab. I
<table><caption>TABLE I. Stellar disk parameters.</caption><thead><tr><th></th><th>Σ<sub>0</sub><br>[10<sup>8</sup> M<sub>☉</sub>/kpc<sup>2</sup>]</th><th>R<sub>d</sub><br>[kpc]</th><th>z<sub>d</sub><br>[kpc]</th></tr></thead><tbody><tr><td>thin</td><td>8.96</td><td>2.5</td><td>0.3</td></tr><tr><td>thick</td><td>1.83</td><td>3.02</td><td>0.9</td></tr></tbody></table>

<figure><img src="image_16.png" /></figure>

FIG. 15. Number density of stars encountered with an impact parameter $b$ divided by $b$ (to ease the comparison between the approximation and the full expression) vs the cosine of the angle between the subhalo trajectory and the normal to the galactic plane $\theta$. The lowest $\cos\theta$ is the closer the subhalo passes to the stellar disk. The approximate computation is given by the black solid line and is independent of $b$ while the full expression gives the colored solid lined for different values of $b$. There is a change in behaviour around $\cos\theta \sim b/R_d$ where $R_d \sim 3$ kpc is the approximate typical length scale of the stellar disk.

## Appendix D: Distribution of impact parameters

Let us define the galactic frame with origin the center of the galaxy $(G, \hat{e}_x, \hat{e}_y, \hat{e}_z)$ where the $z$-axis is perpendicular to the galactic plane (in the following, the galactic plane refers the middle of the stellar disk) and $\hat{e}_x$ and $\hat{e}_y$ are arbitrary orthonormal vectors. We assume that during the crossing the subhalo keeps a rectilinear trajectory with a velocity $v_c$ that intercepts the galactic plane at one point, at position $R$ at time $t = 0$ in the galactic frame. Consider now a star tagged by the letter $i$ that has a position $R_i$ at $t = 0$ in the same frame. We assume that on the interval of time the subhalo crosses the disk, this star keeps a linear trajectory with a velocity $v_i$ such that $v_i \cdot R_i \sim 0$. For simplicity we also introduce $S_i = R_i - R$. Then, at an arbitrary time $t$ the distance between the star and the subhalo is given by

$$d_i(t) = |S_i - v_{r,i}t| \qquad (D1)$$

where we introduce the relative velocity between the star and the subhalo $v_{r,i} = v_c - v_i$. The impact parameter for this specific star is defined as $b_i \equiv \min\{d_i(t)\}$, it yields

$$b_i = \frac{|S_i \times v_{r,i}|}{|v_{r,i}|} \qquad (D2)$$

Now we want to know what is the probability for this star to be at position $R_i$ and have a mass $m_i$. We call the associated probability distribution quantity $p_{(m_i, R_i)}(m_i, R_i) = p_{m_i}(m_i | R_i) p_{R_i}(R_i)$. In our model we make the usual approximation that $p_{m_i}(m_i | R_i) = p_{m_i}(m_i)$ so that it does not depend on the position. Moreover the probability distribution of positions is given by the mass density as

$$p_{R_i}(R_i) = \rho_{\star}(R_i) \left[ \int d^3R \rho_{\star}(R) \right]^{-1} = \frac{\rho_{\star}(R_i)}{m_{\star}^{\text{tot}}} \qquad (D3)$$
where $m_{\star}^{\text{tot}}$ is the total mass of stars and $\rho_{\star}(\mathbf{R}_i)$ their mass density at position $\mathbf{R}_i$. All stars being independent the joint PDF for their mass and position is

$$p_{(m_{\star}, \mathbf{R}_{\star})} (\{m_i\}_{i}, \{\mathbf{R}_i\}_{i}) = \prod_{i=1}^{N_{\star}} \left[ p_{m_{\star}}(m_i) p_{\mathbf{R}_{\star}}(\mathbf{R}_i) \right] \qquad (\text{D4})$$

with $N_{\star} = m_{\star}^{\text{tot}}/\bar{m}_{\star}$ the total number of stars. With all these ingredients it is possible to evaluate the number density of stars that are encountered with impact parameter $b$ and mass $M_{\star}$, knowing the trajectory of the subhalo. It is

$$\frac{d^2 \mathcal{N}}{db dm_{\star}} = \int \prod_{i=1}^{N_{\star}} dm_i d^3 \mathbf{R}_i p_{m_{\star}}(m_i) p_{\mathbf{R}_{\star}}(\mathbf{R}_i) \left[ \sum_{i=1}^{N_{\star}} \delta_{\mathrm{D}}(b_i - b) \delta_{\mathrm{D}}(m_i - m_{\star}) \right]. \qquad (\text{D5})$$

Massaging this expression, it is straightforward to prove that the mass probability distribution can be factored out. This leaves us with the number of density of stars with a given impact parameter that can be written under the compact form

$$\frac{d\mathcal{N}}{db} = \frac{1}{\bar{m}_{\star}} \int d^3\mathbf{R}_{\star} \rho_{\star}(\mathbf{R}_{\star}) \delta_{\mathrm{D}}(b_{\star} - b). \qquad (\text{D6})$$

Form the expression of the impact parameter it is convenient to make the change of variables in this integral $\mathbf{R}_{\star} \rightarrow \mathbf{S}_{\star} = \mathbf{R}_{\star} - \mathbf{R}$ and define the Dirac distribution on the squared value of the impact parameter. Eventually, in order to simplify the computation we will assume that all stars lay within an infinitely thin axi-symmetric disk of surface density $\Sigma_{\star}$ and the integration over the entire 3D space reduces to the integration on the Galactic plane. We therefore write

$$\frac{d\mathcal{N}}{db} = \frac{2b}{\bar{m}_{\star}} \int d^2\mathbf{S}_{\star} \Sigma_{\star} [R_{\star}(\mathbf{S}_{\star})] \delta_{\mathrm{D}}(b_{\star}^2 - b^2). \qquad (\text{D7})$$

Now in order to continue the computation we can parameterize the problem, without loss of generality, by choosing the convenient orientation of the basis $(\hat{e}_x, \hat{e}_y)$ such that $\mathbf{R} = (R, 0, 0)$. Moreover we define the relative velocity direction as $\mathbf{v}_{\mathrm{r},i}/|\mathbf{v}_{\mathrm{r},i}| = (\sin \theta \cos \varphi, \sin \theta \sin \varphi, \cos \theta)$ and $\mathbf{S}_{\star} = (S \cos \phi, S \sin \phi, 0)$ such that $d^2\mathbf{S}_{\star} = sdsd\phi$. Then, the expression of the impact parameter becomes

$$b_{\star}^{2} = s^{2} [1 - \sin^{2} \theta \cos^{2}(\phi - \varphi)]. \qquad (\text{D8})$$

The integration over the Dirac delta distribution in Eq. (D6) after the change of variable can then be done analytically by solving the delta for the angle $\phi$. It gives four distinct solutions in $[\varphi - \pi, \varphi + \pi]$ when $Sc_{\theta} < b < S$. These solutions can be written under the form

$$\phi_j = \varphi + \eta_j \arccos \left[ \frac{\chi_j}{\sin \theta} \sqrt{1 - \frac{b^2}{s^2}} \right], \qquad (\text{D9})$$

where we have introduced $\chi_j = (+1, +1, -1, -1)$ and $\eta_j = (+1, -1, +1, -1)$. The arccos function image being in the range $[0, \pi]$ only, the $\chi_j$ factor parameterizes the two solutions in the interval $[\varphi, \varphi + \pi]$ while the $\eta_j$ gives the two symmetric solutions in the interval $[\varphi - \pi, \varphi]$. Henceforth, in order to perform the integration over the Dirac distribution it is also necessary to provide the absolute value for the derivative of $b_{\star}^2$ with the variable $\phi$ evaluated in the four solution points. Using the properties of these solutions, this takes a simple form

$$\left. \frac{db_{\star}^2}{d\phi} \right|_{\phi=\phi_j} = 2s^2 \sin^2 \theta |\cos(\phi_j - \varphi)| |\sin(\phi_j - \varphi)| = 2\sqrt{s^2 - b^2} \sqrt{b^2 - s^2 \cos^2 \theta}. \qquad (\text{D10})$$

Then, we need to relate the value of $R_{\star}$ to $s$ in the four solution points. We denote $R_j$ these four quantities and introduce a new variable $y = s^2/b^2$. Using the simple relation $\mathbf{R}_{\star} = \mathbf{S}_{\star} + \mathbf{R}$ – that is the definition of $\mathbf{S}_{\star}$ – it is straightforward to show that $R_j^2 = R^2 + b^2y + 2h_j(y, c_{\theta}, \varphi)$ with the shorter notation $c_{\theta} = \cos \theta$ and the functions $h_j$ being

$$h_j(y, c_{\theta}, \varphi) \equiv \frac{\chi_j R b}{\sqrt{1 - c_{\theta}^2}} \left[ \cos \varphi \sqrt{y-1} - \frac{\eta_j}{\chi_j} \sin \varphi \sqrt{1 - y c_{\theta}^2} \right] \qquad (\text{D11})$$

In the end, the number density of stars with impact parameter $b$ can be written

$$\frac{d\mathcal{N}}{db} = \frac{b}{2} \int_{1}^{\frac{1}{c_{\theta}^2}} \frac{dy}{\sqrt{1 - y c_{\theta}^2 \sqrt{y-1}}} \sum_{j=1}^{4} \frac{\Sigma_{\star}(R_j)}{\bar{m}_{\star}} \qquad (\text{D12})$$
and the computation is complete. Eventually we also need relate the angles $\varphi$ and $\theta$ to physical quantities (i.e. the angle marking the direction of the stars and the subhalo). In practice this is difficult in a general case. We therefore consider that stars are motionless during all subhalo crossing so that $\mathbf{v}_{r,i} = \mathbf{v}_c$.

Formally the functions $h_j$ can be rewritten in a more convenient way as

$$h_j = \chi_j R b \sqrt{y} \cos \left( \varphi + \frac{\eta_j}{\chi_j} \arctan \left[ \sqrt{\frac{1 - y c_{\theta}^2}{y - 1}} \right] \right) \quad (D13)$$

Under this form it becomes easier to see that the functions $h_j$ and the associated radii $R_j$ are bounded

$$|h_j| \le R b \sqrt{y} \quad \text{and} \quad |R - b \sqrt{y}| < R_j < R + b \sqrt{y}. \quad (D14)$$

In particular, since $y < 1/c_{\theta}^2$, this previous inequality implies that $|R_j - R| < b/c_{\theta}$. If we consider that the variations of $\Sigma_{\star}$ on R are of typical length $R_d$ it yields that $\Sigma_{\star}(R_j) = \Sigma_{\star}(R)$ as long as $c_{\theta} \gg b/R_d$. Therefore we can simply rewrite the density of encountered stars as

$$\frac{d\mathcal{N}}{db} = 2b \frac{\Sigma_{\star}(R)}{\bar{m}_{\star}} \int_{1}^{\frac{1}{c_{\theta}^2}} \frac{dy}{\sqrt{1 - y c_{\theta}^2 \sqrt{y - 1}}} = \frac{\Sigma_{\star}(R)}{\bar{m}_{\star}} \frac{2\pi b db}{c_{\theta}}. \quad (D15)$$

## Appendix E: Energy and velocity distributions in stellar encounters

In this appendix we detail the PDF of several functions of interest and we prove several properties claimed in the main text. We start by introducing the PDF of final velocity in a given shell and give an analytical derivation. With the results we prove the relation between the median and the average kinetic energy kick of Eq. (40) on general grounds. In a third part, we illustrate this formal derivation with the example of the initial velocity following a Maxwell-Boltzmann distribution. Eventually we conclude by studying the impact of stellar encounters on the density profile using a simple criterion in order to justify the use of the median as the typical kinetic energy kick felts by all particles in a given shell during one crossing of the disk.

### 1. Probability distribution for the final velocity

We suppose that we know the PDF, $p_{\mathbf{v}}(\mathbf{v} | r)$, for the initial velocity $\mathbf{v}$ of particles at position $r$ in a subhalo and the probability for particles in that shell to receive a velocity kick $\Delta\mathbf{v}$ then the probability distribution of final velocity $v_f$ at position $r$ can be written under the form

$$p_{v_f}(v_f | r) = \int d^3\mathbf{v} p_{\mathbf{v}}(\mathbf{v} | r) \int d^d\Delta\mathbf{v} p_{\Delta\mathbf{v}}(\Delta\mathbf{v} | r) \delta_D[v_f - |\mathbf{v} + \Delta\mathbf{v}|]. \quad (E1)$$

We leave room here for the possibility of $\Delta\mathbf{v}$ being a 3D or 2D random vector with the dimension parameter $d$. Let us assume an isotropic initial velocity distribution so that $p_{\mathbf{v}}(\mathbf{v} | r) = p_{\mathbf{v}}(v | r)$ with $v = |\mathbf{v}|$. Then, it is possible to integrate first on the angular distribution of $\mathbf{v}$ in order to remove unnecessary degrees of freedom and get rid of the Dirac delta term. It yields

$$p_{v_f}(v_f | r) = \int d^3\mathbf{v} p_{\mathbf{v}}(v | r) p_{v_f}(v_f | v, r) \quad (E2)$$

with the definition of the PDF of $v_f$ knowing $v$ and $r$ being

$$p_{v_f}(v_f | v, r) = \frac{v_f}{v} \int d^d \Delta \mathbf{v} \frac{p_{\Delta \mathbf{v}}(\Delta \mathbf{v} | r)}{\Delta v} \Theta[\Delta v - |v + v_f|] \Theta[(v + v_f) - \Delta v] \quad (E3)$$

and where we introduced the velocity kick norm $\Delta v = |\Delta\mathbf{v}|$. For an isotropic distribution of $\Delta\mathbf{v}$ such that $p_{\Delta\mathbf{v}}(\Delta\mathbf{v} | r) = p_{\Delta\mathbf{v}}(\Delta v | r)$ it further simplifies to

$$p_{v_f}(v_f | v, r) = \frac{v_f}{v} \frac{\pi^{d/2}}{\Gamma(d/2)} \int_{|v - v_f|}^{|v + v_f|} d\Delta v (\Delta v)^{d-2} p_{\Delta\mathbf{v}}(\Delta v | r). \quad (E4)$$
<figure><img src="image_17.png" /><figcaption>FIG. 16. The inverse CCDF of kinetic energy knowing the velocity at three values $q = 0.25$, $q = 0.5$ and $q = 0.75$ with respect to the initial velocity $v$</figcaption></figure>

Let us now assume that $\Delta v$ follows a Gaussian distribution according to the result of the central limit theorem. Then the distribution on $\Delta v$ takes a simple form

$$p_{\Delta v}(\Delta v | r) = \frac{\pi^{-d/2}}{u^d} e^{-\frac{(\Delta v)^2}{u^2}} \quad \text{with} \quad u = u(r) = \sqrt{\frac{2\mathcal{N}(\delta v)^2}{d}} \qquad (E5)$$

which yields, with the change of variable $\Delta v \to u\sqrt{t}$ in the second line, and the introduction of the incomplete Gamma function in the third line,

$$p_{v_f}(v_f | v, r) = \frac{v_f}{v} \frac{1}{\Gamma(d/2)} \int_{|v-v_f|}^{|v+v_f|} d\Delta v \frac{(\Delta v)^{d-2}}{u^d} e^{-\frac{(\Delta v)^2}{u^2}} \qquad (E6)$$

$$= \frac{v_f}{uv} \frac{1}{2\Gamma(d/2)} \int_{\left(\frac{v-v_f}{u}\right)^2}^{\left(\frac{v+v_f}{u}\right)^2} dt \, t^{\frac{d-3}{2}} e^{-t} \qquad (E7)$$

$$= \frac{v_f}{uv} \frac{1}{2\Gamma(d/2)} \left[ \Gamma\left(\frac{d-1}{2}\right) \left(\frac{v-v_f}{u}\right)^2 - \Gamma\left(\frac{d-1}{2}\right) \left(\frac{v+v_f}{u}\right)^2 \right]. \qquad (E8)$$

In the end, here we have found a generic expression for the PDF of $v_f$ for any initial velocity distribution. Another interesting quantity is the associated CDF of $v_f$ that is defined as

$$F_{v_f}(\langle v_f | r \rangle) \equiv \int_0^{v_f} p_{v_f}(v'_f | r) dv'_f = \int d^3 v p_v(v | r) \left\{ F_{v_f}(\langle v_f | v, r \rangle) \equiv \int_0^{v_f} p_{v_f}(v'_f | v, r) dv'_f \right\}. \qquad (E9)$$

The CDF knowing $v$ and $r$ can therefore be computed with the expression of the PDF derived above. We will see the utility of this expression in the following paragraphs and especially when discussing the density profile modifications due to the encounter. However let us first inspect the properties of the median energy kick in light of this derivation.

## 2. Properties of the median energy kick

In the main text we have introduced the median value of the energy kick received by particles in a subhalo during the encounter with stars when it crosses the stellar disk. Moreover we have based our computation of the tidal radius
on the fact that the typical energy received in a shell is defined as the median value that can be approximated by the average value modulated by a coefficient between 0.5 and 0.7 (in the case of a 2D random walk in the velocity space). We prove this property here and show that it is independent of the initial velocity distribution. Similarly to the CDF for $v_f$ it is possible to introduce a related CDF for the kinetic energy kick $\Delta E$ as follows

$$F_{\Delta E}(\langle \Delta E | r \rangle) = \int d^3 v p_{\mathbf{v}}(v|r) F_{v_f}(\langle \sqrt{v^2 + 2\Delta E} | v, r \rangle). \qquad (E10)$$

Henceforth we can define a median value for $\Delta E$ knowing $v$ and $r$ as $\text{Med}(\Delta E | v, r)$. In Fig. 16 we show that value of $\widehat{\Delta E}(q, v)$ that is given by the implicit equation

$$F_{v_f} \left( \left\langle \sqrt{v^2 + 2\widehat{\Delta E}(q, v)} \right\rangle | v, r \right) = 1 - q, \qquad (E11)$$

for three different values of $q$. In particular, the median $\text{Med}(\Delta E | v, r) = \widehat{\Delta E}(0.5, v)$ is bounded by the asymptotes in $v = 0$ and $v \to \infty$. Series expansions in these two regime then show that the value of the boundaries are such that

$$\frac{d-1}{d} < \frac{\text{Med}(\Delta E | v, r)}{\widehat{\Delta E}} < x \quad \text{with } x \text{ solution of } \frac{\Gamma(d/2, xd/2)}{\Gamma(d/2)} = \frac{1}{2}. \qquad (E12)$$

In practice for $d = 2$, it yields $x = \ln(2)$ and for $d = 3$ it yields $x = 0.789$. Eventually, even though the total median value $\text{Med}(\Delta E)$ cannot be easily computed from $\text{Med}(\Delta E | v, r)$, the properties of the boundaries have to be conserved. Therefore we have shown that whatever the initial velocity distribution of velocity the median kinetic energy kick is always equivalent to the average up to an $\mathcal{O}(1)$ pre-factor.

### 3. The example of a Maxwellian initial velocity

In order to illustrate the theoretical development above and to connect with the main text we consider now that the PDF of the initial velocity is a Maxwell-Boltzmann distribution such that

$$p_{\mathbf{v}}(v|r) = \frac{1}{(2\pi\sigma_{\text{sub}}^2(r))^{3/2}} e^{-\frac{v^2}{2\sigma_{\text{sub}}^2(r)}}. \qquad (E13)$$

In order to compute the median of $\Delta E$ we can now compute the exact PDF and CDF. More precisely, similarly than for the final velocity in Eq. (E8), the PDF for $\Delta E$ is given by

$$p_{\Delta E}(\Delta E | r) = \int d^3 \mathbf{v} p_{\mathbf{v}}(v|r) \int d^d \Delta \mathbf{v} p_{\Delta \mathbf{v}}(\Delta v | r) \delta_D \left[ \Delta E - \frac{(\Delta v)^2}{2} + \mathbf{v} \cdot \Delta \mathbf{v} \right] \qquad (E14)$$

$$= \int d^d \Delta \mathbf{v} p_{\Delta \mathbf{v}}(\Delta v | r) \frac{2\pi}{\Delta v} \int \frac{|2\Delta E - (\Delta v)^2|}{2\Delta v} dv p_{\mathbf{v}}(v|r) \qquad (E15)$$

$$= \int d^d \Delta \mathbf{v} p_{\Delta \mathbf{v}}(\Delta v | r) \frac{1}{\sqrt{2\pi\sigma_{\text{sub}}^2(\Delta v)^2}} \exp \left( \frac{\left(\Delta E - \frac{(\Delta v)^2}{2}\right)^2}{2\sigma_{\text{sub}}^2 |\Delta \mathbf{v}|^2} \right) \qquad (E16)$$

Note that here we integrate over $v$ before integrating over $\Delta v$ for simplicity. According to the dimensionality of the random walk it yields

$$p_{\Delta E}(\Delta E | r) = \begin{cases} \frac{\exp\left(\frac{\Delta E}{2\sigma_{\text{sub}}^2} - \frac{|\Delta E|}{2\sigma_{\text{sub}}^2} \frac{\sqrt{1+s^2}}{s}\right)}{4\sigma_{\text{sub}}^2 s \sqrt{1+s^2}} & \text{if } d=2 \\ \frac{|\Delta E| e^{\frac{\Delta E}{2\sigma_{\text{sub}}^2}}}{4\pi\sigma_{\text{sub}}^4 s^2 \sqrt{1+s^2}} K_1\left[\frac{|\Delta E|}{2\sigma_{\text{sub}}^2} \frac{\sqrt{1+s^2}}{s}\right] & \text{if } d=3 \end{cases} \qquad (E17)$$

where $s^2 \equiv u^2/(8\sigma_{\text{sub}}^2) = \mathcal{N}\overline{\delta v^2}/(4d\sigma_{\text{sub}}^2) = \overline{\Delta E}/(2d\sigma_{\text{sub}}^2)$ and $K_1$ is the modified Bessel function of the second kind of order 1. Now we focus on the $d=2$ case and write down the CDF as

$$F(\langle \Delta E | r \rangle) = \begin{cases} 1 - \frac{1+\xi}{2\xi} e^{-\frac{\Delta E}{2\sigma^2}(\xi-1)} & \text{if } \Delta E \ge 0 \\ \frac{\xi-1}{2\xi} e^{\frac{\Delta E}{2\sigma^2}(1+\xi)} & \text{else.} \end{cases} \qquad (E18)$$
with $\xi = \sqrt{1+s^2}/s$. The complementary CDF, i.e. CCDF, introduced in the main text is denoted $\bar{F}(> \Delta E) = 1-F(< \Delta E)$. Eventually the energy $\Delta E(q)$ is defined implicitly through the CCDF as the solution of $\bar{F}(> \Delta E(q)) = q$. Therefore according to the value of $q$ the energy $\Delta E(q)$ can be written

$$ \frac{\Delta E(q)}{2\sigma^2} = \begin{cases} \frac{1}{1-\xi} \ln\left(\frac{2q\xi}{1+\xi}\right) & \text{if } q < \frac{1+\xi}{2\xi} \\ \frac{1}{1+\xi} \ln\left(\frac{2(1-q)\xi}{\xi-1}\right) & \text{else.} \end{cases} \qquad (E19) $$

As the average kinetic energy kick is such that $\bar{\Delta E}/(2\sigma_{\text{sub}}^2) = 2s^2 = 2(\xi^2 - 1)$, when $q \le 0.5$, it yields

$$ \frac{\text{Med}(\Delta E)}{\bar{\Delta E}} = \frac{\xi + 1}{2} \ln \left( \frac{1 + \xi}{2q\xi} \right) \ge \frac{\xi(1 - 2q) + 1}{2} \ge \frac{1}{2} \qquad (E20) $$

with $\text{Med}(\Delta E) = \Delta E(0.5)$ and where we used the inequality $\ln(1/x) \ge 1-x$. Eventually, when $q \ne 0.5$ the ratio diverges as $\xi \to \infty$ (corresponding to $s \to 0$). When $q = 0.5$ the ratio is a decreasing function of $\xi$ and

$$ \frac{\text{Med}(\Delta E)}{\bar{\Delta E}} \le \lim_{\xi \to 1} \frac{\text{Med}(\Delta E)}{\bar{\Delta E}} = \ln(2). \qquad (E21) $$

which proves the result of Eq. (E12).

### 4. Impact of stellar encounters on the phase space distribution function and on the density profile

In this last part we develop a method to evaluate the new phase space distribution function (PSDF) after one crossing of the disk assuming that the system is initially isotropic and that isotropy is conserved. If we know the initial PSDF $f(v, r)$, then trading $p_v(v|r)$ for $f$ in Eq. (E2) allows to define

$$ \tilde{f}(v_f, r) = \frac{1}{4\pi v_f^2} \int d^3\mathbf{v} f(v, r) p_{v_f}(v_f | v, r). \qquad (E22) $$

We can determine a new profile right after the disk crossing (before any relaxation effect that can drastically change the profile afterward) from the initial one by removing all particles with a final velocity higher than the escape velocity. This amounts to integrate $\tilde{f}$ on $v_f$ between 0 and the escape velocity

$$ v_{es}(r) = \sqrt{2|\Phi(r)|} \qquad (E23) $$

where $\Phi(r)$ is calculated with the initial profile before any possible reorganisation/relaxation<fnref n="3" />:

$$ \tilde{\rho}(r) = \int_{v_f < v_{es}(r)} d^3\mathbf{v}_f \tilde{f}(v_f, r) = \int d^3\mathbf{v} f(v, r) F_{v_f}(v_{es}(r) | v, r). \qquad (E25) $$

The next step is to evaluate the initial PSDF $f$. The common method when isotropy of the system is assumed is to use the Eddington formalism. However, this formalism shows some issues when the initial density profile is truncated at a finite radius. Finding a fully consistent PDF for a NFW profile with a sharp truncation is non trivial [105]. In practice one introduces the inverse gravitational potential $\Psi = |\Phi|$ and the inverse energy $\mathcal{E} = \Psi - v^2/2$. For an isotropic system the PSDF is only dependent on $\mathcal{E}$ and the Eddington formalism provides

$$ f(\mathcal{E}) = \frac{1}{\sqrt{8\pi^2}} \left\{ \frac{1}{\sqrt{\mathcal{E}}} \frac{d\rho}{d\Psi} \bigg|_{\Psi=0} + \int_0^\mathcal{E} \frac{d\Psi}{\sqrt{\mathcal{E} - \Psi}} \frac{d^2\rho}{d\Psi^2} \right\} \qquad (E26) $$

with the relation between $\Psi$ and $\rho$ given by Poisson's equation: $\Delta\Psi = -4\pi G\rho$. When $\Psi = 0$ at finite radius the derivative of $\rho$ with respect to $\Psi$ does not vanish and the first term $\propto 1/\sqrt{\mathcal{E}}$ is divergent and non physical. In order to keep things simple in the following we naively remove this term and compute the associated profile as

$$ \rho(r) = 4\pi\sqrt{2} \int_{0}^{\Psi(r)} \left\{ F(\mathcal{E}) \equiv \int_{0}^{\mathcal{E}} \frac{d\Psi'}{\sqrt{\mathcal{E} - \Psi'}} \frac{d^2\rho}{d\Psi'^2} \right\} \sqrt{\Psi(r) - \mathcal{E}} d\mathcal{E}. \qquad (E27) $$

<footnote n="3">³Note that this approach encompasses a similar method introduced in Ref. [51] and in which case it is possible to properly define a given velocity kick $\delta v$ and $v_f = v + \delta v$. One then recover the same expression by setting</footnote>

<footnote>$$ F_{v_f}(v_{es}(r) | v, r) = \Theta[v_{es}(r) - (v + \delta v)] \qquad (E24) $$</footnote>
Two major caveats can be pointed out here. First the reconstructed density is defined up to a constant that we consider to be 0 here. Moreover when modifying the PSDF it modifies the profile and one should also solve for a new gravitational potential in Poisson's equation. Therefore the initial potential is no longer an NFW profile but goes to 0 at the truncation radius while the impact of stars is evaluated with an exact NFW. Nevertheless, the PSDF in terms of velocity can be evaluated as

$$f(v, r) = F\left(\mathcal{E} = \left\{\Psi - \frac{v^2}{2}\right\}\right). \qquad (E28)$$

In the following, in order to parameterize our uncertainty due to the previously mentioned caveats, we introduce two other simple PSDF. First we consider the Maxwell-Boltzmann distribution Eq. (E13) used in the main text and we simply set

$$f(v, r) = \frac{\rho(r)}{(2\pi\sigma_{\text{sub}}^2(r))^{3/2}} e^{-\frac{v^2}{2\sigma_{\text{sub}}^2(r)}}. \qquad (E29)$$

This distribution is realistic but its main downside is that the velocity of particles can be higher than the escape velocity especially in the outskirts of the structure where the velocity dispersion given by the Jean's equation tends to the gravitational potential. Consequently, close to the truncation radius, the PSDF is not correctly normalised and under-evaluate the density. The second option we investigate is to set a cut-off in the velocity-space and renormalize

$$f(v, r) = \frac{\rho(r)}{K(r)} \left[ e^{-\frac{v^2}{2\sigma_{\text{sub}}^2(r)}} - e^{-\frac{v^2}{2\sigma_{\text{sub}}^2(r)}} \right] \qquad (E30)$$

where we set the normalisation factor $K(r)$ such that

$$\int d^3\mathbf{v} f(v, r) = \rho(r) \qquad (E31)$$

is ensured. Nevertheless, the main issue of this distribution is that the velocity dispersion is no longer $\sigma_{\text{sub}}$ (the real velocity dispersion being lower). Unfortunately trading $\sigma_{\text{sub}} \to \sigma \neq \sigma_{\text{sub}}$ and trying to recover the correct velocity dispersion does not work either in the outskirts of the halo: no good values of $\sigma$ can be found. Nonetheless, let us point out that the latter two configuration give a theoretical uncertainty in the outskirts.

In the left panel of Fig. 17 we represent the evolution of the PSDF with respect to the velocity at a fixed position in a typical subhalo after one disk crossing. We observe, as expected, that the stellar encounters naturally shift the distributions to higher values of velocity. In the right panel we show the corresponding new profile density in the top part and the comparison to the initial profile at the bottom for the different initial PSDF introduced above. We plotted the result for a subhalo crossing the disk at three different distances from the GC. We can remark that even in the most conservative case of the initial truncated Maxwell Boltzmann distribution of Eq. (E31) the density decreases toward the outskirts. At $R=8$ kpc and $R=4.5$ kpc the vertical dash-dotted lines represent the tidal radius obtains from the SL17 recipe using a typical kinetic energy kick equal to the median. In both cases they correspond to positions where the new densities are already below 50% of the initial one. At $R=1$ kpc the tidal radius is evaluated to be 0 as a sizeable part of the central particles are ejected (the blue curves being below the 50% threshold in the bottom panel). In order to determine the true final profile one also should integrate the effect of relaxation after the shock. Nevertheless, this simple analysis is enough to justify the SL17 recipe and the use of the median as a typical energy kick and to compare it to more realistic effects on the density profile.

## Appendix F: Probability distributions of the total energy kick (stars + disk shocking)

As mentioned in the main text we define a total energy kick as $\Delta E_{\text{tot}} = \Delta E + \Delta E_{\text{d}}$ where $\Delta E_{\text{d}}$ is the energy kick due to the disk shocking and $\Delta E$ is the energy kick due to the encounters with stars. In the following we assume that $\Delta E$ is distributed along as in Eq. (E17) in the $d=2$ case and $\Delta E_{\text{d}}$ is distributed along a Gaussian according to

$$p_{\Delta E_{\text{d}}}(\Delta E_{\text{d}}) = \frac{1}{\sqrt{2\pi\sigma_{\text{sub}}^2(\Delta \mathbf{v}_{\text{d}})^2 A_2(\eta)}} e^{-\frac{(\Delta E_{\text{d}} - A_1(\eta)\frac{(\Delta \mathbf{v}_{\text{d}})^2}{2})^2}{2\sigma_{\text{sub}}^2(\Delta \mathbf{v}_{\text{d}})^2 A_2(\eta)}} \qquad (F1)$$

where $\Delta \mathbf{v}_{\text{d}}$ is given in Eq. (4) and the adiabatic correction $A_1$ is introduced in Eq. (6). We also introduced a new adiabatic corrective factor: $A_2(\eta)$ for the dispersion (according to [72]). In the following we also use the parameter
<figure><img src="image_18.png" /><figcaption>FIG. 17. Left panel. Final PSDF (in red) in terms of the velocity at a given radius $r$ inside a subhalo of mass $m_{200} = 8 \times 10^{-6}$ M$_{\odot}$ with median concentration that crosses the disk at 8 kpc from the GC for three different initial PSDF (in blue). The value of $r$ is chosen as 75% of the tidal radius fixed by the virial radius $r_{200}$ here. Eddington is given by Eq. (E28), MB by Eq. (E29) and MB truncated by Eq. (E31). The normalisation was chosen such that the integral between $v = 0$ and $v = v_{es}$ of these curve should be 1 to recover the correct value of the density profile at position $r$. The only correctly normalised distribution is MB truncated. The MB simple PSDF is normalised on the range $v \in [0, \infty[$ and the Eddington PSDF is not normalised at all. The later two therefore under-predict the value of the density. Right top panel. The initial profile in cyan and the new profiles computed with the different initial PSDF for a NFW subhalo crossing the disk at three distances from the GC. Right bottom panel. Ratio of the new profile over the initial profile. The dash-dotted lines correspond to the tidal radius evaluated in the SL17 recipe by when choosing a typical kinetic energy kick of every shell at the median value.</figcaption></figure>

$s_d \equiv (\Delta\mathbf{v}_d)^2/(4\sigma_{\text{sub}}^2)$, similarly to $s = \overline{\Delta E}/(2d\sigma_{\text{sub}}^2)$ introduced above. Then, it is possible to give a PDF for $\Delta E_{\text{tot}}$. In order to do so we further introduce two variables

$$v_{\pm} \equiv \frac{1}{2} \left( \frac{\sqrt{1+s^2}}{s} \pm 1 \right) \sqrt{2A_2 s_d^2} \quad (\text{F2})$$

and we define a pseudo centred reduced variable corresponding to $\Delta E_{\text{tot}}$ of the form

$$\varepsilon \equiv \frac{1}{\sqrt{2A_2 s_d^2}} \left( \frac{\Delta E_{\text{tot}}}{2\sigma_{\text{sub}}^2} - s_d^2 A_1 \right) \quad (\text{F3})$$

in order to simplify the expressions. The main goal of this analysis is to quantify the asymmetry of the PDF around the average value in order to evaluate whether considering an average value for $\varepsilon$ is relevant and if not what should be the better choice. Therefore using an affine shift to define this new variable we do not lose in generality in that sense. With these two definitions we can evaluate a PDF for $\varepsilon$ under the form

$$p_{\varepsilon}(\varepsilon | \nu_{-}, \nu_{+}) = \frac{\nu_{+}\nu_{-}}{\nu_{+} + \nu_{-}} e^{-2\varepsilon\nu_{-} + \nu_{+}^2} \text{erfc}(\nu_{-} - \varepsilon) + \frac{\nu_{+}\nu_{-}}{\nu_{+} + \nu_{-}} e^{2\varepsilon\nu_{+} + \nu_{-}^2} \text{erfc}(\nu_{+} + \varepsilon). \quad (\text{F4})$$

The average value of $\varepsilon$ can be rewritten $\bar{\varepsilon} = (\nu_{+} - \nu_{-})/(2\nu_{+}\nu_{-})$. Moreover from this PDF we can also derive a CCDF for $\varepsilon$. It yields the following expression

$$\bar{F}_{\varepsilon}(\varepsilon | \nu_{-}, \nu_{+}) = -\frac{\nu_{+}\nu_{-}}{\nu_{+} + \nu_{-}} \frac{e^{2\varepsilon\nu_{+} + \nu_{+}^2}}{2\nu_{+}} \text{erfc}(\nu_{+} + \varepsilon) + \frac{\nu_{+}\nu_{-}}{\nu_{+} + \nu_{-}} \frac{e^{-2\varepsilon\nu_{-} + \nu_{-}^2}}{2\nu_{-}} \text{erfc}(\nu_{-} - \varepsilon) + \frac{1}{2} \text{erfc}(\varepsilon). \quad (\text{F5})$$

The main advantage of this parameterization is to be only depends on two parameters $\nu_+$ and $\nu_-$ and therefore to be easy to compute numerically. Moreover it can be used to show that the CCDF evaluated at $\varepsilon = \bar{\varepsilon}$ satisfies

$$\frac{1}{e} \sim 0.368 < \bar{F}_{\varepsilon}(\bar{\varepsilon} | \nu_{+}, \nu_{-}) < \frac{1}{2}. \quad (\text{F6})$$
<figure><img src="image_19.png" /><figcaption>FIG. 18. $\Delta E_{\text{tot}}(q = 0.5)/\overline{\Delta E_{\text{tot}}}$ vs the parameters $s$ and $s_d$ under the assumption that $A_2 = A_1$ for simplicity. The black curves represent contours of constant $\overline{\Delta E_{\text{tot}}}: \overline{\Delta E_{\text{tot}}} = 0.2\sigma^2$ (dash-dotted), $\overline{\Delta E_{\text{tot}}} = 2\sigma^2$ (solid), $\overline{\Delta E_{\text{tot}}} = 20\sigma^2$ (dashed).</figcaption></figure>

Considering the average value is a good way to evaluate that the maximal energy gain for a fraction of at least 37% of the particle. This is understandable as we had, $\text{Med}(\Delta E_{\text{tot}}) \le \overline{\Delta E}$ in the case without disk shocking and because the addition of $\Delta E_d$, which has a symmetric distribution, do not further *un-symmetrise* the PDF for $\Delta E_{\text{tot}}$. Therefore this justify entirely that we can evaluate the total energy as roughly being the sum of the averages. In Fig. 18 we represent $\text{Med}(\Delta E)/\overline{\Delta E_{\text{tot}}}$ vs $s$ and $s_d$ under the assumption $A_1(\eta) = A_2(\eta)$ for simplicity. When $s_d \gg s$ the distribution is symmetrised again with respect to the average and the ratio is close to 1. When $s \gg s_d$ the symmetry is maximally broken and that yields a ratio of $\sim 0.7$ as discussed in the scenario where only stellar encounters impact on the subhalos.