Structure and kinematics of tidally limited satellite galaxies in LCDM Free

ABSTRACT

We use N-body simulations to model the tidal evolution of dark matter-dominated dwarf spheroidal galaxies embedded in cuspy Navarro-Frenk-White subhaloes. Tides gradually peel off stars and dark matter from a subhalo, trimming it down according to their initial binding energy. This process strips preferentially particles with long orbital times, and comes to an end when the remaining bound particles have crossing times shorter than a fraction of the orbital time at pericentre. The properties of the final stellar remnant thus depend on the energy distribution of stars in the progenitor subhalo, which in turn depends on the initial density profile and radial segregation of the initial stellar component. The stellar component may be completely dispersed if its energy distribution does not extend all the way to the subhalo potential minimum, although a bound dark remnant may remain. These results imply that ‘tidally limited’ galaxies, defined as systems whose stellar components have undergone substantial tidal mass-loss, neither converge to a unique structure nor follow a single tidal track. On the other hand, tidally limited dwarfs do have characteristic sizes and velocity dispersions that trace directly the characteristic radius (r_mx) and circular velocity (V_mx) of the subhalo remnant. This result places strong upper limits on the size of satellites whose unusually low velocity dispersions are often ascribed to tidal effects. In particular, the large size of kinematically cold ‘feeble giant’ satellites like Crater 2 or Antlia 2 cannot be explained as due to tidal effects alone in the Lambda Cold Dark Matter scenario.

1 INTRODUCTION

The Lambda Cold Dark Matter (LCDM) scenario makes well-defined and falsifiable predictions for the radial density profile of dark matter haloes, as well as for their abundance as a function of virial¹ mass. These predictions are particularly relevant for the study of dwarf spheroidal (dSph) galaxies, dark matter-dominated systems whose stellar components act as simple kinematic tracers of the structure of their surrounding dark haloes (see e.g. Mateo et al. 1993; Walker et al. 2007).

These studies may be used to test the expected density profiles of cold dark matter haloes, which are well approximated by the Navarro–Frenk–White formula (hereafter, NFW; Navarro, Frenk & White 1996b, 1997). This issue has been debated for decades, albeit with mixed results, with some authors arguing that the structure of dwarf galaxy haloes is consistent with the ‘cuspy’ NFW shape and others claiming that the data suggest density profiles with a sizable constant-density ‘core’ (for reviews see e.g. Gilmore et al. 2007; Bullock & Boylan-Kolchin 2017).

The interpretation of these studies is further complicated by the fact that the assembly of the baryonic component of a galaxy may induce changes in the dark matter density profile, perhaps even erasing the expected cusp and imprinting a core (Navarro, Eke & Frenk 1996a; Read & Gilmore 2005; Governato et al. 2010; Pontzen & Governato 2012). However, such baryon-induced effects should be minimal in very faint dark matter-dominated galaxies, simply because the total fraction of mass in baryonic form is too small to be able to affect gravitationally the dark matter component (see e.g. Peñarrubia et al. 2012; Di Cintio et al. 2014; Benítez-Llambay et al. 2019).

Dwarf galaxies may also be used to probe another robust LCDM prediction, concerning the mass function of low-mass haloes. In LCDM these haloes are so numerous, and their mass function so steep, that accommodating the comparatively scarcer number of dwarfs and their much shallower stellar mass function requires that galaxy formation efficiency should decline steadily with decreasing halo mass, effectively restricting dwarf galaxy formation to haloes in a narrow range of virial mass (Guo et al. 2010; Ferrero et al. 2012; Bullock & Boylan-Kolchin 2017; Sales et al. 2017).

A firm lower limit to that narrow range (⁠|$M_{200}\sim 10^9\, \mathrm{M}_\odot$|⁠) is suggested by the minimum virial temperature needed to allow hydrogen to cool efficiently, after accounting for the effects of an ionizing UV background (the ‘hydrogen cooling limit’, HCL; Efstathiou 1992; Quinn, Katz & Efstathiou 1996; Gnedin 2000; Okamoto & Frenk 2009; Benitez-Llambay & Frenk 2020). In quantitative terms, this implies that essentially all dwarfs with stellar mass, |$M_\star \lt 10^7\, \mathrm{M}_\odot$|⁠, no matter how faint they may be, should form in haloes with characteristic circular velocity² somewhere between |$20\, \mathrm{km\, s^{-1}}$| and |$40\, \mathrm{km\, s^{-1}}$| (Fattahi et al. 2018).

Together with the NFW profile mentioned above, the HCL minimum mass sets a velocity dispersion ‘floor’ for dark matter-dominated dwarfs of given radius, since more massive NFW haloes generally have, at all radii, larger circular velocities than less massive ones. This velocity ‘floor’ is, to first order, simply a fraction of the circular velocity (at that radius) of a halo at the HCL boundary. For example, for an HCL halo with |$V_\mathrm{mx}\approx 20\, \mathrm{km\, s^{-1}}$| (reached at a radius |$r_\mathrm{mx}= 3.4\, \mathrm{kpc}$|⁠), the circular velocity at |$300\, \mathrm{pc}$| is |$\approx 11\, \mathrm{km\, s^{-1}}$|⁠. This implies that a dwarf with 3D half-light radius, r_h ∼ 300 pc should have a line-of-sight (i.e. 1D) velocity dispersion, σ_los, in excess of |$11/\sqrt{3}\approx 6\, \mathrm{km\, s^{-1}}$|⁠. The same argument results in |$\sigma _{\rm los} \gtrsim 9\, \mathrm{km\, s^{-1}}$| for |$r_\mathrm{h}\sim 1\, \mathrm{kpc}$|⁠, and |$\sigma _{\rm los} \gtrsim 4\, \mathrm{km\, s^{-1}}$| for |$r_\mathrm{h}\sim 100\, \mathrm{pc}$|⁠. (These numbers assume a halo of average concentration at z = 0.)

There are a number of dwarfs in the Local Group that appear to violate these limits (see e.g. the compilation maintained by McConnachie 2012), a result which, taken at face value, would call for a review of some of the basic assumptions on which the above predictions are based. However, much of the wealth of available kinematic data on dwarfs concerns satellite galaxies in the Local Group, mainly those orbiting the Milky Way (MW) and the Andromeda (M31) galaxies (for a review, see Simon 2019).

The predictions described above do not apply to satellites, as tides arising from the MW and M31 may strip a significant fraction of the total dark matter content of a dwarf while leaving the stellar component relatively undisturbed (Peñarrubia, Navarro & McConnachie 2008; Errani, Peñarrubia & Tormen 2015; Sanders, Evans & Dehnen 2018). Translating the LCDM predictions described above to the realm of satellite galaxies thus requires a good understanding of how the dark and stellar components of dSphs evolve as a result of tidal effects.

