# Capture of field stars by dark substructures

Jorge Peñarrubia1,2, Raphaël Errani3,4, Matthew G. Walker3, Mark Gieles5,6,
Tjarda C. N. Boekholt7,8
1Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
2Centre for Statistics, University of Edinburgh, School of Mathematics,
Edinburgh EH9 3FD, UK
3McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University, 5000 Forbes Ave., Pittsburgh, PA 15213, United States
4Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR
7550, F-67000 Strasbourg, France
5ICREA, Pg. Lluís Companys 23, E08010 Barcelona, Spain
6Institut de Ciències del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB),
Martí i Franquès 1, E08028 Barcelona, Spain
7Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory,
University of Oxford, Parks Road, Oxford, OX1 3PU, UK
8NASA Ames Research Center, Moffett Field, CA 94035, USA Email:
<EMAIL_ADDRESS>
###### Abstract

We use analytical and $N$-body methods to study the capture of field stars by
gravitating substructures moving across a galactic environment. We find that
the majority of captured stars move on temporarily-bound orbits that escape
from the substructure potential after a few orbital revolutions. In numerical
experiments where a substructure model is immersed in a galaxy on a circular
orbit, we also find particles that remain bound to the substructure potential
for indefinitely-long times. This population is absent from substructure
models initially placed outside the galaxy on an eccentric orbit. We show that
gravitational capture is most efficient in dwarf spheroidal galaxies (dSphs)
on account of their low velocity dispersions and high stellar phase-space
densities. In these galaxies ‘dark’ sub-subhaloes which do not experience in-
situ star formation may capture field stars and become visible as stellar
overdensities with unusual properties: (i) they would have a large size for
their luminosity, (ii) contain stellar populations indistinguishable from the
host galaxy, and (iii) exhibit dark matter (DM)-dominated mass-to-light
ratios. We discuss the nature of several ‘anomalous’ stellar systems reported
as star clusters in the Fornax and Eridanus II dSphs which exhibit some of
these characteristics. A large population of DM sub-subhaloes with a mass
function ${\rm d}N/{\rm d}M_{\bullet}\sim M_{\bullet}^{-\alpha}$ generates
stellar systems with a luminosity function, ${\rm d}N/{\rm d}M_{\star}\sim
M_{\star}^{-\beta}$, where $\beta=(2\alpha+1)/3=1.6$ for $\alpha=1.9$.
Detecting and characterizing these objects in dSphs would provide
unprecedented constraints on the particle mass and cross section of a large
range of DM particle candidates.

###### keywords:

Galaxy: kinematics and dynamics; galaxies: evolution; Cosmology: dark matter.

##  1 Introduction

One of the strongest predictions from Cold Dark Matter (CDM) cosmology is the
existence of self-gravitating haloes devoid of visible matter (i.e., ‘dark’).
Such objects arise because star formation becomes inefficient in haloes with
virial masses below $\sim 10^{7}$–$10^{9}\,\mathrm{M}_{\odot}$ (White & Rees
1978; Bullock et al. 2000; Bovill & Ricotti 2009; Benitez-Llambay & Frenk
2020, Pereira-Wilson et al. 2023), while the minimum subhalo mass associated
with the free-streaming length of ‘cold’ particles with masses above $\sim
1\,\mathrm{GeV/c}^{2}$ is at least $\sim 13$ orders of magnitude smaller,
$<10^{-6}\,\rm M_{\odot}$ (e.g. Schmid et al. 1999; Hofmann et al. 2001; Green
et al. 2005; Loeb & Zaldarriaga 2005; Diemand et al. 2005). Below this scale,
fluctuations of the power spectrum are heavily suppressed, de facto imposing a
truncation at the low end of the halo mass function (e.g. Benson 2017 and
references therein). Hence, detecting this truncation would provide a direct
constraint on the DM particle mass.

The lack of visible matter together with their tiny masses make the detection
of dark subhaloes extremely difficult. Current observational efforts range
from searching for gamma- ray annihilation signals (e.g. Ackermann et al.
2014; Bringmann et al. 2014) to modelling substructure in strongly-lensed
galaxies (Koopmans 2005; Vegetti & Koopmans 2009; Li et al. 2013; Vegetti et
al. 2014). In the Milky Way, encounters with individual subhaloes can induce
significant perturbations in cold tidal streams (Ibata et al. 2002; Johnston
et al. 2002; Yoon et al. 2011; Carlberg 2013; Ngan et al. 2016, Erkal et al.
2016). For example, Bonaca et al. (2019) find that some observed features in
the GD-1 stream, including a gap and an off-stream spur of stars, are best
reproduced by the past encounter with a dark subhalo with a mass
$M_{\bullet}\sim 10^{6}$–$10^{8}\,{\rm M}_{\odot}$ and a scale radius
$c_{\bullet}\lesssim 10\,{\rm pc}$. Puzzlingly, these constraints imply a
matter density comparable to the stellar density in globular clusters,
$\rho_{\bullet}=M_{\bullet}/(2\pi c_{\bullet}^{3})\gtrsim 10^{2}\,{\rm
M}_{\odot}\,{\rm pc}^{-3}$, which is several orders of magnitude denser than
CDM subhaloes with similar masses found in cosmological simulations of
structure formation (e.g. Moliné et al. 2017; Diemer & Joyce 2019).

Inferring the existence of dark subhaloes using dynamical probes is
complicated by the unknown number of baryonic compact objects – such as
stellar black holes, neutron stars, white dwarfs, free-floating planets, giant
molecular clouds etc. – lurking in the Galaxy. In general, it is not
straightforward to isolate perturbations arising from different populations of
gravitating objects (Peñarrubia 2018). In this regard, dwarf spheroidal
galaxies (dSphs) provide relatively clean targets, as their gravitational
potentials appear to be fully dominated by dark matter (e.g. Mateo 1998).

This paper challenges the common assumption that DM subhaloes that do not form
stars in situ remain ‘dark’. Here, we show that dark subhaloes can
capture111In this paper, the word ‘capture’ is broadly used to describe any
dynamical process wherein a field particles undergoes a transition from galaxy
orbit to an orbit around a substructure. baryonic matter as they orbit around
the host galaxy, becoming ‘visible’ as localized substructures of co-moving
bodies with high mass-to-light ratios and extended sizes. Dark subhaloes that
are massive enough to capture field stars would bear resemblance to stellar
clusters, but with atypical properties, e.g. their stellar populations would
be chemically and chronologically identical to the local galactic field, and
they would be DM dominated. As we will show below, these properties are akin
to those of the ‘anomalous’ clusters detected in some Milky Way dSphs,
suggesting the intriguing possibility that these systems may be in fact
agglomerates of field stars captured by dark substructures orbiting in the
host galaxies. This possibility is explored below in some detail.

Capture processes in the classical222Two point masses on circular orbits, plus
a mass-less tracer restricted 3-body problem have been studied for a long
time. For example, the pioneering work of Szebehely (1967) showed that a
finite number of solutions exists where the lightest particle is transferred
from one distinct mode of motion around the most massive point-mass to another
distinct mode around the intermediate-mass one. Hunter (1967), Heppenheimer
(1975) and Heppenheimer & Porco (1976) pointed out that Jupiter’s outer
satellites could have been captured in this way. More recently, Suetsugu &
Ohtsuki (2013) study temporary capture of planetesimals by a giant planet,
while Jílková et al. (2015) consider the scenario where the inner Oort Cloud
was captured from another star during a close encounter in their birth
cluster. Recently, 3-body captures in accretion discs have also gained
attention as a possible source of black-hole binaries. E.g. Li et a. (2022),
Boekholt et al. (2023) and Rowan et al. (2023) show that close encounters
between two black holes orbiting around a supermassive black hole can form
bound pairs.

The dynamics of 3-body capture events is extremely complex. Using a
combination of numerical and analytical methods, Petit & Hénon (1986) showed
that (i) captures in a 3-body system are temporary events that ultimately
induce the dissolution of the bound pair, and (ii) capture only happens for
extremely fined-tuned combinations of impact parameters and relative
velocities which exhibit a self-similar, Cantor-like structure. More recent
work of Boekholt et al. (2023) confirmed these results and found that the
phase space structure that leads to capture resembles a Cantor set with a
fractal dimension of $\simeq 0.4$.

Peñarrubia (2023, hereafter Paper I) studies 3-body captures in a galactic
environment, where the intermediate and lighter bodies (both point masses)
follow orbits in a massive (extended) potential. Using numerical experiments,
Paper I shows that encounters between a massive object and field particles can
be locally described as slingshot manoeuvres, in which the lighter body can
increase/decrease its speed or redirect its path. For a capture event to
happen, the galactic tidal field must decelerate an approaching lighter body
to a degree where it temporarily orbits the intermediate one. Crucially, a
point-mass moving through a sea of lighter particles generates a localized
overdensity – or ‘halo’ – of tidally-trapped particles333A massive perturber
also deflects stars into an overdense ‘wake’ that trails behind it (e.g.
Kalnajs 1971; Weinberg 1986). In contrast, this work focuses on field stars
that become temporarily bound to the perturber. , which reaches a steady state
as the rate of bodies captured from the field becomes comparable to those
being lost to galactic tides.

The reverse process by which stars escape from the intermediate potential in a
restricted 3-body system is not fully understood. For example, Fukushige &
Heggie (2000) show that stars with energies above the energy of escape can
remain inside the tidal radius of the intermediate body for very long times,
and some do not escape at all. These particles are typically dubbed ‘potential
escapers’ (see e.g. Küpper et al. 2010; Daniel et al. 2017a).

In this work we expand the analysis of Paper I to intermediate bodies with an
extended (i.e., not point-like) mass distribution. As an application, we study
capture of field stars by DM subhaloes in a wide range of galactic
environments. The paper is arranged as follows. §2 extends the statistical
method outlined in Paper I to model gravitational captures by point-mass
objects to substructures an extended size. §3 presents numerical experiments
that test the accuracy of the theoretical equations. §4 discusses these
results in the context of dark matter particle physics and outlines future
follow-up work. Finally, in §5, we summarize our main results.

##  2 Statistical theory

This Section summarizes the main techniques applied in this paper to construct
a statistical theory that describes the spatial and kinematical distribution
of tracer particles temporarily bound to self-gravitating substructures moving
in a host galaxy.

Section 2.1 follows the steps outlined in Paper I to compute the average
number of field particles that have negative binding energies
$E=v^{2}/2+\Phi_{\bullet}(r)<0$ within a spherical volume $V=4\pi r^{3}/3$
($N_{b})$, where $r$ and $v$ are measured relative to the substructure, as
well as their number density ($n_{\star}$) and velocity dispersion
($\sigma_{\star}$) profiles in dynamical equilibrium within a generic
potential $\Phi_{\bullet}$. Here, we will focus on capture of field stars,
although this analysis can be extended to other tracers moving in the host
galaxy potential, including gas particles444To model the dynamics of gaseous
particles, a pressure term must be included in the equations of motion
presented in §3.1.

To gain further physical insight, we will inspect the case of dark
substructures sourcing a Hernquist (1990) potential

$\displaystyle\Phi_{\bullet}(r)=-\frac{GM_{\bullet}}{r+c_{\bullet}},$ (1)

which recovers the results derived in Paper I for point-masses if the profile
scale-length is set to $c_{\bullet}=0$. Note that the density profile
associated with the potential (1) has a centrally-divergent cusp, $\rho\sim
r^{-1}$ at $r\ll c_{\bullet}$, as the universal profile found in CDM
simulations of structure formation (e.g. Navarro, Frenk & White 1997).

Section 2.2 applies our results to Milky Way dwarf spheroidal galaxies, which,
due to their intrinsic properties and relatively close distances, appear to be
the most interesting objects to test the existence of dark substructures via
captured field stars.

Section 2.3 analyzes the cosmologically-motivated case of a large population
of extended substructures with a power-law mass function.

###  2.1 Distribution of trapped particles

####  2.1.1 Weak perturbations

