# Disruption of Dark Matter Substructure: Fact or Fiction?

Frank C. van den Bosch, Go Ogiya, Oliver Hahn, Andreas Burkert

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

## To cite this version:

Frank C. van den Bosch, Go Ogiya, Oliver Hahn, Andreas Burkert. Disruption of Dark Matter Substructure: Fact or Fiction?. Monthly Notices of the Royal Astronomical Society, 2018, 474 (3), pp.3043-3066. 10.1093/mnras/stx2956. hal-01703867

HAL Id: hal-01703867

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

Submitted on 21 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Disruption of dark matter substructure: fact or fiction?

Frank C. van den Bosch,<fnref n="1" /> Go Ogiya,<fnref n="2" /> Oliver Hahn<fnref n="2" /> and Andreas Burkert<fnref n="3" />,<fnref n="4" />

Accepted 2017 November 13. Received 2017 November 2; in original form 2017 August 17

## ABSTRACT

Accurately predicting the demographics of dark matter (DM) substructure is of paramount importance for many fields of astrophysics, including gravitational lensing, galaxy evolution, halo occupation modelling, and constraining the nature of DM. Because of its strongly non-linear nature, DM substructure is typically modelled using N-body simulations, which reveal that large fractions of DM subhaloes undergo complete disruption. In this paper, we use both analytical estimates and idealized numerical simulations to investigate whether this disruption is mainly physical, due to tidal heating and stripping, or numerical (i.e. artificial). We show that, contrary to naive expectation, subhaloes that experience a tidal shock $\Delta E$ that exceeds the subhalo's binding energy, $|E_b|$, do not undergo disruption, even when $\Delta E/|E_b|$ is as large as 100. Along the same line, and contrary to existing claims in the literature, instantaneously stripping matter from the outskirts of a DM subhalo also does not result in its complete disruption, even when the instantaneous remnant has positive binding energy. In addition, we show that tidal heating due to high-speed (impulsive) encounters with other subhaloes ('harassment') is negligible compared to the tidal effects due to the host halo. Hence, we conclude that, in the absence of baryonic processes, the complete physical disruption of CDM substructure is extremely rare and that most disruption in numerical simulations therefore must be artificial. We discuss various processes that have been associated with numerical overmerging and conclude that inadequate force softening is the most likely culprit.

**Key words:** gravitation – methods: analytical – methods: numerical – galaxies: haloes – dark matter.

## 1 INTRODUCTION

In the $\Lambda$ cold dark matter ($\Lambda$CDM) concordance cosmology, structure develops as primordial density fluctuations grow to form virialized structures. Because of the negligible thermal velocity of CDM, fluctuations survive the early universe on all scales, and structure develops from the bottom up as small, dense CDM clumps merge together to build-up large dark matter (DM) haloes. During this hierarchical assembly, the inner regions of early virialized haloes often survive accretion on to a larger system, thus giving rise to a population of subhaloes. These are subjected to forces that try to dissolve them: dynamical friction, tides from the host halo, and impulsive encounters with other substructure.

Being able to accurately predict the abundance and demographics of DM substructure is of crucial importance for a wide range of

astrophysics. First of all, it is one of the prime discriminators between different DM models; if DM is warm (WDM), rather than cold, free-streaming and enhanced tidal disruption will cause a significantly reduced abundance of low-mass subhaloes (e.g. Knebe et al. 2008; Lovell et al. 2014; Colín et al. 2015; Bose et al. 2016). If DM has a significant cross-section for self-interaction (SIDM), subhaloes are predicted to have constant-density cores, with significantly lower central densities than their CDM counterparts (e.g. Burkert 2000; Vogelsberger, Zavala & Loeb 2012; Rocha et al. 2013). The main methods that are used to test these different DM models are gravitational lensing, either by using time-delays (e.g. Keeton & Moustakas 2009), flux-ratio anomalies (e.g. Metcalf & Madau 2001; Bradač et al. 2002; Dalal & Kochanek 2002), or perturbations of the surface brightness distribution of lensed arcs and Einstein rings (e.g. Vegetti & Koopmans 2009; Vegetti et al. 2014). Alternatively, low-mass substructure can reveal its presence by creating gaps in tidal (stellar) streams (e.g. Carlberg 2009; Erkal et al. 2016; Sanders, Bovy & Erkal 2016). The detailed structure of

<footnote n="1">¹Department of Astronomy, Yale University, PO Box 208101, New Haven, CT 06520, USA</footnote>

<footnote n="2">²Laboratoire Lagrange, Observatoire de la Côte d'Azur, CNRS, Blvd de l'Observatoire, CS 34229, F-06304 Nice cedex 4, France</footnote>

<footnote n="3">³University Observatory Munich (USM), Scheinerstr. 1, D-81679 Munich, Germany</footnote>

<footnote n="4">⁴Associated Scientist, Max-Planck-Institute for Extraterrestrial Physics (MPE), Giessenbachstr. 1, D-85748 Garching, Germany</footnote>

<footnote>* E-mail: frank.vandenbosch@yale.edu</footnote>
subhaloes can be tested using kinematic data of dwarf galaxies (e.g. Mateo et al. 1993; Łokas 2009; Walker et al. 2009). The latter has given rise to a potential problem for the ΛCDM paradigm in the form of the ‘too-big-to-fail’ problem (Boylan-Kolchin et al. 2010; Boylan-Kolchin, Bullock & Kaplinghat 2011; Garrison-Kimmel et al. 2014; Ogiya & Burkert 2015), which has sparked renewed interest in alternatives such as WDM or SIDM, although baryonic solutions have also been proposed (e.g. Zolotov et al. 2012; Ogiya & Mori 2014; Arraki et al. 2014; Dutton et al. 2016). Finally, DM substructure boosts the expected DM annihilation signal (Bergström et al. 1999), an effect that is typically quantified in terms of a ‘boost-factor’ (e.g. Strigari et al. 2007; Giocoli, Pieri & Tormen 2008b; Kuhlen, Diemand & Madau 2008; Pieri, Bertone & Branchini 2008; Moliné et al. 2017).

In addition to carrying a wealth of information regarding the dark sector, DM substructure is also important for our quest to understand galaxy formation and large-scale structure. DM subhaloes are believed to host satellite galaxies and the demographics of DM substructure is therefore directly related to the (small scale) clustering of galaxies. This idea underlies the popular technique of subhalo abundance matching (e.g. Vale & Ostriker 2004; Conroy, Wechsler & Kravtsov 2006; Guo et al. 2010; Behroozi, Wechsler & Conroy 2013b; Hearin et al. 2013; Moster, Naab & White 2013), in which galaxies are assigned to DM host and sub-haloes in simulations to create mock galaxy samples. This has become a prime tool for interpreting galaxy clustering, galaxy–galaxy lensing, and group multiplicity functions, and is even used to constrain cosmological parameters (Trujillo-Gomez et al. 2011; Reddick et al. 2014; Zentner, Hearin & van den Bosch 2014; Hearin, Watson & van den Bosch 2015; Hearin et al. 2016; Zentner et al. 2016; Lehmann et al. 2017).

The formation and evolution of DM substructure is best studied using numerical *N*-body simulations. Nowadays, large cosmological simulations routinely resolve an entire hierarchy of substructure, with haloes hosting subhaloes, which themselves host sub-subhaloes, etc. (e.g. Ghigna et al. 1998; Tormen, Diaferio & Syer 1998; Diemand et al. 2004a; Gao et al. 2004; Kravtsov et al. 2004; Giocoli, Tormen & van den Bosch 2008a; Giocoli et al. 2010). Since simulations are CPU-expensive and suffer from limiting mass and force resolution, various authors have attempted to construct semi-analytical models for the evolution of DM substructure (e.g. Taylor & Babul 2001; Peñarrubia & Benson 2005; van den Bosch, Tormen & Giocoli 2005; Zentner et al. 2005; Kampakoglou & Benson 2007; Pullen, Benson & Moustakas 2014; Jiang & van den Bosch 2016). However, since we lack accurate, analytical descriptions, based on first principles, for how the mass, structure and orbits of subhaloes evolve with time, these models typically have several free parameters that are tuned to reproduce results from numerical simulations.

Prior to 1997, numerical *N*-body simulations suffered from a serious ‘overmerging’ problem, in that the DM haloes revealed little to no substructure. This was in clear contrast with the wealth of ‘substructure’ (i.e. satellite galaxies) observed in galaxy groups and clusters. While some speculated that baryonic physics would resolve this problem (e.g. Frenk et al. 1988), others argued that overmerging might actually be a numerical artefact arising from insufficient mass and/or force resolution. In particular, Carlberg (1994) and van Kampen (1995, 2000a) identified particle-subhalo two-body heating as the main cause for overmerging, while others blamed inadequate force softening (e.g. Moore, Katz & Lake 1996b; Klypin et al. 1999a). At the close of the last millennium, when simulations started to resolve surviving populations of subhaloes

(e.g. Tormen, Bouchet & White 1997; Brainerd, Goldberg & Verner Villumsen 1998; Moore et al. 1998), the discussion as to what might cause numerical overmerging was quickly eclipsed by a new issue, namely that CDM simulations seem to predict too much substructure (Klypin et al. 1999b; Moore et al. 1999). Over the years, though, it has become clear that this ‘missing satellite problem’ is mainly a manifestation of poorly understood baryonic physics related to galaxy formation (see Bullock & Boylan-Kolchin 2017, for a review).

As for overmerging, it has become clear that modern, state-of-the-art numerical simulations still suffer from numerical overmerging, with important ramifications for small-scale clustering (e.g. Conroy, Wechsler & Kravtsov 2006; Guo & White 2014; Campbell et al. 2017; Moster, Naab & White 2017), semi-analytical models for galaxy formation (e.g. Springel et al. 2001; Kang et al. 2005; Kitzbichler & White 2008), and other studies that rely on the phase-space distribution of subhaloes (e.g. Faltenbacher & Diemand 2006; Wu et al. 2013; Tollet et al. 2017). However, the general consensus seems to be that this numerical overmerging only affects subhaloes below a mass resolution limit of 50–100 particles. This notion is based on the fact that simulations seem to yield consistent, converged results for the mass function of subhaloes above this mass resolution limit (e.g. Springel et al. 2008; Onions et al. 2012; Knebe et al. 2013; Cautun et al. 2014; Griffen et al. 2016; van den Bosch & Jiang 2016).

In a recent study, van den Bosch (2017) showed that subhalo disruption is extremely prevalent in modern simulations, with inferred fractional disruption rates (at $z=0$) of ~13 per cent per Gyr (see also Diemand, Moore & Stadel 2004b). This implies that ~65 (90) per cent of all subhaloes accreted around $z=1$ ($z=2$) are disrupted by $z=0$ (Han et al. 2016; Jiang & van den Bosch 2017). In the simulation studied by van den Bosch (2017), roughly 20 per cent of this disruption occurs above the mass resolution limit of 50 particles, with a mass function of disrupting subhaloes that is indistinguishable from that of the surviving population! What is the dominant cause of this prevalent disruption of subhaloes in numerical simulations? In particular, is it artificial (numerical) or real (physical)? Based on the fact that subhalo mass functions are deemed converged above a mass resolution limit of 50–100 particles, it is tempting to conclude that disruption above this ‘resolution limit’ be physical, rather than numerical. However, convergence is a necessary, but not a sufficient condition, to guarantee that the simulation results are physically correct. Furthermore, there is no consensus as to what physical mechanism dominates, with most studies arguing either for tidal heating or tidal stripping. We therefore initiated a comprehensive study aimed at answering these questions. While researching the existing literature on this topic, we came across a variety of conflicting results and claims, and this paper is largely intended to clarify this confusion, correct a number of misconceptions, and present a comprehensive overview. As such, this paper serves as a compendium to two follow-up papers. In van den Bosch & Ogiya (2018, hereafter Paper II), we use a large suite of idealized numerical simulations and experiments to show that state-of-the-art numerical *N*-body simulations still suffer from an enormous amount of overmerging, mainly driven by inadequate force softening and an amplification of discreteness noise due to the tidal field. In Ogiya et al. (in preparation; hereafter Paper III), we use a large suite of similar idealized simulations, which are properly converged, to probe large parts of parameter space. The results of these simulations can be used to calibrate semi-analytical models of subhalo evolution and to gauge the reliability and accuracy of cosmological simulations. In this paper, we
(i) show that impulsive, tidal heating during peri-centric passage of the host halo can inject as much as 100 times the binding energy. Yet, even under those extreme conditions, a remnant containing ~20 per cent of the mass remains.

(ii) demonstrate that harassment has a negligible impact on the evolution of DM substructure.

(iii) underscore the problems and uncertainties associated with analytical treatments of tidal stripping.

(iv) demonstrate that simply removing matter from the outskirts of a DM subhalo will never result in its complete disruption, even if the remnant has positive binding energy.

(v) correct and elucidate a number of errors and misconceptions regarding the numerical, artificial disruption of substructure.

We start in Section 2 with a discussion of the concept of 'tidal radius', including an overview of the various definitions that are used in the literature. Section 3 describes the idealized, numerical simulations that we use in this paper for comparison with simple analytical estimates of various physical disruption mechanisms (Section 4), and for assessing the impact of instantaneously stripping matter off of DM haloes (Section 5). Section 6 discusses a variety of numerical processes that might potentially give rise to artificial disruption of substructure in N-body simulations, including two-body relaxation, impulsive heating due to interactions with overly massive host halo particles, and inadequate force softening. Finally, Section 7 summarizes our findings and briefly discusses some implications for various areas of astrophysics.

## 2 TIDAL RADIUS

We start our discussion of the tidal evolution of DM substructure by addressing an important concept, namely the tidal radius. There are a number of different definitions that are used in the literature, with varying degrees of approximation. For the sake of completeness, we briefly overview and compare the most commonly used definitions.

Consider two point masses, *m* and *M*, separated by a distance *R*. If we define the tidal radius, *r*<sub>t</sub>, of *m* as the radius from *m* at which the tidal force from *M* exceeds the self-gravity of *m*, then we have that

$$r_t = R \left( \frac{m}{2M} \right)^{1/3} . \qquad (1)$$

This is known as the Roche limit. In general, however, *m* and *M* will be in motion with respect to each other, resulting in an additional centrifugal force. If we assume that the two bodies are on a circular orbit of radius *R*, we can account for this using a derivation based on the restricted three-body problem, resulting in the Jacobi limit

$$r_t = R \left( \frac{m}{3M} \right)^{1/3} . \qquad (2)$$

(e.g. Spitzer 1987; Binney & Tremaine 2008).

Going beyond point particles, and taking account of the extended mass profiles, $m(r)$ and $M(r)$, of the two bodies, the tidal force due to *M* is equal to the gravitational attraction due to *m* at a tidal radius

$$r_{t,1} = R \left[ \frac{m(r_{t,1})/M(R)}{2 - \frac{d \ln M}{d \ln R} |_R} \right]^{1/3} . \qquad (3)$$

(e.g. Tormen et al. 1998). Note that $r_{t,1}$ is only defined when $d \ln M/d \ln R \le 2$. This is a consequence of the fact that (when ignoring angular momentum) the tidal field becomes compressive in all directions (i.e. no mass will be stripped) whenever $d \ln \rho/d \ln R \le -1$ (e.g. Dekel, Devor & Hetzroni 2003).

Although equation (3) is used in several studies (e.g. Klypin et al. 1999a; Hayashi et al. 2003; Taffoni et al. 2003), it ignores the centrifugal force. Assuming that *m* is on a circular orbit of radius *R* around the centre of *M*, this correction results in

$$r_{t,2} = R \left[ \frac{m(r_{t,2})/M(R)}{2 + \frac{\Omega^2 R^3}{GM(R)} - \frac{d \ln M}{d \ln R} |_R} \right]^{1/3} = \left[ \frac{Gm(r_{t,2})}{\Omega^2 - \frac{d^2 \Phi_h}{d \ln^2 R}} \right]^{1/3} . \qquad (4)$$

where $\Omega$ is the angular speed (e.g. King 1962; Tollet et al. 2017). Using that along a circular orbit $\Omega = V_{\text{circ}}(R)/R$, we see that equation (4) reduces to

$$r_{t,2} = R \left[ \frac{m(r_{t,2})/M(R)}{3 - \frac{d \ln M}{d \ln R} |_R} \right]^{1/3} . \qquad (5)$$

Hence, for subhaloes on a circular orbit, the tidal field only becomes compressive in all directions when $d \ln M/d \ln R \ge 3$ ($d \ln \rho/d \ln R \ge 0$). This implies that subhaloes on a circular orbit will always be subject to tidal stripping, unless they happen to find themselves orbiting within the central region of a host halo with a constant density core (see Dekel et al. 2003, for a comprehensive discussion). We emphasize that the above definition for the tidal radius is strictly only valid along circular orbits. However, it is common practice to also apply it to eccentric orbits, in which case one estimates the instantaneous tidal radius by simply using $R$ and $\Omega$ the corresponding instantaneous values, i.e. $\Omega = |\vec{V} \times \vec{R}|/R^2$, where $\vec{V}$ is the instantaneous orbital velocity of the subhalo (see e.g. King 1962; Taylor & Babul 2001; Zentner & Bullock 2003; Peñarrubia & Benson 2005).

Finally, for the sake of completeness, we mention that Klypin et al. (1999a, 2015) have advocated the use of yet another definition of tidal radius. Motivated by the work of Weinberg (1994a,b, 1997), which demonstrates that resonant effects can boost the impact of tidal stripping, Klypin et al. define their tidal radius as the radius where the frequency of the tidal force by the host halo is equal to that of the internal motion within the subhalo. This results in

$$r_{t,3} = R \left[ \frac{m(r_{t,3})}{M(R)} \right]^{1/3} . \qquad (6)$$

Following the suggestion by Klypin et al. (1999a), several studies have taken the tidal radius to be the minimum of $r_{t,1}$ and $r_{t,3}$ (e.g. Oguri & Lee 2004; Arraki et al. 2014; Klypin et al. 2015).

Although the tidal radius is a common concept used to discuss tidal stripping, it is important to realize that all the various definitions discussed above are only approximate (see e.g. Binney & Tremaine 2008; Mo, van den Bosch & White 2010). First of all, the two-dimensional surface along which the radial acceleration of a test particle vanishes is not spherical, even if *m* and *M* are point masses, and can therefore not be characterized by a single radius (see Tollet et al. 2017, for a detailed discussion). Secondly, variance in the orbital motion of the particles within the subhalo gives rise to scatter in the coriolis forces on these particles, which in turn introduces some 'thickness' to the shell of particles for which $\ddot{r} = 0$ (see Read et al. 2006, for a detailed discussion). Thirdly, the various expressions for $r_t$ are all derived under the distant-tide approximation, which implies that the distance *R* between *m* and *M* is large compared to the size of *m*. This condition is not necessarily valid in the case of subhaloes, since the orbit's peri-centre can be comparable to, or even smaller than, the size of the subhalo. Finally, the optimal definition of tidal radius is likely to depend on the orbit in question. For instance, whereas $r_{t,2}$ (equation 4) is expected to be most applicable to circular orbits, $r_{t,1}$ (equation 3) is likely to be more relevant for more radial orbits. In what follows, we will
consider both in order to assess how the ‘choice’ for the definition of tidal radius impacts simple models for tidal stripping. As for $r_{t,3}$ (equation 6), since this definition of tidal radius is basically bracketed by the other two, i.e. one typically has $r_{t,2} < r_{t,3} < r_{t,1}$, we will no longer consider this definition of the tidal radius in what follows.