This is an issue that has been addressed in earlier work, resulting in a number of conclusions and suggestions about how to interpret kinematic and photometric data on Local Group satellites in the context of LCDM. Two highlights of that work include the suggestion that (i) the stellar mass, size, and velocity dispersion of dSphs evolve along well-specified ‘tidal tracks’ which depend only on the total amount of mass lost from within the stellar half-light radius (Peñarrubia et al. 2008), and that (ii) the stellar components of ‘tidally limited’ satellites, defined as those whose stellar mass has been substantially reduced by tides, approach a ‘Plummer-like’ density profile shape roughly independent of the initial distribution of stars (Peñarrubia et al. 2009).

These conclusions, however, were based on simulations exploring a rather limited set of initial conditions, in terms of the assumed initial dSph structure, and also of the assumed radial segregation between stars and dark matter. Considering a broader range of possibilities for either may result in revised predictions that could impact, in particular, the properties of ‘tidally limited’ dwarfs. This is important, especially in light of the recent discovery of a population of dwarfs with unusually large sizes and low velocity dispersions, well below the limits mentioned earlier (see e.g. the Crater 2 and Antlia 2 dSphs; Torrealba et al. 2016; Caldwell et al. 2017; Torrealba et al. 2019).

Although it is tempting to associate such systems with tidally limited dwarfs (Frings et al. 2017; Sanders et al. 2018; Amorisco 2019), it is important to realize that the extreme estimated tidal losses put them in a regime hitherto unexplored in previous work. For example, Fattahi et al. (2018) argue that Crater 2 may be the result of a system that has lost more than |$99{{\ \rm per\ cent}}$| of its stars and has seen its velocity dispersion reduced by a factor of ∼5. As discussed by Errani & Navarro (2021), exploring this regime with N-body simulations requires resolving the progenitor subhalo with over 10⁷ particles, well beyond what has been achieved with cosmological hydrodynamical simulations.

As a result of these limitations, many questions regarding the structure and survival of tidally limited galaxies remain unanswered. Recent work, for example, has argued convincingly that NFW subhaloes are only very rarely fully disrupted by tides, almost always leaving behind self-bound dark remnants that are missed in cosmological simulations of limited resolution (Kazantzidis et al. 2004; Goerdt et al. 2007; Peñarrubia et al. 2010; van den Bosch et al. 2018). How do these results affect the structure and survival of the stellar components of such systems? Do the stellar components of dSphs also survive to some extent (giving rise to ‘micro-galaxies’, as argued in Errani & Peñarrubia 2020), and, if so, what are their properties? Should we expect them all to have similar density profiles? Do they evolve along well-defined ‘tidal tracks’ in luminosity, size, and velocity dispersion? Shoud we expect a large population of tidally limited dwarfs of extremely low surface brightness awaiting discovery?

We explore some of these issues here using a large suite of N-body simulations designed to study the tidal evolution of NFW haloes in the gravitational potential of a much more massive system, with particular emphasis on the regime of extreme tidal mass-loss. These simulations extend our earlier work on the subject (Errani & Navarro 2021), where we focused on the survival and structure of the dark component of the tidal remnant. Our emphasis here is on the evolution of a putative stellar component, assumed to be gravitationally unimportant compared to the dark matter. In the interest of simplicity, we only consider spherical, isotropic models in this contribution, but our approach should be relatively straightforward to extend to include modifications to these assumptions.

2 NUMERICAL METHODS

We use the set of high-resolution N-body simulations of the tidal evolution of NFW dark matter subhaloes introduced in Errani & Navarro (2021) and briefly summarized in this section.

2.1 Subhalo model

As the density profile of equation (1) leads to a diverging cumulative mass for r → ∞, we truncate the profile exponentially at |$10\, r_\mathrm{s}$|⁠. We generate 10⁷-particle equilibrium N-body realizations of NFW haloes with isotropic velocity dispersion, drawn from a distribution function computed numerically using Eddington inversion, following the implementation³ described by Errani & Peñarrubia (2020).

2.2 Host halo

2.3 Orbits

We focus on the evolution of subhaloes on eccentric orbits with peri-to-apocentre ratio of 1: 5. This value is consistent with the average peri-to-apocentre ratio of Milky Way satellite galaxies (Li et al. 2021; table 1), and close to the typical peri-to-apocentre ratio for dark matter substructures determined by van den Bosch et al. (1999) for isotropic distributions. As discussed by Errani & Navarro (2021), the orbital eccentricity affects mainly the rate of tidal stripping, but not the remnant structure.

If we assume an apocentric distance of |$r_\mathrm{apo}= 200\, \mathrm{kpc}$|⁠, a 1: 5 orbit implies a pericentric distance of |$r_\mathrm{peri}= 40\, \mathrm{kpc}$|⁠, and an orbital period of |$T_\mathrm{orb}= 2.5\, \mathrm{Gyr}$|⁠. Note that gravitational effects are scale-free, so all dimensional quantities are just quoted for illustration, and may be rescaled as needed.

We shall only explore the properties of subhaloes at their orbital apocentres, where they spend most of their orbital time, and when they are closest to equilibrium. Large deviations from equilibrium are expected near pericentre, as discussed by Peñarrubia et al. (2009), but we defer their study to future contributions.

We show in Appendix  B that all conclusions derived for the 1:5 orbit also hold for subhaloes on circular orbits, and for more radial orbits with peri-to-apocentre ratios of 1: 10 and 1:20.

2.4 Subhalo tidal evolution

All subhaloes have initial characteristic velocities substantially smaller than that of the host halo, i.e. V_mx0/V₀ ≲ 0.02. This choice allows us to neglect the effects of dynamical friction, and should be appropriate for studying the evolution of the subhaloes of faint and ultrafaint satellites of galaxies like the MW or M31.

The simulations explore a wide range of initial characteristic radii, r_mx0, chosen so that the ratio of initial subhalo crossing time, T_mx0, to the circular orbital time at the pericentre of the orbit, |$T_\mathrm{peri}=2\pi \, r_\mathrm{peri}/V_0$|⁠, lies in the range 2/3 < T_mx0/T_peri < 2.

The structural properties of bound subhalo remnants are determined by first (i) locating the subhalo centre through the shrinking spheres method (Power et al. 2003); then, (ii) assuming spherical symmetry, computing potential, and kinetic energies for N-body particles in a frame of reference comoving with the subhalo; and, then, (iii) removing unbound particles in this comoving frame. We iterate these steps until convergence.

2.5 Simulation code

The N-body models are evolved using the particle-mesh code superbox (Fellhauer et al. 2000). This code employs a high- and a medium-resolution grid of 128³ cells each, co-moving with the subhalo centre of density, with cell size Δx ≈ r_mx/128 and |$\approx 10\, r_\mathrm{mx}/ 128$|⁠, respectively. A third, low-resolution grid is fixed in space with a cell size of |$\approx 500\, \mathrm{kpc}/128$|⁠. The time integration is performed using a leapfrog-scheme with a time-step of Δt = min (T_mx, T_peri)/400. This choice of time-step ensures that, for NFW density profiles, at radii corresponding to the cell size of the highest resolving grid, circular orbits are still resolved by (at least) ≈16 steps.

