DASH: a library of dynamical subhalo evolution

TL;DR

This paper introduces a library of high-resolution N-body simulations for the tidal evolution of dark matter subhaloes, providing a tool to calibrate semi-analytical models and address numerical artefacts in cosmological simulations.

Contribution

It presents a comprehensive library of idealized subhalo evolution simulations and a pre-trained machine learning model for efficient interpolation of subhalo mass evolution.

Findings

Simulation library covers realistic orbital and concentration parameters.

Pre-trained model interpolates simulation data at 0.1 dex accuracy.

Library aids in improving semi-analytical models of subhalo evolution.

Abstract

The abundance and demographics of dark matter substructure is important for many areas in astrophysics and cosmological $N$-body simulations have been the primary tool used to investigate them. However, it has recently become clear that the simulations are subject to numerical artefacts, which hampers a proper treatment of the tidal evolution of subhaloes. Unfortunately, no analytical models that accurately describe subhalo evolution exist either. We therefore present a library of idealized, high resolution $N$-body simulations of the tidal evolution of individual subhaloes that can be used to calibrate semi-analytical models and to complement cosmological simulations. The simulations focus on minor mergers, i.e., the mass of the subhalo is much smaller than that of the host halo, such that the impact of dynamical friction is negligible. This setup allows the adoption of a fixed analytical potential for modelling the host halo. The dynamical evolution of subhaloes is followed with $N$-body computations. In the library, four parameters, two of which characterize the subhalo orbit with respect to the host halo, and the two concentrations of the host- and subhalo, are varied over the ranges encountered in cosmological simulations. We show several representative examples from the library that illustrate the evolution of the subhalo mass and velocity dispersion profiles. Additionally, we make publicly available a pre-trained non-parametric model of the subhalo mass evolution based on random forest regression. This model is able to interpolate the simulation data at the 0.1\,dex level and provides efficient access to the data for further use in modelling.

Figures

Tables

Table 1. Table 1: Summary of the DASH data files available for each simulation

filename	description
subhalo_evo	position, velocity, bound mass fraction, and half-mass radius of subhalo as function of $t$
radprof_rho	subhalo density, $\rho (r,t)$, as function of $r$ and $t$, in units of $200{\rho }_{\mathrm{crit}}\left(0\right)$
radprof_m	enclosed mass of subhalo, , as function of $r$ and $t$, in units of $M_{\mathrm{vir},s}$
radprof_sigmar	radial velocity dispersion, ${\sigma }_r(r,t)$, as function of $r$ and $t$, in units of $V_{\mathrm{vir},s}$
radprof_sigmat	tangential velocity dispersion, ${\sigma }_t(r,t)$, as function of $r$ and $t$, in units of $V_{\mathrm{vir},s}$

Equations

$$\rho (r)={\rho }_0{\left(\frac{r}{r_{\mathrm{s}}}\right)}^{-1}{\left(1+\frac{r}{r_{\mathrm{s}}}\right)}^{-2},$$

$$M_{\mathrm{vir}}=\frac{4\pi }{3}{\Delta }_{\mathrm{vir}}{\rho }_{\mathrm{crit}}{r}_{\mathrm{vir}}^3.$$

$$t_{\mathrm{cross}}=0.978h^{-1}\mathrm{Gyr}{\left(\frac{{\Delta }_{\mathrm{vir}}}{200}\right)}^{-1/2}.$$

$$\Gamma (h,{\Delta }_{\mathrm{vir}})\equiv {\left(\frac{h}{0.678}\right)}^{-1}{\left(\frac{{\Delta }_{\mathrm{vir}}}{200}\right)}^{-1/2}.$$

$$T_{\mathrm{r}}=2{\int }_{r_{\mathrm{p}}}^{r_{\mathrm{a}}}\frac{\mathrm{d}r}{\sqrt{2[E-{\Phi }_{\mathrm{h}}(r)]-L^2/r^2}},$$

$${\Phi }_{\mathrm{h}}(r)=-V_{\mathrm{vir},\mathrm{h}}^2\frac{\ln (1+c_{\mathrm{h}}x)}{f(c_{\mathrm{h}})x}$$

$$f(c)=\ln (1+c)-c/(1+c).$$

$$f_{\mathrm{bound}}(t)>1.79\frac{c_{\mathrm{s}}^2}{f(c_{\mathrm{s}})}\left(\frac{ϵ}{r_{\mathrm{vir},\mathrm{s}}}\right)\left(\frac{r_{\mathrm{h}}(t)}{r_{\mathrm{vir},\mathrm{s}}}\right)$$

$$f_{\mathrm{bound}}(t)>\frac{0.054}{f(c_{\mathrm{s}})}{\left(\frac{c_{\mathrm{s}}}{10}\right)}^2\left(\frac{r_{\mathrm{h}}(t)}{r_{\mathrm{vir},\mathrm{s}}}\right),$$

$$N_{\mathrm{crit}}\equiv 80N^{0.2}.$$

$$f_{\mathrm{bound}}(t)>1.22\times 10^{-3}.$$

$$t=0.12\mathrm{Gyr}\Gamma (h,{\Delta }_{\mathrm{vir}})\times \mathrm{ID}.$$

$$\beta (r)\equiv 1-\frac{{\sigma }_{\mathrm{t}}^2(r)}{2{\sigma }_{\mathrm{r}}^2(r)},$$

$$R^2\equiv 1-\frac{\sum [\log (f_{\mathrm{b},\mathrm{sim}})-\log (f_{\mathrm{b},\mathrm{pred}}){]}^2}{\sum [\log (f_{\mathrm{b},\mathrm{sim}})-⟨\log (f_{\mathrm{b},\mathrm{sim}})⟩{]}^2}.$$

$$⟨\log (f_{\mathrm{b},\mathrm{sim}})⟩=\frac{1}{N_{\mathrm{test}}}\sum \log (f_{\mathrm{b},\mathrm{sim}}).$$

$${\rho }_0=\frac{c^3}{f(c)}{\Delta }_{\mathrm{vir}}{\rho }_{\mathrm{crit}},$$

$$\begin{array}[]{rcl}{\cal M}&\rightarrow&{\cal M}^{\prime}={\cal M}\,({\cal Q}^{3}_{\rm h}/{\cal Q}^{3}_{\rm s})\\ c_{\rm h}&\rightarrow&c^{\prime}_{\rm h}=c_{\rm h}\,{\cal Q}_{\rm h}\\ c_{\rm s}&\rightarrow&c^{\prime}_{\rm s}=c_{\rm s}\,{\cal Q}_{\rm s}\\ x_{\rm c}&\rightarrow&x^{\prime}_{\rm c}=x_{\rm c}\,{\cal Q}^{-1}_{\rm h}\\ \eta&\rightarrow&\eta^{\prime}=\eta\\ f_{\rm b}&\rightarrow&f^{\prime}_{\rm b}=f_{\rm b}\,(\Gamma/{\cal Q}^{3}_{\rm s})\,.\end{array}$$

$${\mathcal{Q}}^3\frac{f(c)}{f(c\mathcal{Q})}\left(\frac{{\Delta }_{\mathrm{vir}}^{\prime }}{200}\right){\left(\frac{h^{\prime }}{0.678}\right)}^2=1.$$

Full text

††pagerange: DASH: a library of dynamical subhalo evolution–A††pubyear: 2013

1 Introduction