## 2.1 Properties of newly accreted subhaloes

As a subhalo moves towards the central region of its host halo, its tidal radius typically shrinks,<fnref n="1" /> reaching a minimum value during peri-centric passage. Hence, it is commonly assumed that the amount of material stripped off a subhalo during a radial orbital period depends on its tidal radius at peri-centre. In what follows, we refer to this as the ‘peri-tidal radius’, for which we use the symbol $r_{t|p}$.

What are realistic values for $r_{t|p}$ for DM subhaloes during their first radial orbit? In order to address this question, we consider DM subhaloes in the cosmological Bolshoi simulation (Klypin, Trujillo-Gomez & Primack 2011), which follows the evolution of 2048³ DM particles using the Adaptive Refinement Tree (ART) code (Kravtsov, Klypin & Khokhlov 1997) in a flat $\Lambda$CDM model with parameters $\Omega_{m,0} = 1 - \Omega_{\Lambda,0} = 0.27$, $\Omega_{b,0} = 0.0469$, $h = H_0/(100 \text{ km s}^{-1} \text{ Mpc}^{-1}) = 0.7$, $\sigma_8 = 0.82$ and $n_s = 0.95$ (hereafter ‘Bolshoi cosmology’). The box size of the Bolshoi simulation is $L_{\text{box}} = 250h^{-1} \text{ Mpc}$, resulting in a particle mass of $m_p = 1.35 \times 10^8 h^{-1} \text{ M}_\odot$.

We use the publicly available halo catalogues<fnref n="2" /> obtained using the phase-space halo finder ROCKSTAR (Behroozi, Wechsler & Wu 2013a), which uses adaptive, hierarchical refinement of friends-of-friends groups in six phase-space dimensions and one time dimension. ROCKSTAR haloes are defined as spheres with an average density equal to $\Delta_{vir}(z)\rho_{crit}(z)$. Here, $\rho_{crit}(z) = 3H^2(z)/8\pi G$ is the critical density for closure at redshift $z$, and $\Delta_{vir}(z)$ is given by the fitting function of Bryan & Norman (1998).

Following van den Bosch (2017, hereafter ‘vdB17’), we use the 19 simulation outputs at $z \le 0.0605$, from which we select a random subset of 5000 newly accreted subhaloes (labelled ‘A’ subhaloes in vdB17) by identifying those subhaloes in the halo catalogues that in the previous output were still considered host haloes and which had never before been subhaloes. In order to focus on well-resolved systems, each subhalo is requested to have at least 250 particles and to be accreted in a host halo with mass $M_h \ge 3 \times 10^{13} h^{-1} \text{ M}_\odot$. For each of the subhaloes thus selected, we register

(i) $x_c \equiv R_c(E)/R_h$, the radius of the circular orbit corresponding to the orbital energy, E, expressed in terms of the virial radius of the host halo, $R_h$.

(ii) $\eta \equiv L/L_c(E)$, the orbital circularity, defined as the ratio of the orbital angular momentum, L, and the angular momentum $L_c(E)$ corresponding to a circular orbit of energy E. Radial and circular orbits have $\eta = 0$ and 1, respectively.

(iii) $M_s/M_h$, the ratio of subhalo mass to host halo mass.

(iv) $c_s$ the concentration parameter of the subhalo.

(v) $c_h$ the concentration parameter of the host halo.

The concentration parameters are defined as the ratios of the virial radius and the scale parameter of the NFW profile (see below) that best fits the density distribution of the (sub)halo (see Behroozi

et al. 2013a, for details on how these quantities are computed in ROCKSTAR).

For each subhalo thus selected we then compute the orbit’s apo-centre, $R_a$, and peri-centre, $R_p$, which are given by the two roots for

$$ \frac{1}{R^2} + \frac{2[\Phi_h(R) - E]}{L^2} = 0, \qquad (7) $$

as well as the radial, orbital period, which is given by

$$ T_r = 2 \int_{R_a}^{R_p} \frac{dR}{\sqrt{2[E - \Phi_h(R)] - L^2/R^2}}, \qquad (8) $$

and the increase in azimuthal angle during a radial period

$$ \Delta\Psi = 2L \int_{R_a}^{R_p} \frac{dR}{R^2 \sqrt{2[E - \Phi_h(R)] - L^2/R^2}} \qquad (9) $$

(Binney & Tremaine 2008). Throughout, we assume that the host halo is an NFW halo (Navarro, Frenk & White 1997), for which

$$ \Phi_h(R) = -V_h^2 \frac{\ln(1 + c_h x)}{f(c_h) x}, \qquad (10) $$

where $V_h = \sqrt{GM_h/R_h}$ is the host halo’s virial velocity, $x = R/R_h$, and $f(x) = \ln(1+x) - x/(1+x)$. Finally, we compute the subhalo’s peri-tidal radius, $r_{t|p}$, in units of the subhalo’s virial radius, $r_s$, using both definitions of the tidal radius; $r_{t,1}$ and $r_{t,2}$. Note that all these quantities are computed using the instantaneous (i.e. at the epoch of accretion) properties of the host halo and subhalo. We thus ignore the fact that $\Phi_h(R)$ is time dependent and may change appreciably between accretion and first peri-centric passage.

The upper left and middle panels of Fig. 1 show the distributions of $x_c$ and $\eta$, characterizing the orbits of newly accreted subhaloes. The median $x_c$ is 1.26, and the distribution is skewed towards higher values, indicating that at accretion the orbits of subhaloes are only weakly bound (cf. vdB17). The $\eta$ distribution is fairly symmetric around $\eta \sim 0.5$, indicating that purely radial and circular orbits are rare (cf. Tormen 1997; Zentner et al. 2005; Wetzel 2011; Jiang et al. 2015). The mass ratio, $M_s/M_h$, is peaked around $10^{-3}$ (upper right-hand panel),<fnref n="3" /> while the concentration ratio follows a roughly lognormal distribution centred around $c_s/c_h \approx 2$. Hence, subhaloes are, on average, more concentrated than their host haloes, which is a natural outcome of the concentration–mass relation of CDM haloes. The lower, middle panel depicts the ratio of the subhalo’s peri-centre $R_p$ to the host halo’s virial ratio, $R_h$. The median $R_p/R_h = 0.37$, while 7.3 per cent (0.9 per cent) of subhaloes have a first peri-centric passage $R_p < 0.1 R_h$ (0.03 $R_h$). These numbers are in good agreement with (Wetzel 2011, see his appendix A). Finally, the lower-right-hand panel plots the distributions of the peri-tidal radius, $r_{t|p}$, in units of the subhalo’s original (at accretion) virial radius, $r_s$. Red and blue histograms correspond to the two different definitions of the tidal radius, $r_{t,1}$ and $r_{t,2}$, respectively. Note that $r_{t,1}$, which ignores the centrifugal force, is significantly larger than $r_{t,2}$; the median $r_{t,1}$ is 0.43$r_s$, with only 0.7 per cent of subhaloes having $r_{t,1}/r_s < 0.1$. In comparison, the median $r_{t,2}$ is 0.21$r_s$, with 21.5 per cent (2.0 per cent) of subhaloes having $r_{t,2}/r_s < 0.1$ (<0.01). Thus, we see that the definition of tidal radius can have a big impact on the expected stripping rate (see Section 4.3).

<footnote n="1">¹ Since the tidal radius depends on the density profile of the host halo, this is not always the case though.</footnote>

<footnote n="2">² http://hipacc.ucsc.edu/Bolshoi/MergerTrees.html</footnote>

<footnote n="3">³ The drop-off at $M_s/M_h < 10^{-3}$ is an artefact of our selection criteria and the limiting mass resolution of the simulation (see van den Bosch et al. 2016).</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1. Properties of subhaloes in the Bolshoi simulation at the moment they are accreted by their host halo. From left to right and top to bottom, the different panels show distributions of the orbital energy, characterized by $x_c = R_c(E)/R_h$, orbital circularity, $\eta = L/L_c(E)$, the mass ratio of subhalo to host halo, $M_s/M_h$, the ratio of concentration parameters of subhalo and host halo, $c_s/c_h$, the orbit’s peri-centre in units of the virial radius of the host halo, $R_p/R_h$, and finally the peri-tidal radius, $r_{t|p}$ (defined as the tidal radius at peri-centre), normalized by the size of the subhalo. For the latter, we have used the two different definitions of the tidal radius, as indicated. See the text for a detailed discussion.</figcaption></figure>

## 3 NUMERICAL SIMULATIONS

Although this paper mainly focuses on simple (semi)-analytical treatments of subhalo evolution, we also present some simulation results for comparison. This section briefly describes these numerical simulations, and we refer the interested reader to Papers II and III for more details (and many more simulation results). Our (idealized) simulations follow the evolution of a single N-body subhalo orbiting in the fixed, analytical potential of a host halo. Our goal is to investigate how much mass is stripped from the subhalo during its first radial orbital period and to compare the results to expectations based on the various processes discussed in Section 4.

Each halo (host and sub) is assumed to (initially) be spherical and to have a NFW density profile

$$ \rho(r) = \rho_0 \left(\frac{r}{r_0}\right)^{-1} \left(1 + \frac{r}{r_0}\right)^{-2}, \qquad (11) $$

where $r_0$ is the characteristic scale radius. If we define the virial radius as enclosing an average density that is 97 times the critical density (which is valid for haloes at $z=0$), the crossing time for such a halo is

$$ t_{\text{cross}} = \frac{r_{\text{vir}}}{V_{\text{vir}}} = 2.006 \text{ Gyr}, \qquad (12) $$

where we have adopted a Hubble constant of $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

Throughout we assume that, prior to stripping, the halo has an isotropic velocity distribution, such that its distribution function (DF) depends only on energy, i.e. $f=f(E)$. We use the method of Widrow (2000) to sample particles from the DF using the standard acceptance-rejection technique (Press et al. 1992; Kuijken & Dubinski 1994; Drakos, Taylor & Benson 2017) and truncate the halo at its virial radius. Since the DF that we use to generate the initial conditions (ICs) is computed using the Eddington (1916) inversion equation, which assumes that the halo extends to infinity, our initial system is not going to be in perfect equilibrium. However, as we demonstrate and discuss in Paper II, this has no significant impact on our results.

Throughout we adopt model units in which the gravitational constant, G, the scale radius, $r_0$, and the initial (virial) mass of the subhalo, $M_s$, are all unity. We restrict ourselves to subhaloes with a concentration parameter $c_s = 10$, for which $t_{\text{cross}} = c_s^{3/2} = 31.6$. Based on equation (12), we thus have that a time interval of $\Delta t = 1$ (model units) corresponds to 63.4 Myr.

The simulations are carried out using a treecode designed for graphic processing unit clusters. This code, which is optimized for speed, uses a second-order Runge-Kutta integrator and is described in Ogiya et al. (2013). We have compared the results of this simulation code with those of the treecode described in Section 5.1, with uses a second-order leap-frog integrator, and find the results to be in excellent agreement.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Results for a subset of the idealized numerical simulations described in Section 3 and listed in Table 1. From left to right, the panels show the time-evolution (where $t/T_r$ is the time normalized by the radial period) of the logarithm of the bound mass fraction, $\log[f_{\text{bound}}]$, the orbital radius in units of the virial radius of the host, $r/r_{\text{vir},h}$, and the orbital speed, in units of the virial velocity of the host, $v/V_{\text{vir},h}$. Different colours correspond to different orbital circularities, as indicated. Note how the subhalo loses between ~40 ($\eta = 1$) and ~90 ($\eta = 0$) per cent of its mass during its first radial period.</figcaption></figure>

We integrate the subhalo in the external tidal field of an NFW host halo whose virial mass (radius) is 1000 (10) times larger than that of the subhalo. For such a large mass ratio, dynamical friction, which is not accounted for in our idealized simulations, can be safely ignored. In the $\Lambda$CDM cosmology, to good approximation, concentration scales with halo mass as $c \propto M^{-0.1}$ (e.g. Dutton & Macciò 2014). Hence, for a mass ratio of 1000, the ratio in concentration parameters is roughly 2, in agreement with the results shown in Fig. 1, and we therefore adopt a concentration for the host halo of $c_h = 5$. Each subhalo is modelled using $N = 10^7$ particles and a softening length of $\epsilon = 0.003$. All simulations adopt a tree opening angle of $\theta = 0.7$, and a fixed time step of $\Delta t = 0.02$ (corresponding to ~1.3 Myr). The latter implies that we resolve the minimal orbital time of the system, defined as $\tau_{\min} \simeq [3\pi/G\bar{\rho}(<\epsilon)]^{1/2}$, with $\bar{\rho}(<\epsilon)$ the average density enclosed by the softening length, with 30 time steps. As discussed in great detail in Paper II, for these parameters we obtain results that are both converged (i.e. stable to an increase in $N$) and reproducible (i.e. different random realizations yield indistinguishable results). We initially position the subhalo at the apo-centre of an orbit with $x_c = 1$ and $\eta = 0.0, 0.1, ..., 1.0$ and integrate the system for a total of 10 Gyr (corresponding to ~7890 time steps). This is slightly longer than a radial orbital period, which varies from $T_r = 9.18$ Gyr for $\eta = 0.0$ to $T_r = 9.53$ Gyr for $\eta = 1.0$. We use an extremely robust, iterative approach to determine the bound fraction, $f_{\text{bound}}$, as function of time, which is described in detail in Appendix A (see also Paper II).

## 3.1 Results

Fig. 2 shows some of the simulation results. The left-hand panel shows the logarithm of the bound fraction as a function of time, normalized by the radial period, $T_r$, while the middle and right-hand panels show the radius and velocity of the subhalo with respect to the host halo, again as functions of $t/T_r$. Different colours correspond to different values of the orbital circularity, $\eta$, as indicated, and for clarity, we only show the results for a subset of six simulations. Table 1 lists the bound fractions at five different epochs for the full set of 11 simulations. Note how, for fixed orbital energy, subhaloes on more radial orbits are stripped of a larger fraction of their mass per orbital period. For the purely radial orbit ($\eta = 0$), as much as 89 per cent of the original subhalo mass is stripped away during its first radial orbit.

<table><caption>**Table 1.** Simulation results.</caption><thead><tr><th>ID</th><th>η</th><th>f<sub>bound</sub><br>0.2T<sub>r</sub></th><th>f<sub>bound</sub><br>0.4T<sub>r</sub></th><th>f<sub>bound</sub><br>0.6T<sub>r</sub></th><th>f<sub>bound</sub><br>0.8T<sub>r</sub></th><th>f<sub>bound</sub><br>T<sub>r</sub></th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th></tr></thead><tbody><tr><td>S0</td><td>0.0</td><td>0.982</td><td>0.909</td><td>0.120</td><td>0.111</td><td>0.110</td></tr><tr><td>S1</td><td>0.1</td><td>0.982</td><td>0.908</td><td>0.262</td><td>0.244</td><td>0.230</td></tr><tr><td>S2</td><td>0.2</td><td>0.982</td><td>0.905</td><td>0.367</td><td>0.347</td><td>0.323</td></tr><tr><td>S3</td><td>0.3</td><td>0.981</td><td>0.901</td><td>0.429</td><td>0.409</td><td>0.378</td></tr><tr><td>S4</td><td>0.4</td><td>0.981</td><td>0.894</td><td>0.479</td><td>0.454</td><td>0.419</td></tr><tr><td>S5</td><td>0.5</td><td>0.980</td><td>0.887</td><td>0.522</td><td>0.489</td><td>0.453</td></tr><tr><td>S6</td><td>0.6</td><td>0.979</td><td>0.877</td><td>0.559</td><td>0.520</td><td>0.484</td></tr><tr><td>S7</td><td>0.7</td><td>0.977</td><td>0.867</td><td>0.594</td><td>0.550</td><td>0.514</td></tr><tr><td>S8</td><td>0.8</td><td>0.975</td><td>0.858</td><td>0.631</td><td>0.580</td><td>0.543</td></tr><tr><td>S9</td><td>0.9</td><td>0.969</td><td>0.849</td><td>0.672</td><td>0.612</td><td>0.573</td></tr><tr><td>S10</td><td>1.0</td><td>0.942</td><td>0.837</td><td>0.732</td><td>0.664</td><td>0.613</td></tr></tbody></table>

*Note.* The bound fractions at five different epochs, $t/T_r = 0.2, 0.4, ..., 1.0$ (Columns 3–7), inferred from our idealized simulations of subhaloes orbiting in the fixed potential of an NFW host. Each simulation has $M_{h}/M_{s} = 1000$, $c_{h} = 5$, $c_{s} = 10$ and $x_{c} = 1.0$, and only differ in their value for the orbital circularity, $\eta$ (Column 2). In each simulation, the subhalo, which is simulated with $N = 10^7$ particles, starts out at its own apo-centre (cf. Fig. 2).

During peri-centric passage, the bound mass fraction can fluctuate rapidly, especially for the more radial orbits (see also Diemand, Kuhlen & Madau 2007; Han et al. 2012). This is a manifestation of the re-binding of some particles as a consequence of the subhalo undergoing re-virialization. The tidal shock causes a stretching of the subhalo in the direction towards the centre of the host halo and a compression in the direction perpendicular to the orbital plane. Hence, immediately following the tidal shock the subhalo finds itself out of virial equilibrium. While some matter has received a large enough boost in kinetic energy that it leaves the subhalo indefinitely, some of the particles that are unbound immediately following the tidal shock get rebound as the subhalo's potential re-adjusts as part of the re-virialization process.

What is most relevant for the discussion in this paper is that none of the simulated subhaloes undergo complete disruption. All orbits simulated have $x_c = 1$, which is at the low end of the distribution (cf. the upper-left-hand panel of Fig. 1). Hence, these are among the most bound orbits among newly accreted subhaloes. And the simulations suggest that at most about 90 per cent of the subhalo mass will be stripped during the first radial period (which lasts
between 9 and 9.5 Gyr). It is intriguing to contrast this with results from the cosmological Bolshoi simulation, in which a significant fraction of subhaloes are disrupted during, or shortly after, their first peri-centric passage (vdB17). This suggests that much of this disruption is likely to be a numerical artefact. We will elaborate on this in the following sections, as well as in Papers II and III.