The convergence tests of Errani & Navarro (2021, appendix A) for the set of simulations used in this work suggest that the structural parameters {r_mx, V_mx} of the N-body models may be considered unaffected by resolution artefacts as long as r_mx > 8Δx (grid resolution being the main limitation for this set of simulations), which translates to a minimum remnant bound mass fraction of M_mx/M_mx0 ≈ 1/300.

2.6 Stellar tracers

3 RESULTS

3.1 Tidal evolution of NFW subhaloes

To illustrate some of the general features of the evolution of the stellar and dark matter components of NFW subhaloes, we discuss here the results of a simulation of a system on an orbit with a pericentre-to-apocentre ratio of 1:5. The example subhalo chosen has an initial mass of |$M_\mathrm{mx0}=10^6\, \mathrm{M_{\odot }}$|⁠, an initial scale radius of |$r_\mathrm{mx0}=0.47\, \mathrm{kpc}$|⁠, and an initial maximum circular velocity of |$V_\mathrm{mx0}= 3.0\, \mathrm{km\, s^{-1}}$|⁠. We will refer to this model as halo h-i. The apocentre and pericentre radii are 200 and |$40\, \mathrm{kpc}$|⁠, respectively, and T_mx0/T_peri ∼ 0.9.

We use the same simulation to follow the evolution of two independent stellar tracers with different degrees of radial segregation relative to the dark matter. In one tracer, the initial star-to-dark matter radial segregation is R_h0/r_mx0 = 3/4; the second tracer is more deeply embedded inside the halo, with R_h0/r_mx0 = 1/4. The projected density profiles of these two systems, normalized to the same central density value, are shown in Fig. 1 with solid red and green lines, respectively. The dark matter is shown as well, with a blue line. (Note that the relative normalization of the density is arbitrary, as the stellar components are assumed to be ‘mass-less’.) We will refer to these two example stellar tracers as s-i and s-ii, for the extended and the deeply embedded tracer, respectively.

Fig. 2 shows snapshots of the dark matter (top row) and of each of the two embedded stellar tracers (middle and bottom rows) in the initial conditions (left-most column), and at different apocentric passages, chosen after 5, 10, 20, and 30 orbits.

As discussed in Section 2.4, tides gradually strip the subhalo, causing V_mx and r_mx to decrease, quite rapidly at first, but slowing down as the subhalo approaches the asymptotic remnant state, where T_mx ≈ T_peri/4. The effects of tidal stripping on this orbit are quite dramatic: the mass within r_mx declines to |$13{{\ \rm per\ cent}}$|⁠, |$6.3{{\ \rm per\ cent}}$|⁠, |$2.6{{\ \rm per\ cent}}$|⁠, and |$1.4{{\ \rm per\ cent}}$| of the initial value at each of the selected snapshots.

The evolution of the stellar tracers is qualitatively similar; their half-light radii decrease gradually as a result of tidal stripping, and the reduction in size decelerates as the tidal remnant stage is approached. Interestingly, the half-light radii of both stellar tracers appear to converge to the same value, despite the fact that they initially differed by a factor of 3. After 30 orbital periods, the difference in half-light radius is just |$\sim 10{{\ \rm per\ cent}}$|⁠. Furthermore, these half-light radii (shown by solid circles in the middle and bottom rows of Fig. 2) become comparable to the characteristic radius, r_mx, of the dark matter (solid circles in the top row of Fig. 2).

This is apparent in Fig. 1, where the dashed lines show the density profiles after 10 orbital periods. Solid circles on each curve mark the current value of the half-light radius. At the final time the radii are quite similar (as highlighted by the grey band), as discussed above. The total mass lost is, however, quite different. Only |$\sim 6{{\ \rm per\ cent}}$| of the initial characteristic dark matter mass remains bound at that time, compared with L/L₀ ≈ 0.24 and 0.01 for the more and less deeply embedded stellar component, respectively.

This discussion foretells some of our main results: ‘tidally limited’ systems, regardless of initial size, converge to a size set by the characteristic radius of the remnant subhalo. As we shall see below, these remnants also converge to the same velocity dispersion, which itself traces the characteristic velocity of the remnant. This implies that, in the ‘tidally limited’ regime, the properties of the stellar component are poor indicators of the initial mass/size/velocity of its progenitor: stars, unlike the dark matter, follow no unique ‘tidal tracks’ as they are stripped.

3.2 Dark matter evolution

As discussed in the previous section, the evolution of stars and dark matter are closely coupled. We study the evolution of the dark matter in a bit more detail next, before discussing the evolution of the stellar components in Section 3.3.

3.2.1 Circular velocity profile

As discussed by Errani & Navarro (2021), tides lead to a gradual change in the shape of the mass profile of an NFW subhalo as well as to a reduction in its characteristic size and circular velocity. We show this in Fig. 3, where we plot the evolution of the circular velocity profile V_c(r) = [GM(< r)/r]^1/2 of the subhalo shown in Fig. 2.

The distance between adjacent circular velocity curves decreases with decreasing remnant mass, as the subhalo asymptotically approaches the asymptotic remnant state, where its crossing time is roughly a quarter of the host halo crossing time at pericentre, T_mx ≈ T_peri/4.

3.2.2 Binding energy

Similar results are obtained when using initial orbital times instead of binding energy, as discussed in Appendix  C.

The energy distribution of an isolated, exponentially truncated NFW cusp with isotropic velocity dispersion may be computed analytically, and is shown by a black curve in Fig. 6. As is clear from this figure, the remnant final energy distribution quickly approaches that of an exponentially truncated cusp. In contrast to the initial NFW profile (shown, for reference, as a dashed curve with arbitrary normalization), which has most mass at low binding energies, the exponentially truncated cusp has a well-defined peak energy, beyond which the energy distribution drops steeply.

Once the subhalo has been stripped to less than 10 per cent of its initial mass, the energy distribution corresponding to an exponentially truncated density cusp provides a good description for the energy distribution measured for the N-body remnant. Deviations at small values of |${\mathcal {E}}_\mathrm{f}$| are expected due to limitations in numerical resolution introduced by the finite code grid size, Δx; energy values corresponding to |$E<\Phi(8\Delta x)$| are shown in grey in Fig. 6.

3.2.3 Initial-to-final energy mapping

The initial and final binding energies of particles that remain bound are closely related. We illustrate this dependence in Fig. 7, which shows final energies as a function of their initial values for the subhalo highlighted in Fig. 2, after |$10\, T_\mathrm{orb}$|⁠. The relation is almost linear for most energy values, except for particles initially less bound than the ‘peak’ energy |${\mathcal {E}}_\mathrm{mx,t}$|⁠.

3.3 Evolution of the stellar component

The results of the previous subsection enable a thorough description of the tidal evolution of the dark matter component of an NFW subhalo. Under the assumption that stars may be treated as tracers of the gravitational potential, their evolution would just follow that of a subset of the dark matter, defined mainly by their initial binding energy distribution.

Let us consider again the example of Fig. 2. The initial energy distributions of stars and dark matter are shown with solid curves in Fig. 9. For exponential profiles, the most bound regions are well approximated by power laws that peak at |${\mathcal {E}}_\star$| (a characteristic value which depends on the radial segregation of stars and dark matter) and are sharply truncated beyond |${\mathcal {E}}_\star$|⁠.

