# Capture of field stars by dark substructures

Jorge Peñarrubia, Raphaël Errani, Matthew G Walker, Mark Gieles, Tjarda C N Boekholt

► To cite this version:

Jorge Peñarrubia, Raphaël Errani, Matthew G Walker, Mark Gieles, Tjarda C N Boekholt. Capture of field stars by dark substructures. Monthly Notices of the Royal Astronomical Society, 2024, 533 (3), pp.3263-3282. 10.1093/mnras/stae1961. insu-04702443

HAL Id: insu-04702443

https://insu.hal.science/insu-04702443v1

Submitted on 19 Sep 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Capture of field stars by dark substructures

Jorge Peñarrubia,<fnref n="1" />,<fnref n="2" />* Raphaël Errani,<fnref n="3" />,<fnref n="4" /> Matthew G. Walker,<fnref n="5" />,<fnref n="6" /> and Tjarda C. N. Boekholt,<fnref n="7" />,<fnref n="8" />

Accepted 2024 August 8. Received 2024 July 22; in original form 2023 August 4

## ABSTRACT

We use analytical and N-body methods to study the capture of field stars by gravitating substructures moving across a galactic environment. The majority of stars captured by a substructure move on temporarily bound orbits that are lost to galactic tides after a few orbital revolutions. In numerical experiments where a substructure model is immersed into a sea of field particles on a circular orbit, we find a population of 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 that exhibit some of these characteristics. DM sub-subhaloes with a mass function $dN/dM_\bullet \sim M_\bullet^{-\alpha}$ are expected to generate stellar systems with a luminosity function, $dN/dM_\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.