## 4 PHYSICAL DISRUPTION PROCESSES

As discussed in Section 1, DM subhaloes are subjected to tidal forces that strip the subhaloes of their mass and which may ultimately result in their complete disruption. Depending on the rate with which the tidal forces change, the subhalo will respond differently. In the limit where the external tidal field changes slowly (i.e. along close-to-circular orbits), the subhalo will lose mass from beyond some limiting tidal radius. We refer to this as ‘tidal stripping’. If the tidal field changes rapidly, particles in the subhalo experience an impulsive kick, resulting in a net heating of the subhalo. We refer to this as tidal heating or tidal shocking. Throughout we distinguish between tidal heating due to the host halo (associated with the fast peri-centric passage along highly eccentric orbits), and due to high-speed encounters with other subhaloes, to which we refer as ‘harassment’. In this section, we describe each of these processes, in turn, and use simple analytical estimates together with the simulation results presented in the previous section to compare their relative importance.

Another process that is relevant for the evolution of DM substructure is dynamical friction, which transfers orbital energy and angular momentum of the subhalo to the particles of the host halo, thereby causing the subhalo to ‘sink’ towards the centre of the host, where the tidal forces are stronger. Dynamical friction is only expected to be significant for the most massive subhaloes, with a mass that exceeds a few per cent of that of the host halo (e.g. Binney & Tremaine 2008; Mo et al. 2010), and various studies have pointed out that indeed dynamical friction makes a negligible contribution to the evolution of the ensemble of subhaloes (e.g. Zhao 2004; Peñarrubia & Benson 2005; Ogiya & Burkert 2016). In this paper, we mainly focus on less massive subhaloes, which vastly outnumber their more massive counter-parts, and therefore ignore dynamical friction in what follows.

### 4.1 Impulsive heating due to the host halo

During the high-speed peri-centric passage, the rapid, impulsive change in the external potential causes a transfer of orbital energy to subhalo internal energy. The resulting increase in subhalo internal energy, $\Delta E$, can be computed using the impulse approximation (Spitzer 1958), and in what follows we closely follow the treatment of Gnedin, Hernquist & Ostriker (1999, hereafter GH099).

Consider a subhalo orbiting a spherical, NFW host halo of mass $M_h$ and concentration parameter $c_h$. Let $\vec{R}_{\text{orb}}(t)$ be the trajectory of the centre of the subhalo with respect to that of the host halo, and $r$ be the distance of a subhalo particle from the centre of the subhalo. The position vector of that particle is then indicated by $\vec{R} = \vec{R}_{\text{orb}} + \vec{r}$. The subhalo’s orbit is confined to a plane and, in the absence of dynamical friction, characterized by two independent parameters; the orbital energy $E$ and the orbital angular momentum $L$, or equivalently, the parameters $x_c = R_c(E)/R_h$, and $\eta = L/L_c(E)$ (see Section 2.1). Following GH099, we assume that the orbit is in the X-Y plane, so that $\vec{R}_{\text{orb}}(t) = \cos\theta(t)\hat{X} + \sin\theta(t)\hat{Y}$, where $\theta = 0$ for $R_{\text{orb}} = R_p$ (i.e. at peri-centric passage). Then, the impulsive

change in the velocity of a subhalo particle located at $\vec{r} = (x, y, z)$ can be written as

$$ \Delta \vec{v} = V_h \frac{Q}{R_h} \{(B_1 - B_3)x, (B_2 - B_3)y, -B_3 z\}, \quad (13) $$

where

$$ Q = \frac{1}{\eta x_c^{1/2}} \left[ \frac{f(c_h)}{f(x_c c_h)} \right]^{1/2}, \quad (14) $$

and<fnref n="4" />

$$ B_1 = \int_{-\theta_m}^{\theta_m} \frac{3\mu_0(\zeta) - \mu_1(\zeta)}{\zeta} \cos^2 \theta \, d\theta, \quad (15) $$

$$ B_2 = \int_{-\theta_m}^{\theta_m} \frac{3\mu_0(\zeta) - \mu_1(\zeta)}{\zeta} \sin^2 \theta \, d\theta, \quad (16) $$

$$ B_3 = \int_{-\theta_m}^{\theta_m} \frac{\mu_0(\zeta)}{\zeta} d\theta. \quad (17) $$

Here, $\zeta = \zeta(\theta) \equiv R_{\text{orb}}/R_h$ is the orbital radius of the subhalo in units of the host halo’s virial radius, $\theta_m = \Delta\Psi/2$ (see equation 9), and the functions $\mu_0(x)$ and $\mu_1(x)$ for an NFW host halo are given by

$$ \mu_0(x) = \frac{1}{f(c_h)} \left[ \ln(1 + c_h x) - \frac{c_h x}{1 + c_h x} \right], \quad (18) $$

and

$$ \mu_1(x) = \frac{1}{f(c_h)} \left[ \frac{(c_h x)^2}{(1 + c_h x)^2} \right]. \quad (19) $$

Note that the above integrations are along one orbit $R_{\text{orb}}(\theta)$ of the subhalo, running from apo-centre to apo-centre, as parametrized by the position angle $\theta$. We compute $R_{\text{orb}}(\theta)$ by integrating the subhalo’s orbit in the host halo potential using a fifth-order Cash–Karp Runge–Kutta method with adaptive step-size control (e.g. Press et al. 1992). When computing the integrals for $B_1$, $B_2$ and $B_3$, we use bi-cubic spline interpolation of the resulting $R_{\text{orb}}(\theta)$.

Equation (13) implies that subhalo particles at subhalo-centric radius $r$ experience an average change in specific (kinetic) energy, per radial period, equal to

$$ \langle \Delta E \rangle (r) = \frac{1}{2} \langle \Delta \vec{v} \cdot \Delta \vec{v} \rangle (r) = \frac{1}{6} Q^2 V_s^2 \chi_{\text{orb}} \frac{r^2}{r_s^2}. \quad (20) $$

Here, $\chi_{\text{orb}} \equiv (B_1 - B_3)^2 + (B_2 - B_3)^2 + B_3^2$, and we have used that, by virtue of the definition of ‘virial radius’, $V_h/R_h = V_s/r_s$. In addition, we have assumed an isotropic velocity distribution, so that $\langle \vec{v} \cdot \Delta \vec{v} \rangle = 0$.

The above derivation is based on the impulse approximation, wherein the subhalo particles are assumed to move little during the tidal shock. This is, however, only appropriate for the least bound particles in the outer regions of the subhalo. The more bound particles will have dynamical times that can be substantially shorter than the duration of the tidal shock, $\tau_{\text{sh}} \simeq R_p/V_p$, where $V_p = \sqrt{2[E - \Phi_h(R_p)]}$ is the subhalo speed at peri-centre. For such particles, the peri-centric passage of the host halo is adiabatic, rather than impulsive, and the conservation of adiabatic invariants prevents these particles from gaining any energy (Spitzer 1987; Weinberg 1994a,b). Gnedin & Ostriker (1999) show that one can

<footnote n="4">⁴ Note that the integrals for $B_k$ ($k = 1, 2, 3$) defined here differ from those in GH099 by a factor $c_h$.</footnote>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Left-hand panel: the distribution of $\log[\Delta E/|E_b|]$ for the same set of newly accreted subhaloes as in Fig. 1. Here, $\Delta E$ is the energy injected into the subhalo during its first radial period, due to the impulsive shock associated with peri-centric passage (equation 22), while $E_b$ is the subhalo's binding energy (equation 24). Ochre and green histograms indicate the results obtained with and without the adiabatic correction (equation 21), respectively, which has little impact. The grey-shaded region indicates where $\Delta E > |E_b|$, and thus where the energy injected due to the impulsive shock exceeds the subhalo's original binding energy. The fractions of subhaloes that meet this criterion are indicated. Right-hand panel: the corresponding fraction of mass lost by the subhalo as a consequence of the tidal shock, calculated under the instantaneous mass-loss approximation as described in the text.</figcaption></figure>

correct for this ‘adiabatic shielding’ by simply multiplying $\langle \Delta E \rangle(r)$ of equation (20) with an adiabatic correction of the form

$$A_{\text{corr}}(r) = [1 + \omega^2(r) \tau_{\text{sh}}^2]^{-3/2} \quad (21)$$

Here, $\omega(r)$ is the orbital frequency of subhalo particles at radius $r$. Throughout, we adopt $\omega(r) = \sigma_r(r)/r$, where $\sigma_r(r)$ is the radial velocity dispersion at subhalo-centric radius $r$, obtained by solving the Jeans equation for a spherical, isotropic NFW subhalo (e.g. equation 11 in van den Bosch et al. 2004).

Finally, we can use the above to compute the total energy injected into the subhalo per radial orbit using

$$\begin{aligned} \Delta E &= 4\pi \int_0^{r_s} \rho_s(r) \langle \Delta E \rangle(r) A_{\text{corr}}(r) r^2 dr \\ &= \frac{1}{2} M_s V_s^2 \frac{\chi_{\text{orb}} Q^2}{3 c_s^2 f(c_s)} C_{\text{corr}}, \end{aligned} \quad (22)$$

where

$$C_{\text{corr}} = \int_{0}^{c_s} \frac{x^3 A_{\text{corr}}(x r_0)}{(1+x)^2} dx, \quad (23)$$

with $r_0 = r_s/c_s$ the scale radius of the subhalo.

### 4.1.1 Application to Bolshoi subhaloes

In order to assess the potential impact of tidal heating, we compute $\Delta E/|E_b|$ for each of the 5000 newly accreted subhaloes from the Bolshoi simulation (see Section 2.1). Here, $E_b$ is the total binding energy of the initial subhalo (i.e. prior to the tidal shock). Since we consider the initial subhalo to be truncated at its virial radius, we have that

$$E_b = -\frac{1}{2} M_s V_s^2 f_E(c_s), \quad (24)$$

where

$$f_E(c) = \frac{c}{2 f^2(c)} \left[ 1 - \frac{1}{(1+c)^2} - \frac{2 \ln(1+c)}{(1+c)} \right] \quad (25)$$

(Mo, Mao & White 1998). The ochre histogram in the left-hand panel of Fig. 3 shows the distribution of $\Delta E/|E_b|$ values for the

newly accreted subhaloes in the Bolshoi simulation, computed using equations (22)–(25). The median is 1.9, indicating that a typical subhalo, during its first-pericentric passage, experiences a tidal, impulsive shock that injects an amount of energy that is roughly twice its original binding energy. If we were to associate $\Delta E > |E_b|$ with the complete disruption of the subhalo, as is done for example in Gonzalez-Casado, Mamon & Salvador-Sole (1994), then the inference would be that 72 per cent of subhaloes do not survive their first peri-centric passage. For comparison, the green histogram shows the distribution obtained when ignoring the adiabatic correction (i.e. by setting $A_{\text{corr}} = 1$). This boosts the distributions to slightly larger values of $\Delta E/|E_b|$ and increases the ‘disruption fraction’ to 0.78.

There are two important points to be made here. First of all, the overall impact of adiabatic shielding is fairly modest, and one does not make a large error if one simply ignores the adiabatic correction. Secondly, if indeed subhaloes for which $\Delta E > |E_b|$ were to experience complete disruption, the inference would be that the vast majority of subhaloes do not survive first peri-centric passage. However, this interpretation of $\Delta E/|E_b|$ is extremely naive and dramatically overestimates the efficiency of subhalo disruption. The reason is that $\Delta E$ does not specify how that energy is distributed over the particles. The following reduction ad absurdum makes it clear why this distribution is important; if all the energy were given to only a single particle, the system clearly would not disrupt, even if $\Delta E \gg |E_b|$. Rather, only that single particle would escape and leave the remainder of the halo virtually intact. As is evident from equation (20), the average change in kinetic energy due to an impulsive tidal shock scales with $r^2$. Hence, particles at the outskirts of the subhalo, which need less energy to escape, receive far more energy than the more bound particles in the central regions.

In order to properly assess the impact of $\Delta E$ on the system as a whole, one needs to compare, for each individual particle $i$, the particle’s individual $(\Delta E)_i$ to its individual binding energy, $|E_i|$. One can then express the impact of the tidal shock in terms of the fraction of particles with $(\Delta E)_i > |E_i|$. We will refer to this fraction as $f_{\text{strip}}$. Note that this ‘instantaneous mass-loss approximation’ (cf. Aguilar & White 1985) ignores potential re-binding and unbinding of particles during the re-virialization process immediately following the tidal shock. However, these processes are neither calculable
<figure><img src="image_8.png" /><figcaption>**Figure 4.** The logarithm of the ratio $\Delta E/|E_b|$ (upper panels) and the corresponding stripped mass fraction, $f_{\text{strip}}$ (lower panels) as functions of the orbital circularity, $\eta$. All results assume a host halo to sub-halo mass ratio of 1000. The solid curves indicate the results obtained using the Monte Carlo method described in the text. Results are shown for different values of the orbital energy, as parametrized via $x_c$ (left-hand column), the concentration of the subhalo, $c_s$ (middle column), and the concentration of the host halo, $c_h$ (right-hand column); corresponding values are indicated. The black line corresponds to the fiducial model, which has $x_c = 1$, $c_s = 10$, and $c_h = 5$. Note that $\Delta E/|E_b| = 1$ corresponds to $f_{\text{strip}} \sim 0.2$; hence even if the energy injected impulsively exceeds the original subhalo binding energy, only a modest fraction of subhalo mass will be stripped. The dotted curves in the lower panels indicate the predictions from equation (27).</figcaption></figure>

nor do they appear to have a significant impact, as demonstrated in Section 4.1.2.

In order to compute $f_{\text{strip}}$ for individual subhaloes, we proceed as follows. We construct an N-body realization, consisting of $N_p = 50000$ particles, of each individual subhalo, using the method described in Section 3. We do not evolve these realizations with an actual simulation code; rather for each particle $i$, we analytically compute its binding energy, $E_i = v_i^2/2 + \Phi_s(r_i)$, and

$$ (\Delta E)_i = A_{\text{corr}}(r_i) \left[ \vec{v}_i \cdot \Delta \vec{v}_i + \frac{1}{2} (\Delta v_i)^2 \right], \quad (26) $$

where $\Delta \vec{v}$ is given by equation (13), and the factor $A_{\text{corr}}$ (equation 21) accounts for the adiabatic shielding of the most bound particles. Finally, we compute the fraction $f_{\text{strip}}$ of particles for which $(\Delta E)_i > |E_i|$.

Before showing the results for the subhaloes in the Bolshoi simulation, we first illustrate how $f_{\text{strip}}$ relates to $\Delta E/|E_b|$, by computing both values for a variety of orbits and concentration parameters of both host- and sub-halo. Fig. 4 shows $\Delta E/|E_b|$ (upper panels) and $f_{\text{strip}}$ (lower panels) as functions of the orbital circularity, $\eta$. The left-hand panels show the dependence on the orbital energy, as parametrized by $x_c = R_c(E)/R_h$, while the middle and right-hand

panels show the dependences on the concentrations of the subhalo, $c_s$, and the host halo, $c_h$, respectively. In all cases, we adopt a host halo-to-subhalo mass ratio of $M_h/M_s = 1000$. As expected, the impact of tidal heating increases if the orbit is more radial (smaller $\eta$) and/or more bound (smaller $x_c$), if the subhalo is less concentrated (smaller $c_s$), or if the host halo is more concentrated (larger $c_h$). It is not uncommon for $\Delta E$ to exceed $|E_b|$ by more than an order of magnitude, especially when $\eta \lesssim 0.5$. However, the stripped mass fraction rarely exceeds 80 per cent. We find a tight relation between $\Delta E/|E_b|$ and $f_{\text{strip}}$, which is reasonably well fit by

$$ f_{\text{strip}} = \exp \left[ -1.65 \left( \frac{\Delta E}{|E_b|} \right)^{-0.44} \right]. \quad (27) $$

The stripped mass fractions predicted from $\Delta E/|E_b|$ using equation (27) are indicated in the lower panels of Fig. 4 as dotted curves, and (except for some of the more extreme concentration parameters) typically agree with the values inferred using the Monte Carlo method described above to better than a few per cent. Using this equation, which holds for NFW subhaloes covering a wide range of concentration parameters and orbital parameters, we see that an impulsive tidal shock with $\Delta E/|E_b| = 1 (10)$ only unbinds roughly 20 (55) per cent of the subhalo mass.
<figure><img src="image_9.png" /><figcaption>Figure 5. The subhalo mass fraction that is stripped off during the first radial period, $f_{\text{strip}}$, as function of the orbital circularity, $\eta$. The subhalo is modelled as a spherical NFW halo with concentration parameter $c_s = 10$, orbiting in a host halo of mass $M_h = 1000 M_s$ and concentration parameter $c_h = 5$, on an orbit with $x_c = 1$. Solid dots indicate the results from our numerical simulations (see Section 3 and Table 1). The solid, ochre curve indicates the fraction of subhalo particles that receive an impulsive, tidal shock, $\Delta E$, that is larger than its binding energy, computed using the method described in Section 4.1.1. Note the good agreement with the simulation results for the more radial orbits with $\eta \lesssim 0.2$. The orange-shaded region indicates the mass fractions stripped during one radial period due to harassment, calculated using the method described in Section 4.2.</figcaption></figure>

The right-hand panel of Fig. 3 shows the distribution of $f_{\text{strip}}$, computed using the Monte Carlo method described above, for the newly accreted subhaloes in the Bolshoi simulation. Ochre and green histograms correspond to the bound fractions obtained with and without the adiabatic correction, respectively, which has little impact. Under the instantaneous mass-loss, approximation tidal heating due to the first peri-centric passage is expected to strip on average about 20–30 per cent of the mass of a subhalo. Note that there are no subhaloes for which $f_{\text{strip}} > 0.9$.

### 4.1.2 Comparison with idealized simulations

We now compare our semi-analytical predictions to the results of our idealized numerical simulations described in Section 3. We thus consider subhaloes with mass $M_s$ and concentration parameter $c_s = 10$ orbiting in a host halo of mass $M_h = 1000 M_s$ and concentration parameter $c_h = 5$. We adopt an orbital energy characterized by $x_c = 1$ and vary the orbital circularity $\eta$. The solid ochre curve in Fig. 5 shows the mass fraction that is stripped due to tidal heating during the subhalo’s first radial orbital period, computed under the instantaneous mass-loss approximation as the fraction of particles with $(\Delta E)_i > |E_i|$. For comparison, the solid dots show the stripped mass fraction in the simulations after one radial period.