In the cold dark matter (CDM) paradigm of cosmic structure formation, smaller perturbations collapse first to form virialized dark matter haloes, leading to a hierarchical assembly of haloes. When dark matter haloes assemble their mass by accreting smaller haloes, they build up a hierarchy of substructure, with subhaloes hosting sub-subhaloes, hosting sub-sub-subhaloes, etc. As these subhaloes orbit their hosts, they experience mass loss due to the combined effect of dynamical friction, tidal stripping and impulsive (tidal) heating (e.g., Mo et al., 2010). The resulting abundance and demographics of substructure depends on the microscopic properties of the dark matter particles, most importantly the free-streaming scale and the strength of dark matter self-interaction (see e.g., Knebe et al., 2008; Rocha et al., 2013; Bose et al., 2016). This is why many efforts are underway to quantify the amount of dark matter substructure using, among others, gravitational lensing (e.g. Vegetti et al., 2012; Shu et al., 2015; Hezaveh et al., 2016), gaps in stellar streams (e.g. Carlberg, 2012; Ngan & Carlberg, 2014; Erkal et al., 2016), and annihilation or decay signals of dark matter particles (e.g. Strigari et al., 2007; Pieri et al., 2008; Hayashi et al., 2016; Hiroshima et al., 2018). In addition, substructure is also directly related to the abundance and properties of satellite galaxies (e.g., Vale & Ostriker, 2006; Newton et al., 2018), and thus to the clustering amplitude of galaxies on small scales (see e.g., Benson et al., 2001; Berlind et al., 2003; Kravtsov et al., 2004). Hence, it is important that we are able to make accurate predictions for the abundance of dark matter substructure for a given cosmological model. Given the highly non-linear nature of the processes involved, this is ideally done using $N$-body simulations.

Modern, state-of-the-art cosmological $N$-body simulations predict that roughly 5-10 percent of a halo’s mass is bound up in substructure, with more massive host haloes having a larger subhalo mass fraction (e.g., Gao et al., 2004; Giocoli et al., 2010). In addition, subhaloes are found to be spatially anti-biased with respect to the dark matter, in that the radial number density profile, $n_{\rm sub}(r)$, is less centrally concentrated than the halo density profile (e.g., Diemand et al., 2004; Nagai & Kravtsov, 2005; Springel et al., 2008). It has also been concluded that the subhalo mass function has a universal form (Jiang & van den Bosch, 2016), albeit with a significant halo-to-halo variance at the massive end (Jiang & van den Bosch, 2017).

Although these trends are well understood (van den Bosch et al., 2005; Zentner et al., 2005; Jiang & van den Bosch, 2016), and the results from numerical simulations seem to be well converged (e.g., Springel et al., 2008; Onions et al., 2012), some issues remain. Foremost among these is the fact that subhaloes in numerical simulations typically experience complete disruption some time after accretion (e.g., Han et al., 2016; van den Bosch, 2017). Although it is often argued that this disruption is a physical consequence of either tidal stripping or tidal heating (Hayashi et al., 2003; Taylor & Babul, 2004; Klypin et al., 2015), others have argued that in a collisionless dark matter simulation, subhaloes should rarely ever completely disrupt. In particular, in van den Bosch et al. (2018, hereafter Paper I), we have demonstrated that neither tidal stripping nor tidal heating is expected to be able to completely unbind the central cusps of CDM substructure. The same conclusion was reached by Peñarrubia et al. (2010) using idealized, high-resolution numerical simulations. In van den Bosch & Ogiya (2018, hereafter Paper II), we used a large suite of similar, idealized simulations to demonstrate that the disruption of subhaloes in $N$-body simulations is predominantly numerical, and triggered by two independent aspects: discreteness noise and inadequate force softening.

An important finding of Paper II is that this artificial, numerical disruption may elude standard convergence tests, in that the method that is typically used to scale the force softening with the mass resolution is inadequate to overcome these problems. Artificial disruption can potentially have far-reaching consequences. After all, unless we can make accurate, and above all reliable, predictions regarding the abundance and structure of dark matter subhaloes, we will forfeit one of the main handles we have on learning about the nature of dark matter. In addition, artificial disruption is also a serious road-block for the small-scale cosmology program, which often relies on numerical simulations to predict the clustering strength of galaxies on small scales. A prime example of this is subhalo abundance matching, which assigns ‘mock’ galaxies to subhaloes identified in numerical simulations in order to predict galaxy-galaxy correlation functions (e.g., Vale & Ostriker, 2006; Conroy et al., 2006; Guo et al., 2010; Hearin et al., 2013). Although one may overcome the implications of artificial disruption by including ‘orphan’ galaxies (i.e., mock galaxies without an associated subhalo in the simulation), this seriously diminishes the information content of small-scale clustering, unless it is well understood how many orphans to add and where.

Unfortunately, it is not clear how to improve $N$-body codes such that the issue of artificial disruption can be avoided. As a consequence, it is difficult to gauge its potential impact on our predictions for the abundance of substructure. However, some insight can be gained from the semi-analytical models constructed by Jiang & van den Bosch (2016). Using accurate halo merger trees, this model uses a simple, orbit-averaged prescription for mass stripping to predict the evolved subhalo mass function. The overall normalization for the efficiency of mass stripping is calibrated by matching the model predictions to those from a high-resolution cosmological $N$-body simulation. In addition, the model includes a treatment of subhalo disruption, which is also calibrated to accurately reproduce the disruption in the simulation. If we use this model, but turn off the disruption (rather than disrupting the subhaloes, we continue to strip their mass), the resulting subhalo mass function is roughly a factor of two higher than with disruption (Green et al., in prep). Hence, if all disruption is indeed artificial, and if the mass stripping model used by Jiang & van den Bosch (2016) is roughly correct, numerical simulations may have been underpredicting the amount of surviving substructure by a factor of two. This would have far-reaching consequences for many areas of astrophysics. For instance, this factor of two is exactly what is needed to solve the ‘galaxy clustering crisis’ in subhalo abundance matching discussed in Campbell et al. (2018).

In order to make more reliable predictions, we need to develop more sophisticated semi-analytical models for the evolution of dark matter substructure. Numerous studies in the past have been devoted to this (e.g., Taylor & Babul, 2001; Peñarrubia & Benson, 2005; Zentner et al., 2005; Diemand et al., 2007; Kampakoglou & Benson, 2007; Gan et al., 2010; Pullen et al., 2014), but they all have one shortcoming in common: they all use the outcome of cosmological $N$-body simulations in order to calibrate one or more ‘fudge’ parameters in their model. And in doing so, their semi-analytical models inherit the shortcomings of the simulations; put differently, by construction the models are only as accurate as the simulations used for their calibration.

In an attempt to bypass this shortcoming, this paper presents a large database (called DASH, for Dynamical Aspects of SubHaloes), of more than 2,000 idealized, high-resolution simulations of the tidal evolution of individual subhaloes. The simulations cover most of the parameter space (mass ratios, orbital parameters and halo concentrations) relevant for modelling the tidal evolution of subhaloes as they orbit their hosts, and each simulation is evolved with sufficient numerical resolution that discreteness noise and force softening do not adversely affect their outcome. The primary goal of DASH is to enable a more accurate calibration and validation of analytical treatments of tidal stripping and heating, thereby allowing for the construction of new and improved semi-analytical models for dark matter substructure that are not hampered by artificial disruption. In addition, we provide a non-parametric model of subhalo mass evolution, using random forest regression, along with the simulation data. This model describes the simulation results at the $\sim$ 0.1 dex level and can be readily used in further modelling.

This paper is organized as follows: §2 describes the simulation setup and gives an overview of the DASH library. §3 presents a few examples, highlighting the type of data that is available. In §4, we demonstrate and validate the performance of a random forest regression model trained on the simulation data to predict bound mass fractions as a function of time since accretion for given orbital parameters and given concentrations of the sub- and host haloes. Finally, §5 summarizes the paper and discusses the future outlook.

2 Overview of the simulations

This section presents an overview of the DASH library of idealized, collisionless $N$-body simulations of halo minor mergers. After describing the initial conditions (§2.1.1), the simulation code (§2.1.2), and the analysis and products (§2.1.3), we describe the method used to sample the parameter space (§2.2) and the data format of the DASH library (§2.3).

2.1 Simulation setup

2.1.1 Halo profiles and merger set-up

Each of the simulations run as part of the DASH library follows an individual $N$-body subhalo as it orbits the fixed, external potential of a host halo. Both the host halo and the initial (prior to the onset of tidal stripping) subhalo are assumed to be spherical, and to have a Navarro-Frenk-White (NFW; Navarro et al., 1997) density profile

[TABLE]