**Key words:** 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 M_\odot$ (White & Rees 1978; Bullock, Kravtsov & Weinberg 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$ GeV/c<fnref n="2" /> is at least $\sim 13$ orders of magnitude smaller, $< 10^{-6} M_\odot$ (e.g. Schmid, Schwarz & Widerin 1999; Hofmann, Schwarz & Stöcker 2001; Green, Hofmann & Schwarz 2005; Loeb & Zaldarriaga 2005; Diemand, Moore & Stadel 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, Johnston & Hogg 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 M_\odot$ and a scale radius $c_\bullet \lesssim 10$ 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 M_\odot \text{pc}^{-3}$, which is several orders of magnitude denser than CDM subhaloes with similar masses found in

<footnote n="1">¹Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="2">²Centre for Statistics, University of Edinburgh, School of Mathematics, Edinburgh EH9 3FD, UK</footnote>

<footnote n="3">³McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, 5000 Forbes Ave., Pittsburgh, PA 15213, USA</footnote>

<footnote n="4">⁴Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="5">⁵ICREA, Pg. Lluís Companys 23, E-08010 Barcelona, Spain</footnote>

<footnote n="6">⁶Institut de Ciències del Cosmos (ICCB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, E-08028 Barcelona, Spain</footnote>

<footnote n="7">⁷Rudolf Peierls Centre for Theoretical Physics, Clarendon Laboratory, University of Oxford, Parks Road, Oxford OX1 3PU, UK</footnote>

<footnote n="8">⁸NASA Ames Research Center, Moffett Field, CA 94035, USA</footnote>

<footnote>*E-mail: jorpega@roe.ac.uk</footnote>
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 capture<fnref n="1" /> 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 dom-inated. 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 classical<fnref n="2" /> restricted three-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 (1977) 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, three-body captures in accretion discs have also gained attention as a possible source of black-hole binaries. E.g. Li, Lai & Rodet (2022), Boekholt, Rowan & Kocsis (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 three-body capture events is extremely complex. Using a combination of numerical and analytical methods, Petit & Hénon (1986) showed that (i) captures in a three-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 ≃ 0.4.

Peñarrubia (2023, hereafter Paper I) studies three-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 sling-shot 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 particles,<fnref n="3" /> 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 three-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, Heggie & Varri 2017).

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. Section 2 extends the statistical method outlined in Paper I to model gravitational captures by point-mass objects to substructures an extended size. Section 3 presents numerical experiments that test the accuracy of the theoretical equations. Section 4 discusses these results in the context of dark matter particle physics and outlines future follow-up work. Finally, in Section 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_*(r) < 0$ within a spherical volume $V = 4\pi r^3/3$, where $r$ and $v$ are measured relative to the substructure, as well as their number density ($n_*$) and velocity dispersion ($\sigma_*$) profiles in dynamical equilibrium within a generic potential $\Phi_*$ . 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 particles.<fnref n="4" />

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

$$ \Phi_*(r) = -\frac{GM_*}{r+c_*} \qquad (1) $$

which recovers the results derived in Paper I for point-masses if the profile scale-length is set to $c_* = 0$. Note that the density profile associated with the potential (1) has a centrally divergent cusp, $\rho \sim r^{-1}$ at $r \ll c_*$ , 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

<footnote n="1">¹In 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.</footnote>

<footnote n="2">²Two point masses on circular orbits, plus a mass-less tracer</footnote>

<footnote n="3">³A 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.</footnote>

<footnote n="4">⁴To model the dynamics of gaseous particles, a pressure term must be included in the equations of motion presented in Section 3.1.</footnote>
distances, appear to be the most interesting objects to test the existence of dark substructures via captured field stars.

Section 2.3 analyses 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(R_{\bullet} + r) \approx n(R_{\bullet}) = n = N/V$, where $R_{\bullet}$ is the 3D position vector of the substructure within the host galaxy, and $r = |R - R_{\bullet}|$ is the relative distance of a star from the centre of the potential $\Phi_{\bullet}$. 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(v) = (2\pi\sigma^2)^{-3/2} \exp[-(v + V_{\bullet})^2/(2\sigma^2)]$ where $\sigma = \sigma(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 $V_{\bullet}$ is the velocity vector of the substructure with respect to the host galaxy centre.

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

$$ 
\begin{aligned}
N_b(r) &= \int_V d^3 r n(r) \int_{E<0} d^3 v p(v) \\
&= \frac{1}{3} \sqrt{\frac{2}{\pi}} z e^{-V_{\bullet}^2/(2\sigma^2)} \frac{n}{\sigma^3} \int_V d^3 r v_e^3(r) + O(v_e/\sigma)^5 \\
&= \frac{32\sqrt{\pi}}{9} (GM_{\bullet})^{3/2} e^{-V_{\bullet}^2/(2\sigma^2)} \frac{n}{\sigma^3} \\
&\quad \times \frac{r^2 - 4c_{\bullet}r + 8c_{\bullet}^{3/2}(r+c_{\bullet})^{1/2} - 8c_{\bullet}^2}{(r+c_{\bullet})^{1/2}},
\end{aligned} 
\quad (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$).

The approximation that field stars move on uncorrelated trajectories is accurate insofar as the number of stars perturbed by $\Phi_{\bullet}$ represents a small fraction of the total number of stars within the volume element, $N = nV$. From equation (2), it is straightforward to show that the fraction of bound stars within the volume $V$ scales as $N_b/N \sim (v_e/\sigma)^3 e^{-V_{\bullet}^2/(2\sigma^2)}$ at $r \gtrsim c_{\bullet}$, hence field particles can be treated independently from each other when the escape velocity is low, $v_e \lesssim \sigma$, or the substructure velocity is high, $V_{\bullet} \gg \sigma$.

The orbits of field stars captured by the potential $\Phi_{\bullet}$ exhibit chaotic fluctuations of energy and angular momentum, as shown in Section 3.2. Statistically, one 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). Petit & Hénon (1986) found that captured particles moving on 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'

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

which roughly corresponds to the time that a particle moving on a straightline 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_{ss} = \alpha N_b$, where $\alpha$ is the so-called abundance parameter, which is set by the 'dynamical survival time'<fnref n="5" /> 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$.

### 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 close to the substructure, where the potential well $|\Phi_{\bullet}|$ may become 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 Section 3): (i) the phase-space density of captured stars becomes approximately constant in the vicinity of the substructure, such that $f_{\bullet}(\mathbf{r}, \mathbf{v}) \approx f_0$, 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, we do not 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. Under this condition,

$$ 
\begin{aligned}
n_{\bullet}[\Phi_{\bullet}(r)] &= \int_{E<0} d^3 v f_{\bullet}(\mathbf{r}, \mathbf{v}) \\
&= 4\pi f_0 \int_0^{|\Phi_{\bullet}|} dE [2(E - \Phi_{\bullet})]^{1/2} \\
&= \frac{8\sqrt{2\pi}}{3} f_0 |\Phi_{\bullet}|^{3/2}.
\end{aligned} 
\quad (4) $$

To find the normalization of the local phase-space density, $f_0$, we match the number density (4) against the steady-state profile derived from equation (2) at low escape velocities ($v_e \lesssim \sigma$), such that $n_{\bullet} =$

<footnote n="5">⁵Defined as the time over which a trapped particle continuously has $E < 0$ before being lost to galactic tides.</footnote>
d N<sub>ss</sub>/d<sup>3</sup>r = α d N<sub>b</sub>/d<sup>3</sup>r. This yields

$$
f_0 = \alpha \frac{n}{(2\pi\sigma^2)^{3/2}} e^{-V_{\bullet}^2/(2\sigma^2)}. \qquad (5)
$$

which corresponds to a local Maxwellian distribution function multiplied by abundance parameter α. The density enhancement induced by the population of energetically bound stars within the volume V can be found by inserting (5) into (4), which yields

$$
\begin{align}
\delta_{\star}(r) &\equiv \frac{n_{\star}(r)}{n} \nonumber \\
&= \alpha \frac{4}{3\sqrt{\pi}} e^{-V_{\star}^2/(2\sigma^2)} \frac{1}{\sigma^3} |\Phi_{\star}|^{3/2} \nonumber \\
&= \alpha \frac{4}{3\sqrt{\pi}} \frac{(GM_{\star})^{3/2}}{\sigma^3} e^{-V_{\star}^2/(2\sigma^2)} \frac{1}{(r+c_{\star})^{3/2}}. \tag{6}
\end{align}
$$

Notice that the density profile 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 haloes ($c_{\bullet} > 0$) the profile becomes
flat for $r \lesssim c_{\bullet}$. For generic potentials, we find that the density en-
hancement (6) scales as $\delta_{\star} \sim |\Phi_{\star}|^{3/2}$, which allows a straightforward
analysis of a wide range of substructure models.

The isotropic 1D velocity dispersion can be derived from the local
distribution function as

$$
\begin{align}
\sigma_{\star}^{2}[\Phi_{\bullet}(r)] &= \frac{1}{3n_{\star}(r)} \int_{E<0} d^{3}v v^{2} f_{*}(\mathbf{r}, \mathbf{v}) \nonumber \\
&= \frac{4\pi f_{0}}{3n_{\star}(r)} \int_{0}^{\Phi_{\bullet}} dE [2(E - \Phi_{\bullet})]^{3/2} \nonumber \\
&= \frac{2}{5} |\Phi_{\bullet}|, \tag{7}
\end{align}
$$

which is a constant fraction of the escape velocity at all radii, $\sigma_*/v_e = (1/5)^{1/2} \approx 0.45$. It is important to emphasize that equation (7) is independent of the speed of the dark object across the host galaxy ($V_\bullet$) as well as of the normalization of the distribution function ($f_0$).

By construction, combining equations (6) and (7) returns a mean
phase-space density of bound stars that remains constant across the
volume element

$$
Q_{\star} \equiv \frac{n_{\star}(r)}{\sigma_{\star}^{3}(r)} = \alpha \frac{5}{3} \sqrt{\frac{10}{\pi}} e^{-V_{\star}^{2}/(3\sigma^{2})} Q. \quad (8)
$$

equation (8) highlights two interesting results. First, Q<sub>*</sub> solely
depends on the potential Φ<sub>•</sub> through the value of the abundance
parameter, α, and second the phase-space density of captured stars
drops exponentially for a potential Φ<sub>•</sub> that is not at rest with the
background (V<sub>•</sub> > 0). In Section 3.1, we carry a number of numerical
experiments that inspect the accuracy of the theoretical profiles
derived above.

Of particular relevance for this paper is the distance at which the
density of bounds stars is equal to that of the field. For reasons that
become clear below, it is convenient to examine the point-mass limit
first. Solving δ(r<sub>ε</sub>) = 1 in equation (6) with c<sub>•</sub> → 0 and α = 1 yields

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

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 (9). 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_*(r_\epsilon^{\text{ext}}) = 1$ happens at $r_\epsilon^{\text{ext}} = r_\epsilon - c_\bullet = \kappa r_\epsilon$, where $\kappa \equiv 1 - c_\bullet/r_\epsilon$ is the com-

pactness 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
overdense 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 overdensity-
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_{\text{min}}$. Using (6) with $c_\bullet = 0$
and (9) and setting $\alpha = 1$, we find

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

here D = (2π n)<sup>-1/3</sup> is a measure of the average separation of stars
in the volume element V (Peñarrubia 2018). Notice again that the
formation of stellar overdensities is strongly suppressed in dark
objects moving at high speed (V<sub>•</sub> ≫ σ) with respect to the galactic
background.

## 2.2 Application: Milky Way dwarf spheroidals

The above analytical expressions can be used to inspect the effi-
ciency 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/σ³ =
(2π D³ σ³)⁻¹ 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
ultrafaint 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 5.

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

Dwarf galaxies with multiple chemo-dynamical components are
particularly interesting in this context. In these systems, the proba-
bility 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 [Fe/H] = -1.8
and -0.65 and ages > 10 Gyr, and ~ 0.2-2 Gyr, respectively
(Rusakov et al. 2021). The metal-rich (MR) component contains
more stars $N_*^{\text{MR}} \sim 3 \times 10^7$, is more centrally concentrated, $r_h^{\text{MR}} \approx$
530 pc and has a colder velocity dispersion, $\sigma^{\text{MR}} \approx 10 \text{ km s}^{-1}$,
than the metal-poor (MP) population, which has $N_*^{\text{MP}} \sim 2 \times 10^7$,
$r_h^{\text{MP}} \approx 1070$ pc and $\sigma^{\text{MP}} \approx 14.4 \text{ km s}^{-1}$ (Walker & Peñarrubia 2011).
Thus, MR stars have much higher phase-space densities than MP
stars,

$$
\frac{Q^{\text{MR}}}{Q^{\text{MP}}} = \frac{n^{\text{MR}}}{n^{\text{MP}}} \left(\frac{\sigma^{\text{MP}}}{\sigma^{\text{MR}}}\right)^3 = \frac{N_{*}^{\text{MR}}}{N_{*}^{\text{MP}}} \left(\frac{r_{h}^{\text{MP}}}{r_{h}^{\text{MR}}} \frac{\sigma^{\text{MP}}}{\sigma^{\text{MR}}}\right)^3 \sim 33. \quad (11)
$$

Middle panel of Fig. 1 shows estimates of $M_{\text{min}}$ derived from
equation (10) for the MW satellites in our sample. To simplify the
analysis, substructures are placed at rest within the host galaxies
<figure><img src="image_2.png" /><figcaption>**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 5. Notice that the phase-space density of the MW stellar halo at $R = 20$ kpc is $\gtrsim 6$ orders of magnitude lower than in dSphs. This difference grows at larger Galactocentric radii. Middle panel: Minimum substructure mass ($M_{\min}$) derived from equation (10). Substructures with $M_{\bullet} > M_{\min}$ are expected to contain stellar overdensities ($\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_{\min}$ in the MW halo is relatively large on account of its high velocity dispersion. Lower panel: Minimum thermal critical radius $r_{\min}^{\text{dSph}}$ of substructures with a mass $M_{\bullet} = M_{\min}$. Stellar overdensities made of capture stars have a physical size comparable to the radii of globular clusters.</figcaption></figure>

($V_{\bullet} = 0$). For reference, we also plot the value of $M_{\min}$ associated with the stellar halo of the Milky Way at a Galactocentric distance of $R = 20$ kpc using the stellar density and velocity dispersion measured by Deason, Belokurov & Evans (2011). This panel reveals a number of interesting points. Notice first that the mass threshold ($M_{\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 ultrafaint dSphs with half-light radii $r_{\min}^{\text{dSph}} \sim 30$ pc are the best targets, as they may be sensitive to the presence of dark substructures with masses above $M_{\bullet} \gtrsim 10^4 M_{\odot}$. In contrast, the mass threshold increases by an order of magnitude, $M_{\bullet} \gtrsim 10^5 M_{\odot}$, in the ‘classical’ dwarf galaxies with $r_{\min}^{\text{dSph}} \gtrsim 300$ pc. In the case of the Fornax dSph, the MR stellar component a mass threshold $M_{\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$ 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 Section 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 overdensity 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 (9) returns a dimensionless mass-to-light ratio

$$ \frac{M}{L} \equiv \frac{M(<r_{\epsilon})}{N_b(r_{\epsilon}) \langle m_{\star} \rangle} = \frac{M_{\min}^3}{M_{\bullet}^2 \langle m_{\star} \rangle}, \quad (12) $$

where $M_{\min}$ is given by equation (10). According to the estimates plotted in Fig. 1, field stars bound to substructures with a mass $M_{\bullet} \sim M_{\min}$ will appear as extremely DM-dominated objects, $M/L = M_{\min}/\langle m_{\star} \rangle \sim 10^4–10^5$ for an average stellar mass $\langle m_{\star} \rangle = 1 M_{\odot}$. In substructures with larger masses, $M_{\bullet} > M_{\min}$, the mass-to-light ratio of captured stars drops as $M/L \sim M_{\bullet}^{-2}$. Note that equation (12) should not be applied to substructures, where $M/L \lesssim 1$, as the theory outlined in Section 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^{\text{MW}} \sim 124 \text{ km s}^{-1}$ and $n^{\text{MW}} = 10^{-4} \text{ pc}^{-3}$ at $R = 20$ 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_{\min} \gtrsim 10^8 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} M_{\odot}$ (e.g. Peñarrubia et al. 2008a; Errani, Peñarrubia & Walker 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$ pc (e.g. Diemar & Joyce 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$). 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.

## 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 $dN/dM_{\bullet} \sim M_{\bullet}^{-\alpha}$ and a slope $\alpha \approx 1.9$ (e.g. Giocoli, Tormen & van den Bosch 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_{\min}$ is expected to arise on scales comparable to the particle free-streaming length. For WIMPS with masses above $\sim 1 \text{GeV/c}^2$ the truncation lies below planet-mass scales, $M_{\min} \lesssim 10^{-6} 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_{\text{out}})$ within a volume size larger than the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left-hand panel: mass $M_*$ and scale radius $c_*$ of DM sub-subhaloes colour-coded according to the number of captured dwarf galaxy stars, $N_* = N_b(r_e)$, equation (2), where $r_e$ is the thermal critical radius, equation (9). In this model, the host galaxy has properties similar to the Fornax dSph: it follows a Plummer profile with a total luminosity of $N^{\text{dSph}} = 5 \times 10^7$, a half-light radius $r_h^{\text{dSph}} = 0.89$ kpc and it has an average velocity dispersion $\sigma^{\text{dSph}} = 12.2$ 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_*/r_e$ of the left-and panel models. Only models with $\kappa > 0$ lead to overdensities of captured stars with $\delta_* > 1$. Right-hand panel: Mass-to-light ratios of the left-hand and middle-panel models derived from equation (12). Stellar systems trapped in subhaloes with masses $M_* \lesssim 10^7$ M$_\odot$ will exhibit dark-matter dominated kinematics.</figcaption></figure>

substructure scale radius, $r_{\text{out}} \gg c_*$. From equation (2)

$$N_* = \frac{32\sqrt{\pi}}{9} (GM_*)^{3/2} e^{-V_*^2/(2\sigma^2)} \frac{n}{\sigma^3} r_{\text{out}}^{3/2}, \quad \text{for } r_{\text{out}} \gg c_* \quad (13)$$

To combine equation (13) with the subhalo mass function, we apply the chain rule $(dN/dM_*)dM_* \sim N_*^{-(2\alpha+1)/3}dN_*$, which leads to a luminosity function

$$\frac{dN}{dM_*} = A_0 \left( \frac{M_*}{\langle m_* \rangle} \right)^{-(2\alpha+1)/3}, \quad (14)$$

where $M_* = N_*\langle m_* \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 $dN/dM_* \sim M_*^{-\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 these objects can be detected as localized stellar overdensities mainly depends on the *compactness* of DM haloes that capture them (see Section 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 Section 3.1. Namely, field stars follow a Plummer profile with a total luminosity of $N^{\text{dSph}} = 5 \times 10^7$, a half-light radius of $r_h^{\text{dSph}} = 0.89$ kpc and a mean velocity dispersion $\sigma^{\text{dSph}} = 12.2$ km s$^{-1}$, similar to the overall properties of the Fornax dSph (see Section 2.2). For simplicity, we assume that dark substructures are at rest with the field ($V_* = 0$).

Left-hand panel of Fig. 2 shows models colour-coded according to the number of bound stars enclosed within their thermal critical radius, $N_b(r_e)$ calculated from equations (2) and (9). As shown in Fig. 1, point-mass objects with $M_* \gtrsim 10^5$ M$_\odot$ contain $N_* \gtrsim 1$ stars trapped within their thermal critical radius. In extended substructures, the mass threshold for capture increases mildly.

Only substructures with positive *compactness* ($\kappa \equiv 1 - c_*/r_e$) can become visible as localized stellar overdensities ($\delta_* > 1$). Middle panel shows that the condition $\kappa > 0$ (or $c_* < r_e$) translates into a linear relation $c_* \lesssim 0.83 GM_*/\sigma^2$, which is steeper than the shallow mass-scale radius relation of Aquarius subhaloes, $c_* \sim M_*^{0.46}$. As a result, Aquarius subhaloes with masses $M_* \lesssim 5 \times 10^6$ M$_\odot$ would be ‘fluffy’ ($\kappa < 0$). In these objects, trapped stars would have sub-dominant densities with respect to the field ($\delta_* < 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_* < \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 overdensities 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 CDM 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 overdensities of captured field stars. We will discuss this issue in Section 4.1.

Right-hand panel of Fig. 2 shows the mass-to-light ratio (12) for the models shown in previous panels. As expected, we find that field stars trapped in low-mass ($M_* \lesssim 10^7$ M$_\odot$) substructures exhibit DM-dominated kinematics ($M/L \gg 1$). As a note of caution, recall that the analytical equations derived in Section 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 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 Section 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:

Section 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.

Section 3.2 injects a static substructure potential $\Phi_{\bullet}$ in a pre-existing field of particles in dynamical equilibrium within $\Phi^{\text{dSph}}$. This set-up is similar to the analytical conditions adopted in Section 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 indefinitely long times, and (ii) temporary captures, which as the name indicates only remain bound for a finite amount of time.

Section 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 |\Phi_{\bullet}/\Phi_{\bullet}|^{-1}$ the phase-space distribution of particles with $E<0$ become indistinguishable from the static case explored in Section 3.2.

In Section 3.4 we analyse 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 innermost regions of the dSph.

Section 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^{\text{dSph}}$ with a total mass of $M^{\text{dSph}} = 3 \times 10^9 M_{\odot}$ and a scale radius $c^{\text{dSph}} = 2$ 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 L_{\odot}$ (e.g. Rusakov et al. 2021), the particle luminosity in our models is $10 L_{\odot}$.

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

$$n(R) = \frac{n_0}{(R/R_0)^{\gamma_f} [1 + (R/R_0)^{\alpha_f}]^{(\beta_f - \gamma_f)/\alpha_f}}, \quad (15)$$

with $n_0$ chosen such that $4\pi \int_0^\infty dr r^2 n(r) = N$. We consider two tracer models: (i) a spherical Plummer-like (1911) profile<fnref n="6" /> with a scale radius $R_0 = 690$ 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$ 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 \text{ 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 \text{ 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

$$\ddot{R}_{\bullet} = -\nabla\Phi^{\text{dSph}}(R_{\bullet}). \quad (16)$$

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

$$\ddot{R} = -\nabla\Phi^{\text{dSph}}(R) - \nabla\Phi_{\bullet}(R - R_{\bullet}) + F_{\text{coll}}, \quad (17)$$

where $R - R_{\bullet}$ is the relative distance of the particle to the substructure, and $F_{\text{coll}} = \sum_{i=1}^{N_{\text{clump}}} f_i$ is the net force induced by $N_{\text{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_{\text{coll}} = 0$. We will analyse 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$ kpc with a tangential velocity $V_{\bullet} = V_c(R_{\bullet}) = 14.35 \text{ km s}^{-1}$. The local dynamical time is therefore $\Omega^{-1} = R_{\bullet}/V_{\bullet} = 500 \text{ pc}/14.35 \text{ km s}^{-1} \simeq 34 \text{ Myr}$. Field stars follow a Plummer (1911) profile, have a local velocity dispersion of $\sigma(R_{\bullet}) = 11.2 \text{ km s}^{-1}$, and are separated by an average distance $D(R_{\bullet}) = 4.8 \text{ pc}$.

For illustration, we choose a substructure mass $M_{\bullet} = 10^6 M_{\odot}$. To estimate the thermal critical radius we adopt $\sigma = \langle v^2 \rangle^{1/2}/3 = 12.2 \text{ km s}^{-1}$ in equation (9), which returns $r_{\epsilon} = 18.9 \text{ 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 \text{ pc}$ (notice that this value is smaller than the typical scale radius of field haloes with a vital mass $M_{\text{vir}}(z=0) \sim 10^6 M_{\odot}$ in the Aquarius simulation, see Fig. 2). The time-scale needed to cross the overdensity size by a random field particle is much shorter than the orbital time, $T(r_{\epsilon}) = r_{\epsilon}/\sigma = 18.9 \text{ pc}/(12.2 \text{ km s}^{-1}) \simeq 1.5 \text{ 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 (20), $r_{\text{out}} = r_t \simeq 190 \text{ pc}$. The volume $V = 4\pi r_{\text{out}}^3/3$ contains approximately $N = nV = 42530$ field stars. From equation (2), we expect that approximately $N_b(r_{\text{out}}) = 2000$ stars will have negative energies.

<footnote n="6">⁶We 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.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Distribution of capture times ($t_c = t - t_{\text{bound}}$) defined as the time spanned since a particle becomes bound ($t_{\text{bound}}$), i.e. the energy *E* flips from positive to negative, until the present integration time (*t*). Time is measured in units Smoluchowski's time scale (*T*) for convenience. Here, we only include particles within a volume size equal to the local tidal radius, $r_{\text{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 strong peak in the distribution at $t_c \approx t$ shifting towards the right as the integration time increases.</figcaption></figure>

## 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 Section 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_{\text{out}}^3/3$ centred at $\Phi_{\bullet}$ with a relative velocity below the escape velocity, $v < v_{\text{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 & Hénon 1986), which take 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 (16) and (17). 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_{\text{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_{\text{bound}} = t$. This spike is populated by particles that were

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Examples of permanent and temporary captures. Left-hand 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 M_{\odot}$ and a compactness $\kappa = 0.2$ (see Section 2.1 for details), which leads to a scale radius $c_{\bullet} = (1-\kappa)r_{\epsilon} = 3.78$ 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 $|L| = |\mathbf{r} \times \mathbf{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.</figcaption></figure>

already bound at $t=0$ ($t_{\text{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_{\text{bound}} > 0$. Their distribution of capture times peaks at $t_c \sim 2T$, although in this experiment we also find a significant number of temporary captures that remain bound for as long as $t_c \sim 50T$. In our models, all particles captured with $t_{\text{bound}} > 0$ are temporary events that ultimately lead to the tracer particle being released back to the galactic potential.

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 \approx 2T$ and remains negative until $t \approx 47T$, which yields a survival time $t_{\text{surv}} \sim 45T$. After $t \gtrsim 47T$, 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_{\text{surv}} \sim T$, until it finally becomes fully unbound and escapes from the volume element $V$ at $t \sim 60T$. 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 ($L = \mathbf{r} \times \mathbf{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
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Upper panel: Radial phase-space location of bound particles in a co-rotating frame, $x' - v_x'$ measured at the end of the simulation ($t = t_f$). Only particles within a narrow plane $|\Delta y'| < 0.2 r_e$ and $|\Delta z'| < 0.2 r_e$ are shown here. The dwarf galaxy centre is located at $x'/r_e \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_{bound} = 0$) and plot them with red dots, whereas temporary captures become bound at $t_{bound} > 0$ and are plotted in blue. Note that permanent captures have negative binding energies measured in a co-rotating frame, $E_{rot} < 0$, and are confined within the tidal radius $|x'| < r_t \approx 10 r_e$, with $r_t$ given by equation (20). Notice also the lack of permanent captures at very small distances from the substructure, $r \ll r_e$. Lower panel: Energy and angular momentum of the particles plotted above given in the same units as in the right-hand 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$.</figcaption></figure>

(such as peri, apo-centres and orbital eccentricity) that change in an apparently random fashion over short time scales.

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'$-axis pointing at the host galaxy centre, and the $z'$-axis aligned with the angular vector $\Omega = (0, 0, \Omega)$. In this frame, the effective potential can be written as (cf. equation 12

of Renaud, Gieles & Boily 2011)

$$ \Phi_{\text{eff}} = \Phi_{\bullet} - \frac{1}{2}(\lambda_1 x'^2 + \lambda_2 y'^2 + \lambda_3 z'^2), \qquad (18) $$

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 equation (17) of Renaud et al. (2011) as

$$ (\lambda_1, \lambda_2, \lambda_3) = \frac{GM^{\text{dSph}}}{(R_{\bullet} + c^{\text{dSph}})^3} \left( \frac{3R_{\bullet}}{R_{\bullet} + c^{\text{dSph}}}, 0, -1 \right). \qquad (19) $$

The tidal radius can be estimated from (19) as

$$ r_t = \left( \frac{GM_{\bullet}}{\lambda_1} \right)^{1/3} = \left[ \frac{M_{\bullet}}{3M^{\text{dSph}}} \left( 1 + \frac{c^{\text{dSph}}}{R_{\bullet}} \right) \right]^{1/3} R_{\bullet}. \qquad (20) $$

which recovers the academic point-mass case for $c^{\text{dSph}} = 0$. Inserting the parameters of the experiments run in this section into (20) yields $r_t \approx 190 \text{ pc} \approx 10 r_e$.

Dotted and solid lines in the upper panel of Fig. 5 show phase-space surfaces that obey the conditions $E = v_x'^2/2 + \Phi_{\bullet}(x', 0, 0) = 0$ and $E_{rot} = v_x'^2/2 + \Phi_{\text{eff}}(x', 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_{rot} < 0$ and can only be found at distances $r < r_t \approx 10 r_e$ from the substructure. Notice also the apparent lack of permanent captures at distances smaller than the thermal critical radius, $|x| \lesssim r_e$. 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 = |\mathbf{r} \times \mathbf{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.

Fig. 6 shows the density enhancement of bound particles, $\delta_* = n_*/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_* = 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 overdensities, i.e. $\delta_* < 1$. We also observe a sharp truncation of the profile at $r \approx 10 r_e$, which roughly corresponds
<figure><img src="image_7.png" /><figcaption>**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 permanent captures beyond the tidal radius $r_{t}$, equation (20), 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.</figcaption></figure>

to the value estimated from (20) (marked with vertical arrows). At small distances, $r \lesssim r_{\epsilon}$, we find no permanent captures. This is likely a numerical artefact 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 overdensity 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 overdensity ($\delta_{\star} > 1$) at small distances from the substructure. We also find that the enhancement profile closely follows the theoretical curve (6) (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 (6).

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 Section 2.1.2, the theoretical profile (6) 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 (6) 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.

### 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 Section 4). To this aim, we re-run the ‘compact’ N-body models presented in Section 3.2 adopting a time-dependent substructure mass

$$M_{\bullet}(t) = M_{\bullet} \left[1 - \exp(-t/\tau)\right], \qquad (21)$$

with $M_{\bullet} = 10^{6} M_{\odot}$ and a fixed scale radius $c_{\bullet} = 3.78$ pc. Here, $\tau$ is a parameter that controls the mass growth rate, $\tau = |\dot{M}_{\bullet}/M_{\bullet}|_{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 Section 3.2.

Fig. 7 shows the distribution of capture times, $t_{c} = t - t_{\text{bound}}$, where $t_{\text{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
<figure><img src="image_8.png" /><figcaption>**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_{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_{out} = 190$ 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 (21). At time-scales $t \gg \tau$, the three models show the same fraction of permanent particles independently of $\tau$.</figcaption></figure>

$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 the start of the experiment, $t_{bound} = 0$, this is not the case in substructures with a mass growth (21), which capture particles onto permanent orbits over an extended period of time, $t_{bound} \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_{rot}$) is not a conserved quantity in time-dependent systems. The results shown in Fig. 7 call for a heuristic classification of orbits according to 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_{bound} \le 2\tau$, and whose energy $E$ remains negative uninterruptedly until the end of the simulation.<fnref n="7" /> 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 model,<fnref n="8" /> $f_{perm} = N_{perm}/(N_{perm} + N_{temp}) = N_{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.

Fig. 8 plots the number of bound particles located within a fixed spherical radius $r_{out} = 190$ 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_{out}$ matches the local tidal radius (20). 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 (21; 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

<footnote n="7">⁷Notice 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_{bound} \le 2\tau$ is $\lesssim 10$ per cent</footnote>

<footnote n="8">⁸The 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.</footnote>
<figure><img src="image_9.png" /><figcaption>**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 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_{\text{out}} = 100$ 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 ratio becomes approximately constant at $t \gtrsim 5$, $T$, where $T$ is Smoluchowski's time scale.</figcaption></figure>

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_{\text{bound}} \approx \Delta N_{\text{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 (21) 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_{\text{bound}} > \Delta N_{\text{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 $|\Phi_\bullet/\Phi_\bullet|^{-1} = |\dot{M}_\bullet/M_\bullet|^{-1}_{\tau=0} = \tau \sim T$. Further theoretical work is need to understand how this capture mechanism works.

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(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 (8) (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 (6) (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_* \sim |\Phi_\bullet|^{3/2}$. Similarly, the lower panel of Fig. 9 shows that the velocity dispersion profile goes as $\sigma_* \sim |\Phi_\bullet|^{1/2}$, which in the case of a Hernquist substructure potential $\Phi_\bullet$ leads to the velocity dispersion profile (7) plotted with a black-dashed line.

## 3.4 Accreted substructures on an eccentric orbit

In the previous sections, we carry out experiments where a sub-structure 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 experimental setup invariably leads to a population of field particles that are captured on to 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)
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Upper panel: Average phase-space density $Q_* = n_*/\sigma_*^3 = n\delta_*/\sigma_*^3$ as a function of distance from the potential $\Phi_*$ measured from the models shown in Fig. 8. Distances are measured in units of the thermal critical radius $r_c$, equation (9). Note that $Q_*$ remains approximately constant within the volume element, and that its value agrees well with the theoretical expectation (8; black-dashed line). Middle panel: Density enhancement profile $\delta_*(r) = n_*(r)/n$. The black-dashed line shows the theoretical curve (6). For reference, we mark with a dotted line the background density value $\delta_* = 1$. Lower panel: Velocity dispersion profiles of the models shown above normalized by the local velocity dispersion of the field, $\sigma_*(r)/\sigma$. Black-dashed line shows the theoretical profile (7).</figcaption></figure>

potential ($\gamma = 1$) with a mass $M^{\text{dSph}} = 3 \times 10^9 \, \text{M}_\odot$ and a scale radius $c^{\text{dSph}} = 2 \, \text{kpc}$. For illustration, we show two experiments where the substructure is inserted in the galactic potential at apocentre, $R_{\text{apo}} = 1380 \, \text{pc}$. The initial velocity is set to $V_*=7.9 \, \text{km s}^{-1}$, which leads to a pericentre of $R_{\text{peri}} = 165 \, \text{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_{\text{apo}}) = R_{\text{apo}}/V_c(R_{\text{apo}}) = 34 \, \text{Myr}$, where $V_c(R_{\text{apo}}) = 39.4 \, \text{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 Section 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 \, \text{pc}$, which falls off very steeply beyond $R \gtrsim 600 \, \text{pc}$. The luminosity-averaged velocity dispersion of these particles is $(v^2)^{1/2}/3 = 11.7 \, \text{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 \, \text{pc}$,

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Left-hand panels: Initial distribution of field particles (grey dots) in dynamical equilibrium within a Dehnen potential with a total mass $M^{\text{dSph}} = 3 \times 10^9 \, \text{M}_\odot$, a scale radius $c^{\text{dSph}} = 2 \, \text{kpc}$ and a central density slope $\gamma^{\text{dSph}} = 1$. Right-hand panels: Distribution of particles that become bound to a Hernquist potential with $M_*=3 \times 10^7 \, \text{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_{\text{out}} = 300 \, \text{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 kpc. Notice that particles captured from a truncated galaxy exhibit a flattened shape roughly oriented along the orbital motion.</figcaption></figure>

which leads to a slightly higher luminosity-averaged velocity dispersion $(v^2)^{1/2}/3 = 19.9 \, \text{km s}^{-1}$.

The substructure sources a static Hernquist potential with a mass $M_*=3 \times 10^7 \, \text{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_*= (1-\kappa)r_c = 0.1 r_c$. To find the thermal critical radii of the two models we insert $\sigma = (v^2)^{1/2}/3$ and $V_*$ in equation (9), which yields $r_c = 390 \, \text{pc}$ and $980 \, \text{pc}$ for the extended and truncated galaxy models, respectively.

Left-hand 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-hand 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 overdensity of bound particles has a close-to-spherical shape and a relatively large size. In contrast, the overdensity in the
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Number of bound particles located within a volume size $r_{\text{out}} = 0.3$ kpc as a function of time for the models plotted in Fig. 10. Dotted lines show the total (bound plus unbound) number of field particles within the volume, $r < r_{\text{out}}$. Black-dashed lines show the theoretical prediction from equation (2), with a 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.</figcaption></figure>

bottom panel has a smaller size and a flattened shape, which appears to be roughly oriented along the orbital motion.

The rate at which field particles are captured in these two models is also strikingly different. Solid lines in Fig. 11 shows the number of bound particles within a volume size $r_{\text{out}} = 300$ pc (marked with black circles in Fig. 10 for reference) as a function of time. In the top panel, we can see that the substructure orbiting in the extended galaxy contains $N_b \sim 3 \times 10^3$ at $t = 0$. By definition, these are *immediate* captures. As the substructure falls into the inner regions of the galaxy, the total number of field particles within the volume (dotted lines) increases systematically until pericentre is reached at $\Omega t \simeq 1.4$. After the pericentre passage, the number of bound particles (solid lines) grows slightly, reaching a plateau at apocentre, $\Omega t \simeq 2.8$.

This situation repeats throughout the next orbital revolution. At the last snapshot of the simulation, the number of bound particles is $N_b(t_f) \sim 4 \times 10^4$, which is a factor $\sim 10$ increases with respect to the value at $t = 0$. Remarkably, the theoretical expectation given by equation (2; black-dashed line) barely changes as a function of orbital phase in spite of the large variation of galactocentric distances and velocities. This 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/(2\sigma^2(R))]$. As the substructure plunges into the inner regions of the galaxy the orbital velocity increases, making capture less efficient. Simultaneously, the number of field particles around the substructure grows, with both effects nearly balancing out in equation (2). The reverse process takes place from pericentre to apocentre.

The substructure model orbiting within a truncated galaxy experiences a very different evolution. In this case, bottom panel of Fig. 11 shows that capture does not occur until the substructure is close to its first orbital pericentre. As expected, the number of immediate captures is zero by construction. After each pericentric passage, we observe a very rapid increase in the population of bound particles, whose size can grow by $\sim 3$ orders of magnitude. However, a large fraction of these particles are quickly lost to galactic tides. At first apogalacticon, the number of bound particles has stabilized at $N_b \sim 10^2$. As it reaches its second pericentre, it grows again up to $N_b \sim 5 \times 10^3$, falling back to $N_b(t_f) \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 Section 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$ pc. Furthermore, it 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 model shown in the bottom panel. The latter approximation is particularly poor after the first pericentric passage. Here, equation (2) predicts $N_b = 0$ because the density of field particles at this location is $n(R_{\text{apo}}) = 0$ (see bottom-left panel of Fig. 10). In contrast, our numerical experiment shows that the substructure has captured $N_b \sim 10^3$ particles after the first orbital revolution. This mismatch arises again at the second apocentre, which corresponds to the final snapshot of the simulation. Here, all particles within the volume under observation (dotted line) are bound to the substructure (solid line). It is worth stressing that region of the galaxy was empty from field particles at $t = 0$. Looking the bottom-right panel of Fig. 10, we find field particles beyond the original truncation radius of the galaxy are either still bound to the substructure (red dots), or were ‘scooped’ from the inner regions and released back to the galaxy by tidal forces (black dots).

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_*(r=0)$, and measure distances in units of the substructure scale radius ($c_*$). 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_{\text{apo}} = 1380$ pc and a tangential velocity $V_* = 7.9 \text{ km s}^{-1}$ at $t = 0$, such that all individual realizations have a common pericentre of $R_{\text{peri}} = 165$ pc.
<figure><img src="image_13.png" /><figcaption>**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_*$). 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 (6). In contrast, substructures evolving in a truncated galaxy generate a local overdensity that falls off steeply at large distances $\delta_* \sim r^{-3}$ at $r \gg c_*$.</figcaption></figure>

Particles captured by the substructure model orbiting in an extended galaxy exhibit a profile shape that roughly scales as $(1 + r/c_*)^{-3/2}$, which matches equation (6). 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. However, the normalization of the profile is not so well matched by the analytical formula. In particular, Fig. 11 shows that equation (6) tends to underpredict the number of field particles bound to the substructure at all points of the orbit.

This mismatch is considerably stronger in substructures orbiting a in a truncated galaxy. In these objects, captured field particles show a profile that falls off more steeply than predicted by equation (6) at large distances, $r \gtrsim c_*$, whereas at small distances, $r \ll c_*$, it converges to the profile of the extended galaxy model. The different behaviour may be traced back to the lack of permanent captures in the truncated galaxy models, which tend to smooth out the overall density enhancement profile (see Fig. 6).

In general, the numerical experiments in this Section indicate that the distribution of field stars captured by substructures moving on eccentric orbits vary as a function of time in different ways depending on whether these objects are form in a pre-existing sea of field particles, or are accreted from outside the galaxy.

## 3.5 Compact versus 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 Section 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_* = 10^6 M_\odot$ moving on a circular orbit $R_* = 0.5$ kpc in the Fornax dSph-like galaxy model introduced in Section 3.1. Recall that the tidal radius given by equation (20) is $r_t \approx 190$ 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_* < 1$) at any distance from the substructure. Application of equation (9) yields a thermal critical radius $r_c = 18.9$ pc, which translates into scale-radii $c_* = 3.78$ pc and $c_* = 28.3$ pc for the ‘compact’ and ‘fluffy’ models, respectively. Notice that the time-scale needed to cross the overdensity size by a random field particle is much shorter than the orbital time around the dSph galaxy potential, $T_c = r_c/\sigma^{\text{dSph}} = 18.9$ pc/(12.2 km s$^{-1}$) $\approx 1.5$ 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 (20), $r_{\text{out}} = r_t = 190$ pc. We run both models from $t=0$ until $t=50$ T.

Fig. 13 highlights the impact of the substructure *compactness* ($\kappa = 1 - c_*/r_c$) on the spatial and kinematic distribution of stars with $E < 0$. As predicted in Section 2.1, at a fixed mass only compact ($\kappa > 0$) objects lead to the formation of overdense 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 (6) and velocity dispersion (7) profiles derived analytically (black-solid lines) against the numerical values (orange dots). For reference, Fig. 13 marks with vertical arrows the location of $c_*$ (in blue) and $r_c$ (in red), respectively. At large distances $r \gg c_*$, we find that the density enhancement approaches the Keplerian behaviour $\delta_* \sim r^{-3/2}$ and the velocity dispersion drops off as $\sigma_* \sim r^{-1/2}$, whereas at small distances $r \ll c_*$ both profiles become roughly constant. However, while in the compact model ($\kappa = +0.8$) trapped stars reach a relatively high density, $\delta_{*,0} = \delta_*(r=0) \approx 9.0$, and a hot velocity dispersion, $\sigma_{*,0}/\sigma = \sigma_*(r=0)/\sigma \simeq 1.9$, stars captured by the fluffy counterpart with a negative compactness ($\kappa = -0.5$) show subdominant densities, $\delta_{*,0} = 0.44$, and a cold dispersion, $\sigma_{*,0}/\sigma \simeq 0.69$, with respect to the field. We discuss the relevance of this result for the detection of dark substructures in Section 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 (~$10^5 M_\odot$) globular clusters (e.g. Larsen, Brodie & Strader 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
<figure><img src="image_14.png" /><figcaption>**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 M_{\odot}$ moving on a circular orbit at a fixed distance $R = 0.5$ kpc from the centre of the Fornax dSph (see Section 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$ pc is the thermal critical radius measured from (9) (black-dashed circle). *Middle left panels:* Velocity dispersion map for the models shown in the left-hand 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 (6) is shown with a black line. *Right-hand panels:* Velocity dispersion profiles in spherical coordinates of the above models. Black lines show the profiles predicted by equation (7). Captured particles have isotropic velocities on account of their chaotic motions in the dwarf potential (see text).</figcaption></figure>

composed of stars or background galaxies (e.g. Stetson, Hesser & Smecker-Hane 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 M_{\odot}$, and a half-light radius of $r_{\tilde{h}}^{F6} = 11 \pm 1.4$ pc (Wang et al. 2019). Crucially, the metallicity $[Fe/H] = -0.71 \pm 0.05$ and age $\sim 2$ 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 five globular clusters, which have much lower metallicities that range between $-2.5 < [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$ 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 ongoing 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 overdensity 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 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 overdensity, it must be compact enough, i.e. $\kappa > 0$ or $c_{\bullet} \lesssim r_{\epsilon} \sim 20$ pc (for illustrative examples, see the numerical models

presented Section 3.1). Interestingly, this size is smaller than the average peak-velocity radius of field CDM haloes of comparable masses, i.e. $r_{\text{max}} \sim 100$ pc for $M_{200} \sim 10^6 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 ($[Fe/H] \simeq -0.6$) and young stellar ages ($\sim 1.5$ 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}^{dSph} \sim (6 \pm 1) \times 10^4$ (Crnojević et al. 2016), a half-light radius of $r_{\tilde{h}}^{dSph} = 299 \pm 12$ pc (Simon et al. 2021), and low velocity dispersion $\sigma^{dSph} = 6.9 \pm 1.1$ km s$^{-1}$. Surprisingly for such a faint galaxy, it contains a single stellar cluster located at $R = 23 \pm 3$ 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$ Gyr, and extremely
metal-poor [Fe/H] = −2.75 ± 0.2 (Weisz, Savino & Dolphin 2023). It has a remarkably large size, $r_h = 15 \pm 1$ 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, Irwin & Kunkel 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 ~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 Section 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 Section 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, Giocoli & Tormen 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 Section 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 Section 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 $|\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 galaxy models, numerical experiments in Section 3.4) show that the steady-state approximation may be reasonably accurate away from orbital pericentre. In truncated galaxy models, it is the local approximation that fails away from pericentre (see Fig. 10).

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 Section 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 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_{vir} \sim 10^9$ M$_\odot$ (e.g. Peñarrubia, McConnachie & Navarro 2008a; Errani et al. 2018). Adopting a mass ratio between the Fornax dSph and the MW of $\sim 10^9$ M$_\odot/10^{12}$ 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$ M$_\odot$ enclosed within its virial radius, $r_{vir}^{dSph} \sim 30$ 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, Navarro & McConnachie 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$ 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, Taibi & Thomas 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 abovementioned 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_{\odot}/M_{\odot} < 10^8$ (Erkal et al. 2016), while strong-lens observations are sensitive to subhaloes with masses $M_{\odot} \gtrsim 6 \times 10^8 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_{\odot} \sim 10^4-10^5 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_{hm} = 3 \times 10^8 M_{\odot} \left( \frac{m_{DM}}{3.3 \text{ keV}} \right)^{-3.33}.$$

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

In addition, the internal structure of these systems may also reveal whether DM particles experience self-interactions. For example, dark substructures with $M_{\odot} \sim 10^6 M_{\odot}$ must be sufficiently compact, $c_{\odot} \lesssim 20$ 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_{\odot} = M_{\odot}/(2\pi c_{\odot}^2) \gtrsim 20 M_{\odot} \text{ pc}^{-3}$, which is higher than CDM subhaloes with similar masses. More precisely, field haloes with a mass of $M_{200} = 10^6 M_{\odot}$ at redshift $z=0$ have a mean peak velocity radius of $r_{\text{max}} \sim 100$ pc, equating this radius to the substructure scale radius, $c_{\odot} = r_{\text{max}}$, yields a considerably lower characteristic density of $\rho_{\odot} \sim 0.1 M_{\odot} \text{ 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-subhaloes in satellite galaxies are expected to be tidally processed, and more concentrated than field haloes 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 Section 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 (14). 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. (2018) 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 Zoutendijk et al. 2020), whereas Marsh & Niemeyer (2019) point out that the cluster could also be used to test quantum fluctuations of ultralight 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_{\odot}$ when their gravitational energy $E = v^2/2 + \Phi_{\odot}(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_{\odot}$ 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_{\odot}/r_{\odot} > 0$, where $r_{\odot}$ and $c_{\odot}$ are the thermal critical radius (9) and scale radius of $\Phi_{\odot}$, 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_{\odot}$ for arbitrarily long times. In static models where $\Phi_{\odot}$ 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 Fukushima & 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|_{t=0}^{-1}$ 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_* = n_*/\sigma_*^3 \approx \text{const.}$ In steady state, these stars follow number density and dispersion profiles that scale as $n_* \sim |\Phi_*|^{3/2}$, and $\sigma_* \sim |\Phi_*|^{1/2}$, respectively (see Section 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 $dN/dM_\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 (14), which scales as $dN/dM_* \sim M_*^{-\beta}$, where $\beta = (2\alpha + 1)/3 \approx 1.6$. However, only subhaloes that are sufficiently 'compact' ($\kappa > 0$) may materialize as localized stellar overdensities in the galactic background ($\delta_* > 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 'fluffy' substructures ($\kappa < 0$) have densities below the background ($\delta_* < 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 (($\sigma_\sigma^2)^{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-hand 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 Section 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 may not 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. Follow-up work is needed to inspect the detectability of 'dark' subhaloes in dSphs with masses below the star formation threshold, as well as observational campaigns aimed at detecting & characterizing those objects. In Section 4.3 we discuss how the detection and characterization of dark subhaloes may provide unique & unprecedented constraints on the particle mass and cross section for a large range of DM candidates.

## ACKNOWLEDGEMENTS

We are grateful to the referee for the constructive criticism and the suggestions that helped to improve the paper. We thank Sergey Koposov, Alan McConnachie and Frank van den Bosch for their valuable comments and suggestions. MGW acknowledges support from National Science Foundation (NSF) grants AST-1813881, AST-1909584 and AST-2206046. MG 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.

In memory of Teo Peñarrubia Reihl (16-3-2018, 6-1-2024). V.A.E.J.H.

## DATA AVAILABILITY

No data were generated for this study.

## REFERENCES

* Ackermann M. et al., 2014, Phys. Rev. D, 89, 042001

* Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781

* Amorisco N. C., Evans N. W., 2012, ApJ, 756, L2

* Bate N. F., McMonigal B., Lewis G. F., Irwin M. J., Gonzalez-Solares E., Shanks T., Metcalfe N., 2015, MNRAS, 453, 690

* Battaglia G. et al., 2006, A&A, 459, 423

* Battaglia G., Taibi S., Thomas G. F., Fritz T. K., 2022, A&A, 657, A54

* Benitez-Llambay A., Frenk C., 2020, MNRAS, 498, 4887

* Benson A. J., 2017, MNRAS, 467, 3454

* Boekholt T. C. N., Rowan C., Kocsis B., 2023, MNRAS, 518, 5653

* Bonaca A., Hogg D. W., Price-Whelan A. M., Conroy C., 2019, ApJ, 880, 38
Bovill M. S., Ricotti M., 2009, ApJ, 693, 1859

Brandt T. D., 2016, ApJ, 824, L31

Bringmann T., Calore F., Di Mauro M., Donato F., 2014, Phys. Rev. D, 89, 023012

Bullock J. S., Kravtsov A. V., Weinberg D. H., 2000, ApJ, 539, 517

Carlberg R. G., 2013, ApJ, 775, 90

Coleman M., Da Costa G. S., Martínez-Delgado D., Bland-Hawthorn J., 2004, ASPC, 327, 173

Contenta F. et al., 2018, MNRAS, 476, 3124

Crnojević D., Sand D. J., Zaritsky D., Spekkens K., Willman B., Hargis J. R., 2016, ApJ, 824, L14

Daniel K. J., Heggie D. C., Varri A. L., 2017, MNRAS, 468, 1453

de Boer T. J. L., Tolstoy E., Saha A., Olszewski E. W., 2013, A&A, 551, A103

Deason A. J., Belokurov V., Evans N. W., 2011, MNRAS, 416, 2903

Dehnen W., 1993, MNRAS, 265, 250

Demers S., Irwin M. J., Kunkel W. E., 1994, AJ, 108, 1648

Despali G., Giocoli C., Tormen G., 2014, MNRAS, 443, 3208

Diemand J., Moore B., Stadel J., 2005, Nature, 433, 389

Diemer B., Joyce M., 2019, ApJ, 871, 168

Eddington A. S., 1916, MNRAS, 76, 572

Enzi W. et al., 2021, MNRAS, 506, 5848

Erkal D., Belokurov V., Bovy J., Sanders J. L., 2016, MNRAS, 463, 102

Errani R., Navarro J. F., 2021, MNRAS, 505, 18

Errani R., Peñarrubia J., 2020, MNRAS, 491, 4591

Errani R., Peñarrubia J., Laporte C. F. P., Gómez F. A., 2017, MNRAS, 465, L59

Errani R., Peñarrubia J., Walker M. G., 2018, MNRAS, 481, 5073

Fukushige T., Heggie D. C., 2000, MNRAS, 318, 753

Gilman D., Birrer S., Nierenberg A., Treu T., Du X., Benson A., 2020, MNRAS, 491, 6077

Giocoli C., Tormen G., van den Bosch F. C., 2008, MNRAS, 386, 2135

Gondolo P., Silk J., 1999, Phys. Rev. Lett., 83, 1719

Green A. M., Hofmann S., Schwarz D. J., 2005, J. Cosmol. Astropart. Phys., 2005, 003

Heppenheimer T. A., 1975, Icarus, 24, 172

Heppenheimer T. A., Porco C., 1977, Icar, 30, 385

Hernquist L., 1990, ApJ, 356, 359

Hofmann S., Schwarz D. J., Stöcker H., 2001, Phys. Rev. D, 64, 083507

Hunter R. B., 1967, MNRAS, 136, 245

Ibata R. A., Lewis G. F., Irwin M. J., Quinn T., 2002, MNRAS, 332, 915

Jílková L., Portegies Zwart S., Pijlloo T., Hammer M., 2015, MNRAS, 453, 3157

Johnston K. V., Spergel D. N., Haydn C., 2002, ApJ, 570, 656

Kahlhoefer F., Kaplinghat M., Slatyer T. R., Wu C.-L., 2019, J Cosmol. Astropart. Phys., 2019, 010

Kalnajs A. J., 1971, Ap&SS, 13, 279

Koopmans L. V. E., 2005, MNRAS, 363, 1136

Koposov S. E., Belokurov V., Torrealba G., Evans N. W., 2015, ApJ, 805, 130

Küpper A. H. W., Kroupa P., Baumgardt H., Heggie D. C., 2010, MNRAS, 407, 2241

Larsen S. S., Brodie J. P., Strader J., 2012, A&A, 546, A53

Li R., Mo H. J., Fan Z., Yang X., Bosch F. C. v. d., 2013, MNRAS, 430, 3359

Li J., Lai D., Rodet L., 2022, ApJ, 934, 154

Loeb A., Zaldarriaga M., 2005, Phys. Rev. D, 71, 103520

Ludlow A. D., Bose S., Angulo R. E., Wang L., Hellwing W. A., Navarro J. F., Cole S., Frenk C. S., 2016, MNRAS, 460, 1214

Marsh D. J. E., Niemeyer J. C., 2019, Phys. Rev. Lett., 123, 051103

Mateo M. L., 1998, ARA&A, 36, 435

McKernan B., Ford K. E. S., Lyra W., Perets H. B., 2012, MNRAS, 425, 460

Moliné Á., Sánchez-Conde M. A., Palomares-Ruiz S., Prada F., 2017, MNRAS, 466, 4974

Nadler E. O. et al., 2021a, Phys. Rev. Lett., 126, 091101

Nadler E. O., Birrer S., Gilman D., Wechsler R. H., Du X., Benson A., Nierenberg A. M., Treu T., 2021b, ApJ, 917, 7

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Ngan W., Carlberg R. G., Bozek B., Wyse R. F. G., Szalay A. S., Madau P., 2016, ApJ, 818, 194

O'Riordan C. M., Despali G., Vegetti S., Lovell M. R., Moliné Á., 2023, MNRAS, 521, 2342

Orkney M. D. A. et al., 2022, MNRAS, 515, 185

Pace A. B., Martinez G. D., Kaplinghat M., Muñoz R. R., 2014, MNRAS, 442, 1718

Pace A. B., Walker M. G., Koposov S. E., Caldwell N., Mateo M., Olszewski E. W., Bailey J. I., Wang I M.-Y., 2021, ApJ, 923, 77

Peñarrubia J., 2018, MNRAS, 474, 1482

Peñarrubia J., 2019, MNRAS, 484, 5409

Peñarrubia J., 2023, MNRAS, 519, 1955

Peñarrubia J., Benson A. J., 2005, MNRAS, 364, 977

Peñarrubia J., McConnachie A. W., Navarro J. F., 2008a, ApJ, 672, 904

Peñarrubia J., Navarro J. F., McConnachie A. W., 2008b, ApJ, 673, 226

Peñarrubia J., Benson A. J., Walker M. G., Gilmore G., McConnachie A. W., Mayer L., 2010, MNRAS, 406, 1290

Pereira-Wilson M., Navarro J. F., Benítez-Llambay A., Santos-Santos I., 2023, MNRAS, 519, 1425

Petit J.-M., Hénon M., 1986, Icarus, 66, 536

Renaud F., Gieles M., Boily C. M., 2011, MNRAS, 418, 759

Rowan C., Boekholt T., Kocsis B., Haiman Z., 2023, MNRAS, 524, 2770

Rusakov V., Monelli M., Gallart C., Fritz T. K., Ruiz-Lara T., Bernard E. J., Cassisi S., 2021, MNRAS, 502, 642

Scheuer P. A. G., 1957, Proc. Cambridge Phil. Soc., 53, 764

Schneider A., Smith R. E., Macciò A. V., Moore B., 2012, MNRAS, 424, 684

Schmid C., Schwarz D. J., Widerin P., 1999, Phys. Rev. D, 59, 043517

Shapley H., 1939, Proc. Natl. Acad. Sci., 25, 565

Simon J. D., 2019, ARA&A, 57, 375

Simon J. D. et al., 2021, ApJ, 908, 18

Smoluchowski M. V., 1916, Zeitschrift fur Physik, 17, 557

Springel V. et al., 2008, MNRAS, 391, 1685

Stetson P. B., Hesser J. E., Smecker-Hane T. A., 1998, PASP, 110, 533

Suetsugu R., Ohtsuki K., 2013, MNRAS, 431, 1709

Szebehely V., 1967, Theory of Orbits: The Restricted Problem of Three Bodies. Academic Press Inc., New York

van den Bosch F. C., Ogiya G., 2018, MNRAS, 475, 4066

Vegetti S., Koopmans L. V. E., 2009, MNRAS, 400, 1583

Vegetti S., Koopmans L. V. E., Auger M. W., Treu T., Bolton A. S., 2014, MNRAS, 442, 2017

Verner G., Demers S., Hardy E., Kunkel W. E., 1981, AJ, 86, 357

Villasenor B., Robertson B., Madau P., Schneider E., 2022, preprint (arXiv:2209.14220)

Walker M. G., Peñarrubia J., 2011, ApJ, 742, 20

Walker M. G., Mateo M., Olszewski E. W., Gnedin O. Y., Wang X., Sen B., Woodroofe M., 2007, ApJ, 667, L53

Wang M. Y. et al., 2019, ApJ, 881, 118

Wang M. et al., 2020, BAAS, 235, 264.06

Weerasooriya S., Bovill M. S., Benson A., Musick A. M., Ricotti M., 2023, ApJ, 948, 87

Weinberg M. D., 1986, ApJ, 300, 93

Weinberg M. D., 1994, AJ, 108, 1398

Weisz D. R., Savino A., Dolphin A. E., 2023, ApJ, 948, 50

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341

Yoon J. H., Johnston K. V., Hogg D. W., 2011, ApJ, 731, 58

Zhao H., 1996, MNRAS, 278, 488

Zoutendijk S. L. et al., 2020, A&A, 635, A107

This paper has been typeset from a TeX/LA TeX file prepared by the author.