For $\eta \lesssim 0.2$, the semi-analytical predictions are in excellent agreement with the simulation results, indicating that the evolution of the bound fraction for subhaloes on relatively radial orbits is governed by tidal shocking and that the impulse approximation combined with the instantaneous mass-loss approximation accurately describes the resulting mass evolution of subhaloes. For less radial orbits, tidal heating by itself underpredicts the stripped mass fraction. This is expected, since the changes in the tidal field become less ‘impulsive’ and the impact of tidal stripping is expected to become more important.

In Paper III in this series (Ogiya et al. in preparation), we compare our analytical estimates presented here to a large suite of (properly converged) simulations, covering a much larger range in orbital energies, orbital angular momentum, and halo concentrations.

## 4.2 Impulsive heating due to subhalo–subhalo encounters

In addition to experiencing tidal shocks during peri-centric passages, subhaloes also undergo tidal heating due to high-speed (impulsive) encounters with other subhaloes. The cumulative impact of many such high-speed impulsive encounters is sometimes called ‘harassment’ (Moore et al. 1996a).

In this section, we estimate how the tidal heating due to harassment compares to that due to a peri-centric passage. Consider a target subhalo (hereafter the subject) of mass $M_s$ and size $r_s$ on a radial orbit ($\eta = 0$) in a host halo of mass $M_h$ and size $R_h$. Assume for simplicity that the host halo has a uniform distribution of subhaloes with mass function $dn/dm$, which are orbiting the host with an orbital velocity dispersion $\sigma_h \simeq (GM_h/2R_h)^{1/2}$ (a good approximation for the global 1D velocity dispersion of an isotropic NFW halo).

In order to calculate the total tidal heating $\Delta E$ of our target subhalo per radial period due to harassment, we first consider a single encounter with a perturbing subhalo of mass $m$ and size $r_m$. Let $b$ and $V$ be the impact parameter and relative velocity of the encounter, which we assume to follow a straight path $\vec{R}_{\text{orb}}(t) = (b, Vt, 0)$ with respect to the centre of the subject mass. As shown by GHO99, in the distant tide approximation ($b \gg r_s$), the velocity kick of a particle located at $\vec{r} = (x, y, z)$ is given by

$$ \Delta \vec{v} = \frac{2 G m}{b^2 V} \times \{(3J_0 - J_1 - I_0)x, (2I_0 - I_1 - 3J_0 + J_1)y, -I_0 z\} \quad (28) $$

where

$$ I_k(b) = \int_{1}^{\infty} \mu_k(b\xi) \frac{d\xi}{\xi^2(\xi^2 - 1)^{1/2}} \quad (29) $$

and

$$ J_k(b) = \int_{1}^{\infty} \mu_k(b\xi) \frac{d\xi}{\xi^4(\xi^2 - 1)^{1/2}} \quad (30) $$

with $\mu_0(x)$ and $\mu_1(x)$ given by equations (18) and (19), respectively. Integrating over the spherically symmetric subject mass yields a total tidal heating of

$$ \Delta E_{\text{dt}}(b, m) = \frac{4G^2m^2}{b^4V^2} M_s \frac{\langle r^2 \rangle_s}{3} \chi_{\text{st}}(b). \quad (31) $$

where the subscript indicates that this is only valid under the distant tide approximation,

$$ \chi_{\text{st}} = \frac{1}{2} \left[ (3J_0 - J_1 - I_0)^2 + (2I_0 - I_1 - 3J_0 + J_1)^2 + I_0^2 \right], \quad (32) $$

and

$$ \langle r^2 \rangle_s = \frac{4\pi}{M_s} \int_0^{r_s} \rho_s(r) r^4 dr, \quad (33) $$

For a head-on collision (i.e. $b = 0$), the cylindrical symmetry assures that the velocity kick for a particle only depends on, and is pointing in the direction of, the cylindrical radius $R = \sqrt{x^2 + z^2}$. In particular,

$$ \Delta v_R = 2 \int_0^\infty \frac{Gm(r)}{r^2} \frac{R}{r} dt = \frac{2 Gm}{V} \frac{I_0(R)}{R}, \quad (34) $$
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Left-hand panel: the tidal heating per crossing time due to impulsive subhalo–subhalo encounters (harassment). Different curves plot the logarithm of $(\Delta E)_{\text{cross}}/|E_b|$ (equation 40) as function of the logarithm of the ratio of the mass of the target subhalo, $M_s$, to the mass of the host halo, $M_h$, for three different values of the assumed concentration of DM subhaloes, as indicated. Right-hand panel: the logarithm of the corresponding stripped mass fraction, $f_{\text{strip}}$, computed under two extreme conditions, where $\Delta E$ is dominated by a single encounter (solid lines) and by many independent encounters (dashed lines). Note how tidal harassment due to subhalo–subhalo encounters is expected to only strip between 1 and 10 per cent of the subhalo’s mass per crossing time. This is significantly less than the mass fraction stripped due to the first peri-centric passage of the host halo (cf. Figs 3 and 4).</figcaption></figure>

which implies a tidal heating

$$ \Delta E_{\text{ho}}(m) = \frac{4G^2m^2}{V^2} \pi \int_0^{r_s} I_0^2(R) \Sigma_s(R) \frac{dR}{R} \quad (35) $$

with $\Sigma_s(R)$ the projected surface density profile of the subject mass, and the subscript indicates that this is the heating for a head-on collision. Throughout, we assume that subhaloes have NFW density profiles,<fnref n="5" /> for which

$$ \Sigma_s(R) = \frac{M_s}{2\pi r_0^2} \frac{1}{f(c_s)} g(R/r_0) \quad (36) $$

with $r_0 = r_s/c_s$ the scale radius, and

$$ g(x) = \begin{cases} \frac{1}{x^{2-1}} \left[ 1 - \frac{2}{\sqrt{1-x^2}} \operatorname{arctanh} \sqrt{\frac{1-x}{1+x}} \right] & x < 1 \\ \frac{1}{3} & x = 1 \\ \frac{1}{x^{2-1}} \left[ 1 - \frac{2}{\sqrt{x^2-1}} \operatorname{arctan} \sqrt{\frac{x-1}{x+1}} \right] & x > 1. \end{cases} \quad (37) $$

We may assume that $\Delta E(b, m)$ for any impact parameter $b$ is a smooth interpolation between $\Delta E_{\text{dt}}(b, m)$, which only holds for $b \gg r_s$ and $\Delta E_{\text{ho}}(m)$, which holds for $b = 0$. In fact, we won’t make a big error if we simply set $\Delta E(b, m)$ to be the minimum of $\Delta E_{\text{dt}}(b, m)$ and $\Delta E_{\text{ho}}(m)$, i.e.

$$ \langle \Delta E \rangle(b, m) = \frac{4G^2m^2M_s}{3V^2} \langle r^2 \rangle_s \begin{cases} \frac{\chi_{\text{st}}(b)}{b^4} & \text{if } b > b_0 \\ \frac{\chi_{\text{st}}(b_0)}{b_0^4} & \text{if } b \le b_0. \end{cases} \quad (38) $$

Here, $b_0$ is defined by $\Delta E_{\text{dt}}(b_0, m) = \Delta E_{\text{ho}}(m)$; i.e. $b_0$ is the impact parameter for which the heating computed using equation (31) is equal to that of a head-on collision (equation 35) and can be solved for using a simple root finder.

Under the assumption that the subhaloes are uniformly distributed within the host halo, the total heating experienced by our target subhalo during one passage through the host halo is given by

$$ (\Delta E)_{\text{cross}} = 2\pi \int_{0}^{M_h} \frac{dm}{m} \frac{dn}{dm} \int_{0}^{R_h} \Delta E(b, m) \frac{b \, db}{\pi R_h^2}. \quad (39) $$

Substituting equation (38) and using that $V^2 \simeq 2\sigma_h^2$ yield

$$ \frac{(\Delta E)_{\text{cross}}}{|E_b|} = \frac{16}{3} \left(\frac{M_h}{M_s}\right)^{2/3} \frac{\langle r^2 \rangle_s}{f_E(c_s)} \times \int_0^1 d\psi \, \psi \frac{dn}{d\ln\psi} \left[ \frac{\chi_{\text{st}}(b_0)}{2b_0^2} + \int_{b_0}^{R_h} \frac{\chi_{\text{st}}(b)}{b^3} db \right]. \quad (40) $$

Here $\psi = m/M_h$, and we have made the additional assumption that $V_h/V_s = (M_h/M_s)^{1/3}$.

The left-hand panel of Fig. 6 shows $(\Delta E)_{\text{cross}}/|E_b|$ as function of the mass ratio $M_h/M_s$, computed with equation (40) using the universal fitting function for the (evolved) subhalo mass function of Jiang & van den Bosch (2016), and adopting a typical subhalo mass fraction of 10 per cent. For simplicity, we assume that all subhaloes (subject as well as perturbers) have the same NFW concentration parameter, $c$. Results are shown for three values of $c$, as indicated. Note that the expected amount of tidal heating due to harassment experienced by a subhalo is typically between 10 and 50 per cent of its binding energy, per crossing time. Typically, as one might expect, more massive subhaloes experience less harassment, but the mass dependence is fairly weak.

As discussed in Section 4.1.1, it is difficult to gauge the impact of $\Delta E$ on a subhalo without addressing how this energy is distributed over the $N$ constituent particles. In the case of harassment, this is difficult to address with accuracy. The reason is that particles in the subject mass move in between separate collisions with individual perturbing subhaloes. However, we may get some insight using two limiting cases. If $\Delta E$ is dominated by the impulsive shock from a single encounter (with a massive subhalo and/or a small impact parameter), then we are in the limit where the energy gain for particle $i$ scales as $(\Delta E)_i \propto r_i^2$ and we can use equation (27) to estimate the stripped mass fraction. The results are indicated as solid lines in the right-hand panel of Fig. 6. If, on the other hand, $\Delta E$ is

<footnote n="5">⁵ We acknowledge that this is not a particularly accurate assumption (e.g. Hayashi et al. 2003), but the detailed density profile of subhaloes has little impact on what is anyways a fairly crude estimate.</footnote>
dominated by the cumulative effect of many encounters (with small subhaloes, covering the entire range of impact parameters), then we are in the regime where $(\Delta E)_i = (\Delta E)_{cross}/N$; i.e. the energy gain is independent of location within the subhalo (see Section 6.3). In this limit, we can compute the stripped mass fraction using a similar Monte Carlo approach as in Section 4.1.1 and simply counting the fraction of particles for which $(\Delta E)_i > E_i$. This results in the stripped mass fractions indicated by dashed lines in the right-hand panel. Interestingly, the two extreme cases yield stripped mass fractions that are relatively similar, spanning the range 0.01 $\lesssim f_{strip} \lesssim 0.1$. This range is indicated as an orange-shaded region in Fig. 5.

Although the above estimate for the tidal heating due to harassment is fairly crude, it is clear that the impact of harassment is subdominant to that due to the subhalo's first peri-centric passage. We therefore conclude that the tidal heating due to subhalo-subhalo encounters can safely be neglected when trying to assess the survivability of DM substructure.

This conclusion is at odds with a study by Knebe et al. (2006), who claim that harassment may contribute as much as 40 per cent to the total mass-loss rate of a subhalo. However, their estimate is based on comparing the force on the subhalo due to the host halo to that due to all other subhaloes, without taking account of the actual tides (relevant for tidal stripping) or the relative velocities (relevant for tidal heating). Although we believe our estimate to be more accurate, a more detailed study is required to settle this disagreement, which is beyond the scope of this paper.

## 4.3 Tidal stripping

In the previous two subsections, we focused on tidal shocking, which manifests itself whenever the tidal field changes rapidly. We now turn to the impact of a slowly varying tidal field. This is the situation one encounters for subhaloes on close to circular orbits. In this limit, one expects the material outside of the tidal radius to be stripped off. Modelling this tidal stripping, however, is far from trivial. First of all, as discussed in Section 2, the tidal radius is a poorly defined concept, with significantly different definitions in use. Secondly, it is unclear a priori at what *rate* the material that is located outside of the (instantaneous) tidal radius is going to be stripped off.

Effectively, all semi-analytical models for the tidal evolution of subhaloes adopt a tidal stripping rate given by

$$ \frac{dm}{dt} = \frac{m(>r_t)}{\tau_{\text{orb}}/\alpha}. \qquad (41) $$

Here, $r_t$ is the instantaneous tidal radius, $\tau_{orb} = 2\pi/\omega$ with $\omega$ the instantaneous angular velocity of the subhalo, and $\alpha$ is a 'free parameter' that differs substantially from author to author. Whereas Taylor & Babul (2001, 2004), Taffoni et al. (2003), and Zentner & Bullock (2003) all adopt $\alpha = 1$, Zentner et al. (2005) and Pullen et al. (2014) tune $\alpha$ by matching their predicted subhalo mass function to simulation results. This yields $\alpha = 3.5$ and 2.5, respectively.<fnref n="6" /> For comparison, Diemand et al. (2007) find that the mass evolution of subhaloes in their 'Via Lactea' simulation is best modelled using equation (41) with $\alpha \sim 6$. Finally, some authors have assumed that the stripping is instantaneous (e.g. Oguri & Lee 2004; Peñarrubia & Benson 2005), which effectively implies $\alpha = \infty$. Since the tidal radius is typically minimal at the orbit's peri-centre, this implies

that per radial orbital period, all the mass outside of the peri-tidal radius, $r_{t,lp}$, will be stripped. This is the same assumption that underlies the orbit-averaged models of van den Bosch et al. (2005) and Jiang & van den Bosch (2016, 2017), which therefore effectively have adopted equation (41) with $\alpha = \infty$.

In this section, we investigate what all these different assumptions imply for the amount of mass that is stripped from a typical DM subhalo during its first orbit after infall. Practically, the way we compute this is as follows. Using orbit-integration, we discretize the subhalo's orbit in 2000 equal time steps $\Delta t$. Let $M_n$, $r_{t,n}$, and $\omega_n$ correspond to the instantaneous bound subhalo mass, the instantaneous tidal radius, and the instantaneous angular velocity of the subhalo at time step $n$. If we ignore for the moment that the subhalo will respond to tidal stripping by re-virialization and simply assume that the mass distribution of the bound material remains identical to the initial mass distribution (prior to any stripping), then

$$ M_{n+1} = M_n - \alpha \dot{M}_n \Delta t, \qquad (42) $$

where

$$ \dot{M}_n = \begin{cases} \frac{[M_n - M(<r_{t,n})]}{(2\pi/\omega_n)} & \text{if } M_n > M(<r_{t,n}) \\ 0 & \text{otherwise.} \end{cases} \qquad (43) $$

Using this method, and the properties and orbits of the 5000 subhaloes at accretion from the Bolshoi simulation depicted in Fig. 1, we obtain the stripped mass fractions shown in Fig. 7. The three panels differ in the value of the parameter $\alpha$, as indicated, while the red and blue histograms correspond to the results obtained using $r_{t,1}$ (equation 3) and $r_{t,2}$ (equation 4), respectively. Note that $\alpha = \infty$ (right-hand panel) corresponds to instantaneous stripping, in which case the stripped mass fraction is simply given by $f_{strip} = M_s(>r_{t,lp})/M_s$.

In general, since $r_{t,2} < r_{t,1}$, accounting for the centrifugal force in the definition of the tidal radius results in predictions for the stripped mass fraction that are significantly larger. Also, $f_{strip}$ depends strongly on the value of $\alpha$: for $\alpha = 1$, the median $f_{strip}$ for subhaloes during their first radial orbit is only 0.10 (0.16) for $r_t = r_{t,1}$ ($r_{t,2}$). For $\alpha = \infty$, this increases to 0.38 (0.63). Note that the assumption of instantaneous stripping combined with $r_t = r_{t,2}$ results in an appreciable fraction of subhaloes with extremely large stripped mass fractions; among our sample of 5000 newly accreted subhaloes 3.4 (1.2) per cent have $f_{strip} > 0.99$ ($>0.999$). Clearly, depending on how one defines the tidal radius and what one assumes regarding the rate at which material beyond the tidal radius is stripped, the predicted impact of tidal stripping varies dramatically. This is also evident from Fig. 8, which compares our semi-analytical predictions for $f_{strip}$ (for the first radial orbital period), to the results from our idealized numerical simulations. The two panels differ in the definition of the tidal radius, as indicated, while different curves in each panel correspond to different values of $\alpha$.

Clearly, no single combination of tidal radius definition and value for $\alpha$ results in predictions for $f_{strip}$ that are in good agreement with the simulation results over the entire range in $\eta$. It is not surprising that tidal stripping is unable to capture the tidal evolution of subhaloes on highly radial (i.e. low $\eta$) orbits, as their evolution is clearly governed by impulsive, tidal heating (cf. Fig. 5). However, it may seem surprising that the tidal stripping model also fails to describe the results along close-to-circular orbits. After all, this is the regime where the tidal field varies slowly, and thus where the tidal stripping model is expected to be most accurate. However, there are several problems with the tidal stripping model that are ultimately responsible for its failure.

<footnote n="6">⁶ As indicated in footnote 13 of Diemand et al. (2007), there is a factor $2\pi$ missing in the formulation of Zentner et al. (2005).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Distributions of the mass fraction $f_{\text{strip}}$ that is stripped off during a subhalo's first radial orbital period due to tidal (non-impulsive) stripping for the same set of newly accreted subhaloes as in Fig. 1. Red and blue histograms are based on the different definitions of the tidal radius, as indicated, while different panels correspond to different stripping rates, as parametrized by $\alpha$ (cf. equation 41). The case where $\alpha = \infty$ (right-hand panel) corresponds to instantaneous stripping, for which $f_{\text{strip}} = M_s(r_{\text{ip}})/M_s$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Same as Fig. 5, except that this time the curves reflect the predicted mass fractions that are stripped off subhaloes during their first radial orbit due to (non-impulsive) tidal stripping. The two panels differ in the adopted definition for the tidal radius, as indicated, while solid, dashed, and dot-dashed curves correspond to assumed stripping rates with $\alpha = \infty$, 3.5, and 1.0, respectively. As in Fig. 5, the solid dots indicate the results from the idealized numerical simulations (see Section 3 and Table 1). See the text for a detailed discussion.</figcaption></figure>

First of all, the notion that only the matter that is initially located outside of $r_t$ will be stripped off is clearly incorrect. After all, some of the particles inside of the tidal radius will be on orbits whose apo-centre $r_{\text{apo}} > r_t$; Hence, if the tidal radius does not evolve much with time (i.e. if the orbit is close to circular), one expects those particles to be stripped as well. And since the dynamical time of particles in the subhalo is comparable to the orbital time of the subhalo within its host, the time-scale over which this stripping occurs should be comparable to the radial orbital period $T_r$. Using the same Monte Carlo realization used to compute the stripped mass fraction due to tidal heating (see Section 4.1.1), we compute the apo-centre for each of the 50 000 subhalo particles. For $\eta = 1$ (i.e. circular orbit), we find that 41 per cent of all subhalo particles have an apo-centre $r_a > r_{t, 2}$. This is in excellent agreement with the stripped mass fraction in the simulation, which has $f_{\text{strip}} = 0.394$ (after one radial period) for $\eta = 1$ (cf. Table 1), and suggests that along orbits that are close to circular, one may estimate

