The phase-space structure of tidally stripped halos

TL;DR

This paper introduces a new method for creating stable, finite, equilibrium models of dark matter halos that resemble tidally truncated profiles, useful for simulations and theoretical analysis.

Contribution

The authors develop a novel iterative truncation method for NFW halos, providing stable models and an analytic approximation for tidally stripped systems.

Findings

Models closely match tidally truncated density profiles

Models are highly stable for N-body simulations

Analytic approximation accurately describes tidally stripped halos

Abstract

We propose a new method for generating equilibrium models of spherical systems of collisionless particles that are finite in extent, but whose central regions resemble dark matter halos from cosmological simulations. This method involves iteratively removing unbound particles from a Navarro-Frenk-White profile truncated sharply at some radius. The resulting models are extremely stable, and thus provide a good starting point for N-body simulations of isolated halos. We provide a code to generate such models for NFW and a variety of other common density profiles. We then develop an analytic approximation to this truncated distribution function. Our method proceeds by analogy with the King model, truncating and shifting the original distribution function of an infinitely extended Navarro-Frenk-White profile in energy space. We show that the density profiles of our models closely resemble the tidally truncated density profiles seen previously in studies of satellite evolution. Pursuing this analogy further with a series of simulations of tidal mass loss, we find that our models provide a good approximation to the full distribution function of tidally stripped systems, thus allowing theoretically motivated phase-space calculations for such systems.

Figures

Tables

Table 1. Table 1: Summary of simulation parameters. Columns give (1) the simulation number (2) the mass ratio between the host and satellite halo (3) the virial radius of the host (4) the apocentric distance (5) the pericentric distance \color black (6) the tangential velocity at apocenter (7) the (radial) orbital period, (8) the circularity of the orbit, (9) the relative energy (defined as the energy divided by the energy of a circular orbit at the virial radius) (10) the radius of a circular orbit with the same energy divided by the virial radius.

Simulation	$M_{host}/M_{sat}$	$R_{vir}/r_s$	$r_a/r_s$	$r_p/r_s$	$v_a/v_{unit}$	$t_{orb}/t_{unit}$	${ϵ}_c$	$\eta$	$R_c/R_{vir}$
1	100	46.4	100	10	0.34	206.8	0.42	0.85	1.26
2	100	46.4	100	50	0.90	299.4	0.92	0.71	1.63
3	300	66.9	100	10	0.51	129.7	0.40	1.09	0.88
4	300	66.9	100	50	1.42	185.4	0.92	0.92	1.13
5	100	46.4	500	50	0.23	1778.5	0.47	0.24	6.14
6	300	66.9	25	10	1.50	31.48	0.82	2.29	0.26
7	50	36.8	80	5	0.19	201.6	0.30	0.86	1.24
8	50	36.8	90	15	0.37	259.7	0.58	0.76	1.48
9	10	21.5	40	10	0.30	196.8	0.71	0.88	1.19
10	10	21.5	25	10	0.42	123.2	0.85	1.14	0.82

Equations

$$\rho (r)=\frac{{\rho }_0}{r/r_s(1+r/r_s{)}^2},$$

$$\rho (r)={\rho }_{-2}\exp \left(-\frac{2}{\alpha }\left[{\left(\frac{r}{r_{-2}}\right)}^{\alpha }-1\right]\right),$$

$$\rho (r)=4\pi {\int }_0^{\Psi (r)}f(\mathcal{E})\sqrt{2(\Psi (r)-\mathcal{E})}d\mathcal{E}.$$

$$f(\mathcal{E})=\frac{1}{\sqrt{8}{\pi }^2}\left[{\int }_0^{\mathcal{E}}\frac{1}{\sqrt{\mathcal{E}-\Psi }}\frac{d^2\rho }{d{\Psi }^2}d\Psi +\frac{1}{\sqrt{\mathcal{E}}}{\left(\frac{d\rho }{d\Psi }\right)}_{\Psi =0}\right].$$