where $r_{\rm s}$ and $\rho_{0}$ are the characteristic scale radius and density, respectively. We define the virial radius, $r_{\rm vir}$, as the radius inside of which the average density is $\Delta_{\rm vir}=200$ times the critical density given by $\rho_{\rm crit}=(3\,H^{2}_{0}/8\pi G)$, where $H_{0}$ and $G$ are the Hubble constant and gravitational constant, respectively. The virial mass of the halo is given by

[TABLE]

We emphasize, though, that the DASH simulations also apply to other values of $\Delta_{\rm vir}$, as detailed in §2.2 and Appendix A. The halo concentration is defined as $c\equiv r_{\rm vir}/r_{\rm s}$, and the virial velocity is defined as the circular velocity at the virial radius, $V_{\rm vir}=\sqrt{GM_{\rm vir}/r_{\rm vir}}$. The mass ratio between the host halo and the initial subhalo is specified by ${\cal M}\equiv M_{\rm vir,h}/M_{\rm vir,s}$, where, as throughout this paper, subscripts ‘h’ and ‘s’ indicate properties of the host- and subhaloes, respectively.

Initial conditions are generated under the assumption that the NFW subhalo is isolated and has an isotropic velocity distribution, such that its phase-space 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). When computing the DF, we follow Kazantzidis et al. (2006) and assume that the initial NFW subhalo has an exponentially decaying density profile for $r>r_{\rm vir,s}$. We assume that the total mass in this exponential extension is $0.05M_{\rm vir,s}$, where $M_{\rm vir,s}=M(<r_{\rm vir,s})$. Requiring a smooth transition in the density profile at $r=r_{\rm vir,s}$ then determines the scale radius of the exponential decay. Note, though, that when we sample particles from the DF thus computed, we apply a hard truncation at $r=r_{\rm vir,s}$ (i.e., no particles are sampled beyond the subhalo’s virial radius). Consequently, the system will deviate somewhat from perfect equilibrium near the truncation radius. However, since we embed the subhalo in an external tidal field, with a corresponding tidal radius that typically lies well inside of $r_{\rm vir,s}$, there is little virtue to having a subhalo whose outskirts are in perfect equilibrium. In fact, one might argue that it is more realistic to truncate the subhalo at the tidal radius corresponding to its initial position within the host halo. However, in addition to the tidal radius being ill-defined (see Paper I for discussion), we have demonstrated in Paper II that none of this matters; the simulation outcome is insensitive to whether we truncate the subhalo at the virial radius or at the initial tidal radius.

Throughout we adopt model units in which the gravitational constant, $G$, the initial virial radius of the subhalo, $r_{\rm vir,s}$, and the initial virial mass of the subhalo, $M_{\rm vir,s}$, are all unity. With this choice, the initial virial velocity of the subhalo, $V_{\rm vir,s}$ is unity, while the host halo has $V_{\rm vir,h}={\cal M}^{1/3}$. Both have the same crossing time, $t_{\rm cross}\equiv r_{\rm vir}/V_{\rm vir}=1$. In physical units, the crossing time is

[TABLE]

where $h=H_{0}/(100\>{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1})$. In what follows, whenever we quote time scales in physical units, we adopt $h=0.678$ (Planck Collaboration et al., 2016) and $\Delta_{\mathrm{vir}}=200$, which implies that a time interval of $(\Delta t)_{\rm model}=1$ corresponds to $1.44\>{\rm Gyr}$. Note, though, that one can scale all the physical time scales quoted in this paper to other values of $H_{0}$ and $\Delta_{\rm vir}$ by simply multiplying the values quoted by the factor

[TABLE]

Put differently, the mapping between time scales in DASH model units and physical units, is given by: $(\Delta t)_{\rm physical}=1.44\>{\rm Gyr}\,\Gamma(h,\Delta_{\mathrm{vir}})\,(\Delta t)_{\rm model}$.

The DASH simulations span a wide range in orbital energy, $E$, and angular momentum, $L$. For convenience, we characterize the orbits using the following two dimensionless quantities:

• •

$x_{\rm c}\equiv r_{\rm c}(E)/r_{\rm vir,h}$, the radius of the circular orbit corresponding to the orbital energy, $E$, expressed in terms of the virial radius of the host halo.

• •

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

We initially position the subhalo at the apocentre of its orbit, and follow its dynamical evolution for a period of $T_{\rm sim}=36\>{\rm Gyr}$. The orbit’s radial period is given by

[TABLE]

(e.g., Binney & Tremaine, 2008), with $r_{\mathrm{p}}$ and $r_{\mathrm{a}}$ the pericentric and apocentric radii of the orbit, respectively, and $\Phi_{\mathrm{h}}(r)$ the gravitational potential due to the host halo. The latter is given by

[TABLE]

where $x=r/r_{\rm vir,h}$ and

[TABLE]

Fig. 1 plots $T_{\mathrm{r}}$ as function of $x_{\mathrm{c}}$ for several values of $\eta$ and $c_{\mathrm{h}}$. Note that, to good approximation, $T_{\mathrm{r}}\sim 6.7\>{\rm Gyr}\,x_{\rm c}^{1.15}$, with only a very weak dependence on $\eta$ or $c_{\mathrm{h}}$. Hence, for the range of orbits covered by DASH, which have $x_{\rm c}\in[0.5,2.0]$ (see §2.2 below), $T_{\rm sim}=36\>{\rm Gyr}$ corresponds to between 2.5 and 12 radial periods. The initial (Cartesian) vectors of the subhalo position and velocity with respect to the host halo are given by ${\bf X}=(r_{\rm a},0,0)$ and ${\bf V}=(0,L/r_{\rm a},0)$, respectively, such that the subhalo orbit is confined to the $x-y$ plane. For each simulation we output a total of 301 snapshots, with the time interval between snapshots fixed at 0.12 Gyr.

2.1.2 Numerical techniques

All DASH simulations have been carried out using a tree code (Barnes & Hut, 1986) specifically developed for graphics processing unit (GPU) clusters (Ogiya et al., 2013). The code uses CPU cores to construct octree structures of the $N$-body particles, while GPU cards are used to compute gravitational accelerations through tree traversal. We employ 1,048,576 particles in each simulation. Forces between particles are softened using a Plummer softening length $\epsilon=0.0003r_{\rm vir,s}$ (see Paper II) and the opening angle of the tree algorithm is set to $\theta=0.7$. Orbit integration uses a leapfrog scheme with a global, adaptive time step $\Delta t=0.2\sqrt{\epsilon/a_{\mathrm{max}}}$ (Power et al., 2003). Here $a_{\mathrm{max}}$ is the maximum, absolute value of acceleration among all particles at that time222Typically the time step is fairly constant throughout the entire simulation, except for a small (typically factor $\sim 2$) decrease in $\Delta t$ during a high-speed pericentric passage of the host halo.. As we have demonstrated in Paper II, these numerical parameters are adequate to properly resolve the tidal evolution of subhaloes. In order to verify this, we have run a subset of our simulations at ten times better mass resolution (using $N\sim 10^{7}$ particles), which yields results that are indistinguishable from our nominal mass resolution.

2.1.3 Data analysis and products

For each simulation output we compute the bound mass fraction of the subhalo, $f_{\rm bound}(t)$, using the iterative method described in detail in Appendix A of Paper I. Briefly, the centre-of-mass position, ${\bf r}_{\rm com}$, and velocity, ${\bf v}_{\rm com}$, are computed using the five percent most bound particles. These quantities are subsequently used to compute the binding energy of each particle. This is iterated until the changes in ${\bf r}_{\rm com}$ and ${\bf v}_{\rm com}$ are smaller than $10^{-4}r_{\rm vir,s}$ and $10^{-4}V_{\rm vir,s}$, respectively, which typically requires 3-10 iterations. The bound fraction, $f_{\rm bound}(t)$, is then defined as the fraction of the original particles that at time $t$ have a negative binding energy.