$f_{\text{strip}}$ as the fraction of particles with $r_a > r_{t, 2}$. Note, though, that this approximation rapidly deteriorates for decreasing $\eta$.<fnref n="7" /> Furthermore, there is an additional problem with the stripping model that is even more difficult to overcome. After the outer layers of a subhalo are (instantaneously) removed, the remaining remnant is no longer in virial equilibrium. As we show in Appendix B, the remnant has a virial ratio $K/|W| > 0.5$. Hence, the only way the system can re-virialize is by converting kinetic to potential energy. In doing so, the system ‘puffs up’; i.e. its characteristic radius increases while its overall density decreases. This, in turn, results in a decrease of the tidal radius, and hence in additional mass-loss (Kampakoglou & Benson 2007). As a result, the mass-loss is a continuous process, even in the static field experienced along a circular orbit. This is evident from the left-hand panel of Fig. 2 (red curve), which shows the

<footnote n="7">⁷ For $\eta = 0.9$, it already overpredicts $f_{\text{strip}}$ by more than 30 per cent.</footnote>
evolution of the bound fraction for the duration of one radial period, $T_r = 9.53$ Gyr. In Paper II, we show that the mass-loss continuous well beyond that, and that subhaloes on circular orbits continue to loose mass even after 60 Gyr of evolution.

## 5 CAN INSTANTANEOUS STRIPPING RESULT IN DISRUPTION?

In an influential paper, Hayashi et al. (2003) pointed out that instantaneously removing all material from a halo down to some truncation radius, $r_{tr}$, may leave a remnant with positive binding energy, as long as $r_{tr} < r_{crit}$. For a spherical, isotropic NFW halo, $r_{crit}$ is 0.77 times the NFW scale radius. This is demonstrated in Appendix B, where it is also shown how the critical radius, $r_{crit}$, depends on the halo's anisotropy. Typically, $r_{crit}$ becomes larger for haloes that are more radially anisotropic, and in general $r_{crit} > 0$, unless haloes are strongly tangentially anisotropic.

This seems to suggest, as was advocated by Hayashi et al., that a subhalo whose tidal radius is smaller than $r_{crit}$ would 'spontaneously' disrupt, even without the need to invoke tidal shock heating. And since a typical NFW halo has roughly between 5 and 10 per cent of its mass enclosed within this critical radius, this suggests that a typical subhalo should physically disrupt once it has lost between 90 and 95 per cent of its original (i.e. at accretion) mass. Interestingly, this is not inconsistent with the results from numerical simulations, which indeed reveal very few subhaloes whose mass is less than 5 per cent of that at accretion (e.g. Han et al. 2016; Jiang & van den Bosch 2016; van den Bosch 2017), and various authors have implemented a treatment of subhalo disruption based on this notion (e.g. Zentner & Bullock 2003; Taylor & Babul 2004; Klypin et al. 2015).

However, as we now demonstrate, this notion is seriously flawed, and, if simulated with sufficient resolution, the instantaneous removal of outer layers always leaves a bound remnant, even when these outer layers contain more than 99.9 per cent of the total virial mass. The reason is once again that the total binding energy is not sufficiently informative. What matters is the distribution of binding energies of the constituent particles, not its sum. Note that we already encountered an example of subhaloes with positive binding energy, namely subhaloes immediately following an impulsive shock with $\Delta E/|E_b| > 1$. As we saw in Section 4.1 neither of these disrupted, even when $\Delta E > 100|E_b|$. In order to test how subhaloes respond to instantaneous mass stripping down to radii below $r_{crit}$, we now use numerical N-body simulations to evolve instantaneously stripped, isotropic NFW haloes and track their bound mass fraction over time.

### 5.1 Numerical simulations

We simulate the evolution of isolated, spherical NFW host haloes that are instantaneously truncated. We use the same method as in Section 3 to draw the initial positions and velocities of the particles: we assume that, prior to truncation, the halo has an isotropic velocity distribution, and we truncate the halo at a radius $r_{tr}$, which is a free parameter, while initializing the velocities as if the halo extends out to infinity (this is the expected result if the stripping is instantaneous). All our simulations have $N = 10^5$ particles inside the truncation radius, and are carried out using a modified version of the hierarchical N-body code treecode, written by Joshua Barnes with some modifications due to John Dubinski. treecode uses a Barnes & Hut (1986) octree to compute accelerations based on a

<table><caption>**Table 2.** Parameters of simulations of instantaneous truncation.</caption><thead><tr><th>ID</th><th>r<sub>tr</sub>/r<sub>0</sub></th><th>f(r<sub>tr</sub>)</th><th>ε/r<sub>0</sub></th><th>f<sub>b, rem</sub><br>(t = 0)</th><th>f<sub>b, rem</sub><br>(50 Gyr)</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>(5)</th><th>(6)</th></tr></thead><tbody><tr><td>A01</td><td>7.11</td><td>0.817</td><td>0.05</td><td>0.98</td><td>0.98</td></tr><tr><td>A02</td><td>6.46</td><td>0.768</td><td>0.05</td><td>0.98</td><td>0.98</td></tr><tr><td>A03</td><td>5.81</td><td>0.715</td><td>0.05</td><td>0.98</td><td>0.97</td></tr><tr><td>A04</td><td>5.16</td><td>0.659</td><td>0.05</td><td>0.98</td><td>0.97</td></tr><tr><td>A05</td><td>4.51</td><td>0.596</td><td>0.05</td><td>0.97</td><td>0.97</td></tr><tr><td>A06</td><td>3.85</td><td>0.528</td><td>0.05</td><td>0.97</td><td>0.96</td></tr><tr><td>A07</td><td>3.19</td><td>0.451</td><td>0.05</td><td>0.96</td><td>0.95</td></tr><tr><td>A08</td><td>2.52</td><td>0.365</td><td>0.05</td><td>0.94</td><td>0.93</td></tr><tr><td>A09</td><td>1.83</td><td>0.265</td><td>0.05</td><td>0.91</td><td>0.90</td></tr><tr><td>A10</td><td>1.48</td><td>0.209</td><td>0.05</td><td>0.89</td><td>0.88</td></tr><tr><td>A11</td><td>1.10</td><td>0.147</td><td>0.05</td><td>0.83</td><td>0.82</td></tr><tr><td>A12</td><td>0.68</td><td>0.076</td><td>0.05</td><td>0.71</td><td>0.71</td></tr><tr><td>A13</td><td>0.48</td><td>0.046</td><td>0.05</td><td>0.61</td><td>0.62</td></tr><tr><td>A14</td><td>0.35</td><td>0.027</td><td>0.05</td><td>0.46</td><td>0.49</td></tr><tr><td>A15</td><td>0.20</td><td>0.011</td><td>0.05</td><td>0.23</td><td>0.31</td></tr><tr><td>A16</td><td>0.10</td><td>0.003</td><td>0.05</td><td>0.15</td><td>0.18</td></tr><tr><td>A17</td><td>0.08</td><td>0.001</td><td>0.05</td><td>0.00</td><td>Disrupt</td></tr><tr><td>B01</td><td>2.52</td><td>0.365</td><td>0.0125</td><td>0.94</td><td>0.93</td></tr><tr><td>B02</td><td>1.83</td><td>0.265</td><td>0.0092</td><td>0.91</td><td>0.90</td></tr><tr><td>B03</td><td>1.48</td><td>0.209</td><td>0.0074</td><td>0.89</td><td>0.88</td></tr><tr><td>B04</td><td>1.10</td><td>0.147</td><td>0.0051</td><td>0.83</td><td>0.82</td></tr><tr><td>B05</td><td>0.68</td><td>0.076</td><td>0.0034</td><td>0.73</td><td>0.72</td></tr><tr><td>B06</td><td>0.35</td><td>0.027</td><td>0.0017</td><td>0.55</td><td>0.57</td></tr><tr><td>B07</td><td>0.20</td><td>0.011</td><td>0.0010</td><td>0.43</td><td>0.46</td></tr><tr><td>B08</td><td>0.10</td><td>0.003</td><td>0.0005</td><td>0.34</td><td>0.37</td></tr><tr><td>B09</td><td>0.08</td><td>0.001</td><td>0.0004</td><td>0.31</td><td>0.35</td></tr></tbody></table>

Note. Simulations of instantaneously truncated NFW haloes. Listed are the simulation ID (column 1), the ratio of the truncation radius, $r_{tr}$, to the scale radius $r_0$ (column 2), the fraction $f(r_{tr}) = M(r_{tr})/M_{vir}$ of the virial mass located inside the truncation radius (column 3), the softening length in units of the scale radius (column 4), and the bound fractions of the remnants, $f_{b, rem}$, at $t = 0$ (column 5) and $t = 50$ Gyr (column 6). Simulations A01–A17 belong to TruncA, and all adopt a softening length of $\epsilon = 0.05$, while the TruncB simulations B01–B09 adopt a softening length that scales linearly with the truncation radius according to $\epsilon = 0.05(r_{tr}/r_{vir})$.

multipole expansion up to quadrupole order, and uses a straightforward second-order leap-frog integration scheme to solve the equations of motion. Since we use fixed time steps, our integration scheme is fully symplectic. Forces between particles are softened using a simple Plummer softening.

As in Section 3, we adopt model units in which the gravitational constant, $G$, the scale radius, $r_0$, and the virial mass of the halo, $M_{vir}$, are all unity. Using numerical simulations of NFW haloes in isolation (see Paper II), we infer an optimal softening length of $\epsilon = 0.05$ (in model units) for simulations with $N = 10^5$ particles inside the virial radius. In order to gauge the sensitivity to the specific choice of the softening length, we run two different sets of simulations: in set TruncA we adopt a softening length of $\epsilon = 0.05$, independent of $r_{tr}$, while in set TruncB, we scale the softening length linearly with the size of the simulated system according to $\epsilon = 0.05(r_{tr}/r_{vir})$; such a scaling is consistent with simple expectations for the optimal softening length (e.g. van Kampen 2000a; Power et al. 2003). We run each truncated halo in isolation for 40 000 time steps of $\Delta t = 0.02$, corresponding to a total integration time of $\sim 50$ Gyr. Finally, we adopt a tree opening angle of $\theta = 0.7$, and we have verified that our results are extremely stable to changes in $\Delta t$ and/or $\theta$. The parameters of the simulations in TruncA and TruncB are listed in Table 2.
<figure><img src="image_13.png" /><figcaption>Figure 9. The bound fraction, $f_{b,tot}$, for NFW haloes with concentration parameter $c = 10$ as a function of their truncation radius, $r_{tr}$, expressed in units of the scale radius $r_0$. The green solid line indicates $M(r_{tr})/M_{vir}$ and reflects the initial mass fraction enclosed by $r_{tr}$. The open circles are the mass fractions that are found to be bound in numerical simulations with $N = 10^5$ particles after 50 Gyr of evolution. Blue and red circles correspond to simulations from TruncA and TruncB, respectively, that only differ in their choice of softening length, as indicated. The dashed, vertical line corresponds to $r_{tr} = 0.77r_0$, where the binding energy of the remnant transits from being negative to being positive. The blue, downward arrow indicates that the TruncA simulation with $r_{tr} = 0.078r_0$ disrupts (i.e. has $f_{bound} \to 0$).</figcaption></figure>

## 5.2 Results

We use the method described in Appendix A to investigate how the bound-mass of the truncated halo evolves over time. We distinguish two bound fractions:

$$f_{b,rem}(t) \equiv \frac{M(t)}{M(r_{tr})} = \frac{N_{bound}}{N_{tot}}, \qquad (44)$$

and

$$f_{b,tot}(t) \equiv \frac{M(t)}{M_{vir}} = f_{b,rem}(t) \times \frac{M(r_{tr})}{M_{vir}}. \qquad (45)$$

Here, $M(t)$ is the bound mass of the halo at time $t$, $M(r)$ is the initial (prior to truncation) mass of the subhalo inside radius $r$, $M_{vir} = M(r_{vir})$ is the subhalo's initial virial mass, $N_{bound}$ is the number of bound particles at time $t$, and $N_{tot} = 10^5$ is the total number of particles used in the simulation. Thus, $f_{b,rem}$ is the fraction of the mass initially located inside the truncation radius that remains bound, while $f_{b,tot}$ is the fraction of mass that remains bound compared to the original virial mass of the halo in question.

Fig. 9 summarizes the simulation results (see also Table 2). The open circles indicate $f_{b,tot}$ at the end of each simulation (i.e. after 50 Gyr of evolution). Blue and red circles correspond to simulations from TruncA and TruncB, respectively, and the results are shown as a function of the truncation radius, $r_{tr}$, in units of the halo's (original) scale radius, $r_0$. The green, solid line indicates $M(r_{tr})/M_{vir}$, and reflects the initial mass fraction enclosed by $r_{tr}$. The dashed, vertical line marks $r_{tr} = 0.77r_0$, which is the critical value for the truncation radius, to the left of which the remnants (immediately after the instantaneous truncation) have positive binding energy. Clearly, there is no indication of anything special happening around this scale. In fact, in all simulations we find that part of the remnant remains bound, and quickly resettles into a new, virialized halo. As

<figure><img src="image_14.png" /><figcaption>Figure 10. The density profiles (upper panel), enclosed mass profiles (middle panel), and radial velocity dispersion profiles (lower panel) at five different evolutionary stages for an NFW halo that is instantaneously truncated (at $t = 0$) at $r_{tr} = 0.68r_0$ (simulation B05 in Table 2). Note that the profiles are for the bound particles only.</figcaption></figure>

an example, Fig. 10 shows the evolution in the density profile (upper panel), enclosed mass profile (middle panel), and radial velocity dispersion profile (lower panel) for an NFW halo that is instantaneously truncated at $r_{tr} = 0.68r_0$. The inner region quickly (within ~2 Gyr) re-virializes, while it takes more than a Hubble time for the weakly bound particles to build a virialized structure that extends to $> 10r_0$ (i.e. beyond the original virial radius of the subhalo).

In general, a smaller truncation radius implies a larger fraction of particles inside this truncation radius that are unbound (and escape), but in each case a bound subset remains. The only exception is simulation A17 from TruncA, which has $r_{tr} = 0.078r_0$ and which experiences complete disruption immediately after the onset of the simulation (indicated by a blue, downward pointing arrow in Fig. 9). We point out, though, that in this case $r_{tr}$ is only ~1.5 times the softening length; clearly, the potential is ‘over-softened’ and has little to no resemblance to the central region of an NFW profile. The same simulation in TruncB (simulation B09), which uses a much smaller softening length (see Table 2), results in ~35 per cent of the initial particles of the remnant surviving as a self-bound entity.

In all simulations, we find that the bound fraction does not evolve significantly with time: the fraction of particles that remain bound
together is basically the same after more than 50 Gyr of evolution, as it is at $t = 0$ (cf. columns (5) and (6) of Table 2). Only when $r_{tr} \lesssim r_0/2$ do we find that the bound fraction evolves somewhat (typically $f_b, \text{rem}$ increases by a few per cent), during the first $\sim 0.2$ Gyr after the onset of the simulation. This is a manifestation of (violent) relaxation in response to the instantaneous removal of the halo's outer layers; as shown in Appendix B, the remnant can be far from virial equilibrium.

To summarize, unless DM haloes are strongly tangentially anisotropic, instantaneously stripping matter from the outskirts may leave a remnant with total, positive binding energy. However, contrary to naive expectations (e.g. Hayashi et al. 2003; Taylor & Babul 2004; Klypin et al. 2015), a system with positive binding energy does not necessarily disrupt entirely. After all, the system will typically have a broad distribution of binding energies, and a significant fraction of the particles can still be bound, even if the total binding energy is positive. In the N-body experiments described above, we only find that a truncated, isotropic NFW halo disrupts if the truncation radius is less than, or comparable to, the softening length used. Although we have only tested this for isotropic NFW haloes, we emphasize that (the central parts) of NFW haloes in cosmological simulations are indeed close to isotropic (e.g. Navarro et al. 2010). We therefore conclude that instantaneous truncation of an NFW halo does not lead to halo disruption. Artificial disruption in simulations may occur, however, if the system is simulated with insufficient force resolution (see Section 6.4 and Paper II).

# 6 NUMERICAL DISRUPTION PROCESSES

As discussed in Section 1, it was not until 1997 that numerical N-body simulations started to resolve surviving populations of DM substructure. Prior to this, the simulations suffered from a serious overmerging problem, mainly due to inadequate force softening (Moore et al. 1996b). However, numerical overmerging (i.e. the artificial disruption of substructure in numerical simulations) is still present in modern state-of-the-art simulations. This is the reason why methods to populate DM-only simulations with galaxies, such as semi-analytical models for galaxy formation (e.g. Springel et al. 2001; Kang et al. 2005; Kitzbichler & White 2008), subhalo abundance matching models (e.g. Conroy et al. 2006; Guo & White 2014; Campbell et al. 2017), and empirical models (e.g. Moster et al. 2010, 2017; Tollet et al. 2017), often include ‘orphan’ galaxies, i.e. mock galaxies without an associated subhalo in the simulation. Although it is therefore recognized that simulations continue to experience artificial disruption, it is generally assumed that this happens only for subhaloes below a mass resolution of 50–100 particles. This notion is based on the fact that subhalo mass functions are typically converged down to this mass resolution limit (e.g. Springel et al. 2008; Onions et al. 2012; Knebe et al. 2013; van den Bosch & Jiang 2016). The general consensus therefore is that any disruption of subhaloes above this ‘resolution limit’ must be physical in origin (see Diemand et al. 2004b, for a detailed discussion).

However, as we demonstrate in Paper II, state-of-the-art numerical simulations of CDM structure formation still suffer from abundant artificial disruption well above this mass resolution limit, and it is therefore imperative that we revisit the issue of numerical overmerging. There are a number of processes that may (potentially) give rise to numerical, artificial disruption of DM substructure in numerical N-body simulations. These are (i) evaporation resulting from two-body relaxation of the subhalo, (ii) evaporation due to two-body encounters with particles from the host halo, (iii) tidal

heating due to impulsive encounters with particles from the host halo, and (iv) disruption due to issues with force softening.

Each of these processes have been discussed in more or less detail in various previous studies, including Carlberg (1994), van Kampen (1995, 2000a,b), Moore et al. (1996b), and Klypin et al. (1999a). For the sake of completeness, and in order to correct and elucidate inconsistencies and errors in some of these previous studies, we discuss each of these processes in detail in what follows.

