# Micro galaxies as a falsifiable prediction of LCDM cosmology

Raphaël Errani, Rodrigo Ibata, Julio F Navarro, Jorge Peñarrubia, Matthew

G Walker

► To cite this version:

Raphaël Errani, Rodrigo Ibata, Julio F Navarro, Jorge Peñarrubia, Matthew G Walker. Micro galaxies as a falsifiable prediction of LCDM cosmology. The Astrophysical Journal, 2024, 968 (2), pp.89. 10.3847/1538-4357/ad402d. hal-04329360

HAL Id: hal-04329360

https://hal.science/hal-04329360v1

Submitted on 18 Jul 2025

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
<figure><img src="image_2.png" /></figure>

# Microgalaxies in LCDM

Raphaël Errani<fnref n="1" />,<fnref n="2" />, Rodrigo Ibata<fnref n="2" /> [ID], Julio F. Navarro<fnref n="3" /> [ID], Jorge Peñarrubia<fnref n="4" />, and Matthew G. Walker<fnref n="1" /> [ID]

Received 2023 November 24; revised 2024 March 18; accepted 2024 April 16; published 2024 June 13

## Abstract

A fundamental prediction of the Lambda cold dark matter cosmology is the centrally divergent cuspy density profile of dark matter haloes. Density cusps render cold dark matter haloes resilient to tides, and protect dwarf galaxies embedded in them from full tidal disruption. The hierarchical assembly history of the Milky Way may therefore give rise to a population of “microgalaxies”; i.e., heavily stripped remnants of early accreted satellites, which can reach arbitrarily low luminosity. Assuming that the progenitor systems are dark matter dominated, we use an empirical formalism for tidal stripping to predict the evolution of the luminosity, size, and velocity dispersion of such remnants, tracing their tidal evolution across multiple orders of magnitude in mass and size. The evolutionary tracks depend sensitively on the progenitor distribution of stellar binding energies. We explore three cases that likely bracket most realistic models of dwarf galaxies: one where the energy distribution of the most tightly bound stars follows that of the dark matter, and two where stars are defined by either an exponential density or surface brightness profile. The tidal evolution in the size-velocity dispersion plane is quite similar for these three models, although their remnants may differ widely in luminosity. Microgalaxies are therefore best distinguished from globular clusters by the presence of dark matter; either directly, by measuring their velocity dispersion, or indirectly, by examining their tidal resilience. Our work highlights the need for further theoretical and observational constraints on the stellar energy distribution in dwarf galaxies.

*Unified Astronomy Thesaurus concepts: Cold dark matter (265); Dwarf spheroidal galaxies (420); Low surface brightness galaxies (940); the Milky Way (1054); N-body simulations (1083); Star clusters (1567); Tidal disruption (1696)*

## 1. Introduction

In the Lambda cold dark matter (LCDM) cosmology, galaxies are predicted to form in the potential wells of dark matter overdensities (White & Rees 1978). Through a history of accretion and merger events, these overdensities give rise to a complex clustering hierarchy of *haloes* and *subhaloes* (for reviews, see, e.g., Frenk & White 2012; Zavala & Frenk 2019).

Computer simulations suggest that the internal mass distribution of haloes is well approximated by the Navarro–Frenk–White (NFW) profile (Navarro et al. 1996, 1997). This means that cold dark matter (CDM) haloes are predicted to have remarkably high central densities: for NFW profiles, the density formally diverges as $d \ln \rho_{\text{NFW}}/d \ln r \to -1$ for $r \to 0$. In contrast, some other theories of dark matter, such as self-interacting dark matter (SIDM), predict haloes with constant-density cores, i.e., $d \ln \rho_{\text{core}}/d \ln r \to 0$ for $r \to 0$ (Burkert 2000; Spergel & Steinhardt 2000; Colín et al. 2002).

Dwarf spheroidal (dSph) galaxies are promising observational probes of the properties of galactic dark matter subhaloes: their dynamical mass-to-light ratios, inferred from stellar kinematics, suggest that most dSphs are dark matter-dominated objects (Mateo et al. 1993; Walker et al. 2007; McConnachie 2012; for reviews, see Mateo 1998; Simon 2019; Battaglia & Nipoti 2022).

For the Fornax and Sculptor dSphs, several studies find evidence for cored dark matter density profiles (Walker & Peñarrubia 2011; Amorisco & Evans 2012; Amorisco et al. 2013; Diakogiannis et al. 2017; Pascale et al. 2018; Read et al. 2019). Instead, for the Draco dSph, most studies favor a cuspy density profile (Jardel et al. 2013; Read et al. 2018; Massari et al. 2020; Hayashi et al. 2020). Whether this “diversity” (Oman et al. 2015) is due to the intrinsic properties of dark matter, or driven by baryonic effects (see, e.g., Santos-Santos et al. 2020), remains a matter of debate.

An independent approach to infer properties of dark matter substructures on galactic scales is to study their response to tides. Dwarf galaxies that have been accreted onto the Milky Way are subject to tidal forces, which induce mass loss and structural changes (Peñarrubia et al. 2008). In the case of LCDM, the centrally divergent density profile of subhaloes renders them very resilient to the effects of tides. Indeed, numerical simulations suggest that smooth tidal fields do not fully disrupt NFW subhaloes (see, e.g., Peñarrubia et al. 2010; van den Bosch et al. 2018; Errani & Peñarrubia 2020, hereafter EP20) but rather lead them to asymptotically approach a stable remnant state (Errani & Navarro 2021, hereafter EN21).

On the other hand, for cored dark matter substructures, tides may trigger a runaway process that leads to their full disruption (Peñarrubia et al. 2010; Errani et al. 2023a). The mere existence of heavily stripped substructures can hence help to constrain the density structure of haloes and thereby the nature of dark matter. Crucially, stars embedded inside cuspy dark matter subhaloes will be protected from full tidal disruption if their

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA; errani@cmu.edu</footnote>

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

<footnote n="3">³Department of Physics and Astronomy, University of Victoria, Victoria, BC V8P 5C2, Canada</footnote>

<footnote n="4">⁴Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>
distribution of binding energies within the subhalo extends all the way to the most-bound states (Errani et al. 2022, hereafter E +22). Heavily stripped dwarf galaxies may thereby give rise to a population of “micro galaxies” (EP20), i.e., comoving groups of stars embedded in heavily stripped dark matter subhaloes.

Detailed modeling of the structural changes that dwarf galaxies undergo when subject to strong tidal fields is complicated by the limited particle number and spatial resolution of *N*-body simulations: insufficient resolution leads to an inaccurate evolution of structural parameters (EP20; EN21) and may result in the artificial disruption of substructures (van den Bosch et al. 2018). This is of particular relevance for dwarf galaxies that pass through the inner regions of the Milky Way, where the combined tidal field of disk, bulge, and halo is strongest, significantly increasing the rate of tidal stripping (D’Onghia et al. 2010; Errani et al. 2017; Kelley et al. 2019).

In this work, we model the observable properties of the faintest Milky Way satellites, assuming that they are embedded in LCDM subhaloes, and accurately following their tidal evolution over many orders of magnitude in mass, size, and luminosity. Our models use the empirical tidal stripping framework introduced in E+22, which allows us to make predictions for spatial and mass scales inaccessible to current cosmological simulations.

We discuss our findings in the light of recently discovered Milky Way satellites (Torrealba et al. 2019; Mau et al. 2020; Cerny et al. 2023, 2023; Smith et al. 2024), which have sizes and luminosities at the interface between the globular cluster and dwarf galaxy regime. Could these satellites have their origin in accreted dwarfs and be among the first microgalaxy candidates?

The paper is structured as follows. We show how resolution limits of cosmological *N*-body simulations affect the tidal survival of DM subhaloes in Section 2. In Section 3, we summarize the empirical tidal stripping framework used in this study, and apply it to model the tidal evolution of dwarf galaxies over many orders of magnitude in mass loss. In Section 4, we compare the model predictions against observed properties of Milky Way satellites, paying particular attention to systems at the interface of the globular cluster and dwarf galaxy regime. Finally, we summarize our main results and conclusions in Section 5.

## 2. Artificial Tidal Disruption in Cosmological N-body Simulations

Convergence in *N*-body simulations depends on a combination of the spatial resolution (force softening length $\epsilon$, or particle mesh cell size $\Delta x$), particle number *N*, and time step $\Delta t$ (see, e.g., Power et al. 2003; Springel et al. 2008; van den Bosch & Ogiya 2018, EN21). In this section, we perform controlled simulations to illustrate that the resolution of current cosmological *N*-body simulations is not enough to reliably trace the tidal evolution of most subhaloes that pass through the inner regions of the Milky Way.

### 2.1. Pericentric/Apocentric Radii and Accretion Redshift

Cosmological simulations indicate that the present-day pericentric and apocentric radii of subhaloes correlate strongly with the time when they were first accreted onto the main halo. Figure 1 illustrates this result for the case of the Aquarius-A2

simulation (Springel et al. 2008), where we have tracked the orbits of all subhaloes with virial<fnref n="5" /> mass $M_{200} \gtrsim 10^8 M_\odot$ accreted onto the A2 main halo (as in Errani et al. 2017). Subhalo orbits have been computed as point-masses in a time-evolving analytical potential fitted to the main halo (Equations (22) and (23) in Buist & Helmi 2016). This setup allows us to follow the orbits of all subhaloes resolved at accretion until $z=0$, without losing subhaloes to artificial disruption.

The integration includes the adiabatic orbital contraction in the growing main halo potential, but neglects the effects of dynamical friction, which generally would affect only the few most massive satellites (Peñarrubia & Benson 2005), reducing their orbital radii even further. Clearly, subhaloes that orbit the inner regions of the main halo have been on average accreted earlier and have been therefore subjected to tidal effects for longer.

Our results show<fnref n="6" /> that, on average, subhaloes on orbits with pericenters $r_{\text{peri}} \lesssim 10$ kpc were accreted more than $\gtrsim 10$ Gyr ago. A subhalo on an orbit with $r_{\text{peri}} = 10$ kpc and $r_{\text{apo}} = 50$ kpc has completed more than $\gtrsim 15$ orbital periods since accretion; subhaloes in the inner regions of the Milky Way have been subjected to strong tidal fields for extended periods of time.

### 2.2. Tidal Disruption and Numerical Resolution

To illustrate the effects of insufficient numerical resolution on tidal evolution, we evolve four different *N*-body realizations of an NFW subhalo in a static isothermal host potential with a circular velocity $V_{\text{host}} = 220$ km s⁻¹ (see Equation (5)). The *N*-body realizations differ only in the number of *N*-body particles used to simulate the subhalo.

#### 2.2.1. Initial Conditions

As an example, we consider an NFW subhalo with a virial mass of $M_{200} \approx 7 \times 10^9 M_\odot$, close to the redshift $z=0$ hydrogen cooling limit (Benitez-Llambay & Frenk 2020). For an average concentration of $c = r_{200}/r_s \approx 13$ (Ludlow et al. 2016), this corresponds to an initial characteristic mass and size of $M_{\text{mx0}} \equiv M(<r_{\text{mx0}}) = 2 \times 10^9 M_\odot$ and $r_{\text{mx0}} = 7$ kpc, respectively (see Equation (1) for a definition).

Using the code introduced in EP20, available online,<fnref n="7" /> we generate four *N*-body realizations of this NFW subhalo: For the first realization, labeled A1, we choose an *N*-body particle mass identical to that used in the Aquarius-A1 simulation, $m_p = 1.71 \times 10^3 M_\odot$. Similarly, for the second realization, labeled A2, we use $m_p = 1.37 \times 10^4 M_\odot$, matching the particle mass of the Aquarius-A2 simulation. The models A3 and A4 adopt particle masses that match those of the Aquarius-A3 and A4 resolution levels, respectively. The parameters of these

<footnote n="5">⁵ The virial radius $r_{200}$ encloses a mass $M_{200}$ so that the mean density $M_{200}/(4/3 \pi r_{200}^3)$ is 200 × larger than the critical density for closure, $\rho_{\text{crit}}$. At redshift $z=0$, $\rho_{\text{crit}} = 3H_0^2/(8\pi G)$.</footnote>

<footnote n="6">⁶ Comparison against the orbits of Milky Way (MW) satellite galaxies (Li et al. 2021; Battaglia et al. 2022) suggests that the subhaloes in our controlled Aq-A2 setup are on average on more radial orbits than MW satellites. This resembles the “tangential velocity excess” of MW satellites discussed in Cautun & Frenk (2017). Note, however, that our setup assumes a spherically symmetric halo and does not include a disk, which may play an important role in shaping the orbital anisotropies of satellites in the inner regions of our Galaxy (see Riley et al. 2019).</footnote>

<footnote n="7">⁷ The code used to generate *N*-body models as described in EP20, as well as an implementation of the EN21 tidal evolution model, is available at https://github.com/rrerrani.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Distribution of apocenter $r_{apo}$ (■) and pericenter $r_{peri}$ (■) for subhaloes in a Milky Way–like host halo as a function of accretion redshift $z_{acc}$ (assuming Aquarius cosmological parameters, $H_0 = 73$ km s$^{-1}$Mpc$^{-1}$, $\Omega_m = 0.25$, $\Omega_\Lambda = 0.75$). The distributions are computed from subhaloes with peak virial masses of $M_{200} \ge 10^8 M_\odot$ in the Aquarius-A2 simulation. Median radii are shown as solid lines, with shaded regions corresponding to the 16th–84th percentiles. On average, subhaloes accreted at earlier $z_{acc}$ have smaller apocenters and pericenters than those accreted more recently. Subhaloes with $r_{peri} < 20$ kpc have $z_{acc} \gtrsim 2$ (i.e., they were accreted $\gtrsim 10$ Gyr ago, see top axis). For reference, the evolution of the host halo virial radius $r_{200}$ is shown as a black solid curve.</figcaption></figure>

