Microgalaxies in LCDM

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.


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(Navarro et al. , 1997)).This means that CDM haloes are predicted to have remarkably high central densities: for NFW profiles, the density formally diverges as d ln  NFW /d ln  → −1 for  → 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  core /d ln  → 0 for  → 0 (Burkert 2000;Spergel & Steinhardt 2000;Colín et al. 2002).
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. 2018and 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. 2023).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 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., co-moving 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 -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. 2023a,b;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 micro galaxy candidates?
The paper is structured as follows.We show how resolution limits of cosmological -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.
Convergence in -body simulations depends on a combination of the spatial resolution (force softening length , or particle mesh cell size Δ), particle number  and time step Δ (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 -body simulations is not enough to reliably trace the tidal evolution of most subhaloes that pass through the inner regions of the Milky Way.

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 virial1 mass  200 ≳ 10 8 M ⊙ 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 (eq.22 and 23 in Buist & Helmi 2016).This setup allows us to follow the orbits of all subhaloes resolved at accretion until  = 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 show2 that, on average, subhaloes on orbits with pericentres  peri ≲ 10 kpc were accreted more than ≳ 10 Gyr ago.A subhalo on an orbit with  peri = 10 kpc and  apo = 50 kpc has completed more than ≳ 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.

Tidal Disruption and Numerical Resolution
To illustrate the effects of insufficient numerical resolution on tidal evolution, we evolve four different -body realizations of an NFW subhalo in a static isothermal host potential with a circular velocity  host = 220 km s −1 (see Eq. 5).The -body 1 The virial radius  200 encloses a mass  200 so that the mean density  200 /(4/3  3 200 ) is 200× larger than the critical density for closure,  crit .At redshift  = 0,  crit = 32 0 /(8 ). 2 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 disc, 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).) and pericentre  peri ( ) for subhaloes in a Milky Way-like host halo as a function of accretion redshift  acc (assuming Aquarius cosmological parameters,  0 = 73 km s −1 Mpc −1 , Ω m = 0.25, Ω Λ = 0.75).The distributions are computed from subhaloes with peak virial masses of  200 ≥ 10 8 M ⊙ in the Aquarius-A2 simulation.Median radii are shown as solid lines, with shaded regions corresponding to the 16 th − 84 th percentiles.On average, subhaloes accreted at earlier  acc have smaller apo-and pericentres than those accreted more recently.Subhaloes with  peri < 20 kpc have  acc ≳ 2 (i.e., they were accreted ≳ 10 Gyr ago, see top axis).For reference, the evolution of the host halo virial radius  200 is shown as a black solid curve.realizations differ only in the number of -body particles used to simulate the subhalo.
Using the code introduced in EP20, available online3, we generate four -body realizations of this NFW subhalo: For the first realization, labelled A1, we choose an -body particle mass identical to that used in the Aquarius-A1 simulation,  p = 1.71 × 103 M ⊙ .Similarly, for the second realization, labelled A2, we use  p = 1.37×104M ⊙ , 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 simulations are summarized in Table 1.As the NFW profile has a divergent total mass, we taper the profile exponentially beyond 10  s (where by  s we denote the NFW scale radius as in Eq. 1).Consequently, the A1, A2, A3 and A4 realizations of our example subhalo consist of  = 4.4 × 10 6 , 5.6 × 10 5 , 1.5 × 10 5 and 1.9 × 10 4 particles, respectively.is chosen to approximately match the (redshift  = 0) hydrogen cooling limit (Benítez-Llambay et al. 2019), and the initial subhalo characteristic size  mx0 is that of a  = 0 average-concentration halo (Ludlow et al. 2016).Listed also are the host halo and orbit properties underlying the  -body and empirical models.

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 co-moving with the subhalo, with a cell size of Δ ≈  mx0 /256 ≈ 27 pc.Superbox uses two additional grids, one co-moving of medium resolution (10  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  = 22 pc (Springel et al. 2008), similar to the particle-mesh cell size Δ ≈ 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 -body particle mass  p , and thereby the total number of -body particles in each simulation.

Simulation Results
The NFW subhalo models are placed on an orbit with pericentric and apocentric distance  peri = 10 kpc and  apo = 50 kpc, respectively, and evolved for 13 Gyr.Figure 2 shows the evolution of  mx , the bound4 mass within the radius,  mx , where the circular velocity peaks, as a function of time.The four -body realizations are shown in different shades of blue.Over the first two decades in mass loss, all four -body 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 -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  mx drops below  mx < 3000, resolution artifacts result in subhalo densities being systematically underestimated.Evidently, a tool different from classical cosmological -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 (see footnote 3).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.

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

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 density profile with a scale radius  s and a scale density  s , (1) Instead of referring to  s and  s directly, we characterize the subhalo using two equivalent parameters: the peak velocity  mx ≈ 1.65  s √︁   s of the subhalo circular velocity profile  c = √︁   (< )/, and the radius  mx ≈ 2.16  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), with  ≈ 0.3 obtained from fits to -body simulations.For  cut / s → ∞, the above equation reduces to the initial NFW profile, whereas for strong tidal truncation,  cut / s → 0, the profile converges to an exponentially truncated cusp, For the exponentially truncated cusp,  mx ≈ 1.79  cut ,  mx ≈ 1.94  cut √︁   cut .The truncated cusp has a total mass of  cut = 4 cut  3 cut , roughly twice as large as the mass  mx ≈ 0.54  cut enclosed within  mx .As tides strip the subhalo, its characteristic velocity  mx and size  mx decrease, following a tidal track (see Peñarrubia et al. 2008) that couples the evolution of  mx to that of  mx relative to their initial values (here, we use the track of EN21; see their equation 5): with  ≈ 0.4 and  ≈ 0.65 obtained from fits to -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 pericentre  peri and apocentre  apo within an isothermal host potential, where by  host we denote the host's constant circular velocity, and  0 is an arbitrary scale radius.The rate of tidal stripping depends on the density contrast between subhalo and host, measured at the pericentre (EN21 equations 12 and 15).For a given pericentre, 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 pericentre, the tidal evolution gradually decelerates until a final remnant state is reached where the subhalo mean density within  mx is ≈ 16 times higher than the mean host density at the pericentre.In terms of time scales, heavily stripped subhaloes converge towards a characteristic circular time ) and stars (red ), corresponding to snapshot 0 shown in Fig. 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  mx / mx0 = 0, −1, ... , −5 (snapshots 1 , ... , 5 ).The stars initially follow a 2D exponential surface brightness profile (model exp2D, see Sec. 3.3.1)with a half-light radius of  h0 =  mx0 /16.The stellar energy distribution is normalized so that max(d /dE ) = 1; the normalization for the dark matter energy distribution is arbitrary.
peri /4 set by the circular time of the host at the pericentre  peri = 2 peri / 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 -body example of Sec.2.2.3 (see Tab. 1), and place it on the same orbit.The time evolution of the dark matter subhalo is shown in Fig. 2: for the range of remnant masses resolved in the -body simulations, there is excellent agreement between the empirical model and the simulation results.During ∼ 13 Gyr of evolution, the subhalo is stripped to a remnant mass ∼ 10 6 times smaller than its initial mass.

