Milking the spherical cow – on aspherical dynamics in spherical coordinates

Pontzen, Andrew; Read, Justin I.; Teyssier, Romain; Governato, Fabio; Gualandris, Alessia; Roth, Nina; Devriendt, Julien

doi:10.1093/mnras/stv1032

Abstract

Galaxies and the dark matter haloes that host them are not spherically symmetric, yet spherical symmetry is a helpful simplifying approximation for idealized calculations and analysis of observational data. The assumption leads to an exact conservation of angular momentum for every particle, making the dynamics unrealistic. But how much does that inaccuracy matter in practice for analyses of stellar distribution functions, collisionless relaxation, or dark matter core-creation? We provide a general answer to this question for a wide class of aspherical systems; specifically, we consider distribution functions that are ‘maximally stable’, i.e. that do not evolve at first order when external potentials (which arise from baryons, large-scale tidal fields or infalling substructure) are applied. We show that a spherically symmetric analysis of such systems gives rise to the false conclusion that the density of particles in phase space is ergodic (a function of energy alone). Using this idea we are able to demonstrate that: (a) observational analyses that falsely assume spherical symmetry are made more accurate by imposing a strong prior preference for near-isotropic velocity dispersions in the centre of spheroids; (b) numerical simulations that use an idealized spherically symmetric setup can yield misleading results and should be avoided where possible; and (c) triaxial dark matter haloes (formed in collisionless cosmological simulations) nearly attain our maximally stable limit, but their evolution freezes out before reaching it.

galaxies: evolution, galaxies: formation, galaxies: haloes

1 INTRODUCTION

Spherical symmetry is a foundational assumption of many dynamical analyses. The primary motivation is simplicity, since few astronomical objects are actually spherical. For example, observations and simulations both suggest that gravitational potential wells generated by dark matter haloes are typically triaxial (e.g. Dubinski & Carlberg 1991; Cole & Lacey 1996; Jing & Suto 2002; Kasun & Evrard 2005; Hayashi, Navarro & Springel 2007; Loebman et al. 2012; Schneider, Frenk & Cole 2012). Characterizing dark matter haloes by spherically averaged densities and velocities (e.g. Dubinski & Carlberg 1991; Navarro, Frenk & White 1996; Taylor & Navarro 2001; Stadel et al. 2009) at best tells only part of the story. At worst, it could be severely misleading.

The question of whether baryonic processes can convert dark matter cusps into cores (Pontzen & Governato 2014) provides one motivation for a detailed study of the relationship between spherical and near-spherical dynamics. To explain why, we need to look forward to some of our results. Later in this paper, we cut a dark matter halo out of a cosmological simulation, then match it to an exactly spherical halo with an identical density and velocity anisotropy profile. This gives us two easy-to-compare equilibrium structures – the first triaxial, the second spherical – to perform a dynamical comparison. We expose each to the same time-varying gravitational potential, mimicking the effects of stellar feedback (there are no actual baryons in these runs). After 1 Gyr, the triaxial halo's averaged density profile flattens into a convincing dark matter core, but the spherical halo maintains its cusp (see Fig. 1).

Figure 1.

One motivation for studying the relationship between spherical and aspherical dynamics is that the conversion of a dark matter cusp into a core by baryonic processes is qualitatively different in the two cases. Here we show the inner log density slope from a numerical experiment on two haloes. One is spherical (dashed line) and the other triaxial (solid line) but their spherically averaged properties are initially identical. An external potential has been added at the centre during the times indicated by the grey bands, with the fluctuations mimicking stellar feedback. The triaxial halo develops a clear core, whereas the spherical halo almost maintains its central density cusp. A complete description and analysis is given in Section 3.5.

Open in new tab Download slide

This example, which is fully explored in Section 3.5, illustrates how it is dangerous to use spherically symmetric simulations to infer anything about dynamics – even spherically averaged dynamics – in the real universe. A spherical system does not evolve in the same way as the spherical averages of a triaxial system.

Ignoring asphericity can also lead to observational biases (e.g. Hayashi & Navarro 2006; Corless & King 2007). From a dynamical standpoint, the nature of orbits in triaxial potentials is fundamentally different from those in spherical potentials: although the total angular momentum of any self-gravitating system must always be conserved, it is only in the spherical case that this conservation holds for individual particles. Conversely, a large fraction of dark matter particles near the centre of cosmological haloes will be on box orbits which do not conserve their individual angular momenta (de Zeeuw & Merritt 1983; Merritt & Valluri 1996; Holley-Bockelmann et al. 2001; Adams et al. 2007). One practical consequence is that asphericity may be responsible for filling the loss cones of supermassive black holes at the centre of the corresponding galaxies (Merritt & Poon 2004).

Finally, it is known that asphericity plays a fundamental role in setting the equilibrium density profile during gravitational cold collapse (e.g. Huss, Jain & Steinmetz 1999). The underlying process is known as the radial orbit instability (ROI; Henon 1973; Saha 1991; Huss et al. 1999; MacMillan, Widrow & Henriksen 2006; Bellovary et al. 2008; Barnes, Lanzel & Williams 2009; Maréchal & Perez 2009); a related effect was discussed by Adams et al. (2007). The name arises because particles on radial orbits are perturbed on to more circular trajectories. At the same time, the density distribution becomes triaxial. Even in the case of a uniform spherical collapse, this symmetry-breaking process is still triggered, presumably by numerical noise; the tangential component of forces must be unphysically suppressed for the system to remain spherical (Huss et al. 1999; MacMillan et al. 2006).

Despite all this, assuming spherical symmetry is very tempting because it makes life so much easier. Defining spherically averaged quantities is a well-defined and sensible procedure even if we actually have the full distribution function in hand (as in simulations): departures from spherical symmetry are sufficiently small that different averaging procedures lead to consistent results (Saha & Read 2009). Additionally, when an aspherical halo is in equilibrium, we have shown numerically that a ‘sphericalized’ version of it is also in equilibrium (see appendix B of Pontzen & Governato 2013). This is helpful because it allows one to make a meaningful analysis in spherical coordinates, even when the system is aspherical. But it breaks down when out-of-equilibrium processes are included, as in the stellar-feedback-driven core-creation example above.

The present paper formalizes the idea of spherical analysis performed on aspherical systems and follows through the consequences. We will study equilibrium distribution functions in nearly, but not exactly, spherically symmetric potentials, and focus on maximally stable systems (which we define as being stable against all possible external linear perturbations). We will find that, in spherical coordinates, such systems appear to be ‘ergodic’ (meaning that their distribution functions depend on energy alone, following the terminology of Binney & Tremaine 2008) because the individual particles move randomly in angular momentum while maintaining a near-constant energy. It is important to emphasize that this describes the appearance of the system when analysed in a spherical coordinate system and the true system need not be chaotic for the result to hold, provided any isolating integrals are not closely related to angular momentum.

The formal statement of this idea is derived in Section 2. A brief overview of the required background is given in Appendix A, and a second-order derivation of the evolution is given in Appendix B. Section 3 outlines the practical consequences, starting by recasting and extending the phenomenology of the ROI. We describe an immediate implication for observational studies of aspherical systems, namely a new way to break the anisotropy degeneracy. Finally we return to the motivating problem above and explain why triaxial systems can undergo cusp–core transitions more easily than spherical systems. Section 4 concludes and points to open questions and future work.

2 ASPHERICAL DYNAMICS IN SPHERICAL COORDINATES

In this section we consider an aspherical system which is maximally stable against external linear perturbations. We assume that an observer of this system analyses it assuming spherical symmetry. We will show that this observer (falsely) concludes that the system is ergodic, i.e. that the density of particles in phase space is a function of energy alone. The derivation requires the use of action-angle coordinates; a crash course is provided in Appendix A.

2.1 Single particles

Given any near-spherical system, the Hamiltonian in the spherical action-angle variables is

\begin{equation} H(\boldsymbol {J}, \boldsymbol {\Theta }) = H_0(\boldsymbol {J}) + \delta H(\boldsymbol {J}, \boldsymbol {\Theta })\textrm {,} \end{equation}

(1)

where H₀ is the sum of kinetic and potential energies in the spherical background, |$\boldsymbol {J}=(J_{\rm r},j,j_z)$| is the vector of spherical actions (see Appendix A), |$\boldsymbol {\Theta }$| is the vector of spherical angles and δH contains the perturbation (which includes the aspherical correction to the potential). The orbit of a particle in exact spherical symmetry, δH = 0, is described by Hamilton's equations:

\begin{equation} \dot{\boldsymbol {J}}_0 = -\frac{\mathrm{\partial} H_0}{\mathrm{\partial} \boldsymbol {\Theta }} = 0 \textrm {;} \hspace{28.45274pt} \dot{\boldsymbol {\Theta }}_0 = \left. \frac{\mathrm{\partial} H_0}{\mathrm{\partial} \boldsymbol {J}} \right|_{\boldsymbol {J}=\boldsymbol {J}_0} \equiv \Omega _0(\boldsymbol {J}_0)\textrm {,} \end{equation}

(2)

which defines the constant background orbital frequencies |$\Omega _0(\boldsymbol {J}_0)$|⁠. The expressions |$\boldsymbol {J}_0$|⁠, |$\boldsymbol {\Theta }_0$| and |$\Omega _0(\boldsymbol {J}_0)$| will be used throughout to refer to the background (δH = 0) solution. This algebraically simple form of the equations of motion is the reason for using action-angle variables, since it immediately integrates to

\begin{equation} \boldsymbol {J}_0(t) = \boldsymbol {J}_0 = \textrm {constant; }\hspace{14.22636pt} \boldsymbol {\Theta }_0(t) = \boldsymbol {\Theta }_0(0) + \boldsymbol {\Omega }_0(\boldsymbol {J}_0) t\textrm {,} \end{equation}

(3)

where |$\boldsymbol {J}_0$| and |$\boldsymbol {\Theta }_0(0)$| specify the initial action and angle coordinates of the orbit.

We now consider the effect of the aspherical correction to the potential encoded in δH, using standard Hamiltonian perturbation theory (e.g. Lichtenberg & Lieberman 1992; Binney & Tremaine 2008). First, taking advantage of the angle coordinates |$\boldsymbol {\Theta }$| being periodic in 2π, δH is expressed as

\begin{equation} \delta H(\boldsymbol {J},\boldsymbol {\Theta }) = \sum _{\boldsymbol {n}} \delta H_{\boldsymbol {n}}(\boldsymbol {J}) {\rm e}^{{\rm i}\boldsymbol {n}\cdot \boldsymbol {\Theta }}\textrm {.} \end{equation}

(4)

This equation states that, at any fixed |$\boldsymbol {J}$|⁠, one can expand the periodic Θ dependence in a Fourier series without loss of generality.

We are interested in the evolution of |$\boldsymbol {J}$| at first order in the perturbation. Hamilton's relevant equation now reads:

\begin{equation} \dot{\boldsymbol {J}} = - \frac{\mathrm{\partial} H}{\mathrm{\partial} \boldsymbol {\Theta }} = - \sum _{\boldsymbol {n}} {\rm i} \boldsymbol {n} \delta H_{\boldsymbol {n}}(\boldsymbol {J}) {\rm e}^{{\rm i} \boldsymbol {n} \cdot \boldsymbol {\Theta }}\textrm {.} \end{equation}

(5)

Because δH is small, the result to first-order accuracy is given by substituting the zero-order solution (3) into equation (5) and integrating to give

\begin{equation} \boldsymbol {J}(t) = \boldsymbol {J}_0 - \sum _{\boldsymbol {n}} \frac{\boldsymbol {n}}{\boldsymbol {n} \cdot \boldsymbol {\Omega }_0} \delta H_{\boldsymbol {n}}(\boldsymbol {J}_0) {\rm e}^{{\rm i}\boldsymbol{n} \cdot \boldsymbol {\Theta }_0(t)} + \cdots \textrm {.} \end{equation}