simulations are summarized in Table 1. As the NFW profile has a divergent total mass, we taper the profile exponentially beyond $10 r_s$ (where by $r_s$ we denote the NFW scale radius as in Equation (1)). Consequently, the A1, A2, A3, and A4 realizations of our example subhalo consist of $N = 4.4 \times 10^6$, $5.6 \times 10^5$, $1.5 \times 10^5$, and $1.9 \times 10^4$ particles, respectively.

### 2.2.2. Simulation Code

The simulations are carried out using the particle-mesh code SUPERBOX (Fellhauer et al. 2000). The code employs a high-resolution cubic grid comoving with the subhalo, with a cell size $\Delta x \approx r_{mx0}/256 \approx 27$ pc. SUPERBOX uses two additional grids, one comoving of medium resolution ($10 r_{mx0}/256$), and a fixed low-resolution ($500$ kpc/$256$) grid containing the entire simulation box. For reference, the force softening length of the Aquarius-A1 simulation equals $\epsilon = 22$ pc (Springel et al. 2008), similar to the particle-mesh cell size $\Delta x \approx 27$ pc of the highest-resolving grid used in our numerical experiments. Note that the only parameter varied in the simulations discussed in this section is the N-body particle mass $m_p$, and thereby the total number of N-body particles in each simulation.

### 2.2.3. Simulation Results

The NFW subhalo models are placed on an orbit with pericentric and apocentric distance $r_{peri} = 10$ kpc and $r_{apo} = 50$ kpc, respectively, and evolved for 13 Gyr. Figure 2 shows the evolution of $M_{mx}$, the bound<fnref n="8" /> mass within the radius, $r_{mx}$, where the circular velocity peaks, as a function of time. The four N-body realizations are shown in different shades of blue. Over the first two decades in mass loss, all four N-body

<table><caption>**Table 1**
Parameters of the Example Subhalo Model Considered in the N-body Experiments of Section 2 and the Empirical Models of Section 3</caption><tbody><tr><td rowspan="3">SUBHALO</td><td>Profile</td><td>NFW (Equation (1))</td></tr><tr><td>$M_{mx0}$, $V_{mx0}$, $r_{mx0}$</td><td>$2 \times 10^9 M_\odot$, 35 km s<sup>-1</sup>, 7 kpc</td></tr><tr><td>N-body $\Delta x$</td><td>$r_{mx0}/256 \approx 27$ pc</td></tr><tr><td rowspan="4">HOST</td><td rowspan="4">N-body $m_p$</td><td>A1: $1.71 \times 10^3 M_\odot$</td></tr><tr><td>A2: $1.37 \times 10^4 M_\odot$</td></tr><tr><td>A3: $4.91 \times 10^4 M_\odot$</td></tr><tr><td>A4: $3.93 \times 10^5 M_\odot$</td></tr><tr><td>Profile</td><td>Isothermal (Equation (5))</td></tr><tr><td>V<sub>host</sub></td><td>220 km s<sup>-1</sup></td></tr><tr><td>r<sub>peri</sub></td><td>10 kpc</td></tr><tr><td>r<sub>peri</sub>: r<sub>apo</sub></td><td>1: 5</td></tr></tbody></table>

**Note.** The initial subhalo mass $M_{mx0}$ is chosen to approximately match the (redshift $z=0$) hydrogen cooling limit (Benítez-Llambay et al. 2019), and the initial subhalo characteristic size $r_{mx0}$ is that of a $z=0$ average-concentration halo (Ludlow et al. 2016). Listed also are the host halo and orbit properties underlying the N-body and empirical models.

realizations show a similar mass evolution. Beyond that point, however, the evolution differs between the models. The A4 model fully disrupts after only ~1 Gyr; the A2 model fully disrupts after ~2.5 Gyr. The A1 model survives slightly longer, but eventually also disrupts after a total of ~3.5 Gyr of evolution.

The increase in N-body particle number has hence delayed, but not prevented, the full disruption of the subhalo. This numerical experiment suggests that even simulations with a resolution as high as that of the Aquarius-A1 simulation will suffer from the artificial depletion of substructures on orbits that reach the innermost regions of the Galaxy. Even before the subhalo is (artificially) disrupted, its structural evolution is compromised by insufficient resolution: the convergence study listed in Appendix A of EN21 shows that, once the number of bound particles within $r_{mx}$ drops below $N_{mx} < 3000$, resolution artifacts result in subhalo densities being systematically underestimated. Evidently, a tool different from classical cosmological N-body simulations is needed to model the tidal remnants of heavily stripped dwarf galaxies.

To address this issue, we use an empirical model for the tidal evolution of subhaloes (and embedded dwarf galaxies), which allows us to study their tidal evolution over many orders of magnitude in mass loss. A black solid curve in Figure 2 shows the evolution of the subhalo computed using the empirical EN21 model, available online.⁷ The model suggests that mass loss keeps decelerating, and that a stable remnant state is asymptotically approached. The details of this model, as well as its extension to model the evolution of dwarf galaxies, are summarized in the next chapter.

## 3. Tidal Evolution Model

We summarize now the empirical model for tidal stripping used in this work.

### 3.1. Model for the Evolution of the Dark Matter Component

For the evolution of the dark matter component, we rely on the empirical model introduced and tested in EN21. In the initial conditions, the subhalo is assumed to follow an NFW

<footnote n="8">⁸ The bound mass is computed by (1) defining the subhalo center through the shrinking spheres method (Power et al. 2003), (2) computing the subhalo potential under the assumption of spherical symmetry, (3) removing particles unbound in this potential. These steps are iterated until convergence.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Tidal stripping of a subhalo with an initial mass of $M_{mx0} = 2 \times 10^9 M_\odot$ in a Milky Way-like host (see Table 1), computed using N-body simulations (blue curves) and the empirical EN21 model (black solid curve). The four N-body simulations have particle masses $m_p$ chosen to match those of the Aquarius-A1(=), A2(=), A3(=), and A4(=) resolution levels, respectively. The evolution of the subhalo mass $M_{mx}$ is shown as a function of time t. Numbered circles @, ..., @ mark remnant masses of $log_{10} M_{mx}/M_{mx0} = 0, -1, ..., -6$. Shaded regions correspond to initial conditions spanning a range of circular velocities $20 \le V_{mx0}/\text{km s}^{-1} \le 40$. Insufficient numerical resolution causes the N-body models disrupt artificially, whereas the empirical model predicts an evolution toward a stable asymptotic remnant state.</figcaption></figure>

density profile with a scale radius $r_s$ and a scale density $\rho_s$,

$$\rho_{\text{NFW}}(r) = \rho_s (r/r_s)^{-1} (1 + r/r_s)^{-2}. \quad (1)$$

Instead of referring to $r_s$ and $\rho_s$ directly, we characterize the subhalo using two equivalent parameters: the peak velocity $V_{mx} \approx 1.65 r_s \sqrt{G\rho_s}$ of the subhalo circular velocity profile $V_c = \sqrt{GM(<r)/r}$, and the radius $r_{mx} \approx 2.16 r_s$ where the peak velocity is reached. The effect of tides on the subhalo is modeled as an exponential truncation and a renormalization of the density profile (see EN21, Equation (7)),

$$\rho(r) = \rho_{\text{NFW}}(r) \times \exp(-r/r_{\text{cut}})/(1 + r_s/r_{\text{cut}})^{\kappa}, \quad (2)$$

with $\kappa \approx 0.3$ obtained from fits to N-body simulations. For $r_{cut}/r_s \to \infty$, the above equation reduces to the initial NFW profile, whereas for strong tidal truncation, $r_{cut}/r_s \to 0$, the profile converges to an exponentially truncated cusp,

$$\rho_{\text{asy}}(r) = \rho_{\text{cut}} \exp(-r/r_{\text{cut}}) (r/r_{\text{cut}})^{-1}. \quad (3)$$

For the exponentially truncated cusp, $r_{mx} \approx 1.79 r_{cut}$, $V_{mx} \approx 1.94 r_{cut} \sqrt{G\rho_{cut}}$. The truncated cusp has a total mass of $M_{cut} = 4\pi\rho_{cut}r_{cut}^3$, roughly twice as large as the mass $M_{mx} \approx 0.54 M_{cut}$ enclosed within $r_{mx}$. As tides strip the subhalo, its characteristic velocity $V_{mx}$ and size $r_{mx}$ decrease, following a tidal track (see Peñarrubia et al. 2008) that couples the evolution of $V_{mx}$ to that of $r_{mx}$ relative to their initial values (here, we use the track of EN21; see their Equation (5)):

$$V_{\text{mx}}/V_{\text{mx0}} = 2^{\alpha} (r_{\text{mx}}/r_{\text{mx0}})^{\beta} [1 + (r_{\text{mx}}/r_{\text{mx0}})^2]^{-\alpha}, \quad (4)$$

with $\alpha \approx 0.4$ and $\beta \approx 0.65$ obtained from fits to N-body simulations.

We also use the EN21 model to estimate the time it takes for a subhalo to be stripped to a given remnant mass, assuming that the subhalo is on an orbit with a constant pericenter $r_{peri}$ and

apocenter $r_{apo}$ within an isothermal host potential,

$$\Phi_{\text{host}} = V_{\text{host}}^2 \ln(r/r_0), \quad (5)$$

where by $V_{host}$ we denote the host's constant circular velocity, and $r_0$ is an arbitrary scale radius.

The rate of tidal stripping depends on the density contrast between subhalo and host, measured at the pericenter (EN21 Equations (12) and (15)). For a given pericenter, the only effect of orbital eccentricity is to delay tidal evolution with respect to the evolution on a circular orbit. For subhaloes that are underdense with respect to the mean host density at the pericenter, the tidal evolution gradually decelerates until a final remnant state is reached where the subhalo mean density within $r_{mx}$ is $\approx 16$ times higher than the mean host density at the pericenter. In terms of timescales, heavily stripped subhaloes converge toward a characteristic circular time $T_{mx} = 2\pi r_{mx}/V_{mx} \approx T_{peri}/4$ set by the circular time of the host at the pericenter $T_{peri} = 2\pi r_{peri}/V_{host}$.

As an example, we use the empirical method outlined here to model the tidal evolution of a dark matter subhalo with the same characteristic mass and size as in the N-body example of Section 2.2.3 (see Table 1), and place it on the same orbit. The time evolution of the dark matter subhalo is shown in Figure 2: for the range of remnant masses resolved in the N-body simulations, there is excellent agreement between the empirical model and the simulation results. During $\sim 13$ Gyr of evolution, the subhalo is stripped to a remnant mass $\sim 10^6$ times smaller than its initial mass.

## 3.2. Model for the Evolution of Stellar Tracers

We now turn our attention to modeling the effect of tides on a stellar system embedded within a dark matter subhalo. We use the model outlined and tested in E+22, where stars are assumed to be massless and collisionless tracers within the subhalo potential. The effect of tides on a dwarf galaxy can then be described through subsequent truncations in the energy distributions of dark matter and stars.

Under the assumptions of spherical symmetry and an isotropic velocity dispersion, for a given subhalo potential $\Phi(r)$, the initial stellar component is fully defined by its energy distribution

$$dN_{\star}/d\mathcal{E}|_{i} = (4\pi)^{2} f_{\star}(\mathcal{E}) p(\mathcal{E}), \quad (6)$$

where by $f_*(\mathcal{E})$ we denote the initial stellar distribution function, and by $p(\mathcal{E})$, we denote the density of states in the initial NFW potential. Both are functions of energy $\mathcal{E} = 1 - E/\Phi_0$, where $E = v^2/2 + \Phi(r)$, and $\Phi_0$ denotes the potential minimum. For a given initial stellar density profile $\rho_*(r)$, the stellar distribution function can be obtained by Eddington inversion. We compute $f_*(\mathcal{E})$ and $p(\mathcal{E})$ using the implementation described in EP20, available online (see footnote 7). We will discuss our choices of initial stellar density profiles separately in Section 3.3, and now proceed assuming that an appropriate initial $dN_*/d\mathcal{E}$ has been defined.

The effect of tides on the stellar energy distribution $dN_*/d\mathcal{E}$ can be approximated by a tapered truncation. We model the truncation empirically through (E+22, Equation (9)),

$$dN_{\star}/d\mathcal{E}|_{i,t} = \frac{dN_{\star}/d\mathcal{E}|_{i}}{1 + (a\mathcal{E}/\mathcal{E}_{\text{mx},t})^{b}}, \quad (7)$$
<figure><img src="image_6.png" /><figcaption>Figure 3. Solid curves show the initial energy distribution $dN/d\mathcal{E}$ of dark matter (gray —) and stars (red—), corresponding to snapshot @ shown in Figure 4. Dashed curves show the tidally truncated energy distributions of dark matter (—) and stars (—) in the initial conditions for remnant subhalo masses of $\log_{10} M_{\text{mx}}/M_{\text{mx0}} = 0, -1, \dots, -5$ (snapshots ①, \dots, ⑤). The stars initially follow a 2D exponential surface brightness profile (model exp2D, see Section 3.3.1) with a half-light radius of $R_{h0} = r_{\text{mx0}}/16$. The stellar energy distribution is normalized so that $\max(dN/d\mathcal{E}) = 1$; the normalization for the dark matter energy distribution is arbitrary.</figcaption></figure>