Model for the Evolution of Stellar Tracers
We now turn our attention to modelling 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 Φ(), the initial stellar component is fully defined by its energy distribution where by  ★ (E) we denote the initial stellar distribution function, and by  (E) we denote the density of states in the initial NFW potential.Both are functions of energy E = 1 − /Φ 0 , where  =  2 /2 + Φ(), and Φ 0 denotes the potential minimum.For a given initial stellar density profile  ★ (), the stellar distribution function can be obtained by Eddington inversion.We compute  ★ (E) and  (E) using the implementation described in EP20, available online (see footnote 3).We will discuss our choices of initial stellar density profiles separately in Sec.3.3, and now proceed assuming that an appropriate initial d ★ /dE has been defined.The effect of tides on the stellar energy distribution d ★ /dE can be approximated by a tapered truncation.We model the truncation empirically through (E+22 eq.9), where by d ★ /dE | i,t we denote the truncated energy distribution in the initial conditions, and E mx,t is the energy scale beyond which the energy distribution is truncated.The parameters  ≈ 0.85,  ≈ 12 are obtained through fits to -body simulations.
As an example, Fig. 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 ( h / mx = 1/16).Energies E = 1 − /Φ 0 are measured relative to the potential minimum Φ 0 .In this notation, the most-bound energy state (the ground state) is E = 0, and the boundary between bound and unbound particles lies at E = 1.
The energy-truncated system is initially out of equilibrium.We model the return to virial equilibrium empirically using fits to -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 pre-virialization energies E to final energies E f after virialization follows, on average, the empirical relation (E+22 eq.12) with  ≈ 0.8 and  ≈ −3.Note that the pre-virialization energies E are defined in the initial NFW potential, whereas the energies E f after virialization 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  h ≈ 440 pc, and an initial luminosity of  ≈ 10 6 L ⊙ .We embed the dwarf galaxy in a subhalo with the same characteristic mass and size as in the -body example of Sec.2.2.3 (see Tab. 1).This results in a line-of-sight velocity dispersion of ⟨ 2 los ⟩ 1/2 ≈ 12 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 Sec.2.2.3.
After ∼ 6 Gyr of tidal evolution, snapshot 5 , the dwarf galaxy has been stripped to the size of a micro galaxy, with an enclosed dark matter mass of ∼ 10 4 M ⊙ , a luminosity of only ∼ 50 L ⊙ and a size of several tens of parsecs, well beyond the resolution limits of -body simulations like Aquarius-A1 (see Figure 2).After ∼ 12 Gyr of evolution, snapshot 6 , the stellar system has been stripped to roughly solar luminosity, but still encloses ∼ 10 3 M ⊙ of dark matter.

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 , projected halflight radius  h and line-of-sight velocity dispersion  los .In Sec.3.3.1,3.3.2and 3.3.3,we define the three different stellar density profiles that we adopt in our modelling.In Sec.3.3.4,we briefly discuss the similarities and differences between these three models.Finally, in Sec.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).