(6)

Consequently, as |$\boldsymbol {n} \cdot \boldsymbol {\Omega }_0 \rightarrow 0$|⁠, the linear-order correction to the orbit of a particle can become large even if the aspherical correction to the potential (δH) is small, an effect known as resonance (Binney & Tremaine 2008). Consider now the evolution of the background energy along the perturbed trajectory, given by

\begin{equation} H_0(\boldsymbol {J}(t)) \simeq H_0(\boldsymbol {J}_0) + \frac{\mathrm{\partial} H_0}{\mathrm{\partial} \boldsymbol {J}} \cdot \left(\boldsymbol {J}(t) - \boldsymbol {J}_0\right) + \cdots \textrm {,} \end{equation}

(7)

where we have Taylor-expanded to first order around |$\boldsymbol {J}_0$|⁠. Substituting equation 6 for J(t), there is a cancellation between numerator and denominator:

\begin{equation} H_0(\boldsymbol {J}(t)) \simeq H_0(\boldsymbol {J}_0) - \sum _{\boldsymbol {n}}\delta H_{\boldsymbol {n}}(\boldsymbol {J}_0) {\rm e}^{{\rm i} \boldsymbol {n} \cdot \boldsymbol {\Theta }_0(t)} +\cdots \end{equation}

(8)

and so the fractional variation in H₀ remains small, even if |$\boldsymbol {J}$| changes significantly over time. In other words, according to linear perturbation theory, particles migrate large distances in |$\boldsymbol {J}$| along surfaces of approximately constant background energy H₀. One can verify this constrained-migration prediction in numerical simulations of dark matter haloes – an explicit demonstration is given in Appendix A3. This gives us some intuition for the result to come: a population of particles will seem to ‘randomize’ their actions (including angular momentum), but not their energy distribution.

The extent of the migration will depend on the nature of the potential in which a particle orbits. To quantify this requires going beyond linear perturbation theory and is the subject of ‘KAM theory’ after Kolmogorov (1954), Arnold (1963) and Moser (1962); see e.g. Binney & Tremaine (2008), Goldstein, Poole & Safko (2002), Lichtenberg & Lieberman (1992) for introductions. Colloquially the result is that for any given small perturbation the migration of typical orbits is also small. Arnold diffusion offers the most famous route to more significant diffusion through action space (see Lichtenberg & Lieberman 1992); but in our case, there is a more immediate reason why the KAM result does not in fact hold. Specifically, KAM theory relies on the frequencies |$\boldsymbol {\Omega }(\boldsymbol {J})$| being non-degenerate – i.e. that any change in the action leads also to a change in the frequencies, thus shutting off resonant migration. In smooth potentials, Ω is almost a function of energy alone (see Appendix A3) and so the migration can be substantial.

Overall, we informally expect particles to redistribute themselves randomly within the action shell of fixed background energy until they are evenly spread, implying a distribution function that appears ergodic in a spherical analysis. This does not imply the orbits are chaotic in the traditional sense; it is only because we are analysing an aspherical system in spherical coordinates that the phenomenon arises. With this in mind, we now turn to a more formal demonstration of the result.

2.2 The distribution function

So far, we have discussed how a single particle orbiting in a mildly aspherical potential does not conserve its spherical actions (e.g. angular momentum). We informally suggested that a population will appear to ‘randomize’ the spherical actions at fixed energy. We now show more formally that a distribution function of particles subject to aspherical perturbations will be most stable when it is spread evenly on surfaces of constant H₀.

We start by decomposing the true distribution function of particles in phase space, f, in terms of a spherical background f₀ and a perturbation δf. To make sure the split between spherical background and aspherical perturbation is uniquely defined, we take f₀ as the distribution function obtained when we perform a naive analysis averaging out the aspherical contribution:

\begin{equation} f_0(\boldsymbol {J}) \equiv \frac{1}{(2\pi)^3} \int \mathrm{d}^3\Theta f(\boldsymbol {J},\boldsymbol {\Theta })\textrm {.} \end{equation}

(9)

By Jeans’ theorem, f₀ is an equilibrium distribution function in the spherical background because it is constructed from spherical invariants |$\boldsymbol {J}$| alone. Analogous to equations (1) and (4) one can write the full distribution function f as

\begin{eqnarray} f(\boldsymbol {J},\boldsymbol {\Theta }) & =& f_0(\boldsymbol {J}) + \delta f(\boldsymbol {J}, \boldsymbol {\Theta }) \nonumber \\ & = &f_0(\boldsymbol {J}) + \sum _{\boldsymbol {n}} \delta f_{\boldsymbol {n}}(\boldsymbol {J}) {\rm e}^{{\rm i} \boldsymbol {n} \cdot \boldsymbol {\Theta }}\textrm {.} \end{eqnarray}

(10)

The whole f is to be in equilibrium in the true system, which implies that ∂f/∂t = 0. We can turn this into an explicit condition on f using the collisionless Boltzmann equation,

\begin{equation} 0=\frac{\mathrm{\partial} f}{\mathrm{\partial} t} = [H,f] \equiv \frac{\mathrm{\partial} H}{\mathrm{\partial} \boldsymbol {\Theta }} \cdot \frac{\mathrm{\partial} f}{\mathrm{\partial} \boldsymbol {J}} - \frac{\mathrm{\partial} H}{\mathrm{\partial} \boldsymbol {J}} \cdot \frac{\mathrm{\partial} f}{\mathrm{\partial} \boldsymbol {\Theta }}\textrm {.} \end{equation}

(11)

Expanding to linear order in H and f gives the condition

\begin{equation} \sum _{\boldsymbol {n}} \left( \boldsymbol {\Omega }_0(\boldsymbol {J}) \cdot \boldsymbol {n} \delta f_{\boldsymbol {n}} (\boldsymbol {J}) - \frac{\mathrm{\partial} f_0}{\mathrm{\partial} \boldsymbol {J}} \cdot \boldsymbol {n} \,\delta H_{\boldsymbol {n}}(\boldsymbol {J}) \right) {\rm e}^{{\rm i} \boldsymbol {n} \cdot \boldsymbol {\Theta }} = 0\textrm {.} \end{equation}

(12)

The different |$\boldsymbol {\Theta }$| dependence of each term in the sum means that the term in brackets must be zero for each |$\boldsymbol {n}$|⁠. We can read most of these separate terms as determining which |$\delta f_{\boldsymbol {n}}$| are required for a stable distribution function under a given aspherical potential correction |$\delta H_{\boldsymbol {n}}$|⁠. However, for the resonant contributions (⁠|$\boldsymbol {n}=\boldsymbol {n}_{\perp }$| where |$\boldsymbol {\Omega }_0 \cdot \boldsymbol {n}_{\perp }=0$|⁠) the |$\delta f_{\boldsymbol {n}}$| term vanishes to leave the condition

\begin{equation} \left(\frac{\mathrm{\partial} f_0}{\mathrm{\partial} \boldsymbol {J}}\cdot \boldsymbol {n}_{\perp }\right) \,\delta H_{\boldsymbol {n}_{\perp }}(\boldsymbol {J})=0\textrm {.} \end{equation}

(13)

To be able to find a stable |$\delta f_{\boldsymbol {n}}$| for any given |$\delta H_{\boldsymbol {n}}$|⁠, one therefore requires that f₀ is a function only of H₀, since then |$\mathrm{\partial} f_0 / \mathrm{\partial} \boldsymbol {J} = \boldsymbol {\Omega }_0\, \mathrm{d}f_0 / \mathrm{d}H_0$| and consequently the dot product in equation (13) vanishes.

This is the core result claimed at the start of the section: f₀ = f₀(H₀), i.e. the distribution function implied by a spherical analysis appears to be ergodic. It is not a necessary condition for achieving equilibrium, since for any given aspherical system certain |$\delta H_{\boldsymbol {n}_{\perp }}$| will be zero. Rather, the result should be read as applying to the maximally stable distribution function – a distribution function that can self-stabilize under any linear perturbation to its potential.

We again emphasize that the distribution function f₀ is fictional. There is no sense in which the true distribution function, f, is actually ergodic. The statement is about how the system appears to be when it is analysed using spherically averaged quantities, equation (9). Yet, it establishes a way in which we can understand these spherically averaged quantities in a systematic way – the system is most stable if it appears ergodic, regardless of what the underlying dynamics is really up to. In the remainder of this paper we will refer to such systems as ‘spherically ergodic’ (SE).

2.3 Testable predictions

We have established that systems which appear to be ergodic in a spherical analysis are maximally stable. Now we need to devise a connection to observable or numerically measurable quantities.

A distribution function f(H₀) that is truly a function of energy alone has an isotropic velocity distribution (Binney & Tremaine 2008). To test for isotropy, one calculates β(r) according to the usual spherically averaged definition

\begin{equation} \beta (r) = 1 - \frac{\langle v_{\rm t}^2\rangle (r)}{2 \langle v_{\rm r}^2 \rangle (r)} \end{equation}

(14)

where v_r is a particle's radial velocity, v_t its tangential velocity and the angle-bracket averages are taken in radial shells. For a population on radial orbits, β(r) = 1; conversely for purely circular motion, β(r) = −∞. Between these two extremes, an ergodic population has β(r) = 0 (Binney & Tremaine 2008).

Intuitively, a SE system (in the sense defined in the previous section) should therefore be approximately isotropic. However one has to handle that expectation with a little care because the true population f is not ergodic and the measured velocity dispersions, even in spherical polar coordinates, may inherit information from f that is not present in f₀.

Instead, we will now construct a more rigorously justifiable, slightly different statement that still connects SE populations to velocity isotropy. Measuring the mean of any function of the spherical actions |$q(\boldsymbol {J})$|⁠, we obtain

\begin{equation} \int \mathrm{d}^3 J\, \mathrm{d}^3 \Theta \, f(\boldsymbol {J},\boldsymbol {\Theta }) q(\boldsymbol {J}) = \int \mathrm{d}^3 J \, f_0\left(H_0(\boldsymbol {J})\right) q(\boldsymbol {J})\textrm {,} \end{equation}

(15)

an exact result. Therefore any statement about averages over spherical actions carries information only about f₀ – the spherical part of the distribution function. This allows us to derive unambiguous implications of a SE population.

The most familiar action is the specific scalar angular momentum j. Because it is a scalar for each particle, averages over this quantity do not express anything about a net spin of the halo but rather about the mix of circular and radial orbits, just like the traditional velocity anisotropy. Radially biased populations have 〈j〉 ≃ 0 whereas populations on circular orbits have 〈j〉 = j_c, where j_c is the maximum angular momentum available at a given energy. So, velocity anisotropy can be conveniently represented in terms of the mean scalar angular momentum.

We can go further and calculate a function, 〈j〉(E), where the average is taken only over particles at a particular specific energy. This quantity can be represented in terms of the ratio of two integrals of the form (15):

\begin{eqnarray} \langle j \rangle (E) &=& \frac{\iiint \mathrm\;{d}J_{\rm r}\, \mathrm{d}j\, \mathrm{d}j_z\, f_0(H_0) j \delta (H_0-E)}{\iiint \mathrm\;{d}J_{\rm r}\, \mathrm{d}j\, \mathrm{d}j_z\, f_0(H_0) \delta (H_0-E)}\textrm {.} \end{eqnarray}

(16)

The triple integral ranges over the physical phase space coordinates: 0 ≤ J_r < ∞, 0 ≤ j < ∞, −j ≤ j_z ≤ j. One can immediately perform the j_z integrals; then the J_r integral can be completed by changing variables to H₀ (recalling ∂H₀/∂J_r ≡ Ω_r) and consuming the δ function. After this manipulation j can only range between 0 and j_c(E) where j_c(E) is the specific angular momentum corresponding to a circular orbit with specific energy E; there are no physical orbits with more angular momentum at the specified E. The final, exact result is:

\begin{eqnarray} \langle j \rangle (E) = \int _0^{j_{\rm c}(E)} \mathrm{d}j\, \Omega _{\rm r}(E,j)^{-1}\, j^2 \Big/ \int _0^{j_{\rm c}(E)} \mathrm{d}j\, \Omega _{\rm r}(E,j)^{-1}\, j. \end{eqnarray}

(17)

We now have a firm prediction for SE populations. Namely, if we bin particles in E and measure 〈j〉(E) in each bin, the results should be predicted by equation (17), which is a function only of the potential (through Ω_r). Equation (17) does not exhaust the possible tests for spherical ergodicity, but it is sufficient for our present exploratory purposes.

For smooth potentials, Ω_r(E, j) varies very little between j = 0 and j = j_c, and one can approximate it very well as a function of E alone (Appendix A3). In this case, the integrals follow analytically and one has the result

\begin{equation} \langle j \rangle (E) \simeq \frac{2}{3} j_{\rm c}(E)\textrm {.} \end{equation}

(18)

This is a helpful simplification to set expectations, but throughout this paper when showing the SE limit, we will use the exact expression given by equation (17).

3 EXAMPLE APPLICATIONS

So far, we have motivated and derived a formal result that aspherical systems are most stable when they appear ergodic in spherical coordinates. We derived one practical consequence for the angular momentum distribution, equation (17), which in an approximate sense states that the velocity distribution will appear isotropic. We are now in a position to test whether numerical simulations actually tend towards this maximally stable limit in a variety of situations, beginning with cosmological collisionless dark matter haloes.

3.1 Cosmological dark matter haloes

Let us re-examine the three high-resolution, dark-matter-only zoom cosmological simulations used in the analysis of Pontzen & Governato (2013). The three each have several million particles in their z = 0 haloes which correspond in turn to a dwarf irregular, L_⋆ galaxy and cluster. The force softening ϵ, virial radius r₂₀₀ (at which the mean density enclosed is 200 times the critical density) and virial masses M₂₀₀ are 65, 170, 690 pc; 98, 301, 1430 kpc and 2.8 × 10¹⁰, 8.0 × 10¹¹, 8.7 × 10¹³ M_⊙ respectively. For further details see Pontzen & Governato (2013).

The upper panel of Fig. 2 shows the anisotropy for our cosmological haloes. To compare the three directly, we scale the radius by r_max (respectively 27, 57 and 340 kpc), the radius at which the circular velocity (GM(< r)/r)^1/2 reaches its maximum, v_max (=56, 150 and 610 km s⁻¹). We restrict attention to the region well within the virial radius; here, the anisotropy β(r) typically lies between the purely radial and isotropic cases (e.g. Bellovary et al. 2008; Navarro et al. 2010).

Figure 2.

The velocity anisotropy of the inner parts of three sample high-resolution cosmological dark matter haloes (simulated without baryons), plotted as a function of radius (upper panel) and in energy shells (lower panel). The upper panel shows the classic velocity anisotropy β(r) defined in the text, for which a purely radial population has β = 1 and a population on circular orbits β = −∞. The lower panel shows our alternative in energy space which can be more precisely related to the theoretical arguments presented in Section 2; here 〈j〉(E)/j_circ = 0 for radial orbits and 1 for circular orbits, and approximately 2/3 for an isotropic population. Both panels show that the haloes have near-isotropic orbits with a slight radial bias. The range of the two plots is roughly comparable, but we caution that the mapping from r to E is not unique (see Fig. 3).

Open in new tab Download slide

We now want to link this relatively familiar velocity anisotropy to the alternative 〈j〉(E) statistic that was directly predicted by the SE property in Section 2. For each particle we calculate the specific energy |$E = \dot{\boldsymbol {x}}^2/2 + \Phi \left(|\boldsymbol {x}|\right)$|⁠, where |$\boldsymbol {x}$| is the vector displacement from the halo centre. To make this quantity agree exactly with H₀ in the terminology of Section 2, we ignore asphericity when calculating the potential energy, defining it as

\begin{equation} \Phi (r) \equiv \int _0^r \frac{G M(<r^{\prime })}{r^{\prime 2}} \mathrm{d}r^{\prime } \end{equation}

(19)

where M( < r′) is the mass enclosed inside a sphere of radius r′. (The numerical integration is performed by binning particles in shells of fixed width ϵ, chosen to coincide with the force softening in the simulation; within these bins the density is taken to be constant.) The physics is invariant if a constant is added to the potential; we have chosen to fix its scale by setting Φ(0) = 0.

For each particle we also calculate the specific angular momentum |$j=\left| \boldsymbol {x} \times {\dot{\boldsymbol {x}}} \right|$|⁠, and the specific angular momentum of a circular orbit at the same energy, j_circ(E) which is given by simultaneously solving

\begin{eqnarray} E &=& \Phi (r)+\frac{j_{\mathrm{circ}}^2}{2 r^2} \hspace{14.22636pt}\textrm { and }\hspace{14.22636pt} \Phi ^{\prime }(r) = \frac{j_{\mathrm{circ}}^2}{r^3}\textrm {,} \end{eqnarray}

(20)

to eliminate r in favour of j_circ.

We plot 〈j〉(E)/j_circ(E) in bins containing 1000 particles each in the lower panel of Fig. 2. To facilitate comparison with the top panel, a population on purely radial orbits would have β = 1 and 〈j〉/j_circ = 0, whereas a purely circular distribution function corresponds to β = −∞ or 〈j〉/j_circ = 1. Isotropic, purely spherical populations have β(r) = 0 and 〈j〉/j_circ ≃ 2/3 as discussed at the end of Section 2.3. When compared against each other in this way, the two panels agree well: the populations are on near-isotropic orbits with a slight radial bias.

Quantitatively, how well do these results agree? As a guideline, we can compare results for models with constant anisotropy. From Binney & Tremaine (2008), a distribution function f(j, E) = j^−2βf₁(E) generates a constant anisotropy β(r) = β. We can calculate the connection to the new statistic by generalizing the reasoning of Section 2.3, with the result that

\begin{equation} \langle j \rangle (E) = \frac{\int _0^{j_{\mathrm{circ}}} \mathrm{d}j \, \Omega _{\rm r}^{-1} \, j^{2-2\beta }}{\int _0^{j_{\mathrm{circ}}} \mathrm{d}j \, \Omega _{\rm r}^{-1} j^{1-2\beta }} \simeq \frac{2-2\beta }{3-2\beta } j_{\mathrm{circ}}\textrm {,} \end{equation}

(21)

where the first result is exact and the second follows from assuming Ω_r is independent of j (which is an excellent approximation; see Appendix A3). Consistent with equation (18), β = 0 gives 〈j〉(E)/j_circ ≃ 2/3. But as the system becomes more radially biased we can now calculate that, for example, β = 0.2 corresponds to 〈j〉(E)/j_circ ≃ 0.62. Despite the various approximations involved, these values therefore correctly relate the values of β in the top panel of Fig. 2 with the 〈j〉 results in the lower panel.

That said, a detailed comparison as a function of radius is hard because particles at a given radius r have a wide spread of energies E. Fig. 3 illustrates this relationship for the ‘MW’ halo. The density shows the probability of a particle at radius r also having specific energy E, p(E|r). The minimum E at each r is set by Φ(r), which gives the energy of a particle at apocentre (solid line). A more typical E is given by the energy of a circular orbit at r (dashed line), and this gives some intuition for mapping results from the top panel of Fig. 2 on to the bottom panel. However, any E exceeding Φ(r) is theoretically possible. So β at any given radius actually represents an average over particles of many different energies.

The relationship between r/rmax and E/vmax2 for circular orbits (dashed line), radial orbits at apocentre (solid line) and particles in the ‘MW’ run (density shows the number of particles with energy E at each radius r). The mapping between energy and radius is fuzzy, so that anisotropy at high E can easily contaminate β(r) at small r.

Figure 3.

The relationship between r/r_max and E/v_max² for circular orbits (dashed line), radial orbits at apocentre (solid line) and particles in the ‘MW’ run (density shows the number of particles with energy E at each radius r). The mapping between energy and radius is fuzzy, so that anisotropy at high E can easily contaminate β(r) at small r.

Open in new tab Download slide

Before moving on, we should remind ourselves that 〈j〉(E) was just one of many possible tests for SE populations that follow from equation (15); we chose it because it ties most closely to the better-known β(r) statistic. Therefore, a population close to the dashed lines in Fig. 4 could in principle look far from being ergodic by other measures. For example, instead of being evenly spread in angular momentum, the individual particles could all cluster tightly around the mean value j ≃ 2j_c/3. We have checked for all results in this paper that the particles really do spread evenly through j/j_c between zero and one [noting that there is a degeneracy factor arising from the j_z direction, so that p(j) ∝ j]. However for simplicity, all our figures focus on the 〈j〉(E) relation.

Figure 4.

The β(r) velocity anisotropy (upper panel) and its equivalent in energy space 〈j〉(E) (lower panel, as Fig. 2) for the evolution over time of a halo that is initially in equilibrium, but unstable against the ROI. The initial conditions are represented by the dark blue line, t = 0; there is a delay of 12t_dyn before any significant evolution can be seen (light blue line). Subsequent lines are shown every 4t_dyn. Once the instability kicks in and generates aspherical perturbations, there is a rapid evolution towards the predicted SE limit (lower panel).

Open in new tab Download slide

3.2 The classical radial orbit instability

Having established that, loosely speaking, 〈j〉(E) represents the anisotropy in energy shells in the same way that β(r) does in radial shells, we can return to our prediction (Section 2.3) for the former quantity, which is shown by the dashed line in the lower panel of Fig. 2. The prediction is almost, but not quite, satisfied in cosmological dark matter haloes. Since the condition is only reached in a maximally stable object, approximate agreement is an acceptable situation.

Because cosmological haloes initially form from near-cold collapse, the ROI (van Albada 1982; Barnes, Hut & Goodman 1986; Saha 1991; Weinberg 1991) is invoked to explain how the radially infalling material gets scattered on to a wider variety of orbits (MacMillan et al. 2006; Bellovary et al. 2008; Barnes et al. 2009), isotropizing the velocity dispersion. Most tellingly, numerical experiments by Huss et al. (1999); MacMillan et al. (2006) show that suppressing the instability (by switching off non-radial forces) results in a qualitatively different density profile as the end-point of collapse. We will now show that the isotropization of velocity dispersion during the ROI can be interpreted in terms of an evolution towards stability in the terms of Section 2. A similar analysis can be applied to populations on initially tangentially biased orbits (Barnes et al. 1986), but we will focus here on the radial case since it is of more cosmological relevance (see Introduction).

Consider what happens to a halo that is intentionally designed to be unstable. First, we will show the classic ROI at work by constructing a spherical halo with particles that are on radially biased orbits. We initialize our particles such that they solve the Boltzmann equation and so are stable in exact spherical symmetry. In practice, however, the strong radial bias means that any slight numerical noise will trigger the ROI. By initializing an unstable equilibrium in this way, we avoid confusion from violent relaxation processes associated with out-of-equilibrium collapse (Lynden-Bell 1967).

Specifically, the initial conditions are set up in a similar fashion to Read et al. (2006), with particle positions drawn from a generalized Hernquist density profile (Hernquist 1990; Dehnen 1993):

\begin{equation} \rho (r) = \frac{M(3-\gamma )}{4\pi a^3} \left(\frac{r}{a}\right)^{-\gamma }\left(1+\frac{r}{a}\right)^{\gamma -4}\textrm {,} \end{equation}

(22)

which has a circular velocity reaching a peak at r_max = (2 − γ)a and implies the gravitational potential

\begin{equation} \Phi (r) = \frac{GM}{2-\gamma }\left[\left(1+a/r\right)^{\gamma -2}-1\right]. \end{equation}

(23)

We choose γ = 1 to roughly mimic an NFW halo (Navarro, Frenk & White 1997) in the innermost parts of interest. The velocities are sampled (using an accept–reject algorithm) from a numerically calculated Osipkov–Merritt distribution function (Binney & Tremaine 2008):