All simulations in the DASH library initially have sufficient numerical resolution to properly resolve the dynamics of the dark matter subhalo. However, as highlighted in Paper II, simulation results can become unreliable once the subhalo has lost a significant fraction of its initial mass due to tidal evolution. Paper II presented the two criteria that a subhalo in a numerical simulation needs to satisfy in order for its dynamical evolution to be reliable. The first criterion, which is motivated by the work of Power et al. (2003), tests whether the softening length is sufficiently small to resolve the relevant accelerations. It translates into a requirement for the bound mass fraction given by

[TABLE]

with $r_{\rm h}(t)$ the half-mass radius of the (bound part of) the subhalo at time $t$. Using that all DASH simulations adopt $\varepsilon=0.0003r_{\rm vir,s}$, we thus have that the simulation results need to satisfy

[TABLE]

for the results to be deemed reliable. The second criterion is related to discreteness noise, and puts a constraint on the number of bound particles in the subhalo. In particular, it states that the number of bound particles needs to exceed

[TABLE]

with $N$ the number of particles in the initial subhalo. Once the number of bound particles falls below this critical value, the subhalo starts to experience a runaway instability, triggered by discreteness noise, which quickly leads to its demise (i.e., artificial disruption). Since all DASH simulations have $N=1,048,576$, we have that $N_{\mathrm{crit}}=1,280$, which implies that the DASH simulations are only reliable for

[TABLE]

As discussed in §2.3 below, the DASH library contains, for each output, the bound mass fraction and the half-mass radius of the subhalo, which the user can use to test whether the output satisfies both criteria. More than 99.5% of all the simulation outputs available in the DASH library satisfy both criteria (9) and (11).

2.2 The DASH parameter space

The simulations described above are characterized by six parameters: the mass and concentration of the host halo, $M_{\rm vir,h}$ and $c_{\rm h}$, respectively, the mass and concentration of the subhalo, $M_{\rm vir,s}$ and $c_{\rm s}$, respectively, and the orbital parameters $x_{\rm c}$ and $\eta$. In order to limit the numerical cost of sampling this six-dimensional parameter space without sacrificing the volume sampled, we adopt a strategy that devotes more computational efforts to the regions of parameter space with higher probability. In doing so, we are guided by physical considerations and results from cosmological simulations.

2.2.1 Initial mass ratio

Since dynamical timescales and the impact of tidal forces depend only on density, and since all haloes, by definition, have the same virial density, our simulation results should be independent of $M_{\rm vir,h}$ and $M_{\rm vir,s}$, significantly reducing the dimensionality of our parameter space. However, there is one important caveat. In reality, a subhalo orbiting a host halo experiences a dynamical drag force caused by the (perturbed) matter in the host halo (‘dynamical friction’; e.g., Chandrasekhar, 1943) and by its own stripped material (‘self-friction’; e.g., Fujii et al., 2006; Fellhauer & Lin, 2007). The impact of these drag forces, which result in orbital decay, increase strongly with decreasing mass ratio ${\cal M}=M_{\rm vir,h}/M_{\rm vir,s}$. When {\cal M}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}100, though, their impact can safely be ignored (e.g., Mo et al., 2010). Throughout, we therefore restrict ourselves to simulations with ${\cal M}=1000$. As we specifically demonstrate in §3.1, these simulations accurately capture the tidal evolution of subhaloes for any {\cal M}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}100, independent of the absolute value of $M_{\rm vir,h}$.

2.2.2 Sampling orbital parameters

Using a state-of-the-art cosmological simulation, Jiang et al. (2015) studied the orbital parameters of dark matter subhaloes at their moment of accretion into their host halo (see also e.g., Tormen, 1997; Zentner et al., 2005; Khochfar & Burkert, 2006; Wetzel, 2011; van den Bosch, 2017). In particular, Jiang et al. (2015) measured the radial and tangential components, $V_{\mathrm{r}}$ and $V_{\mathrm{\theta}}$, of the relative velocity vector between the host- and subhaloes at infall. We convert their bivariate distribution of $V_{\mathrm{r}}/V_{\rm vir,h}$ and $V_{\mathrm{\theta}}/V_{\rm vir,h}$ to the corresponding bivariate distribution of $x_{\rm c}$ and $\eta$, assuming an NFW density profile for the host halo with a concentration parameter $c_{\mathrm{h}}=5$. The resulting PDF $P(x_{\rm c},\eta)$ is shown in Fig. 2, where the small panels show the corresponding, marginalized distributions of $x_{\rm c}$ (side panel) and $\eta$ (bottom panel).

In order to ensure that our sampling of $x_{\rm c}$ and $\eta$ covers the entire range relevant for modelling the assembly and evolution of dark matter substructure, we proceed as follows. We sample $\eta$ linearly over the entire range from $\eta=0$ (purely radial orbit) to $1.0$ (circular orbit) in 10 steps of $\Delta\eta=0.1$, and $x_{\rm c}$ logarithmically over the range from $x_{\rm c}=0.5$ to $2.0$, in 10 steps of $\Delta\log x_{\rm c}\simeq 0.06$. Hence, there are a total of 121 combinations of $(x_{\rm c},\eta)$. Together with 121 combinations of $(c_{\rm h},c_{\rm s})$ (see below), this would imply that we need to run more than 14,500 simulations. In addition to being prohibitively expensive, this is also unwarranted, as we will demonstrate in §4. Instead, we significantly reduce the number of simulations by focusing most of our efforts on the more likely parts of parameter space. The red crosses in Fig. 2 indicate the 45 most likely combination of $x_{\rm c}$ and $\eta$. For each of these we run the 45 most likely combinations of $(c_{\rm h},c_{\rm s})$, using the sampling strategy detailed below. The black crosses, on the other hand, indicate the combinations of $(x_{\rm c},\eta)$ for which we only run one realization of halo concentrations, namely the one for $(c_{\rm h},c_{\rm s})=(5.0,10.0)$.

2.2.3 Sampling halo concentrations

In deciding on how to sample $c_{\rm h}$ and $c_{\rm s}$, we are again guided by cosmological simulations. These reveal that haloes of fixed mass have a log-normal distribution of halo concentrations, with a scatter of $\sim 0.12$ dex, and with a median that decreases with halo mass roughly as $c\propto M^{-0.1}$ (e.g., Bullock et al., 2001; Prada et al., 2012; Dutton & Macciò, 2014; Diemer & Kravtsov, 2015; Ludlow et al., 2016). For DASH we sample both $c_{\rm s}$ and $c_{\rm h}$ using $\log c=0.5,0.6,0.7,...,1.5$. This covers the range 3.1\lower 3.01385pt\hbox{;\buildrel<\over{\sim};}c_{\rm vir}\lower 3.01385pt\hbox{;\buildrel<\over{\sim};}31.5, which is adequate for the vast majority of all haloes in our mass range of interest (roughly $10^{7}<M_{\rm vir}/(\>h^{-1}\rm M_{\odot})<10^{15}$). This sampling, though, yields a total of 121 unique pairs of $(c_{\rm h},c_{\rm s})$. Similar to our sampling of the orbital parameters, we focus most of our efforts on the ‘most likely’ combinations of $c_{\rm h}$ and $c_{\rm s}$. First, using the fitting function of Fakhouri et al. (2010), we compute the halo merger rate as a function of the mass ratio of the merging haloes. We only consider minor mergers with a mass ratio ${\cal M}>100$, representative of our simulations. Next, we draw concentrations for each of the merging haloes using the concentration-mass relation of Ludlow et al. (2016), and account for log-normal scatter with $\sigma_{\log c}=0.12$. Using this Monte-Carlo procedure, we compute the PDF for $\log[c_{\rm s}/c_{\rm h}]$. Results for different values of $M_{\rm vir,h}$ are shown in Fig. 3, and are in good agreement with the results obtained from cosmological simulations (e.g., Paper I). Note that more massive host haloes tend to have somewhat higher concentration ratios on average. However, the strength of this mass dependence is weak compared to the width of the individual PDFs, and we therefore ignore it in what follows. Fig. 4 shows the joint distribution of $c_{\mathrm{h}}$ and $c_{\mathrm{s}}$, derived using the PDF for $\log[c_{\rm s}/c_{\rm h}]$ for a host halo of mass $M_{\rm vir,h}=10^{12}\>h^{-1}\rm M_{\odot}$ and adopting a uniform distribution333Although halo concentrations follow a log-normal distribution at fixed mass, we adopt (for simplicity) a uniform distribution such that the simulations are relevant for a wide range in halo masses. in $c_{\rm h}$. Black crosses indicate combinations of $(c_{\rm h},c_{\rm s})$ for which we have run simulations for all 45 most likely orbital parameter combinations $(x_{\rm c},\eta)$ (the red crosses in Fig. 2). White crosses, on the other hand, indicate combinations of $(c_{\rm h},c_{\rm s})$ for which we have run only one orbital configuration, namely with $(x_{\rm c},\eta)=(1.0,0.5)$.