The initial energy distribution of particles that remain bound after 10 orbital periods are shown by the lines connecting circles of the corresponding colour. The black dashed curves overlapping each of these lines indicate the result of applying the ‘filter function’ introduced in Section 3.2.2 (equation 9). Note that the same filter, applied to either stars or dark matter, yields an excellent description of the initial energy distribution of the particles that remain bound after |$10\, T_\mathrm{orb}$|⁠.

This has a number of important implications. It suggests that (i) in practice the only parameter that matters when determining the probability that a given particle will remain bound (or be stripped) is its initial binding energy. It also implies that (ii) stellar components whose initial energy distributions peak beyond the truncation energy, |${\mathcal {E}}_\mathrm{mx,t}$|⁠, will share the same outermost energy distribution as the dark matter remnant. Finally, (iii) it suggests that the final structure of the bound stellar remnant will depend on how stars populate energies more bound than |${\mathcal {E}}_\mathrm{mx,t}$| in the initial system; in other words, the ‘tidal tracks’ and final structure of the stellar remnant should depend sensitively on the stars’ initial density profile and radial segregation relative to the dark matter. We explore these implications in more detail next.

3.3.1 Initial and final stellar remnant structure

Fig. 10 shows two examples, one with α = β = 3 (shown in magenta) and another with α = β = 6 (shown in orange). The value of |${\mathcal {E}}_\mathrm{s}$| is chosen so that both models have the same initial 2D half-light radius, R_h0/r_mx0 = 1/4 (isotropic solutions for different initial half-light radii are discussed in Appendix  E). As α = β, the scale energy and peak energy coincide; i.e. |${\mathcal {E}}_\star = {\mathcal {E}}_\mathrm{s} \approx 1/3$|⁠. The α = β = 3 model spans a much wider range of energies than the α = β = 6 model (see top panel), and thereby also samples more energies that probe deeper into the dark matter density cusp. As a consequence, the resulting surface brightness profile of the α = β = 3 model has a smaller core radius (defined as the radius where the projected stellar density drops by a factor of 2 from its central value; i.e. Σ(R_c) = Σ(0)/2) than the α = β = 6 model. The corresponding 2D and 3D density profiles are shown in the bottom panels of Fig. 10. For reference, the ratio R_c/R_h between core and half-light radius in the initial conditions is 0.5 and 0.9 for α = β = 3 and α = β = 6, respectively.

Many different density profiles may be approximated by varying α and β in equation (13). For example, a Plummer model may be approximated quite well by α = 14, β = 0.45, whereas exponential profiles are well fit by α = 3.3, β = 4.0. Indeed, the density profile for the α = β = 3 case is almost indistinguishable from an exponential profile over two decades in radius, as may be seen in the bottom panels of Fig. 10.

As expected, the final and initial energy distribution of the stars are highly dependent on each other. For example, we show in Fig. 11 the final energy distributions for the α = β = 3 (magenta) and α = β = 6 (orange) stellar tracers in the ‘tidally limited’ regime; i.e. after substantial mass-loss. Various snapshots are shown at apocentric passages and span a large range of bound mass fraction: 1/100 < M_mx/M_mx0 < 1/10. Energies |${\mathcal {E}}_\mathrm{f} = 1-E/\Phi _{0}$| are normalized by their instantaneous potential minimum (approximated as |$\Phi _0 \approx \Phi _\mathrm{0,asy} \approx -3.35\, V_\mathrm{mx}^2$|⁠). The shape of the energy distributions converges rapidly (as shown in Fig. 6 for dark matter alone), and little evolution in the energy distribution shape is seen for remnant masses smaller than |$0.1\, M_\mathrm{mx0}$|⁠.

The final energy distributions of the remnants may also be computed using the ‘filter’ function of equation (9) to select bound particles from the initial conditions, combined with the initial-to-final energy mapping given by equation (12). The result of this calculation is shown by the solid black lines in Fig. 11 and are seen to approximate well the results of the simulations. A step-by-step description of how to implement this calculation is given in Appendix  G.

The different initial energy distributions give rise to different density profiles for the stellar remnants, depending on the values of α and β chosen. We show this in Fig. 12, where we plot the 3D density profiles corresponding to the energy distributions shown in Fig. 11. Radii are scaled to the instantaneous characteristic size of the underlying dark matter halo, r_mx(t). All profiles are sharply truncated beyond r_mx, and have similar outer slopes (consistent with the findings of Peñarrubia et al. (2009), who studied tidal stripping of stellar King models deeply embedded within NFW subhaloes). The inner regions, however, differ, retaining clean memory of their initial profiles. As in the initial conditions, the α = β = 6 remnant profile has a larger ratio of core to half-light radius (R_c/R_h ∼ 0.8) than the α = β = 3 profile (R_c/R_h ∼ 0.6).

3.3.2 Tidal tracks

The results presented in Fig. 9 suggest that luminosity and size of the stellar components should remain largely unaffected by tides until the truncation in energy starts to overlap the energy distribution of the stars. Using the notation illustrated in that figure, this occurs when |${\mathcal {E}}_\mathrm{mx,t}$| becomes comparable or smaller than |${\mathcal {E}}_\star$|⁠. We show this in Fig. 13, where we plot the evolution of the stellar projected half-light radius R_h, as well as of the central line-of-sight velocity dispersion σ_los and luminosity (i.e. stellar mass), L, of the two stellar populations highlighted in the example shown in Fig. 2.

The evolution is plotted as a function of the tidal truncation energy, |${\mathcal {E}}_\mathrm{mx,t}$|⁠, which is a proxy for the total mass-loss due to tides (see Fig. 5). The smaller |${\mathcal {E}}_\mathrm{mx,t}$| the higher the mass-loss, so tidal evolution runs from right to left in Fig. 9. The location of the peak |${\mathcal {E}}_\star$| for each of the two initial stellar energy distributions is indicated by the vertical arrows.

As expected, if the tidal truncation energy does not overlap the stellar energy distribution (i.e. if |${\mathcal {E}}_\mathrm{mx,t}\gt {\mathcal {E}}_\star$|⁠), the size of the stellar component remains largely unaffected by tides. However, once |${\mathcal {E}}_\mathrm{mx,t}$| drops below |${\mathcal {E}}_\star$|⁠, the stellar components rapidly lose mass, decrease in size, and gradually become kinematically colder. As expected, the more extended stellar component (R_h0/r_mx0 = 3/4) is more readily affected than the more deeply segregated component with R_h0/r_mx0 = 1/4.

Interestingly, once |${\mathcal {E}}_\mathrm{mx,t}$| becomes much smaller than |${\mathcal {E}}_\star$|⁠, the stellar half-light radius converges to the dark matter characteristic radius, r_mx, regardless of its initial value. In addition, the velocity dispersion becomes directly proportional to the dark matter characteristic velocity, V_mx. In other words, the stellar components of these ‘tidally limited systems’ become direct tracers of the characteristic parameters of the tidal remnant, almost regardless of their initial radial segregation (see also Kravtsov 2010).