## 6.1 Evaporation due to two-body self-relaxation

In a gravitational N-body system, encounters (also called ‘collisions’) among the particles drive the system towards equipartition of kinetic energy. The time-scale for this two-body relaxation process is roughly

$$t_{\text{relax}} = \frac{N}{8 \ln \Lambda} t_{\text{cross}} \qquad (46)$$

(e.g. Binney & Tremaine 2008). Here, $t_{\text{cross}} \approx r/\sigma$ is the crossing time of the system of size $r$ and velocity dispersion $\sigma$, and $\ln \Lambda = \ln(b_{\text{max}}/b_{\text{min}})$ is the Coulomb logarithm, with $b_{\text{max}} \sim r$ and $b_{\text{min}} = 2Gm_p/\sigma^2$ the maximum and minimum impact parameters, and $m_p$ the particle mass. Using the virial theorem, one then finds that $\ln \Lambda = \ln(N/2) \sim \ln N$.

Two-body relaxation drives the local velocity distribution of the particles towards a Maxwellian, which has a tail that exceeds the local escape speed. Hence, from time to time the gravitational encounters can give enough energy to a particle so that it can escape. This causes the N-body system to evaporate, which ultimately leaves a final system with two particles on a Keplerian orbit. The evaporation time-scale is roughly 140 times the two-body relaxation time (Binney & Tremaine 2008), which for a system of point particles implies

$$t_{\text{evap}} \approx \frac{15 N}{\ln N} t_{\text{cross}}. \qquad (47)$$

In N-body simulations of DM substructure, the number of particles used is vastly smaller (typically by more than 50 orders of magnitude) than the actual number of constituents of the physical system being simulated. Hence, the relaxation and evaporation times in the simulation are much, much smaller than they are in reality, which in principle could give rise to artificial, excessive evaporation.

When considering the impact of two-body relaxation in numerical N-body simulations, one has to modify the above estimate to take account of force softening, which suppresses encounters with large deflection angles. As shown in Farouki & Salpeter (1982), the net effect of force softening is to boost the minimum impact parameter, $b_{\text{min}}$, such that the Coulomb logarithm becomes

$$\Lambda = \min\{N, r/(4\epsilon)\}, \qquad (48)$$

with $\epsilon$ the softening length. Hence, softening increases the two-body relaxation time (and thus also the evaporation time), but by how much depends somewhat on the exact value of $\epsilon$. This is demonstrated in Fig. 11, where the green lines plot $\log [t_{\text{evap}}/t_{\text{cross}}]$ as function of $\log[N]$ for an NFW subhalo with concentration parameter $c = 10$. The grey-shaded region marks time-scales that are less than 30 crossing times, which is a conservative upper limit on the age of the Universe. Hence, time-scales above this grey region are not expected to be of much relevance. In the left-hand panel, we adopt $\epsilon/r_0 = 1.7 \times 10^{-2} (N/10^5)^{-0.23}$, where $r_0 = r_s/c$ is the scale radius of the subhalo. As shown in Dehnen (2001), this scaling for the (Plummer) softening length optimizes between systematic force
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Time-scales, in units of the crossing time, for the disruption of subhaloes due to numerical effects as a function of the number of particles, *N*, used in simulating the subhalo (which is assumed to have a concentration parameter *c* = 10). In particular, green, blue, and red lines correspond to the times scales for evaporation due to two-body relaxation, evaporation due to collisional encounters with host halo particles, and disruption due to impulsive heating by host halo particles. For the latter, the red dashed line shows the proper, corrected time-scale. See the text for details. The grey-shaded region corresponds to $t/t_{cross} \le 30$, which roughly corresponds to time-scales shorter than a Hubble time. Left- and right-hand panels correspond to different assumptions regarding the softening length, as indicated.</figcaption></figure>

errors (which dominate when $\varepsilon$ is large) and force noise (which dominates when $\varepsilon$ is small).<fnref n="8" /> In the right-hand panel, we adopt $\varepsilon/r_0 = 0.126 (N/10^5)^{-0.5}$, which corresponds to the optimal softening length advocated by Power et al. (2003), based on a lower limit required to prevent discreteness effects. For $N \lesssim 2000$, the Power et al. softening results in a much longer evaporation time; the dramatic upturn in $\log[t_{evap}/t_{cross}]$ at small $N$ occurs when the softening length becomes comparable to the size of the subhalo itself; obviously, in such an extreme case, two-body relaxation is completely irrelevant, but at the same time, such a simulation clearly is not able to capture the relevant dynamics. The softening advocated by Dehnen (2001) always has $\varepsilon/r_0 \ll 1$, which results in much smaller evaporation times for small $N$. However, even for $N = 10$, the evaporation time exceeds the Hubble time, and it therefore seems safe to conclude that artificial evaporation of substructure due to two-body relaxation should be of little concern (see also Moore et al. 1996b).

There is one important caveat, though: the above estimate for $t_{relax}$ is based on the classical two-body treatment of Chandrasekhar (1943), in which all collisions are assumed to be independent. This ignores resonant effects and self-gravity ('collective relaxation'), both of which can significantly boost relaxation with respect to the two-body rate (Weinberg 1993). In particular, as Weinberg demonstrates, the dominant contribution to relaxation of a gravitational system comes from large-scale modes, with wavelengths comparable to the size of the system. As a consequence, the treatment of force softening has little to no effect on the actual relaxation rate. This explains, among others, why relaxation rates in expansion-type codes, such as the self-consistent field code, reveal relaxation at a level that is similar to that of tree-based codes (Hernquist & Ostriker 1992), and why the amount of spurious (artificial) fragmentation apparent in warm DM simulations is virtually independent

of the softening length (Power et al. 2016). In Paper II, we demonstrate that these same large-scale fluctuations, in the presence of a tidal field, cause a run-away instability in the mass evolution of DM subhaloes in N-body simulations.

To summarize, we caution that although the two-body relaxation rate (equation 46) may exceed the Hubble time, the actual relaxation time in numerical simulations may be significantly shorter, and the overall impact of relaxation in N-body simulations remains a contentious topic (see e.g. Melott 2007; Hahn & Angulo 2016; Power et al. 2016, and references therein).

## 6.2 Evaporation due to two-body collisions with host particles

In an attempt to explain the origin of numerical ‘overmerging’, Carlberg (1994) suggested that subhaloes in numerical simulations may artificially disrupt due to two-body encounters between particles in the subhalo, and those in the host halo. He argued that, since the host halo particles are so much hotter (i.e. have higher velocities) than the subhalo particles, the former will transfer kinetic energy to the latter, thereby ‘boiling’ the subhalo into dissolution. This idea was further promoted by van Kampen (1995, 2000a,b) who suggested that it is the dominant numerical process to explain the overmerging of substructure in cosmological N-body simulations. Unfortunately, Carlberg (1994) does not give a derivation of his expression for the heating time-scale, which in addition contains a typo, in that the indices ‘c’ and ‘h’ in his equation (11) are swapped. van Kampen (1995) does present a correct derivation for the relaxation time, but then incorrectly assumes that this is comparable to the *disruption* time, thereby underestimating the latter by two orders of magnitude. In order to set the record straight, we give a detailed derivation in what follows.

Consider a dynamically cold system of size $r_c$, consisting of $N_c$ particles of mass *m* and with a velocity dispersion $\sigma_c^2 \simeq G N_c m / r_c$, moving through a hot system of size $r_h$ consisting of $N_h \gg N_c$ particles of identical mass *m*. The velocity dispersion of the hot

<footnote n="8">⁸ This relation, which is taken from table 2 in Dehnen (2001), corresponds to a Hernquist (1990) profile, which is similar, though, to an NFW profile.</footnote>
system is $\sigma_h^2 \simeq G N_h m / r_h$. Each encounter between a hot particle and a cold particle changes the velocity of the cold particle such that

$$ (\Delta v^2)_{\text{enc}} = \frac{4G^2m^2}{b^2V^2}, \qquad (49) $$

(Binney & Tremaine 2008). Here, b is the impact parameter of the encounter and V is the relative speed between the particles. Since the hot particles are moving much faster than the cold ones, we have that $\langle V^2 \rangle = \sigma_h^2 + \sigma_c^2 \simeq \sigma_h^2$. Per crossing through the hot system, each particle in the subhalo has on average about $2(N_h/r_h^2)b \, db$ encounters with impact parameters in the range b, b + db. Integrating over b yields the total $\Delta v^2$ per crossing time, which is

$$ (\Delta v^2)_{\text{cross}} = \frac{8 \ln \Lambda_h}{N_h} \sigma_h^2. \qquad (50) $$

Hence, the relaxation time for the cold system is

$$ t_{\text{relax}}^{\text{ch}} = \frac{\sigma_c^2}{(\Delta v^2)_{\text{cross}}} t_{\text{cross,h}} = \frac{N_h}{8 \ln \Lambda_h} \left( \frac{\sigma_c^2}{\sigma_h^2} \right) t_{\text{cross,h}}, \qquad (51) $$

where the superscript ‘ch’ indicates that this is two-body relaxation of cold particles due to hot particles. Using that the densities of collapsed haloes are independent of mass, we have that $t_{\text{cross,c}} = t_{\text{cross,h}}$, which allows us to rewrite equation (51) as

$$ t_{\text{relax}}^{\text{ch}} = \frac{\sigma_h}{\sigma_c} \frac{\ln \Lambda_c}{\ln \Lambda_h} t_{\text{relax}}^{\text{cc}}, \qquad (52) $$

where $t_{\text{relax}}^{\text{cc}}$ is the two-body relaxation time of the cold system in isolation, given by equation (46), and we have used that $N_h/N_s = (\sigma_h/\sigma_c)^3$, which follows from the virial scaling relations for DM haloes, according to which $\sigma \propto M^{1/3}$.

We thus see that the relaxation time due to encounters with the hot particles is significantly longer than that due to encounters with its own particles. This has its origin in the fact that encounters with host halo particles have, on average, a larger impact parameter and a higher relative velocity (cf. equation 49). This is partially compensated by the fact that there are more host halo particles than subhalo particles, but the combined effect is such that self-interactions are more important than those between host and subhalo particles.

The relaxation time derived here is identical to that in van Kampen (1995). However, van Kampen then argued that the evaporation time due to these particle–subhalo interactions is identical to the relaxation time (i.e. no multiplication with a factor 140 is required). Van Kampen incorrectly assumes that relaxation heats the cold system to the ‘temperature’ (velocity dispersion) of the hot system. If this were true, then after a single relaxation time a very large fraction of subhalo particles would evaporate, such that the evaporation time-scale is only slightly longer than the relaxation time-scale. But this is incorrect. After all, the relaxation time is defined as the time required to heat the velocity distribution of the cold system to a Maxwellian with a velocity dispersion $\sigma_c$, rather than $\sigma_h$ (see equation 51). Hence, one still requires of the order of 140 relaxation times for the cold system to evaporate, which implies that the collisional heating due to interactions with host halo particles is less important than that due to self-interactions.

This is illustrated in Fig. 11, where the blue lines correspond to $t_{\text{evap}}^{\text{ch}}$ for a subhalo embedded in a host halo with a mass 30 times larger than that of the subhalo. Since $t_{\text{evap}}^{\text{ch}}$ is roughly proportional to $(M_h/M_c)^{1/3}$ and most subhaloes have $M_h/M_c \gg 30$, this estimate for $t_{\text{evap}}^{\text{ch}}$ may be considered a conservative lower limit. Note that indeed $t_{\text{evap}}^{\text{ch}} > t_{\text{evap}}^{\text{cc}}$, except when the softening length becomes comparable

to the size of the subhalo in question (which is not a meaningful situation to consider).

To summarize, contrary to claims by Carlberg (1994) and van Kampen (1995, 2000a), evaporation due to particle–subhalo interactions is always subdominant to that due to self-interactions among the subhalo particles. And since the latter is always longer than the Hubble time, evaporation due to collisional encounters with host halo particles cannot be an important numerical effect.

## 6.3 Impulsive heating due to encounters with host halo particles

A subtly different numerical disruption mechanism is impulsive heating of the subhalo due to encounters with the artificially massive particles of the host halo. This is different from the two-body particle–subhalo heating picture outlined above, which is a collisional process. The heating mechanism considered here, and first discussed in Moore et al. (1996b), considers the subhalo as a collisionless system being tidally heated due to the high-speed, impulsive encounters with the particles of the host halo.

Consider a particle from the hot host halo, which is modelled as a Plummer sphere of mass m and size $\epsilon$ (the softening length), that has an impulsive encounter with a (cold) subhalo of mass $M_c = N_c m$ and size $r_c$. Let b and V = V $\hat{\epsilon}_z$ be the impact parameter and relative velocity (taken to be along the z-direction) of the encounter.

As long as the impact parameter $b > r_c$, we can use the distant tide approximation, according to which a subhalo particle increases its kinetic energy as follows:

$$ \Delta v^2 = \frac{4G^2m^2}{V^2} \frac{R^2}{b^4}. \qquad (53) $$

Here, $R^2 = x^2 + y^2$ is the cylindrical radius of the particle with respect to the centre of the subhalo. For a penetrating encounter with $b = 0$, one finds that

$$ \Delta v^2 = \frac{4G^2m^2}{V^2} \frac{R^2}{(R^2 + \epsilon^2)^2}. \qquad (54) $$

(Binney & Tremaine 1987). By smoothly interpolating between these two cases, we obtain an expression for $\Delta v^2(R)$ that is valid for any impact parameter:

$$ \Delta v^2(R) = \frac{4G^2m^2}{V^2} \frac{R^2}{b^4 + (R^2 + \epsilon^2)^2}. \qquad (55) $$

By far the largest increase of kinetic energy is due to encounters with $b \sim 0$ and for particles with $R \sim \epsilon$. Using the same approach as in Section 6.2, we now integrate over impact parameter to obtain the $\Delta v^2(R)$ per crossing time $t_{\text{cross,h}}$:

$$ \begin{aligned} (\Delta v^2)_{\text{cross}}(R) &= 32 \frac{\sigma_h^2}{N_h} R^2 \int_0^{r_h} \frac{b \, db}{b^4 + (R^2 + \epsilon^2)^2} \\ &\approx 8\pi \frac{\sigma_h^2}{N_h} \frac{R^2}{R^2 + \epsilon^2}, \end{aligned} \qquad (56) $$

where we have used that $V^2 \simeq \sigma_h^2$ and that $r_h \gg r_c \ge R$. If we assume that the subhalo follows an NFW density profile with scale parameter $r_0 = r_c/c$, with c the subhalo’s concentration parameter, then the total increase in energy, per crossing time, due to impulsive shocks from high-speed encounters with the host halo particles is

$$ \begin{aligned} (\Delta E)_{\text{cross}} &= \frac{1}{2} \int_{0}^{r_c} \Sigma_c(R) (\Delta v^2)_{\text{cross}}(R) 2\pi R \, dR \\ &= 4\pi \frac{|E_c|}{N_c} \frac{\sigma_c}{\sigma_h} \frac{h_c(\epsilon/r_0)}{f_E(c)}, \end{aligned} \qquad (57) $$
where $E_c$ is the subhalo's binding energy (equation 24), and

$$h_c(y) = \frac{1}{f(c)} \int_0^c \frac{x^3 g(x) \, dx}{x^2 + y^2} \qquad (58)$$

with $g(x)$ given by equation (37).

The function $h(y)$ transits from $\mathcal{O}(1)$ for $y \ll 1$ to zero for $y \gg 1$. Hence, as long as the softening length used is small compared to the scale radius of the subhalo, one can simply set $h_c(\varepsilon/r_s) = 1$ in the above expression for $(\Delta E)_{cross}$. This is immediately evident from equation (56), which shows that every subhalo particle, independent of its cylindrical radius $R$, experiences roughly the same increase in kinetic energy per unit time. The only exception are the particles with $R < \varepsilon$, for which $(\Delta v^2)_{cross}(R) \propto R^2$. As long as $\varepsilon \ll r_0$, one therefore does not make a large error by simply writing the total increase in kinetic energy of the subhalo, per crossing time, as $(\Delta E)_{cross} = \frac{1}{2} N_c m (\Delta v^2)_{cross}(R)$. This yields equation (57) with $h_c(\varepsilon/r_s) = 1$.

If we follow Moore et al. (1996b) and van Kampen (2000a) and define the characteristic time-scale for this heating process as $t_{\text{heat}}^{\text{ch}} = (|E_c|/(\Delta E)_{cross}) t_{\text{cross}}$, then

$$t_{\text{heat}}^{\text{ch}} = \frac{2}{\pi} \ln \Lambda_h \frac{f_{\text{E}}(c)}{h_c(\varepsilon/r_0)} t_{\text{relax}}^{\text{ch}} = \zeta(c) \ln \Lambda_h t_{\text{relax}}^{\text{ch}}, \qquad (59)$$

where the second equality is only valid in the limit $\varepsilon \ll r_0$, and $\zeta(c)$ is a function that depends weakly on the concentration of the subhalo, increasing from ~0.5 for $c = 5$ to ~1.3 for $c = 50$. Hence, in a typical N-body simulation with sufficiently small softening length, the time-scale for tidal heating due to impulsive encounters with the host halo particles is roughly an order of magnitude longer than the particle–subhalo two-body relaxation time, and thus roughly an order of magnitude shorter than the corresponding evaporation time-scale.

This is illustrated in the left-hand panel of Fig. 11, where the red, solid line indicates $t_{\text{heat}}^{\text{ch}}$, computed using equation (59), once again assuming that $M_h/M_c = 30$. Note that this time-scale falls in the grey zone (i.e. becomes less than the Hubble time) for $N_c \lesssim 100$ particles. When using the softening length advocated by Power et al. (2003), $t_{\text{heat}}^{\text{ch}}$ is somewhat larger for small $N_c$, but still of order the Hubble time for $N_c \lesssim 100$. As face-value, this seems to suggest that particle–subhalo heating might actually be an important mechanism causing artificial disruption of subhaloes in numerical simulations. However, this inference is based on the assumption that $t_{\text{heat}}^{\text{ch}}$ corresponds to the actual disruption time, which, as we now argue, is incorrect.

As discussed in Sections 4.1.1 and 4.2, the problem is that the ratio between the total energy gain, $\Delta E$, and the total binding energy, $|E_c|$, is of little use to gauge whether or not the system will disrupt. It is important to account for how the energy $\Delta E$ is distributed over the particles. Contrary to the case of tidal heating by a single encounter, where the energy gain of particle $i$ is proportional to the square of the particle's distance from the centre of the subhalo, $(\Delta E)_i \propto r_i^2$, in the case considered here, every subhalo particle roughly receives the same amount of energy (see equation 60); the only exception are the particles with $R < \varepsilon$, for which $(\Delta E)_i \propto R_i^2$.