where by $dN_*/d\mathcal{E}|_{l,t}$ we denote the truncated energy distribution in the initial conditions, and $\mathcal{E}_{\text{mx},t}$ is the energy scale beyond which the energy distribution is truncated. The parameters $a \approx 0.85$, $b \approx 12$ are obtained through fits to N-body simulations.

As an example, Figure 3 shows the initial stellar and dark matter energy distributions for a stellar tracer with a 2D exponential surface brightness profile, deeply embedded in an NFW halo ($R_h/r_{\text{mx}} = 1/16$). Energies $\mathcal{E} = 1 - E/\Phi_0$ are measured relative to the potential minimum $\Phi_0$. In this notation, the most-bound energy state (the ground state) is $\mathcal{E} = 0$, and the boundary between bound and unbound particles lies at $\mathcal{E} = 1$.

The energy-truncated system is initially out of equilibrium. We model the return to virial equilibrium empirically using fits to N-body simulations. The virialization process preserves (on average) the order of energies: the most-bound particles prior to virialization are also (on average) the most-bound particles after virialization. The mapping of initial virialization energies $\mathcal{E}$ to final energies $\mathcal{E}_f$ after virialization follows on average the empirical relation (E+22, Equation (12))

$$ \bar{\mathcal{E}}_f = [1 + (c\mathcal{E}/\mathcal{E}_{\text{mx},t})^d]^{1/d}, \quad (8) $$

with $c \approx 0.8$ and $d \approx -3$. Note that the virialization energies $\mathcal{E}$ are defined in the initial NFW potential, whereas the energies after virialization $\mathcal{E}_f$ are defined in the tidally stripped subhalo potential.

Under the assumption of isotropic velocity dispersion and spherical symmetry, we then reconstruct the distribution function from the virialized energy distribution in the evolved subhalo potential. All observable properties of the stellar component can be derived from the evolved distribution function.

The simple empirical method outlined here allows us to model the tidal evolution of accreted substructures down to tiny remnant masses and sizes. In Figure 4, we apply the empirical tidal stripping model to follow the evolution of a dwarf galaxy with a (2D) exponential surface brightness profile with an

initial projected half-light radius of $R_h \approx 440$ pc, and an initial luminosity of $L \approx 10^6 L_\odot$. We embed the dwarf galaxy in a subhalo with the same characteristic mass and size as in the N-body example of Section 2.2.3 (see Table 1). This results in a line-of-sight velocity dispersion of $(\sigma_{\text{los}}^2)^{1/2} \approx 12 \text{ km s}^{-1}$. The luminosity, size, and velocity dispersion of the dwarf galaxy model are therefore roughly comparable to those of classical Milky Way satellites like Ursa Minor, Sextans, or Sculptor. We model the tidal evolution of the dwarf galaxy assuming the same orbit as in Section 2.2.3.

After ~6 Gyr of tidal evolution, snapshot ⑤, the dwarf galaxy has been stripped to the size of a microgalaxy, with an enclosed dark matter mass of $\sim 10^4 M_\odot$, a luminosity of only $\sim 50 L_\odot$, and a size of several tens of parsecs, well beyond the resolution limits of N-body simulations like Aquarius-A1 (see Figure 2). After ~12 Gyr of evolution, snapshot ⑥, the stellar system has been stripped to roughly solar luminosity, but still encloses $\sim 10^3 M_\odot$ of dark matter.

### 3.3. Stellar Component Models

Tidally stripped dwarf galaxies retain some memory of the properties of their initial stellar distribution E+22. We aim to explore how different initial stellar distributions affect the tidal evolution of the remnant luminosity $L$, projected half-light radius $R_h$, and line-of-sight velocity dispersion $\sigma_{\text{los}}$. In Sections 3.3.1, 3.3.2, and 3.3.3, we define the three different stellar density profiles that we adopt in our modeling. In Section 3.3.4, we briefly discuss the similarities and differences between these three models. Finally, in Section 3.4, we describe the time evolution of the surface brightness and line-of-sight velocity dispersion profiles during tidal stripping. All stellar models discussed here provide a plausible description of the observed surface brightness profiles of faint Milky Way satellites (see Appendix C for a comparison against available observational data).

#### 3.3.1. 2D exponential—exp2D

The surface brightness profiles of many Milky Way dSph galaxies are well approximated by 2D exponential profiles (see, e.g., Irwin & Hatzidimitriou 1995 for fits to the Carina, Draco, Fornax, and Leo I dwarf galaxies; or more recently Wang et al. 2019 for Fornax; Sestito et al. 2023 for Ursa Minor; and Jensen et al. 2024 for multiple Milky Way satellites including Sculptor and several ultrafaint systems). Motivated by this observational evidence, we adopt a 2D exponential as “fiducial” profile for our models, parameterized through

$$ \Sigma_*(R) = \Sigma_0 \exp(-R/R_*) \quad (9) $$

where $\Sigma_0 = L/(2\pi R_*^2)$, for a (2D) half-light radius of $R_h \approx 1.68 R_*$. Under the assumption of spherical symmetry, we use the inverse Abel transform to reconstruct the corresponding 3D stellar density profile (see, e.g., Binney & Tremaine 1987, Appendix 1B(4)),

$$ \rho_*(r) = -\frac{1}{\pi} \int_r^\infty \frac{d\Sigma_*(R)/dR}{\sqrt{R^2 - r^2}} dR = \rho_{s*} K_0(r/R_*), \quad (10) $$

where $K_0$ denotes the (zeroth-order) modified Bessel function of second type, and $\rho_{s*} = \Sigma_0/(\pi R_*)$, for a (3D) half-light radius of $r_h \approx 2.22 R_*$. Asymptotically, for $r \to 0$ (see Abramowitz &
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Tidal evolution of a dwarf galaxy embedded in a dark matter subhalo, computed using the empirical energy-truncation model of E+22. The dark matter surface-density $\rho_{2D}$ is shown in gray, and the surface-brightness $\Sigma_{*}$ of the embedded dwarf is shown in color. Snapshots, labeled @, ..., @, are shown for fixed subhalo remnant masses of $\log_{10} M_{\text{mx}}/M_{\text{mx0}} = 0, -1, \dots, -6$. The scale is kept constant between the first two panels and the smaller top left panels. The stars initially follow a (2D) exponential surface brightness profile with half-light radius $R_{h0} = r_{\text{mx0}}/16$ (model exp2D, see Section 3.3.1). Initially, the stars are hardly affected by tides, and mainly dark matter is lost (first three panels). Once the tidal energy truncation reaches the stellar component, the characteristic size of the dark matter halo $r_{\text{mx}}$ (shown as black circles) and the stellar half-light radius $R_{h}$ (red circles) evolve in parallel (last three panels).</figcaption></figure>

Stegun 1972, Equation (9.6.8)),

$$ \rho_{*}(r) \sim \rho_{s*} \ln(R_{*}/r), \qquad (11) $$

i.e., the 3D density distribution that generates a 2D exponential surface brightness profile diverges logarithmically at the center.

### 3.3.2. 3D Exponential—exp3D

We furthermore consider a stellar tracer that follows a 3D exponential density profile,

$$ \rho_{*}(r) = \rho_{0} \exp(-r/r_{*}), \qquad (12) $$

where $\rho_{0} = L/(8\pi r_{*}^{3})$, for a 3D half-light radius of $r_{h} \approx 2.67 r_{*}$. The 2D projection of this profile gives

$$ \Sigma_{*}(R) = 2 \int_{R}^{\infty} \frac{\rho_{*}(r) r}{\sqrt{r^{2} - R^{2}}} dr = \Sigma_{0} \frac{R}{r_{*}} K_{1}(R/r_{*}), \qquad (13) $$

where $K_{1}$ denotes the (first-order) modified Bessel function of second type, and $\Sigma_{0} = 2\rho_{0}r_{*}$, for a 2D half-light radius of $R_{h} \approx 2.03 r_{*}$ (see Abramowitz & Stegun 1972, Equation (9.6.9)).

### 3.3.3. 2D Logarithmic Cusp—log2D

Finally, we consider a stellar tracer, which, in 3D, diverges at the center like an NFW dark matter profile, $\rho_{*}(r) \sim r^{-1}$ for $r \ll r_{*}$. For this model, we choose a stellar density profile formally identical to Equation (3),

$$ \rho_{*}(r) = \rho_{s*} \frac{r_{*}}{r} \exp(-r/r_{*}), \qquad (14) $$

where $\rho_{s*} = L/(4\pi r_{*}^{3})$, with a 3D half-light radius of $r_{h} \approx 1.68 r_{*}$. The 2D projection gives

$$ \Sigma_{*}(R) = \Sigma_{s} K_{0}(R/r_{*}), \qquad (15) $$

where $K_{0}$ denotes the (zeroth-order) modified Bessel function of second type, and $\Sigma_{s} = 2\rho_{s*}r_{*}$, for a 2D half-light radius of $R_{h} \approx 1.26 r_{*}$. For $R \to 0$,

$$ \Sigma_{*}(R) \sim \Sigma_{s} \ln(r_{*}/R), \qquad (16) $$

i.e., the (2D) surface brightness profile of an exponentially truncated cuspy stellar (3D) density profile logarithmically diverges at the center (see Abramowitz & Stegun 1972, Equation (9.6.8)).

<caption>**Table 2**</caption>

Central Asymptotics for the 3D Stellar Density Profiles $\rho_{*}(r)$ and the 2D Surface Brightness Profiles $\Sigma_{*}(R)$ Corresponding to the exp2D, exp3D, and log2D Models (for Definitions, See Sections 3.3.1, 3.3.2, and 3.3.3, Respectively)

<table><thead><tr><th>Model</th><th>Equation</th><th>3D Asymptotics<br>r &lt;&lt; r<sub>h</sub></th><th>2D Asymptotics<br>R &lt;&lt; R<sub>h</sub></th></tr></thead><tbody><tr><td>exp2D</td><td>(9), (10)</td><td>ρ<sub>*</sub>(r) ~ ln(R<sub>*</sub>/r)</td><td>Σ<sub>*</sub>(R) → const</td></tr><tr><td>exp3D</td><td>(12), (13)</td><td>ρ<sub>*</sub>(r) → const</td><td>Σ<sub>*</sub>(R) → const</td></tr><tr><td>log2D</td><td>(14), (15)</td><td>ρ<sub>*</sub>(r) ~ 1/r</td><td>Σ<sub>*</sub>(R) ~ ln(r<sub>*</sub>/R)</td></tr></tbody></table>

### 3.3.4. Comparison of Initial Profiles

Here, we discuss the observable differences between the exp2D, exp3D, and log2D stellar models defined in the previous three subsections. For reference, Table 2 summarizes the central asymptotic behavior of the three stellar models.

In projection, the exp2D, exp3D models are cored, i.e., $d \ln \Sigma_{*}/d \ln R \to 0$ for $R \to 0$. In contrast, the surface brightness profile of the log2D model has a very shallow, logarithmically diverging central cusp. The surface brightness profiles are shown in the top right panel of Figure 5, scaled to a half-light radius of $R_{h} = 440$ pc.

The exp2D model, shown in red, is identical to the example stellar tracer shown previously in Figures 3 and 4. (Projected) 3D stellar density profiles are shown in the bottom left panel of Figure 5.

In 3D, only the exp3D model has a centrally finite density: the exp2D diverges logarithmically, and the log2D diverges with a power-law slope of $d \ln \rho_{*}/d \ln r \to -1$ for $r \to 0$, like an NFW profile. Consequently, the stellar energy distribution of the log2D model has the same inner asymptotics as the dark matter when embedded in an NFW halo. This is shown in the top left panel of Figure 5 for a stellar tracer deeply embedded in an NFW halo ($R_{h}/r_{\text{mx}} = 1/16$). Note that the energy distributions of the exp2D and exp3D models drop faster toward the most-bound state than the log2D model.

Finally, in the bottom right panel of Figure 5, we show the line-of-sight velocity dispersion profiles of the three different stellar models, embedded in an NFW (subhalo) potential with parameters as in Table 1. The luminosity-averaged line-of-sight velocity dispersion,

$$ \langle \sigma_{\text{los}}^{2} \rangle = \frac{2\pi}{L} \int_{0}^{\infty} \Sigma_{*}(R) \sigma_{\text{los}}^{2}(R) R \, dR, \qquad (17) $$

$$ = \frac{4\pi G}{3} \int_{0}^{\infty} r \nu_{*}(r) M(<r) dr, \qquad (18) $$
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Top left: Initial energy distributions for the exp2D (—), exp3D (—), and log2D (—) stellar models (Sections 3.3.1, 3.3.2, 3.3.3) embedded in an NFW dark matter halo with parameters as in Table 1. All three models share the same projected half-light radius of $R_{h0}/r_{mx0} = 1/16$. For comparison, the energy distribution of the underlying NFW halo is shown as a gray dashed curve (—). Bottom left: corresponding 3D stellar density profiles $\rho_*(r)$. Top right: 2D stellar surface brightness profiles $\Sigma_*(R)$. Bottom right: Line-of-sight velocity dispersion profiles $\sigma_{los}(R)$. The luminosity-averaged line-of-sight velocity dispersion (Equation (17)) is virtually identical between the models and is shown as a gray band.</figcaption></figure>