This is a general result, which applies to all systems in the ‘tidally limited’ regime. We illustrate this in Fig. 14 where we plot the ‘tidal tracks’ of α = β = 3 systems with different initial radial segregation; i.e. R_h0/r_mx0 = 1/2, 1/4, 1/8, and 1/16, each shown with different symbols (Appendix  F shows the same tidal tracks in units of remnant bound mass fraction M_mx/M_mx0). As the subhalo gets stripped, its characteristic radius and velocity (r_mx and V_mx) follow the track indicated by the dashed line. The 3D half-light radius and velocity dispersion of the stellar components (the latter multiplied by |$\sqrt{3}$| to convert it to a circular velocity) evolve as well, generally getting smaller and colder as stripping progresses.

Only the most embedded tracers ‘puff up’ initially, increasing their size⁶ as they get kinematically colder, but the net increase is small. Interestingly, once the stellar component parameters reach the tidal track of the dark matter remnant, they follow it closely. Qualitatively, this result applies to all the values of α and β we have tried, and it therefore appears of general applicability: the size and velocity dispersion of a tidally limited stellar system traces closely the r_mx and V_mx of the dark remnant.⁷

This result has a couple of interesting implications. One is that there is no unique ‘tidal track’ describing the evolution of a stellar system in the size-velocity plane, as had been suggested by Peñarrubia et al. (2008), i.e. the tidal track depends on the initial distribution of stars within the subhalo. With hindsight, it is easy to trace that conclusion to the fact that those authors tested a single stellar model (a King profile) and a rather limited range of radial segregation. In fact, Fig. 14 shows clearly that the same stellar tidal remnant may be produced from very different initial radial segregations, each following different tracks.

A second implication is that the tidal track of a subhalo (i.e. the dashed line in Fig. 14) sets a firm lower limit on the velocity dispersion of an embedded stellar remnant of given size. (Or, equivalently, an upper limit to the size of an embedded system of given velocity dispersion.)

We discuss next an application of this argument to the interpretation of the sizes and velocities of dwarf galaxies in the Local Group.

3.4 Application to Local Group dwarfs

We now use the results of the previous subsections to interpret available structural data on Local Group dwarfs, in the context of LCDM. As discussed in Section 1, the assumption of a minimum mass for galaxy formation, such as that expected from the hydrogen cooling limit, together with the mass–concentration relation expected for NFW haloes in LCDM, yield strong constraints on the size and velocity dispersion of dwarf galaxies.

Following Fattahi et al. (2018), the grey band in Fig. 15 indicates the circular velocity profiles of NFW haloes with V_mx in the range 20–40 km s⁻¹, the virial mass range expected to host all isolated dwarfs with |$M_\star \lt 10^7\, \mathrm{M}_\odot$|⁠. The solid orange line indicates the r_mx–V_mx relation expected from the LCDM mass–concentration relation (Ludlow et al. 2016). The accompanying error band delineates the one- and two-sigma scatter inferred from cosmological N-body simulations.

For dSphs, the velocity dispersion of the stars may be used to estimate the circular velocity, V_h, at the 3D stellar half-light radius, r_h, using the Wolf et al. (2010) mass estimator.⁸ While not exact, we expect these estimates to be accurate to about 0.1 dex (see also Walker et al. 2009; Campbell et al. 2017; Errani et al. 2018). Dwarfs unaffected by tides should lie on the grey area if they follow LCDM predictions. Encouragingly, Local Group isolated (field) dwarfs, shown by magenta symbols in Fig. 15, seem to be in reasonable agreement with that prediction.

Systems below the grey band are, in the LCDM context, typically interpreted as satellites whose velocity dispersions and sizes have been affected by tides, bringing them below the grey band. Indeed, no dwarf with |$M_\star \lt 10^7\, \mathrm{M}_\odot$| should fall above the grey band, according to LCDM. The ones that seem to deviate from that precept, such as Bootes 2 and Carina 3, have fairly large uncertainties in their circular velocity estimates, mainly due to the extremely small number of stars with measured radial velocities. Furthermore, Ji et al. (2016) argue that the velocity dispersion of Bootes 2 is likely overestimated due to the presence of binaries, and that the reported velocity dispersion is best seen as an upper limit. It remains to be seen whether these galaxies are actually in conflict with LCDM once their velocity dispersions are better constrained.

Galaxies below the grey band are all satellites, a requisite for a tidal interpretation of their location in the size-velocity plane. However, tides cannot shift them to arbitrarily low velocities without affecting their size, as shown by our discussion of the tracks in Fig. 14. Indeed, the tidal track of the halo with the minimum mass needed to form a galaxy should provide firm limits on the size and velocity of equilibrium ‘tidally limited’ dwarfs. This may be seen as a lower limit on the velocity at given size, or, equivalently, as an upper limit to the size given a velocity.

This limit is indicated by the blue line in Fig. 15, which outlines the tidal track of a V_mx = 20 km s⁻¹ halo of average concentration. This restricts the loci of tidally limited dwarfs to the region shaded in blue: no dwarf should populate the region below the blue line.

Interestingly, there are a number of satellites that seem to violate this condition. These are galaxies of unusually low velocity dispersion for their size, like the MW ‘feeble giant’ satellites Crater 2 and Antlia 2, as well as the M31 satellites And 19, And 21 and And 25. Our results imply that tidal effects are not a viable explanation for these galaxies in the context of LCDM.

These conclusions are consistent with those of Sanders et al. (2018), who reported ‘tension’ when trying to account for the size of Crater 2 in tidally limited NFW haloes, and argued that the tension may be alleviated if the progenitor subhalo had a sizable constant-density ‘core’. Sanders et al. (2018) reached these conclusions considering also flattened models in their analysis, which suggests that angular momentum plays a subdominant role in determining the remnant structure of dispersion-supported tidally stripped systems, see Appendix  A.

On the other hand, Frings et al. (2017) reports simulations that can apparently produce satellites as large and as kinematically cold as Crater 2. However, these results are likely affected by numerical limitations. For example, the closest ‘Crater 2 analogue’ that they report corresponds to a halo whose maximum circular velocity has been reduced to nearly 1/10th of its original value. As discussed by Errani & Navarro (2021), adequately resolving the tidal remnant in that case is extremely difficult, and would require many more particles and much better spatial resolution than used by Frings et al. (2017). In addition, some of the Frings et al. (2017) satellites have density profile cusps somewhat shallower than NFW, hindering a detailed comparison with our results.

Finally, Amorisco (2019) argued that transients induced by tides may be the cause of the large radii and low velocity dispersions of these objects. This interpretation could in principle be checked by looking for other signs of ongoing tidal disturbance, such as a velocity gradient across the system (Li et al. 2021), extended tidal tails, or ‘bumps’ in the surface density profile caused by escaping stars (Peñarrubia et al. 2009). We have not considered such transients in the analysis presented here, which only considers the (equilibrium) structure of galaxies near apocentre.