2.3 Data structure

Using the sampling strategy outlined above, there are a total of 2,177 idealized $N$-body simulations that we have run. It is hierarchically structured, with different layers of directories corresponding to the various adopted orbital or concentration parameters. Each directory presents the results of the simulation in the form of five simple text files (see Table 1), and each file contains a header listing the parameter set used in the simulation, as well as a brief description of the data file. Rather than accessing the results for individual simulations, the user can also download a single compressed archive file ($\sim 0.6$ GB) containing the results for all 2,177 simulations as well as a Python routine that uses a non-parametric model based on random forest regression (see §4), to predict the bound mass fraction of a subhalo as a function of time, for any configuration captured by our parameter space. In what follows we briefly describe the information that is provided for each simulation.

The file subhalo_evo lists the temporal evolution of several subhalo properties. The first column lists the ID of each simulation snapshot, which can be used to compute the corresponding physical time as

[TABLE]

Columns 2-4 and 5-7 list the Cartesian components of the position and velocity vectors of the centre-of-mass of the subhalo, respectively. The former and latter are normalized by the virial radius and velocity of the host halo, $r_{\rm vir,h}$ and $V_{\rm vir,h}$, respectively. Finally, columns 8 and 9 list the subhalo bound mass fraction, $f_{\rm bound}(t)$, and the subhalo’s half-mass radius, $r_{\rm h}(t)$. The latter is normalized by the virial radius of the initial subhalo, $r_{\rm vir,s}$. Note that these outputs describe the tidal evolution of a (spherical NFW) subhalo of any mass, $M_{\rm vir,s}$, in a (spherical NFW) host halo of any mass, $M_{\rm vir,h}$, as long as {\cal M}=M_{\rm vir,h}/M_{\rm vir,s}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}100. For smaller mass ratios dynamical friction, which is not accounted for, is important.

The other data available for each simulation in the DASH library are four text files that list the time evolution of four different radial profiles: the normalized density, $\rho(r,t)/[200\,\rho_{\rm crit}(0)]$, the enclosed, normalized mass profile, $M(<r,t)/M_{\rm vir,s}$, and the normalized radial and tangential velocity dispersion profiles, $\sigma_{\mathrm{r}}(r,t)/V_{\rm vir,s}$ and $\sigma_{\mathrm{t}}(r,t)/V_{\rm vir,s}$, respectively. The positions and velocities of all particles are defined with respect to the centre-of-mass position and velocity of the subhalo, and only bound particles are taken into account. The first row of each of the radprof files lists the radial bins, which span the range $-2.95\leq\log(r/r_{\rm vir,s})\leq 0.95$ in equally spaced bins of width $\Delta\log(r/r_{\rm vir,s})=0.1$. The subsequent rows list the radial profiles for each of the 301 simulation snapshots, with row $j$ corresponding to snapshot ID=$j-2$. The corresponding physical time is given by equation (12).

3 Some examples from the simulations

While classical dynamical friction is absent in the simulations, another source of friction caused by the mass stripped from the subhalo may alter the mass and orbital evolution of the subhalo in the tidal interactions. §3.1 studies the validity condition that must be satisfied in order to neglect this friction force. The subsequent subsections present some examples of the DASH simulations, namely the mass evolution (§3.2) and radial profiles (§3.3) of the tidally stripped subhaloes.

3.1 Dependence on the initial mass ratio

As mentioned in §2.1, the host halo is modelled as an analytical NFW potential and the initial mass ratio between the host- and subhaloes, ${\cal M}\equiv M_{\mathrm{vir,h}}/M_{\mathrm{vir,s}}$, is fixed at $10^{3}$ for all simulations. In this subsection, we discuss the validity condition for these assumptions.

In actual mergers, friction forces may alter the orbit of subhaloes. The classical dynamical friction (Chandrasekhar, 1943) is caused by wakes formed in the density field of the larger host system due to the gravitational force of smaller objects, such as supermassive black holes, galactic bars, and dark matter subhaloes (e.g., Weinberg & Katz, 2007; Antonini & Merritt, 2012; Ogiya & Burkert, 2016). This kind of friction force is absent in our simulations because we employ an analytical host potential and the response in the density field of the host system due to the gravitational force of the subhalo is therefore not taken into account. The timescale of orbital decay due to this dynamical friction is roughly equal to ${\cal M}\tau_{\rm ff}$ (e.g., Mo et al., 2010), where $\tau_{\rm ff}$ is the free-fall time of the host halo, which is of the order of 2-3 Gyr at $z=0$. Since all our simulations adopt ${\cal M}=10^{3}$, is it clear that we can safely neglect the effects of dynamical friction, and thus use an analytical potential to model the host halo.

However, in addition to this ‘classical’ dynamical friction, which is not accounted for in our simulations, there is an additional friction force that may contribute, and which is included in our simulations. This force is due to the mass stripped from the subhalo itself, and we therefore refer to this force as ‘self-friction’. Fig. missing 5 shows an example of the tidal evolution of a subhalo in the DASH simulation with $x_{\mathrm{c}}=1.0$, $\eta=0.5$, $c_{\mathrm{h}}=5$, $c_{\mathrm{s}}=10$ and ${\cal M}=10^{3}$. Time progresses from the top left to the bottom right, as indicated, while colors indicate column density (in arbitrary units). As is evident, the tidal field of the host halo strips mass from the subhalo, giving rise to a leading and a trailing tidal arm that roughly trace out the subhalo’s orbit. The gravitational force from this stripped material on the bound remnant, can result in a net deceleration, thus giving rise to a self-friction (Fujii et al. 2006; Fellhauer & Lin 2007; Paper II).

In order to estimate the impact of self-friction, we have performed a number of simulations in which we vary the mass ratio, ${\cal M}$, while keeping all other parameters fixed to $x_{\mathrm{c}}=1.0$, $\eta=0.5$, $c_{\mathrm{h}}=5$ and $c_{\mathrm{s}}=10$. Each subhalo is simulated with 1,048,576 particles and uses our fiducial softening and time-stepping. The results are shown in Fig. missing 6, which plots the orbital radius (upper panel) and bound mass fraction (lower panel) as functions of time. The black dotted line corresponds to the results from a ${\cal M}=10^{3}$ simulation with an order of magnitude more particles ($N=10,485,760$), and is shown to demonstrate that these results are well converged; the black dotted line lies exactly on top of the results from our fiducial resolution simulation, shown as a purple solid line. In the absence of any friction, the results from all these simulations should all be identical. This follows from the scale-free nature of gravity, the fact that (sub)haloes with the same concentration parameter but different mass, have identical density profiles as function of the normalized radius $r/r_{\rm vir}$, and the fact that the tidal radius scales with the ratio of the densities of host and subhalo444A comprehensive review is found in Paper I (see also e.g., Binney & Tremaine 2008).