shown as a gray-shaded band, is nearly identical for the three models. We denote by $\nu_*(r)$ the (3D) stellar tracer density, normalized so that $4\pi \int_0^\infty r^2 \nu_*(r) dr = 1$, and $L = 2\pi \int_0^\infty R \Sigma_*(R) dR$ is the total luminosity. $M(<r)$ is the total mass enclosed within radius $r$.

Equation (18) follows from the projected virial theorem (see, e.g., Amorisco & Evans 2012; Errani et al. 2018), and guarantees that $\langle \sigma_{los}^2 \rangle$ is independent of anisotropies in the velocity dispersion. Note that a stochastic realization of the luminosity-averaged $\langle \sigma_{los}^2 \rangle$ is the only kinematic observable accessible for faint stellar systems that contain only a few stars sufficiently bright to accurately measure their radial velocities, and that this estimate may be artificially inflated if some of these stars are in binary systems. Distinguishing between these three models on kinematic grounds hence poses a considerable observational challenge.

### 3.4. Tidal Evolution of Stellar Components

As tides strip a dwarf galaxy, its surface brightness and velocity dispersion profiles evolve. For a given subhalo and orbit, the evolution depends on the shape of the initial stellar density profile, and on how deeply embedded the stellar profile is within its surrounding dark matter subhalo.

The top panel of Figure 6 shows the tidal evolution of the surface brightness profiles for stellar tracers with initial conditions as in Figure 5. Profiles are shown for subhalo remnant masses of $\log_{10} M_{\text{mx}}/M_{\text{mx0}} = 0, -2, -4, -6$, as labeled in the left-hand panel. The location of the (2D) half-light radius is indicated by the label location in the left-hand panel, and by filled circles in the central and right-hand panel.

For all three models shown, the surface brightness profile is hardly affected by tides in the early stages of tidal evolution. Once the underlying subhalo has lost ~99% of its initial characteristic mass, the tidal energy truncation reaches the stellar component (snapshot ②, compare with Figure 3). From that point on, the surface brightness of the stellar component decreases with each further tidal energy truncation, and the half-light radius of the stellar component decreases.

For the exp2D stellar model (left-hand panel), this evolution is illustrated also in Figure 4. The surface brightness of the stellar component is shown in color, whereas the dark matter surface density is shown in gray. After a slight initial expansion of the stellar half-light radius $R_h$ (shown as a red circle), once the tidal energy truncation has reached the stellar component, dark matter and stars are truncated approximately at the same radius, and the stellar half-light radius $R_h$ evolves thereafter in sync (Kravtsov 2010) with the characteristic size $r_{mx}$ (black circle) of the underlying dark matter subhalo.

Returning to Figure 6, we see that the surface brightness of the exp2D (red) and exp3D (orange) models drops much more rapidly during tidal stripping than the surface brightness of the log2D (purple) model. This disparate evolution is easily understood when looking at the underlying stellar energy distributions (Figure 5).

As tides truncate the dwarf galaxy in energy, the log2D model loses stars and dark matter at a similar pace, as its stellar energy distribution follows that of the dark matter toward the most-bound energy states. On the other hand, the most-bound energy states are hardly populated in the exp3D stellar model: its stellar energy distribution drops much more rapidly toward the most-bound states than that of the dark matter, making the stellar component vulnerable to full tidal disruption. The
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Top panels: Time evolution of the surface brightness profiles of the exp2D (—), exp3D (—), and log2D (—) stellar models, with initial conditions as shown in Figure 5. The exp2D tracer is identical to the model of Figure 4, with selected snapshots, labeled @; ..., @, shown for fixed subhalo remnant masses of $log_{10} M_{\text{mx}}/M_{\text{mx0}} = 0, -2, -4, -6$. The label location coincides with the projected half-light radius. For the exp3D and log2D models, the half-light radius is marked using a filled circle. Note that the time evolution of the surface brightness profiles differs substantially between the three stellar models shown. Bottom panels: Time evolution of the corresponding line-of-sight velocity dispersion profiles. Note that, in the tidally limited regime, the velocity dispersion profile drops steeply beyond $R_h$.</figcaption></figure>

prospect of observing microgalaxies hence depends critically on how stars are distributed energetically within the dark matter subhalo.

The bottom panel of Figure 6 shows the time evolution of the corresponding line-of-sight velocity dispersion profiles. The average velocity dispersion of all three stellar models drops monotonously with tidal mass loss: the evolution of the line-of-sight velocity dispersion depends only weakly on the shape of the initial stellar profile.

## 4. Application to Local Group Satellites

We are now ready to apply the model outlined in Section 3 to predict the luminosity, structure, and kinematics of heavily stripped microgalaxies, and to compare predicted properties with those of observed stellar systems in the Local Group.

In particular, we aim to address the nature of recently discovered stellar systems with half-light radii $1 \lesssim R_h/\text{pc} \lesssim 30$ and luminosities of $10 \lesssim L/L_\odot \lesssim 10^3$ (Torrealba et al. 2019; Mau et al. 2020; Cerny et al. 2023, 2023; Smith et al. 2024). These systems populate a parameter space at the boundary between the globular cluster and dwarf galaxy regimes. If gravitationally bound, then, in principle, they could be one of the following three possibilities:

(i) self-gravitating star clusters devoid of dark matter, like globular clusters,

(ii) dark matter-dominated galaxies so deeply embedded within their dark matter haloes that their sizes and luminosities have not been affected by tides since formation, or

(iii) dark matter-dominated microgalaxies of tidal origin, formed through stripping of larger and more luminous progenitors.

We will in the following discuss these potential formation scenarios guided by the tidal evolutionary tracks computed for the exp2D, exp3D, and log2D stellar models.

### 4.1. Half-light Radii and Luminosities

Figure 7(a) shows a compilation<fnref n="9" /> of projected half-light radii, $R_h$, and luminosities, $L$, of Local Group dwarf galaxies (blue circles) and globular clusters (yellow triangles). On average, at equal luminosity, globular clusters are significantly more compact than dwarf galaxies. For example, at $L = 10^6 L_\odot$, globular clusters span a range in sizes of roughly $1 \sim 10$ pc, more than an order of magnitude more compact than dwarfs of the same luminosity, $0.1 \sim 1$ kpc.

This clear separation in size becomes ambiguous at lower luminosities. Ultrafaint dwarfs like Bootes 2, Carina 3, and Draco 2 have half-light radii of 20–30 pc, not too different from extended Milky Way globular clusters like Palomar 5 or

<footnote n="9">⁹ The dwarf galaxy properties are taken from McConnachie (2012; version 2021 January, with updates for Antlia 2, Ji et al. 2021; Bootes 2, Bruce et al. 2023; Crater 2, Ji et al. 2021; Draco 2, Martin et al. 2016; Longeard et al. 2018; Tucana, Taibi et al. 2020; Tucana 2, Chiti et al. 2021; And 19, Collins et al. 2020; and And 21, Collins et al. 2021). For globular clusters, the data shown is from Harris (1996; version 2010 December; with updated half-light radii and velocity dispersions for Pal-5 from Kuzma et al. 2015; Gieles et al. 2021; for NGC 2419 from Baumgardt et al. 2009, and for Pal-14 from Hilker 2006; Jordi et al. 2009). Objects marked as “unidentified” at the boundary of the globular cluster and dwarf galaxy regimes are as compiled in Cerny et al. (2023, 2023). Data for Ursa Major 3/Unions 1 are from Smith et al. (2024).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a) Luminosities *L* and half-light radii *R*<sub>h</sub> of Local Group dwarf galaxies (●) and globular clusters (♥). Ambiguous objects at the boundary between the dwarf and globular cluster regime are marked using gray crosses (x). For references, see footnote 9. Solid curves show tidal evolutionary tracks for the exp2D (—), exp3D (■), and log2D (—) stellar models (with initial conditions as in Figure 5). Snapshots corresponding to the exp2D model, labeled ●, ..., ⊕, are highlighted for fixed subhalo remnant masses of log<sub>10</sub> *M*<sub>max</sub>/*M*<sub>max0</sub> = 0, -1, ..., -5. (b) Line-of-sight velocity dispersion (*σ*<sub>los</sub><sup>2</sup>)<sup>1/2</sup> (Equation (17)) and projected half-light radius *R*<sub>h</sub>. The top black curve labeled “initial” shows the line-of-sight velocity dispersion expected for a stellar tracer of half-light radius *R*<sub>h</sub> embedded in a ~10<sup>9</sup> *M*<sub>☉</sub> dark matter halo. Yellow crosses (■) show velocity dispersions computed under the assumption of self-gravity for the ambiguous objects shown as gray crosses in the top left panel. If instead these objects are dark matter-supported galaxies with 1 ≲ *R*<sub>h</sub>/pc ≲ 30, their dispersion should fall in the region highlighted in blue (—). The region is computed assuming a dwarf galaxy progenitor embedded in a subhalo with a halo mass close to the *z* = 0 hydrogen cooling limit. For earlier formation redshifts, the lower bound shifts upwards; see Appendix B. (c) Luminosities *L* and dynamical masses *M*<sub>1.8</sub> (Equation (21)). Diagonal lines mark constant dynamical mass-to-light ratios of *M*<sub>1.8</sub>/*L* = 1 *M*<sub>☉</sub>/*L*<sub>☉</sub> and 100 *M*<sub>☉</sub>/*L*<sub>☉</sub>. (d) Projected half-light radii *R*<sub>h</sub> and mean densities *ρ*<sub>1.8</sub> = *ρ*(<1.8 *R*<sub>h</sub>) (Eq. 22) enclosed within a spherical radius of 1.8 *R*<sub>h</sub>.</figcaption></figure>

Palomar 14 (albeit at lower luminosity). Recently discovered
objects, whose nature is still unclear, are labeled “unidentified”
and depicted as gray crosses in Figure 7(a). They populate a
region of parameter space at the boundary between the globular
cluster and dwarf galaxy regimes. As an example, we
show error bars for the faint Milky Way satellite Ursa
Major 3/Unions 1 (UMa3/U1 for short), with a half-light

radius of only (3 ± 1)pc, a total stellar mass of $M_* = 16_{-5}^{+6} M_\odot$, corresponding to a luminosity of $M_V \approx +2.2_{-0.3}^{+0.4}$ mag (Smith et al. 2024).

To address whether these ambiguous objects could be the
remnants of tidally stripped dwarf galaxies, we show the
evolution of three dwarf galaxy models in Figure 7(a),
computed using the method outlined in Section 3. As an
example, we assume that the progenitor dwarf galaxy is initially embedded in a dark matter subhalo with structural properties as listed in Table 1, with a virial mass close to the redshift $z=0$ hydrogen cooling limit. The systematics arising from earlier formation redshifts and lower halo masses are studied in Appendix B. The example progenitor considered here has an initial luminosity of $L = 10^6 L_\odot$, a half-light radius of $R_h = 440$ pc, and a velocity dispersion $\sigma_{los} \approx 12$ km s$^{-1}$, chosen to roughly resemble classical Milky Way satellites like Ursa Minor, Sextans, or Sculptor.

The evolution of size and luminosity for a dwarf galaxy with an exp2D surface brightness profile is shown in red (identical to the model shown in Figures 3 and 4, with circled numbers @, ..., @ corresponding to fixed subhalo remnant masses of $log_{10} M_{mx}/M_{mx0} = 0, -1, ..., -6$). As tides remove stars from the exp2D model, its size and luminosity drop. The evolution moves the dwarf galaxy model roughly along the observed mass–size relation of dwarf galaxies. The evolution of a dwarf galaxy with an initial half-light radius and luminosity different from the example model shown in Figure 7(a) can be roughly estimated by shifting the tidal evolutionary track within the size–luminosity plane.

As discussed in Section 3.3, the evolution of size and luminosity crucially depends on the shape of the initial stellar density profile, which in turn depends on the degree to which stellar binding energies extend to the most-bound energy states in the subhalo. The evolution of the exp3D model is shown in orange. For this model, the luminosity drops faster with radius than in the case of the exp3D model.

On the other hand, for the log2D model (shown in purple), the drop in luminosity is less steep$^{10}$: tidal evolution shifts the dwarf galaxy model roughly along the observed mass–size relation of dwarf galaxies, and, for the example initial conditions chosen here, the stripped model roughly approaches the region of parameter space occupied by the unidentified systems.

It is clear that the majority of the unidentified systems depicted as gray crosses in Figure 7(a) cannot be reached by tidally stripping progenitor systems with initial luminosities and sizes comparable to those assumed in the example above. A progenitor whose tidal descendent is consistent with the unidentified systems can be found by shifting the tidal tracks upwards and to the left (Leo 1 could be one example). We may conclude that, if the unidentified systems are indeed microgalaxies, then their progenitors must have been systems of higher luminosity and/or surface brightness than the typical dSphs, and their stellar density must be closer to that of the log2D model than exp3D.

The tidal tracks shown in Figure 7(a) highlight the fact that tidal evolution in the size/luminosity plane is very sensitive to the underlying distribution of stellar binding energies. Further constraints are needed to either rule out or to confirm conclusively the possibility that many of the unidentified systems are indeed microgalaxies.

## 4.2. Half-light Radii and Velocity Dispersions

We now discuss to what extent measurements of the line-of-sight velocity dispersion $\langle\sigma_{los}^2\rangle^{1/2}$ (Equation (17)) can help to constrain the nature of the unidentified stellar systems. Figure 7(b) shows the half-light radii $R_h$ and line-of-sight velocity dispersions $\langle\sigma_{los}^2\rangle^{1/2}$ for the same dwarf galaxies and globular clusters as in Figure 7(a). Dwarf galaxies, on average, are more extended than globular clusters at equal velocity dispersion.

### 4.2.1. Self-gravitating Star Clusters

For most of the objects marked as “unidentified” in Figure 7(a), no robust measurement of velocity dispersion is available to date. If these objects were self-gravitating and devoid of dark matter, their velocity dispersion could be estimated from their half-light radius $R_h$ and their total stellar mass $M_*$. For a 3D exponential stellar density profile (for definition, see Section 3.3.2) with $R_h \approx 2.03 r_*$ and total stellar mass $M_*$, Equation (18) gives

$$ \langle\sigma_{\text{los}}^2\rangle = (5/96) GM_* / r_* \qquad (19) $$

For each “unidentified” system shown in Figure 7(a), we use Equation (19) to estimate the (virial) velocity dispersion the object would have in the absence of dark matter. For this calculation, we assume a mean stellar mass-to-light ratio of $\Upsilon \approx 1.6$ (Woo et al. 2008). The results of this exercise$^{11}$ are marked as yellow crosses in Figure 7(b). For reference, dashed diagonal lines in Figure 7(b) are also computed from Equation (19), for constant stellar masses of $M_* = \Upsilon L = 10^2$, $10^3$, and $10^4 M_\odot$.

### 4.2.2. Dark Matter-dominated Objects

For the case of dark matter-dominated dwarf galaxies, the gravitational potential is sourced primarily by the dark matter component. Their stellar velocity dispersions depend mainly on the mass of the halo, and the degree of spatial segregation between stars and dark matter.

In the LCDM cosmology, dwarf galaxies are expected to form within a relatively narrow range of halo masses (see, e.g., Guo et al. 2010, 2011; Fattahi et al. 2018). An estimate of the characteristic minimum halo mass may be obtained from the constraint that, to form stars, gas must be able cool in presence of the cosmic UV background. Only haloes with masses above the hydrogen cooling limit provide the necessary conditions for the onset of star formation.

At redshift $z=0$, this limiting virial mass is roughly $\sim 5 \times 10^9 M_\odot$ (see, e.g., Benitez-Llambay & Frenk 2020). The velocity dispersion of an exponential stellar profile embedded in an average-concentration NFW halo with a halo mass close to the $z=0$ hydrogen cooling limit ($r_{mx0} = 7$ kpc, $V_{mx0} = 35$ km s$^{-1}$, see Table 1) is shown as a black solid curve in Figure 7(b), labeled “initial.” This curve is computed from Equation (18), with the integral evaluated numerically. We

<footnote n="10">¹⁰ Some simple analytical insight in the disparate evolution of stellar models with exp2D, exp3D, and log2D density profiles can be gained from the tidally limited regime, when the stellar component and the subhalo are trimmed down to similar sizes and evolve in sync ($R_h \propto r_{mx}$ and $M_h \propto M_{mx}$). Calling $\alpha_*$ the slope of the stellar energy distribution $dN_*/dE \propto E^{\alpha_*}$ for $E \ll 1$, then $L \propto M_h^{(\alpha_*+1)/2}$ (assuming NFW haloes, where $dM/dE \propto E$ for $E \ll 1$). Combining this with the tidal track of Equation (4) with $M_h \propto R_h^{2.9+1}$ gives $d \ln L/d \ln R_h = (2\beta+1)(\alpha_*+1)/2$. Hence, the evolution of luminosity and size of a stellar tracer component is directly related to how stars populate the most-bound energy states of the subhalo, parameterized through the log-slope $d \ln(dN/dE)/d \ln E = \alpha_*$ for $E \ll 1$.</footnote>