To summarize, galaxies below the blue tidal track shown in Fig. 14 present a serious challenge to the tidal interpretation suggested by models based on the LCDM scenario. Possible resolutions of this challenge include: (i) revisions to the observed photometric or kinematic parameters of these galaxies, many of which are extremely challenging to study due to their extremely low surface brightness; (ii) the possibility that these galaxies are not bound, equilibrium systems, but are in the process of being tidally stripped, (iii) that these systems formed in haloes with masses substantially below the hydrogen cooling limit, or (iv) that the dark matter halo of these galaxies deviates from the NFW shape, perhaps signalling the effects of baryons on the inner cusp of the halo even in ultrafaint galaxies,⁹ or, perhaps more intriguingly, effects associated with the intimate nature of dark matter, such as finite self-interactions, or other such deviations from the canonical LCDM paradigm.

4 SUMMARY AND CONCLUSIONS

We have studied the tidal evolution of dark matter-dominated satellites embedded in LCDM dark matter subhaloes as they orbit in the potential of a massive galaxy. Our study assumes that subhaloes may be approximated by spherically symmetric, isotropic, cuspy Navarro–Frenk–White density profiles. We further assume that subhalo masses are much smaller than the host halo mass, so as to neglect the effects of dynamical friction. As angular momentum plays only a subordinate role in determining if a dark matter particle gets stripped or remains bound (see Appendix  A), we expect our results to hold not only for the isotropic models studied in this work, but also for mildly anisotropic models. Our main findings may be summarized as follows.

• Tides lead NFW subhaloes to evolve following ‘tidal tracks’ that describe the changes in characteristic size and velocity as tidal mass-losses accumulate. These tidal tracks lead to a well-defined asymptotic tidal remnant with characteristic crossing time set by the orbital time at pericentre, T_mx ≈ T_peri/4. The asymptotic structure of the tidal remnant is well approximated by an exponentially truncated cusp.

• Tidal evolution can be conveniently modelled in terms of binding energy; tides gradually truncate subhaloes energetically, systematically removing particles with low binding energies. The bound remnant consists of particles initially more strongly bound than the energy threshold imposed by tides.

• For the majority of particles, the initial binding energy alone is sufficient to determine when and whether a particle will be stripped or not. Initial and final binding energies are strongly correlated, in a way that allows the structure of the remnant to be inferred as a function of tidal mass-loss. These relations hold regardless of the orbital eccentricity of the subhalo.

• If gravitationally unimportant, the stellar components of tidally stripped NFW subhaloes evolve according to how stars populate the initial energy distribution of the subhalo. The same gradual energy truncation applies to dark matter and stars, independent of the initial density structure or radial segregation of the stellar component.

• The structure of stellar remnants retain memory of their initial structure, and, in particular, of how stars populate the most bound energy levels of the inner dark matter cusp. Their evolution is thus sensitively dependent on their initial structure and radial segregation. Unlike the dark matter, stellar components do not follow unique ‘tidal tracks’ in stellar mass, size, and velocity dispersion.

• The stellar component may be completely dispersed by tides if the stars do not populate the most bound energy levels of the dark matter cusp, making the potential existence of ‘micro-galaxies’ critically dependent on how stars populate the most bound energy levels of an NFW subhalo.

• ‘Tidally limited’ satellites, defined as those which have lost a substantial fraction of their initial stellar mass, have radii and velocity dispersions that trace directly the characteristic radius and velocity of the subhalo remnant.

Most Local Group dwarfs have structural parameters consistent with them being embedded in cuspy NFW haloes, but there are also a number whose dynamical masses are below what is expected from LCDM. The properties of such systems, all of which are M31 or MW satellites, are usually assumed to result from the effects of tidal stripping. Our work sheds further insight into this interpretation.

Indeed, coupled with the assumption of a minimum halo mass, such as that suggested by hydrogen cooling limit considerations, our results place strong constraints on the size and velocity dispersion of tidally limited systems embedded in the remnants of cuspy NFW subhaloes. These constraints may be expressed as a firm lower limit on the velocity dispersion of an embedded stellar remnant of given size, or, equivalently, as an upper limit to the size of an embedded system of given velocity dispersion. Inspection of available data for dwarf galaxies in the Local Group, however, reveals a number of ultrafaint satellites that breach these limits.

The findings we report here imply that such systems cannot be understood as equilibrium systems embedded in the tidal remnants of cuspy NFW subhaloes. This presents a serious challenge to our understanding of the formation and evolution of these unusual galaxies in LCDM. Reconciling systems like Crater 2, Antlia 2, or And 19 with LCDM seems to require substantial revision to at least one of the assumptions of this work. Either the reported observational parameters of such systems are in substantial error, or the systems are far from dynamical equilibrium, or they inhabit subhaloes whose density profiles differ significantly from NFW.

Neither alternative is particularly attractive in the context of LCDM. It may indicate that baryons may affect the inner halo cusp even in extremely faint dwarfs or, more intriguingly, may indicate effects associated with the intimate nature of the dark matter, such as finite self-interactions or other such deviations from the canonical LCDM paradigm. Resolving this challenge will likely require a concerted numerical and observational effort designed to improve our understanding of the formation and evolution of some of the faintest galaxies known.

ACKNOWLEDGEMENTS

RE acknowledges support provided by a CITA National Fellowship. Both RE and RI acknowledge funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No. 834148). This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES

APPENDIX A: STRIPPING AND ANGULAR MOMENTUM

Tides gradually truncate NFW subhaloes in energy, stripping the least-bound particles first, as shown in Fig. 4. The truncation in energy is relatively sharp, and may be approximated using the ‘filter function’ of equation (9). This function gradually truncates the initial energy distribution beyond some energy |${\mathcal {E}}_\mathrm{mx,t}$|⁠. While in Fig. 4 almost all particles that fall ‘to the left’ of the truncation energy |${\mathcal {E}}_\mathrm{mx,t}$| remain bound, there is a tail of particles towards less-bound energies where energy alone is not sufficient to tell whether these particles will be stripped, or remain bound. This is shown in the left-hand panel of Fig. A1.

In the highlighted range of energies (black dashed selection), |$0.25 \le {\mathcal {E}}\le 0.3$|⁠, half of particles get stripped, while the other half remains bound. The right-hand panel of Fig. A1 shows the (initial) angular momentum distributions of particles in the selected energy range, expressed in units of the angular momentum of a circular orbit |$L_\mathrm{c}({\mathcal {E}})$| at energy |${\mathcal {E}}$|⁠. The (initial) angular momentum distributions are remarkably similar. Only particles on nearly circular orbits exhibit a mild tendency to remain bound. Indeed, whereas |$\sim 48{{\ \rm per\ cent}}$| of particles with L/L_c < 0.9 remain bound (and |$\sim 52{{\ \rm per\ cent}}$| are stripped), |$\sim 61{{\ \rm per\ cent}}$| of particles with L/L_c ≥ 0.9 remain bound (and |$\sim 39{{\ \rm per\ cent}}$| are stripped). This difference, however, affects fewer than |$\sim 20{{\ \rm per\ cent}}$| of all particles in that energy bin.