While the simulations with {\cal M}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}100 are indeed indistinguishable, those with smaller mass ratios clearly deviate, both in terms of their orbit and in terms of the bound mass fraction. Notably, the apocentric distance reached after the first pericentric passage becomes smaller for smaller ${\cal M}$, indicating that the subhalo is experiencing self-friction due to its own stripped material. Over time this self-friction causes the orbit to shrink, which reduces the pericentre of the orbit, thereby exposing the subhalo to a stronger tidal tidal, which in turn results in a reduced bound mass fraction. As is evident from Fig. missing 6, the impact of self-friction is more pronounced for smaller values of ${\cal M}$. This is easy to understand; as ${\cal M}$ decreases, the force from the stripped material on the bound remnant relative to that from the host halo increases, causing a more pronounced deceleration along its orbit. We present a more detailed study of self-friction in a forthcoming paper (Miller et al. 2019, in preparation). Relevant for this study, though, is that self-friction is negligible for {\cal M}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}100, and that the evolution of the bound mass fraction does not depend on ${\cal M}$ as long as this is the case. Hence, the DASH simulations presented here, which all have ${\cal M}=10^{3}$, are valid for any mass ratio {\cal M}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}100.

3.2 Bound mass fraction

Fig. missing 7 demonstrates the evolution of the bound mass fraction of the subhalo, $f_{\mathrm{b}}$, obtained in some example simulations. Fixing the structural parameters of the host- and subhaloes (upper panels), the subhalo mass is more significantly decreased on more tightly bound (smaller $x_{\mathrm{c}}$) or radial (smaller $\eta$) orbits, since the orbits have smaller pericentres and the subhalo feels a stronger tidal force from the host potential. The lower panels illustrate that when the subhaloes are on the same orbit, a larger $c_{\mathrm{h}}$ (smaller $c_{\mathrm{s}}$) leads to more significant mass loss because of the stronger tidal force of the host halo (more loosely bound subhalo structure). These results, which only represent a tiny fraction of the entire DASH library, are consistent with intuition.

Less intuitive is the behavior in the bound mass fraction on a highly radial orbit (i.e., small $\eta$) close to a pericentric passage. In some cases, the bound mass fraction fluctuates wildly, dropping steeply, only to increase again immediately thereafter (see for instance the purely radial, $\eta=0$ orbit in the upper-left panel). This arises i) because the method to compute the bound mass and orbit of the subhalo allows re-binding of particles (see Appendix A of Paper I for details), and ii) because the subhalo is impulsively heated by tidal shocking, especially on radial orbits (e.g., Spitzer, 1987; Gnedin et al., 1999), which leads to the strong reduction of the bound mass. A fraction of the temporarily evaporated subhalo mass is re-bound during the subsequent re-virialization process of the subhalo.

Another interesting feature apparent in Fig. missing 7 is that, except for the step-like behavior near pericentre, $f_{\mathrm{b}}$ roughly behaves like a power-law function for $t/T_{\mathrm{r}}>1$. The slope of this power-law function clearly depends on the orbital and concentration parameters, and the DASH library, which contains simulation data for a large parameter space, can be used to calibrate these scaling relations, which in turn can be used to model the tidal evolution of substructure in (semi-)analytical models of structure formation. An alternative method, which we explore in §4 below, is to use machine learning algorithms, such as random forest regression, to process and distill the huge amount of data available in the DASH library. As we demonstrate, this allows for reasonably accurate predictions for the bound mass fraction of subhaloes as a function of their time since accretion.

3.3 Radial profiles

The DASH library also includes density and velocity dispersion profiles of tidally stripped subhaloes. An example is shown in Fig. missing 8, which plots the profiles for a subhalo for three snapshots of the simulation with $c_{\mathrm{h}}=5.0$, $c_{\mathrm{s}}=10.0$, $x_{\mathrm{c}}=1.00$ and $\eta=0.5$. Each of these snapshots corresponds to an epoch at which the subhalo is near its apocentre, at which point the subhalos is in a fairly relaxed state (e.g., Aguilar & White, 1986; Peñarrubia et al., 2009), and the corresponding bound mass fractions are indicated in the middle panel. The top panel of Fig. missing 8 shows the subhalo’s density profile, normalized by $200\rho_{\rm crit}$. Note how tidal stripping mainly removes mass from the outskirts, while leaving the central densities almost unaffected. For comparison, the dashed and dotted curves are the model predictions of Hayashi et al. (2003) and Peñarrubia et al. (2010), respectively. Both models suggest that the density profile of a stripped subhalo only depends on the initial density profile (prior to stripping) and the present bound mass fraction. Whereas the model by Hayashi et al. (2003) fits the profiles extremely well, the model by Peñarrubia et al. (2010) predicts a shallower outer density profile at later times. However, as we demonstrate in Green et al. (in preparation), neither the Hayashi et al. (2003) nor the Peñarrubia et al. (2010) model can adequately describe the evolution of the subhalo density profile under all conditions encountered in the DASH library, and we therefore develop a new and improved model based on the entire set of over 2,000 DASH simulations. The middle and bottom panels of Fig. missing 8 show the radial velocity dispersion profiles and the corresponding profiles of the velocity anisotropy parameter

[TABLE]

(Binney & Tremaine, 2008). Note that, by construction, the initial subhalo is isotropic ($\beta=0$) throughout. At later times, as its bound mass fraction decreases, the subhalo becomes more and more tangentially anisotropic ($\beta<0$) in its outskirts, while the radial velocity dispersion profile decreases on all scales. Hence, the bound remnant becomes colder and colder as more and more mass is stripped, and since particles on more radial orbits reach larger apocentric distances, they are more likely stripped, thereby causing the remnant to become tangentially anisotropic. Finally, the black, dashed lines show the results from a simulation with 10 times higher mass resolution. The fact that the resulting profiles are indistinguishable from those of the nominal resolution simulation indicates that the DASH simulations are well converged.

4 An application of machine learning to the DASH library

The DASH library is a homogeneously structured, large dataset that can easily be explored with machine learning (ML). As a first example, we apply a commonly-used ML regression method, random forests (RF; Breiman, 2001) as implemented in scikit-learn (Pedregosa et al., 2011)555https://scikit-learn.org, to predict $f_{\mathrm{b}}$ as a function of five features: $t/T_{\mathrm{r}}$, $x_{\mathrm{c}}$, $\eta$, $c_{\mathrm{h}}$ and $c_{\mathrm{s}}$. Because $c_{\mathrm{h}}$, $c_{\mathrm{s}}$ and $x_{\mathrm{c}}$ are equally binned in logarithmic space in the DASH library (see §2.2), we adopt $\log{(c_{\mathrm{h}})}$, $\log{(c_{\mathrm{s}})}$ and $\log{(x_{\mathrm{c}})}$ as the actual features. Furthermore, we train the model to predict $\log{(f_{\mathrm{b}})}$ as the target since $f_{\mathrm{b}}$ varies over several orders of magnitude. In the following analysis, only data points satisfying the numerical criteria of Paper II are included (a brief summary is given in §2.3).

The ML algorithm we adopt, RF, is based on decision trees (Quinlan, 1986). While the decision tree method is intuitive and useful, trained models tend to be overfitted, i.e., the training data set is very accurately reproduced while poor predictions are made for untrained cases. In order to avoid this overfitting issue, RF constructs an ensemble of decision trees and adopts the mean prediction of individual decision trees as the final prediction. The ensemble consists of 20 decision trees with a maximum depth of 20. Here, the depth of the tree is the number of layers from a root to a leaf. The other hyperparameters are set to their default values in the scikit-learn implementation.

In order to increase the confidence in our trained model, we adopt group $k$-fold cross validation (e.g., Browne, 2000, and references therein). First, the full data set is divided into $k$ subsets composed of data points from randomly selected simulations. Note that all time steps from each simulation are assigned to one group such that they are all placed within the same subset and thus are not split between multiple subsets. Then, the RF is trained with $k-1$ subsets and the trained model is tested on the remaining subset. This training and test procedure is iterated $k$ times and the performance ($R^{2}$ score, see below for details) is measured as the average of the $k$ models. We set $k=5$ and the number of simulations in each subset is almost the same (435 or 436). To further verify the robustness of our trained model, we perform a test with a reduced number of $t/T_{\mathrm{r}}$ bins. While the other four parameters ($c_{\mathrm{h}}$, $c_{\mathrm{s}}$, $x_{\mathrm{c}}$, and $\eta$) have 11 bins for each, $t/T_{\mathrm{r}}$ has 301 bins in the DASH library. This might artificially weight $t/T_{\mathrm{r}}$ more strongly relative to the other features. To assuage this concern, we also train a model in which we use only 11 of the 301 epochs, equally spaced in time, such that all five features have 11 bins.