<footnote n="11">¹¹ The velocity dispersion estimates for self-gravitating clusters are fairly insensitive to the choice of the underlying density profile: for Plummer spheres, we find $\langle\sigma_{los}^2\rangle = (\pi/32) GM_* / R_h$, whereas before by $M_*$ we denote the total stellar mass, and by $R_h$ the projected half-light radius. For a given $M_*$ and $R_h$, the dispersion $\langle\sigma_{los}^2\rangle^{1/2}$ differs by $\sim 4\%$ between a Plummer sphere and a 3D exponential profile.</footnote>
discuss the redshift-dependence of these initial conditions in Appendix B.

More generally, cosmological simulations suggest that satellites of Milky Way–like host galaxies formed within dark matter haloes of the peak circular velocity $20 \lesssim V_{\text{mx}}/\text{km s}^{-1} \lesssim 40$ (e.g., Fattahi et al. 2018). This range of potential initial haloes is indicated by the gray band in Figure 7(b), which includes a $\pm 0.15$ dex variation in the concentration around the mean. Note that many dSphs of the Local Group fall within this band, which suggests that they have experienced only minor tidal perturbations.

The present-day sizes and dispersions of dwarf galaxies and their remnants depend on their properties at formation, and on their tidal history. In particular, dwarf galaxies that experienced little tidal mass loss would be located within the gray band of Figure 7(b).

Tidal stripping reduces the stellar velocity dispersions, providing a potential explanation as to why some observed dwarfs are found below the gray band. To illustrate this, we show the tidal evolution of half-light radius $R_h$ and velocity dispersion $\langle\sigma_{\text{los}}^2\rangle^{1/2}$ of the three example dwarf galaxy models with exp2D, exp3D, and log2D stellar tracers using red, orange, and purple curves in Figure 7(b). Both models are initially located at the position marked ⊙, and are embedded in the same dark matter halo. As previous, circled numbers ⊙, ..., ⊙ correspond to the snapshots shown in Figure 4 at fixed subhalo remnant masses of $\log_{10} M_{\text{mx}}/M_{\text{mx}0} = 0, -1, \dots, -6$.

In the early stages of tidal evolution, the tidal energy truncation does not reach the binding energies of the stars (compare with Figure 3). Dark matter on weakly-bound orbits is removed by tides. Some of these orbits do contribute to the density within the stellar radii: after their removal, the velocity dispersion of an embedded stellar tracer drops. This causes the initial near-vertical evolution in the size–velocity dispersion plane.

Once the tidal energy truncation has reached the stellar component (in the example shown, this corresponds roughly to snapshot ②), the size and velocity dispersion of stars and dark matter evolve in sync, both decreasing monotonously.<fnref n="12" />

This is the regime we refer to as *tidally limited* in E+22. In this regime, for the exp2D stellar tracer, we find $R_h/r_{\text{mx}} \approx 0.8$. Similarly, we find $R_h/r_{\text{mx}} \approx 0.9$ for the exp3D stellar tracer, and $R_h/r_{\text{mx}} \approx 0.6$ for the log2D model. For all three models, in the tidally limited regime,<fnref n="13" /> $\langle\sigma_{\text{los}}^2\rangle^{1/2}/V_{\text{mx}} \approx 0.5$.

The initial shape of the stellar density profile appears to play only a secondary role in the evolution of half-light radius and

luminosity-averaged line-of-sight velocity dispersion; the model predictions for the exp2D, exp3D, and log2D stellar tracers in Figure 7(b) differ very little.

We can now return to the objects marked “unidentified” in Figure 7(a). If these objects were dark matter dominated and only weakly affected by tides, they would have velocity dispersions that fall in the gray-shaded band of Figure 7(b). If instead these objects were subjected to strong tides, their velocity dispersions should fall within the region highlighted in blue in Figure 7(b). In either of the two cases, for many of the “unidentified” objects, the predicted velocity dispersions for the dark matter-dominated formation scenario are higher than those expected for self-gravitating objects devoid of dark matter (yellow crosses). As an example, for the case of the faint Milky Way satellite Ursa Major 3/Unions 1, in Figure 7(b), we show a velocity dispersion estimate of $\sigma_{\text{los}} = 1.9_{-1.1}^{+1.4}$ (obtained after removing the furthest outlier in velocity as in Smith et al. (2024)). Note that this dispersion estimate may still be inflated by the presence of binary stars. Taking for now the available dispersion estimate at face value, it appears that UMa3/U1 is consistent with being a dark matter-dominated system. If in dynamical equilibrium and devoid of dark matter, UMa3/U1 would be expected to have a velocity dispersion of only $\sim 50 \text{ m s}^{-1}$ (Errani et al. 2023b), below the range plotted in Figure 7(b).

These results suggest that accurate velocity dispersion estimates of microgalaxy candidates can help to distinguish them from globular cluster remnants: If devoid of dark matter, each “ambiguous” object should have a velocity dispersion that is roughly consistent with that expected for a self-gravitating object of its luminosity. On the other hand, if dark matter dominated, the velocity dispersion should be in excess of that expected for a self-gravitating object of its luminosity. Crucially, our model predictions for the size–velocity dispersion relation of microgalaxies are relatively insensitive to the (unknown) initial energy distribution of stars within the dark matter halo. The required level of measurement accuracy is, however, challenging, with (“virial,” Equation (18)) velocity dispersions in the dark matter-dominated case ranging from several hundred meters per second to a few kilometers per second. Accurate multiepoch spectroscopy will be required to overcome uncertainties due to the potential contribution of binaries to the observed dispersions (see, e.g., McConnachie & Côté 2010; Koposov et al. 2011; Minor et al. 2019).

### 4.3. Dynamical Mass-to-light Ratio

As discussed in the previous section, microgalaxies would differ from globular clusters in their dark matter content. The virial theorem can be used to estimate the dynamical mass enclosed within the luminous radius of a stellar system through the combined measurement of half-light radius $R_h$ and luminosity-averaged velocity dispersion $\langle\sigma_{\text{los}}^2\rangle$ (Illingworth 1976; Merritt 1987; Amorisco & Evans 2012). We will in the following adopt the mass estimator derived in Errani et al. (2018), which is insensitive to anisotropies in the velocity dispersion, and minimizes uncertainties introduced by the (unknown) dark matter density profile shape and extent. The estimator yields the mass enclosed within a spherical radius of $\sim 1.8 R_h$,

$$M_{1.8} \equiv M(<1.8 R_h) \approx 3.5 \times 1.8 R_h \sigma_{\text{los}}^2 G^{-1}. \quad (21)$$

$$\langle\sigma_{\text{los}}^2\rangle = \frac{GM_{\text{cut}}(r_*/r_{\text{cut}})(3 + r_*/r_{\text{cut}})}{6 r_{\text{cut}}(1 + r_*/r_{\text{cut}})^3}. \quad (20)$$

<footnote n="12">¹² For heavily stripped stellar systems in LCDM, a decrease in the line-of-sight velocity dispersion is always accompanied by a decrease in size of the bound component. The velocity dispersions and sizes of the “feeble giant” (Torrealba et al. 2016) satellites Ant 2, Cra 2, And 19, And 25 can therefore not be reached through tides from initial conditions within the gray band of Figure 7(b) (assuming they are bound objects in dynamical equilibrium); see Borukhovetskaya et al. (2022), E+22.</footnote>

<footnote n="13">¹³ A simple analytical estimate for the line-of-sight velocity dispersion in the tidally limited regime can be obtained from Equation (18), assuming as stellar profile a 3D exponential sphere (Equation (12), where $r_* \approx R_h/2.03$), and as underlying potential a truncated NFW cusp (Equation (3), where $r_{\text{cut}} \approx r_{\text{mx}}/1.79$, and $\sqrt{GM_{\text{cut}}/r_{\text{cut}}} \approx V_{\text{mx}}/0.546$). For this choice of stellar profile and dark matter potential, the integral in Equation (18) has an analytical solution, and</footnote>

<footnote>For a tidally limited galaxy, $R_h \approx r_{\text{mx}}$, and Equation (20) yields $\langle\sigma_{\text{los}}^2\rangle^{1/2} \approx 0.5 V_{\text{mx}}$.</footnote>
In Figure 7(c), we compare the dynamical mass $M_{1.8}$ of Local Group dwarf galaxies and globular clusters against their total luminosity $L$. The globular cluster luminosity and dynamical mass are distributed with little scatter around a line of constant dynamical mass-to-light ratio $M_{1.8}/L \approx 1 M_{\odot}/L_{\odot}$. Dwarf galaxies, on the other hand, span a wide range of dynamical mass-to-light ratios, $10 \lesssim (M_{1.8}/L)(M_{\odot}/L_{\odot})^{-1} \lesssim 10^4$. For the case of the faint stellar system UMa3/U1, taking the velocity dispersion estimate of $\sigma_{10\text{os}} \approx 1.9 \text{ km s}^{-1}$ at face value, we find a dynamical mass-to-light ratio of order $\sim 10^3$, suggesting that UMa3/U1 is consistent with being heavily dark matter dominated.

The tidal evolution of the dynamical mass-to-light ratio depends on the shape of the stellar density profile. The dynamical mass-to-light ratio of the exp2D and exp3D models (shown in Figure 7(c) as red and orange curves, respective) gradually increases in the regime of heavy mass loss, as the stellar energy distribution of these two models drops more steeply toward the most-bound states than that of the surrounding dark matter halo (see Figure 5). As the tidal mass loss progresses, the exp2D and exp3D models become more and more dark matter dominated. In contrast, the mass-to-light ratio of the log2D stellar model (with a stellar energy distribution that traces that of the dark matter) converges during tidal stripping to a constant value. This asymptotic mass-to-light ratio is lower than the mass-to-light ratio at accretion. Yet, for the examples considered here, the dwarf galaxies do remain dark matter dominated throughout their evolution. Microgalaxies can therefore be distinguished from globular clusters by their dark matter content, either directly, through measurements of their stellar velocity dispersions (see Section 4.2), or indirectly, by studying their susceptibility to tidal forces along their orbit (see Section 4.4).

### 4.4. Enclosed Mean Density