2D Exponential -exp2D
The surface brightness profiles of many Milky Way dwarf spheroidal 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. 2019for 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 where Σ 0 = /(2 2 ★ ), for a (2D) half-light radius of  h ≈ 1.68  ★ .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)), where  0 denotes the (zeroth-order) modified Bessel function of second type, and Asymptotically, for  → 0 (see Abramowitz & Stegun 1972 Eq. 9.6.8), i.e., the 3D density distribution that generates a 2D exponential surface brightness profile diverges logarithmically at the centre.

2D Logarithmic Cusp -log2D
Finally, we consider a stellar tracer which, in 3D, diverges at the centre like an NFW dark matter profile,  ★ () ∼  −1 for  ≪  ★ .For this model, we choose a stellar density profile formally identical to Eq. 3, where Initial energy distributions for the exp2D ( ), exp3D ( ) and log2D ( ) stellar models (Sec. 3.3.1,3.3.2,3.3.3)embedded in an NFW dark matter halo with parameters as in Tab 1.All three models share the same projected half-light radius of  h0 / mx0 = 1/16.For comparison, the energy distribution of the underlying NFW halo is shown as a grey dashed curve ( ).Bottom left.Corresponding 3D stellar density profiles  ★ ( ).Top right.2D stellar surface brightness profiles Σ ★ ().Bottom right.Line-of-sight velocity dispersion profiles  los ().The luminosity-averaged line-of-sight velocity dispersion (Eq.17) is virtually identical between the models and is shown as a grey band.

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 behaviour of the three stellar models.
In projection, the exp2D, exp3D models are cored, i.e., d ln Σ ★ /d ln  → 0 for  → 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 Fig. 5, scaled to a half-light radius of  h = 440 pc.
The exp2D model, shown in red, is identical to the example stellar tracer shown previously in Fig. 3 and 4. (Deprojected) 3D stellar density profiles are shown in the bottom-left panel of Fig. 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  ★ /d ln  → −1 for  → 0, like an NFW profile.Consequently, the stellar energy distribution

Model
Eq. 3D Asymptotics 2D Asymptotics 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 . Top panels.Time evolution of the surface brightness profiles of the exp2D ( ), exp3D ( ) and log2D ( ) stellar models, with initial conditions as shown in Fig. 5.The exp2D tracer is identical to the model of Fig. 4, with selected snapshots, labelled 0 , ... , 6 , shown for fixed subhalo remnant masses of log 10  mx / 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  h .sight velocity dispersion, shown as a grey-shaded band, is nearly identical for the three models.We denote by  ★ () the (3D) stellar tracer density, normalized so that 4 ∫ ∞ 0  2  ★ () d = 1, and  = 2 ∫ ∞ 0 Σ ★ d is the total luminosity. (< ) is the total mass enclosed within radius .
Eq. 18 follows from the projected virial theorem (see, e.g.Amorisco & Evans 2012;Errani et al. 2018), and guarantees that ⟨ 2 los ⟩ is independent of anisotropies in the velocity dispersion.Note that a stochastic realization of the luminosityaveraged ⟨ 2 los ⟩ 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.