Fig. missing 9 compares the evolution of $f_{\mathrm{b}}$ in the $N$-body simulations with the predictions made by the trained RF model and depicts that the model is predictive at the $\sim$ 0.1 dex level. Note that the model predictions shown in Fig. missing 9 are of simulations included only in the test fold for this particular trained model. In our case, a measure of accuracy of the trained model, the $R^{2}$ score (coefficient of determination) is defined as follows:

[TABLE]

The summation runs over all data points in the test set. The mean value in the test set is simply given as

[TABLE]

Here, $N_{\mathrm{test}}=N_{\mathrm{snap}}\times N_{\mathrm{sim}}/5\approx 131,000$ is the number of data points in the test set and $N_{\mathrm{sim}}(=2,177)$ and $N_{\mathrm{snap}}(=301)$ represent the number of simulations in the DASH library and the number of snapshots in each simulation, respectively. We emphasize again that only the data points fulfilling the numerical convergence criteria (see Paper II and §2.1.3) are included in this analysis. The trained model yields $R^{2}>0.98$ for the test set in all five cross-validation cases, indicating that the trained model works well for untrained cases within the covered parameter space. The $R^{2}$ score does not change even if we reduce the number of $t/T_{\mathrm{r}}$ bins by a factor of 30, verifying that the higher number of bins in $t/T_{\mathrm{r}}$ does not matter for our model.

scikit-learn also reports the importance of features, i.e., how much the model depends on each feature for predicting $f_{\mathrm{b}}$. The derived importance for each feature is [$t/T_{\mathrm{r}}$, $\log{(c_{\mathrm{h}})}$, $\log{(c_{\mathrm{s}})}$, $\log{(x_{\mathrm{c}})}$, $\eta$] = $[0.46,0.05,0.21,0.09,0.19]$, meaning that $f_{\mathrm{b}}$ depends strongly on time while $c_{\mathrm{s}}$ and $\eta$ play the most dominant roles among the four given parameters. This information may be useful in constructing semi-analytical models and in determining the parameter sets to be studied in the subsequent expansion of the DASH library.

The model trained with the full data set is saved in the file named dash_fb_rf.joblib and available within the Jupyter notebook, dash_fb_predict.ipynb, in the DASH library. Inputting the four parameters we vary, $x_{\mathrm{c}}$, $\eta$, $c_{\mathrm{h}}$, and $c_{\mathrm{s}}$, as well as the time of interest, one obtains the expected trajectories of the mass evolution of tidally stripped subhaloes. While the model is accurate and easy to use, it is not robust for extrapolation beyond the sampled region of parameter space because of the nature of decision tree-based algorithms, which can only interpolate between the data points in the training set. The prediction made by RF corresponds to the mean prediction of individual decision trees and hence it does not work in the parameter space where no data points are found. Other ML algorithms, e.g., support vector machines and others based on neural networks, are needed to construct models that can be extrapolated, but these frameworks generally require more complicated data preprocessing treatment and tuning of hyperparameters than those for RF to obtain good models.

Subsequent studies will use the other types of data available in DASH, such as the radial profiles of mass density and velocity dispersion, in order to investigate the dynamical evolution of tidally stripped subhaloes in more detail (e.g., Green et al., in prep.).

5 Summary and discussion

Cosmological $N$-body simulations are the prime tool used to study the hierarchical assembly of dark matter haloes. They reveal that virialized dark matter haloes have a universal density profile (e.g., Navarro et al., 1997), and that roughly 10 percent of their mass is bound up in distinct subhaloes (e.g., Ghigna et al., 1998; Gao et al., 2004; Giocoli et al., 2010). According to the same simulations, a large fraction of these subhaloes completely disrupt after a few orbital periods (Han et al., 2016; van den Bosch, 2017). It has recently been argued that the majority of this disruption is artificial (Peñarrubia et al., 2010; van den Bosch et al., 2018; van den Bosch & Ogiya, 2018), and thus that state-of-the-art cosmological simulations still suffer from an appreciable amount of ‘over-merging’. Most importantly, van den Bosch & Ogiya (2018) argued that this problem may go unnoticed in standard numerical ‘convergence’ tests. Hence, it is prudent that we consider alternative methods to predict the abundance and demographics of dark matter substructure, which is a potentially powerful Rosetta stone for deciphering the nature of dark matter.

One alternative to numerical simulations is a semi-analytical approach that combines halo merger trees, constructed using the framework of extended Press-Schechter theory (Bond et al., 1991), with a treatment of the tidal evolution of subhaloes as they orbit their host. These models are not hampered by discreteness issues or limiting force resolution responsible for artificial disruption. In addition, these models are far less CPU-intensive than actual $N$-body simulations, thus allowing for an extensive exploration of parameter space. Numerous models along this line have been constructed in the past (Taylor & Babul, 2001, 2004; van den Bosch et al., 2005; Peñarrubia & Benson, 2005; Zentner et al., 2005; Diemand et al., 2007; Kampakoglou & Benson, 2007; Gan et al., 2010; Pullen et al., 2014; Jiang & van den Bosch, 2016). Unfortunately, since we lack a purely analytical treatment of tidal stripping and heating, these models typically contain one or more ‘fudge’ parameters. These are tuned by requiring the model to reproduce the subhalo mass functions taken from cosmological $N$-body simulations. The obvious downside of this approach is that the models thereby inherit the over-merging problems of the simulations.

In order to overcome this dilemma, we need idealized simulations that are (i) well resolved and free from artificial disruption, and (ii) optimized to allow for calibration of semi-analytical treatments of tidal stripping and heating. This paper presents the DASH library, consisting of 2,177 idealized, high-resolution ($N=1,048,576$), collisionless $N$-body simulations of individual dark matter subhaloes orbiting in the potential of a static, analytical host halo. The simulations have sufficient mass and force resolution to overcome artificial disruption (i.e., they satisfy the numerical reliability criteria given by equations [9] and [11]), and sample the entire parameter space of orbital energies, orbital angular momenta, and halo concentrations relevant for dark matter substructure. All simulations adopt a host halo-to-subhalo mass ratio of ${\cal M}=M_{\rm vir,h}/M_{\rm vir,s}=1000$, for which dynamical friction, which is not accounted for in the DASH simulations, is negligible. Because of the scale-free nature of the tidal evolution of subhaloes (see §3.1), the DASH simulations are valid for any mass ratio large enough such that dynamical friction is negligble (i.e., {\cal M}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}100). Each simulation is evolved for a period of roughly 36 Gyr, during which the subhalo undergoes anywhere between 2 and 12 radial orbits. For each simulation, the DASH library, which is publicly available, contains simple text files that present, among others, the temporal evolution of the subhalo’s bound mass fraction, and the density and velocity dispersion profiles of the bound particles of the subhalo at 301 outputs equally spaced in time. The library also contains a Python code, trained on the DASH simulation data, that uses random forest regression to predict the bound mass fraction of subhaloes as a function of time for given halo concentrations and orbital parameters. This code, which is accurate at the 0.1 dex level, conveniently summarizes the main results from our large suite of simulations.

In a forthcoming paper (Green et al., in prep.), we use the DASH library to calibrate a new and improved semi-analytical model for the tidal evolution of subhaloes, which we will subsequently use in combination with accurate halo merger trees (e.g., Parkinson et al., 2008; Jiang & van den Bosch, 2014) to predict the subhalo mass function of CDM haloes, unhindered by artificial disruption. This will shed new light on the level of reliability of the subhalo demographics that have been extracted from cosmological $N$-body simulations.