The number of field stars temporarily bound to a moving substructure sourcing
a potential $\Phi_{\bullet}$ can be estimated statistically from the (local)
distribution function under the following simplifying assumptions: (i) field
stars move on uncorrelated (random) trajectories within a small volume element
$V$, (ii) their number density is roughly constant within the volume $V$, such
that $n({\textbf{{R}}}_{\bullet}+{\textbf{{r}}})\approx
n({\textbf{{R}}}_{\bullet})=n=N/V$, where ${\textbf{{R}}}_{\bullet}$ is the 3D
position vector of the substructure within the host galaxy. This is known as
the local approximation, and holds insofar as the (local) density profile
rolls slowly, i.e. $r\ll d=|\nabla n/n|^{-1}$. (iii) The relative velocity
distribution of particles within the volume element $V=4\pi r^{3}/3$ follows a
Maxwellian distribution displaced by the reflex velocity of the point-mass
$p(\mbox{\boldmath$v$})=(2\pi\sigma^{2})^{-3/2}\exp[-(\mbox{\boldmath$v$}+\mbox{\boldmath$V$}_{\bullet})^{2}/(2\sigma^{2})$
where $\sigma=\sigma(\mbox{\boldmath$R$}_{\bullet})$ is the local, one-
dimensional velocity dispersion of field stars. Under the Maxwellian
approximation, the mean-squared (relative) velocity between the background
particles and the substructure is $\langle
v^{2}\rangle=3\sigma^{2}+V_{\bullet}^{2}$, where ${\textbf{{V}}}_{\bullet}$ is
the velocity vector of the substructure with respect to the host galaxy
centre.

The approximation that field stars move on uncorrelated trajectories is
accurate insofar as the number of stars bound to $\Phi_{\bullet}$ represents a
small fraction of the total number of stars within the volume element,
$N=n\,V$. Paper I shows that this condition requires an escape velocity
$v_{e}\lesssim\sigma$, where
$v_{e}(\mbox{\boldmath$r$})=\sqrt{2|\Phi_{\bullet}(\mbox{\boldmath$r$})|}$ is
calculated at a distance
$r=|\mbox{\boldmath$R$}-\mbox{\boldmath$R$}_{\bullet}|$ from the centre of the
potential $\Phi_{\bullet}$.

Under these conditions, the average number of stars in the volume $V$ with
negative specific energy $E=v^{2}/2+\Phi_{\bullet}<0$ can be derived from the
local phase space density of field stars as

$\displaystyle N_{b}(r)$ $\displaystyle=\int_{V}{\rm
d}^{3}r\,n(\mbox{\boldmath$r$})\int_{E<0}{\rm d}^{3}v\,p(\mbox{\boldmath$v$})$
(2)
$\displaystyle=\frac{1}{3}\sqrt{\frac{2}{\pi}}ze^{-V_{\bullet}^{2}/(2\sigma^{2})}\frac{n}{\sigma^{3}}\int_{V}{\rm
d}^{3}r\,v_{e}^{3}(\mbox{\boldmath$r$})+\mathcal{O}(v_{e}/\sigma)^{5}$
$\displaystyle=\frac{32\sqrt{\pi}}{9}(GM_{\bullet})^{3/2}\,e^{-V_{\bullet}^{2}/(2\sigma^{2})}\frac{n}{\sigma^{3}}$
$\displaystyle~{}~{}~{}\times\frac{r^{2}-4c_{\bullet}r+8c_{\bullet}^{3/2}(r+c_{\bullet})^{1/2}-8c_{\bullet}^{2}}{(r+c_{\bullet})^{1/2}},$

for $r\ll d$ and $v_{e}/\sigma\lesssim 1$. To gain physical intuition, the
last expression adopts a Hernquist potential (1) with a escape speed
$v_{e}(r)=\sqrt{2|\Phi_{\bullet}|}=\sqrt{2GM_{\bullet}/(r+c_{\bullet})}$. It
is trivial to show Equation (2) recovers Equation (11) of Paper I for a point-
mass with $c_{\bullet}=0$. Note that the number of stars with negative
energies within the volume $V$ is proportional to the local mean phase-space
density of the field, $Q\equiv n/\sigma^{3}$, and that $N_{b}$ drops for
substructures that are not at rest with the background ($V_{\bullet}>0$).

It is worth emphasizing that captured stars move on chaotic orbits that do not
conserve energy or angular momentum, as shown in §3.2. However, from a
statistical point of view, we can count how many stars are bound to the
substructure potential (i.e. $E<0$) at any given time and compare that number
against the value estimated in Equation (2). Hénon & Petit (1986) found that
chaotic orbits only remain trapped within the potential $\Phi_{\bullet}$ over
a finite amount of time before being lost to galactic tides. As shown in Paper
I, this leads to a population of bound field stars that reaches a steady state
as the capture rate (defined as the net number of field particles with $E$
flipping from positive to negative values) equals the loss rate (i.e the net
number of field particles with $E$ flipping from negative to positive values).
Paper I finds that steady state is typically reached over Smoluchowski’s
(1916) ‘fluctuation mean life’

$\displaystyle T(r)=\sqrt{\frac{2\pi}{3}}\frac{r}{\langle
v^{2}\rangle^{1/2}},$ (3)

which roughly corresponds to the time that a particle moving on a straight-
line trajectory takes to cross the volume element, i.e. the so-called
‘crossing’ time. Typically, this time-scale is much shorter than the orbital
time of the substructure around the host galaxy. Paper I shows that in steady-
state the average number of field stars with $E<0$ can be estimated from
Equation (2) as $N_{\rm ss}=\alpha\,N_{b}$, where $\alpha$ is the so-called
abundance parameter, which is set by the ‘dynamical survival time’ 555Defined
as the time over which a trapped particle continuously has $E<0$ before being
lost to galactic tides.  that trapped objects remain bound to the potential
$\Phi_{\bullet}$. Numerical experiments shown in Paper I show that $\alpha$ is
close to unity. The estimates shown below set the abundance parameter
$\alpha=1$.

The density enhancement due to the presence of trapped interstellar stars in
steady state can be directly derived from (2) as

$\displaystyle\delta_{\star}(r)$
$\displaystyle\equiv\frac{1}{4\pi\,n\,r^{2}}\frac{{\rm d}N_{\rm ss}}{{\rm
d}r}$ (4)
$\displaystyle=\alpha\frac{4}{3\sqrt{\pi}}e^{-V_{\bullet}^{2}/(2\sigma^{2})}\frac{1}{\sigma^{3}}|\Phi_{\bullet}|^{3/2}$
$\displaystyle=\alpha\frac{4}{3\sqrt{\pi}}\frac{(GM_{\bullet})^{3/2}}{\sigma^{3}}e^{-V_{\bullet}^{2}/(2\sigma^{2})}\frac{1}{(r+c_{\bullet})^{3/2}},$

which converges to a centrally-divergent ‘density spike’ $\delta_{\star}\sim
r^{-3/2}$ (Gondolo & Silk 1999) in the point-mass limit $c_{\bullet}\to 0$.
For extended halos ($c_{\bullet}>0$) the density profile is flat for
$r\lesssim c_{\bullet}$. Note that the density enhancement (4) scales as
$\delta_{\star}\sim|\Phi_{\bullet}|^{3/2}$, which allows a straightforward
analysis of a wide range of substructure potentials.

The velocity dispersion of stars trapped in a potential $\Phi_{\bullet}$ in
steady-state can be found from (4) by solving the isotropic Jeans equations

$\displaystyle\sigma_{\star}^{2}(r)$
$\displaystyle=\frac{1}{\delta_{\star}(r)}\int_{r}^{\infty}{\rm
d}r^{\prime}\delta_{\star}(r^{\prime})\bigg{|}\frac{{\rm
d}\Phi_{\bullet}}{{\rm d}r}\bigg{|}_{r^{\prime}}$ (5)
$\displaystyle=\frac{2}{5}|\Phi_{\bullet}|$
$\displaystyle=\frac{2}{5}\frac{GM_{\bullet}}{r+c_{\bullet}}.$

which shows that $\sigma_{\star}/v_{e}=(1/5)^{1/2}\approx 0.45$ at all radii.
It is important to emphasize that Equation (5) is independent of the speed of
the dark object across the host galaxy ($V_{\bullet}$).

Of particular relevance is the distance at which the density of stars equals
that of the galactic background. It is convenient to examine the point-mass
limit first. Solving $\delta(r_{\epsilon})=1$ in Equation (4) with
$c_{\bullet}\to 0$ and $\alpha=1$ yields

$\displaystyle
r_{\epsilon}=\bigg{(}\frac{16}{9\pi}\bigg{)}^{1/3}\,e^{-V_{\bullet}^{2}/(3\sigma^{2})}\frac{GM_{\bullet}}{\sigma^{2}},$
(6)

which is dubbed the ‘thermal’ critical radius in Paper I owing to the
Maxwellian velocity dependence $\exp[-V_{\bullet}^{2}/(3\sigma^{2})]$
multiplying the critical radius $r_{0}=2GM_{\bullet}/\sigma^{2}$ in (6). Thus,
the density of field stars bound to an object with mass $M_{\bullet}$ will
exceed that of the galactic background on scales below the thermal critical
radius, i.e. $\delta_{\star}>1$ at $r<r_{\epsilon}$.

For objects with an extended mass distribution,
$\delta_{\star}(r_{\epsilon}^{\rm ext})=1$ happens at $r_{\epsilon}^{\rm
ext}=r_{\epsilon}-c_{\bullet}=\kappa\,r_{\epsilon}$, where $\kappa\equiv
1-c_{\bullet}/r_{\epsilon}$ is the compactness of the substructure. It is
straightforward to show that the condition $\delta_{\star}>1$ requires
$\kappa>0$ (or $c_{\bullet}<r_{\epsilon}$), which means that only
substructures that are sufficiently compact can contain an over-dense region
of captured field stars.

We inspect this issue in next Section with the aid of $N$-body tools.

In what follows, we assume that a dark substructure becomes ‘visible’ –i.e. it
can be detected as a localized stellar over-density– if the volume
$V_{\epsilon}=4\pi r_{\epsilon}^{3}/3>0$ contains at least one bound star,
i.e. $N_{b}(r_{\epsilon})>1$. This condition can only be satisfied by dark
objects above a certain mass threshold, $M_{\bullet}>M_{\rm min}$. Using (4)
with $c_{\bullet}=0$ and (6) and setting $\alpha=1$ we find

$\displaystyle M_{\rm
min}=\frac{3\pi^{1/3}}{4}\,e^{V_{\bullet}^{2}/(3\sigma^{2})}\frac{D\,\sigma^{2}}{G},$
(7)

here $D=(2\pi\,n)^{-1/3}$ is a measure of the average separation of stars in
the volume element $V$ (Peñarrubia 2018). Notice again that the formation of
stellar over-densities is strongly suppressed in dark objects moving at high
speed ($V_{\bullet}\gg\sigma$) with respect to the galactic background.

####  2.1.2 Strong perturbations

The statistical theory presented above assumes that the substructure potential
$\Phi_{\bullet}$ induces weak perturbations on the trajectories of field stars
as they move across the volume $V=4\pi r^{3}/3$. One can easily show that this
approximation breaks down on scales $r\lesssim r_{\epsilon}$, where the
potential well $|\Phi_{\bullet}|$ becomes deeper than the specific kinetic
energy of field particles, $K=3\sigma^{2}/2$.

To derive the distribution of captured stars at small distances from the
substructure we use two empirical results found in Paper I (see also §3): (i)
the phase-space density of captured stars becomes approximately constant in
the vicinity of the substructure, such that $f_{\star}({\bf r},{\bf v})\simeq
f_{0}$ on scales $r\lesssim r_{\epsilon}$, and (ii) the population of captured
stars reaches a steady state on time scales comparable to Smoluchowski’s
(1916) ‘fluctuation mean life’ defined by Equation (3).

The density profile of an equilibrium ensemble of stars moving in the
potential $\Phi_{\bullet}$ can be derived from the local distribution function
without a priori knowledge of their orbital trajectories. That means, there is
no need to specify whether orbits are regular (permanently bound) or chaotic
(temporarily bound), as long as the distribution of these objects is
homogeneous in phase-space. Hence,

$\displaystyle n_{\star}[\Phi_{\bullet}(r)]$ $\displaystyle=\int_{E<0}{\rm
d}^{3}vf_{\star}({\bf r},{\bf v})$ (8) $\displaystyle=4\pi
f_{0}\,\int_{0}^{|\Phi_{\bullet}|}{\rm d}E[2(E-\Phi_{\bullet})]^{1/2}$
$\displaystyle=\frac{8\sqrt{2}\pi}{3}f_{0}|\Phi_{\bullet}|^{3/2},$

whereas the isotropic 1D velocity dispersion is

$\displaystyle\sigma^{2}_{\star}[\Phi_{\bullet}(r)]$
$\displaystyle=\frac{1}{3\,n_{\star}(r)}\int_{E<0}{\rm
d}^{3}v\,v^{2}f_{\star}({\bf r},{\bf v})$ (9)
$\displaystyle=\frac{4\pi\,f_{0}}{3\,n_{\star}(r)}\,\int_{0}^{|\Phi_{\bullet}|}{\rm
d}E[2(E-\Phi_{\bullet})]^{3/2}$ $\displaystyle=\frac{2}{5}|\Phi_{\bullet}|.$

Notice that the velocity dispersion profile (9) matches Equation (5), whereas
the density profile (8) recovers Equation (4) for a local Maxwellian
distribution function

$\displaystyle
f_{0}=\alpha\frac{n}{(2\pi\sigma^{2})^{3/2}}\,e^{-V_{\bullet}^{2}/(2\sigma^{2})}.$
(10)

Hence, the profiles derived at distances $r\gtrsim r_{\epsilon}$ under the
condition that $\Phi_{\bullet}$ induces weak perturbations also extend to
regions where $\Phi_{\bullet}$ dominates the local gravitational field,
$r\lesssim r_{\epsilon}$. We have numerically checked that captured particles
have a constant phase-space density in dark substructures with an extended
mass profile (see Section 3), which agrees with the numerical results shown in
Paper I in which the intermediate object is a point-mass.

By construction, combination of Equations (8), (9) and (10) returns a mean
phase-space density of bound stars that remains constant across the volume
element

$\displaystyle
Q_{\star}\equiv\frac{n_{\star}(r)}{\sigma^{3}_{\star}(r)}=\alpha\frac{5}{3}\sqrt{\frac{10}{\pi}}e^{-V^{2}_{\bullet}/(2\sigma^{2})}\,Q.$
(11)

Equation (11) highlights two interesting results. First, $Q_{\star}$ solely
depends on the potential $\Phi_{\bullet}$ through the value of the abundance
parameter, $\alpha$, and second the phase-space density of captured stars
drops exponentially for a potential $\Phi_{\bullet}$ that is not at rest with
the background ($V_{\bullet}>0$). Numerical tests in §3.1 (see also Appendix
3) discuss the accuracy of these theoretical profiles.

Figure 1: Upper panel: Stellar phase-space density $Q=n/\sigma^{3}=1/(2\pi
D^{3}\sigma^{3})$ of MW dSphs as a function of their half-light radii.
References for the data shown are given in footnote 2.2. Notice that the
phase-space density of the MW stellar halo at $R=20\,{\rm kpc}$ is $\gtrsim 6$
orders of magnitude lower than in dSphs. This difference grows at larger
Galactocentric radii. Middle panel: Minimum substructure mass ($M_{\rm min}$)
derived from Equation (7). Substructures with $M_{\bullet}>M_{\rm min}$ are
expected to contain stellar over-densities ($\delta_{\star}>1$) with at least
one bound field star ($N_{b}>1$), see text. The metal-rich (MR) and metal-poor
(MP) components of the Fornax dSph are highlighted with red and blue filled
dots, respectively. Notice that the value of $M_{\rm min}$ in the MW halo is
relatively large on account of its high velocity dispersion. Lower panel:
Minimum thermal critical radius $r_{\epsilon}^{\rm min}$ of substructures with
a mass $M_{\bullet}=M_{\rm min}$. Stellar over-densities made of capture stars
have a physical size comparable to the radii of globular clusters.

###  2.2 Application: Milky Way dwarf spheroidals

The above analytical expressions can be used to inspect the efficiency of
gravitational captures in different galactic environments. According to our
statistical theory, capture of field stars is most likely in regions where the
stellar phase-space density $Q=n/\sigma^{3}=(2\pi D^{3}\,\sigma^{3})^{-1}$ is
high, that is the intra-stellar separation $D$ is small and the velocity
dispersion is low (i.e. ‘cold’). Upper panel of Fig. 1 shows that the galaxies
with the highest phase-space densities are the ultra-faint dwarf spheroidals,
which also are the smallest, coldest and most metal-poor galaxies in the known
Universe (e.g. Simon 2019). References for the data shown are given in
footnote 2.255footnotetext: Half-light radii and velocity dispersion for the
dSphs shown in Fig. 1 are from McConnachie (2012) (version January 2021), with
updates for Antlia 2, Crater 2 (Ji et al. 2021) and Tucana (Taibi et al.
2020). .

In what follows, we assume that the stellar component of satellite galaxies
follow a Plummer profile, $n^{\rm dSph}(r)=n^{\rm
dSph}_{0}[1+(r/a)^{2}]^{-5/2}$, with a central density $n^{\rm
dSph}_{0}=N^{\rm dSph}_{\star}/(4\pi a^{3}/3)$, and a (3d) half-light radius
$r^{\rm dSph}_{h}=1.305\,a$. The average inter-stellar separation is estimated
as $D=(2\pi n_{0}^{\rm dSph})^{-1/3}$. The number of stars in a dSph galaxy is
calculated from the total stellar mass as $N^{\rm dSph}_{\star}=M^{\rm
dSph}_{\star}/\langle m_{\star}\rangle$, where, $\langle m_{\star}\rangle$ is
the average mass of a single star. Given that dSphs have approximately flat
velocity dispersion profiles, we set $\sigma=\sigma^{\rm dSph}$ (Walker et
al., 2007).

Dwarf galaxies with multiple chemo-dynamical components are particularly
interesting in this context. In these systems, the probability to capture
field stars strongly depends on stellar metallicity. For example, the Fornax
dSph contains two prominent chemo-dynamical populations (Battaglia et al.
2006; Walker & Peñarrubia 2011; Amorisco & Evans 2012) with metallicities
$\mathrm{[Fe/H]}=-1.8$ and $-0.65$ and ages $>10\,{\rm Gyr}$, and $\sim
0.2$–$2\,{\rm Gyr}$, respectively (Rusakov et al. 2021). The metal-rich (MR)
component contains more stars $N_{\star}^{\rm MR}\sim 3\times 10^{7}$, is more
centrally concentrated, $r_{h}^{\rm MR}\approx 530\,{\rm pc}$ and has a colder
velocity dispersion, $\sigma^{\rm MR}\approx 10\,{\rm km\,s^{-1}}$, than the
metal-poor (MP) population, which has $N_{\star}^{\rm MP}\sim 2\times 10^{7}$,
$r_{h}^{\rm MP}\approx 1070\,{\rm pc}$ and $\sigma^{\rm MP}\approx 14.4\,{\rm
km\,s^{-1}}$ (Walker & Peñarrubia 2011). Thus, MR stars have much higher
phase-space densities than MP stars,

$\displaystyle\frac{Q^{\rm MR}}{Q^{\rm MP}}=\frac{n^{\rm MR}}{n^{\rm
MP}}\bigg{(}\frac{\sigma^{\rm MP}}{\sigma^{\rm
MR}}\bigg{)}^{3}=\frac{N_{\star}^{\rm MR}}{N_{\star}^{\rm
MP}}\bigg{(}\frac{r_{h}^{\rm MP}}{r_{h}^{\rm MR}}\frac{\sigma^{\rm
MP}}{\sigma^{\rm MR}}\bigg{)}^{3}\sim 33.$ (12)

Middle panel of Fig. 1 shows estimates of $M_{\rm min}$ derived from Equation
(7) for the MW satellites in our sample. To simplify the analysis,
substructures are placed at rest within the host galaxies ($V_{\bullet}=0$).
For reference, we also plot the value of $M_{\rm min}$ associated with the
stellar halo of the Milky Way at a Galactocentric distance of $R=20\,{\rm
kpc}$ using the stellar density and velocity dispersion measured by Deason et
al. (2011). This panel reveals a number of interesting points. Notice first
that the mass threshold ($M_{\rm min}$) is lowest in the smallest dSphs, which
simply reflects the tight correlation between half-light radius and phase-
space density shown in the upper panel. If the goal is to detect ‘dark’
substructures that contain a stellar population of captured field stars, then
ultra-faint dSphs with half-light radii $r_{h}^{\rm dSph}\sim 30\,{\rm pc}$
are the best targets, as they may be sensitive to the presence of dark
substructures with masses above $M_{\bullet}\gtrsim 10^{4}\,{\rm M}_{\odot}$.
In contrast, the mass threshold increases by an order of magnitude,
$M_{\bullet}\gtrsim 10^{5}\,{\rm M}_{\odot}$, in the ‘classical’ dwarf
galaxies with $r_{h}^{\rm dSph}\gtrsim 300\,{\rm pc}$. In the case of the
Fornax dSph, the MR stellar component a mass threshold $M_{\rm min}$ that is
approximately an order of magnitude lower than for the MP population.

Stars trapped in dark substructures can lead to the formation of localized
stellar overdensities ($\delta_{\star}>1$) in the host galaxy. Bottom panel of
Fig. 1 shows that the size of the overdense regions lies in the range
$r_{\epsilon}\gtrsim 1$–$10\,{\rm pc}$, thus being comparable to the size of
stellar clusters. This leads to the intriguing possibility that some of the
‘anomalous’ stellar systems detected in MW dSphs may actually correspond to
stars temporarily trapped in dark substructures. We discuss this scenario in
§4.1.

Overdensities of trapped stars are expected to have high mass-to-light ratios.
This can be shown by measuring the substructure mass within the over-density
volume as
$M=M_{\bullet}(<r_{\epsilon})=M_{\bullet}r_{\epsilon}^{2}/(r_{\epsilon}+c_{\bullet})^{2}$,
and comparing it to the bound stellar mass $L=N_{b}(r_{\epsilon})\,\langle
m_{\star}\rangle$. For the point-mass case $c_{\bullet}=0$, applying (2) and
(6) returns a dimension-less mass-to-light ratio

$\displaystyle\frac{M}{L}\equiv\frac{M(<r_{\epsilon})}{N_{b}(r_{\epsilon})\,\langle
m_{\star}\rangle}=\frac{M_{\rm min}^{3}}{M^{2}_{\bullet}\,\langle
m_{\star}\rangle},$ (13)

where $M_{\rm min}$ is given by equation (7). According to the estimates
plotted in Fig. 1, field stars bound to susbtructures with a mass
$M_{\bullet}\sim M_{\rm min}$ will appear as extremely DM-dominated objects,
$M/L=M_{\rm min}/\langle m_{\star}\rangle\sim 10^{4}$–$10^{5}$ for an average
stellar mass $\langle m_{\star}\rangle=1\,{\rm M}_{\odot}$. In substructures
with larger masses, $M_{\bullet}>M_{\rm min}$, the mass-to-light ratio of
captured stars drops as $M/L\sim M_{\bullet}^{-2}$. Note that Equation (13)
should not be applied to substructures where $M/L\lesssim 1$, as the theory
outlined in §2.1 ignores the contribution of trapped stars to the underlying
potential.

In the MW halo, capture of field stars is inefficient on account of their high
velocity dispersion. For illustration, here we adopt $\sigma^{\rm MW}\sim
124\,{\rm km\,s^{-1}}$ and $n^{\rm MW}=10^{-4}\,{\rm pc}^{-3}$ at $R=20\,{\rm
kpc}$ from the MW centre (e.g. Deason et al. 2011), noting in passing that
capture of field stars in this particular stellar halo model becomes
systematically less efficient at larger distances. Note that the velocity
dispersion of the MW stellar halo is approximately one order of magnitude
higher than in the classical dSphs. As a result, the minimum substructure mass
needed to capture stars from the field increases up to $M_{\rm min}\gtrsim
10^{8}\,{\rm M}_{\odot}$, which is $\gtrsim 3$ orders of magnitude larger than
in the MW dwarf spheroidals. Crucially, this mass threshold is lower than the
virial mass of the classical dSphs, which appear to be embedded in DM haloes
with virial masses $10^{9}$–$10^{10}\,{\rm M}_{\odot}$ (e.g. Peñarrubia et al.
2008; Errani et al. 2018), suggesting that dSphs may contain a population of
captured MW halo stars. However, field haloes in this mass range have scale
radii that are a factor $\gtrsim 100$ larger than the thermal critical radius
plotted in the bottom panel of Fig. 1, $c_{\bullet}\gtrsim 10^{3}\,{\rm pc}$
(e.g. Joyce & Diemand 2019), which implies a negative compactness parameter,
$\kappa<0$. Therefore, field stars trapped in dSphs will appear as a diffuse
envelope of co-moving stars with subdominant densities with respect to the
local MW background ($\delta_{\star}\ll 1$). Therefore, the detection of this
stellar population will be challenging, and may require simultaneous modelling
of the kinematic and chemical composition of a large sample of stars at the
locations of dSphs. This problem will be explored in a separate contribution.

Figure 2: Left panel: mass $M_{\bullet}$ and scale radius $c_{\bullet}$ of DM
sub-subhaloes colour-coded according to the number of captured dwarf galaxy
stars, $N_{\star}=N_{b}(r_{\epsilon})$, Equation (2), where $r_{\epsilon}$ is
the thermal critical radius, Equation (6). In this model, the host galaxy has
properties similar to the Fornax dSph: it follows a Plummer profile with a
total luminosity of $N^{\rm dSph}=5\times 10^{7}$, a half-light radius
$r_{h}^{\rm dSph}=0.89\,{\rm kpc}$ and it has an average velocity dispersion
$\sigma^{\rm dSph}=12.2\,{\rm km\,s^{-1}}$. White solid and dotted lines
respectively show the virial mass-size relation of haloes and subhaloes taken
from the Aquarius simulations (Springel et al. 2008, fig. 26) extrapolated
down to the mass scales of interest (see text). Middle panel Compactness
$\kappa=1-c_{\bullet}/r_{\epsilon}$ of the left-panel models. Only models with
$\kappa>0$ lead to over-densities of captured stars with $\delta_{\star}>1$.
Right panel: Mass-to-light rations of the left and middle-panel models derived
from Equation (13). Stellar systems trapped in subhaloes with masses
$M_{\bullet}\lesssim 10^{7}\,{\rm M}_{\odot}$ will exhibit dark-matter
dominated kinematics.

###  2.3 Sub-subhalo populations in a CDM framework

According to CDM models of structure formation, all galaxies are expected to
contain a large population of DM subhaloes with a power-law mass function
${\rm d}N/{\rm d}M_{\bullet}\sim M_{\bullet}^{-\alpha}$ and a slope
$\alpha\approx 1.9$ (e.g. Giocoli et al. 2008; Springel et al. 2008). In the
classical fluid limit, the number of subhaloes diverges at low masses. In
contrast, if DM is made of cold massive particles a truncation at a minimum
subhalo mass $M_{\rm min}$ is expected to arise on scales comparable to the
particle free-streaming length. For WIMPS with masses above $\sim
1\mathrm{GeV/c}^{2}$ the truncation lies below planet-mass scales, $M_{\rm
min}\lesssim 10^{-6}\,\rm M_{\odot}$ (e.g. Schmid et al. 1999; Hofmann et al.
2001).

Within a CDM context, dark subhaloes hosting a population of captured field
stars are expected to follow a luminosity function dictated by the underlying
halo mass function. To illustrate this, let us define $N_{\star}=N_{b}(r_{\rm
out})$ within a volume size larger than the substructure scale radius, $r_{\rm
out}\gg c_{\bullet}$. From Equation (2)

$\displaystyle
N_{\star}=\frac{32\sqrt{\pi}}{9}(GM_{\bullet})^{3/2}\,e^{-V_{\bullet}^{2}/(2\sigma^{2})}\frac{n}{\sigma^{3}}\,r_{\rm
out}^{3/2},~{}~{}~{}~{}{\rm for}~{}~{}~{}r_{\rm out}\gg c_{\bullet}.$ (14)

To combine Equation (14) with the subhalo mass function we apply the chain
rule $({\rm d}N/{\rm d}M_{\bullet}){\rm d}M_{\bullet}\sim
N_{\star}^{-(2\alpha+1)/3}{\rm d}N_{\star}$, which leads to a luminosity
function

$\displaystyle\frac{{\rm d}N}{{\rm
d}M_{\star}}=A_{0}\bigg{(}\frac{M_{\star}}{\langle
m_{\star}\rangle}\bigg{)}^{-(2\alpha+1)/3},$ (15)

where $M_{\star}=N_{\star}\langle m_{\star}\rangle$ is the mass of a stellar
clump, and $A_{0}$ is an arbitrary normalization. Thus, the luminosity
function of stellar substructures also follows a power-law ${\rm d}N/{\rm
d}M_{\star}\sim M_{\star}^{-\beta}$, albeit with a shallower index than the
underlying subhalo mass function, $\beta=(2\alpha+1)/3=1.6$ for $\alpha=1.9$.

Whether or not the these objects can be detected as localized stellar over-
densities mainly depends on the compactness of DM haloes that capture them
(see §2). To illustrate this point, Fig. 2 shows the mean properties of field
stars trapped in dark substructures covering a wide range of masses and scale
radii. For this plot, we adopt the same dwarf galaxy model as outlined in
§3.1. Namely, field stars follow a Plummer profile with a total luminosity of
$N^{\rm dSph}=5\times 10^{7}$, a half-light radius of $r_{h}^{\rm
dSph}=0.89\,{\rm kpc}$ and a mean velocity dispersion $\sigma^{\rm
dSph}=12.2\,{\rm km\,s^{-1}}$, similar to the overall properties of the Fornax
dSph (see §2.2). For simplicity, we assume that dark substructures are at rest
with the field ($V_{\bullet}=0$).

Left panel of Fig. 2 shows models colour-coded according to the number of
bound stars enclosed within their thermal critical radius,
$N_{b}(r_{\epsilon})$ calculated from Equations (2) and (6). As shown in Fig.
1, point-mass objects with $M_{\bullet}\gtrsim 10^{5}\,{\rm M}_{\odot}$
contain $N_{\star}\gtrsim 1$ stars trapped within their thermal critical
radius. In extended substructures, the mass threshold for capture increases
mildly.

It is interesting to extrapolate the virial mass-scale radius relation of CDM
field haloes (white dotted line) and subhaloes (white solid line) found in the
Aquarius simulations (Springel et al. 2008) down to the relevant mass
scales666More precisely, we adopt the power-law fits of the $r_{\rm
max}-V_{\rm max}$ measurements plotted in fig. 26 of Springel et al. (2008).
As proxy for total halo mass, we use the halo virial mass $M_{200}$ derived
from $r_{\rm max}-V_{\rm max}$ under the assumption that field haloes as well
as subhaloes follow NFW profiles. For Hernquist spheres the circular velocity
peaks at the scale radius, i.e. $r_{\rm max}=c_{\bullet}$.. The first
noteworthy result is that subhaloes are more compact than field haloes. This
is because substructures that lose a significant amount of mass to tides tend
to be more concentrated than field halos of the same mass (e.g. Errani &
Navarro 2021). In particular, extrapolation of the Aquarius subhaloes down to
$M_{\bullet}\gtrsim 10^{5}\,{\rm M}_{\odot}$ leads to scale-radii
$c_{\bullet}\gtrsim 25\,{\rm pc}$, whereas field haloes of the same mass are a
factor $\sim 2$ larger, $c_{\bullet}\gtrsim 50\,{\rm pc}$ (Ludlow et al. 2016;
see also Diemer & Joyce 2019).

Only substructures with positive compactness ($\kappa\equiv
1-c_{\bullet}/r_{\epsilon}$) can become visible as localized stellar over-
densities ($\delta_{\star}>1$). Middle panel shows that the condition
$\kappa>0$ (or $c_{\bullet}<r_{\epsilon}$) translates into a linear relation
$c_{\bullet}\lesssim 0.83\,GM_{\bullet}/\sigma^{2}$, which is steeper than the
shallow mass-scale radius relation of Aquarius subhaloes, $c_{\bullet}\sim
M_{\bullet}^{0.46}$. As a result, Aquarius subhaloes with masses
$M_{\bullet}\lesssim 5\times 10^{6}\,{\rm M}_{\odot}$ would be ‘fluffy’
($\kappa<0$). In these objects, trapped stars would have sub-dominant
densities with respect to the field ($\delta_{\star}<1$), complicating their
detection in photometric surveys. However, in spectroscopic surveys these
objects may appear as localized regions with distinct kinematics, i.e. a low
velocity dispersion $\sigma_{\star}<\sigma$, and/or a significant velocity
offset $\Delta v$.

It is important to stress that the mass-size relation of field haloes and
subhaloes contain a significant amount of scatter that is not plotted in Fig.
2 (e.g. the concentration of field haloes has a standard deviation of $\sim
0.15\,$dex at a fixed mass, Ludlow et al. 2016). Crucially, a low-probability,
high-density tail in the mass-size distribution would lead to a population of
low-luminosity over-densities with the size of stellar clusters. This is of
particular relevance as the existence of such a high-density, low-probability
tail of substructures is a strong prediction from the cold dark matter
paradigm (Peñarrubia et al. 2010; Errani & Peñarrubia 2020; Errani & Navarro
2021). Unfortunately, given the poor theoretical understanding of the mass
function and the density profile of sub-subhaloes in dSphs, at present it is
not possible to make predictions on the number of stellar over-densities of
captured field stars. We will discuss this issue in §4.1.

Right panel of Fig. 2 shows the mass-to-light ratio (13) for the models shown
in previous panels. As expected, we find that field stars trapped in low-mass
($M_{\bullet}\lesssim 10^{7}\,{\rm M}_{\odot}$) substructures exhibit DM-
dominated kinematics ($M/L\gg 1$). As a note of caution, recall that the
analytical equations derived in §2.1 neglect the self-gravity of captured
stars, which is not a valid assumption in systems with mass-to-light ratios
$M/L\lesssim 1$.

In summary, Fig. 2 suggests that field stars captured by dark sub-subhaloes
may resemble stellar clusters with anomalous properties: (i) they would have a
large size for their luminosity, (ii) contain stellar populations with ages
and metallicities indistinguishable from the host galaxy, and (iii) exhibit
DM-dominated mass-to-light ratios.

##  3 Statistical Experiments

This Section carries a suite of numerical experiments that help us to test the
accuracy of the theoretical equations derived in §2.1, and shed light on the
dynamics of stars trapped by a substructure orbiting in a Fornax-like dSph
galaxy. The Section goes as follows

§3.1 provides a brief overview of the properties of the galactic field
particles and the numerical set-up used to solve the equations of motion. We
refer interested readers to Appendix A of Paper I for a detailed description
of the integration tools.

§3.2 injects a static substructure potential $\Phi_{\bullet}$ in a pre-
existing field of particles in dynamical equilibrium within $\Phi^{\rm dSph}$.
This set-up is similar to the analytical conditions adopted in §2.1, and we
can therefore anticipate a close match between theoretical predictions and the
numerical results. Interestingly, we will see that the dynamics of field
particles trapped in $\Phi_{\bullet}$ can be broadly separated in two
families: (i) permanent captures, which become bound immediately after the
substructure is placed in the field at $t=0$, and remain bound for arbitrarily
long times, and (ii) temporary captures, which as the name indicates only
remain bound for a finite amount of time.

§3.3 explores a physically-motivated case where a substructure mass grows
while moving along a circular orbit, thus sourcing a time-dependent potential
$\Phi_{\bullet}(t)$. We will see that these models also capture particles
trapped on permanent orbits around $\Phi_{\bullet}$, and that one time-scales
$t\gtrsim|\dot{\Phi}_{\bullet}/\Phi_{\bullet}|^{-1}$ the phase-space
distribution of particles with $E<0$ become indistinguishable from the static
case explored in §3.2.

In §3.4 we analyze capture of field stars by a substructure accreted onto a
dwarf galaxy on an eccentric orbit. The potential $\Phi_{\bullet}$ is
initially placed at orbital apocentre in a region populated by no field
particles (therefore it contains no permanent captures by construction), with
a small orbital pericentre that reaches the inner-most regions of the dSph.

§3.5 considers substructure models with a fixed mass and different scale
radii. The numerical results stress that only substructures that are
sufficiently compact generate overdensities of captured field stars.

###  3.1 Initial conditions & set-up

We generate realizations of $N=5\times 10^{6}$ stellar tracer particles in
dynamical equilibrium within a Dehnen (1993) potential $\Phi^{\rm dSph}$ with
a total mass of $M^{\rm dSph}=3\times 10^{9}\,{\rm M}_{\odot}$ and a scale
radius $c^{\rm dSph}=2\,{\rm kpc}$. We run experiments with cuspy ($\gamma=1$)
and cored ($\gamma=0$) profiles. Given that the stellar luminosity of the
Fornax dSph is $\sim 5\times 10^{7}\,{\rm L}_{\odot}$ (e.g. Rusakov et al.
2021), the particle luminosity in our models is $10\,{\rm L}_{\odot}$.

Field particles in this potential follow an $\alpha$-$\beta$-$\gamma$ profile
(Zhao 1996)

$\displaystyle
n(R)=\frac{n_{0}}{(R/R_{0})^{\gamma_{f}}[1+(R/R_{0})^{\alpha_{f}}]^{(\beta_{f}-\gamma_{f})/\alpha_{f}}},$
(16)

with $n_{0}$ chosen such that $4\pi\int_{0}^{\infty}{\rm d}r\,r^{2}\,n(r)=N$.
We consider two tracer models: (i) a spherical Plummer-like (1911)
profile777We notice in order to generate a positive distribution function via
an Eddington inversion, the inner slope $\gamma_{f}$ can be small but not
exactly zero. with a scale radius $R_{0}=690\,{\rm pc}$ and slopes
$(\alpha_{f},\beta_{f},\gamma_{f})=(2,5,0.1)$, and (ii) a truncated model with
$(\alpha_{f},\beta_{f},\gamma_{f})=(2,30,0.2)$ and a scale radius
$R_{0}=460\,{\rm pc}$. To guarantee dynamical equilibrium, orbital velocities
are assigned using Eddington (1916) inversion (see Errani & Peñarrubia 2020).
For the Plummer profile, this returns a luminosity-averaged 1D velocity
dispersion $\langle v^{2}\rangle^{1/2}/3=12.2\,{\rm km\,s^{-1}}$. This is our
reference model, which approximately matches the overall phase-space density,
size and luminosity of the Fornax dSph plotted in Fig. 1. The luminosity-
averaged velocity dispersion of the truncated model is $\langle
v^{2}\rangle^{1/2}/3=11.7\,{\rm km\,s^{-1}}$.

The motion of individual tracer particles are solutions to two sets of
differential equations. The first set describes the orbit of the substructure
in the dSph potential

$\displaystyle\ddot{\mbox{\boldmath$R$}_{\bullet}}=-\nabla\Phi^{\rm
dSph}({\mbox{\boldmath$R$}_{\bullet}}).$ (17)

The second set describes the trajectories of tracer particles in the
galactocentric frame

$\displaystyle\ddot{\mbox{\boldmath$R$}}=-\nabla\Phi^{\rm
dSph}(\mbox{\boldmath$R$})-\nabla\Phi_{\bullet}(\mbox{\boldmath$R$}-\mbox{\boldmath$R$}_{\bullet})+{\mbox{\boldmath$F$}}_{\rm
coll},$ (18)

where $\mbox{\boldmath$R$}-\mbox{\boldmath$R$}_{\bullet}$ is the relative
distance of the particle to the substructure, and ${\mbox{\boldmath$F$}}_{\rm
coll}=\sum_{i=1}^{N_{\rm clump}}{\mbox{\boldmath$f$}}_{i}$ is the net force
induced by$N_{\rm clump}$ clumps orbiting in the substructure potential
$\Phi_{\bullet}$ (for example, these could be planets orbiting the Sun in
Paper I, or sub-subhaloes in the DM halo of dSphs in the current work). For
simplicity, here we assume that the substructure potential is ‘smooth’ by
setting $F_{\rm coll}=0$. We will analyze the effect of random encounters with
clumps on the dynamics of captured particles in follow-up work.

In Sections 3.2 and 3.3, substructures are placed on circular orbits around a
cored DM halo ($\gamma=0$) at a galactocentric radius $R_{\bullet}=0.5\,{\rm
kpc}$ with a tangential velocity $V_{\bullet}=V_{c}(R_{\bullet})=14.35\,{\rm
km\,s^{-1}}$. The local dynamical time is therefore
$\Omega^{-1}=R_{\bullet}/V_{\bullet}=500\,{\rm pc}/14.35\,{\rm
km\,s^{-1}}\simeq 34\,{\rm Myr}$. Field stars follow a Plummer (1911) profile,
have a local velocity dispersion of $\sigma(R_{\bullet})=11.2\,{\rm
km\,s^{-1}}$, and are separated by an average distance
$D(R_{\bullet})=4.8\,{\rm pc}$.

For illustration, we choose a substructure mass $M_{\bullet}=10^{6}\,{\rm
M}_{\odot}$. To estimate the thermal critical radius we adopt $\sigma=\langle
v^{2}\rangle^{1/2}/3=12.2\,{\rm km\,s^{-1}}$ in Equation (6), which returns
$r_{\epsilon}=18.9\,{\rm pc}$. Note that this size is comparable to that of
the extended stellar clusters discussed in the Introduction. We consider a
‘compact’ substructure with $\kappa=+0.8$, which translates into a scale
radius $c_{\bullet}=(1-\kappa)\,r_{\epsilon}=3.78\,{\rm pc}$ (notice that this
value is smaller than the typical scale radius of field haloes with a vitial
mass $M_{\rm vir}(z=0)\sim 10^{6}M_{\odot}$ in the Aquarius simulation, see
Fig. 2). The time-scale needed to cross the over-density size by a random
field particle is much shorter than the orbital time,
$T_{\epsilon}=r_{\epsilon}/\sigma=18.9\,{\rm pc}/(12.2\,{\rm
km\,s^{-1}})\simeq 1.5\,{\rm Myr}$. As a result, in these models the
population of captured particles reaches steady state quickly.

We choose a volume size around $\Phi_{\bullet}$ which matches the local tidal
radius derived below from Equation (21), $r_{\rm out}=r_{t}\simeq 190\,{\rm
pc}$. The volume $V=4\pi r^{3}_{\rm out}/3$ contains approximately
$N=n\,V=42530$ field stars. From Equation (2), we expect that approximately
$N_{b}(r_{\rm out})=2000$ stars will have negative energies.

Figure 3:  Distribution of capture times ($t_{c}=t-t_{\rm bound}$) defined as
the time span since a particle becomes bound ($t_{\rm bound}$), i.e. the
energy $E$ flips from positive to negative, to the current simulation time
($t$). Time is measured in units Smoluchowski’s time scale ($T$). Here, we
only include particles within a volume size equal to the local tidal radius,
$r_{\rm out}=r_{t}$ (see text). Notice the presence of a substantial number of
particles that become bound at $t\approx 0$, which result in a a strong peak
in the distribution at $t_{c}\approx t$ shifting towards the right as the
integration time increases.

###  3.2 Static substructures potential on circular orbits

In this Section we place a substructure potential $\Phi_{\bullet}$ at a
galactocentric radius $R_{\bullet}$ with a circular velocity
$V_{\bullet}=V_{c}(R_{\bullet})$ in a sea of tracer particles with an extended
profile (see §3.1) in dynamical equilibrium. By chance, at the initial
snapshot ($t=0$) a number of tracer particles are found within a volume
element $V=4\pi r^{3}_{\rm out}/3$ centred at $\Phi_{\bullet}$ with a relative
velocity below the escape velocity, $v<v_{e}$, which implies a negative
energy, $E<0$. In what follows, we refer to these particles as ‘immediate
captures’ in order to distinguish them from ‘chaotic 3-body captures’ (e.g.
Petit & Henon 1986) taking place during the dynamical integration of their
orbits as a result of the interplay between the tracer particle and the
substructure and host galaxy potentials.

Next, we follow the motions of the substructure and the field particles by
solving Equations (17) and (18). At each time-step, we identify particles
located within the volume $V$ whose energy flips from positive to negative,
label them as ‘captured’, and mark their ID’s and the time when this happens
($t_{\rm bound})$.

As the simulation proceeds, we notice that orbits with $E<0$ can be broadly
separated into two families. The first group is made of ‘permanent’ captures,
which correspond to particles that remain bound for arbitrarily long times and
were already bound at $t=0$ (i.e. these are immediate captures that survive
until the final snapshot of the simulation). Here, it is worth stressing that
not all immediate captures are permanently bound, as some of them become
tidally stripped by the galactic potential. We will come back to this point
below. The second group corresponds to ‘temporary’ captures. As the name
indicates, these particles were unbound ($E>0$) at the initial snapshot $t=0$.
At some point during the integration their energy sign flips from positive to
negative, and back to positive, as these particles only remain bound for a
finite amount of time.

The presence of particles on permanently-bound orbits can be easily identified
in the distribution of capture times plotted in Fig. 3 at three different
snapshots of the simulation. Independently of the time of the measurement, we
find a very prominent spike centred at $t_{c}=t-t_{\rm bound}=t$. This spike
is populated by particles that were already bound at $t=0$ ($t_{\rm
bound}=0$). Away from the spike, we observe a broad distribution of capture
times, which is populated by particles that were initially unbound and become
captured at a later time of the simulation, such that $t_{\rm bound}>0$. Their
distribution of capture times peaks at $t_{c}\sim 2\,T$, although in this
experiment we also find a significant number of temporary captures that remain
bound for as long as $t_{c}\sim 50\,T$. In our models, all particles captured
with $t_{\rm bound}>0$ are temporary events that ultimately lead to the tracer
particle being released back to the galactic potential.

Figure 4: Examples of permanent and temporary captures. Left panels: Motion of
captured particles in a coordinate frame centred at the potential
$\Phi_{\bullet}$ of a ‘compact’ Hernquist substructure with
$M_{\bullet}=10^{6}\,{\rm M}_{\odot}$ and a compactness $\kappa=0.2$ (see §2.1
for details), which leads to a scale radius
$c_{\bullet}=(1-\kappa)\,r_{\epsilon}=3.78\,{\rm pc}$. Notice how the
temporary capture eventually leaves the volume of observation as it becomes
gravitationally unbound. Upper and lower right panels: Time-evolution of the
specific energy $E=v^{2}/2+\Phi_{\bullet}$ and angular momentum
$|\mbox{\boldmath$L$}|=|\mbox{\boldmath$r$}\times\mbox{\boldmath$v$}|$,
respectively. Energy is measured in units of the mean kinetic energy of field
stars, $K=3\sigma^{2}/2$, and angular momentum in units of
$L_{0}=r_{\epsilon}\,\sigma$. Time is measured in units of Smoluchowski’s
(1916) fluctuation mean-life (3). In contrast to temporary captures, particles
on permanent orbits exhibit negative binding energies ($E<0$) during the
entire time span of the integration.

Since we know the particle ID’s, we can follow the trajectories of permanent
and temporary captures individually. For illustration, Fig. 4 shows the motion
of two field particles captured by the potential $\Phi_{\bullet}$, with orange
and blue lines showing temporary and permanent captures, respectively.
Temporary captures are transient events that inevitably result in the field
particle escaping from the potential $\Phi_{\bullet}$ after a finite amount of
time. In the example plotted in Fig. 4, the specific energy of the field
particle $E=v^{2}/2+\Phi_{\bullet}$ flips from positive to negative $t\simeq
2\,T$ and remains negative until $t\simeq 47\,T$, which yields a survival time
$t_{\rm surv}\sim 45\,T$. After $t\gtrsim 47\,T$, the particle goes through a
short period of rapid energy oscillations in which the orbit becomes bound and
unbound repeatedly over short time-scales, $t_{\rm surv}\sim T$, until it
finally becomes fully unbound and escapes from the volume element $V$ at
$t\sim 60\,T$. In contrast, permanent captures have $E<0$ at all times and
show milder variations of energy. It is important to emphasize that none of
these particles conserve energy or angular momentum
($\mbox{\boldmath$L$}=\mbox{\boldmath$r$}\times\mbox{\boldmath$v$}$), which
greatly complicates an analytical description of their motion. For example,
their orbits are not typically restricted to a constant orbital plane –even
when $\Phi_{\bullet}$ is spherical & static– and exhibit time-varying orbital
parameters (such as peri, apo-centres and orbital eccentricity) that change in
an apparently random fashion over short time scales.

Figure 5: Upper panel: Radial phase-space location of bound particles in a co-
rotating frame, $x^{\prime}$ – $v_{x}^{\prime}$ measured at the end of the
simulation ($t=t_{f}$). Only particles within a narrow plane $|\Delta
y^{\prime}|<0.2\,r_{\epsilon}$ and $|\Delta z^{\prime}|<0.2\,r_{\epsilon}$ are
shown here. The dwarf galaxy centre is located at
$x^{\prime}/r_{\epsilon}\approx-25$. We identify ‘permanent’ captures as those
particles that remain bound for the whole span of the simulation $t_{c}=t_{f}$
(or $t_{\rm bound}=0$) and plot them with red dots, whereas temporary captures
become bound at $t_{\rm bound}>0$ and are plotted in blue. Note that permanent
captures have negative binding energies measured in a co-rotating frame,
$E_{\rm rot}<0$, and are confined within the tidal radius
$|x^{\prime}|<r_{t}\simeq 10\,r_{\epsilon}$, with $r_{t}$ given by Equation
(21). Notice also the lack of permanent captures at very small distances from
the substructure, $r\ll r_{\epsilon}$. Lower panel: Energy and angular
momentum of the particles plotted above given in the same units as in the
right panels of Fig. 4. Black-solid line shows the circular velocity of the
Hernquist substructure. Notice that permanent and temporary particles appear
segregated in the integral-of-motion space. In particular, temporary captures
can be rarely found on circular orbits. Note also that this division blurs
close to the ’fringe’ region, $E\approx 0$.

As noted above, only a fraction of particles with $E<0$ at $t=0$ (the so-
called immediate captures) remain bound for an indefinite amount of time. Fig.
5 shows that only immediate captures that were initially located within the
local tidal radius become permanently bound. This plot shows particles in a
co-rotating coordinate system centred at the substructure, with the
$x^{\prime}$-axis pointing at the host galaxy centre, and the
$z^{\prime}$-axis aligned with the angular vector
${\mbox{\boldmath$\Omega$}}=(0,0,\Omega)$. In this frame, the effective
potential can be written as (c.f. Eq. 12 of Renaud et al. 2011)

$\displaystyle\Phi_{\rm eff}=\Phi_{\bullet}-\frac{1}{2}(\lambda_{1}\,x^{\prime
2}+\lambda_{2}\,y^{\prime 2}+\lambda_{3}\,z^{\prime 2}),$ (19)

where $\lambda_{i}$ are the 3 eigen-values of the local tidal tensor. For a
Dehnen (1993) sphere with $\gamma=0$ they can be expressed analytically using
Eq. (17) of Renaud et al. (2011) as

$\displaystyle(\lambda_{1},\lambda_{2},\lambda_{3})=\frac{GM^{\rm
dSph}}{(R_{\bullet}+c^{\rm
dSph})^{3}}\bigg{(}\frac{3R_{\bullet}}{R_{\bullet}+c^{\rm
dSph}},0,-1\bigg{)}.$ (20)

The tidal radius can be estimated from (20) as

$\displaystyle
r_{t}=\bigg{(}\frac{GM_{\bullet}}{\lambda_{1}}\bigg{)}^{1/3}=\bigg{[}\frac{M_{\bullet}}{3M^{\rm
dSph}}\bigg{(}1+\frac{c^{\rm
dSph}}{R_{\bullet}}\bigg{)}\bigg{]}^{1/3}R_{\bullet},$ (21)

which recovers the academic point-mass case for $c^{\rm dSph}=0$. Inserting
the parameters of the experiments run in this Section into (21) yields
$r_{t}\simeq 190\,{\rm pc}\approx 10\,r_{\epsilon}$.

Dotted and solid lines in the upper panel of Fig. 5 show phase-space surfaces
that obey the conditions $E=v_{x}^{\prime
2}/2+\Phi_{\bullet}(x^{\prime},0,0)=0$ and $E_{\rm rot}=v_{x}^{\prime
2}/2+\Phi_{\rm eff}(x^{\prime},0,0)=0$, respectively, with primes denoting
quantities measured in a co-rotating frame. By definition, temporary captures
(shown in blue) populate a phase-space region where $E<0$ over a wide range of
distances. Interestingly, permanent captures are more confined in phase space.
In particular, we find that all permanent captures obey $E_{\rm rot}<0$ and
can only be found at distances $r<r_{t}\simeq 10\,r_{\epsilon}$ from the
substructure. Notice also the apparent lack of permanent captures at distances
smaller than the thermal critical radius, $|x|\lesssim r_{\epsilon}$. As
discussed below, this is likely due to the finite number of particles in our
models.

Lower panel of Fig. 5 shows the energy ($E$) and angular momentum modulus
($L=|\mbox{\boldmath$r$}\times\mbox{\boldmath$v$}|$) of the particles plotted
in the upper panel. For reference, the angular momentum of circular orbits at
a given energy is shown with a black-solid line. The first noteworthy result
is that permanent and temporary captures occupy different regions of the
integral-of-motion space. In particular, it is rare to find permanent orbits
with low angular momentum, which suggests that orbits with high-eccentricity
do not remain bound to the substructure potential for an indefinite amount of
time. In contrast, temporary captures avoid orbits with low eccentricity. This
result agrees with the numerical experiments published in Paper I, which show
that particles trapped in the potential $\Phi_{\bullet}$ exhibit a ‘super-
thermal’ eccentricity distribution with an excess of particles moving on very
eccentric orbits. Notice also that the division between temporary and
permanent captures blurs in the energy ’fringe’, $|E|\lesssim 0.2\,K$, where
$K=3\sigma^{2}/2$ is the average kinetic energy of field particles.

Given that captured particles do not conserve $E$ or $L$, as shown in Fig. 4,
it is remarkable that temporary and permanent orbits remain locked within
their respective regions as they drift in $E$–$L$ space. This empirical result
may serve as a starting point for follow up theoretical work to better
understand the dynamical mechanisms that lead to gravitational capture.

Figure 6: Density enhancement ($\delta_{\star}=n_{\star}/n$) of particles with
negative binding energies $E<0$ as a function of distance from the
substructure. Top, middle and bottom panels split bound particles according to
whether they are permanent, temporary or any kind of captures (see text).
Horizontal dotted lines mark the local density of the field
$\delta_{\star}=1$. Green dots show the density enhancement derived from the
equilibrium Initial Conditions (ICs) by placing a substructure Hernquist
potential $\Phi_{\bullet}$ at 100 random locations in the galaxy and
identifying particles with $E<0$. Note that lack of permament captures beyond
the the tidal radius $r_{t}$, Equation (21), marked with vertical arrows.
Notice also that the density enhancement profiles do not evolve with time, and
that stellar overdensities ($\delta_{\star}>1$) of field particles mainly
originate from temporary captures. Bottom panel shows that the theoretical
profiles accurately describe the distribution of all bound particles, that is
the superposition of permanent and temporary populations.

Fig. 6 shows the density enhancement of bound particles,
$\delta_{\star}=n_{\star}/n$, as a function of distance from the substructure
for (i) permanent captures (top panel), (ii) temporary captures (middle panel)
and all bound particles (bottom panel) at different snapshots of the
simulation. For ease of reference, the background density $\delta_{\star}=1$
is marked with horizontal dotted lines. As expected, in the static experiments
we find a very weak temporal variation of the profiles, suggesting that the
distribution of bound particles reaches steady state quickly. Focusing first
on the top panel, we find that permanent captures do not generate over-
densities, i.e. $\delta_{\star}<1$. We also observe a sharp truncation of the
profile at $r\approx 10\,r_{\epsilon}$, which roughly corresponds to the value
estimated from (21) (marked with vertical arrows). At small distances,
$r\lesssim r_{\epsilon}$, we find no permanent captures. This is likely a
numerical artifact resulting from finite sampling. To prove this point, we
place a substructure potential $\Phi_{\bullet}$ at 100 random locations within
the galaxy and identify bound particles in the unperturbed models at $t=0$.
The averaged density enhancement profile is shown with green dots. Recall that
permanent captures are already bound at $t=0$, and given that their spatial
distribution shows no temporal variation it is not surprising that these
objects follow a similar profile as the one derived from random sampling the
Initial Conditions (ICs). However, there are two visible discrepancies: at
large radii, $r\gtrsim r_{t}$, the profile derived from the ICs is not
truncated, which indicates that the truncation exhibited by permanent captures
originates from dynamical evolution. A second mismatch can be observed at
small radii, $r\lesssim r_{\epsilon}$, where we find that the profile derived
from the ICs extends well below that of permanent captures. This is because
the IC profile is generated by placing the substructure potential
$\Phi_{\bullet}$ at 100 random locations in the galaxy, which enlarges the
statistical sample of bound field stars, allowing us to measure their density
at smaller radii. Crucially, the density enhancement converges slowly to
unity, $\delta_{\star}\to 1$ in the limit $r/r_{\epsilon}\to 0$, which again
indicates that the presence of permanent captures does not generate a stellar
over-density at the substructure location.

The density enhancement generated by temporary captures is plotted in the
middle panel Fig. 6. In contrast to permanent captures, temporary captures
generate an over-density ($\delta_{\star}>1$) at small distances from the
substructure. We also find that the enhancement profile closely follows the
theoretical curve (4) (black-dashed lines) on scales $r\lesssim r_{\epsilon}$.
However, on larger scales $r\gtrsim r_{\epsilon}$ the density of temporary
captures falls off more quickly than predicted by the statistical theory.
Beyond the tidal radius, $r\gtrsim r_{t}$, the profile of temporary captures
approaches again the theoretical curve (4).

Bottom panel show the profiles generated by all particles with $E<0$ at three
different snapshots of the simulation. In agreement with the theoretical
predictions from §2.1.2, the theoretical profile (4) matches the numerical
result after we consider the entire population of bound particles without
specifying whether captures are temporarily or permanently bound. Notice also
that there is no particular feature that marks the location of the tidal
radius in the distribution of bound stars, and that Equation (4) is accurate
in the vicinity and beyond the tidal radius, which suggests that the
superposition of permanent and temporary orbits conspires to generate a
population of bound stars with a constant phase-space density across the
volume under observation. We come back to this point below.

Figure 7: Upper panel: Distribution of capture times ($t_{c}$) defined as the
time since a particle becomes bound (energy flips from positive to negative)
to the current simulation time ($t$). In substructure models with different
growth rates ($\tau$), we find a substantial number of particles captured at
early times ($t_{\rm bound}\lesssim 2\tau$) that remain bound for the rest of
the integration, the so-called ‘permanent’ captures, which materialize as a
peak in the distribution at long capture times $t_{c}\sim t$, which shifts
towards the right of this plot as the integration time increases. Lower panel:
Fraction of permanent captures within a volume element $r_{\rm out}=190\,{\rm
pc}$ as a function of integration time. The fraction of permanent captures is
roughly constant in the static case ($\tau=0$), and peaks around $t\sim\tau$
for substructures with a time-dependent mass (22). At time-scales $t\gg\tau$,
the three models show the same fraction of permanent particles independently
of $\tau$.

###  3.3 Growing substructures on a circular orbit

In the previous Section, a static substructure potential $\Phi_{\bullet}$ is
instantly injected in a sea of field particles in dynamical equilibrium. This
experimental set up is mathematically convenient because it removes any
explicit time dependence from the analytical estimates, but it has limited
physical applications.

In this Section, we inspect a physically-motivated scenario in which a
substructure has a mass that grows as it moves along a fixed orbit across a
sea of field particles (see discussion in §4). To this aim, we re-run the
’compact’ $N$-body models presented in §3.2 adopting a time-dependent
substructure mass

$\displaystyle M_{\bullet}(t)=M_{\bullet}\,[1-\exp(-t/\tau)],$ (22)

with $M_{\bullet}=10^{6}\,{\rm M}_{\odot}$ and a fixed scale radius
$c_{\bullet}=3.78\,{\rm pc}$. Here, $\tau$ is a parameter that controls the
mass growth rate, $\tau=|\dot{M}_{\bullet}/M_{\bullet}|^{-1}_{t=0}$. The
substructure potential $\Phi_{\bullet}$ vanishes in the limit $t\to 0$ for
$\tau>0$, which allows us to easily modify the number of immediate captures by
choosing different values of the time-scale $\tau$. By construction, in the
limit $t\gg\tau$ the time-dependent model approach the static case explored in
§3.2.

Fig. 7 shows the distribution of capture times, $t_{c}=t-t_{\rm bound}$, where
$t_{\rm bound}$ is the integration time at which a particle becomes bound, in
models with three different mass-growth rates, $\tau/T=0,1$ and $3$. Orange
lines correspond to the same distribution plotted in Fig. 3 for static models
($\tau=0)$. Comparison between the panels from top to bottom shows the
emergence of prominent spike of particles with $t_{c}\approx t$, which is the
telltale of the build-up of a permanent captures trapped in the time-evolving
substructure potential, $\Phi_{\bullet}(t)$. Focusing on the models with a
slow growth rate ($\tau=3\,T$, red lines) we find that the spike is absent
early on during the simulation ($t/T=5$, top panel), and that it only starts
to arise at $t/T\gtrsim 20$ (second panel from the top). At the end of the
simulation at $t=50\,T$ (third panel from the top), we find that all three
models show a clear excess of captures with long capture times, $t_{c}\sim t$,
with the peak of the distribution broadening up for models with slow growth
rates ($\tau>0$). This suggests that while in a static potential ($\tau=0$)
permanent captures are trapped immediately at $t=0$, this is not the case in
substructures with a mass growth (22), which capture particles onto permanent
orbits over an extended period of time, $t\lesssim 2\tau$.

Unfortunately, the theoretical conditions that determine whether a field
particle becomes permanently or temporarily bound to a time-dependent
potential $\Phi_{\bullet}(t)$ are more difficult to study than in the static
case, because the energy measured in a co-rotating frame ($E_{\rm rot})$ is
not a conserved quantity in time-dependent systems. The results shown in Fig.
7 call for a heuristic classification of orbits according the the amount of
time that they remain bound. In particular, in what follows we identify
permanent captures as particles that become bound early on, $t_{\rm bound}\leq
2\tau$, and whose energy $E$ remains negative uninterruptedly until the end of
the simulation888Notice that some temporary captures may also remain bound for
relatively long times. However, Fig. 3 suggests that the fraction of temporary
captures that form within $t_{\rm bound}\leq 2\tau$ is $\lesssim 10\%$. These
are the particles that populate the spikes of the distributions plotted in
Fig. 3. With this empirical definition at hand, it is straightforward to
measure the fraction of bound particles that are captured on permanent orbits
as a function of time, which is shown in the bottom panel. The first
noteworthy result is that the fraction of permanent captures is approximately
constant in the static model999The potential $\Phi_{\bullet}$ moving across
the galaxy heats up the population of field particles. This translates into a
capture rate that decreases with time (see also Fig. 8), and to a slightly
rising fraction of permanent captures at $t\gtrsim 10\,T$., $f_{\rm
perm}=N_{\rm perm}/(N_{\rm perm}+N_{\rm temp})=N_{\rm perm}/N_{b}\simeq 0.35$
As expected, increasing the time-scale $\tau$ systematically decreases the
number of permanent captures at $t=0$. In particular, the model with the
slowest growth ($\tau=3\,T$, red line) contains no permanent captures
initially. As time proceeds, the number of permanent captures grows, peaking
at $t\approx\tau$. Once the substructure stops growing, the fraction of
permanent captures approaches the static value at $t\gg\tau$, independently of
the choice of $\tau$. This suggests that the exact manner in which
$\Phi_{\bullet}$ is inserted into the field does not affect the final
distribution of bound particles, insofar as the captured population has enough
time to reach steady state.

Figure 8: Upper panel: Distribution of capture times ($t_{c}$) defined as the
time since a particle becomes bound (energy flips from positive to negative)
to the current simulation time ($t$). Notice the presence of a substantial
number of particles that become bound at $t\approx 0$, which result in a a
strong peak in the distribution at $t_{c}\approx t$ shifting towards the right
as the integration time increases. Lower panel: Permanent-to-temporary ratio
(or ptt in short) of particles captured within a volume element $r_{\rm
out}=100\,{\rm pc}$ as a function of integration time. Permanent captures are
defined as particles with a capture time $t_{c}\approx t$, i.e. they remain
bound during the whole integration time. In contrast, temporary captures only
remain bound for a fraction of the integration time, i.e. $t_{c}\lesssim t$.
Notice that the ptt ration becomes approximately constant at $t\gtrsim 5,T$,
where $T$ is Smoluchowski’s time scale.

Fig. 8 plots the number of bound particles located within a fixed spherical
radius $r_{\rm out}=190\,{\rm pc}$ from the centre of $\Phi_{\bullet}(t)$ as a
function of time (red lines) for three different values of the time-scale
$\tau$. Recall that in the static case, this choice of $r_{\rm out}$ matches
the local tidal radius (21). We find that the total number of particles with
$E=v^{2}/2+\Phi_{\bullet}(t)<0$ within the volume ($N_{b}$) grows in
proportion to the substructure mass. The experimental values are in excellent
agreement with the values derived from Equation (2) using a time-dependent
mass (22) (black-dashed lines).

Fig. 8 also illustrates the process by which a population of captured
particles reaches steady state. Orange and green lines respectively show the
number of particles in the volume $V$ that become bound ($E$ flips from
positive to negative) and unbound (either $E$ flips from negative to positive,
or the particle leaves the volume under observation) within a time interval
$\Delta t=0.1\,T$. In models where $\Phi_{\bullet}$ is introduced suddenly
($\tau=0$), we find that these two values become approximately equal after a
few snapshots. At this point, steady state has been reached. In models where
the potential $\Phi_{\bullet}$ grows slowly ($\tau>T$) we observe two distinct
regimes: (i) at $t\lesssim\tau$, particles are captured at a rate that
overcomes the number of unbinding events, which leads to a net growth of the
population of captured particles. (ii) On longer time-scales, $t\gtrsim\tau$,
the number of bound particles reaches a steady state, wherein the rate of
capture becomes roughly equal to the loss rate, $\Delta N_{\rm
bound}\approx\Delta N_{\rm unbound}$. From this time on, all captures are
temporarily bound, which leads to no net variation in the average number of
stars with $E<0$. It is worth to highlight that in all models the capture rate
and the loss rate conspire to yield an average number of bound particles
within the volume $V$ that is very close to the value derived from Equation
(2) with a time-dependent mass (22) at all times.

Blue-dotted lines in Fig. 8 show the number of permanent captures as a
function of time. Notice that all particles trapped onto permanent orbits are
captured early on, $t\lesssim 2\tau$, and that at the end of the simulation
$t_{f}=50\,T$ all models contain a similar number of permanent captures
independently of the choice of $\tau$. At early times, the build up of
permanently-bound orbits leads to a capture rate that exceeds the loss rate,
$\Delta N_{\rm bound}>\Delta N_{\rm unbound}$, and thus to a net growth of the
average number of bound particles, $N_{b}(t)$. This may suggest that permanent
captures result from an impulsive response of the galactic field to the growth
of the substructure potential on short time-scales
$|\dot{\Phi}_{\bullet}/\Phi_{\bullet}|^{-1}=|\dot{M}_{\bullet}/M_{\bullet}|^{-1}_{t=0}=\tau\sim
T$. Further theoretical work is need to understand how this capture mechanism
works.

Figure 9: Upper panel: Average phase-space density
$Q_{\star}=n_{\star}/\sigma^{3}_{\star}=n\,\delta_{\star}/\sigma^{3}_{\star}$
as a function of distance from the potential $\Phi_{\bullet}$ measured from
the models shown in Fig. 8. Distances are measured in units of the thermal
critical radius $r_{\epsilon}$, Equation (6). Note that $Q_{\star}$ remains
approximately constant within the volume element, and that its value agrees
well with the theoretical expectation (11) (black-dashed line). Middle panel:
Density enhancement profile $\delta_{\star}(r)=n_{\star}(r)/n$. The black-
dashed line shows the theoretical curve (4). For reference, we mark with a
dotted line the background density value $\delta_{\star}=1$. Lower panel:
Velocity dispersion profiles of the models shown above normalized by the local
velocity dispersion of the field, $\sigma_{\star}(r)/\sigma$. Black-dashed
line shows the theoretical profile (5).

Fig. 9 shows that the final distribution of bound particles around the
potential $\Phi_{\bullet}$ at $t=t_{f}$ is not sensitive to the rate of growth
of the substructure potential ($\tau$). This result holds insofar as the
simulation time is longer than the time-scale on which the potential grows,
i.e. $t_{f}\gg\tau$. Let us focus first on the upper panel, which plots the
average phase-space density of captured particles as a function of distance to
the substructure, $Q_{\star}(r)\equiv n_{\star}(r)/\sigma^{3}_{\star}(r)$. As
expected, we find that $Q_{\star}$ is approximately constant across the volume
$V$, and that the measured value of $Q_{\star}$ approximately matches Equation
(11) (marked with a black-dashed line) independently of the choice of $\tau$.
Section 2.1.2 predicts that bound particles with a constant phase-space
density in steady state follow the density profile (4) (black-dashed line).
The middle panel shows that this expectation is largely correct. Regardless of
whether $\Phi_{\bullet}$ is introduced gradually ($\tau>0)$ or suddenly
($\tau=0$), we find that the number density of captured particles scales as
$\delta_{\star}\sim|\Phi_{\bullet}|^{3/2}$. Similarly, the lower panel of Fig.
9 shows that the velocity dispersion profile goes as
$\sigma_{\star}\sim|\Phi_{\bullet}|^{1/2}$, which in the case of a Hernquist
substructure potential $\Phi_{\bullet}$ leads to the velocity dispersion
profile (5) plotted with a black-dashed line.

Figure 10:  Left panels: Initial distribution of field particles (grey dots)
in dynamical equilibrium within a Dehnen potential with a total mass $M^{\rm
dSph}=3\times 10^{9}\,{\rm M}_{\odot}$, a scale radius $c^{\rm dSph}=2\,{\rm
kpc}$ and a central density slope $\gamma^{\rm dSph}=1$. Right panels:
Distribution of particles that become bound to a Hernquist potential with
$M_{\bullet}=3\times 10^{7}\,{\rm M}_{\odot}$ and a compactness $\kappa=+0.9$
(see text). Particles that become captured during the integration time and
have $E<0$ ($E>0$) at the final snapshot ($t=t_{f}$) are highlighted with red
(black) dots. The trajectory of the substructure across the galaxy is shown
with a solid line. Black circles mark a volume size $r_{\rm out}=300\,{\rm
pc}$ around the substructure potential. Inlets in the upper-right corner zooms
in at the substructure location (marked with a cyan dot for ease of
reference). Distances are given in $\,{\rm kpc}$. Notice that particles
captured from a truncated galaxy exhibit a flattened shape roughly oriented
along the orbital motion.

###  3.4 Accreted substructures on an eccentric orbit

In the previous Sections, we carry out experiments where a substructure is
injected on a circular orbit either suddenly ($\tau=0$), fast ($\tau=T$), or
slowly ($\tau>T)$ into a pre-existing field of field particles in dynamical
equilibrium. This experiminental setup invariably leads to a population of
field particles that are captured onto permanent orbits around the potential
$\Phi_{\bullet}(t)$. In this Section, we present models where a static
substructure potential is initially placed in a region of the galaxy devoid of
field stars on an eccentric orbit. By design, these experiments contain no
immediate or permanent captures.

The numerical set-up is slightly different from previous Sections. Here, we
consider a dwarf galaxy sourcing a cuspy Dehnen (1993) potential ($\gamma=1$)
with a mass $M^{\rm dSph}=3\times 10^{9}\,{\rm M}_{\odot}$ and a scale radius
$c^{\rm dSph}=2\,{\rm kpc}$. For illustration, we show two experiments where
the substructure is inserted in the galactic potential at apocentre, $R_{\rm
apo}=1380\,{\rm pc}$. The initial velocity is set to $V_{\bullet}=7.9\,{\rm
km\,s^{-1}}$, which leads to a pericentre of $R_{\rm peri}=165\,{\rm pc}$.
Thus, this orbit penetrates the inner-most regions of the galaxy. We then
follow the motion of the substructure for two dynamical times measured at
apocentre, $\Omega^{-1}(R_{\rm apo})=R_{\rm apo}/V_{c}(R_{\rm apo})=34\,{\rm
Myr}$, where $V_{c}(R_{\rm apo})=39.4\,{\rm km\,s^{-1}}$ is the circular
velocity. Before showing the outcome of this experiment, it should be pointed
out that we have explored several combinations of apo- and peri-centres and
number of orbital periods, finding similar results as discussed below.

We consider two populations of tracer particles in equilibrium within the
galactic potential (see §3.1 for details). The first population follows a
‘truncated’ density profile $(\alpha_{f},\beta_{f},\gamma_{f})=(2,30,0.2)$
with a scale radius $R_{0}=460\,{\rm pc}$, which falls off very steeply beyond
$R\gtrsim 600\,{\rm pc}$. The luminosity-averaged velocity dispersion of these
particles is $\langle v^{2}\rangle^{1/2}/3=11.7\,{\rm km\,s^{-1}}$. The second
population follows an ‘extended’ Plummer-like profile with
$(\alpha_{f},\beta_{f},\gamma_{f})=(2,5,0.1)$ and the same scale radius,
$R_{0}=460\,{\rm pc}$, which leads to a slightly higher luminosity-averaged
velocity dispersion $\langle v^{2}\rangle^{1/2}/3=19.9\,{\rm km\,s^{-1}}$.

The substructure sources a static Hernquist potential with a mass
$M_{\bullet}=3\times 10^{7}\,{\rm M}_{\odot}$. To set the scale radii of the
models we use the thermal critical radius at $t=0$ as a reference. For
illustration, we set the compactness to $\kappa=+0.9$ in both galaxy models,
which leads to a scale radius
$c_{\bullet}=(1-\kappa)r_{\epsilon}=0.1\,r_{\epsilon}$. To find the thermal
critical radii of the two models we insert $\sigma=\langle
v^{2}\rangle^{1/2}/3$ and $V_{\bullet}$ in Equation (6), which yields
$r_{\epsilon}=390\,{\rm pc}$ and $980\,{\rm pc}$ for the extended and
truncated galaxy models, respectively.

Left panels of Fig. 10 show the projected locations of tracer particles with
extended (top panel) and truncated (bottom panel) profiles in equilibrium
($t=0$). In the right panels, we overplot the trajectory of the substructure
as it falls into the galaxy on an eccentric orbit (solid lines), and mark its
final location with a cyan dot. Red dots show particles with $E<0$ at the
final snapshot $(t_{f})$. Black dots correspond to particles that have been
captured for a finite amount of time before being lost to galactic tides.

As expected, zooming in at the substructure location we find that particles
with $E<0$ generate a local overdensity of field particles. Interestingly,
their observed spatial distribution depends strongly on the initial profile of
the field. In the extended galaxy model (top panel), the over-density of bound
particles has a close-to-spherical shape and a relatively large size. In
contrast, the over-density in the bottom panel has a smaller size and a
flattened shape, which appears to be roughly oriented along the orbital
motion.

Figure 11: Number of bound particles located within a volume size $r_{\rm
out}=0.3\,{\rm kpc}$ as a function of time for the models plotted in Fig. 10.
Black-dashed lines show the theoretical prediction from Equation (2), with the
time-dependent number density and velocity dispersion of field particles
measured at the location of the substructure, $n(t)=n[R(t)]$ and
$\sigma(t)=\sigma[R(t)]$. Pericentric passages are marked with vertical dotted
lines.

The rate at which field particles are captured in these two models is also
strikingly different. Fig. 11 shows the number of bound particles within a
volume size $r_{\rm out}=300\,{\rm pc}$ (marked with black circles in Fig. 10
for reference) as a function time. The top panel shows that the substructure
potential orbiting in the extended galaxy contains $N_{b}\sim 3\times 10^{3}$
at $t=0$. These are immediate captures by definition. As the substructure
falls into the inner regions of the galaxy, the number of bound particles
grows slightly, suffering a significant increase after the first pericentric
passage. At the last snapshot of the simulation, the number of bound particles
has increased by a factor $\sim 10$ up to $N_{b}(t_{f})\sim 4\times 10^{4}$.
Remarkably, the theoretical expectation given by Equation (2) (black-dashed
line) barely changes along the orbital phase in spite of the large variation
of galactocentric distances and velocities. This result reflects the counter
balance between the variation of the mean phase-space density of field
particles along the orbit, $Q(R)=n(R)/\sigma^{3}(R)$, and the exponential
suppression of captures at high velocities,
$\xi=\exp[-V^{2}_{\bullet}/(2\sigma^{2}(R))]$. As the substructure plunges
into the inner regions of the galaxy the orbital velocity increases, making
capture less efficient. Simultaneously, it encounters field particles with
higher phase-space densities, with both effects nearly balancing out in
Equation (2). The reverse process takes place from pericentre to apocentre.

In contrast, the substructure orbiting in the truncated galaxy does not
capture field particles until it is close to its first orbital pericentre at
$\Omega\,t\simeq 1.2$. Bottom panel of Fig. 11 shows a rapid increase in the
population of bound particles, which grows $\sim 3$ orders of magnitude after
the pericentric passage. A large fraction of these particles are quickly lost
to galactic tides, such that after its first apogalacticon the number of bound
particles plateaus at $N_{b}\sim 10^{2}$. As it reaches its second pericentre,
it grows up again to $N_{b}\sim 5\times 10^{3}$, falling back to $N_{b}\sim
10^{2}$ at its second apocentre. Comparison with the theoretical expectation
(black-dashed lines) shows that Equation (2) cannot accurately describe this
behaviour. There are two main reasons for the mismatch. In the first place,
the statistical theory derived in §2.1 relies on the local approximation,
which assumes that the properties of the galactic field do not vary strongly
as a function of distance from the substructure. This is clearly not the case
in this particular galaxy model, which is strongly truncated beyond $R\gtrsim
600\,{\rm pc}$. Second, it also assumes that the properties of the field are
time-invariant, and that the population of captured particles is in steady
state at all times. None of these conditions apply to the these particular
models, which are strongly time dependent. The latter approximation is
particularly poor after the first pericentric passage. Here, the numerical
experiments show that the substructure has captured $N_{b}\sim 10^{3}$
particles after the first orbital revolution. In contrast, Equation (2)
predicts $N_{b}=0$ because the density of field particles at this location is
$n(R_{\rm apo})=0$ (see bottom-left panel of Fig. 10). This mismatch repeats
again at the second apocentre. Clearly, we observe field particles in the
vicinity of the substructure in the bottom-right panel of Fig. 10, which are
not predicted by Equation (2. These particles are either still bound to the
substructure (red dots), or were ‘scooped’ from the inner regions of the
galaxy and released from the substructure by galactic tides (black dots).

Figure 12: Number density profile of field particles bound to the substructure
models shown in Fig. 10. Distances are measured relative to the final
substructure location and given in units of the substructure scale radius
($c_{\bullet}$). The profiles are averaged over 500 realizations of the model
to reduce statistical noise (see text). Notice that the model associated with
an extended galaxy model roughly follows the theoretical profile (4). In
contrast, substructures evolving in a truncated galaxy generate a local
overdensity that falls off steeply at large distances $\delta_{\star}\sim
r^{-3}$ at $r\gg c_{\bullet}$.

Fig. 12 shows that distribution of bound particles in the extended and
truncated galaxies plotted in Fig. 10 follow very different profiles. For ease
of comparison, we normalize the profiles to its central value,
$n_{0}=n_{\star}(r=0)$, and measure distances in units of the substructure
scale radius ($c_{\bullet}$). Given the relative small number of particles
with $E<0$ at the final snapshot, we reduce statistical noise by generating
200 and 600 realizations of the extended and truncated galaxy models,
respectively. This is done placing substructures at random positions over the
surface of a sphere with a radius $R_{\rm apo}=1380\,{\rm pc}$ and a
tangential velocity $V_{\bullet}=7.9\,{\rm km\,s^{-1}}$ at $t=0$, such that
all individual realizations have a common pericentre of $R_{\rm
peri}=165\,{\rm pc}$.

The substructure orbiting in an extended galaxy contains a population of
captured field particles with a profile that roughly follows Equation (4).
This is a remarkable result considering the strong strong temporal evolution
of the field along the substructure orbit, and that the statistical theory is
derived for static systems. The good agreement suggests the (local) steady-
state expressions are accurate at leading order insofar as the orbital period
around the host galaxy is longer than the internal crossing time of the
substructure. In our models, the crossing time
($T_{\epsilon}=r_{\epsilon}/\langle v^{2}\rangle^{1/2}\simeq 1.5\,{\rm Myr}$)
is much shorter than their orbital time in the dSph potential ($t_{\rm
dyn}=R_{\bullet}/V_{\bullet}\sim 30\,{\rm Myr}$), which suggests that the
steady-state equations are reasonably accurate far from orbital pericentre.

In contrast, field particles bound to the substructure orbiting in an extended
galaxy show a profile that falls off more steeply than predicted by Equation
(4) at large distances, whereas at small distances, $r\ll c_{\bullet}$, it
matches the profile of the extended galaxy model.

In general, the numerical experiments shown in this Section indicate that the
distribution of stars captured by moving substructures depend on whether these
objects are form in a pre-existing sea of field particles, or are accreted
from outside the galaxy.

Figure 13: Upper left panels: Locations projected onto the orbital plane of
stellar tracer-particles temporarily-bound to a substructure with a mass
$M_{\bullet}=10^{6}\,{\rm M}_{\odot}$ moving on a circular orbit at a fixed
distance $R=0.5\,{\rm kpc}$ from the centre of the Fornax dSph (see §3.1).
These models adopt two different scale-radii (marked with black-solid circles
for reference): ’compact’ ($\kappa=+0.8$, or $c_{\bullet}=0.2\,r_{\epsilon}$
(upper panels) and ‘fluffy’ ($\kappa=-0.5$, or $c_{\bullet}=1.5\,r_{\epsilon}$
( lower panels), where $r_{\epsilon}=18.9\,{\rm pc}$ is the thermal critical
radius measured from (6) (black-dashed circle). Middle left panels: Velocity
dispersion map for the models shown in the left panels. Middle right panels:
Density enhancement as a function of distance to the substructure centre. Red
and blue arrows mark the location of the thermal critical radius and the scale
length of the substructure, respectively. The theoretical profile predicted by
Equation (4) is shown with a black line. Right panels: Velocity dispersion
profiles in spherical coordinates of the above models. Black lines show the
profiles predicted by Equation (5). Captured particles have isotropic
velocities on account of their chaotic motions in the dwarf potential (see
text).

###  3.5 Compact vs. fluffy substructures

The statistical theory outlined Section 2.1 predicts that only dark subhaloes
that are sufficiently compact ($\kappa>0$) can become visible as localized
overdensities of field stars. To illustrate this result, this Section presents
similar experiments as shown in §3.2 with positive and negative compactness
values.

Fig. 13 shows the projection onto the orbital plane of particles bound to a
Hernquist substructure potential with a mass $M_{\bullet}=10^{6}\,{\rm
M}_{\odot}$ moving on a circular orbit $R_{\bullet}=0.5\,{\rm kpc}$ in the
Fornax dSph-like galaxy model introduced in §3.1. Recall that the tidal radius
given by Equation (21) is $r_{t}\approx 190\,{\rm pc}$. We consider two models
with $\kappa=+0.8$ (upper panels) and $\kappa=-0.5$ (lower panels), with the
negative sign of $\kappa$ implying that the number density of captured field
particle does not exceed the background ($\delta_{\star}<1)$ at any distance
from the substructure. Application of Equation (6) yields a thermal critical
radius $r_{\epsilon}=18.9\,{\rm pc}$, which translates into scale-radii
$c_{\bullet}=3.78\,{\rm pc}$ and $c_{\bullet}=28.3\,{\rm pc}$ for the
‘compact’ and ‘fluffy’ models, respectively. Notice that the time-scale needed
to cross the over-density size by a random field particle is much shorter than
the orbital time around the dSph galaxy potential,
$T_{\epsilon}=r_{\epsilon}/\sigma^{\rm dSph}=18.9\,{\rm pc}/(12.2\,{\rm
km\,s^{-1}})\approx 1.5\,{\rm Myr}$. As a result, the population of captured
particles reaches steady state quickly. We choose a volume size around the
substructure equal to its tidal radius (21), $r_{\rm out}=r_{t}=190\,{\rm
pc}$. We run both models from $t=0$ until $t=50\,T$.

Fig. 13 highlights the impact of the substructure compactness
($\kappa=1-c_{\bullet}/r_{\epsilon}$) on the spatial and kinematic
distribution of stars with $E<0$. As predicted in §2.1, at a fixed mass only
compact ($\kappa>0$) objects lead to the formation of over-dense regions of
trapped stars with a velocity dispersion that is systematically hotter than in
fluffy ($\kappa<0$) counterparts. These differences can be better quantified
in the middle-right and right-most panels, which compare the density
enhancement (4) and velocity dispersion (5) profiles derived analytically
(black-solid lines) against the numerical values (orange dots). For reference,
Fig. 13 marks with vertical arrows the location of $c_{\bullet}$ (in blue) and
$r_{\epsilon}$ (in red), respectively. At large distances $r\gg c_{\bullet}$,
we find that the density enhancement approaches the Keplerian behaviour
$\delta_{\star}\sim r^{-3/2}$ and the velocity dispersion drops off as
$\sigma_{\star}\sim r^{-1/2}$, whereas at small distances $r\ll c_{\bullet}$
both profiles become roughly constant. However, while in the compact model
($\kappa=+0.8$) trapped stars reach a relatively high density,
$\delta_{\star,0}=\delta_{\star}(r=0)\approx 9.0$, and a hot velocity
dispersion, $\sigma_{\star,0}/\sigma=\sigma_{\star}(r=0)/\sigma\simeq 1.9$,
stars captured by the fluffy counterpart with a negative compactness
($\kappa=-0.5$) show subdominant densities, $\delta_{\star,0}=0.44$, and a
cold dispersion, $\sigma_{\star,0}/\sigma\simeq 0.69$, with respect to the
field. We discuss the relevance of this result for the detection of dark
substructures in §4.

An important property of the models plotted in Fig. 13 is that they appear
approximately spherical and show no signature of tidal tails. This is in stark
contrast with the substructures moving on eccentric orbits plotted in Fig. 10,
which are elongated and exhibit prominent tidal tails roughly aligned with the
orbit of the substructure. This emphasizes the need to model the motion of the
substructure in the host potential and the distribution of field particles
simultaneously, as discussed in the following Section.

##  4 Discussion

###  4.1 Anomalous stellar clusters in dSphs

Several Milky Way dSphs contain stellar substructures that may be good
candidates for being composed of field stars captured by dark substructures,
i.e. they have a large size for their luminosity, their stellar populations
are indistinguishable from those of the host galaxy, and they exhibit DM-
dominated mass-to-light ratios.

For example, the Fornax dSph is known to contain 5 massive ($\sim 10^{5}\,{\rm
M}_{\odot}$) globular clusters (e.g. Larsen et al. 2012). There has been a
long-standing debate in the literature about a possible sixth cluster, named
Fornax 6 (F6), which was first noted by Shapley (1939). While early studies
debated whether it was composed of stars or background galaxies (e.g. Stetson
et al. 1998), recent ground-based and Gaia data show that F6 is clearly an
overdensity of stars (Wang et al. 2019). F6 has total stellar mass of
$M_{\star}^{F6}=(7.2\pm 2.2)\times 10^{3}\,{\rm M}_{\odot}$, and a half-light
radius of $r_{h}^{F6}=11\pm 1.4\,{\rm pc}$ (Wang et al. 2019). Crucially, the
metallicity $\mathrm{[Fe/H]}=-0.71\pm 0.05$ and age $\sim 2\,{\rm Gyr}$ of F6
are very similar to the average metallicity and age of Fornax’s MR stars (Pace
et al. 2021). This is in stark contrast with the other 5 globular clusters,
which have much lower metallicities that range between
$-2.5<\mathrm{[Fe/H]}<-1.4$. With this data at hand F6 members are therefore
indistinguishable from MR stars in the Fornax dSph. Importantly, the velocity
dispersion of F6 is unexpectedly high, $\sigma^{F6}=5.6\pm 2.0\,{\rm
km\,s^{-1}}$, which translates into an ‘anomalous’ mass-to-light ratio of
$15<M/L<258$ when virial equilibrium is assumed (Pace et al. 2021). The
inflated mass-to-light ratio has been interpreted as a result of on-going
tidal disruption, although no tidal tails originating from the cluster have
been found.

Alternatively, the unusual properties of F6 can be explained if this stellar
system is made of field stars temporarily captured by a dark substructure
orbiting in the Fornax dSph potential. According to Fig. 2 a stellar over-
density with $N_{\star}\sim 10^{4}$ members and a mass-to-light ratio of
$M/L\sim 100$ would require a minimum subhalo mass of $M_{\bullet}\sim
10^{6}\,{\rm M}_{\odot}$, which lies significantly below the minimum subhalo
mass that can trigger star formation, and is consistent with our working
assumption that this object would be ‘dark’. To be visible as an over-density,
it must be compact enough, i.e. $\kappa>0$ or $c_{\bullet}\lesssim
r_{\epsilon}\sim 20\,{\rm pc}$ (for illustrative examples, see the numerical
models presented §3.1). Interestingly, this size is smaller than the average
peak-velocity radius of field CDM haloes of comparable masses, i.e. $r_{\rm
max}\sim 100\,{\rm pc}$ for $M_{200}\sim 10^{6}\,{\rm M}_{\odot}$ at redshift
$z=0$. We discuss the implications of this estimate in Section 4.3.

In addition, de Boer et al. (2013) report the presence of two stellar
overdensities in the Fornax dSph with lower luminosities than F6. These
stellar systems also have elongated shapes (see their Fig. 2a), high
metallicities ($\mathrm{[Fe/H]}\simeq-0.6$) and young stellar ages ($\sim
1.5\,{\rm Gyr}$), similar to those of the Metal Rich component of the Fornax
dSph. VST observations confirmed the presence of those substructures, and
revealed two additional overdensities that share similar properties (Bate et
al. 2015). More recently, using DES data Wang et al. (2020) find several high-
density regions of high-metallicity & young stellar ages, some of which were
new and some of which were previously known (Coleman et al. 2004). These
substructures are therefore good candidates for being made of MR field stars
captured by a population of massive dark objects. Spectroscopic measurements
of their systemic velocity and proper motions, together with their internal
velocity dispersion are needed in order to test this scenario.

The Eridanus II dSph is a dwarf spheroidal galaxy with low luminosity,
$N_{\star}^{\rm dSph}\sim(6\pm 1)\times 10^{4}$ (Crnojević et al. 2016), a
half-light radius of $r_{h}^{\rm dSph}=299\pm 12\,{\rm pc}$ (Simon et al.
2021), and low velocity dispersion $\sigma^{\rm dSph}=6.9\pm 1.1\,{\rm
km\,s^{-1}}$ (Li et al. 2017). Surprisingly for such a faint galaxy, it
contains a single stellar cluster located at $R=23\pm 3\,{\rm pc}$ from its
centre (Koposov et al. 2015). The cluster has very unusual properties, e.g. it
is ancient $\sim 13.5\pm 0.3\,{\rm Gyr}$, and extremely metal-poor
$\mathrm{[Fe/H]}=-2.75\pm 0.2$ (Weisz et al. 2023). It has a remarkably large
size, $r_{h}=15\pm 1\,{\rm pc}$, for a cluster with luminosity $M_{V}=-2.7\pm
0.3$ (Simon et al. 2021) (note that Crnojević et al. 2016 derive slightly
different values). Like F6, Eri II cluster is elongated, $\epsilon\approx
0.31\pm 0.05$, with a shape that is remarkably well aligned with the
ellipticity of the Eri II galaxy (Simon et al. 2021). An elongated shape could
be the telltale signature of these systems being close to full tidal
destruction and therefore in a disequilibrium state. Yet, $N$-body models that
adopt this scenario struggle to reproduce the extended size of these clusters,
and predict bright tidal tails that should have already been observed in the
data currently available (Orkney et al. 2022). Crucially in the context of
this paper, the age and metallicity of the Eri II cluster are statistically
indistinguishable from the dSph itself (Crnojević et al. 2016; Simon et al.
2021; Weisz et al. 2023), which makes this system a good candidate for being
composed of field stars captured by a dark substructure. Again, spectroscopic
follow-up measurements are needed to test the main prediction from this
scenario, i.e whether the Eridanus II cluster exhibits anomalously high
velocity dispersion compatible with being a DM-dominated system.

As a cautionary note, it must be pointed out that the nature of some of those
stellar substructures is still debated. For example, there were early claims
that F6 is not a stellar cluster, but actually an overdensity of unresolved
background galaxies (Verner et al. 1981; Demers et al. 1994; Stetson et al.
1998). Only recently, DECam imaging and Gaia astrometric data have shown that
F6 is clearly an overdensity of stars (Wang et al. 2019), a conclusion later
supported by Magellan/M2FS spectroscopy, which helped to identify $\sim
15$–$17$ likely stellar members (Pace et al. 2021). Yet, the relatively small
size of the data set may be still prone to systematics arising from low-number
statistics.

The unknown physical separation between these substructures and their host
galaxies introduces additional uncertainty in the models. For simplicity, it
is common to assume that these objects are located at a galactocentric
distance equal to the projected separation on the sky, which is the minimum
distance allowed by the data. Given that the phase-space density of field
stars decreases with distance to the host galaxy centre, this choice can
potentially bias future model constraints. For example, were these
substructures located in the outskirts of the dwarf galaxies, substructures
would need to have systematically larger masses in order to capture the same
number of field stars (see §2).

The results enclosed in this paper call for follow-up observational efforts to
measure the spatial distribution, kinematics & chemical composition of stellar
substructures in dSphs with better accuracy, as they might provide unique
constraints on the population of dark objects orbiting in these galaxies, as
well as on the particle nature of dark matter, as briefly discussed in §4.3.

###  4.2 Theory: simplifications and current uncertainties

The theoretical models presented in this work rely on a number of assumptions
that are worth discussing here.

First, our substructure models source a spherical potential. Yet, anomalous
stellar systems detected in dSphs, such as F6 and the lone cluster in Eridanus
II, tend to exhibit elongated shapes that are well aligned with the morphology
of the host galaxy (Wang et al. 2019; Simon et al. 2021). In the scenario
proposed in this paper, stars trapped in a dark substructure may appear
aspherical due to an intrinsic triaxial shape of the objects that captured
them, and/or a triaxial dwarf galaxy potential (e.g. Allgood et al. 2006;
Despali et al. 2014). Interestingly, the numerical experiments plotted in Fig.
10 show that spherical substructures accreted onto the galaxy on eccentric
orbits may also generate flattend spatial overdensities that appear to be
aligned with the orbital motion. Further theoretical work is needed to
understand the possible connection the elongated shape of these stellar
clusters and the expected triaxiality of dark matter haloes as well as their
formation history.

The statistical experiments in §3 adopt substructure models with masses that
are either constant, or grow with time. The latter models can be applied for
example to study the population of field particles trapped around intermediate
mass black holes (IMBHs) that grow in disks surrounding supermassive black
holes (e.g. McKernan et al. 2012). However, our models do not cover the case
of self-gravitating subhaloes orbiting in a parent halo, which experience
tidal stripping and mass loss after each pericentre passage (e.g. Peñarrubia
et al. 2010). This shortcoming may be addressed by running live $N$-body
simulations where substructures are modelled as self-gravitating objects.
Unfortunately, the numerical tools required to explore this scenario are more
complex and computationally expensive than the ones used for this work. The
results enclosed in this paper suggest that the analytical expressions
outlined in §2 can be applied to time-dependent substructures,
$\Phi_{\bullet}(t)$, insofar as its time-evolution is not impulsive, i.e.t it
occurs on time scales $|\dot{\Phi}_{\bullet}/\Phi_{\bullet}|^{-1}\gtrsim T$,
and the local approximation is reasonably accurate. Under these conditions,
the population of bound field stars can be assumed to be in steady state in a
time-varying potential $\Phi_{\bullet}(t)$, see Fig. 8. For substructures
moving on eccentric orbits around extended galaxies, numerical experiments in
§3.4) show that the steady-state approximation may be reasonably accurate away
from orbital pericentre. In galaxies with steep density profiles, it is the
local approximation that fails away from orbital pericentre.

Our experiments do not explore a cosmologically-motivated scenario in which
field stars form in a dark matter halo that contains a pre-existing population
of dark matter clumps. By definition, dark substructures orbiting in these
galaxies would not host a population of permanent captures, yet they may be
able to trap field stars onto temporary orbits. Given the numerical
experiments shown in §3, we expect a different distribution of captured field
stars compared to models where substructures grow in a pre-existing sea of
particles. Cosmological hydrodynamical simulations that incorporate star
formation/feedback are needed to tackle this issue.

The number, distribution and density profiles of dark matter sub-subhaloes in
the satellite galaxies of MW-like haloes are notoriously uncertain. For
illustrative purposes, it is useful to estimate how many subhaloes that one
would expect in a Fornax-like dwarf spheroidal with a virial mass of $M_{\rm
vir}\sim 10^{9}\,{\rm M}_{\odot}$ (e.g. Peñarrubia et al. 2008a; Errani et al.
2018). Adopting a mass ratio between the Fornax dSph and the MW of $\sim
10^{9}\,{\rm M}_{\odot}/10^{12}\,{\rm M}_{\odot}=10^{-3}$, and re-scaling CDM
haloes of MW-like galaxies down three orders of magnitude suggests that the
Fornax dSph should contain of the order of $\sim 100$ satellite sub-subhaloes
with $M_{\bullet}\gtrsim 10^{6}\,{\rm M}_{\odot}$ enclosed within its virial
radius, $r^{\rm dSph}_{\rm vir}\sim 30\,{\rm kpc}$ (e.g. Weerasooriya et al.
2023 and references therein).

However, this naive scaling is only valid for field haloes. As pointed out
above, dwarf spheroidals tend to be accreted early onto larger galaxies and
lose a large fraction of their mass to Galactic tides (e.g. Peñarrubia &
Benson 2005; Errani et al. 2017), which possibly also removes a large fraction
of the sub-subhalo population that fell in embedded with the dSph,
particularly those with large orbital apocentres and long orbital periods (e.g
Peñarrubia et al. 2008b; Errani & Navarro 2021).

At present, it is not possible to turn to cosmological simulations of
structure formation for accurate predictions on the number of dark sub-
subhaloes in dSphs. Current $N$-body methods struggle to resolve objects with
internal crossing times as short as $T\sim 1\,{\rm Myr}$ in Milky Way-like
haloes. Furthermore, these simulations suffer from well-known numerical
artifacts (such as self-heating and artificial disruption) on scales
comparable to the resolution of the simulation (e.g. van den Bosch & Ogiya
2018; Errani & Peñarrubia 2020). These issues call for dedicated high-
resolution $N$-body simulations that address these numerical shortcomings and
provide reliable statistics of ensembles of sub-subhaloes in dwarf spheroidal
galaxies.

Alternatively, since the orbits of MW dSphs are relatively well known thanks
to the astrometry provided by the Gaia mission (e.g. Battaglia et al. 2022 and
references therein), it is now feasible to model the properties of the
surviving population of dark substructures in each individual MW dSph using
constrained $N$-body simulations that do not suffer from the above-mentioned
numerical artifacts. We leave these questions to future work.

###  4.3 Dark Matter particle constraints

The detection and characterization of dark sub-subhaloes in dSphs may provide
unprecedented constraints on the mass and self-interacting cross-section of DM
particle models. For example, the vast majority of gaps in cold tidal streams
around MW-like galaxies are sensitive to subhaloes with masses in the range
$10^{6}<M_{\bullet}/\,{\rm M}_{\odot}<10^{8}$ (Erkal et al. 2016), while
strong-lens observations are sensitive to subhaloes with masses
$M_{\bullet}\gtrsim 6\times 10^{8}\,{\rm M}_{\odot}$ (e.g. O’Riordan et al.
2023). In contrast, Figs. 1 and 2 show that dSphs may be able to probe subhalo
masses down to $M_{\bullet}\sim 10^{4}$–$10^{5}\,{\rm M}_{\odot}$.

Measuring the masses of dark sub-subhaloes down to those scales could
significantly improve current bounds on the free-streaming length of DM
particles. This is best illustrated by quantifying the half-mode mass
$m_{hm}$, defined as the mass scale where the WDM power spectrum is suppressed
by half with respect to CDM models (Schneider et al. 2012)

$m_{\rm hm}=3\times 10^{8}\,{\rm M}_{\odot}\bigg{(}\frac{m_{\rm DM}}{3.3\,{\rm
keV}}\bigg{)}^{-3.33}.$

Detecting field stars trapped in a dark subhalo with a mass of
$M_{\bullet}\simeq 10^{6}\,{\rm M}_{\odot}\gtrsim m_{\rm hm}$ would imply a
half-mass mode of $m_{\rm hm}\gtrsim 18\,{\rm keV}$, significantly tightening
existing constraints from gravitational lensing ($m_{\rm DM}\gtrsim 5\,{\rm
keV}$; Gilman et al. 2020), Lyman-alpha forest ($m_{\rm DM}\gtrsim 3\,{\rm
keV}$; Villasenor et al. 2022), or the combination of strong gravitational
lensing, the Ly-$\alpha$ forest, and the number of luminous satellites in the
Milky Way, which put a lower particle mass limit of $m_{\rm DM}\gtrsim 6\,{\rm
keV}$ (Enzi et al. 2021). Detecting objects with masses $M_{\bullet}\simeq
10^{5}\,{\rm M}_{\odot}\gtrsim m_{\rm hm}$ would push the constraints up to
$m_{\rm hm}\gtrsim 37\,{\rm keV}$.

In addition, the internal structure of these systems may also reveal whether
DM particles experience self-interactions. For example, dark substructures
with $M_{\bullet}\sim 10^{6}\,{\rm M}_{\odot}$ must be sufficiently compact,
$c_{\bullet}\lesssim 20\,{\rm pc}$, in order to generate a stellar overdensity
of field stars in a Fornax-like dSph (see middle panel of Fig. 2). This
condition translates into a characteristic density
$\rho_{\bullet}=M_{\bullet}/(2\pi c_{\bullet}^{3})\gtrsim 20\,{\rm
M}_{\odot}\,{\rm pc}^{-3}$, which is higher than CDM subhaloes with similar
masses. More precisely, field halos with a mass of $M_{200}=10^{6}\,{\rm
M}_{\odot}$ at redshift $z=0$ have a mean peak velocity radius of $r_{\rm
max}\sim 100\,{\rm pc}$. Equating this radius to the substructure scale
radius, $c_{\bullet}=r_{\rm max}$, yields a considerably lower characteristic
density of $\rho_{\bullet}\sim 0.1\,{\rm M}_{\odot}\,{\rm pc}^{-3}$ (Ludlow et
al. 2016). On the one hand, the existence of DM sub-subhaloes with abnormally-
high densities may point to interactions in the dark sector (e.g. Kahlhoefer
et al. 2019). On the other, it is not immediately clear whether such a high
density would necessarily be in conflict with CDM predictions. First, because
sub-subhalos in satellite galaxies are expected to be tidally processed, and
more concentrated than field halos of the same mass (e.g. Peñarrubia et al.
2010; Errani & Navarro 2021). And second, because it is possible that the
detection of dense sub-subhaloes may be due to observational biases. Indeed,
compact substructures produce higher density enhancements in the field (see
§2), which can therefore be more easily detected as localized overdensities
($\delta_{\star}>1$). In this picture, anomalous objects like F6 may be
sampling the high-density, low-probability tail of the sub-subhalo population.
That is, they would represent the proverbial ‘tip of the iceberg’, pointing to
the presence of a much larger population of diffuse ($\delta_{\star}<1$)
substructures with a power-law luminosity function (15). There are several
ways to test this prediction, e.g. by modelling statistical fluctuations of
number counts in photometric surveys (e.g. Scheuer 1957), or by searching for
clumps with low-velocity dispersion and/or significant velocity offsets in
kinematic surveys of dSphs (e.g. Pace et al. 2014).

It is also worth noting that the presence of a large DM component in clusters
such as F6 and the Eri II would have important implications for their survival
in the host galaxy. For example, using collisional $N$-body simulations
Contenta et al. (2019) show that the Eri II cluster (with no DM in the star
cluster) quickly disrupt in a cuspy DM halo, favouring dSph models with a
cored DM profile. In addition, Brandt (2016) shows that the survival of a
stellar cluster near the centre of a dwarf galaxy depends on the number of
massive compact halo objects populating the DM halo (see also Li et al. 2017;
Zoutendijk et al. 2020), whereas Marsh & Niemeyer (2019) point out that the
cluster could also be used to test quantum fluctuations of ultra-light DM
models. However, if these objects are embedded in dense DM haloes, they would
become resilient to tidal stripping and would also be protected against
collisions with nearby compact objects and rapid fluctuations of the local
tidal field by dynamical invariance (e.g. Weinberg 1994; Peñarrubia 2019),
weakening the constraints derived from their very survival.

##  5 Summary

Substructures orbiting a larger galaxy can capture field stars that pass
nearby with low relative velocities. In this work, stars become captured by
the substructure potential $\Phi_{\bullet}$ when their gravitational energy
$E=v^{2}/2+\Phi_{\bullet}(r)$ flips from positive to negative.

Our numerical experiments show that the orbits of stars captured through this
mechanism can be separated into two families: (i) Temporary captures. The
majority of field stars trapped in the potential $\Phi_{\bullet}$ move on
chaotic orbits that orbit around the substructure for a finite amount of time
before being released back to the galactic potential. This leads to a net
number of temporary captures that converges to a steady state as the rate of
capture becomes roughly equal to the escape rate. Stars captured on temporary
orbits generate stellar overdensities ($\delta_{\star}>1)$ at the location of
substructures that are sufficiently compact,
$\kappa=1-c_{\bullet}/r_{\epsilon}>0$, where $r_{\epsilon}$ and $c_{\bullet}$
are the thermal critical radius (6) and scale radius of $\Phi_{\bullet}$,
respectively. (ii) Permanent captures. In addition, we find that immersing a
substructure potential in a galaxy that is in dynamical equilibrium at $t=0$
invariably results in a population of field particles that remain bound to
$\Phi_{\bullet}$ for arbitrarily-long times. In static models where
$\Phi_{\bullet}$ does not vary with time and moves on a circular orbit around
the host galaxy, we find that particles on permanent orbits were already bound
to the substructure potential at $t=0$, and remain within its tidal radius
indefinitely. These particles are therefore reminiscent of the ‘potential
escapers’ discussed by Fukushige & Heggie (2000). Our experiments also show
that substructures with a time-growing mass, $M_{\bullet}(t)$, can also
capture stars on permanent orbits at $t>0$, with a rate that peaks at
$t\sim\tau$, where $\tau=|\dot{M}_{\bullet}/M_{\bullet}|^{-1}_{t=0}$ is the
time-scale that determines how fast the substructure mass varies. Importantly,
we find that stars on permanent orbits have a density that lies below the host
galaxy background ($\delta_{\star}<1$), and therefore do not generate visible
overdensities.

The superposition of permanent and temporary orbits results in a population of
captured field stars with an homogeneous phase-space distribution around the
substructure potential, $Q_{\star}=n_{\star}/\sigma_{\star}^{3}\approx{\rm
const}$. In steady state, these stars follow number density and dispersion
profiles that scale as $n_{\star}\sim|\Phi_{\bullet}|^{3/2}$, and
$\sigma_{\star}\sim|\Phi_{\bullet}|^{1/2}$, respectively (see §2.1).
Experimental results show that these theoretical predictions are accurate for
substructures that source a static potential and move on a circular orbit
across a sea of field particles initially in equilibrium. They become less
accurate when applied to permanent & temporary captures individually, or to
substructures on very eccentric orbits moving through a rapidly-changing
background.

In the current cosmological paradigm, galaxies are expected to host a large
number of DM subhaloes, many of which may contain captured field stars. Here,
we show that these objects would have a luminosity function governed by the
subhalo mass function. Adopting a CDM-motivated mass function ${\rm d}N/{\rm
d}M_{\bullet}\sim M_{\bullet}^{-\alpha}$ with $\alpha\simeq 1.9$ (e.g.
Springel et al. 2008) leads to stellar substructures with a power-law
luminosity function (15), which scales as ${\rm d}N/{\rm d}M_{\star}\sim
M_{\star}^{-\beta}$, where $\beta=(2\alpha+1)/3\approx 1.6$. However, only
subhaloes that are sufficiently ‘compact’ ($\kappa>0$) may materialize as
localized stellar over-densities in the galactic background
($\delta_{\star}>1)$.

In dwarf spheroidal galaxies, field stars captured by compact substructures
may resemble stellar ‘clusters’ with anomalous properties, e.g. they have
extended sizes for their luminosity, contain stellar populations
indistinguishable from the field, and exhibit DM-dominated mass-to-light
ratios ($M/L\gg 1$). In contrast, field stars captured by ‘fluffly’
substructures ($\kappa<0$) have densities below the background
($\delta_{\star}<1)$, and are therefore more difficult to detect. In
spectroscopic surveys they may appear as clumps of co-moving field stars with
low velocity dispersion ($\langle\sigma_{\star}^{2}\rangle^{1/2}<\sigma$)
and/or a significant velocity offset ($V_{\bullet}\neq 0$). In cluster models
that contain no DM, the presence of stellar systems with extended sizes and
high mass-to-light ratios is typically attributed to on-going tidal
disruption, which can be tested by searching for the associated tidal tails
(Orkney et al. 2022). Yet, our experiments show that stars captured by
substructures moving on eccentric orbits are released back to the galactic
potential along tidal tails (see right panels Fig. 10), which complicates a
clear-cut distinction between the two scenarios.

Capture models make a number of unique testable predictions. E.g. stellar
overdensities composed of captured field stars have stellar ages and
metallicities indistinguishable from those of the host galaxy, and high mass-
to-light ratios $M/L\gg 1$ indicative of the presence of a substantial DM
component. In galaxies with multiple chemo-dynamical populations, like the
Fornax dSph, these models also make predictions on the ratio of captured stars
as a function of metallicity. In addition, the density and velocity dispersion
profiles of captured stars provide simultaneous constraints on the
substructure mass profile, ($M_{\bullet},c_{\bullet})$, as well as on its
systemic velocity with respect to the host galaxy, $V_{\bullet}$, which can be
potentially measured with a combination of accurate astrometric &
spectroscopic data. The numerical experiments outlined in §3 show that the
distribution of captured field stars also depends on the formation mechanism
of the substructure, e.g. whether it grew within a pre-existing sea of field
particles, or it formed outside of the galaxy and was accreted at a later time
(see Fig. 12). This suggests that it may be possible to constrain the
formation history, mass profile and systemic motion of the substructure
simultaneously by fitting the distribution and kinematics of stellar systems
with anomalous properties.

We have given numerical and theoretical arguments that indicate that ‘dark’
subhaloes that were not sufficiently massive to trigger in-situ star formation
may become ‘visible’ by capturing baryonic particles from the host galaxy
field. The important implication of this result is that dark subhaloes cannot
be completely invisible. If they contain gravitationally-bound baryonic
matter, they must emit and absorb radiation, which opens up a new avenue to
test CDM predictions on halo mass scales that have not been probed to date.
This calls for observational campaigns aimed at detecting & characterizing
‘dark’ subhaloes in dSphs with masses below the star formation threshold. In
§4.3 we discuss how the detection and characterization of these objects may
provide unique constraints on the particle mass and cross section for a large
range of DM candidates.

## Acknowledgements

We thank Sergey Koposov, Alan McConnachie and Frank van den Bosch for their
valuable comments and suggestions. M.G.W. acknowledges support from National
Science Foundation (NSF) grants AST-1813881, AST-1909584 and AST-2206046. M.G.
acknowledges financial support from the grants PID2021-125485NB-C22,
EUR2020-112157, CEX2019-000918-M funded by MCIN/AEI/10.13039/501100011033
(State Agency for Research of the Spanish Ministry of Science and Innovation)
and SGR-2021-01069 (AGAUR). RE acknowledges support from the European Research
Council (ERC) under the European Unions Horizon 2020 research and innovation
programme (grant agreement number 834148), and from the National Science
Foundation (NSF) grant AST-2206046. TB’s research was supported by an
appointment to the NASA Postdoctoral Program at the NASA Ames Research Center,
administered by Oak Ridge Associated Universities under contract with NASA.

## Data Availability

No data were generated for this study.

## References

  * Ackermann et al. (2014) Ackermann, M., Albert, A., Anderson, B., et al. 2014, Phys. Rev. D, 89, 042001
  * Amorisco & Evans (2012) Amorisco N. C., Evans N. W., 2012, ApJL, 756, L2. doi:10.1088/2041-8205/756/1/L2
  * Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781. doi:10.1111/j.1365-2966.2006.10094.x
  * Bate et al. (2015) Bate N. F., McMonigal B., Lewis G. F., Irwin M. J., Gonzalez-Solares E., Shanks T., Metcalfe N., 2015, MNRAS, 453, 690. doi:10.1093/mnras/stv1684
  * Battaglia et al. (2006) Battaglia G., Tolstoy E., Helmi A., Irwin M. J., Letarte B., Jablonka P., Hill V., et al., 2006, A&A, 459, 423. doi:10.1051/0004-6361:20065720
  * Benitez-Llambay & Frenk (2020) Benitez-Llambay A., Frenk C., 2020, MNRAS, 498, 4887. doi:10.1093/mnras/staa2698
  * Benson (2017) Benson, A. J. 2017, MNRAS, 467, 3454
  * Boekholt, Rowan, & Kocsis (2023) Boekholt T. C. N., Rowan C., Kocsis B., 2023, MNRAS, 518, 5653. doi:10.1093/mnras/stac3495
  * Bovill & Ricotti (2009) Bovill M. S., Ricotti M., 2009, ApJ, 693, 1859. doi:10.1088/0004-637X/693/2/1859
  * Brandt (2016) Brandt T. D., 2016, ApJL, 824, L31. doi:10.3847/2041-8205/824/2/L31
  * Bringmann et al. (2014) Bringmann, T., Calore, F., Di Mauro, M., & Donato, F. 2014, Phys. Rev. D, 89, 023012
  * Bullock et al. (2000) Bullock J. S., Kravtsov A. V., Weinberg D. H., 2000, ApJ, 539, 517
  * Carlberg (2013) Carlberg, R. G. 2013, ApJ, 775, 90
  * Coleman et al. (2004) Coleman M., Da Costa G. S., Martínez-Delgado D., Bland-Hawthorn J., 2004, ASPC, 327, 173
  * Contenta et al. (2018) Contenta F., Balbinot E., Petts J. A., Read J. I., Gieles M., Collins M. L. M., Peñarrubia J., et al., 2018, MNRAS, 476, 3124. doi:10.1093/mnras/sty424
  * Crnojević et al. (2016) Crnojević D., Sand D. J., Zaritsky D., Spekkens K., Willman B., Hargis J. R., 2016, ApJL, 824, L14. doi:10.3847/2041-8205/824/1/L14
  * Daniel, Heggie, & Varri (2017) Daniel K. J., Heggie D. C., Varri A. L., 2017, MNRAS, 468, 1453. doi:10.1093/mnras/stx571
  * de Boer et al. (2013) de Boer T. J. L., Tolstoy E., Saha A., Olszewski E. W., 2013, A&A, 551, A103. doi:10.1051/0004-6361/201220855
  * Deason, Belokurov, & Evans (2011) Deason A. J., Belokurov V., Evans N. W., 2011, MNRAS, 416, 2903. doi:10.1111/j.1365-2966.2011.19237.x
  * Despali, Giocoli, & Tormen (2014) Despali G., Giocoli C., Tormen G., 2014, MNRAS, 443, 3208. doi:10.1093/mnras/stu1393
  * Diemand et al. (2005) Diemand, J., Moore, B., & Stadel, J. 2005, Nature, 433, 389
  * Diemer & Joyce (2019) Diemer B., Joyce M., 2019, ApJ, 871, 168. doi:10.3847/1538-4357/aafad6