Tidal Evolution of the 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 Fig. 5. Profiles are shown for subhalo remnant masses of log 10  mx / mx0 = 0, −2, −4 , −6, as labelled 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 per cent of its initial characteristic mass, the tidal energy truncation reaches the stellar component (snapshot 2 , compare with Fig. 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 Fig. 4. The surface brightness of the stellar component is shown in colour, whereas the dark matter surface density is shown in grey.After a slight initial expansion of the stellar half-light radius  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  h evolves thereafter in sync (Kravtsov 2010) with the characteristic size  mx (black circle) of the underlying dark matter subhalo.
Returning to Fig. 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 (Fig. 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 towards 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 towards the most-bound states than that of the dark matter, making the stellar component vulnerable to full tidal disruption.The prospect of observing micro galaxies hence depends critically on how stars are distributed energetically within the dark matter subhalo.
The bottom panel of Fig. 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.

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 micro galaxies, and to compare predicted properties with those of observed stellar systems in the Local Group.
In particular, we aim to address the nature of recentlydiscovered stellar systems with half-light radii 1 ≲  h /pc ≲ 30 and luminosities of 10 ≲ /L ⊙ ≲ 10 3 (Torrealba et al. 2019;Mau et al. 2020;Cerny et al. 2023a,b;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 micro galaxies 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.

Half-light Radii and Luminosities
Figure 7a shows a compilation5 of projected half-light radii,  h , and luminosities, , 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  = 10 6 L ⊙ , globular clusters span a range in sizes of roughly 1 ∼ 10 pc, more than an order of magnitude more compact than dwarfs of the same luminosity, 0.1 ∼ 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 to 30 pc, not too different from extended Milky Way globular clusters like Palomar 5 or Palomar 14 (albeit at lower luminosity).Recently discovered objects, whose nature is still unclear, are labelled "unidentified" and depicted as grey crosses in Fig. 7a.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  ★ = 16 + 6 −5 M ⊙ , corresponding to a luminosity of  V ≈ +2.2 + 0.4 −0.3 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 Fig. 7a, computed using the method outlined in Sec. 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 Tab 1, with a viral mass close to the redshift  = 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  = 10 6 L ⊙ , a half-light radius of  h = 440 pc and a velocity dispersion  los ≈ 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 Fig. 3 and 4, with circled numbers 0 , ... , 6 corresponding to fixed subhalo remnant masses of log 10  mx / 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 Fig. 7a 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 steep6: tidal evolution shifts the dwarf galaxy model roughly along the observed masssize 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 grey crosses in Fig. 7a cannot be reached by tidally stripping progenitor systems with initial luminosites 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 micro galaxies, 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 Fig. 7a 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 micro galaxies.

Half-light Radii and Velocity Dispersions
We now discuss to what extent measurements of the line-ofsight velocity dispersion ⟨ 2 los ⟩ 1/2 (Eq.17) can help to constrain the nature of the unidentified stellar systems.Fig. 7b shows the half-light radii  h and line-of-sight velocity dispersions ⟨ 2 los ⟩ 1/2 for the same dwarf galaxies and globular clusters as in Fig. 7a.Dwarf galaxies, on average, are more extended than globular clusters at equal velocity dispersion.

Self-gravitating Star Clusters
For most of the objects marked as "unidentified" in Fig. 7a, 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  h and their total stellar mass  ★ .For a 3D exponential stellar density profile (for definition see Sec. 3.3.2) with  h ≈ 2.03  ★ and total stellar mass  ★ , Eq. 18 gives 6 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 ( h ∝  mx and  h ∝  mx ).Calling  ★ the slope of the stellar energy distribution d ★ /dE ∝ E ★ for E ≪ 1, then  ∝  ( ★+1) /2 h (assuming NFW haloes, where d/dE ∝ E for E ≪ 1).Combining this with the tidal track of Eq. 4 with  h ∝  2+1 h gives d ln /d ln  h = (2 + 1) (  ★ + 1)/2.Hence, the evolution of luminosity and size of a stellar tracer component is directly related to how stars populate the mostbound energy states of the subhalo, parameterized through the log-slope d ln(d /dE )/d ln E =  ★ for E ≪ 1.
For each "unidentified" system shown in Fig. 7a, we use Eq.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 Υ ≈ 1.6 (Woo et al. 2008).The results of this exercise7 are marked as yellow crosses in Fig. 7b.For reference, dashed diagonal lines in Fig. 7b are also computed from Eq. 19, for constant stellar masses of  ★ = Υ = 10 2 , 10 3 and 10 4 M ⊙ .

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. 2010Guo et al. , 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  = 0, this limiting virial mass is roughly ∼ 5 × 10 9 M ⊙ (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  = 0 hydrogen cooling limit ( mx0 = 7 kpc,  mx0 = 35 km s −1 , see Table 1) is shown as a black solid curve in Fig. 7b, labelled "initial".This curve is computed from Eq. 18, with the integral evaluated numerically.We 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 peak circular velocity 20 ≲  mx /km s −1 ≲ 40 (e.g., Fattahi et al. 2018).This range of potential initial haloes is indicated by the grey band in Fig. 7b, which includes ±0.15 dex variation in 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 grey band of Fig. 7b.
Tidal stripping reduces the stellar velocity dispersion, providing a potential explanation as to why some observed dwarfs are found below the grey band.To illustrate this, we show the tidal evolution of half-light radius  h and velocity dispersion ⟨ 2 los ⟩ 1/2 of the three example dwarf galaxy models with exp2D, exp3D and log2D stellar tracers using red, orange and   Luminosities  and half-light radii  h of Local Group dwarf galaxies ( ) and globular clusters ( ).Ambiguous objects at the boundary between the dwarf-and globular cluster regime are marked using grey crosses ( ).For references, see footnote 5. Solid curves show tidal evolutionary tracks for the exp2D ( ), exp3D ( ) and log2D ( ) stellar models (with initial conditions as in Fig. 5).Snapshots corresponding to the exp2D model, labelled 0 , ... , 5 , are highlighted for fixed subhalo remnant masses of log 10  mx / mx0 = 0, −1, ... , −5.(b) Line-of-sight velocity dispersion ⟨  2 los ⟩ 1/2 (Eq.17) and projected half-light radius  h .The top black curve labelled "initial" shows the line-of-sight velocity dispersion expected for a stellar tracer of half-light radius  h embedded in a ∼ 10 9 M ⊙ dark matter halo.Yellow crosses ( ) show velocity dispersions computed under the assumption of self-gravity for the ambiguous objects shown as grey crosses in the top-left panel.If instead these objects are dark matter-supported galaxies with 1 ≲  h /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 to the  = 0 hydrogen cooling limit.For earlier formation redshifts, the lower bound shifts upwards, see Appendix B. (c) Luminosities  and dynamical masses  1.8 (Eq.21).Diagonal lines mark constant dynamical mass-to-light ratios of  1.8 / = 1 M ⊙ /L ⊙ and 100 M ⊙ /L ⊙ .(d) Projected half-light radii  h and mean densities ρ1.8 (Eq.22) enclosed within a spherical radius of 1.8  h .purple curves in Fig. 7b.Both models are initially located at the position marked 0 , and are embedded in the same dark matter halo.As previous, circled numbers 0 , ... , 6 correspond to the snapshots shown in Fig. 4 at fixed subhalo remnant masses of log 10  mx / mx0 = 0, −1, ... , −6.
In the early stages of tidal evolution, the tidal energy truncation does not reach the binding energies of the stars (compare with Fig. 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 2 ), size and velocity dispersion of stars and dark matter evolve in sync, both decreasing monotonously8.This is the regime we refer to as tidally limited in E+22.In this regime, for the exp2D stellar tracer, we find  h / mx ≈ 0.8.Similarly, we find  h / mx ≈ 0.9 for the exp3D stellar tracer, and  h / mx ≈ 0.6 for the log2D model.For all three models, in the tidally limited regime9, ⟨ 2 los ⟩ 1/2 / mx ≈ 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 Fig. 7b differ very little.
We can now return to the objects marked "unidentified" in Fig. 7a.If these objects were dark matter-dominated and only weakly affected by tides, they would have velocity dispersions that fall in the grey-shaded band of Fig. 7b.If instead these objects were subjected to strong tides, their velocity dispersions should fall within the region highlighted in blue in Fig. 7b.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 Fig. 7b, we show a velocity dispersion estimate of  los = 1.9 + 1.4 −1.1 (obtained after removing the furthest outlier from the sample of measured 8 For heavily stripped stellar systems in LCDM, a decrease in the line-ofsight 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 grey band of Fig. 7b (assuming that they are bound objects in dynamical equilibrium), see Borukhovetskaya et al. ( 2022), E+22.9 A simple analytical estimate for the line-of-sight velocity dispersion in the tidally limited regime can be obtained from Eq. 18, assuming as stellar profile a 3D exponential sphere (Eq.12, where  ★ ≈  h /2.03), and as underlying potential a truncated NFW cusp (Eq.3, where  cut ≈  mx /1.79, and √︁   cut / cut ≈  mx /0.546).For this choice of stellar profile and dark matter potential, the integral in Eq. 18 has an analytical solution, and For a tidally limited galaxy,  h ≈  mx , and Eq.20 yields ⟨  2 los ⟩ 1/2 ≈ 0.5  mx .velocities as described 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 ∼50 m/s (Errani et al. 2024), below the range plotted in Fig. 7b.
These results suggest that accurate velocity dispersion estimates of micro galaxy 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 matterdominated, 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 micro galaxies 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", Eq. 18) velocity dispersions in the dark matter-dominated case ranging from several hundred ms −1 to a few km s −1 .Accurate multi-epoch 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).

Dynamical Mass-to-light Ratio
As discussed in the previous section, micro galaxies 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  h and luminosity-averaged velocity dispersion ⟨ 2 los ⟩ (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 ∼ 1.8  h , In Fig. 7c, we compare the dynamical mass  1.8 of Local Group dwarf galaxies and globular clusters against their total luminosity .The globular cluster luminosity and dynamical mass are distributed with little scatter around a line of constant dynamical mass-to-light ratio  1.8 / ≈ 1 M ⊙ /L ⊙ .Dwarf galaxies, on the other hand, span a wide range of dynamical mass-to-light ratios, 10 ≲ ( 1.8 /) (M ⊙ /L ⊙ ) −1 ≲ 10 4 .For the case of the faint stellar system UMa3/U1, taking the velocity dispersion estimate of  los ≈ 1.9 km s −1 at face value, we find a dynamical mass-to-light ratio of order ∼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 Fig. 7c 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 towards the most-bound states than that of the surrounding dark matter halo (see Fig. 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 which 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.Micro galaxies can therefore be distinguished from globular clusters by their dark matter content, either directly, through measurements of their stellar velocity dispersions (see Sec. 4.2), or indirectly, by studying their susceptibility to tidal forces along their orbit (see Sec. 4.4).

Enclosed Mean Density
Using the enclosed dynamical mass  1.8 defined by Eq.21, we can estimate the density of a stellar system averaged within a spherical radius of 1.8  h , Fig. 7d shows the enclosed mean densities ρ1.8 of Milky Way satellites as well as their half-light radii  h .All elements of Fig. 7d are analogous to Fig. 7b, 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, 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 1 ).Once the tidal energy truncation affects the stellar component, its size decreases (snapshot 2 ).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 Fig. 7a, 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 matterdominated micro galaxies (blue region).For the case of the faint satellite UMa3/U1, the mean density estimated for a velocity dispersion of  los ≈ 1.9 km s −1 is roughly consistent with the LCDM prediction.If devoid of dark matter, the low stellar mass of UMa3/U1 ( ★ ≈ 16 M ⊙ ) would result in a mean density well below the LCDM prediction, falling roughly on the lowermost dashed curve in Fig. 7d.
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, ρMW ( ⊙ ) ≈ 5 × 10 7 M ⊙ kpc −3 (assuming  ⊙ = 8.3 kpc and a circular velocity of 240 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. (2024) 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.

How Many?
The model outlined in Sec. 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 semi-analytical 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 ∼300 and ∼600 luminous satellites within 400 kpc of the galactic centre.Similarly, Manwadkar & Kravtsov (2022) compute a total of roughly ∼ 440 satellites with   < 0 (i.e.log 10 /L ⊙ ≳ 1.9) and half-light radii  h ≳ 10 pc within 300 kpc.Ahvazi et al. ( 2024) estimate a total of ∼ 300 satellites with   < 0 within 300 kpc, while Newton et al. (2018) argue for a somewhat lower value of ∼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  sub /d ln  sub ≈ −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 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 star formation which, at  = 0, equals ∼ 5 × 10 9 M ⊙ (see Appendix B for a discussion on the redshift dependence of our results).Galaxies which 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 8 M ⊙ , while Manwadkar & Kravtsov (2022) suggest that half of the haloes with a peak mass of only ∼4 × 10 7 M ⊙ may host a galaxy.
For a numerical example, we use as before the Aquarius-A2 main halo as a guide (see Sec. 2).Within the  = 0 virial radius ( 200 ≈ 246 kpc), we count a total of  sub = 16, 65, 152 and 347 subhaloes with peak masses above 5×10 9 M ⊙ , 1×10 9 M ⊙ , 5 × 10 8 M ⊙ and 2 × 10 8 M ⊙ , 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 Sec. 2, we treat subhaloes as point-masses, and integrate their orbits from their respective redshift of accretion  acc till  = 0. We assume a time-evolving, analytical potential fitted to the A2 main halo (see Sec. 2.1 for details).This allows us to follow the orbits of all subhaloes resolved at accretion till  = 0, without losing subhaloes to artificial disruption.We find that within a Galactocentric radius of  200 ≈ 246 kpc there are 31 subhaloes with peak virial masses  peak ≥ 1 × 10 9 M ⊙ on orbits with pericentric distances  peri ≤ 20 kpc, and 16 subhaloes have  peri ≤ 10 kpc.Decreasing the mass threshold to  peak ≥ 2 × 10 8 M ⊙ increases the number of subhaloes with  peri ≤ 20 kpc to 155, and the number of those with  peri ≤ 10 kpc 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 micro galaxies 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.

Summary and Conclusions
In the LCDM cosmology, dark matter haloes are predicted to have remarkably dense centres where the density profiles formally diverge as d ln /d ln  = −1 for  → 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 "micro galaxies", i.e., co-moving 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", FDM), or self-interacting dark matter (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.

Caveats
We make several modelling assumptions in this study, which we summarize here as caveats to our conclusions.i) We assume that dark matter haloes prior to accretion are wellapproximated by NFW profiles (Navarro et al. 1996(Navarro et al. , 1997)), which centrally diverge as () ∼  −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 () ∼  −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 "micro galaxies" 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 modelled with the current framework.
v) We base our analysis on tidal evolutionary tracks derived from high-resolution -body simulations, assuming a powerlaw 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.

Conclusions
Our main conclusions are summarized below.
i) Current cosmological simulations cannot accurately predict the structure, kinematicsm and abundance of subhaloes and embedded dwarf galaxies in the innermost regions of the Galaxy.Insufficient -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 micro galaxies, co-moving 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 micro galaxies.For micro galaxies with half-light radii of 1 ≲  h /pc ≲ 100, we predict (virial) velocity dispersions between several hundred m s −1 and a few km s −1 .Resolving such dispersions for faint stellar systems represents a considerable challenge.Multi-epoch spectroscopy will be required to accurately account for the contribution of binary stars to the observed dispersions.
vi) Dark matter-dominated micro galaxies 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 micro galaxies.
The discovery of dark matter-dominated micro galaxies 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, micro galaxies would offer equally insight in the physical properties of dark matter subhaloes, and the baryonic processes that shape the galaxies embedded therein.for haloes with peak circular velocities 32 <  mx /km s −1 < 45 and 45 <  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 10 th to the 90 th percentile of the underlying distribution.Grey shaded bands are taken from Oman et al. (2015, their Fig. 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.(Fattahi et al. 2016) simulations.Note that for halo masses 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  200 from a mass function d/d 200 ∼  −1.9 200 (see e.g.Springel et al. 2008;Benson 2020) (2) drawing halo concentrations  from the Ludlow et al. (2016) mass-concentration relation at redshift  = 0 assuming a log-normal scatter of 0.15 dex, (3) computing NFW circular velocity profiles for the sampled values of  200 and  and (4) binning the circular velocity profiles using the same bins in  mx as adopted in Oman et al. (2015).We find excellent agreement between the -body circular velocity profiles of Oman et al. (2015) and our model, shown respectively as grey bands and filled blue circles with error bars in Fig. 8.

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-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 Fig. 9, we show the typical threshold (virial) mass for star formation computed using the Benitez-Llambay & Frenk (2020) critical mass model, as well 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 theshold 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  mx =  mx /(4/3  3 mx ) of the least-massive star forming haloes.This evolution is shown in the bottom panel of Fig. 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 band corresponds to a scatter in concentration of ±0.15 dex.
This increase of the characteristic density has direct observable consequences for tidally limited systems.Fig. 10 is identical in structure to Fig. 7 (b) and (d) and shows the average density  1.8 (top panel) and the luminosity-averaged velocity dispersion ⟨ 2 los ⟩ 1/2 (bottom panel) of exponential (exp3D) stellar tracers.Here, we choose halo masses corresponding to the minimum mass for star formation at redshifts  = 0, . . ., 20, with matching average concentration for the respective redshift.
Solid curves show  1.8 (top panel) and ⟨ 2 los ⟩ 1/2 (bottom panel) for stellar tracers with a half-light radius  h embedded in the (initial) NFW halo at different redshifts.Dashed curves of matching colour show the tidal evolution for a tidally limited stellar tracer ( h =  mx ), delimiting the minimum density (and velocity dispersion) of a tidally stripped system for a given half-light radius  h and formation redshift .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  h than their counterparts that formed at lower redshift.

C. Asymptotic Slopes of dSph Surface Brightness Profiles
In Sec.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 (Sec.C.1), and compare the surface brightness profiles underlying the present work against a selection of observed profiles of Milky Way satellites (Sec.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 mono-energetic stellar distribution function, centred on  ★ .The (spherical) stellar density profile is readily obtained by integrating over velocity space, assuming isotropic kinematics, for  so that  ★ ≥ Φ(), otherwise  ★ () = 0.This monoenergetic density profile is cored ( ★ () → const for  → 0) for all potentials Φ() 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 Fig. 11 shows the distribution of binding energies of a stellar   7b, 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.
tracer with a log2D surface brightness profile (see Sec. profile, colour-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 slope of d ln Σ ★ /d ln  = −0.4 is shown as a grey dashed curve: In the inner regions of the dwarf, a surface brightness profile where stars follow the dark matter energetically remains very shallow.

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 (Eq.C3) of the surface brightness profiles of selected Milky Way satellites, consistent with the enclosed luminosity  (< ) and the local surface brightness Σ ★ () at a given radius.The available data for the dwarfs shown above constrain the inner slope to Γ max ≲ 0.4, leaving room for a "cuspy" (sub)component to the surface brightness profile.Solid curves show the constraints on Γ max that could be obtained from surface brightness profiles that (exactly) follow the exp2D, exp3D and log2D models.
cores in turn render the dwarfs vulnerable to tidal mass loss and facilitate their full tidal disruption (Peñarrubia et al. 2010;Errani et al. 2023), 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.
Fig. 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  by the respective half-light radius  h of a fitted 2D exponential (exp2D) model.Similarly, the surface brightness of each dwarf is normalized by Σ ★ ( h ) of the fitted 2D exponential.Coloured 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.For the other objects shown, all models provide a reasonable description of the data at those radii where data is available.
A necessary condition for tidally stripped "micro galaxies" 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 Γ max ≡ −d ln Σ ★ /d ln  for  → 0. We find that the value Γ max that is compatible with (1) the observed surface brightness Σ ★ () at radius , (2) the enclosed luminosity (< ), and (3) a surface brightness profile that is either a power law, or flattens off monotonically towards the centre, is given by where Σ★ (< ) ≡  (< )/( 2 ) is the average surface brightness enclosed within the radius .In Fig. 13, we show the value of Γ 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 Γ max ≈ 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.

Figure 1 .
Figure 1.Distribution of apocentre  apo () and pericentre  peri ( ) for subhaloes in a Milky Way-like host halo as a function of accretion redshift  acc (assuming Aquarius cosmological parameters,  0 = 73 km s −1 Mpc −1 , Ω m = 0.25, Ω Λ = 0.75).The distributions are computed from subhaloes with peak virial masses of  200 ≥ 10 8 M ⊙ in the Aquarius-A2 simulation.Median radii are shown as solid lines, with shaded regions corresponding to the 16 th − 84 th percentiles.On average, subhaloes accreted at earlier  acc have smaller apo-and pericentres than those accreted more recently.Subhaloes with  peri < 20 kpc have  acc ≳ 2 (i.e., they were accreted ≳ 10 Gyr ago, see top axis).For reference, the evolution of the host halo virial radius  200 is shown as a black solid curve.

Figure 2 .
Figure2.Tidal stripping of a subhalo with an initial mass of  mx0 = 2 × 10 9 M ⊙ in a Milky Way-like host (see Table1), computed using  -body simulations (blue curves) and the empirical EN21 model (black solid curve).The four  -body simulations have particle masses  p chosen to match those of the Aquarius-A1 ( ), A2 ( ), A3 ( ) and A4 ( ) resolution levels, respectively.The evolution of the subhalo mass  mx is shown as a function of time .Numbered circles 0 , ... , 6 mark remnant masses of log 10  mx / mx0 = 0, −1, ... , −6.Shaded regions correspond to initial conditions spanning a range of circular velocities 20 ≤  mx0 /km s −1 ≤ 40.Insufficient numerical resolution causes the  -body models disrupt artificially, whereas the empirical model predicts an evolution towards a stable asymptotic remnant state.

Figure 3 .
Figure 3. Solid curves show the initial energy distribution d /dE of dark matter (grey) and stars (red ), corresponding to snapshot 0 shown in Fig.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  mx / mx0 = 0, −1, ... , −5 (snapshots 1 , ... , 5 ).The stars initially follow a 2D exponential surface brightness profile (model exp2D, see Sec. 3.3.1)with a half-light radius of  h0 =  mx0 /16.The stellar energy distribution is normalized so that max(d /dE ) = 1; the normalization for the dark matter energy distribution is arbitrary.

Figure 4 .
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  2D is shown in grey, and the surface-brightness Σ ★ of the embedded dwarf is shown in colour.Snapshots, labelled 0 , ... , 6 , are shown for fixed subhalo remnant masses of log 10  mx / mx0 = 0, −1, ... , −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  h0 =  mx0 /16 (model exp2D, see Sec. 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  mx (shown as black circles) and the stellar half-light radius  h (red circles) evolve in parallel (last three panels).
Fig 5 for a stellar tracer deeply embedded in an NFW halo ( h / mx = 1/16).Note that the energy distributions of the exp2D and exp3D models drop faster towards the most-bound state than the log2D model.Finally, in the bottom-right panel of Fig 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 Table1.The luminosity-averaged line-of- Half-light radius  h and luminosity Half-light radius  h and velocity dispersion ⟨  2 los ⟩ 1/2 (Eq.17) Luminosity  and dynamical mass  1.8 =  (< 1.8  h ) (Eq. 21)

Figure
Figure 7. (a)Luminosities  and half-light radii  h of Local Group dwarf galaxies ( ) and globular clusters ( ).Ambiguous objects at the boundary between the dwarf-and globular cluster regime are marked using grey crosses ( ).For references, see footnote 5. Solid curves show tidal evolutionary tracks for the exp2D ( ), exp3D ( ) and log2D ( ) stellar models (with initial conditions as in Fig.5).Snapshots corresponding to the exp2D model, labelled 0 , ... , 5 , are highlighted for fixed subhalo remnant masses of log 10  mx / mx0 = 0, −1, ... , −5.(b) Line-of-sight velocity dispersion ⟨  2 los ⟩ 1/2 (Eq.17) and projected half-light radius  h .The top black curve labelled "initial" shows the line-of-sight velocity dispersion expected for a stellar tracer of half-light radius  h embedded in a ∼ 10 9 M ⊙ dark matter halo.Yellow crosses ( ) show velocity dispersions computed under the assumption of self-gravity for the ambiguous objects shown as grey crosses in the top-left panel.If instead these objects are dark matter-supported galaxies with 1 ≲  h /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 to the  = 0 hydrogen cooling limit.For earlier formation redshifts, the lower bound shifts upwards, see Appendix B. (c) Luminosities  and dynamical masses  1.8 (Eq.21).Diagonal lines mark constant dynamical mass-to-light ratios of  1.8 / = 1 M ⊙ /L ⊙ and 100 M ⊙ /L ⊙ .(d) Projected half-light radii  h and mean densities ρ1.8 (Eq.22) enclosed within a spherical radius of 1.8  h .

Figure 8 .
Figure 8.Average circular velocity profiles  c = [  (<  )/ ] 1/2for haloes with peak circular velocities 32 <  mx /km s −1 < 45 and 45 <  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 10 th to the 90 th percentile of the underlying distribution.Grey shaded bands are taken fromOman et al. (2015, their  Fig. 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 thePower et al. (2003) convergence radius.

Figure 9 .
Figure 9. Top panel: Virial mass threshold for star formation as a function of redshift, adapted from Pereira-Wilson et al.Halo masses computed using the Benitez-Llambay & Frenk (2020) critical mass are shown as filled circles.For redshifts  ≥ 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  mx =  mx /(4 /3  3 mx ) for haloes with virial mass as in the top panel, using the Ludlow et al. (2016) concentration-mass-redshift relation (±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× the critical density of the universe (see footnote 1).

Figure 10 .
Figure10.Top panel: like Fig.7d, showing the average density  1.8 of a stellar tracer (exp3D, see Sec. 3.3.2) embedded in a dark matter halo.The halo masses are chosen to match the threshold mass for star formation as shown in Fig.9.Solid curves show  1.8 for stellar tracers embedded in NFW haloes at redshfits  = 0, . . ., 20.Dashed curves show the tidal evolution of a tidally limited stellar tracer ( h =  mx ), which defines a lower bound on  1.8 for a given half-light radius  h and redshift .Bottom panel: like Fig.7b, 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.

Figure 11 .
Figure 11.Top panel: Binding energy distribution d /dE for the log2D stellar tracer (Sec.3.3.3)embedded in an NFW halo, for  h / mx = 1/16.The distribution is normalized to max(d /dE ) = 1.Each shaded energy band contributes 20 per cent 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.

Figure 12 .Figure 13 .
Figure12.Comparison of the stacked and normalized surface brightness profiles of selected Milky Way satellites against the exp2D ( ), exp3D ( ) and log2D ( ) models (Sec.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 inMoskowitz & 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.

Table 1 .
Parameters of the example subhalo model considered in the body experiments of Sec. 2 and the empirical models of Sec. 3. The initial subhalo mass  mx0 Table 1), computed using  -body simulations (blue curves) and the empirical EN21 model (black solid curve).
The four  -body simulations have particle masses  p chosen to match those of the Aquarius-A1 ( ), A2 ( ), A3 ( ) and A4 ( ) resolution levels, respectively.The evolution of the subhalo mass  mx is shown as a function of time .Numbered circles 0 , ... , 6 mark remnant masses of log 10  mx / mx0 = 0, −1, ... , −6.Shaded regions correspond to initial conditions spanning a range of circular velocities 20 ≤  mx0 /km s −1 ≤ 40.Insufficient numerical resolution causes the  -body models disrupt artificially, whereas the empirical model predicts an evolution towards a stable asymptotic remnant state.