Finally, we emphasize that although the parameter space covered by the DASH library is vast, it is by no means exhaustive. One obvious shortcoming, as discussed above, is that the DASH simulations are inadequate to describe major mergers with {\cal M}\lower 3.01385pt\hbox{;\buildrel<\over{\sim};}100. In those cases, dynamical friction due to the host, and self-friction due to tidally stripped material, cause the orbit of the subhalo to decay, exposing it to stronger tides. Another degree of freedom not covered here is the inner density slope of dark matter haloes. It is well known that observations of dwarf galaxies often suggest that their haloes have constant density cores, rather than the steep $r^{-1}$-cusps predicted by dark matter-only simulations (e.g., Burkert, 1995; Gentile et al., 2004; Oh et al., 2011; Hayashi & Chiba, 2015). Such cores can be created within the CDM paradigm by a variety of baryonic processes (e.g. El-Zant et al., 2001; Inoue & Saitoh, 2011; Pontzen & Governato, 2012; Ogiya & Mori, 2014), and have a dramatic impact on the tidal evolution of subhaloes (Peñarrubia et al., 2010; Errani et al., 2015; Ogiya, 2018). In addition, baryons modify the potentials of host- and subhaloes through the bulges and discs that they form at the halo centres, and these also strongly impact the tidal fields (Errani et al., 2017; Garrison-Kimmel et al., 2017). Finally, in the DASH simulations presented here, the host halo is assumed to be spherically symmetric, which allows us to completely specify each orbit with only two parameters (energy and angular momentum). Cosmological simulations, though, indicate that dark matter haloes are expected to be triaxial systems (e.g., Jing & Suto, 2002; Allgood et al., 2006; Hayashi et al., 2007), consistent with the shapes of the gravitational potentials of galaxies and clusters as inferred from a variety of observations (e.g., Oguri et al., 2005; Corless & King, 2007; Law & Majewski, 2010). Triaxial systems have a much richer variety of orbits, which is likely to impact the tidal evolution of subhaloes.

In the near future, we therefore anticipate augmenting the DASH library with a suite of simulations that probe some of this extended parameter space. These additional simulations will be particularly useful for informing the semi-analytical treatments mentioned above. Ideally, any such semi-analytical treatment should capture the actual physics of tidal stripping and heating, and should thus be able to correctly predict the tidal evolution of subhaloes in triaxial potentials, in the presence of orbital decay due to dynamical friction, or in the case where the potential of the host- and/or subhalo has been modified due to the impact of baryonic processes. It remains to be seen to what extent the models can meet this challenge, and it is our hope that the DASH library presented here, as well as its future extensions, will play an important role in this process.

Acknowledgements

We are grateful to the developers of scikit-learn for making their code publicly available. GO and OH acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 679145, project ‘COSMO-SIMS’). FvdB is supported by the National Aeronautics and Space Administration through Grant No. 17-ATP17-0028 issued as part of the Astrophysics Theory Program, and by the US National Science Foundation through grant AST 1516962.

Appendix A Simulation Invariance

As discussed in §2.1.1, the DASH simulations have been run in model units for which $G=M_{\rm vir,s}=r_{\rm vir,s}=1$. In these units, the subhalo has a crossing time $t_{\rm cross}\equiv r_{\rm vir,s}/V_{\rm vir,s}$ of unity. Converting to physical time units, simply requires multiplying the model time units by $1.44\,\Gamma\>{\rm Gyr}$ with $\Gamma=\Gamma(h,\Delta_{\rm vir})$ given by Equation (4). Hence, each DASH simulation is applicable to any combination of $\Delta_{\rm vir}$ and $H_{0}=100\,h\>{\rm km}\,{\rm s}^{-1}\,{\rm Mpc}^{-1}$, each of which with its own scaling between time in model units and time in physical units.

Changing $\Delta_{\rm vir}$ and/or $h$, while keeping the simulation parameters (${\cal M}$, $x_{\rm c}$, $\eta$, $c_{\rm h}$, $c_{\rm s}$) fixed, corresponds to changing the actual physical densities. For example, increasing $\Delta_{\rm vir}$ implies that the virial radius of the host halo decreases, such that a given value of $x_{\rm c}$ corresponds to a smaller physical radius, where the density of the host halo is larger. But, since the densities of the subhalo change similarly, and since the tidal evolution only depends on the ratio of densities (i.e., gravity is scale-free), the outcome of the simulation is entirely invariant to these changes in $\Delta_{\rm vir}$ and/or $h$.

However, when changing $\Delta_{\rm vir}$ and/or $h$ one can also re-scale the DASH simulations in another way, one that keeps the physical densities, and hence the mapping between model time and physical time, invariant. This scaling, however, requires a mapping between the parameter set $\{{\cal M},x_{\rm c},\eta,c_{\rm h},c_{\rm s}\}$ for the DASH simulation ($\Delta_{\rm vir}=200$ and $h=0.678$) and another parameter set $\{{\cal M}^{\prime},x^{\prime}_{\rm c},\eta^{\prime},c^{\prime}_{\rm h},c^{\prime}_{\rm s}\}$ corresponding to $\Delta^{\prime}_{\rm vir}$ and $h^{\prime}$. This scaling keeps the characteristic density and scale-radius, i.e., $\rho_{0}$ and $r_{\rm s}$ in equation (1), invariant. Consequently, a different value of $\Delta_{\rm vir}$, which corresponds to a different virial radius, now implies a different value for the mass and concentration parameter of the halo. And since the orbital radius remains invariant, also the dimensionless parameter $x_{\rm c}$ will change its value. Using that the characteristic density of an NFW halo with concentration parameter $c$ is given by

[TABLE]

it is straightforward to show that a DASH simulation for $\{{\cal M},c_{\rm h},c_{\rm s},x_{\rm c},\eta\}$ can be used to represent the evolution of a subhalo with parameters $\{{\cal M}^{\prime},c^{\prime}_{\rm h},c^{\prime}_{\rm s},x^{\prime}_{\rm c},\eta^{\prime}\}$ for any other combination of $\Delta^{\prime}_{\rm vir}$ and $\rho^{\prime}_{\rm crit}$ using the following mapping:

[TABLE]

Here $\Gamma=\Gamma(h^{\prime},\Delta^{\prime}_{\rm vir})$ is given by Equation (4), and ${\cal Q}={\cal Q}(c,\Delta^{\prime}_{\rm vir})$ is the root of

[TABLE]

Note that, whereas the $(\Delta_{\rm vir},h)$-dependent time-scaling is exact, this density-invariant mapping is only approximate. One of the reasons is that, because of dynamical friction, our simulations are only invariant to changes in ${\cal M}$ as long as the mass ratio M_{\rm vir,h}/M_{\rm vir,s}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}100. Hence, one can only use this mapping as long as ({\cal Q}_{\rm h}/{\cal Q}_{\rm s})^{3}\lower 3.01385pt\hbox{;\buildrel>\over{\sim};}0.1. In addition, the initial subhalo is only initialized out to a truncation radius $r_{\rm trunc}$ (see §2.1.1), which is equal to the virial radius, but only for $\Delta_{\rm vir}=200$ and $h=0.678$. For other values of $\Delta_{\rm vir}$ and $h$, the subhalo is effectively truncated at a radius $r_{\rm trunc}=r^{\prime}_{\rm vir}/{\cal Q}_{\rm s}$. As we have demonstrated in Paper II, the simulation outcome depends only very weakly on where exactly the initial subhalo is truncated, as long as it is outside of the initial tidal radius. The weak dependence mainly originates from ‘self-friction’ (see §3.1), which depends on the amount of mass that is stripped from the subhalo, which is larger if the initial truncation radius is larger. Hence, as long as the impact of self-friction is weak, which is almost always the case, we expect the density-invariant mapping of Equation (17) to be reasonably accurate. We have verified that this is indeed the case by running a few simulations for different values of $\Delta_{\rm vir}$, and comparing the resulting $f_{\rm b}(t)$ to predictions from the re-scaled DASH simulations based on the mapping of Equation (17). Except for the time period prior to the first pericentric passage, we find this mapping to be accurate at the few percent level. The reason why the mapping fails prior to the first-pericentric passage is simply that the mapping predicts that the initial subhalo starts out with a bound mass fraction $f^{\prime}_{\rm b}$ that is not equal to unity. However, after first pericentric passage the subhalo is basically stripped down to the same physical radius as in the fiducial case, and the mapping of Equation (17) is reliable.