While all models studied in this work are based on spherical, isotropic initial conditions that can be expressed as |$f({\mathcal {E}})$| distribution functions, because of the small effect of angular momentum on determining if a particle gets stripped, we expect that our conclusions hold also for mildly anisotropic models (cast, e.g. as |$f({\mathcal {E}},L_z)$| distribution functions).

APPENDIX B: ORBITAL ECCENTRICITY

Section 3.2 discusses tidal stripping of NFW subhaloes on an orbit with pericentre-to-apocentre ratio of |$1{\, :\, }5$|⁠. This Appendix extends the analysis to circular (⁠|$1{\, :\, }1$|⁠) and more radial orbits (⁠|$1{\, :\, }10$|⁠, |$1{\, :\, }20$|⁠). The circular orbit has a radius of |$40\, \mathrm{kpc}$| and an orbital period of |$T_\mathrm{orb}= 1.12 \, \mathrm{Gyr}$|⁠. All eccentric orbits have an apocentre distance of |$r_\mathrm{apo}=200\, \mathrm{kpc}$|⁠, allowing the subhalo to reach equilibrium between subsequent pericentre passages. The |$1{\, :\, }10$| and |$1{\, :\, }20$| orbits consequently have pericentre distances of |$r_\mathrm{peri}= 20\,$| and |$10\, \mathrm{kpc}$|⁠, and orbital periods of 2.31 and |$2.25\, \mathrm{Gyr}$|⁠, respectively. Subhalo models were chosen from the Errani & Navarro (2021) set of simulations and, on the selected orbit, were stripped to less than one per cent of their initial mass. For the circular orbit, we used the same subhalo as for the |$1{\, :\, }5$| of Section 3.2, with |$r_\mathrm{mx0}= 0.48\, \mathrm{kpc}$| and |$V_\mathrm{mx0}= 3.0\, \mathrm{km\, s^{-1}}$|⁠. For the |$1{\, :\, }10$| and |$1{\, :\, }20$| radial orbits, we used subhaloes with |$r_\mathrm{mx0}= 0.39\, \mathrm{kpc}$|⁠, |$V_\mathrm{mx0}= 3.3\, \mathrm{km\, s^{-1}}$| and |$r_\mathrm{mx0}= 0.29\, \mathrm{kpc}$|⁠, |$V_\mathrm{mx0}= 3.8\, \mathrm{km\, s^{-1}}$|⁠.

The left-hand column of Fig. B1 shows the truncation in energy of subhaloes on these three different orbits (as shown in Fig. 4 for a |$1{\, :\, }5$| orbit). The filter function of equation (9) (two examples are shown as black dashed curves in each panel) provides an approximative description for the shape of the energy truncation regardless of orbital eccentricity. The small deviations seen may contribute to the minor eccentricity dependence of tidal evolutionary tracks as shown in fig. 6 of Errani & Navarro (2021).

The right-hand column of Fig. B1 shows that the initial-to-final energy mapping (as shown in Fig. 8 for a |$1{\, :\, }5$| orbit) is well-described by equation (12), regardless of eccentricity.

APPENDIX C: ORBITAL PERIODS

As M_mx/M_mx0 decreases, the bound remnant consists of particles which in the initial conditions were located at decreasing distance to the subhalo centre. Initial radii alone, however, do not determine whether a particle will remain bound or not. Indeed, at fixed radius in the left-hand panel of Fig. C1, some particles get stripped, but others remain bound.

In contrast, the truncation in the initial radial orbital period T_r imposed by tides is steep (see right-hand panel of Fig. C1), selecting most particles with radial orbital periods shorter than the truncation period. Tidal evolution stalls once the relaxed remnant has a crossing time of roughly one quarter of the host halo crossing time at pericentre, 2πr_mx/V_mx ≈ T_peri/4. In the initial conditions, this corresponds to a selection of particles with orbital periods shorter than ≈T_peri/16.

APPENDIX D: RADIUS-DEPENDENT ENERGY DISTRIBUTION

No particles located at a radius r may have an energy lower than the potential ϕ(r) = 1 − Φ(r)/Φ₀ at that radius, shown as a solid black line. The dark matter (NFW) particles span a wide range in radii and energies, with substantial mass located at high binding energies within the central dark matter density cusp. In contrast, the initial energy distribution of stellar tracers is more localized, i.e. the range of radii and energies where most particles are located is narrower.

The average energy at a fixed radius for dark matter particles is shown as a dashed curve in all panels, while the average energy of each stellar tracer is shown using a dotted line. For the more deeply embedded stellar system, at the half-light radius R_h (corresponding approximately to the radius where most particles are located, see colour coding), stellar particles are on average more strongly bound than dark matter. Hence tidal stripping will initially preferentially remove dark matter. This causes the mass enclosed within the stellar half-light radius to drop (and thereby the stellar velocity dispersion to decrease), and the half-light radius to expand during relaxation. On the other hand, for the more extended stellar tracer, at the half-light radius, stars and dark matter have similar binding energies. As tides strip the system, both dark matter and stars are hence lost from the beginning, and R_h drops, consistent with the models shown in Fig. F1.

APPENDIX E: ISOTROPIC STELLAR MODELS

This appendix presents the structural properties of spherical stellar systems with isotropic velocity dispersion, embedded in NFW dark matter haloes, with energy distribution |$\mathrm{d}N_\star / \mathrm{d}{\mathcal {E}}$| given by equation (13). Fig. E1 shows the energy distributions (top panel), logarithmic slopes of the surface brightness profiles (2nd row), surface brightness profiles (3rd row), and line-of-sight velocity dispersion profiles (bottom row) for models with α = β = 3 (left-hand column), and α = β = 6 (right-hand column). Profiles for four different segregations, R_h/r_mx = 1/2 (red), 1/4 (orange), 1/8 (light blue), and 1/16 (dark blue), are shown.

Structural properties are computed directly from the distribution function |$f({\mathcal {E}}) \propto (\mathrm{d}N_\star / \mathrm{d}{\mathcal {E}})/ p({\mathcal {E}}) ~$|⁠, where |$p({\mathcal {E}})$| is the phase-space volume accessible in an NFW potential to a particle with energy |${\mathcal {E}}$| (i.e. the density of states).

The more deeply embedded the stellar profile is within its dark matter halo, the further the maximum of the energy distribution is shifted towards high binding energies. The range of energies selected by the α = β = 3 model is much wider than that of the α = β = 6 model. For equal values of R_h/r_mx, the energy |${\mathcal {E}}_\star$| at which the distribution peaks is near identical between the two models.

The exact shape of the stellar density profile depends on where it is embedded inside the dark matter halo. α = β = 3 models approximately resemble exponential density profiles, while α = β = 6 profiles have larger core radius (i.e. a radius R_c so that Σ(R_c) = Σ(0)/2) at equal half-light radius.