Using the enclosed dynamical mass $M_{1.8}$ defined by Equation (21), we can estimate the density of a stellar system averaged within a spherical radius of $1.8 R_h$,

$$ \bar{\rho}_{1.8} = M_{1.8} (1.8 R_h)^{-3} (4\pi/3)^{-1}. \qquad (22) $$

Figure 7(d) shows the enclosed mean densities $\bar{\rho}_{1.8}$ of Milky Way satellites as well as their half-light radii $R_h$. All elements of Figure 7(d) are analogous to Figure 7(b), but expressed in terms of density. This allows for an intuitive description of the tidal evolution of dwarf galaxies. In the early stages of tidal evolution, the luminosity and extent of the stellar component are hardly affected. Yet, the mean density enclosed within the luminous radii drops, as dark matter on weakly-bound orbits is stripped (snapshot ①). Once the tidal energy truncation affects the stellar component, its size decreases (snapshot ②). At the same time, its mean density gradually increases: the smaller the remaining stellar system, the higher the average density of the underlying cuspy dark matter halo enclosed within the luminous radii.

Returning to the objects marked as “unidentified” in Figure 7(a), we note that, for many of these objects, in the absence of dark matter their mean densities (yellow crosses) would be substantially lower than the model prediction for dark matter-dominated microgalaxies (blue region). For the case of the faint satellite UMa3/U1, the mean density estimated for a

velocity dispersion of $\sigma_{10\text{os}} \approx 1.9 \text{ km s}^{-1}$ is roughly consistent with the LCDM prediction. If devoid of dark matter, the low stellar mass of UMa3/U1 ($M_* \approx 16 M_{\odot}$) would result in a mean density well below the LCDM prediction, falling roughly on the lowermost dashed curve in Figure 7(d).

Remarkably, if self-gravitating, the mean densities of several of the “unidentified” systems are comparable to the mean density of the Milky Way at the solar circle, $\bar{\rho}_{\text{MW}}(R_{\odot}) \approx 5 \times 10^7 M_{\odot} \text{ kpc}^{-3}$ (assuming $R_{\odot} = 8.3 \text{ kpc}$ and a circular velocity of $240 \text{ km s}^{-1}$). This in turn implies that these objects are likely vulnerable to Galactic tides, which offers an alternative probe to constrain their potential dark matter content.

Studying the vulnerability to tidal disruption of stellar systems is of particular relevance for those objects where accurate estimates of the internal velocity dispersion are unavailable, or deemed unreliable because of the potential contribution of binary stars to the measured dispersion. As long as the Galactocentric velocity of a stellar system can be measured and its orbit constrained, the strength of the tidal field experienced on that orbit can serve to inform whether the presence of dark matter is necessary to avoid full tidal disruption over a time span comparable to the stellar ages in the system. We explore this method with more detail in Errani et al. (2023b) for the example of the Milky Way satellite Ursa Major 3/Unions 1 (Smith et al. 2024), concluding that only if stabilized by the presence of dark matter could the system survive for more than two radial orbital periods on its current orbit.

### 4.5. How Many?

The model outlined in Section 3 enables predictions of the luminosities, sizes, and velocity dispersions of heavily stripped dwarf galaxies—but not their expected abundance in the Milky Way. Simulation studies (e.g., van den Bosch et al. 2005; Springel et al. 2008; Libeskind et al. 2010; Errani et al. 2017) and semianalytical frameworks (e.g., Han et al. 2016; Green et al. 2021) have been used to estimate the subhalo mass function, and/or the radial distribution of subhaloes within a larger host. Tollerud et al. (2008), e.g., suggest that the Milky Way could host between $\sim 300$ and $\sim 600$ luminous satellites within 400 kpc of the galactic center. Similarly, Manwadkar & Kravtsov (2022) compute a total of roughly $\sim 440$ satellites with $M_V < 0$ (i.e., $\log_{10} L/L_{\odot} \gtrsim 1.9$) and half-light radii $R_h \gtrsim 10$ pc within 300 kpc. Ahvazi et al. (2024) estimate a total of $\sim 300$ satellites with $M_V < 0$ within 300 kpc, while Newton et al. (2018) argue for a somewhat lower value of $\sim 120$ satellites. The actual number of heavily stripped satellites in the inner regions of the Milky Way is shown to be sensitive to its detailed mass assembly history (see, e.g., Bose et al. 2020).

The number of luminous satellites is determined by the interplay of the subhalo mass function, and the lowest halo mass that allows star formation. For the Aquarius A halo, Springel et al. (2008) measure $d \ln N_{\text{sub}}/d \ln M_{\text{sub}} \approx -1.9$, i.e., roughly speaking, the number of subhaloes per log-spaced mass bin is inversely proportional to the subhalo mass. Only haloes above some threshold mass allow hydrogen gas to cool in the presence of the cosmic UV background, and in turn, to collapse and form stars (see, e.g., Bullock et al. 2000; Gnedin 2000). The critical mass model discussed in Benitez-Llambay & Frenk (2020), for example, implies a redshift-dependent minimum halo mass for the star formation, which, at $z=0$, equals $\sim 5 \times 10^9 M_{\odot}$ (see Appendix B for a discussion
on the redshift dependence of our results). Galaxies that formed their stars prior to reionization may have formed in haloes of substantially smaller masses: e.g., Bovill & Ricotti (2009) suggest a threshold of ~2 × 10⁸ M☉, while Manwadkar & Kravtsov (2022) suggest that half of the haloes with a peak mass of only ~4 × 10⁷ M☉ may host a galaxy.

For a numerical example, we use as before the Aquarius-A2 main halo as a guide (see Section 2). Within the $z=0$ virial radius ($r_{200} \approx 246$ kpc), we count a total of $N_{\text{sub}} = 16, 65, 152$, and 347 subhaloes with peak masses above $5 \times 10^9 M_{\odot}$, $1 \times 10^9 M_{\odot}$, $5 \times 10^8 M_{\odot}$, and $2 \times 10^8 M_{\odot}$, respectively. To roughly estimate the number of subhaloes that reach the inner regions of a Milky Way-like dark matter halo, we return to study the orbits of subhaloes in the Aquarius-A2 merger tree. Using the same setup as in Section 2, we treat subhaloes as point-masses, and integrate their orbits from their respective redshift of accretion $z_{\text{acc}}$ until $z=0$. We assume a time-evolving, analytical potential fitted to the A2 main halo (see Section 2.1 for details). This allows us to follow the orbits of all subhaloes resolved at accretion until $z=0$, without losing subhaloes to artificial disruption. We find that within a Galactocentric radius of $r_{200} \approx 246$ kpc there are 31 subhaloes with peak virial masses $M_{\text{peak}} \ge 1 \times 10^9 M_{\odot}$ on orbits with pericentric distances $r_{\text{peri}} \le 20$ kpc, and 16 subhaloes have $r_{\text{peri}} \le 10$ kpc. Decreasing the mass threshold to $M_{\text{peak}} \ge 2 \times 10^8 M_{\odot}$ increases the number of subhaloes with $r_{\text{peri}} \le 20$ kpc to 155, and the number of those with $r_{\text{peri}} \le 10$ kpc increases to 82. Note that cosmological variance will introduce a substantial halo-to-halo scatter of these numbers, in particular at the higher-mass end (see Figure 8 in Springel et al. 2008 comparing the subhalo abundance in the Aquarius A–F main haloes).

The strong dependence of the subhalo abundance on the choice of threshold mass is a direct consequence of the steep underlying subhalo mass function. Observational constraints on the abundance of microgalaxies would inform us equally about the number of dark matter subhaloes in the inner Milky Way, and about the minimum mass of star-forming subhaloes.

# 5. Summary and Conclusions

In the LCDM cosmology, dark matter haloes are predicted to have remarkably dense centers where the density profiles formally diverge as $d \ln \rho/d \ln r = -1$ for $r \to 0$. These density cusps render LCDM subhaloes resilient to the effect of tides, and prevent the full tidal disruption of subhaloes and embedded dwarf galaxies. The hierarchical accretion history of the Milky Way may therefore give rise to a population of heavily stripped “microgalaxies,” i.e., comoving groups of stars, held together and shielded from Galactic tides by a surrounding dark matter subhalo.

The resilience to full tidal disruption distinguishes LCDM from other dark matter theories, many of which predict haloes with central constant-density cores, and/or lower mean central densities, such as ultralight-particle DM models (“fuzzy dark matter”), or SIDM models (prior to the onset of core collapse).

In this study, we have used an empirical framework to model the evolution of size, luminosity, and velocity dispersion of stellar systems embedded in LCDM subhaloes, following their tidal evolution over many orders of magnitude in mass, luminosity, and size. We discuss our findings in the context of recently discovered faint stellar systems with structural properties at the boundary of the globular cluster and dwarf galaxy regimes.

## 5.1. Caveats

We make several modeling assumptions in this study, which we summarize here as caveats to our conclusions.

(i) We assume that dark matter haloes prior to accretion are well approximated by NFW profiles (Navarro et al. 1996, 1997), which centrally diverge as $\rho(r) \sim r^{-1}$. High-resolution cosmological simulations suggest that Einasto profiles, with a centrally decreasing power-law slope, provide a slightly better description of the inner regions of CDM haloes (Navarro et al. 2010; Ludlow et al. 2013; Wang et al. 2020), whereas “prompt cusps” (Delos & White 2023a) are argued to cause a steep $\rho(r) \sim r^{-3/2}$ divergence in the innermost regions (Ishiyama et al. 2010; Ogiya & Hahn 2018; Delos & White 2023b).

(ii) Baryonic feedback processes may erode dark matter cusps, rendering CDM subhaloes more vulnerable to tidal mass loss and disruption. The existence of “microgalaxies” hinges on the dark matter cusps being intact, and stars populating the most-bound energy states within them.

(iii) We model stars as tracers of the underlying dark matter potential. Stellar systems that are gravitationally dominant over a dark matter component, or fully self-gravitating, cannot be described by our framework.

(iv) The evolution of both dark matter and stars assumes spherical, dispersion-supported systems with isotropic kinematics. The tidal evolution of rotating or strongly anisotropic systems cannot be modeled with the current framework.

(v) We base our analysis on tidal evolutionary tracks derived from high-resolution N-body simulations, assuming a power-law extrapolation for remnant masses that are unresolved in the simulations. Recent (semi)analytical models propose tidal tracks with (slightly) different asymptotics (Amorisco 2021; Stücker et al. 2023).

(vi) The concentration–mass–redshift relation underlying this study is calibrated to dark matter-only cosmological simulations, ignoring the effect of baryonic feedback on halo concentrations. Note however that, for the halo masses of interest in this study, comparison of the scatter in halo concentration in our initial conditions against hydrodynamical simulations shows good agreement; see Appendix A.

(vii) Finally, we assume a redshift-dependent minimum halo mass for star formation following the Benitez-Llambay & Frenk (2020) model (after reionization) and Tegmark et al. (1997) model (prior to reionization). We discuss the systematics arising from lowering the threshold halo mass for star formation in Appendix B.

## 5.2. Conclusions

Our main conclusions are summarized below.

(i) Current cosmological simulations cannot accurately predict the structure, kinematics, and abundance of subhaloes and embedded dwarf galaxies in the innermost regions of the Galaxy. An insufficient N-body particle number and spatial resolution drive the artificial disruption of substructures in the inner regions of simulated Milky Way–like haloes. A route to circumvent such
resolution limits are analytical models like the one discussed in this work.

(ii) Consistent with earlier work, we show that, if a stellar tracer populates the energy states of a CDM halo all the way down to the most-bound state, then smooth tidal fields cannot fully disrupt that stellar tracer. This suggests the possible existence of microgalaxies, comoving stellar systems of low total luminosity, stabilized against Galactic tides by a surrounding dark matter subhalo.

(iii) The evolution of luminosity and size of such objects is very sensitive to the (unknown) initial binding energy distribution of a stellar tracer within its surrounding dark matter halo. Stellar systems with near-identical surface brightness profiles may show substantial differences in their binding energy distribution.

(iv) In contrast, the tidal evolution of half-light radius and velocity dispersion only weakly depends on the initial stellar energy distribution.

(v) Velocity dispersion measurements provide a robust criterion to distinguish self-gravitating stellar systems from dark matter-dominated microgalaxies. For microgalaxies with a half-light radii of $1 \lesssim r_h/pc \lesssim 100$, we predict (virial) velocity dispersions between several hundred meters per second and a few kilometers per second. Resolving such dispersions for faint stellar systems represents a considerable challenge. Multiepoch spectroscopy will be required to accurately account for the contribution of binary stars to the observed dispersions.

(vi) Dark matter-dominated microgalaxies have larger mean densities than their self-gravitating counterparts of equal luminosity. For systems where the available velocity measurements are sufficiently accurate to constrain their orbits in the Milky way, but not sufficiently precise to resolve their internal dynamics, the systems' susceptibility to tidal mass loss can serve as an alternative criterion to distinguish globular clusters from microgalaxies.

The discovery of dark matter-dominated microgalaxies would inform us about the ability of dark matter subhaloes to survive in the Milky Way tidal field, and at the same time, suggest that stars can populate the most-bound energy states of a halo. Hence, microgalaxies would offer equal insight in the physical properties of dark matter subhaloes, and the baryonic processes that shape the galaxies embedded therein.

## Acknowledgments

R.E. and R.I. acknowledge funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 834148). R.E. and M.W. acknowledge support from the National Science Foundation (NSF) grants AST-2206046 and AST-1909584.

## Appendix A

### Comparison of the Scatter in the Initial Conditions against Cosmological Simulations