$$\rho (r)=\{\begin{array}{ll}\frac{{\rho }_0{r}_s^3}{r(r_s+r{)}^2}& \text{mbox}ifr<r_{vir}\\ \frac{{\rho }_0}{c(1+c{)}^2}{\left(\frac{r}{r_{vir}}\right)}^{ϵ}\exp \left(-\frac{r-r_{vir}}{r_d}\right)& \text{mbox}ifr>r_{vir}\end{array},$$

$$ϵ=-\frac{r_s+3r_{vir}}{r_s+r_{vir}}+\frac{r_{vir}}{r_d}.$$

$$t_{rel}(r)$$

$$t_{evap}(r)$$

$$F_{iso}(Z)=F_0\exp [Z],$$

$$F_{Wooley}(Z)=\{\begin{array}{ll}{F}_{iso}(Z-Z_t)& Z\ge Z_t\\ 0& Z\le Z_t.\end{array}$$

$$F_{King}(Z)=\{\begin{array}{ll}{F}_{iso}(Z-Z_t)-F_{iso}(0)& Z\ge Z_t\\ 0& Z\le Z_t.\end{array}$$

$$\begin{gathered} F_{NFW}(Z)=F_0{Z}^{3/2}(1-Z{)}^{-5/2}{\left(-\frac{\ln Z}{1-Z}\right)}^q \\ \times \exp (p_{1}Z+p_2{Z}^2+p_3{Z}^3+p_4{Z}^4), \end{gathered}$$

$$F(Z)=\{\begin{array}{ll}{F}_{NFW}(Z)-F_{NFW}(Z_t)& Z\ge Z_t\\ 0& Z\le Z_t.\end{array}$$

$$F(Z^{\prime })=\{\begin{array}{ll}{F}_{NFW}(Z^{\prime }+Z_t)-F_{NFW}(Z_t)& Z^{\prime }\ge 0\\ 0& Z^{\prime }\le 0,\end{array}$$

$$\frac{d^2\Psi }{d{r}^2}+\frac{2}{r}\frac{d\Psi }{dr}=-4\pi G\rho (\Psi ),$$

$$\frac{d^{2}P}{d{R}^2}+\frac{2}{R}\frac{dP}{dR}=-p(P)=-4\pi {\int }_0^{P}F(Z)\sqrt{2(P-Z)}dZ.$$

$$F_M(<r)=\frac{M(<r)}{M(r_{max})}.$$

$$f(\mathcal{E})=\frac{1}{\sqrt{8}{\pi }^2}\left[{\int }_0^{\mathcal{E}}\frac{1}{\sqrt{\mathcal{E}-\Psi }}\frac{d^2\rho }{d{\Psi }^2}d\Psi +\frac{1}{\sqrt{\mathcal{E}}}{\left(\frac{d\rho }{d\Psi }\right)}_{\Psi =0}\right].$$

$$\begin{array}{c}\frac{d^2\rho }{d{\Psi }^2}={\left(\frac{d\Psi }{dr}\right)}^{-2}\left[\frac{d^2\rho }{d{r}^2}-{\left(\frac{d\Psi }{dr}\right)}^{-1}\frac{d^2\Psi }{d{r}^2}\frac{d\rho }{dr}\right]\\ \mathrm{or},\mathrm{since}\frac{d\Psi }{dr}=-\frac{GM}{r^2}\\ \mathrm{and}\frac{d^2\Psi }{d{r}^2}=\frac{2GM}{r^3}-4\pi G\rho \\ \frac{d^2\rho }{d{\Psi }^2}=\left(\frac{r^4}{G^2{M}^2}\right)\left[\frac{d^2\rho }{d{r}^2}+\left(\frac{r^2}{GM}\right)\left[\frac{2GM}{r^3}-4\pi G\rho \right]\frac{d\rho }{dr}\right].\end{array}$$

$$F(<\mathcal{E})=\frac{{\int }_0^{\mathcal{E}}f(\mathcal{E})\sqrt{\Psi -\mathcal{E}}d\mathcal{E}}{{\int }_0^{\Psi }f(\mathcal{E})\sqrt{\Psi -\mathcal{E}}d\mathcal{E}}.$$

$$\mathcal{E}=\Psi -\frac{1}{2}{v}^2.$$

Full text

††pubyear: 2016††pagerange: The phase-space structure of tidally stripped halos–References

1 Introduction

On large scales, a number of independent and complementary tests, including the spectrum of fluctuations in the microwave background (e.g. Planck Collaboration et al., 2016), galaxy clustering (e.g. Alam et al., 2016), and weak gravitational lensing (e.g. Kitching et al., 2014), provide overwhelming support for the existence of dark matter. This component dominates over regular matter by a factor of \sim 5\ (Planck Collaboration et al., 2016), and collapsed, virialized dark matter halos are thought to be the site of all galaxy formation (see Frenk & White, 2012, for a review).

On smaller scales, observational tests including galaxy kinematics (e.g. Ouellette et al. in prep., Battaglia et al. 2013), satellite kinematics (e.g. Prada et al., 2003; Guo et al., 2012), and weak or strong gravitational lensing (e.g. Okabe et al., 2013; Umetsu et al., 2016), amongst others, are beginning to probe the structure of individual dark matter halos, placing direct constraints on their density profile and velocity structure, as well as central concentration, shape, and substructure. For the moment, however, most of our knowledge of these non-linear scales comes from numerical simulations, which have studied the formation and evolution of halos at ever increasing resolution (e.g. Diemand et al., 2007; Springel et al., 2008; Diemand et al., 2008; Stadel et al., 2009; Gao et al., 2012).

One of the main results of these simulations has been the observation that halos have a universal density profile (UDP), when averaged spherically. The classic approximation to this profile is from the work of Navarro et al. (1996, 1997, – NFW hereafter):

[TABLE]

where $\rho_{0}$ is a characteristic density and $r_{s}$ is the scale radius, corresponding to the point where the logarithmic slope is $d\ln\rho/d\ln r=-2$. This profile has no outer limit a priori; as $r$ goes to infinity the central potential remains finite but the mass diverges. In a cosmological context, halos are usually considered out to the virial radius, the radius within which they are in approximate virial equilibrium. In practice, the virial radius can be defined by one of several overdensity criteria, or it can be specified in terms of the concentration parameter $c\equiv r_{vir}/r_{s}$.

More recent work at higher resolution has determined that halos differ slightly but systematically from the original NFW fit, and that a better approximation is the Einasto profile (Navarro et al., 2004; Merritt et al., 2006; Gao et al., 2008):

[TABLE]

which contains an extra shape parameter $\alpha$ that seems to vary systematically with mass and redshift (Gao et al., 2008; Dutton & Macciò, 2014; Klypin et al., 2016). In this paper we will consider only the simpler NFW model, although our results are easily extended to an Einasto profile, and our initial condition code, described in the Appendix, includes the Einasto profile as an option.

Since the discovery of the universal density profile, there has been extensive work characterizing halo properties such as shape, concentration, spin and substructure (see Taylor 2011 for older references, or Klypin et al. 2016 for more recent references). Average values for properties such as concentration are now well determined as a function of mass, redshift and cosmology. On the other hand, given the complexity of halo growth through hierarchical merging in a cosmological context, idealized, isolated simulations of mergers between pairs of halos may be better suited to determining the physical mechanisms by which these average trends are established. For this work, it is convenient to be able to construct models whose density \colorblack is truncated at a finite radius \colorblack and/or whose mass converges to a finite value. Prior work (e.g. Moore et al., 2004; Kazantzidis et al., 2004) has used models with exponentially truncated NFW profiles, but there is no particular theoretical motivation for this choice. Our initial goal at the outset of this work was to provide a better motivated model for the initial conditions for simulations of isolated, truncated systems.

The standard theory of cosmological structure formation also predicts that halos should grow through repeated, hierarchical mergers, and that the cores of smaller merging systems should survive as tidally-stripped, self-bound substructure within galaxy, group and cluster halos. In previous work, Hayashi et al. (2003) studied the evolution of halo substructure by simulating the tidal stripping of a smaller satellite halo by a larger host halo. They found that, independent of orbit, the density profile of the satellite halo changed in a predictable way. This change could be described by an empirical model which only required one additional parameter, \colorblack equivalent to a tidal radius $r_{t}$. Though this paper did a thorough job of describing the profile of tidally stripped halos, \colorblack it did not determine the exact mechanisms behind the changes. Subsequent work by Kampakoglou & Benson (2007) showed that tidal forces acting on individual particles may explain some of the form of the \colorblack density profile of stripped systems, while Choi et al. (2009) showed that particles become unbound based more on their energy than on their angular momentum. Despite this work, the net effect of tidal forces on the full distribution function remains somewhat unclear.

Observational tests of these theoretical predictions based on lensing (e.g. Grillo et al., 2015; Jauzac et al., 2016) or internal kinematics (e.g. Ouellette et al. in prep.) are still in their infancy, but show promise given forthcoming datasets from very large surveys. Other substructure calculations, such as the boost factor for dark matter annihilation (e.g. Anderhalden & Diemand, 2013; Sánchez-Conde & Prada, 2014, and references therein), depend not only on the density profile, but also on the full distribution function (DF) of dark matter in stripped systems. For these applications, and in order to understand the mechanisms of tidal stripping fully, it would be useful to have a simple analytic model for the DF of a tidally stripped system.

In our experiments with \colorblack the truncation of distribution functions (DFs) in energy space, we have found a simple method for realizing spatially truncated halos whose central regions resemble NFW profiles. In this paper we describe the method, but also show that the truncated systems it produces closely resemble tidally stripped subhalos orbiting within a larger potential. Thus our model can be used to represent either stable, isolated systems with a finite extent (e.g.  for merger simulations) or tidally stripped systems (e.g. for substructure calculations).

The outline of the paper is as follows: in Section 2, we describe an iterative algorithm for generating the ICs for an isolated, \colorblack spatially truncated NFW profile. In Section 3, we compare this approach to an analytic model where the NFW DF is truncated in energy and shifted, analogously to a King model. In Section 4 we describe our simulations of tidal stripping of halos on various orbits, and in Section 5 we compare the properties of the stripped remnants to the analytic truncated models. Finally, in Section 6 we discuss the implications of our results and summarize our main conclusions.

2 Creating Truncated Initial Conditions

To simulate how a spherical dark matter halo evolves in isolation, we require a method for generating initial conditions (ICs) for particles which collectively form a stable self-gravitating system. Given an analytic expression for the density profile, we can integrate this to obtain the normalized cumulative mass distribution $M(<r)/M_{tot}$, and select particles randomly from this distribution, mapping enclosed mass fraction to radius. Clearly, profiles whose total mass $M_{tot}$ diverges at large radii need to be truncated in some way. The most common solution for cosmological halos (e.g. Kazantzidis et al., 2004; Moore et al., 2004; Kazantzidis et al., 2006; Peñarrubia et al., 2010) is to use the exponentially truncated NFW profile first introduced by Springel & White (1999), although this form is motivated more by mathematical convenience than by \colorblack any physical \colorblack argument.

There are two approaches to determining particle velocities (assumed in the simplest case to be spherically symmetric and isotropic in velocity space). We can either calculate the velocity dispersion at each radius from the Jeans equation, making the approximation that the velocity distribution is Maxwellian, (e.g. Hernquist, 1993), or, for a more accurate model, we can draw particle energies from the full DF (e.g. Kazantzidis et al., 2004) and convert these to velocities. We review the latter approach below, before introducing and testing our truncation method.

2.1 Generating Initial Conditions from the Distribution Function

We will briefly review the main properties of the DF; a more detailed explanation can be found in Binney & Tremaine (1987). The DF, $f$, describes the mass per phase-space volume. \colorblack For a spherical, isotropic, self-gravitating system, it can be written as a function of a single variable, say a relative energy $\mathcal{E}$. Given a potential $\Phi$, we will define the \colorblack (potential and total) relative energies as $\Psi=-\Phi+\Phi_{0}$ and $\mathcal{E}=\Psi-v^{2}/2$. The parameter $\Phi_{0}$ is free, and is generally set to the value of the potential at the outer boundary of the system, such that the $f>0$ only when $\mathcal{E}>0$, while $f=0$ when $\mathcal{E}\leq 0$.

A relationship between the density profile of a system and the DF can be found by integrating $f(\mathcal{E})$ over all velocities. For a spherically symmetric system:

[TABLE]

Equation (3) can be inverted (Eddington, 1916), and thus the DF can be expressed in terms of the density:

[TABLE]

A particle at radius $r$ then has energy $\mathcal{E}\in(0,\Psi(r))$ with a probability proportional to $f(\mathcal{E})\sqrt{\Psi(r)-\mathcal{E}}$.

\color

black Given a density profile $\rho$, the relative potential $\Psi(r)$ can be calculated using Poisson’s equation, and the derivatives $d\rho/d\Psi$ and $d^{2}\rho/d\Psi^{2}$ can be evaluated analytically or numerically. The distribution function $f(\mathcal{E})$ can than be calculated for a given relative energy $\mathcal{E}$ using Equation (4). To generate a model halo, a radius and a relative energy are selected at random for each particle, in such a way as to reproduce the correct density profile and distribution function. Once \colorblack the radius and relative energy have been assigned, the norm of the velocity of each particle can be calculated as \colorblack $v=\sqrt{2(\Psi(r)-\mathcal{E})}$. Finally, 3D position and velocity \colorblack components can be chosen at random assuming (spatial) spherical symmetry, and isotropy in velocity space, respectively. \colorblack Although these techniques are well known, for convenience we review them in Appendix A.

2.2 Truncating the NFW Profile at a Finite Radius

The exponentially truncated NFW profile was first described in Springel & White (1999), and is identical to the NFW profile within the virial radius, $r_{vir}$, but is truncated exponentially beyond that. How fast this decay occurs depends on the parameter $r_{d}$:

[TABLE]

where the constants $\rho_{0}$ and $r_{s}$ are as in Equation (1). Additionally, to ensure that the logarithmic slope at $r_{vir}$ is continuous, there is the constraint:

[TABLE]

Not only does this solution have little physical motivation, \colorblack however, but it is also discontinuous in the second derivative of \colorblack $\rho(r)$. Since Equation (4) depends on this second derivative through the ${d^{2}\rho}/{d\Psi^{2}}$ term, the DF may not be monotonically increasing for certain choices of $r_{d}$, resulting in unphysical behaviour.

\color

black Instead, we propose a different modification to the NFW profile. Our solution is to consider particles as if they were sampled from an infinitely extended NFW profile, but only choose those within some radius $r_{cut}$. Particle energies are initially assigned from the \colorblack DF corresponding to the full (infinitely extended) NFW \colorblack profile, as described in Section 2.1. We then iteratively remove any unbound particles, using at each step of the iteration the potential defined by the set of particles remaining.111\colorblackNote Choi et al. (2007) propose an alternative, but more complicated method, truncating the initial density profile, calculating a distribution function via Eddington inversion, recalculating the resulting density profile, and iterating over these steps until convergence. For a halo with $r_{cut}=10\,r_{s}$, this process converges within 10 iterations or fewer, leaving approximately 65% of the mass within $r_{cut}$ bound. The final result is a system with a density profile that matches NFW at small radii, but drops off more steeply at large radii, reaching zero density at $r_{cut}$. While this truncated profile is instantaneously self-bound, its long term stability is unclear. We will consider this point in the next section.

2.3 Stability Tests

Before we test the stability of our truncated profile, we need to estimate the possible contribution from collisional effects. These can complicate the interpretation of any N-body simulation. Our method for determining particle ICs assumes that the particles exist in a smooth continuous potential; each particle can be considered a Monte Carlo sampling of this potential. However, individual particle-particle interactions lead to a gradual loss of the initial structure, most noticeably in the dense center of the halo.

There are two time scales on which collisional effects are important. The first, the relaxation time scale, $t_{rel}$, corresponds to the time which it takes a typical particle’s velocity to change by an order of itself. The second, slower timescale is the evaporation time, $t_{evap}$, the amount of time it takes a typical particle to reach escape speed, and thus ‘evaporate’ from the system. Following the arguments from Binney & Tremaine (1987), we calculate:

[TABLE]

where $N(<r)$ is the number of particles within radius $r$. Using these formulae, we can define a ‘(central) relaxation radius’, $r_{relax}(t)$, and a ‘(central) evaporation radius’, $r_{evap}(t)$, such that $t_{rel}(r_{rel})=t$ and $t_{evap}(r_{evap})=t$ respectively, within which relaxation and evaporation are important.

To verify the stability of our initial satellite halo, we evolve it in isolation with the $N$-body code Gadget-2 (Springel, 2005). We generate an initial NFW system with \colorblack$2\times 10^{6}$ particles within a radius $r_{cut}=10r_{s}$. After iteratively removing unbound particles (as described in Section 2.2), the final number of particles in the truncated system is \colorblack$N=1286991$. We then assign each particle a mass $m$ such that the total satellite has mass $M_{sat}$. We simulate the evolution of the system using the softening length proposed in van Kampen (2000), $\epsilon=0.5r_{h}N^{-1/3}$, where $r_{h}$ is the half-mass radius. For our ICs, \colorblack $\epsilon=0.01r_{s}$.

Fig. 1 shows the evolution of the density profile with time. The outer part of the profile appears to be completely stable, with any systematic changes invisible on the scale of the figure. The only visible change is in the region interior to $r_{evap}$. \colorblackThis is consistent with the results \colorblack of Hayashi et al. (2003), who found \colorblack that the internal structure evolves at radii where the evaporation time is close to $t$. \colorblack (We have also confirmed the predicted scaling of $r_{evap}$ with particle number using lower resolution simulations.) We conclude that this central change is a collisional effect due to finite resolution, and that our ICs are otherwise extremely stable.

One disadvantage of our iterative unbinding method is that the resulting truncated DF does not have a simple expression, even for cases where the original, un-truncated DF is simple. Thus, in the next section we will consider a slightly different analytic approach, and show that it produces remarkably similar results.

3 An Analytic Model for the Truncated Distribution Function

Our goal is to derive an analytic approximation for the DF of a truncated system with a NFW or UDP-like central density profile. We will take an approach similar to the derivation of the King model (King, 1966), \colorblack truncating and lowering the DF in energy space, as described below. \colorblack For the NFW profile, this approach has been proposed previously by Widrow & Dubinski (2005), who showed that it leads to sharply truncated density profiles (see their Figure 1).

3.1 Review of the King Model

The King model is a lowered isothermal model based on the (infinitely extended) isothermal sphere, but with a lowered relative energy $\mathcal{E}$. Additionally, it subtracts a constant term to ensure that the DF is continuous at this truncation energy. This last step results in a more extended profile (Hunter, 1977), and presumably increases the stability of the model (e.g. Guo & Li, 2008).

As previously, we will define relative energies $\Psi=-\Phi+\Phi_{0}$ and $\mathcal{E}=\Psi-v^{2}/2$. Suppose initially we choose $\Phi_{0}=0$. The DF for an isothermal sphere of velocity dispersion $\sigma$ can be written:

[TABLE]

where $Z=\mathcal{E}/\sigma^{2}$ \colorblack is a dimensionless, scaled version of the relative energy. As with the NFW profile, this DF corresponds to a profile that extends to infinity and has infinite mass.

A solution to the infinite mass problem of the isothermal sphere is to lower the relative energy, letting $Z\rightarrow Z-Z_{t}$ (Woolley, 1954). This produces a truncated DF of the form:

[TABLE]

If \colorblack one uses the freedom of $\Phi_{0}$ to express \colorblack $F$ in terms of $Z^{\prime}=Z-Z_{t}$, the DF has the same form as \colorblack that of the infinitely extended isothermal sphere, but \colorblack the density now drops to zero at a finite radius.

This form of the DF introduces a new complication, \colorblack however, as it is now discontinuous at $F(Z_{t})$. A solution to this \colorblack problem is to subtract a constant term. DFs of this form are known as King models (Michie & Bodenheimer, 1963; King, 1966):

[TABLE]

3.2 Energy Truncation of a NFW Distribution Function

For mathematical convenience, we introduce the following dimensionless variables; $R=r/r_{s}$, $p=\rho/\rho_{0}$, $P=\Psi/(4\pi G\rho_{0}r_{s}^{2})$, $Z=\mathcal{E}/(4\pi G\rho_{0}r_{s}^{2})$ and $F=(4\pi G)^{3/2}r_{s}^{3}\rho_{0}^{1/2}f$, \colorblack where $\Psi(r)=-\Phi(r)+\Phi_{0}$ and $\mathcal{E}=\Psi(r)-v^{2}/2$, as above. Here, the energies have been normalized by the magnitude of the central potential of a NFW profile $|\Phi_{NFW}(r=0)|=4\pi G\rho_{0}r_{s}^{2}$.

The DF for the NFW profile can be determined numerically from Equation (4). However, in practice, we use the analytic approximation proposed by Widrow (2000):

[TABLE]

where $q=-2.7419$, $p_{1}=0.3620$, $p_{2}=-0.5639$, $p_{3}=-0.0859$, $p_{4}=-0.4912$ and $F_{0}=0.091968$. \colorblack We find that this approximation agrees with the numerically calculated DF to within 2%.

We wish to truncate the NFW DF at some truncation energy $\mathcal{E}_{t}$, or in dimensionless form $Z_{t}=\mathcal{E}_{t}/(4\pi G\rho_{0}r_{s}^{2})$. \colorblackThe modified DF is then given by:

[TABLE]

We can exploit the freedom of $\Phi_{0}$, picking a new value such that the relative energy is zero on the boundary. Given this new value of $\Phi_{0}$, and denoting the new relative energy variable $Z^{\prime}=Z-Z_{t}$, the DF then becomes:

[TABLE]

where $Z_{t}$ the truncation energy defined using the original value of $\Phi_{0}$.

Note that this derivation is slightly different \colorblack from that of the lowered isothermal sphere model outlined in Section 3.1. Comparing Equations (10) and (13), we find that applying the latter to the isothermal sphere will recover the King model, but multiplied by a constant term $e^{Z_{t}}$.

It can be shown that the \colorblackrelative energy of the energy-truncated NFW profile has a maximum \colorblack value $Z=1-Z_{t}$; this corresponds to the central relative potential $P(0)$ of the truncated system. The final system \colorblack has a finite radius, $r_{t}$. Increasing the truncation energy will decrease the central potential, the truncation radius, and the total mass of the system.

We need two more equations to describe this system. The density profile can be recovered from Equation (3), and the relative potential can be determined from Poisson’s equation:

[TABLE]

or, re-expressed in the dimensionless parameters:

[TABLE]

The initial conditions are $P(0)=1-Z_{t}$ and $dP(0)/dR=0$. Equation (15) can be integrated numerically until $P(R_{t})=0$.

Putting all this together, the density of the truncated halo at radius ${\color{black}R}$ can be determined using the following steps:

Specify the dimensionless truncation energy $Z_{t}\in(0,1)$. 2. 2.

\color

blackGiven the distribution function, calculate the potential at radius ${\color{black}R}$ by numerically integrating Poisson’s equation (Equation (15)) with initial conditions $P(0)=1-Z_{t}$, $dP(0)/dR=0$. 3. 3.

Find the density by numerically integrating Equation (3).

The potential $\Phi(r)=-\Psi(r)+\Phi_{0}$ can be recovered once the truncation radius has been determined\colorblack, since $\Phi_{0}=-GM(r_{t})/r_{t}$.

Finally, we note that although we have discussed the energy-truncation method specifically for an NFW profile, it can be used for any \colorblack density profile with a known or calculable DF, including an Einasto profile\colorblack, and could also be extended to the various theoretically motivated models of the intrinsic halo DF (e.g. Hjorth & Williams, 2010; Pontzen & Governato, 2013; Beraldo e Silva et al., 2014). For cored density profiles, energy truncation will generally reduce the central density of the system significantly as the outer radius decreases, as discussed by Widrow & Dubinski (2005) and Peñarrubia et al. (2010).

3.3 Properties of the Truncated Model

As we increase the truncation energy, \colorblack the mass and extent of the NFW profile will decrease progressively. Fig. 2 shows how the density profile changes as a function of the dimensionless truncation energy $Z_{t}$, relative to the original NFW profile. Note that by definition, the density drops to zero outside the truncation radius (indicated by the vertical dashed lines).

The relationships between the truncation energy $Z_{t}$ and the total bound mass and truncation radius are shown in Fig. 3. Both the mass and the truncation radius are smooth, decreasing functions of the truncation energy, as expected. \colorblack The polynomial fits to the data shown (and valid over the range of the plot) are $M_{f}\equiv M(r<r_{t})/M_{NFW}(r<r_{t})=0.35Z_{t}^{2}-1.14Z_{t}+0.83$ and $r_{f}\equiv\log_{10}(r_{t}/r_{s})=-3.70Z_{t}^{3}+5.93Z_{t}^{2}-4.56Z_{t}+2.01$. We can also fit the inverse relations: $Z_{t}=0.6M_{f}^{2}-1.7M_{f}+1.02$ and $Z_{t}=0.2r_{f}^{3}-0.4r_{f}^{2}-0.39r_{f}+0.94$.

Finally, we return to our original goal, to establish \colorblack an analytic approximation to the DF produced by the iterative unbinding procedure we introduced in Section 2. In Fig. 4, we show how the ICs derived in Section 2 compare to the \colorblack analytic energy-truncated NFW model developed in this section. The histogram shows the ICs, the lower \colorblack dashed curve shows the result of truncating the NFW DF at the energy $\mathcal{E}_{t}$ (i.e. $f(\mathcal{E})=f_{NFW}(\mathcal{E}+\mathcal{E}_{t})$), while the upper \colorblack dotted curve shows the full model, including the shift to make $f$ continuous at zero: $f(\mathcal{E})=f_{NFW}(\mathcal{E}+\mathcal{E}_{t})-f_{NFW}(\mathcal{E}_{t})$. As before, the \colorblack solid line shows the original untruncated profile. Both \colorblack models provide a good match to the density profile of the ICs. As with the King model, subtracting the constant term $f_{NFW}(\mathcal{E}_{t})$ produces a more extended profile. We will adopt this version as our final analytic model for the DF, on the assumption that it is slightly more stable than the model where $f$ is discontinuous at zero. Overall, these results suggest that an energy-truncated DF comes close to describing \colorblack ICs obtained using the method outlined in Section 2. In principle, ICs could therefore be generated directly from the energy-truncated DF, although in practice our code implements the iterative unbinding procedure.

4 Simulating Tidally Stripped Halos

The very similar truncated profiles derived either by iterative unbinding in Section 2, or analytically in Section 3, were designed to represent isolated systems of finite mass and radial extent. We note, however, that they \colorblack also look very similar to the \colorblack density profiles of tidally stripped halos orbiting within the potential of a larger system (e.g. Hayashi et al., 2003). To pursue this analogy, we will compare our truncated models \colorblack directly to tidally stripped halos taken from simulations \colorblack of satellite mass loss. We describe the simulations and basic analysis \colorblack below, and then compare the simulated systems to our models in Section 5.

4.1 Simulation Parameters

Simulations of a smaller ‘satellite’ halo orbiting within the potential of a larger static ‘host’ halo were performed using the $N$-body code Gadget-2 (Springel, 2005). This code was modified to contain a fixed background potential corresponding to a host halo with an NFW profile. We use the mass and scale radius of the satellite halo ($M_{sat}$ and $r_{s}$) as the mass and distance units. Time is given in units $t_{unit}=\sqrt{r_{s}^{3}/GM_{sat}}$, and velocity in units of $v_{unit}=\sqrt{GM_{sat}/r_{s}}$. The host and satellite halos were assumed to have the same initial density within their outer, or ‘virial’ radii, as would be the case for a merger between two cosmological halos at a fixed redshift. We generated initial conditions for the satellite using our iterative unbinding algorithm with $r_{cut}=r_{vir}=10\,r_{s}$, as described in Section 2. The virial radius of the main halo scales as $(M_{host}/M_{sat})^{1/3}$, while the scale radius of the host was set assuming $c_{host}=10$.

\color

black

The orbital parameters of infalling satellite halos have been studied extensively in simulations (Tormen, 1997; Ghigna et al., 1998; Vitvitska et al., 2002; Benson, 2005; Wang et al., 2005; Zentner et al., 2005; Khochfar & Burkert, 2006; Wetzel, 2011; Jiang et al., 2015). The parameters used in this paper (see Table 1) cover the full range of energy and angular momentum expected for cosmological mergers (e.g. Jiang et al., 2015), with the exception of Simulations 5 and 6 which have unusually low and high energies, respectively. This allows us to test our tidal-stripping model not only for cosmological orbits, but also for a few more extreme cases.

We considered four different host/satellite mass ratios, $M_{host}/M_{sat}=300,\,100,\,50$ and $10$. Since the host halo was modelled as a fixed background potential, the satellite is not subject to dynamical friction in our simulations, and we expect its specific energy and angular momentum to be roughly conserved, even as it loses mass. The assumption that the host halo is static becomes less physically valid for low mass ratios, where dynamical friction plays a larger role. Nonetheless, we run a few cases at smaller mass ratios to test the effect of a larger satellite on the evolution of the density profile. At large mass ratios, where the satellite is small compared to the scale of the background potential, we expect satellite evolution to become independent of mass ratio. We chose a limiting value of $M_{host}/M_{sat}=300$ for practical reasons, as the ratio of the simulation time step to the orbital period is becoming very long at this point.

4.2 Locating the Satellite Remnant

We identified the satellite remnant at any given time using a method similar to the one outlined in Tormen et al. (1997, 1998). This is also the method used in Hayashi et al. (2003). There are two steps to this method; first the highest density peak of the particles is located approximately by iteratively decreasing the radius $R$ of a sphere, and re-centering it at each step on the center of mass of the particles contained within the sphere. We decreased the sphere by $0.9R$ on each iteration, and repeated until there were fewer than 100 particles within the sphere. We expect our results to be insensitive to this choice of final particle number, as discussed in Tormen et al. (1997). The velocity of the satellite frame was then calculated as the average velocity of all satellite particles within a sphere of radius $r_{cut}$ \colorblack (the original truncation radius of the satellite) centered on the highest density peak. The second step of this process was to identify which subset of particles was self-bound in this frame. This was calculated by iteratively removing unbound particles in the rest frame of the satellite until the algorithm converged.

Finally, we found that in a few cases where the algorithm had trouble locating a self-bound remnant (particularly at late times in Simulation \colorblack 6), we were able to improve the algorithm by first approximating the location of the satellite remnant by integrating \colorblack forward from the previous snapshot (assuming the satellite was a point mass orbiting in the potential of the host), and then only considering particles within $2r_{cut}$ of this predicted location.

Fig. 5 shows how the recovered satellite mass decreases with time for the ten simulations.

Given a bound remnant and associated satellite reference frame in each step, we then fit the remnant’s profile to the analytic model from Section 3, as described below.

5 Results

5.1 Fitting Criteria

To compare our tidally stripped satellites to the energy truncated DF model, we need a way of determining the free parameter, \colorblack $Z_{t}$. We do this by normalizing our DF models, such that their scale radii and central densities before truncation match those of the original NFW profile from which the ICs for the simulations were derived by the iterative method. This accounts for the fact that the initial conditions in the simulations are already truncated at a finite radius. $Z_{t}$ is then chosen so that both the simulations and the models have the same mass. We compare the analytic model to the simulations at apocenter on successive orbits, since these are the times when we expect the satellite to be closest to equilibrium.

5.2 Density Profiles, Enclosed Mass and Circular Velocity Profiles

The top panels of Fig. 6 show a comparison of the satellite density profiles from simulations (points) with the energy-truncated model (solid curves), where $Z_{t}$ has been fixed as described above. The thick dashed curve shows an untruncated NFW profile. \colorblack We demonstrate how the profile changes in time for Simulations 4 and 3 in the two left-most panels, as well as how the different simulations compare in the two right-most panels. Note that Simulation 4 and Simulations 2, 5, 8 and 9 are the orbits with the slowest mass loss, and therefore the ones we might expect to be the most successfully described by our model. \colorblackThe bottom panels of Fig. 6 show the relative residuals in density, ($\rho_{sim}/\rho_{model}-1$). These are generally less than 10% \colorblack for the orbits with slower mass loss, or less than 20% for orbits with faster mass loss. They are largest for Simulation 6, which loses mass extremely rapidly (see Fig. 5), or at late times for Simulation 3. The residuals are generally slightly larger at large radii, where an excess relative to the model is also visible in the top panels. Previous authors (e.g. Peñarrubia et al., 2008, 2009; Peñarrubia et al., 2010, \colorblack and earlier references therein) have noted a \colorblack distinct population of particles in the outer parts of tidally stripped systems, at radii too large to have crossed the satellite in the time since last pericentric passage. This may account for some of the excess seen here.

Fig. 7 and Fig. 8 show a comparison of the enclosed mass and circular velocity profiles, respectively. Linestyles are as in Fig. 6. Here too, residuals are generally at the 10–20% level or less. \colorblackThe largest deviations occur in Simulation 6, which loses mass the fastest. It should also be noted that Simulation 6 corresponds to an extreme, uncosmological orbit, with unusually large energy and extremely rapid mass loss, so it is less representative of realistic halo mergers. We \colorblack also note that our mass profiles and circular velocity profiles are constrained by the condition that the stripped mass fraction in the simulations matches the corresponding fraction in the models; a slightly different choice of $Z_{t}$ would have improved the agreement at some radii (e.g. near the peak of the circular velocity curve in Simulation 1), at the expense of a slightly worse fit close to $r_{cut}$.

\color

black In general, the radius at which evaporation becomes important is not visible on the scales of Figs. 6 – 8 . The exception is Simulation 5, which has a very long orbital time scale. Thus, the deficit of mass at small radii for this simulation can be attributed to collisional effects.

5.3 Distribution Function and Moments of the Phase-Space Distribution

Finally, we can attempt to compare the full DF of our simulated systems to the analytic energy-truncated model. We note \colorblack that the construction of the DF requires a frame in which to define velocities, and a set of particles to integrate over when calculating the potential. Thus the ‘distribution function’ of a subsystem within a larger halo is a slightly problematic concept, relative to the usual definition for an isolated system. Here, to be concrete, we define velocities in the mean center-of-mass frame of the self-bound remnant, and calculate potential energies summing only over those particles that are bound.

Given this convention, the DFs \colorblackfor the first four orbits of Simulation 4 are shown in Fig. 9. Each point represents the number of particles in a bin of normalized \colorblack relative energy $Z$. This was found by binning the particles in $500$ equally sized bins in both radius and velocity. The phase-space density of each bin was then calculated as the number of particles divided by $16\pi^{2}r^{2}v^{2}drdv$. \colorblack As usual, the relative energy is $\mathcal{E}=\Psi(r)-v^{2}/2$, where $\Psi(r)=-\Phi(r)+\Phi(r_{max})$, and $r_{max}$ is the radius of the remnant. Only bins with phase-space volume greater than $10^{-5}$ were plotted, to avoid numerical errors resulting from dividing by small numbers.

The solid lines show the prediction of the analytic model from Section 3 for comparison, while the dashed line shows the original (untruncated) NFW DF. Generally speaking, there is good agreement between the tidally stripped DFs and the analytic models for low \colorblackrelative energies, though it is difficult to compare the simulation to the model for large values of $Z$ (corresponding to particles at small radii and/or low energies), since the phase-space volume becomes very small in this limit.

\color

black In Fig. 10 we also show the evolution of the pseudo-phase-space density for Simulations 3 and 4 (top panel). As first noted by Taylor & Navarro (2001), the spherically averaged pseudo-phase-space density, $\rho/\sigma^{3}$, of isotropic, cosmological halos appears to follow a simple power law as a function of radius. The dashed line shows this power-law for an infinitely extended NFW profile. For the truncated analytic models (solid lines), we can calculate the velocity dispersion from the probability distribution function for velocities, $P(v)\propto f(\mathcal{E})v^{2}/\rho(r)$. The analytic models also follow a power law out to around the truncation radius, but with a flatter slope. The relative increase in pseudo-phase-space density at large radii is expected, since energy truncation reduces the number of particles with large velocities, and thus the velocity dispersion, in these regions. A similar flattening of the slope has also been seen in tidally truncated subhalos from self-consistent cosmological simulations (Vera-Ciro et al., 2014).

Interestingly, while the simulation results (points) match the analytic models extremely well at small radii, they deviate from them systematically at large radii. The bottom panel shows the likely reason for this discrepancy: the tidally stripped simulations are not isotropic in their outer regions. Plotting the anisotropy parameter $\beta=1-(\sigma_{\theta}^{2}+\sigma_{\phi}^{2})/2\sigma_{r}^{2}$ as a function of radius, we see that the analytic model and simulation results in the top panel begin to differ at radii where $\beta$ is becoming significantly different from zero, suggesting that our analytic model fails to match the outer parts of the simulations primarily because the assumption of isotropy is no longer valid there.

To get a better sense of where the simulations differ from the models in \colorblack the full phase space, we can also bin the DF in radius and (total) velocity. Fig. 11 shows this \colorblack 2-D distribution for the first three orbits of Simulation 4 (top three panels), the corresponding analytic models (middle panels), and the differences between the two (lower panels). Generally speaking, there is good agreement between the models and the simulations \colorblack in this projection of the full phase space. The tidally stripped satellites show a slight excess of particles at large radii and large velocities, and this excess may even have a caustic-like structure to it. It seems likely that this is related to the excess component at large radii \colorblack noted by earlier authors and discussed in the previous section. Peñarrubia et al. (2009) found that this outer material, while still bound, is on its way out of the system and is mostly lost on the next orbit. Such non-equilibrium effects would not be captured by our analytic truncation model. There is also a slight deficit of particles at small radii and low velocities. Here, relaxation may play a role in scattering particle energies, and depleting the lowest-energy orbits.

6 Discussion and Conclusions

Dark matter halos play a central role in the current picture of cosmological structure formation, and their properties and evolution determine many of the broad trends in the properties of galaxies, galaxy groups and galaxy clusters. Our understanding of halo formation and evolution is still limited, however; many basic questions, such as the origin of the universal density profile or the concentration-mass relation, have not yet been answered fully. Progress in understanding halo properties will come from a combination of fully realistic, large-volume cosmological simulations, together with simpler, idealized simulations of individual halos. Our original goal, in this work, was to develop a better algorithm for generating initial conditions for the latter.

We have found through experimentation that two different approaches yield similar results: either truncating an NFW profile abruptly at some radius $r_{cut}$, and then iteratively removing unbound particles until convergence is reached, or using an analytic NFW distribution function truncated at some energy and shifted such that the resulting distribution function is continuous at $\mathcal{E}=0$. The latter technique, inspired by the King model (King, 1966) \colorblack and proposed previously by Widrow & Dubinski (2005), produces almost the same density profile as the former, and allows us to construct simple analytic models for spatially finite systems whose central regions resemble NFW profiles. For generating initial conditions in practice, we have used the first method; our tests show that these models are extremely stable, and thus well suited to the study of isolated systems.

For either of these solutions, the density profile drops off steeply, reaching zero at a finite truncation radius. This behaviour is familiar from numerical studies of tidally stripped halos, starting with (Hayashi et al., 2003). Pursuing the suggestion of a connection between the two, we have simulated the tidal evolution of satellite subhalos in the potential of a larger system. We find that the resulting tidally stripped remnants match our theoretical models at the 10–20% level or better in density profile, enclosed mass and/or circular velocity profile. In phase space they also look similar, particularly when we plot phase-space density \colorblack projected in terms of radius and \colorblack total velocity. Some of the minor differences may relate to the presence of unrelaxed material close to the outer edge of the system, as discussed by Peñarrubia et al. (2009), or to relaxation effects near its centre. The origin and full significance of these differences will require further work to clarify. Nonetheless, overall our analytic model provides a good first approximation to the detailed distribution function of tidally stripped systems.

This result is actually slightly surprising. Tidal mass loss is a complex process, and has been investigated by many authors previously, using different approaches and approximations. Taylor & Babul (2001), for instance, developed a 1-D model for tidal stripping, calculating the average effects of tidal heating and expansion, as well as adiabatic cooling, within spherical shells, and removing the mass outside an instantaneous tidal limit over a timescale equal to the instantaneous (angular) orbital period, $t_{orb}=2\pi/\omega$. A limitation of this model, in the present context, is that it did not specify the actual density profile of the satellite at any one time, but only the net effect of heating and cooling on mass loss. Benson et al. (2002) considered a similar model, but with several corrections, notably using the shorter of the angular orbital period and the radial infall time $R/v_{R}$ as the timescale for mass loss. They found this gives better agreement with mass loss rates from simulations for very radial orbits. Analytic models of a type similar have also been developed by several other groups (e.g. Zentner & Bullock, 2003; Oguri & Lee, 2004; van den Bosch et al., 2005; Zentner et al., 2005).

Kampakoglou & Benson (2007) considered the problem of mass loss in more detail, calculating the effects of heating \colorblack on circular orbits more precisely using particles to sample the distribution function in energy and angular momentum, as well as for particle orbits of different inclination. They found that by successively heating and removing unbound particles, they predicted a profile similar to the one seen in their simulations, at least in the limit of weak tidal fields. This may provide the first hint that modelling mass loss in energy space can naturally explain the profile of tidally stripped cosmological halos.

Given that they were trying to establish simpler, semi-analytic prescriptions for tidal mass loss, none of these approaches considered the full complexity of the problem. In reality, the tidal boundary is non-spherical even for a spherical system on a circular orbit; for a general orbit it is also time-varying, and the timescales for relaxation and mass removal in the outer parts of the satellite are unfortunately very close to the orbital timescales. Furthermore, real halos are usually triaxial, and may have their spin partly coupled to their orbital angular momentum. Thus a general, predictive and highly accurate model for tidal mass loss seems a distant prospect.

Nonetheless, some simple patterns do emerge. Hayashi et al. (2003) first pointed out that all tidally truncated NFW profiles look the same, and can be fit by a single additional parameter; they showed that mass loss in their simulations is essentially a 1-dimensional sequence, a result later confirmed for a wider range of profiles by Peñarrubia et al. (2010). Our results now connect this stripping sequence to a simple cutoff and shift in the underlying distribution function, parameterized by an increasing truncation energy $\mathcal{E}_{t}$. Exactly how heating, relaxation, and unbinding conspire to produce this simple result is unclear \colorblack. A hint comes from the work of Choi et al. (2009), who show that individual particles in tidally stripped systems do or don’t become unbound based primarily on their energy, rather than their angular momentum. Thus, a simple cut in energy may provide a fairly accurate description of the mass-loss process.

Producing a more detailed description of tidal mass loss will be the focus of future work. In the interim, the models presented here provide a new, physically motivated way of generating initial conditions for isolated halos or tidally truncated subhalos. \colorblack Applications include studies of tidal stream formation (e.g. Amorisco, 2015, and references therein), disk heating (e.g. Moetazedian & Just, 2016, and references therein), mergers (e.g. Carucci et al., 2014, and references therein), and dwarf galaxy evolution (e.g. Tomozeiu et al., 2016, and references therein). A python code for generating ICs \colorblack using our iterative method is available on-line. Its use is described in the appendix.

\color

black Finally, we note several limitations to this work. We have studied the tidal evolution of idealized, spherical, isotropic halos, with an NFW density profile and corresponding DF. Halos in self-consistent cosmological simulations differ from this idealized case in a number of ways. It is tempting to look to self-consistent simulations to try to understand tidal mass loss in a more realistic situation; Springel et al. (2008), for instance, show the density profiles of tidally truncated subhalos resolved in the Aquarius simulations with $10^{5}-10^{6}$ particles or more. Broadly speaking they resemble our models, with a profile similar to field halos (i.e. close to an NFW or Einasto profile) in the inner parts, and truncated abruptly in the outer parts. Unfortunately, in self-consistent simulations the exact profile and distribution function of a subhalo will depend on the group-finder used. Group-finders use various different criteria for defining the boundary of a subhalo, producing slightly different results (e.g. Muldrew et al. 2011; see Onions et al. 2012 for a general review). Establishing precisely which particles are or aren’t associated with a subhalo, in a region dominated by the background potential of the main system, is essential impossible. Thus, the idealized simulations presented here actually provide more reliable information about the behaviour of stripped systems close to the tidal boundary, and may even help to evaluate different group-finding schemes in realistic simulations.

The isotropic NFW models considered in this work do not necessarily provide the most accurate description of cosmological halos, either. The highest-resolution simulations indicate that an Einasto profile is a better fit to the density profile (e.g. Navarro et al., 2010; Klypin et al., 2016), but real halos are also anisotropic, triaxial, and have more complicated correlations between shape and anisotropy (e.g. Wojtak et al., 2013). There are also several proposed models for the ‘true’ DF of dark matter halos, based on maximizing entropy under various constraints, or other arguments (e.g. Hjorth & Williams 2010; Pontzen & Governato 2013; Beraldo e Silva et al. 2014; see Halle et al. 2017 for a discussion). Our focus here is not on which of these theoretical models is correct, but on the generic effects of tidal mass loss on any DF. Our analysis is general and may be applied to these theoretical models, or to any other model with an explicit DF.

Many real systems of interest also have a separate baryonic component that is important, if only as a tracer of dynamics and mass loss. Various authors have considered the density profiles and/or phase-space distributions of luminous stars within a surrounding, tidally limited dark matter potential, and their evolution through tidal interactions (e.g. Mashchenko & Sills, 2004, 2005a, 2005b; Peñarrubia et al., 2008; Sales et al., 2010; Kazantzidis et al., 2011; Amorisco & Evans, 2011). Our focus here is simply on understanding the phase-space evolution of the dark matter particles themselves. We leave further discussion of the evolution of self-consistent, two-component systems to future work.

Acknowledgements

NED acknowledges support from the government of Ontario, through an OGS award. JET acknowledges financial support from NSERC Canada, through a Discovery Grant. The authors also wish to thank Frank van den Bosch and the anonymous referee for useful comments.

\color

black

Appendix A ICICLE – A Code for Generating Isolated Initial Conditions

The authors have created a \colorblackpublicly available Python package, Icicle (Initial Conditions for Isolated CoLlisionless systEms), which can create stable initial conditions for spherical, isotropic, collisionless systems with various density profiles. The code currently supports NFW (either exponentially truncated, abruptly truncated, or truncated \colorblack using the iterative method described in Section 2), King, Hernquist and Einasto profiles.

A.1 Files

The main part of the program is in the code “ICICLE.py". The content of this code is described thoroughly in Section A.3; \colorblack briefly, given a model it returns positions and velocities of $n$ particles within that distribution.

For each model there is an additional file named “ICs_model.py". These files contain information about the specific model (the density profile, cumulative mass distribution, gravitational potential and distribution function). \colorblackThe user can specify parameters by use of a parameter file.

A.2 Output

\color

black The output is written to the filename specified. The first line of the text file has the following information: the number of particles, the mass of each particle and the value of the gravitational constant. Subsequent lines display the particle number (indexed from 0), the x, y and z positions, and the x, y, and z velocity components.

A.3 Positions and Velocities

In this section we outline the steps needed to select positions and velocities given an isotropic density profile, in a manner similar to Kazantzidis et al. (2004). First, the radial distance \colorblack for each particle is selected using the cumulative mass distribution, and then \colorblack a direction for the position vector is selected assuming spherical symmetry. Next, an energy is selected from the energy distribution. Once the position and energy of the particle have been determined, it is straightforward to calculate the velocity of the particle. Finally, the velocity direction is chosen isotropically.

A.3.1 Positions

\color

black The radius for each particle is chosen in such a way as to reproduce the mass distribution. Consider a density profile where the mass within a given radius $r$ is $M(<r)$, and the largest radius is $r_{max}$. The mass fraction interior to $r$ is then:

[TABLE]

\color

black If we choose a random variate $x$ uniformly on the interval $[0,1]$, we can then set $F_{M}(<r)=x$ and solve for $r$ to find a corresponding radius. In the limit of a large number of particles, the resulting set of radii will reproduce the desired density profile. Given a radius, we can then choose a direction isotropically, by choosing a random point on the surface of a unit sphere.

A.3.2 Calculating The Distribution Function

The distribution function, $f(\mathcal{E})$, is given by Eddington’s formula:

[TABLE]

Here $\mathcal{E}$, $\Psi$ and $\rho$ are the relative energy, the relative potential and the density, respectively. \colorblack The second, boundary term is equal to $f(0)$ in general (since the first term evaluates to zero when $\mathcal{E}=0$), so it vanishes for any system with $f(0)=0$.

Consider the term $d^{2}\rho/d\Psi^{2}$. This can be expressed as:

[TABLE]

\color

black This form is more convenient in cases where $\rho$ is an analytic function of radius, as it avoids having to take potentially noisy numerical derivatives.

A.3.3 Choosing From the Distribution Function

The probability that a particle located at radius $r$ has a relative energy $\mathcal{E}$ is proportional to $f(\mathcal{E})\sqrt{\Psi-\mathcal{E}}$, \colorblack with a maximum energy of $\mathcal{E}_{max}=\Psi(r)$.

Thus, the cumulative distribution function (CDF) $F(<\mathcal{E})$, for a particle at position $r$ with relative potential energy $\Psi(r)$, is:

[TABLE]

\color

black Choosing a random variate $y$ uniformly on the interval $[0,1]$, we can set $F(<\mathcal{E}))=y$ and solve for $\mathcal{E}$ to find the corresponding relative energy. In practice, we do this by linear interpolation of a numerically determined CDF.

A.3.4 Velocities

\color

black Once a position and a relative energy have been chosen for the particle, the velocity magnitude can be calculated from:

[TABLE]

\color

black Finally, a direction for the velocity vector can be chosen isotropically in the same way as described for the positions.

A.4 Profile options and documentation

\color

black The code, a list of supported profile types, and full documentation are available as a package, Icicle, at the URL: https://github.com/ndrakos/ICICLE.