In order to assess the impact of $(\Delta E)_{cross}$ on the system as a whole, we follow the approach of Section 4.1.1 and compare the individual $(\Delta E)_i$ to the individual binding energies, $|E_i|$, and define the disruption time as the number of crossing times required so that $(\Delta E)_i > |E_i|$ for the majority of particles (i.e. 95 per cent). To illustrate how this modifies the inferred heating time-scale, we proceed as follows. We construct an N-body realization of an isotropic NFW subhalo with concentration parameter $c = 10$, using the method

described in Section 3. For each individual particle, we compute the ratio $|E_i|/(\Delta E)_i$. Here, $|E_i| = v_i^2/2 + \Phi(r_i)$ is the specific binding energy of particle $i$, and

$$(\Delta E)_i = \frac{1}{2} (\Delta v^2)_i \approx 4\pi \frac{\sigma_h^2}{N_h} \frac{R^2}{R^2 + \varepsilon^2} \qquad (60)$$

is the amount of specific kinetic energy added per crossing time due to impulsive encounters with the host particles (cf. equation 56). Note that we have assumed that the $\vec{v}_i \cdot \Delta \vec{v}_i$ term (cf. equation 26), summed over all encounters, adds to zero. This is motivated by the fact that the particle is moving in between encounters. Hence, integrated over a crossing time, roughly every angle between $\vec{v}_i$ and $\Delta \vec{v}_i$ is equally represented.

For each particle, we interpret $|E_i|/(\Delta E)_i$ as the number of crossing times required to unbind the particle. This ignores the fact that the system will continuously re-virialize in order to adjust to the newly acquired (impulsively added) energy and the resulting mass-loss, and what follows is therefore only a rough estimate. We define the 'corrected' time-scale for disruption due to particle–subhalo tidal heating as the number of crossing times required for 95 per cent of all subhalo particles to become unbound; that is, we find the 95 percentile of the distribution of $|E_i|/(\Delta E)_i$ and interpret that as $t_{\text{heat,corr}}^{\text{ch}}/t_{\text{cross}}$. The $t_{\text{heat,corr}}^{\text{ch}}$ thus obtained is shown as dashed, red lines in Fig. 11. As long as $\varepsilon \ll r_0$, this correction boosts $t_{\text{heat}}^{\text{ch}}$ by roughly an order of magnitude, making it very similar to the two-body evaporation time, $t_{\text{evap}}^{\text{cc}}$. When the softening length becomes comparable to the size of the subhalo, $t_{\text{heat,corr}}^{\text{ch}}$ becomes independent of $N_c$, and many orders of magnitude larger than the Hubble time.

We conclude that similar to two-body relaxation, artificial tidal heating due to impulsive encounters with the overly massive particles of the host halo cannot be responsible for artificial disruption of subhaloes in numerical simulations.

## 6.4 Force softening

Despite many decades of simulation work, and numerous previous studies, there is still no consensus as to the optimal choice for the softening parameter in N-body simulations. Numerous studies have argued that minimizing the scatter in certain statistical descriptions of the matter field (i.e. power spectrum or two-point correlation function) among an ensemble of N-body realizations requires a softening length $\varepsilon$ that is at least as large as the mean inter-particle separation $l = L_{\text{box}}/N^{1/3}$ (e.g. Melott et al. 1997; Splinter et al. 1998; Romeo et al. 2008; Joyce, Marcos & Baertschiger 2009). Here, $L_{\text{box}}$ is the (periodic) box size of the cosmological simulation and N is the total number of particles used. Most cosmological simulations, though, are run with much smaller softening lengths, typically of the order of $\varepsilon/l \sim 0.01-0.02$. For example, the Millennium simulation (Springel et al. 2005) has $\varepsilon/l \sim 0.022$, while the MDPL and SMDPL simulations of Klypin et al. (2016) adopt $\varepsilon/l \sim 0.019$ and $\varepsilon/l \sim 0.014$, respectively. Whether such simulations produce reliable results on scales below the (comoving) mean inter-particle separation is still a matter of fierce, ongoing debate (see e.g. Knebe et al. 2000; Melott 2007; Joyce et al. 2009; Power et al. 2016; Benhaiem, Joyce & Sylos Labini 2017, and references therein).

Several studies have focused on the optimal softening length for simulating individual objects (i.e. individual DM haloes or galaxies). For example, Dehnen (2001), extending studies by Merritt (1996) and Athanassoula et al. (2000), suggests that one should aim to minimize force errors. This can be achieved by optimizing the softening length such that it simultaneously minimizes the force bias (which increases with increasing $\varepsilon$) and the force noise (which
decreases with increasing $\epsilon$). Based on such considerations, Dehnen (2001) advocates that a collisionless system with a Hernquist (1990) density profile (which has the same $r^{-1}$-cusp as a NFW profile) should be simulated with a Plummer softening $\epsilon/r_0 = 0.017 N_5^{-0.23}$, where $r_0$ is the scale radius of the Hernquist (or NFW) profile and $N_5 = N/10^5$ with $N$ the number of particles used to model the system in question. Power et al. (2003) use a different criterion to optimize the softening length: by minimizing the impact of two-body effects on the (central) density profiles of virialized DM haloes in cosmological simulations, they advocate a much larger softening length of $\epsilon/r_0 = 0.126 (c/10) N_5^{-0.5}$, where $c$ is the NFW halo concentration parameter. Note that both of these criteria have the optimal softening length depend on the size of the system as well as the number of particles that is used to represent the system. If we adopt a concentration-mass relation of the form $c \propto M^{-0.1}$ (e.g. Dutton & Macciò 2014) and use that the virial radii of DM haloes scale according to $r_{vir} \propto M^{1/3} \propto N^{1/3}$, we find that the optimal softening length depends on halo mass according to $\epsilon_{opt} \propto M^{\alpha}$, with $\alpha = +0.2$ and $-0.17$ for the criteria of Dehnen and Power et al., respectively. This scaling is sufficiently weak so that the softening is close to optimal for haloes spanning a relatively large range in halo masses, as required for cosmological simulations.

Of relevance for the discussion in this paper is how softening impacts the dynamical evolution of DM substructure. To our knowledge, there has been no concerted effort to investigate this. We fill this gap in Paper II, with a detailed study of how softening impacts the evolution and survivability of subhaloes. We demonstrate that whenever the tidal field is strong (i.e. close to the centre of the host halo), the survivability of substructure is extremely sensitive to the softening used. Furthermore, as a subhalo experiences mass-loss, both its size and its number of particles change with time, which in turn causes an evolution in the optimal softening length. As we demonstrate in Paper II, typically $\epsilon_{opt}$ drops by an order of magnitude after the first peri-centric passage, after which it remains fairly unchanged. Since simulations do not account for this drastic drop in the optimal softening length, substructure is typically evolved with a softening length that is too large, which results in significant artificial subhalo disruption.<fnref n="9" /> Based on the discussion in this section, and the results presented in Paper II, we therefore conclude that inadequate softening is the dominant driver of subhalo disruption in numerical simulations.

## 7 SUMMARY AND DISCUSSION

Being able to accurately predict the abundance and demographics of DM substructure is of paramount importance for many fields of astrophysics: gravitational lensing, galaxy evolution, halo occupation modelling, and even constraining the nature of the DM. DM substructure is subject to tidal stripping and heating, which are highly non-linear processes and therefore best studied using numerical $N$-body simulations. These reveal prevalent subhalo disruption, with only ~35 per cent of subhaloes accreted at $z=1$ surviving until $z=0$ (Jiang & van den Bosch 2016). Based on a hand-full of convergence studies, it is generally assumed that most of the disruption of subhaloes with more than 50–100 particles is physical, rather than numerical. However, there is no consensus in the literature as to what causes subhaloes to disrupt. In addition, convergence is

only a necessary, but not a sufficient condition, to guarantee that the simulation results are physically correct.

All this warrants a detailed, comprehensive investigation into the disruption of substructure in numerical simulations. This is the first paper in a series in which we use both idealized numerical simulations and (semi)-analytical treatments to study under what conditions DM substructure undergoes physical and/or numerical disruption. In this paper, we have focused mainly on analytical treatments of tidal heating and tidal stripping. We also presented a crude analytical treatment of potential numerical disruption mechanisms. Our main conclusions are as follows:

(i) During their first orbit after accretion, subhaloes pass the centre of their host halo at a median peri-centric distance of ~0.37$R_{vir,h}$, while ~7 per cent come within 0.1$R_{vir,h}$.

(ii) Contrary to naive expectation, subhaloes that experience a tidal shock $\Delta E$ that exceeds the subhalo's binding energy, $|E_b|$, do not necessarily undergo disruption. Rather, they experience mass-loss. We have presented a fitting function (equation 27) that relates the amount of mass stripped to the ratio $\Delta E/|E_b|$, which is valid for subhaloes with NFW density distributions. An NFW subhalo for which $\Delta E/|E_b| = 1 (100)$ is stripped of only ~20 (80) per cent of its mass.

(iii) During first peri-centric passage, tidal heating increases the internal energy of the subhalo by a median $\Delta E \sim 2|E_b|$. This results in the subhalo losing ~30 per cent of its mass. In the rare case of a purely radial orbit, the subhalo experiences a tidal shock that may strip as much as 90 per cent of its mass.

(iv) The impulse approximation, combined with the instantaneous mass-loss approximation of Aguilar & White (1985), accurately predicts the amount of mass-loss experienced by subhaloes on fairly radial orbits ($\eta \lesssim 0.2$).

(v) Tidal heating due to high-speed (impulsive) encounters with other subhaloes, sometimes called 'harassment', is negligible compared to the tidal heating associated with the peri-centric passage of the host halo.

(vi) Modelling tidal stripping in the non-impulsive regime is extremely difficult, mainly because the concept of tidal radius is ill defined and because we lack a theory to describe how a gravitational $N$-body system re-virializes as it undergoes stripping. Depending on which definition one adopts for the tidal radius, and what one assumes regarding the *rate* at which material outside of the (instantaneous) tidal radius is stripped, one can obtain dramatically different results, none of which adequately reproduce simulation results.

(vii) Instantaneously, stripping matter from the outskirts of an NFW halo can leave a remnant with positive binding energy (cf. Hayashi et al. 2003). The remnant’s binding energy depends strongly on the orbital anisotropy of the subhalo, with more radially anisotropic subhaloes having larger (i.e. more positive or less negative) binding energies for a given amount of stripping.

(viii) Instantaneously stripping matter from an NFW halo never leads to subhalo disruption, even if the remnant has positive binding energy.

(ix) Contrary to claims in the literature, two-body relaxation of subhaloes due to collisions between subhalo particles and (the more abundant) particles from the host halo is less efficient than that due to collisions among subhalo particles themselves. Neither is efficient enough, however, to cause significant evaporation of subhaloes in numerical simulations.

(x) Tidal heating of subhaloes due to high-speed, impulsive encounters with host halo particles that are artificially massive (due to

<footnote n="9">⁹ The simulations used in this paper adopted a softening that was carefully calibrated using the detailed convergence studies described in Paper II, and therefore do not suffer from this effect.</footnote>
the limiting mass resolution of the simulation) cannot cause artificial disruption of substructure.

Based on these results, we conclude that it is extremely difficult to physically disrupt CDM subhaloes (in the absence of baryonic effects; see below). Tidal stripping alone never leads to complete disruption, and tidal heating is not effective in the central regions of subhaloes, which are adiabatically ‘shielded’ (Spitzer 1987; Weinberg 1994a,b). In principle, multiple tidal shocks, associated with multiple peri-centric passages, could ultimately cause the subhalo to disrupt. We did not explore such a scenario in this paper, since we lack a reliable analytical treatment for how subhaloes re-virialize after tidally induced mass-loss. However, we do investigate the impact of multiple peri-centric passages in Papers II and III, in which we use hundreds of idealized numerical simulations. As we demonstrate in those papers, physical disruption in CDM-based DM-only simulations is extremely rare (see also Diemand et al. 2007; Peñarrubia et al. 2010). On the other hand, we do find that the simulation results are extremely sensitive to the choice of softening length. In particular, we demonstrate that state-of-the-art cosmological simulations still suffer from serious overmerging, largely as a consequence of inappropriate force softening. In addition, we demonstrate that N-body simulations suffer from a discreteness-driven run-away instability, which makes it impossible to reliably trace the evolution of subhaloes in a strong tidal field unless the subhalo has at least 1000 particles.

Numerical overmerging is a serious road-block for the many astrophysical applications that require accurate characterization of halo substructure. In the coming decade, a number of large galaxy surveys, such as the Dark Energy Spectroscopic Instrument (DESI), the Large Synoptic Survey Telescope (LSST), Euclid, and the Wide Field InfraRed Space Telescope (WFIRST), will provide astrophysicists with an unprecedented wealth of data regarding galaxy evolution and cosmology. In order to optimize the scientific impact of these huge investments, it is prudent that galaxy clustering in these surveys be interpreted in an unbiased, maximally informative manner. A popular method to do so is to compare the observations to mock data, obtained by populating DM (sub)haloes in DM-only simulations with ‘mock’ galaxies, using methods such as subhalo abundance matching (see Section 1 for references). This entire programme faces a severe challenge if those simulations underpredict the abundance of subhaloes due to artificial disruption. In principle, one can correct for this by including ‘orphans’, i.e. mock galaxies without an associated subhalo in the simulation (e.g. Kitzbichler & White 2008; Moster et al. 2013, 2017; Guo & White 2014); however, unless it is known how many orphans to add, and where, this seriously diminishes the information content of small-scale clustering (see Campbell et al. 2017, for a detailed discussion). DM substructure is also an important discriminator between different DM models (cold versus warm versus self-interacting), and is important for translating a potential, observed DM annihilation signal into a mass and annihilation cross-section of the DM particle. Unless we can make accurate, and above all reliable, predictions regarding the abundance and structure of DM subhaloes, we will forfeit one of the main handles we have on learning about the nature of DM.

We end by emphasizing that this work has entirely focused on DM only, without taking account of potential baryonic effects. There is a rapidly expanding literature on how baryons may impact the abundance and demographics of DM substructure (e.g. Macciò et al. 2006; Weinberg et al. 2008; Dolag et al. 2009; Arraki et al. 2014; Brooks & Zolotov 2014; Fiacconi et al. 2016; Wetzel et al. 2016; Despali & Vegetti 2017; Garrison-Kimmel et al. 2017).

Our work, in no way, aims to undermine the importance of baryons. However, before we can make reliable predictions for how baryonic processes modify substructure in the dark sector, we first need to establish a better understanding of how tides impact substructure, and develop tools and criteria to assess the reliability of simulations (be it hydro or N-body). Hence, this work is to be considered a necessary first step towards a more reliable treatment of DM substructure (including satellite galaxies), and not the final answer.

## ACKNOWLEDGEMENTS

We are grateful to the referee, Chris Power, for an insightful referee report, and to Andrew Hearin, Fangzhou Jiang, and Fred Rasio for useful discussion. FvdB is supported by the Klaus Tschira Foundation and by the US National Science Foundation through grant AST 1516962 and is grateful to the Munich Excellence Cluster for its hospitality. Part of this work was performed at the Aspen Center for Physics, which is supported by the National Science Foundation under grant PHY-1066293, and at the Kavli Institute for Theoretical Physics, which is supported by the National Science Foundation under grant PHY-1125915. OH and GO were supported by funding form the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 679145, project ‘COSMO-SIMS’).

## REFERENCES

* Aguilar L. A., White S. D. M., 1985, ApJ, 295, 374

* Arraki K. S., Klypin A., More S., Trujillo-Gomez S., 2014, MNRAS, 438, 1466

* Athanassoula E., Fady E., Lambert J. C., Bosma A., 2000, MNRAS, 314, 475

* Barnes J., Hut P., 1986, Nature, 324, 446

* Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013a, ApJ, 762, 109

* Behroozi P. S., Wechsler R. H., Conroy C., 2013b, ApJ, 770, 57

* Benhaiem D., Joyce M., Sylos Labini F., 2017, MNRAS, 470, 4099

* Bergström L., Edsjö J., Gondolo P., Ullio P., 1999, Phys. Rev. D, 59, 043506

* Binney J., Tremaine S., 1987, Galactic Dynamics. Princeton Univ. Press, NJ

* Binney J., Tremaine S., 2008, Galactic Dynamics, 2nd edn. Princeton Univ. Press, NJ

* Bose S. et al., 2017, MNRAS, 464, 4520

* Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., 2010, MNRAS, 406, 896

* Boylan-Kolchin M., Bullock J. S., Kaplinghat M., 2011, MNRAS, 415, L40

* Bradač M., Schneider P., Steinmetz M., Lombardi M., King L. J., Porcas R., 2002, A&A, 388, 373

* Brainerd T. G., Goldberg D. M., Verner Villumsen J., 1998, ApJ, 502, 505

* Brooks A. M., Zolotov A., 2014, ApJ, 786, 87

* Bryan G. L., Norman M. L., 1998, ApJ, 495, 80

* Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343

* Burkert A., 2000, ApJ, 534, L143

* Campbell D., van den Bosch F. C., Padmanabhan N., Mao Y.-Y., Zentner A. R., Lange J. U., Jiang F., Villarreal A., 2017, MNRAS, preprint (arXiv:1705.06347)

* Carlberg R. G., 1994, ApJ, 433, 468

* Carlberg R. G., 2009, ApJ, 705, L223

* Cautun M., Hellwing W. A., van de Weygaert R., Frenk C. S., Jones B. J. T., Sawala T., 2014, MNRAS, 445, 1820

* Chandrasekhar S., 1943, ApJ, 97, 255

* Colín P., Avila-Reese V., González-Samaniego A., Velázquez H., 2015, ApJ, 803, 28

* Conroy C., Wechsler R. H., Kravtsov A. V., 2006, ApJ, 647, 201

* Dalal N., Kochanek C. S., 2002, ApJ, 572, 25

* Dehnen W., 2001, MNRAS, 324, 273

* Dekel A., Devor J., Hetzroni G., 2003, MNRAS, 341, 326
Despali G., Vegetti S., 2017, MNRAS, 469, 1997

Diemand J., Moore B., Stadel J., Kazantzidis S., 2004a, MNRAS, 348, 977

Diemand J., Moore B., Stadel J., 2004b, MNRAS, 352, 535

Diemand J., Kuhlen M., Madau P., 2007, ApJ, 667, 859

Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497

Drakos N. E., Taylor J. E., Benson A. J., 2017, MNRAS, 468, 2345

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

Dutton A. A. et al., 2016, MNRAS, 461, 2658

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

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

Faltenbacher A., Diemand J., 2006, MNRAS, 369, 1698

Farouki R. T., Salpeter E. E., 1982, ApJ, 253, 512

Fiacconi D., Madau P., Potter D., Stadel J., 2016, ApJ, 824, 144