The more deeply embedded a stellar tracer is within its dark matter halo, the lower is its central velocity dispersion relative to the subhalo characteristic velocity V_mx (bottom panel). A constant density tracer of infinite extent has a 3D velocity dispersion profile equal to |$V_\mathrm{esc}(r)/\sqrt{2}$|⁠, hence, for reference, the escape velocity profile of an NFW subhalo is shown in the bottom panel of Fig. E1 (re-scaled by a factor of |$\sqrt{3}$| to compare against the 1D line-of-sight velocity profiles shown).

APPENDIX F: TIDAL TRACKS

Structural parameters of dwarf spheroidal galaxies evolve differently under tidal stripping depending on how deeply embedded they are within the dark matter halo, and which energies they populate within the halo. Fig. F1 shows the evolution of half-light radius R_h (top panel), line-of-sight velocity dispersion σ_los (2nd panel from top), luminosity L (3rd panel from top) and mass-to-light ratio averaged within the half-light radius Υ = L/[2M(< R_h)], as a function of remnant dark matter mass fraction M_mx/M_mx0. The evolution is shown for stellar tracers with initial half-light radii of R_h0/r_mx0 = 1/2 (red), 1/4 (orange), 1/8 (light blue), 1/16 (dark blue). Each point corresponds to one apocentre snapshot, taken from the set of simulations described in Section 2. The evolution of the subhalo characteristic radius r_mx and velocity V_mx are shown using solid black curves. Stellar tracers with α = β = 3 and α = β = 6 (equation 13) energy distribution are shown in the left-hand column and right-hand column, respectively. The full initial profiles for these tracers are shown in Fig. E1.

In the ‘tidally limited’ regime, the stellar half-light radii R_h trace the dark matter characteristic size r_mx (for the α = β = 3 and α = β = 6 models, we find in the tidally limited regime |$R_\mathrm{h}\approx 0.93\, r_\mathrm{mx}$| and |$R_\mathrm{h}\approx 1.26\, r_\mathrm{mx}$|⁠, respectively). Similarly, the stellar velocity dispersion σ_los evolves parallel to the dark matter characteristic velocity V_mx (for α = β = 3 and α = β = 6, we find |${\sigma _\mathrm{los}}\approx 0.70\, V_\mathrm{mx}$| and |${\sigma _\mathrm{los}}\approx 0.76\, V_\mathrm{mx}$|⁠). Once a stellar tracer has been stripped to the size of R_h ≈ r_mx, it becomes impossible to reconstruct its original extent from structural properties alone.

For order-of-magnitude estimates, a simple power-law approximation to this integral may be used for heavily stripped systems. For NFW (equation 1) density profiles and those of an exponentially truncated cusp (equation 3), the energy distribution |$\mathrm{d}N_\mathrm{DM}/ \mathrm{d}{\mathcal {E}}$| is well-approximated by a power law of slope ≈1 towards the most-bound energies, i.e. |$\mathrm{d}N_\mathrm{DM}/ \mathrm{d}{\mathcal {E}}\propto {\mathcal {E}}$| for |${\mathcal {E}}\ll 1$|⁠. Assuming that the truncation happens sharpy at an energy |${\mathcal {E}}_\mathrm{mx,t} \ll 1$|⁠, and that the density profile has converged to that of an exponentially truncated cusp, this gives |$M_\mathrm{mx}\propto M \approx \int _0^{{\mathcal {E}}_\mathrm{mx,t}} \mathrm{d}{\mathcal {E}}~ \mathrm{d}N_\mathrm{DM}/\mathrm{d}{\mathcal {E}}\propto {\mathcal {E}}_\mathrm{mx,t}^2$|⁠. Similarly, for stellar tracers with an energy distribution parametrized through equation (13) and |${\mathcal {E}}_\mathrm{mx,t} \ll {\mathcal {E}}_\star$|⁠, |$L \approx \int _0^{{\mathcal {E}}_\mathrm{mx,t}} \mathrm{d}{\mathcal {E}}~ {\mathcal {E}}^\alpha \propto {{\mathcal {E}}_\mathrm{mx,t}}^{\alpha +1}$|⁠. Consequently L ∝ M^{(α + 1)/2} for highly stripped stellar systems.

APPENDIX G: EMPIRICAL MODEL FOR EVOLVED STELLAR PARAMETERS

In this appendix, we discuss how to construct the relaxed energy distribution of a tidally stripped (stellar) tracer in a cold dark matter subhalo. For this, we make use of the observation that stellar energy distributions are subject to the same tidal truncation as the underlying NFW dark matter halo (Section 3.3), and of the energy mapping between the initial NFW halo and the stripped, relaxed system discussed in Section 3.2.3.

Given an initial stellar energy distribution, |$\mathrm{d}N_\star / \mathrm{d}{\mathcal {E}}|_\mathrm{i}$|⁠, as well as the desired remnant mass fraction of the underlying NFW dark mater halo, M_mx/M_mx0, the procedure involves the following steps:

• For the desired remnant bound mass fraction M_mx/M_mx0 of the underlying dark matter subhalo, the tidal truncation energy |${\mathcal {E}}_\mathrm{mx,t}$| is computed from equation (10).

• The energy distribution |$\mathrm{d}N_\star / \mathrm{d}{\mathcal {E}}|_\mathrm{i,t}$| of those particles in the initial conditions which remain bound at the given tidal truncation energy |${\mathcal {E}}_\mathrm{mx,t}$| follows from applying the filter of equation (9) to the initial energy distribution.

• Finally, we need to model the relaxation process to pass from the initial truncated energy distribution |$\mathrm{d}N_\star / \mathrm{d}{\mathcal {E}}|_\mathrm{i,t}$| to the relaxed, equilibrium system |$\mathrm{d}N / \mathrm{d}{\mathcal {E}}|_\mathrm{f}$|⁠. To take into account the scatter¹⁰ in the energy mapping shown in the top panel of Fig. 7, we compute the final energy distribution through the convolution

  $$\begin{eqnarray*} \left. \frac{\mathrm{d}N_\star }{\mathrm{d}{\mathcal {E}}} \right|_\mathrm{f} = \int _{0}^{1} \mathrm{d}{\mathcal {E}}_\mathrm{i} ~ \left. \frac{\mathrm{d}N_\star }{\mathrm{d}{\mathcal {E}}} \right|_\mathrm{i,t} ~ \mathrm{Lognormal}\left(\bar{{\mathcal {E}}}_\mathrm{f}({\mathcal {E}}_\mathrm{i}) , \mathrm{0.03\, dex} \right), \end{eqnarray*}$$

  (G1)

where |$\mathrm{Lognormal}\left(\bar{{\mathcal {E}}}_\mathrm{f}({\mathcal {E}}_\mathrm{i}) , \mathrm{0.03\, dex} \right)$| is a base-10 lognormal of width 0.03 dex modelling the scatter around the average relation of equation (12).

An implementation of this procedure is made available online,¹¹ taking as input the (stellar) initial energy distribution (parametrized through equation 13) as well as a desired dark matter remnant mass fraction M_mx/M_mx0, and returning the evolved stellar energy distribution.