In this section, we show that our initial conditions are consistent with the distribution of circular velocity profiles measured for dark matter haloes in cosmological N-body simulations.

As initial conditions to our modeling framework, we assume NFW density profiles (Equation (1)) with halo concentrations

that follow the Ludlow et al. (2016) concentration–mass–redshift relation (see Section 3). We take into account a (logarithmic) scatter of 0.15 dex in concentration $c = r_{200}/r_s$ at a given virial mass $M_{200}$ (for definition, see footnote 5). The value of 0.15 dex used in this study is a conservative value motivated by the scatter measured in cosmological dark matter–only simulations. Dutton & Macciò (2014) find that, at $z=0$, the scatter around the concentration–mass relation is well approximated by a log-normal distribution of width 0.11 dex. Similarly, Moliné et al. (2017, see their Figures 4 and 5) measure a scatter between 0.10–0.15 dex for the ELVIS (Garrison-Kimmel et al. 2014) and VL-II (Diemand et al. 2008) simulations. At higher redshifts, the expected scatter is slightly lower (see, e.g., Ludlow et al. 2016, Figure B3).

At a fixed value of $M_{200}$, assuming NFW profiles, the scatter in the halo characteristic size $r_{mx} \approx 2.16 \times r_{200}/c$ is identical to the scatter in concentration. Note that $r_{mx}$ defines the characteristic size of a halo (and, through $V_{mx}$, it is mean density) only in combination with a functional form for the density profile. Cosmological simulations show that the circular velocity profiles of CDM haloes have very flat maxima (see, e.g., Figure 1 in Navarro et al. 2010). Small deviations from the NFW functional form of an individual halo, as well as N-body discreteness noise, may result in substantial scatter in the exact location of the maximum of the circular velocity profile even in haloes with otherwise very similar density profiles (see, e.g., Knebe et al. 2013).

In Figure 8, we compare the circular velocity profiles of haloes studied in Oman et al. (2015, their Figure 2) against the equivalent profiles modeled using our simple assumptions. The profiles shown in Oman et al. (2015) are obtained from haloes in the EAGLE-HR (Schaye et al. 2015) and APOSTLE-L2 (Fattahi et al. 2016) simulations. Note that, for halo masses

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Average circular velocity profiles $V_c = [GM(<r)/r]^{1/2}$ for haloes with peak circular velocities $32 < V_{mx}/km s^{-1} < 45$ and $45 < V_{mx}/km s^{-1} < 63$. Profiles computed analytically using the same assumptions as adopted in this work are shown as filled blue circles with error bars, spanning from the 10th to the 90th percentile of the underlying distribution. Gray shaded bands are taken from Oman et al. (2015, their Figure 2) and show the same range of circular velocity profiles measured in the EAGLE-HR (Schaye et al. 2015) and APOSTLE-L2 (Fattahi et al. 2016) simulations. The shading is omitted at radii below the Power et al. (2003) convergence radius.</figcaption></figure>
relevant to our study, the profiles discussed in Oman et al. (2015) are virtually identical between the dark matter-only and hydrodynamical runs. We model the concentration-driven scatter in circular velocity profiles by (1) drawing halo virial masses $M_{200}$ from a mass function $dN/dM_{200} \sim M_{200}^{-1.9}$ (see, e.g., Springel et al. 2008; Benson 2020) (2) drawing halo concentrations $c$ from the Ludlow et al. (2016) mass-concentration relation at redshift $z=0$ assuming a log-normal scatter of 0.15 dex, (3) computing NFW circular velocity profiles for the sampled values of $M_{200}$ and $c$, and (4) binning the circular velocity profiles using the same bins in $V_{mx}$ as adopted in Oman et al. (2015). We find excellent agreement between the $N$-body circular velocity profiles of Oman et al. (2015) and our model, shown respectively as gray bands and filled blue circles with error bars in Figure 8.

## Appendix B

### Formation Redshift and Its Effect on the Structure and the Kinematics of Tidally Limited Dwarfs

The characteristic density of subhaloes is related to their collapse redshift: on average, subhaloes that collapsed earlier also have larger characteristic densities, reflecting the larger mean density of the universe at collapse. Similarly, the threshold halo mass for star formation shows a redshift-

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

**Figure 9.** Top panel: Virial mass threshold for star formation as a function of redshift, adapted from Pereira-Wilson et al. (2023). Halo masses computed using the Benitez-Llambay & Frenk (2020) critical mass model are shown as filled circles. For redshifts $z \ge 11$, we also show the halo mass constraints for hydrogen cooling as in Tegmark et al. (1997; filled squares). For both models, the threshold halo mass decreases with increasing redshift. Bottom panel: characteristic density $\rho_{mx} = M_{mx}/(4\pi/3 r_{mx}^3)$ for haloes with virial mass as in the top panel, using the Ludlow et al. (2016) concentration-mass-redshift relation ($\pm 0.15$ dex scatter in concentration, shaded band). The characteristic density of haloes at the threshold mass for star formation monotonously increases with redshift. For reference, a dashed curve shows a density of 200 $\times$ the critical density of the universe (see footnote 5).

dependence, and was lower at higher redshifts than it is today. In this section, we discuss the observational consequences of the interplay between these two redshift-dependent systematics.

In the top panel of Figure 9, we show the typical threshold (virial) mass for star formation computed using the Benitez-Llambay & Frenk (2020) critical mass model, as well as for redshifts prior to reionization using the hydrogen cooling limit as in Tegmark et al. (1997). For both models, the threshold halo mass decreases as redshift increases. Combining this threshold virial mass with halo concentrations computed from the Ludlow et al. (2016) concentration-mass-redshift relation allows to calculate the evolution of the characteristic density

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

**Figure 10.** Top panel: like Figure 7(d), showing the average density $\rho_{1.8}$ of a stellar tracer (exp3D, see Section 3.3.2) embedded in a dark matter halo. The halo masses are chosen to match the redshift-dependent threshold mass for star formation as shown in Figure 9. Solid curves show $\rho_{1.8}$ for stellar tracers embedded in NFW halos at redshifts $z = 0, \dots, 20$. Dashed curves show the tidal evolution of a tidally limited stellar tracer ($R_h = r_{mx}$), which defines a lower bound on $\rho_{1.8}$ for a given half-light radius $R_h$ and redshift $z$. Bottom panel: like Figure 7(b), showing the luminosity-averaged velocity dispersion, but using the same threshold halo masses as in the top panel. For a given half-light radius, a tidally limited stellar system that formed at a higher redshift is expected to have a larger velocity dispersion than a system that formed at a lower redshift.
$\rho_{\text{mx}} = M_{\text{mx}}/(4\pi/3 r_{\text{mx}}^3)$ of the least-massive star-forming haloes. This evolution is shown in the bottom panel of Figure 9: even though the threshold mass for star formation decreases with redshift, the characteristic density of star-forming haloes monotonically increases with redshift. Note that the (rather flat) concentration-mass relation at high redshifts compresses the difference in characteristic density between the Benitez-Llambay & Frenk (2020) and Tegmark et al. (1997) models. The shaded bands corresponds to a scatter in concentration of ±0.15 dex.

This increase of the characteristic density has direct observable consequences for tidally limited systems. Figure 10 is identical in structure to Figures 7(b) and (d) and shows the average density $\rho_{1.8}$ (top panel) and the luminosity-averaged velocity dispersion $\langle\sigma_{\text{los}}^2\rangle^{1/2}$ (bottom panel) of exponential (exp3D) stellar tracers. Here, we choose halo masses corresponding to the minimum mass for star formation at redshifts $z=0,...,20$, with a matching average concentration for the respective redshift.

Solid curves show $\rho_{1.8}$ (top panel) and $\langle\sigma_{\text{los}}^2\rangle^{1/2}$ (bottom panel) for stellar tracers with a half-light radius $R_h$ embedded in the (initial) NFW halo at different redshifts. Dashed curves of matching color show the tidal evolution for a tidally limited stellar tracer ($R_h = r_{\text{mx}}$), delimiting the minimum density (and velocity dispersion) of a tidally stripped system for a given half-light radius $R_h$ and formation redshift $z$. Tidally limited stellar tracers embedded in progenitor haloes that collapsed at earlier (i.e., higher) redshift are expected to have larger average densities and larger average velocity dispersions at fixed $R_h$ than their counterparts that formed at lower redshift.

## Appendix C

### Asymptotic Slopes of DSph Surface Brightness Profiles

In Section 3.3, we discuss the effect of different stellar binding energy distributions on the tidal evolution of dwarf galaxy structural parameters. Assuming a stellar tracer deeply embedded in the power-law cusp of an NFW potential, spherical symmetry, and isotropic kinematics, the binding energy distribution can be estimated from the observed surface brightness profile. In this appendix, we broadly discuss how the stellar binding energies shape the observed surface brightness profile (Section C.1), and compare the surface brightness profiles underlying the present work against a selection of observed profiles of Milky Way satellites (Section C.2).

#### C.1. Illustration of the Connection between Binding Energy Distribution and Surface Brightness Profile

To gain some analytical insight into how the distribution of stellar binding energies shapes the stellar surface brightness profile, consider a monoenergetic stellar distribution function,

$$f_{*}(E) = \delta(E - E_{*}), \qquad (\text{C1})$$

centered on $E_*$.

$$\rho_{*}(r) \propto \int_{\text{all v}} dv \ v^2 \delta\left(\frac{v^2}{2} + \Phi(r) - E_{*}\right) \propto \sqrt{E_{*} - \Phi(r)}, \qquad (\text{C2})$$

for $r$ so that $E_* \ge \Phi(r)$; otherwise, $\rho_*(r) = 0$. This monoenergetic density profile is cored ($\rho_*(r) \to \text{const}$ for $r \to 0$) for all

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Top panel: Binding energy distribution $dN/d\mathcal{E}$ for the log2D stellar tracer (Section 3.3.3) embedded in an NFW halo, for $R_h/r_{\text{mx}} = 1/16$. The distribution is normalized to $\max(dN/d\mathcal{E}) = 1$. Each shaded energy band contributes 20% to the total luminosity. Bottom panel: Surface brightness profile of the log2D stellar tracer, with shaded regions illustrating the contribution to the total luminosity of the energy bands highlighted in the top panel. All energy bands with a lower bound in the binding energy produce a cored surface brightness profile. The logarithmic cusp in the surface brightness profile is caused by the energy band that includes the most-bound state.</figcaption></figure>

potentials $\Phi(r)$ with finite central escape velocity. A centrally divergent surface brightness profile is therefore a clear indication of the binding energy states being populated all the way down to the most-bound state.

We now extend this picture to the contribution of bands in energy to the surface brightness profile. The top panel of Figure 11 shows the distribution of binding energies of a stellar tracer with a log2D surface brightness profile (see Section 3.3.3 for definition) embedded in an NFW halo, for $R_h/r_{\text{mx}} = 1/16$. Bands of a binding energy are color-coded, each band containing 20% of the total luminosity. Note that toward the most-bound energies, the stellar energy distribution is asymptotic to the dark matter (NFW) distribution. The bottom panel of Figure 11 shows the corresponding surface brightness profile, color-coded by the respective contribution of each energy band. All energy bands with a lower bound in the binding energy produce a cored surface brightness profile. The logarithmic cusp of the total surface brightness profile results from the energy band that includes the most-bound state. To guide the eye, a power law with a slope of $d \ln \Sigma_*/d \ln R = -0.4$ is shown as a gray dashed curve: in the inner regions of the dwarf, a surface brightness profile where stars follow the dark matter energetically remains very shallow.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Comparison of the stacked and normalized surface brightness profiles of selected Milky Way satellites against the exp2D (—), exp3D (—), and log2D (—) models (Section 3.3) studied in this work. The surface brightness profiles are computed from star counts available through the DECALS survey (Dey et al. 2019) with photometric selection as in Moskowitz & Walker (2020). The cored exp2D and exp3D appear to better describe the data for Dra 1 and Boo 1 than the logarithmically diverging log2D profile. For CVen 1, UMa 2, Ret 2, all profiles are consistent with the available data.</figcaption></figure>

## C.2. Constraints on the Central Asymptotic Slope from Observations

In this section, we compare the surface brightness profiles used in the present study against available observational data of selected Milky Way dwarf galaxies. The observational situation constitutes a conundrum. For the brightest dwarfs with well-measured surface brightness profiles, stellar feedback has likely affected both the dark matter and stellar distributions, causing central constant-density cores in the dark matter profiles (Pontzen & Governato 2012; Oñorbe et al. 2015). Density cores in turn render the dwarf's vulnerable to tidal mass loss and facilitate their full tidal disruption (Peñarrubia et al. 2010; Errani et al. 2023a), preventing them from becoming tidally stripped “micro galaxies.” Faint dwarfs with few stellar tracers on the other hand may have CDM haloes that are virtually unaffected by their embedded stellar components (Peñarrubia et al. 2012), but the low number of available stars in these objects limits to what degree the central slopes of the surface brightness profiles can be constrained observationally.

As an example, we select five Milky Way satellites spanning a broad range of luminosities: Dra 1, CVen 2 (both with $L \approx 3 \times 10^5 L_\odot$), Boo 1 ($\approx 2 \times 10^4 L_\odot$), UMa 2 ($\approx 4 \times 10^3 L_\odot$) and Ret 2 ($\approx 1 \times 10^3 L_\odot$). We compare the exp2D, exp3D, and log2D stellar distributions as defined in Section 3.3 against star counts available through the DECALS survey (Dey et al. 2019) for a photometric selection as in Moskowitz & Walker (2020).