Frenk C. S., White S. D. M., Davis M., Efstathiou G., 1988, ApJ, 327, 507

Gao L., White S. D. M., Jenkins A., Stoehr F., Springel V., 2004, MNRAS, 355, 819

Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Kirby E. N., 2014, MNRAS, 444, 222

Garrison-Kimmel S. et al., 2017, MNRAS, 471, 1709

Ghigna S., Moore B., Governato F., Lake G., Quinn T., Stadel J., 1998, MNRAS, 300, 146

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

Giocoli C., Pieri L., Tormen G., 2008b, MNRAS, 387, 689

Giocoli C., Tormen G., Sheth R. K., van den Bosch F. C., 2010, MNRAS, 404, 502

Gnedin O. Y., Ostriker J. P., 1999, ApJ, 513, 626

Gnedin O. Y., Hernquist L., Ostriker J. P., 1999, ApJ, 514, 109 (GHO99)

Gonzalez-Casado G., Mamon G. A., Salvador-Sole E., 1994, ApJ, 433, L61

Griffen B. F., Ji A. P., Dooley G. A., Gómez F. A., Vogelsberger M., O'Shea B. W., Frebel A., 2016, ApJ, 818, 10

Guo Q., White S., 2014, MNRAS, 437, 3228

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

Hahn O., Angulo R. E., 2016, MNRAS, 455, 1115

Han J., Jing Y. P., Wang H., Wang W., 2012, MNRAS, 427, 2437

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

Hayashi E., Navarro J. F., Taylor J. E., Stadel J., Quinn T., 2003, ApJ, 584, 541

Hearin A. P., Zentner A. R., Berlind A. A., Newman J. A., 2013, MNRAS, 433, 659

Hearin A. P., Watson D. F., van den Bosch F. C., 2015, MNRAS, 452, 1958

Hearin A. P., Zentner A. R., van den Bosch F. C., Campbell D., Tollerud E., 2016, MNRAS, 460, 2552

Hernquist L., 1990, ApJ, 356, 359

Hernquist L., Ostriker J. P., 1992, ApJ, 386, 375

Jiang F., van den Bosch F. C., 2016, MNRAS, 458, 2848

Jiang F., van den Bosch F. C., 2017, MNRAS, 472, 657

Jiang L., Cole S., Sawala T., Frenk C. S., 2015, MNRAS, 448, 1674

Joyce M., Marcos B., Baertschiger T., 2009, MNRAS, 394, 751

Kampakoglou M., Benson A. J., 2007, MNRAS, 374, 775

Kang X., Jing Y. P., Mo H. J., Börner G., 2005, ApJ, 631, 21

Keeton C. R., Moustakas L. A., 2009, ApJ, 699, 1720

King I., 1962, AJ, 67, 471

Kitzbichler M. G., White S. D. M., 2008, MNRAS, 391, 1489

Klypin A., Gottlöber S., Kravtsov A. V., Khokhlov A. M., 1999a, ApJ, 516, 530

Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999b, ApJ, 522, 82

Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102

Klypin A., Prada F., Yepes G., Heß S., Gottlöber S., 2015, MNRAS, 447, 3693

Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340

Knebe A., Kravtsov A. V., Gottlöber S., Klypin A. A., 2000, MNRAS, 317, 630

Knebe A., Power C., Gill S. P. D., Gibson B. K., 2006, MNRAS, 368, 741

Knebe A., Arnold B., Power C., Gibson B. K., 2008, MNRAS, 386, 1029

Knebe A. et al., 2013, MNRAS, 435, 1618

Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, ApJS, 111, 73

Kravtsov A. V., Berlind A. A., Wechsler R. H., Klypin A. A., Gottlöber S., Allgood B., Primack J. R., 2004, ApJ, 609, 35

Kuhlen M., Diemand J., Madau P., 2008, ApJ, 686, 262

Kuijken K., Dubinski J., 1994, MNRAS, 269, 13

Lehmann B. V., Mao Y.-Y., Becker M. R., Skillman S. W., Wechsler R. H., 2017, ApJ, 834, 37

Łokas E. L., 2009, MNRAS, 394, L102

Lovell M. R., Frenk C. S., Eke V. R., Jenkins A., Gao L., Theuns T., 2014, MNRAS, 439, 300

Macciò A. V., Moore B., Stadel J., Diemand J., 2006, MNRAS, 366, 1529

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

Melott A. L., 2007, preprint (arXiv:0709.0745)

Melott A. L., Shandarin S. F., Splinter R. J., Suto Y., 1997, ApJ, 479, L79

Merritt D., 1996, AJ, 111, 2462

Metcalf R. B., Madau P., 2001, ApJ, 563, 9

Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319

Mo H., van den Bosch F. C., White S., 2010, Galaxy Formation and Evolution. Cambridge Univ. Press, Cambridge

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

Moore B., Katz N., Lake G., Dressler A., Oemler A., 1996a, Nature, 379, 613

Moore B., Katz N., Lake G., 1996b, ApJ, 457, 455

Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, ApJ, 499, L5

Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19

Moster B. P., Somerville R. S., Maulbetsch C., van den Bosch F. C., Macciò A. V., Naab T., Oser L., 2010, ApJ, 710, 903

Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121

Moster B. P., Naab T., White S. D. M., 2017, MNRAS, preprint (arXiv:1705.05373)

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

Navarro J. F. et al., 2010, MNRAS, 402, 21

Ogiya G., Burkert A., 2015, MNRAS, 446, 2363

Ogiya G., Burkert A., 2016, MNRAS, 457, 2164

Ogiya G., Mori M., 2014, ApJ, 793, 46

Ogiya G., Mori M., Miki Y., Boku T., Nakasato N., 2013, J. Phys.: Conf. Ser., 454, 012014

Oguri M., Lee J., 2004, MNRAS, 355, 120

Onions J. et al., 2012, MNRAS, 423, 1200

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

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

Pieri L., Bertone G., Branchini E., 2008, MNRAS, 384, 1627

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

Power C., Robotham A. S. G., Obreschkow D., Hobbs A., Lewis G. F., 2016, MNRAS, 462, 474

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in FORTRAN. The Art of Scientific Computing. Cambridge Univ. Press, Cambridge

Pullen A. R., Benson A. J., Moustakas L. A., 2014, ApJ, 792, 24

Read J. I., Wilkinson M. I., Evans N. W., Gilmore G., Kleyna J. T., 2006, MNRAS, 366, 429

Reddick R. M., Tinker J. L., Wechsler R. H., Lu Y., 2014, ApJ, 783, 118

Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Óñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81

Romeo A. B., Agertz O., Moore B., Stadel J., 2008, ApJ, 686, 1

Sanders J. L., Bovy J., Erkal D., 2016, MNRAS, 457, 3817

Spitzer L. Jr, 1958, ApJ, 127, 17

Spitzer L., 1987, Dynamical Evolution of Globular Clusters. Princeton Univ. Press, NJ

Splinter R. J., Melott A. L., Shandarin S. F., Suto Y., 1998, ApJ, 497, 38

Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, MNRAS, 328, 726

Springel V. et al., 2005, Nature, 435, 629

Springel V. et al., 2008, MNRAS, 391, 1685
Strigari L. E., Koushiappas S. M., Bullock J. S., Kaplinghat M., 2007, Phys. Rev. D, 75, 083526

Taffoni G., Mayer L., Colpi M., Governato F., 2003, MNRAS, 341, 434

Taylor J. E., Babul A., 2001, ApJ, 559, 716

Taylor J. E., Babul A., 2004, MNRAS, 348, 811

Tollet É., Cattaneo A., Mamon G. A., Moutard T., van den Bosch F. C., 2017, MNRAS, 471, 4170

Tormen G., 1997, MNRAS, 290, 411

Tormen G., Bouchet F. R., White S. D. M., 1997, MNRAS, 286, 865

Tormen G., Diaferio A., Syer D., 1998, MNRAS, 299, 728

Trujillo-Gomez S., Klypin A., Primack J., Romanowsky A. J., 2011, ApJ, 742, 16

van den Bosch F. C., 2017, MNRAS, 468, 885 (vdB17)

van den Bosch F. C., Jiang F., 2016, MNRAS, 458, 2870

van den Bosch F. C., Norberg P., Mo H. J., Yang X., 2004, MNRAS, 352, 1302

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

van den Bosch F. C., Jiang F., Campbell D., Behroozi P., 2016, MNRAS, 455, 158

van den Bosch F. C., Ogiya G., 2018, MNRAS, submitted

van Kampen E., 1995, MNRAS, 273, 295

van Kampen E., 2000a, preprint (arXiv:astro-ph/0002027)

van Kampen E., 2000b, preprint (arXiv:astro-ph/0008453)

Vale A., Ostriker J. P., 2004, MNRAS, 353, 189

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

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

Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740

Walker M. G., Mateo M., Olszewski E. W., Peñarrubia J., Wyn Evans N., Gilmore G., 2009, ApJ, 704, 1274

Weinberg M. D., 1993, ApJ, 410, 543

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

Weinberg M. D., 1994b, AJ, 108, 1403

Weinberg M. D., 1997, ApJ, 478, 435

Weinberg D. H., Colombi S., Davé R., Katz N., 2008, ApJ, 678, 6

Wetzel A. R., 2011, MNRAS, 412, 49

Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, ApJ, 827, L23

Widrow L. M., 2000, ApJS, 131, 39

Wu H.-Y., Hahn O., Evrard A. E., Wechsler R. H., Dolag K., 2013, MNRAS, 436, 460

Zentner A. R., Bullock J. S., 2003, ApJ, 598, 49

Zentner A. R., Berlind A. A., Bullock J. S., Kravtsov A. V., Wechsler R. H., 2005, ApJ, 624, 505

Zentner A. R., Hearin A. P., van den Bosch F. C., 2014, MNRAS, 443, 3044

Zentner A. R., Hearin A., van den Bosch F. C., Lange J. U., Villarreal A., 2016, MNRAS, preprint (arXiv:1606.07817)

Zhao H., 2004, MNRAS, 351, 891

Zolotov A. et al., 2012, ApJ, 761, 71

## APPENDIX A: COMPUTING BOUND FRACTIONS OF N-BODY HALOES

This appendix describes the method we use to compute the bound fraction of an N-body (sub)halo. We consider a particle *i* to be bound if

$$E_i \equiv \frac{1}{2} m_i v_i^2 - \sum_{j \neq i} \frac{G m_i m_j Q_j}{(|\mathbf{r}_j - \mathbf{r}_i|^2 + \varepsilon^2)^{1/2}} < 0. \quad (A1)$$

Here, $\varepsilon$ is the softening length and $Q_j$ is equal to 1 (0) if particle $j$ is bound (unbound). Hence, in order to be able to compute $E_i$, one first needs to know which are the bound particles. Obviously, this can only be solved using an iterative scheme. In our analysis, we proceed as follows:

(i) We begin by making an initial guess for the boundness, $Q_i$, of each particle. At $t=0$ (ICs), we assume that each particle is bound

$(Q_i = 1 \forall i = 1, ..., N)$, while at later times we assume that $Q_i$ is the same as in the previous simulation output.

(ii) For each particle, we compute $E_i$ using equation (A1). We update $Q_i$ accordingly and compute the new, bound fraction $f_{\text{bound}} = \frac{1}{N_{\text{tot}}} \sum_i Q_i$.

(iii) Compute the centre-of-mass position and velocity of the halo, $\mathbf{r}_{\text{com}}$ and $\mathbf{v}_{\text{com}}$, as the average position and velocity of the $N = \text{MAX}(N_{\text{min}}, f_{\text{com}} N_{\text{bound}})$ particles that are most bound. Here, $N_{\text{bound}} = f_{\text{bound}} N_{\text{tot}}$, while $N_{\text{min}}$ and $f_{\text{com}}$ are free parameters.

(iv) Update $\mathbf{v}_i$ for each particle using the new centre-of-mass velocity.

(v) Go back to (ii) and iterate until the changes in $\mathbf{r}_{\text{com}}$ and $\mathbf{v}_{\text{com}}$ are smaller than $10^{-4} r_{\text{vir}}$ and $10^{-4} V_{\text{vir}}$, respectively. This typically requires 3–10 iterations.

When computing the gravitational potential term of equation (A1), we use a Barnes & Hut octree with the same opening angle, $\theta$, as in the simulation. We also adopt the same softening. Generally, smaller values of $f_{\text{com}}$ results in a more noisy time evolution of $\mathbf{r}_{\text{com}}$ and $\mathbf{v}_{\text{com}}$, while for larger values of $f_{\text{com}}$ it is more difficult to trace $f_{\text{bound}}(t)$ when it becomes small. We obtain stable results for $N_{\text{min}} = 10$ (we never use fewer than 10 particles to determine the centre-of-mass properties$^{10}$) and $f_{\text{com}} = 0.05$ (centre-of-mass properties are determined using the 5 per cent most bound particles), which are the parameters we use throughout.

## APPENDIX B: THE ENERGY STRUCTURE OF DM HALOES

Consider a spherical DM halo with density distribution $\rho(r)$ and enclosed mass profile

$$M(r) = 4\pi \int_{0}^{r} \rho(r') r'^2 dr'. \quad (B1)$$

The total energy of the DM particles inside radius $r$ is given by $E(r) = K(r) + W(r)$, where $K$ and $W$ are the corresponding kinetic and potential energies, respectively. For a spherical system,

$$W(r) = 2\pi \int_{0}^{r} \rho(r') \Phi(r') r'^2 dr'. \quad (B2)$$

Using that the gravitational potential

$$\Phi(r) = -4\pi G \left[ \frac{1}{r} \int_{0}^{r} \rho(r') r'^2 dr' + \int_{r}^{\infty} \rho(r') r' dr' \right] \quad (B3)$$

(Binney & Tremaine 2008), and using integration by parts, this can be written as

$$W(r) = -4\pi G \int_{0}^{r} \rho(r') M(r') r' dr' - 2\pi G M(r) \int_{r}^{\infty} \rho(r') r' dr'. \quad (B4)$$

The second term describes the contribution to the gravitational potential energy due to matter located at radii larger than $r$ and vanishes if $r \to \infty$. In what follows, we compute the energy of a spherical halo whose mass is instantaneously truncated (e.g. due to tides) at radius $r_i$. In that case, $\rho(r) = 0$ for $r > r_i$ and $W(r_i)$ is simply given by the first term of equation (B4).

<footnote>$^{10}$ Except when $N_{\text{bound}} < 10$, in which case we adopt $N = N_{\text{bound}}$.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure A1.** Energy structure of spherical, instantaneously truncated NFW haloes. The left-hand panel plots the total binding energy, $E(r)$, as a function of the truncation radius, $r$, in units of the NFW scale radius, $r_0$. The middle panel shows the corresponding virial ratio, $K(r)/|W(r)|$. Different colours correspond to different values of the anisotropy parameter; from blue to red $\beta = 1.0, 0.8, ..., -3.0$. The panel on the right shows $\log[r_{\text{crit}}/r_0]$ as function of $\beta$, where $r_{\text{crit}}$ is defined as the truncation radius below which the remnant has positive binding energy. $r_{\text{crit}} \simeq 0.77r_0$ for an isotropic NFW halo ($\beta = 0$), becomes $\sim 2.84r_0$ for an NFW halo with purely radial orbits ($\beta = 1$), and is absent (i.e. the remnant always has negative binding energy) whenever $\beta \lesssim -1.5$.</figcaption></figure>

The Jeans equation for hydrostatic equilibrium in a spherical stellar system is

$$ \frac{d(\rho\sigma_r^2)}{dr} + 2\frac{\beta\rho\sigma_r^2}{r} = -\rho\frac{d\Phi}{dr}, \qquad (B5) $$

where $\beta \equiv 1 - \sigma_\theta^2/\sigma_r^2$ describes the velocity dispersion anisotropy (here assumed constant). Haloes with $\beta = 0$ are isotropic, haloes with $\beta < 0$ are tangentially anisotropic, and haloes with $0 < \beta \le 1$ are radially anisotropic. The solution of equation (B5) is

$$ \sigma_r^2(r) = \frac{r^{-2\beta}}{\rho(r)} \int_r^\infty \rho(r') M(r') r'^{2\beta-2} dr', \qquad (B6) $$

where we have used that $d\Phi/dr = GM(r)/r^2$. The kinetic energy of particles inside radius $r_t$ is given by

$$ K(r_t) = 2\pi \int_0^{r_t} \rho(r) \langle v^2 \rangle(r) r^2 dr. \qquad (B7) $$

If we assume that the halo has no net streaming motion, so that $\sigma_\phi^2 = \sigma_\theta^2$, then $\langle v^2 \rangle(r) = [1 + 2(1 - \beta)]\sigma_r^2(r)$ and we have that

$$ K(r_t) = 2\pi G [1 + 2(1 - \beta)] \int_0^{r_t} dr \, r^{2-2\beta} \\ \times \int_r^\infty \rho(r') M(r') r'^{2\beta-2} dr'. \qquad (B8) $$

The left-hand panel of Fig. A1 plots the total binding energy $E(r_t) = K(r_t) + W(r_t)$ of a spherical NFW halo that is instantaneously

truncated at radius $r_t$. Different lines correspond to different values of the anisotropy parameter $\beta$, as described in the figure caption. Note that these results are independent of halo concentration. The middle panel shows the corresponding virial ratio $K(r_t)/|W(r_t)|$, which asymptotes to 0.5 (corresponding to virial equilibrium) in the limit $r_t \to \infty$. Note that $K(r_t)/|W(r_t)| > 1$ corresponds to $E(r_t) > 0$, and thus a positive, total binding energy. Finally, the right-hand panel plots $\log[r_{\text{crit}}/r_0]$ as function of $\beta$, where $r_{\text{crit}}$ is defined such that $E(r_t) > 0$ for $r_t < r_{\text{crit}}$. For an isotropic NFW profile ($\beta = 0$), we have that $r_{\text{crit}} = 0.77r_0$ in agreement with Hayashi et al. (2003). For a maximally, radially anisotropic NFW halo, $r_{\text{crit}} = 2.84r_0$, while $r_{\text{crit}} = 0$ if $\beta \lesssim -1.6$. In this case, the tidally truncated remnant will always have negative, total binding energy, independent of $r_t$. To see this, recall that the limit $\beta \to -\infty$ corresponds to a halo in which all particles are on circular orbits, so that $\langle v^2 \rangle(r) = v_{\text{circ}}^2(r) = GM(r)/r$. In that case, equation (B7) reduces to

$$ K(r_t) = 2\pi G \int_0^{r_t} \rho(r) M(r) r dr = -\frac{W(r_t)}{2}. \qquad (B9) $$

Hence, we have that $K(r_t)/|W(r_t)| = \frac{1}{2}$, and thus $E(r_t) = -K(r_t) < 0$.

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