\begin{equation} f(Q) = \frac{\sqrt{2}}{4 \pi^2} \frac{\mathrm{d}}{\mathrm{d}Q} \int _Q^0 \frac{\mathrm{d}\Phi }{\sqrt{\Phi - Q}} \frac{\mathrm{d}}{\mathrm{d}\Phi } \left[\left( 1+r(\Phi )^2/r_{\rm a}^2\right) \rho (r(\Phi ) \right] \!\!\!\!\!\! \end{equation}

(24)

where the parameter Q is defined by

\begin{equation} Q \equiv E + \frac{j^2}{2r_{\rm a}^2}\textrm {.} \end{equation}

(25)

The value of r_a is known as the anisotropy radius because the velocity anisotropy is given by

\begin{equation} \beta (r) = \frac{1}{1 + r_{\rm a}^2/r^2}\textrm {,} \end{equation}

(26)

showing that β(r) ≃ 0 (isotropic) for r ≪ r_a and β(r) ≃ 1 (radial) for r ≫ r_a. We used the minimum value of r_a for which the distribution function is everywhere positive, making the orbits as radially biased as possible; for γ = 1, this is r_a ≃ 0.21 (Meza & Zamorano 1997). We draw 10⁶ particles and evolve the system using ramses (Teyssier 2002), with mesh refinement based on the number of particles per cell, resulting in a naturally adaptive force softening reaching a minimum of ϵ = 90 pc.

Our expectation that numerical noise triggers the ROI is borne out by the numerical experiment. The upper panel of Fig. 4 shows the radial anisotropy β(r) over time. We have defined a single dynamical time, t_dyn, at the peak of the velocity curve so that t_dyn ≡ r_max/v_max. The six solid lines show the population at t = 0, 12, 16, 20, 24 and 28 dynamical times. At first, β(r) appears stable, but suddenly after 16 t_dyn it becomes considerably more isotropic. (This time-scale is sensitive to the particle resolution, although quite weakly; with half the number of particles, each double the mass, the instability sets in around 2 t_dyn earlier.) At the same time, asphericity in the potential develops. To demonstrate this, we determine the inertia tensor of the entire density distribution and calculate the ratio of the principal axes; at t = 0, the ratio is 1.0 by construction. By 16 t_dyn the ratio is 1.8. It stabilizes at around 26 t_dyn with a ratio of 4.3.

This is symptomatic of the classic ROI in action. We can follow the same process from our energy/angular-momentum standpoint in the bottom panel of Fig. 4. The initial conditions (Δt = 0) have a kink in them, with most regions of energy space appearing radially biased but some showing a slight circular preference (around E ≃ 1.5v_max²). We verified that this is an artefact of the Opsikov–Merritt construction and is consistent with the radial β(r) being radially biased everywhere; recall Fig. 3 shows that the relationship between r and E is non-trivial.

As time progresses, the 〈j〉(E) distribution correctly tends towards the SE limit, as predicted. The SE limit is attained to good accuracy for energies E < v_max² by Δt = 20 t_dyn. At larger energies, it is likely difficult to achieve SE because of the long time-scales and weak gravitational fields involved. Like 〈j〉, β(r) at 28 t_dyn shows an isotropic distribution at the centre (r → 0); but β(r) > 0 everywhere for r > 0.1 r_max. We verified that this continuing radial bias throughout the halo is produced by a number of high-energy, loosely bound particles plunging through, i.e. the high-E particles from the lower panel in Fig. 4 pervade all radii in the upper panel.

Let us briefly recap: what we have so far discussed provides a new view of an existing phenomenon. We have recast the ROI and its impact on spherically averaged quantities as an evolution towards an analytically derived maximally stable class of distribution functions. Viewed in energy shells instead of radial shells, the distinction between the regions that reach the SE limit and the slowly evolving, loosely bound regions that remain radially biased is considerably clearer. Now, because our underlying analytic result is derived without requiring the population to be self-gravitating, we can now go further and consider a different version of the ROI – and show that it continues to operate even when a system is in a completely stable equilibrium.

3.3 The continuous radial orbit instability

Our next target for investigation is to broaden the conditions: according to the results of Section 2, a subpopulation of particles should undergo something much like the ROI even when the global potential is completely stable. A suitable name for this phenomenon would seem to be the ‘continuous ROI’, since it continues indefinitely after the global potential has stabilized.

We perform another numerical experiment to demonstrate the effect. First, to avoid confusion from cosmological infall and tidal fields, we create a stable, isolated, triaxial cosmological halo by extracting from our cosmological run a region of 3r₂₀₀ around our ‘Dwarf’. We then evolve this isolated region for 2 Gyr to allow any edge effects to die away, and verify that the density profile out to r₂₀₀ is completely stable. As before we define a dynamical time for the system of t_dyn ≡ r_max/v_max = 470 Myr.

After the 2 Gyr ≃ 4t_dyn has elapsed, we select all particles with j < 0.2 j_circ. These particles are, at the moment of selection, on preferentially radial orbits. We trace our particles forward through time, measuring 〈j〉(E) in each snapshot. The results are shown in Fig. 5 for various times between Δt = 0.0 and 4.0 t_dyn after selection. Over this period, the mean angular momentum significantly increases towards the SE limit at every energy. The changes are much faster at low energies where the particles are tightly bound and the local dynamical time is short compared to the globally defined t_dyn.

Figure 5.

The mean angular momentum of particles as a function of energy, like the lower panel of Fig. 4 – but now for a subpopulation in a dark matter halo that is globally in equilibrium. At Δt = 0, we select particles on predominantly radial (j/j_circ < 0.2) orbits; at later times, the subpopulation mean evolves back towards the population mean. The lines from bottom to top show the state at selection and after 0.2, 0.4, 0.6, 1.0 and 4.0 t_dyn, respectively.

Open in new tab Download slide

The evolution is rapid until 4 t_dyn after which the 〈j〉(E) remains near-static (except at E > 2.5v_max² where it continues to slowly rise). Eventually the subpopulation has 〈j〉(E) ≃ 0.5j_circ(E) independent of E. This establishes that particles at low angular momentum within a completely stable, unevolving halo, automatically evolve towards higher angular momentum. After a few dynamical times, their angular momentum becomes comparable to that of a randomly selected particle from the full population, although with a continuing slight radial bias.

One can understand this incompleteness in a number of equivalent ways. Within our formal picture, it arises from the fact that only certain |$\delta H_{\boldsymbol {n}}$| are non-zero in equation (12). More intuitively, the continuing process of subpopulation evolution is being driven by particles on box orbits that change their angular momentum at near-fixed energy – but not all particles are on such orbits, and so some memory of the initial selection persists. We explicitly verified that particles selected on circular orbits also mix back into the population; the global equilibrium is maintained by particles diffusing in both directions in the angular momentum space.

3.4 Observations of dwarf spheroidal galaxies

In the previous section, we established that subpopulations on initially radially biased orbits evolve towards velocity isotropy even if the global potential is stable. We explained this in terms of our earlier calculations, and now turn to how the results might impact on observations.

An important current astrophysical question is how dark matter is distributed within low-mass galaxies: this can discriminate between different particle physics scenarios (Pontzen & Governato 2014). The smallest known galaxies, dwarf spheroidals surrounding our Milky Way, in principle provide a unique laboratory from this perspective. But determining the distribution of dark matter is a degenerate problem because of the unknown transverse velocities of the stellar component (e.g. Charbonnier et al. 2011). Furthermore, the use of spherical analyses seems inappropriate since the underlying systems are known to be triaxial (e.g. Bonnivard et al. 2015).

This is exactly the kind of situation we set up in our initial calculations (Section 2): an aspherical system being analysed in spherical symmetry. We have shown that the results apply both to tracer and self-gravitating populations. So, we can go ahead and apply it to the stars in observed systems: the stellar population of a dwarf spheroidal system will be maximally stable if it appears ergodic in a spherical analysis. This implies a strong prior on what spherical distribution functions are actually acceptable and therefore, in principle, lessens the anisotropy degeneracy.

This idea warrants exploration in a separate paper; here we will briefly test whether the idea is feasible by inspecting some simulated dwarf spheroidal satellite systems. In particular we use a gas-dynamical simulation of ‘MW’ region (see above) using exactly the same resolution and physics as Zolotov et al. (2012); see that work for technical details (although note the actual simulation box is a different realization).

We analyse all satellites with more than 10⁴ star particles, which gives us objects lying in the ranges 5.4 × 10⁷ < M_⋆ < 4.9 × 10⁸ M_⊙, 28 < v_max < 60 km s⁻¹ and 4.6 < r_max < 9.6 kpc. We first verified that these do not host rotationally supported discs; we then calculated β(r) profiles between 3ϵ ≃ 0.5 kpc and 0.6 r_max. The stellar half-light radius is 0.15 < r_1/2/r_max < 0.3, so our calculations extend into the outer edges of each visible system.

The results are shown in the upper panel of Fig. 6, with arrows indicating r_1/2. The outer regions display a range of different behaviours from strongly radial to circular. However, in each case the β(r) is nearly isotropic as r → 0. This is consistent with a view in which the centres have achieved stability while the outer parts are being harassed by the parent tides and stripping.

Figure 6.

The β(r) (upper panel) and 〈j〉(E) (lower panel) relations for the stellar populations of dwarf satellites around a Milky Way-like central. In the upper panel arrows indicate the location of r_1/2, the radius enclosing half the stellar mass. The profile is plotted exterior to 3ϵ where ϵ is the softening length. Because of the tidal interaction with the parent galaxy, the outer parts of each object are out of equilibrium and display biases ranging from strongly radial to circular. The lower panel shows the same populations in energy space. At low energies, tightly bound particles are now seen to be close to the SE limit. At higher energies, a spread is seen but not as large as would naively be expected from the β(r) relation.

Open in new tab Download slide

The picture is reinforced by considering the 〈j〉(E) statistic (lower panel of Fig. 6). Here, stars on tightly bound orbits (small E) are very close to the SE limit. In the less-bound regions, the agreement worsens. However, a quantitative analysis using the approximate relation (21) shows that the 〈j〉(E) relation stays much closer to the SE limit than the β(r) relation naively implies. This is probably because the shape of β(r) is in part determined by out-of-equilibrium processes, coupled to the multivalued relationship between r and E (Fig. 3).

In any case these results suggest that we should adopt priors on spherical Jeans or Schwarzschild analyses that strongly favour near-isotropy in the centre of spheroidal systems – or, more accurately, strongly favour near-ergodicity for tightly bound stars. In future work we will expand on these ideas and apply them to observational data, since they could lead to a substantial weakening of the problematic density–anisotropy degeneracies.

3.5 Dark matter cusp–core transitions

In Pontzen & Governato (2012, hereafter PG12) we established that dark matter can be redistributed when intense, short, repeated bursts of star formation repeatedly clear the central regions of a forming dwarf galaxy of dense gas (see also Read & Gilmore 2005; Mashchenko, Couchman & Wadsley 2006). The resulting time-changing potential imparts net energy to the dark matter in accordance with an impulsive analytic approximation. This type of activity has now been seen or mimicked in a large number of simulations, allowing glimpses of the dependency of the process on galactic mass, feedback type and efficiency (Governato et al. 2012; Martizzi et al. 2012; Peñarrubia et al. 2012; Zolotov et al. 2012; Garrison-Kimmel et al. 2013; Teyssier et al. 2013; Di Cintio et al. 2014; Oñorbe et al. 2015).

We based our PG12 analysis on 3D zoom cosmological simulations, but our analytic model assumed exact spherical symmetry. By definition, in the analytic model all particles conserve their angular momentum at all times. Taken at face value, energy gains coupled to constant angular momentum would leave a radially biased population in the centre of the halo.

The top panel of Fig. 7 shows the measured velocity anisotropy in the inner 5 kpc for the zoom simulations in PG12 at z = 0; the solid line shows the feedback simulation which has developed a core, whereas the dotted line shows the dark-matter-only simulation which maintains its cusp. The difference between the two cases is the opposite to that naively expected: the centre of the cored halo has a more isotropic velocity dispersion than the centre of the cusped halo. The lower panel shows the equivalent picture in energy space [noting that Φ(0) = 0 and Φ(5 kpc) = 3 790 km² s⁻²]. For clarity only a small fraction of the 0 < j/j_circ < 1 interval is now plotted; the differences between the cusped (dotted) and cored (solid) lines are relatively small, but significant. At all energies the cored profile has more angular momentum than the cusped profile; in fact, it lies right on the SE limit (dashed line) whereas the cusped profile is biased by around 0.06 to radial orbits.¹

Figure 7.

The velocity anisotropy [β(r), upper panel] and angular momentum [〈j〉(E), lower panel] for cusped and cored haloes from PG12; note the much-expanded y-axis scales relative to previous figures, which are required to highlight the differences. The cored cases (solid lines) are almost perfectly SE (lower panel) and hence isotropic (upper panel). This contrasts with the cusped case (dotted lines) which has a slight but significant radial bias (seen as high β in the upper panel and low j in the lower panel).

Open in new tab Download slide

While the PG12 model correctly describes the energy shifts, it misses the angular momentum aspect of core-creation, which has been emphasized elsewhere (e.g. Tonini, Lapi & Salucci 2006). We will now show that a complete description of cusp–core transitions involves two components: energy gains consistent with PG12, and re-stabilization of the distribution function consistent with the description of Sections 2.2 and 3.3. We again use the ramses code to follow isolated haloes with particle mass 2 × 10⁴ M_⊙; however we adapted the code to add an external, time-varying potential to the self-gravity. The time-stepping is calculated adaptively for each particle based on the cell crossing and free-fall times and so is automatically adjusted in the modified potential.

We took our stabilized, isolated, triaxial dark matter halo extracted as described in Section 3.3 and created a sphericalized version of it as follows. For each particle we generate a random rotation matrix following an algorithm given by Kirk (1992), then multiply the velocity and position vector by this matrix. Finally, we verified that the final particles are distributed evenly in solid angle, and that the spherically averaged density profile and velocity anisotropy are unchanged. The triaxial and spherical haloes are both NFW-like and stable over more than 2 Gyr when no external potential is applied.

For our science runs, we impose an external potential corresponding to 10⁸ M_⊙ gas in a spherical ball 1 kpc in radius, distributed following ρ ∝ r⁻²; this implies a potential perturbation of ΔΦ = −700 km² s⁻² at 500 pc, for instance. The potential instantaneously switches off at 100 Myr, back on at 200 Myr, off at 300 Myr and so forth until it has accomplished four ‘bursts’. The period and the mass in gas are motivated by figs 1 and 2 of PG12 and fig. 7 of Teyssier et al. (2013). We also tried imposing potentials with different regular periods and with random fluctuations, none of which altered the behaviour described below.

The top two panels of Fig. 8 show the time-dependent behaviour of the central, most-bound 0.1 per cent of particles. The triaxial and spherical halo results are illustrated by solid and dashed lines, respectively. The top panel shows how the coupling of the external potential is similar; in particular, it results in a mean increase in specific energy of ΔE ≃ 200 km² s⁻² for particles in both cases. The final shift in the spherical case is slightly larger than that in the triaxial case, a difference which is unimportant for what follows (it would tend to create a larger core if anything).

Figure 8.

The PG12 mechanism is reproduced by adding an external potential (to mimic ‘gas’) to a DM-only simulation for the time intervals indicated by the grey bands. We measure the response of a completely spherical halo (dotted lines) and a realistically triaxial halo (solid lines) as described in the text. In both cases we track the mean energy (top panel) and angular momentum (middle panel) of the 0.1 per cent of particles that start out most bound. We also measure the log slope of the density profile at each time output (bottom panel). The energy shift (top) predicted by PG12 is found to hold in both cases. In the spherical case, the angular momentum remains constant up to discreteness effects; it therefore drops relative to the circular angular momentum (middle panel). Conversely, in the triaxial case, the continuous ROI causes the angular momentum to rise in proportion to j_circ. Only when the angular momentum is allowed to rise does a significant core form, meaning that the spherical case unphysically suppresses core formation (lower panel).

Open in new tab Download slide

The middle panel displays the mean specific angular momentum j of the same particles as a fraction of the circular angular momentum j_circ(E). In the spherical case (dashed line), this quantity drops significantly because j is fixed but j_circ is rising over time (since it increases with E). By contrast in the triaxial case j/j_circ returns to its original value, meaning that j has risen by the same fraction as j_circ(E). This is dynamic confirmation of the discussion earlier in this section: the spherically symmetric approach predicts an increasing bias to radial orbits – whereas in the realistic aspherical case, the stability requirements quickly erase this bias.

The bottom panel of Fig. 8 shows the measured density slope α = dln ρ/dln r at 500 pc for the two experiments. Both start at α ≃ −1.5; the triaxial case correctly develops a core (with α ≃ −0.1) whereas the completely spherical case maintains a cusp (with α ≃ −0.9). This discrepancy is causally connected to the angular momentum behaviour: the mean radius of a particle increases with increasing angular momentum j, even at fixed energy E, so a typical particle migrates further outwards in the triaxial case compared to the spherical.

We can therefore conclude that asphericity is a pre-requisite for efficient cusp–core transitions. However there is a subsidiary issue worth mentioning. After about 3 Gyr we find that the spherical halo autonomously does start increasing 〈j〉/j_circ and the dark matter density slope becomes shallower. This is because the potential fluctuations have generated a radially biased population which is unstable, and a global ROI is triggered by numerical noise over a sufficiently long period (as in Section 3.2).

In fact, with ‘live’ baryons rather than imposed fluctuations, there are aspherical perturbations which accelerate the re-equilibration process further and renew a global ROI (Section 3.2), encouraging the entire population towards spherical ergodicity. As we saw in Fig. 7, the cored halo from PG12 has an almost perfectly SE population for E < 3 000 km² s⁻² (unlike the cusped case). We verified that this is also true of Teyssier et al. (2013). Generating a core through potential fluctuations seems to complete a relaxation process that otherwise freezes out at an incomplete stage during collisionless collapse.

4 CONCLUSIONS

Let us return to the original question: how much do inaccuracies inherent in the spherical approximation really matter in practical situations? The answer is, unfortunately, that ‘it depends’; but we can now distinguish two cases as a rule of thumb.

For dynamical calculations or simulations, the inaccuracy matters a great deal. Neglecting the aspherical part of the potential unphysically freezes out the ROI and related effects, so can lead to qualitatively incorrect behaviour.
For the analysis of observations or simulations in equilibrium, the assumption is far more benign – it is, in fact, extremely powerful when handled with care. The underlying aspherical system and the fictional spherical system both appear to be in equilibrium; the mapping between the two views yields striking insight into (for example) spheroidal stellar distribution functions and dark matter halo equilibria.

These conclusions are based on the fact that, when aspherical systems are analysed in spherical coordinates, there is an attractor solution for the spherically averaged distribution function f₀ – namely, it tends towards being ergodic (i.e. f₀ is well-approximated as a function of energy alone). We demonstrated this using Hamiltonian perturbation theory (see Section 2 and Appendix B), and subsequently used the term ‘spherically ergodic’ to describe a distribution function f with the property that its spherical average f₀ is ergodic in this way.

The result follows because the orbits do not respect spherical integrals-of-motion such as angular momentum. Note, however, that the physical orbits may still possess invariants in a more appropriate set of coordinates. The apparent chaotic behaviour and tendency of particles to spread evenly at each energy is a helpful illusion caused by adopting coordinates that are only partially appropriate.

In Section 3 we showed that the idea has significant explanatory power. First, we inspected the selection of equilibrium in triaxial dark matter haloes (Section 3.1) which led us to consider the classical ROI (Section 3.2). We demonstrated that the instability naturally terminates very near our SE limit (lower panel, Fig. 4). Particles at high energies have long dynamical times which causes them to freeze out: they evolve towards, but do not reach, the SE limit. Because these high-energy particles often stray into the innermost regions (Fig. 3), the velocity anisotropy β(r) continues to display a significant radial bias at all r after the instability has frozen out. Moreover in the case of self-consistently formed cosmological haloes (Fig. 2), even at low energies there is a slight radial bias. The bias is only erased when a suitable external potential is applied (e.g. during baryonic cusp-core transformations), forcing the system to stabilize itself against a wider class of perturbations than it can self-consistently generate (Fig. 7); we will return to this issue momentarily.

One novel aspect of our analysis compared to previous treatments of the ROI is that it applies as much to tracer particles as to self-gravitating populations. As a first example, in Section 3.3 we demonstrated how a subset of dark matter particles chosen to be on radially biased orbits mix back into the population (Fig. 5). This is the case even for haloes that, as a whole, are in stable equilibrium; therefore we referred to the phenomenon as a ‘continuous’ ROI.

The same argument implies that stars undergo the ROI as easily as dark matter particles. In particular, without a stable disc to enforce extra invariants, dwarf spheroidal galaxies likely have stellar populations with a near-SE distribution (Section 3.4; Fig. 6). This provides a footing on which to base spherical Jeans or Schwarzschild analyses of observed systems: it implies an extra prior which can be formulated loosely as stating that β(r) → 0 as r → 0. Such a prior could be powerful in breaking the degeneracy between density estimates and anisotropy (e.g. Charbonnier et al. 2011), in turn tightening limits on the particle physics of dark matter (Pontzen & Governato 2014).

As a final application, we turned to the question of the baryonic processes that convert a dark matter cusp into a core (Section 3.5). Angular momentum is gained by individual particles during the cusp-flattening process (Tonini et al. 2006) but our earlier work (especially PG12) has focused primarily on the energy gains instead. In Fig. 8 we see that, in spherical or in triaxial dark matter haloes, time-dependent perturbations (corresponding to the behaviour of gas in the presence of bursty star formation feedback) always lead to a rise in energy. However, in the spherical case the mean angular momentum as a fraction of j_circ drops because the angular momentum of each particle is fixed, whereas j_circ rises with E. Only in the triaxial case does 〈j〉/j_circ stay constant, indicating a re-isotropization of the velocities. We tied the increase in 〈j〉 to the continuous ROI which always pushes a population with low 〈j〉/j_circ back towards the SE limit (Section 3.3). The consequence is that – in realistically triaxial haloes – the final dark matter core size will be chiefly determined by the total energy lost from baryons to dark matter, with little sensitivity to details of the coupling mechanism.

Although the PG12 model can be completed neatly in this way, other analytic calculations or simulations based on spherical symmetry will need to be considered on a case-by-case basis. In our exactly spherical test cases (Fig. 8), core development is substantially suppressed. This certainly cautions against taking the results of purely spherical analyses at face value. On the other hand, any slight asphericity is normally sufficient to prevent the unphysical angular-momentum lock-up – in particular, we verified that the simulations in Teyssier et al. (2013) were unaffected by this issue because their baryons settle into a flattened disc. The analytic result is generic, so the exact shape and strength of the asphericity is a secondary effect in determining the final spherically averaged distribution function f₀.

That said, to understand the way in which f₀ actually evolves towards the SE limit (and perhaps freezes out before it gets there) requires going to second order in perturbation theory, as shown in Appendix B. At this point it may also become important to incorporate self-gravity, i.e. the instantaneous connection between δf and δH. This approach has been investigated more fully elsewhere, leading to a different set of insights regarding the onset of the ROI (e.g. Saha 1991) as opposed to its end state. The present work actually suggests that the ROI can be cast largely as a kinematical process, and that the self-gravity is a secondary aspect; it would be interesting to further understand how these two views relate. But of more immediate practical importance is to apply the broader insights about dwarf spheroidal stellar equilibria to observational data, something that we will attempt in the near future.

All simulation analysis made use of the pynbody suite (Pontzen et al. 2013). We are grateful to the anonymous referee for helpful comments. AP acknowledges conversations with James Binney, Simon White, Chervin Laporte and Chiara Tonini as well as support from the Royal Society and, during 2013 when a significant part of this work was undertaken, the Oxford Martin School. JIR acknowledges support from SNF grant PP00P2_128540/1. FG acknowledges support from HST GO-1125, NSF AST-0908499. NR is supported by STFC and the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no 306478-CosmicDawn. JD's research is supported by Adrian Beecroft and the Oxford Martin School. This research used the DiRAC Facility, jointly funded by STFC and the Large Facilities Capital Fund of BIS.

1

The dashed line for the isotropic j/j_circ in the lower panel is calculated using the cored potential. However the dependence on potential is very weak, as discussed in Section 3; calculating it instead with the cusped potential leads to differences of less than one per cent at every energy.

REFERENCES

Adams

F. C.

Bloch

A. M.

Butler

S. C.

Druce

J. M.

Ketchum

J. A.

,

ApJ

,

2007

, vol.

670

pg.

1027

Crossref

Search ADS

Arnold

V. I.

,

Russian Math. Surv.

,

1963

, vol.

18

pg.

9

Crossref

Search ADS

Barnes

J.

Hut

P.

Goodman

J.

,

ApJ

,

1986

, vol.

300

pg.

112

Crossref

Search ADS

Barnes

E. I.

Lanzel

P. A.

Williams

L. L. R.

,

ApJ

,

2009

, vol.

704

pg.

372

Crossref

Search ADS

Bellovary

J. M.

Dalcanton

J. J.

Babul

A.

Quinn

T. R.

Maas

R. W.

Austin

C. G.

Williams

L. L. R.

Barnes

E. I.

,

ApJ

,

2008

, vol.

685

pg.

739

Crossref

Search ADS

Binney

J. J.

Tremaine

S.

,

Galactic Dynamics

,

2008

2nd edn

Princeton, NJ

Princeton University Press

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Bonnivard

V.

Combet

C.

Maurin

D.

Walker

M. G.

,

MNRAS

,

2015

, vol.

446

pg.

3002

Crossref

Search ADS

Bryan

S. E.

Mao

S.

Kay

S. T.

Schaye

J.

Dalla Vecchia

C.

Booth

C. M.

,

MNRAS

,

2012

, vol.

422

pg.

1863

Crossref

Search ADS

Charbonnier

A.

et al. ,

MNRAS

,

2011

, vol.

418

pg.

1526

Crossref

Search ADS

Cole

S.

Lacey

C.

,

MNRAS

,

1996

, vol.

281

pg.

716

Crossref

Search ADS

Corless

V. L.

King

L. J.

,

MNRAS

,

2007

, vol.

380

pg.

149

Crossref

Search ADS

de Zeeuw

T.

Merritt

D.

,

ApJ

,

1983

, vol.

267

pg.

571

Crossref

Search ADS

Dehnen

W.

,

MNRAS

,

1993

, vol.

265

pg.

250

Crossref

Search ADS

Di Cintio

A.

Brook

C. B.

Maccio

A. V.

Stinson

G. S.

Knebe

A.

Dutton

A. A.

Wadsley

J.

,

MNRAS

,

2014

, vol.

437

pg.

415

Crossref

Search ADS

Dubinski

J.

Carlberg

R. G.

,

ApJ

,

1991

, vol.

378

pg.

496

Crossref

Search ADS

Garrison-Kimmel

S.

Rocha

M.

Boylan-Kolchin

M.

Bullock

J.

Lally

J.

,

MNRAS

,

2013

, vol.

433

pg.

3539

Crossref

Search ADS

Goldstein

H.

Poole

C.

Safko

J.

,

Classical Mechanics

,

2002

3rd edn

San Francisco

Addison Wesley

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Governato

F.

et al. ,

MNRAS

,

2012

, vol.

422

pg.

1231

Crossref

Search ADS

Hayashi

E.

Navarro

J. F.

,

MNRAS

,

2006

, vol.

373

pg.

1117

Crossref

Search ADS

Hayashi

E.

Navarro

J. F.

Springel

V.

,

MNRAS

,

2007

, vol.

377

pg.

50

Crossref

Search ADS

Henon

M.

,

A&A

,

1973

, vol.

24

pg.

229

Hernquist

L.

,

ApJ

,

1990

, vol.

356

pg.

359

Crossref

Search ADS

Holley-Bockelmann

K.

Mihos

J. C.

Sigurdsson

S.

Hernquist

L.

,

ApJ

,

2001

, vol.

549

pg.

862

Crossref

Search ADS

Huss

A.

Jain

B.

Steinmetz

M.

,

ApJ

,

1999

, vol.

517

pg.

64

Crossref

Search ADS

Jing

Y. P.

Suto

Y.

,

ApJ

,

2002

, vol.

574

pg.

538

Crossref

Search ADS

Kasun

S. F.

Evrard

A. E.

,

ApJ

,

2005

, vol.

629

pg.

781

Crossref

Search ADS

Kirk

D.

,

Graphics Gems III. AP Professional

,

1992

MA

Cambridge

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Kolmogorov

A. N.

,

Dokl. Akad. Nauk SSSR

,

1954

, vol.

98

pg.

527

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Lichtenberg

A. J.

Lieberman

M. A.

,

Regular and Stochastic Motion

,

1992

2nd edn

New York

Springer

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Loebman

S. R.

Ivezić

Ž.

Quinn

T. R.

Governato

F.

Brooks

A. M.

Christensen

C. R.

Jurić

M.

,

ApJ

,

2012

, vol.

758

pg.

L23

Crossref

Search ADS

Lynden-Bell

D.

,

MNRAS

,

1967

, vol.

136

pg.

101

Crossref

Search ADS

MacMillan

J. D.

Widrow

L. M.

Henriksen

R. N.

,

ApJ

,

2006

, vol.

653

pg.

43

Crossref

Search ADS

Maréchal

L.

Perez

J.

,

Transp. Theory Stat. Phys.

,

2011

, vol.

40

pg.

425

Crossref

Search ADS

Martizzi

D.

Teyssier

R.

Moore

B.

Wentz

T.

,

MNRAS

,

2012

, vol.

422

pg.

3081

Crossref

Search ADS

Mashchenko

S.

Couchman

H. M. P.

Wadsley

J.

,

Nature

,

2006

, vol.

442

pg.

539

Crossref

Search ADS

PubMed

Merritt

D.

Poon

M. Y.

,

ApJ

,

2004

, vol.

606

pg.

788

Crossref

Search ADS

Merritt

D.

Valluri

M.

,

ApJ

,

1996

, vol.

471

pg.

82

Crossref

Search ADS

Meza

A.

Zamorano

N.

,

ApJ

,

1997

, vol.

490

pg.

136

Crossref

Search ADS

Moser

J.

,

Nachr. Akad. Wiss. Gttingen Math.-Phys. Kl. II

,

1962

pg.

1

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Navarro

J. F.

Frenk

C. S.

White

S. D. M.

,

ApJ

,

1996

, vol.

462

pg.

563

Crossref

Search ADS

Navarro

J. F.

Frenk

C. S.

White

S. D. M.

,

ApJ

,

1997

, vol.

490

pg.

493

Crossref

Search ADS

Navarro

J. F.

et al. ,

MNRAS

,

2010

, vol.

402

pg.

21

Crossref

Search ADS

Oñorbe

J.

Boylan-Kolchin

M.

Bullock

J. S.

Hopkins

P. F.

Kerěs

D.

Faucher-Giguère

C.-A.

Quataert

E.

Murray

N.

,

MNRAS

,

2015

submitted arXiv:1502.02036

Peñarrubia

J.

Pontzen

A.

Walker

M. G.

Koposov

S. E.

,

ApJ

,

2012

, vol.

759

pg.

L42

Crossref

Search ADS

Pontzen

A.

Governato

F.

,

MNRAS

,

2012

, vol.

421

pg.

3464

(PG12)

Crossref

Search ADS

Pontzen

A.

Governato

F.

,

MNRAS

,

2013

, vol.

430

pg.

121

Crossref

Search ADS

Pontzen

A.

Governato

F.

,

Nature

,

2014

, vol.

506

pg.

171

Crossref

Search ADS

PubMed

Pontzen

A.

Roskar

R.

Stinson

G.

Woods

R.

,

pynbody: N-Body/SPH Analysis for Python. Astrophysics Source Code Library ascl:1305.002

,

2013

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Read

J. I.

Gilmore

G.

,

MNRAS

,

2005

, vol.

356

pg.

107

Crossref

Search ADS

Read

J. I.

Wilkinson

M. I.

Evans

N. W.

Gilmore

G.

Kleyna

J. T.

,

MNRAS

,

2006

, vol.

367

pg.

387

Crossref

Search ADS

Saha

P.

,

MNRAS

,

1991

, vol.

248

pg.

494

Crossref

Search ADS

Saha

P.

Read

J. I.

,

ApJ

,

2009

, vol.

690

pg.

154

Crossref

Search ADS

Schneider

M. D.

Frenk

C. S.

Cole

S.

,

J. Cosmol. Astropart. Phys.

,

2012

, vol.

5

pg.

30

Crossref

Search ADS

Stadel

J.

Potter

D.

Moore

B.

Diemand

J.

Madau

P.

Zemp

M.

Kuhlen

M.

Quilis

V.

,

MNRAS

,

2009

, vol.

398

pg.

L21

Crossref

Search ADS

Taylor

J. E.

Navarro

J. F.

,

ApJ

,

2001

, vol.

563

pg.

483

Crossref

Search ADS

Teyssier

R.

,

A&A

,

2002

, vol.

385

pg.

337

Crossref

Search ADS

Teyssier

R.

Pontzen

A.

Dubois

Y.

Read

J. I.

,

MNRAS

,

2013

, vol.

429

pg.

3068

Crossref

Search ADS

Tonini

C.

Lapi

A.

Salucci

P.

,

ApJ

,

2006

, vol.

649

pg.

591

Crossref

Search ADS

Valluri

M.

Debattista

V. P.

Quinn

T.

Moore

B.

,

MNRAS

,

2010

, vol.

403

pg.

525

Crossref

Search ADS

van Albada

T. S.

,

MNRAS

,

1982

, vol.

201

pg.

939

Crossref

Search ADS

Weinberg

M. D.

,

ApJ

,

1991

, vol.

368

pg.

66

Crossref

Search ADS

Zolotov

A.

et al. ,

ApJ

,

2012

, vol.

761

pg.

71

Crossref

Search ADS

Zorzi

A. F.

Muzzio

J. C.

,

MNRAS

,

2012

, vol.

423

pg.

1955

Crossref

Search ADS

APPENDIX A: BACKGROUND

A1 A brief review of actions

This Appendix contains a very brief review of actions which are necessary for deriving the main result in the paper. For more complete introductions see Binney & Tremaine (2008) or Goldstein et al. (2002). We consider any mechanical problem described by phase-space coordinates |$\boldsymbol {q}$| and momenta |$\boldsymbol {p}$|⁠, with Hamiltonian |$H(\boldsymbol {p},\boldsymbol {q})$| so that the equations of motion are

\begin{equation} \dot{\boldsymbol {p}} = -\frac{\mathrm{\partial} H}{\mathrm{\partial} \boldsymbol {q}}\textrm {,} \hspace{28.45274pt} \dot{\boldsymbol {q}} = \frac{\mathrm{\partial} H}{\mathrm{\partial} \boldsymbol {p}}\textrm {,} \end{equation}

(A1)

where |$\dot{\boldsymbol {q}}=\mathrm{d}\boldsymbol {q}/\mathrm{d}t$| and t denotes time.

The actions can be defined starting from any coordinate system in which the motion is separately periodic in each dimension (i.e. for each i, q_i and p_i repeat every Δt_i). The actions J_i are then given by

\begin{equation} J_i \equiv \frac{1}{2\pi }\int _0^{\Delta t_i} p_i \dot{q}_i \,\mathrm{d}t\textrm {(no sum over}\; i. \end{equation}

(A2)

By construction the actions J_i do not change under time evolution and are therefore integrals of motion.

We can complete the set of six phase-space coordinates with angles |$\boldsymbol {\Theta }$| in such a way that equations of motion of the form (A1) apply. Since J_i is constant, we must then have

\begin{equation} 0 = \dot{J}_i=-\frac{\mathrm{\partial} H}{\mathrm{\partial} \Theta _i}\textrm {,}\hspace{28.45274pt} \dot{\Theta }_i = \frac{\mathrm{\partial} H}{\mathrm{\partial} J_i} \equiv \Omega _{\rm i}(\boldsymbol {J})\textrm {,} \end{equation}

(A3)

where the first equation establishes that H can have no Θ dependence, and the second that consequently the Θ_i each increase at a constant rate in time specified by the frequency |$\Omega _{\rm i}(\boldsymbol {J})$|⁠. The convenience of this set of coordinates is that all the time evolution of a particle trajectory is represented in a very simple way:

\begin{equation} J_i = \textrm {constant,} \hspace{28.45274pt} \Theta _i = \textrm {constant} + \Omega _{\rm i}\, t\textrm {.} \end{equation}

(A4)

Furthermore the equations of motion are canonical, which is sufficient to demonstrate that the coordinates are canonical – in other words, the measure appearing in phase space integrals is d³J d³Θ.

We now specialize to the spherical case, with polar coordinates |$\boldsymbol {q}=(r,\theta ,\phi )$| and momenta |$\boldsymbol {p}=(\dot{r},r^2 \dot{\theta }, r^2 \sin ^2 \theta \dot{\phi })$|⁠. Consider first J_ϕ,

\begin{equation} J_{\phi } = \frac{1}{2\pi } \int _0^{\Delta t_{\phi }} \mathrm{d}t \dot{\phi }^2 r^2 \sin ^2 \theta = \frac{1}{2\pi} \int _0^{2\pi } \mathrm{d}\phi \dot{\phi } r^2 \sin ^2 \theta = j_z, \!\!\!\!\! \end{equation}

(A5)

where j_z is the z component of the specific angular momentum, |$j_z=\dot{\phi } r^2 \sin ^2 \theta$| which is a constant of motion. Next consider J_θ,

\begin{equation} J_{\theta } = \frac{1}{2\pi } \int _0^{\Delta t_{\theta }} \mathrm{d}t \dot{\theta }^2 r^2 = \frac{1}{\pi } \int _{\theta _{\rm a}}^{\theta _b} \mathrm{d}\theta \sqrt{j^2 - \frac{j_z^2}{\sin ^2 \theta }} \textrm {,} \end{equation}

(A6)

where |$j^2=j_x^2+j_y^2+j_z^2$| is the square of the total specific angular momentum, and θ varies between θ_a and θ_b over the course of an orbit. To evaluate the integral requires a change of variables to relate θ_a and θ_b to the inclination of the orbit and so to j and j_z (e.g. Goldstein et al. 2002, equation 10.135); the final result is that

\begin{equation} J_{\theta } = j-j_z\textrm {.} \end{equation}

(A7)

Because linear combinations of actions are still actions (i.e. they still satisfy equation (A3)) one can take j_z as the first and j as the second action in place of J_ϕ and J_θ.

The last action is the radial action,

\begin{equation} J_{\rm r} = \frac{1}{2 \pi } \int _0^{\Delta t_{\rm r}} \mathrm{d}t \dot{r}^2 = \frac{1}{ \pi} \int _{r_{\rm a}}^{r_b} \mathrm{d}r \sqrt{2E - j^2/r^2 - 2 \Phi (r)}\textrm {,} \end{equation}

(A8)

where |$\dot{r}$| has been evaluated by energy conservation, and r librates between r_a and r_b over the period of an orbit.

Although the complexity of expression (A8) appears to make using the actions cumbersome, the great simplification it brings to the equations of motion in the background (i.e. equation A3) makes the perturbation theory tractable. For that reason we have used the action-angle coordinates in our analytic derivation, Section 2, but avoided them when discussing results from simulations in Section 3.

A2 Perturbed trajectories in the harmonic oscillator case

Section 2 used Hamiltonian perturbation theory to discuss the behaviour of particles in aspherical potentials. To connect this more firmly with the equations of motion, it can be helpful to study orbits in a specific potential and connect the solutions with the more general statements made by the perturbation theory. In this Appendix, we use the harmonic oscillator as such an illustrative example.

Consider the Hamiltonian for a single particle in an anisotropic harmonic oscillator potential,

\begin{equation} H = \frac{1}{2} \left(\pi _x^2 + \pi _y^2 + \pi _z^2\right) + \frac{\omega _0^2}{2} ( x^2+(1+\epsilon ) y^2+(1-\delta ) z^2)\;,\!\!\!\!\!\! \end{equation}

(A9)

where π_x, π_y, π_z denote the momentum in the three Cartesian directions and without loss of generality ϵ ≥ 0 and δ ≥ 0.

We want to compare the motion in this potential to the behaviour in a sphericalized version with Hamiltonian H₀; the change in the Hamiltonian between the true and sphericalized cases is given by

\begin{equation} \delta H = H-H_0 = \frac{1}{2} \omega _0^2 (\epsilon y^2 - \delta z^2 )\textrm {.} \end{equation}

(A10)

We can immediately read off the first result, which is that the magnitude of δH is bounded. The true solution moves around on a fixed H surface, meaning the fractional error takes a maximum value given by

\begin{equation} \frac{\left| \delta H \right|}{H}< \max (\epsilon ,\delta )\textrm {.} \end{equation}

(A11)

This is the equivalent of the general statement that H₀ variations are small, given by equation (8).

On the other hand, the angular momentum does change significantly. We can see this as follows: each of x, y and z undergoes oscillation at the frequencies ω₀, (1 + ϵ)^1/2ω₀ and (1 − δ)^1/2ω₀, respectively. Assuming these new frequencies are not commensurate, the relative phases between the different oscillations slowly shift until at some point all three separated oscillators reach x = 0, y = 0 and z = 0 at the same moment. At this point, since the velocity remains finite, the angular momentum has become zero. It may take a number of dynamical times before this happens, but (for example) at the centre of dark matter haloes the dynamical time is very short compared to the Hubble time so angular momentum conservation is effectively destroyed.

All the above is illustrated in Fig. A1. The left panel contrasts the orbits for a spherical harmonic oscillator (dashed line) and an aspherical oscillator (solid line) projected in the (x, y) plane. The latter obeys equation (A9) with ϵ = δ = 0.1 (and the former with ϵ = δ = 0). The orbits for the spherical case are closed because the frequencies are identical, so the relative phase of the x and y part of the motion remains fixed. Orbits in the aspherical potential are more complex as the relative phase of the Cartesian components gradually changes; in fact, a particle will sometimes plunge through the centre of the potential. This is known as a ‘box orbit’. In more general triaxial potentials, a variety of orbit types are possible (e.g. Merritt & Valluri 1996); the importance of these will be considered momentarily.

Figure A1.

An example orbit in a 3D anisotropic harmonic oscillator. Left: the solid line shows a portion of the orbit projected into the x-y plane. The equivalent orbit in an isotropic potential is closed and given by the blue dashed line. Right: the orbit analysed in terms of J_r and j, the actions for the spherical background Hamiltonian. Because ΔH is bounded (see text) the particle has to stay within a narrow band in the J_r-j plane, corresponding to the linear theory resonance (dotted line). The blue dot shows the action for the orbit in the spherical potential (for which the actions are by definition fixed).

Open in new tab Download slide

The same portion of the orbit is illustrated in the right panel, but now projected into the spherical actions plane (J_r, j). For the spherical case, J_r and j are exactly conserved by construction and the orbit appears as a single point. In the aspherical case (solid line), J_r and j are not even approximately conserved over more than a dynamical time. However, even then, the orbit remains close to the dashed line of constant H₀, as required by equation (A11) and more generally by equation (8). From this figure, it is intuitively plausible that the particle is equally likely to be found anywhere along the constant H₀ contour, which is the essential result of Section 2.2.

Another way to look at the effect is to consider the relationship between the angular momentum and the harmonic oscillator's Cartesian actions, J_x, J_y and J_z which remain constant even in the aspherical case. Specializing for simplicity to the case that J_z = 0, one can show that

\begin{equation} j = 2 \sqrt{J_x J_y} \sin (\Theta _x-\Theta _y)\textrm {,} \end{equation}

(A12)

where Θ_x and Θ_y are the angles conjugate to J_x and J_y. Only if Θ_x − Θ_y remains constant is j an integral of motion; as soon as the oscillator is aspherical, one has

\begin{equation} j = 2 \sqrt{J_x J_y} \sin ( \phi _0 + (\Omega _{\rm x} - \Omega _{\rm y}) t) \end{equation}

(A13)

so that j oscillates on the time-scale 2π(Ω_x − Ω_y)⁻¹ ≃ πϵ⁻¹. The situation in 3D is qualitatively similar.

The harmonic oscillator orbits discussed above are all box orbits. More general triaxial potentials support other orbit types (loop orbits, for example, which are much more tightly constrained; or chaotic orbits, which are even less tightly constrained than box orbits). However the Hamiltonian analysis in the main text is general for all these types of possible orbit. The fraction of different orbit types will determine how fast and how far particles diffuse along the contour. In realistic cosmological dark matter haloes, most orbits in the central regions are indeed of the box or chaotic type (Zorzi & Muzzio 2012). Even with the baryonic contribution to the gravitational force included (which partially sphericalizes the central potential), the large majority of particles remain on the same class of orbit as the dark progenitor (Valluri et al. 2010) so long as feedback is strong enough to prevent long-lived central baryon concentrations developing (Bryan et al. 2012).

A3 The action space of dark matter haloes

The action-angle space of our equilibrium ‘Dwarf’ dark-matter-only halo is illustrated in Fig. A2, along with some particle orbits (solid lines) over 1.5 Gyr ≃ 3 t_dyn. The Hamiltonian is a function of J_r and j only for any spherical potential, and so we have suppressed the third action j_z. As expected, the particles explore the space, approximately running along the H₀ contours, giving rise to the continuous ROI described in Section 3.3. The contours of H₀ give a great deal of dynamical information, because the background frequencies are defined by |$\boldsymbol {\Omega }_0 \equiv \mathrm{\partial} H_0 / \mathrm{\partial} \boldsymbol {J}$|⁠. These frequencies are at the core of perturbation theory through equation (6), but they also determine the higher-order behaviour as follows.

$Orbits in the action space of the equilibrium ‘Dwarf’ cosmological halo. Dotted contours of constant spherical specific energy H0 are spaced at 600 km2 s−2. (For this halo $v_{\mathrm{max}}^2 \simeq 3\,100\,\mathrm{km}^2\,\mathrm{s}^{-2}$.) The straightness of these contours is part of the reason that orbits can efficiently explore the space, as described in the text.$

Figure A2.

Orbits in the action space of the equilibrium ‘Dwarf’ cosmological halo. Dotted contours of constant spherical specific energy H₀ are spaced at 600 km² s⁻². (For this halo |$v_{\mathrm{max}}^2 \simeq 3\,100\,\mathrm{km}^2\,\mathrm{s}^{-2}$|⁠.) The straightness of these contours is part of the reason that orbits can efficiently explore the space, as described in the text.

Open in new tab Download slide

The first step to understanding the behaviour of resonances is to move to secular perturbation theory (e.g. section 2.4 of Lichtenberg & Lieberman 1992). Secular analysis splits resonant orbits into two classes, known as accidentally and intrinsically degenerate. The intrinsic case refers to the situation where the resonance condition applies globally; in other words, that |$\boldsymbol {n} \cdot \boldsymbol {\Omega }_0$| is near-constant along lines of constant H₀. Suppose, conversely, that |$\boldsymbol {\Omega }_0$| did vary along these directions; then, since |$\boldsymbol {\Omega }_0$| is defined by the normal to the H₀ contours (⁠|$\boldsymbol {\Omega }_0 = \mathrm{\partial} H_0/\mathrm{\partial} \boldsymbol {J}$|⁠), this would imply a significant curvature of the dotted lines in Fig. A2. Since there is no such curvature, we can read off that the frequencies do not change and the dynamics is in the intrinsically degenerate regime, giving rise to large-scale migrations. [The same property also means that the approximation Ω_r = Ω_r(H₀) used in reaching equation (18) will be extremely accurate.] We have established using numerical investigations that this intrinsic degeneracy property is generic to any spherical action space with smooth potentials, rather than being specific to cosmological haloes.

APPENDIX B: THE EVOLUTION OF f₀

In the main text (Section 2), we showed that a distribution function f is maximally stable to linear perturbations if its sphericalized part f₀ appears ergodic, f₀ = f₀(H₀). However we did not discuss the actual evolution of f₀ to see whether this limit is likely to be achieved.

This requires time-dependent perturbation theory. We start by writing the collisionless Boltzmann equation,

\begin{eqnarray} \frac{\mathrm{\partial} f}{\mathrm{\partial} t} = [H, f] = [H_0, f_0] + [H_0, \delta f] + [\delta H, f_0] + [\delta H, \delta f]\textrm {,} \nonumber\\ \!\!\!\!\!\!\! \end{eqnarray}

(B1)

which is an exact expression. The first term vanishes identically; the second and third terms are linear order, and the final term is second order. The evolution of f₀ is given by taking the time derivative of equation (9) and interchanging the derivative and integral:

\begin{eqnarray}\frac{\mathrm{\partial} f_0}{\mathrm{\partial} t} &=& \frac{1}{(2\pi )^3}\int \mathrm{d}^3 \, \Theta \, [H,f] \nonumber \\ & =& \frac{1}{(2\pi )^3}\sum _{\boldsymbol {m},\boldsymbol {n}} \int \mathrm{d}^3 \Theta \, \left[\delta H_{\boldsymbol {m}}\,{\rm e}^{{\rm i} \boldsymbol {m} \cdot \boldsymbol {\Theta }},\, \delta f_{\boldsymbol {n}} \,{\rm e}^{{\rm i} \boldsymbol {n} \cdot \boldsymbol {\Theta }} \right]\textrm {,} \end{eqnarray}

(B2)

where the two linear-order terms have vanished after integrating over Θ. Expanding the Poisson bracket and integrating the remaining term give

\begin{equation} \frac{\mathrm{\partial} f_0}{\mathrm{\partial} t} = - {\rm i} \sum _{\boldsymbol {n}} \boldsymbol {n} \cdot \frac{\mathrm{\partial} }{\mathrm{\partial} \boldsymbol {J}} \left( \delta f_{\boldsymbol {n}} \delta H_{-\boldsymbol {n}} \right)\textrm {.} \end{equation}

(B3)

This shows that the evolution of f₀ is a fundamentally non-linear phenomenon. To make further progress, we can eliminate |$\delta f_{\boldsymbol {n}}$|⁠, showing that f₀ evolution depends only on |$\delta H_{\boldsymbol {n}}$|⁠. First, Fourier transform the time-dependence of f₀ and δf, so that

\begin{eqnarray} f_0(\boldsymbol {J},t) & =& \int \mathrm{d}\omega \, {\rm e}^{{\rm i} \omega t}\, \tilde{f}_0(\boldsymbol {J}, \omega )\hspace{11.38092pt} \textrm { and} \end{eqnarray}

(B4)

\begin{eqnarray} \delta f(\boldsymbol {J}, \boldsymbol {\Theta }, t) & =& \sum _{\boldsymbol {n}} \int \mathrm{d}\omega \delta \tilde{f}_{\boldsymbol {n}} (\boldsymbol {J}, \omega )\,{\rm e}^{{\rm i} \omega t + {\rm i} \,\boldsymbol {\Theta } \cdot \boldsymbol {n}}\textrm {.} \end{eqnarray}

(B5)

For simplicity, we will restrict attention to the case where δH is constant in time. Then the evolution of δf is given by the linear-order part of equation (B1), Fourier transformed to give

\begin{equation} \left( \omega + \boldsymbol {\Omega }_0 \cdot \boldsymbol {n} \right) \delta \tilde{f}_{\boldsymbol {n}} = \left(\frac{\mathrm{\partial} \tilde{f}_0}{\mathrm{\partial} \boldsymbol {J}} \cdot \boldsymbol {n}\right) \delta H_{\boldsymbol {n}}\textrm {.} \end{equation}

(B6)

This is just the time-dependent, Fourier-transformed version of equation (12). Note that, for δH and δf to be real, the Fourier coefficients must satisfy

\begin{equation} \delta \tilde{f}_{\boldsymbol {n}}(\omega )^* = \delta \tilde{f}_{-\boldsymbol {n}}(-\omega ) \textrm { and } \delta H_{\boldsymbol {n}}^* = \delta H_{-\boldsymbol {n}}\textrm {.} \end{equation}

(B7)

These requirements are consistent with the relation (B6).

We can now substitute the linear-order solution (B6) for |$\delta f_{\boldsymbol {n}}$| in equation (B3) to get the leading-order evolution equation for f₀:

\begin{equation} \omega \tilde{f}_0 = -\sum _{\boldsymbol {n}} \boldsymbol {n} \cdot \frac{\mathrm{\partial} }{\mathrm{\partial} \boldsymbol {J}} \left[ \frac{\delta H_{\boldsymbol {n}} \delta H_{\boldsymbol {-n}} \boldsymbol {n} \cdot \mathrm{\partial} \tilde{f}_0/\mathrm{\partial} \boldsymbol {J}}{\omega + \boldsymbol {\Omega }_0 \cdot \boldsymbol {n}} \right]\textrm {.} \end{equation}

(B8)

Using the reality condition, equation (B7), we can pair up negative and positive |$\boldsymbol {n}$| modes, giving an alternative version of the expression that is more explicitly symmetric:

\begin{equation} \tilde{f}_0 = \sum _{\boldsymbol {n}} \boldsymbol {n} \cdot \frac{\mathrm{\partial} }{\mathrm{\partial} \boldsymbol {J}} \left[ \frac{\left| \delta H_{\boldsymbol {n}} \right|^2 \boldsymbol {n} \cdot \mathrm{\partial} \tilde{f}_0/\mathrm{\partial} \boldsymbol {J} }{\left| \boldsymbol {\Omega }_0 \cdot \boldsymbol {n} \right|^2 - \omega ^2} \right]\textrm {,} \end{equation}

(B9)

where we have divided both sides by ω and so the result is technically only applicable for ω ≠ 0. Provided the evolution of f₀ is smooth, its Fourier transform |$\tilde{f}_0$| is continuous, so this is not a problem in practice.

Equation (B9) is enough to give some insight into the relaxation process. For simplicity, consider a perturbation consisting of a single, resonant |$\boldsymbol {n}_{\perp }$|-mode |$\delta H_{\boldsymbol {n}_{\perp }}$| with |$\boldsymbol {n}_{\perp } \cdot \boldsymbol {\Omega }_0=0$|⁠. Then we can explicitly invert the Fourier transform to yield

\begin{equation} \frac{\mathrm{\partial} ^2 f_0}{\mathrm{\partial} t^2} = \boldsymbol {n}_{\perp } \cdot \frac{\mathrm{\partial} }{\mathrm{\partial} \boldsymbol {J}} [ \left| \delta H_{\boldsymbol {n}_{\perp }} \right|^2 \boldsymbol {n}_{\perp } \cdot \mathrm{\partial} f_0/\mathrm{\partial} \boldsymbol {J} ]\textrm {.} \end{equation}

(B10)

This is a wave equation in the |$\boldsymbol {n}_{\perp }$| direction with varying speed-of-sound proportional to |$|\delta H_{\boldsymbol {n}_{\perp }}|$|⁠. Any unevenness in the resonant directions will flow away at a speed proportional to the strength of the asphericity, showing explicitly that |$\tilde{f}_0$| evolves towards the SE limit. We hope to investigate the full behaviour of equation (B9) for a variety of regimes in future work.

Download all slides

Month:	Total Views:
January 2017	3
February 2017	1
March 2017	7
June 2017	1
August 2017	3
September 2017	1
October 2017	3
November 2017	1
December 2017	1
February 2018	2
March 2018	2
April 2018	2
May 2018	3
June 2018	1
July 2018	4
August 2018	8
September 2018	8
October 2018	15
November 2018	3
December 2018	13
January 2019	9
February 2019	14
March 2019	12
April 2019	17
May 2019	15
June 2019	10
July 2019	10
August 2019	6
September 2019	2
October 2019	9
November 2019	11
December 2019	3
January 2020	26
February 2020	9
March 2020	4
April 2020	2
May 2020	3
June 2020	2
July 2020	4
August 2020	6
September 2020	25
October 2020	3
November 2020	10
December 2020	7
January 2021	7
February 2021	7
March 2021	16
April 2021	6
May 2021	8
June 2021	1
July 2021	9
August 2021	7
September 2021	9
October 2021	9
November 2021	7
December 2021	6
January 2022	10
February 2022	9
March 2022	5
April 2022	2
May 2022	2
June 2022	6
July 2022	10
August 2022	10
September 2022	9
October 2022	6
November 2022	8
December 2022	19
January 2023	42
February 2023	12
March 2023	20
April 2023	18
May 2023	12
June 2023	21
July 2023	8
August 2023	28
September 2023	31
October 2023	18
November 2023	20
December 2023	30
January 2024	40
February 2024	25
March 2024	14
April 2024	23
May 2024	21
June 2024	19
July 2024	12
August 2024	21
September 2024	45
October 2024	24
November 2024	14
December 2024	14
January 2025	17
February 2025	30
March 2025	37
April 2025	19
May 2025	7

Article Contents

Milking the spherical cow – on aspherical dynamics in spherical coordinates

Abstract

1 INTRODUCTION

2 ASPHERICAL DYNAMICS IN SPHERICAL COORDINATES

2.1 Single particles

2.2 The distribution function

2.3 Testable predictions

3 EXAMPLE APPLICATIONS

3.1 Cosmological dark matter haloes

3.2 The classical radial orbit instability

3.3 The continuous radial orbit instability

3.4 Observations of dwarf spheroidal galaxies

3.5 Dark matter cusp–core transitions

4 CONCLUSIONS

REFERENCES

APPENDIX A: BACKGROUND

A1 A brief review of actions

A2 Perturbed trajectories in the harmonic oscillator case

A3 The action space of dark matter haloes

APPENDIX B: THE EVOLUTION OF f₀

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Milking the spherical cow – on aspherical dynamics in spherical coordinates

Abstract

1 INTRODUCTION

2 ASPHERICAL DYNAMICS IN SPHERICAL COORDINATES

2.1 Single particles

2.2 The distribution function

2.3 Testable predictions

3 EXAMPLE APPLICATIONS

3.1 Cosmological dark matter haloes

3.2 The classical radial orbit instability

3.3 The continuous radial orbit instability

3.4 Observations of dwarf spheroidal galaxies

3.5 Dark matter cusp–core transitions

4 CONCLUSIONS

REFERENCES

APPENDIX A: BACKGROUND

A1 A brief review of actions

A2 Perturbed trajectories in the harmonic oscillator case

A3 The action space of dark matter haloes

APPENDIX B: THE EVOLUTION OF f0

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

APPENDIX B: THE EVOLUTION OF f₀