Figure 12 shows the surface brightness profiles of the five Milky Way satellites listed above, after fitting and subtracting a uniform background component. For each dwarf, we normalize the projected radii $R$ by the respective half-light radius $R_h$ of a fitted 2D exponential (exp2D) model. Similarly, the surface brightness of each dwarf is normalized by $\Sigma_*(R_h)$ of the fitted 2D exponential. Colored lines show the surface brightness profiles corresponding to the exp2D, exp3D, and log2D models. For Draco and Bootes 1, the exp2D, exp3D appear to describe the measured profiles better than the log2D model.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Maximum inner asymptotic slope $\Gamma_{max} = -d \ln \Sigma_*(R)/d \ln R$ (Equation (C3)) of the surface brightness profiles of selected Milky Way satellites, consistent with the enclosed luminosity $L(<R)$ and the local surface brightness $\Sigma_*(R)$ at a given radius. The available data for the dwarfs shown above constrain the inner slope to $\Gamma_{max} \lesssim 0.4$, leaving room for a “cuspy” (sub) component to the surface brightness profile. Solid curves show the constraints on $\Gamma_{max}$ that could be obtained from surface brightness profiles that (exactly) follow the exp2D, exp3D, and log2D models.</figcaption></figure>

For the other objects shown, all models provide a reasonable description of the data at those radii where data are available.

A necessary condition for tidally stripped “microgalaxies” to exist is that at least a (sub)component of the stellar binding energy distribution extends to the most-bound state. For spherical systems with isotropic kinematics, this translates to a cuspy (sub)component to the surface brightness profile. To address to what extent the observed surface brightness profiles can constrain the presence (or absence) of a central cuspy component, we adopt an approach originally introduced in Navarro et al. (2004) to constrain the inner asymptotic slopes of dark matter (3D) density profiles. Following the same reasoning, but using projected quantities, we aim to compute the maximum asymptotic power-law slope of the surface brightness profile $\Gamma_{max} \equiv -d \ln \Sigma_*/d \ln R$ for $R \to 0$. We find that the value $\Gamma_{max}$ that is compatible with (1) the observed surface brightness $\Sigma_*(R)$ at radius $R$, (2) the enclosed luminosity $L(<R)$, and (3) a surface brightness profile that either is a power law or flattens off monotonically toward the center is given by

$$ \Gamma_{max} = 2[1 - \Sigma_*(R) / \bar{\Sigma}_*(<R)], \quad (C3) $$

where $\bar{\Sigma}_*(<R) \equiv L(<R)/(\pi R^2)$ is the average surface brightness enclosed within the radius $R$. In Figure 13, we show the value of $\Gamma_{max}$ as constrained by five example Milky Way dwarf galaxies. This simple analysis suggests that central cusps in the surface brightness profile with an asymptotic power-law slope of $\Gamma_{max} \approx 0.4$ are consistent with the available data, allowing in principle for the existence of a (sub)component of the stellar energy distribution that extends all the way to the most-bound state.

## ORCID iDs

Rodrigo Ibata https://orcid.org/0000-0002-3292-9709

Julio F. Navarro https://orcid.org/0000-0003-3862-5076

Matthew G. Walker https://orcid.org/0000-0003-2496-1925
References

Abramowitz, M., & Stegun, I. A. 1972, Handbook of Mathematical Functions
(New York: Dover)

Ahvazi, N., Benson, A., Sales, L. V., et al. 2024, MNRAS, 529, 3387

Amorisco, N. C. 2021, arXiv:2111.01148

Amorisco, N. C., Agnello, A., & Evans, N. W. 2013, MNRAS, 429, L89

Amorisco, N. C., & Evans, N. W. 2012, MNRAS, 419, 184

Battaglia, G., & Nipoti, C. 2022, NatAs, 6, 659

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

Baumgardt, H., Côté, P., Hilker, M., et al. 2009, MNRAS, 396, 2051

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

Benítez-Llambay, A., Frenk, C. S., Ludlow, A. D., & Navarro, J. F. 2019, MNRAS, 488, 2387

Benson, A. J. 2020, MNRAS, 493, 1268

Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, NJ: Princeton Univ. Press)

Borukhovetskaya, A., Navarro, J. F., Errani, R., & Fattahi, A. 2022, MNRAS, 512, 5247

Bose, S., Deason, A. J., Belokurov, V., & Frenk, C. S. 2020, MNRAS, 495, 743

Bovill, M. S., & Ricotti, M. 2009, ApJ, 693, 1859

Bruce, J., Li, T. S., Pace, A. B., et al. 2023, ApJ, 950, 167

Buist, H. J. T., & Helmi, A. 2016, A&A, 589, C3

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

Burkert, A. 2000, ApJL, 534, L143

Cautun, M., & Frenk, C. S. 2017, MNRAS, 468, L41

Cerny, W., Martínez-Vázquez, C. E., Drlica-Wagner, A., et al. 2023, ApJ, 953, 1

Cerny, W., Simon, J. D., Li, T. S., et al. 2023, ApJ, 942, 111

Chiti, A., Frebel, A., Simon, J. D., et al. 2021, NatAs, 5, 392

Colín, P., Avila-Reese, V., Valenzuela, O., & Firmani, C. 2002, ApJ, 581, 777

Collins, M. L. M., Tollerud, E. J., Rich, R. M., et al. 2020, MNRAS, 491, 3496

Collins, M. L. M., Read, J. I., Ibata, R. A., et al. 2021, MNRAS, 505, 5686

Delos, M. S., & White, S. D. M. 2023a, JCAP, 2023, 008

Delos, M. S., & White, S. D. M. 2023b, MNRAS, 518, 3509

Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168

Diakogiannis, F. I., Lewis, G. F., Ibata, R. A., et al. 2017, MNRAS, 470, 2034

Diemand, J., Kuhlen, M., Madau, P., et al. 2008, Natur, 454, 735

D'Onghia, E., Springel, V., Hernquist, L., & Keres, D. 2010, ApJ, 709, 1138

Dutton, A. A., & Macciò, A. V. 2014, MNRAS, 441, 3359

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

Errani, R., Navarro, J. F., Ibata, R., & Peñarrubia, J. 2022, MNRAS, 511, 6001

Errani, R., Navarro, J. F., Peñarrubia, J., Famaey, B., & Ibata, R. 2023a, MNRAS, 519, 384

Errani, R., Navarro, J. F., Smith, S. E. T., & McConnachie, A. W. 2023b, ApJ, 965, 20

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

Fattahi, A., Navarro, J. F., Frenk, C. S., et al. 2018, MNRAS, 476, 3816

Fattahi, A., Navarro, J. F., Sawala, T., et al. 2016, MNRAS, 457, 844

Fellhauer, M., Kroupa, P., Baumgardt, H., et al. 2000, NA, 5, 305

Frenk, C. S., & White, S. D. M. 2012, AnP, 524, 507

Garrison-Kimmel, S., Boylan-Kolchin, M., Bullock, J. S., & Lee, K. 2014, MNRAS, 438, 2578

Gieles, M., Erkal, D., Antonini, F., Balbinot, E., & Peñarrubia, J. 2021, NatAs, 5, 957

Gnedin, N. Y. 2000, ApJ, 542, 535

Green, S. B., van den Bosch, F. C., & Jiang, F. 2021, MNRAS, 503, 4075

Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010, MNRAS, 404, 1111

Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101

Han, J., Cole, S., Frenk, C. S., & Jing, Y. 2016, MNRAS, 457, 1208

Harris, W. E. 1996, AJ, 112, 1487

Hayashi, K., Chiba, M., & Ishiyama, T. 2020, ApJ, 904, 45

Hilker, M. 2006, A&A, 448, 171

Illingworth, G. 1976, ApJ, 204, 73

Irwin, M., & Hatzidimitriou, D. 1995, MNRAS, 277, 1354

Ishiyama, T., Makino, J., & Ebisuzaki, T. 2010, ApJL, 723, L195

Jardel, J. R., Gebhardt, K., Fabricius, M. H., Drory, N., & Williams, M. J. 2013, ApJ, 763, 91

Jensen, J., Hayes, C. R., Sestito, F., et al. 2024, MNRAS, 527, 4209

Ji, A. P., Koposov, S. E., Li, T. S., et al. 2021, ApJ, 921, 32

Jordi, K., Grebel, E. K., Hilker, M., et al. 2009, AJ, 137, 4586

Kelley, T., Bullock, J. S., Garrison-Kimmel, S., et al. 2019, MNRAS, 487, 4409

Knebe, A., Pearce, F. R., Lux, H., et al. 2013, MNRAS, 435, 1618

Koposov, S. E., Gilmore, G., Walker, M. G., et al. 2011, ApJ, 736, 146

Kravtsov, A. 2010, AdAst, 2010, 281913

Kuzma, P. B., Da Costa, G. S., Keller, S. C., & Maunder, E. 2015, MNRAS, 446, 3297

Li, H., Hammer, F., Babusiaux, C., et al. 2021, ApJ, 916, 8

Libeskind, N. I., Yepes, G., Knebe, A., et al. 2010, MNRAS, 401, 1889

Longeard, N., Martin, N., Starkenburg, E., et al. 2018, MNRAS, 480, 2609

Ludlow, A. D., Bose, S., Angulo, R. E., et al. 2016, MNRAS, 460, 1214

Ludlow, A. D., Navarro, J. F., Boylan-Kolchin, M., et al. 2013, MNRAS, 432, 1103

Manwadkar, V., & Kravtsov, A. V. 2022, MNRAS, 516, 3944

Martin, N. F., Geha, M., Ibata, R. A., et al. 2016, MNRAS, 458, L59

Massari, D., Helmi, A., Mucciarelli, A., et al. 2020, A&A, 633, A36

Mateo, M., Olszewski, E. W., Pryor, C., Welch, D. L., & Fischer, P. 1993, AJ, 105, 510

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

Mau, S., Cerny, W., Pace, A. B., et al. 2020, ApJ, 890, 136

McConnachie, A. W. 2012, AJ, 144, 4

McConnachie, A. W., & Côté, P. 2010, ApJL, 722, L209

Merritt, D. 1987, ApJ, 313, 121

Minor, Q. E., Pace, A. B., Marshall, J. L., & Strigari, L. E. 2019, MNRAS, 487, 2961

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

Moskowitz, A. G., & Walker, M. G. 2020, ApJ, 892, 27

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 462, 563

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

Navarro, J. F., Hayashi, E., Power, C., et al. 2004, MNRAS, 349, 1039

Navarro, J. F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 402, 21

Newton, O., Cautun, M., Jenkins, A., Frenk, C. S., & Helly, J. C. 2018, MNRAS, 479, 2853

Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., et al. 2015, MNRAS, 454, 2092

Ogiya, G., & Hahn, O. 2018, MNRAS, 473, 4339

Oman, K. A., Navarro, J. F., Fattahi, A., et al. 2015, MNRAS, 452, 3650

Pascale, R., Posti, L., Nipoti, C., & Binney, J. 2018, MNRAS, 480, 927

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

Peñarrubia, J., Benson, A. J., Walker, M. G., et al. 2010, MNRAS, 406, 1290

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

Peñarrubia, J., Pontzen, A., Walker, M. G., & Koposov, S. E. 2012, ApJL, 759, L42

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

Pontzen, A., & Governato, F. 2012, MNRAS, 421, 3464

Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14

Read, J. I., Walker, M. G., & Steger, P. 2018, MNRAS, 481, 860

Read, J. I., Walker, M. G., & Steger, P. 2019, MNRAS, 484, 1401

Riley, A. H., Fattahi, A., Pace, A. B., et al. 2019, MNRAS, 486, 2679

Santos-Santos, I. M. E., Navarro, J. F., Robertson, A., et al. 2020, MNRAS, 495, 58

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Sestito, F., Zaremba, D., Venn, K. A., et al. 2023, MNRAS, 525, 2875

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

Smith, S. E. T., Cerny, W., Hayes, C. R., et al. 2024, ApJ, 961, 92

Spergel, D. N., & Steinhardt, P. J. 2000, PhRvL, 84, 3760

Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS, 391, 1685

Stücker, J., Ogiya, G., Angulo, R. E., Aguirre-Santaella, A., & Sánchez-Conde, M. A. 2023, MNRAS, 521, 4432

Taibi, S., Battaglia, G., Rejkuba, M., et al. 2020, A&A, 635, A152

Tegmark, M., Silk, J., Rees, M. J., et al. 1997, ApJ, 474, 1

Tollerud, E. J., Bullock, J. S., Strigari, L. E., & Willman, B. 2008, ApJ, 688, 277

Torrealba, G., Belokurov, V., & Koposov, S. E. 2019, MNRAS, 484, 2181

Torrealba, G., Koposov, S. E., Belokurov, V., & Irwin, M. 2016, MNRAS, 459, 2370

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

van den Bosch, F. C., Ogiya, G., Hahn, O., & Burkert, A. 2018, MNRAS, 474, 3043

van den Bosch, F. C., Tormen, G., & Giocoli, C. 2005, MNRAS, 359, 1029

Walker, M. G., Mateo, M., Olszewski, E. W., et al. 2007, ApJL, 667, L53

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

Wang, J., Bose, S., Frenk, C. S., et al. 2020, Natur, 585, 39

Wang, M. Y., de Boer, T., Pieres, A., et al. 2019, ApJ, 881, 118

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

Woo, J., Courteau, S., & Dekel, A. 2008, MNRAS, 390, 1453

Zavala, J., & Frenk, C. S. 2019, Galax, 7, 81