Conserved actions, maximum entropy and dark matter haloes

Abstract

We use maximum entropy arguments to derive the phase-space distribution of a virialized dark matter halo. Our distribution function gives an improved representation of the end product of violent relaxation. This is achieved by incorporating physically motivated dynamical constraints (specifically on orbital actions) which prevent arbitrary redistribution of energy.

We compare the predictions with three high-resolution dark matter simulations of widely varying mass. The numerical distribution function is accurately predicted by our argument, producing an excellent match for the vast majority of particles.

The remaining particles constitute the central cusp of the halo (≲4 per cent of the dark matter). They can be accounted for within the presented framework once the short dynamical time-scales of the centre are taken into account.

galaxies: haloes, dark matter

1 INTRODUCTION

For over two decades it has been possible to use numerical methods to model systems of cold, collisionless dark matter particles collapsing under gravity to form stable, virialized ‘haloes’ (Frenk et al. 1985; Dubinski & Carlberg 1991; for a review see Frenk & White 2012). The density of these haloes declines with radius following a slowly changing power-law dependence, roughly ρ ∼ r⁻¹ at small radii and ρ ∼ r⁻³ in the outer regions (Navarro, Frenk & White 1996b; Kravtsov, Klypin & Khokhlov 1997; Moore et al. 1998). Despite some early uncertainty, recent simulations with independent computer codes all reproduce this result (e.g. Diemand et al. 2008; Stadel et al. 2009; Navarro et al. 2010). The universal behaviour seems to be independent of the power spectrum of initial linear density fluctuations (Moore et al. 1999; Reed et al. 2005; Wang & White 2009) as well as the mass of the collapsed object and the epoch of collapse, although together these determine a scale radius for the transition from r⁻¹ to r⁻³ behaviour (Cole & Lacey 1996; Navarro, Frenk & White 1997; Bullock et al. 2001b; Eke, Navarro & Steinmetz 2001; Macciò et al. 2007). Even simulations of ‘cold collapse’, for which initial conditions consist of a homogeneous sphere of particles, seem to produce similar universal profiles (Huss, Jain & Steinmetz 1999).

It remains an outstanding question, however, whether this universality can be adequately explained from first principles. Until this question is answered, we do not fully understand what the universality means and must rely on new simulations to predict the effect of changes in the initial conditions or particle properties. The experimental result that monolithic collapse produces the same types of system as hierarchical merging (Huss et al. 1999; Moore et al. 1999; Wang & White 2009) is provocative: it means that any explanation for universality which invokes a specific cosmology (e.g. Syer & White 1998; Dekel et al. 2003; Salvador-Solé et al. 2012) must be describing a special case of a more general process (Manrique et al. 2003).

Attempts to understand collisionless gravitational collapse’s insensitivity to initial conditions were pioneered by Lynden-Bell (1967) in the context of self-gravitating stellar systems. Adopting Boltzmann’s procedure for deriving the thermodynamics of collisional systems, Lynden-Bell maximized the entropy of the systems subject to fixed energy. This implies a density profile obeying ρ(r) ∝ r⁻², so disagrees with the results of numerical experiments. Moreover this approach gives rise to a number of physically questionable conclusions (for reviews see Padmanabhan 1990; Lynden-Bell 1999). The clearest of these is the ‘gravothermal catastrophe’: entropy can be increased without bound by transferring energy from the innermost orbits of a self-gravitating system to the outermost orbits (Lynden-Bell & Wood 1968; Tremaine, Henon & Lynden-Bell 1986). This implies the existence of a runaway physical instability in which the majority of material collapses into an extreme central density cusp or black hole.

Observations and numerics both suggest that something prevents the above catastrophe from occurring on any reasonable time-scale. As a result, alterations to normal statistical mechanical arguments have been proposed (e.g. Tsallis 1988; Hjorth & Williams 2010). However one way to look at the situation is that, in reality, a physical constraint must prevent the arbitrary redistribution of energy (White & Narayan 1987). In the spirit of Jaynes (1957) a general explanation for the stable state can then be based on the standard formulation in the presence of constraints other than energy (Stiavelli & Bertin 1987). The additional constraints will represent the incompleteness of the energy equilibriation effects of violent relaxation.

This is the approach we adopt in the present work. The likely distribution of particles in phase space is selected by maximizing entropy:

\begin{equation} S= -k \int f \ln f\, \mathrm{d}^6 \omega , \end{equation}

(1)

where k is Boltzmann’s constant and f (ω) is the probability of finding a particle in a specified region of phase space ω, subject to relevant constraints on the desired solution. (This approach can be motivated by showing that the vast majority of states consistent with a given set of constraints are to be found near the maximum entropy solution; see Appendix A for further discussion and references.) Although for convenience we will repeatedly refer to f as describing probabilities for finding particles in particular states, we do not need to specify the mass or number of these particles since the arguments are based on purely collisionless dynamics.

The constraints applied arise from the dynamical evolution of collisionless systems. We will argue that, in the late stages of violent relaxation, there is a diffusion of particles in phase space which approximately conserves the sum of orbital actions. This sum is progressively better conserved as equilibrium is approached and therefore its role in establishing that equilibrium cannot be ignored. Applying the maximum entropy recipe, we will show that the phase-space structure of equilibrium haloes is then reproduced over orders of magnitude in probability density. To our knowledge, this is the first instance of maximum entropy reasoning, applied to 6D phase space and subject to physically motivated constraints, producing such success in a collisionless system.

The remainder of this work is organized as follows. First, the conservation of action is discussed (Section 2.1). Then a canonical ensemble constructed on this basis (Section 2.3) yields a phase-space distribution in quantitative agreement with high-resolution numerical experiments (Section 3). A discrepancy affecting a small fraction of particles at low angular momentum is highlighted in Section 3.3. Finally, we discuss the predicted radial density profiles which again highlight the need for special treatment of low angular momentum orbits (Section 3.5). We conclude in Section 4.

2 THE ANALYTIC PHASE-SPACE DISTRIBUTION

2.1 Conservation of action

Our first task is to identify and explain some relevant quantities which should be held fixed when maximizing entropy. This will be central to our argument because such constraints represent the incompleteness of violent relaxation’s tendency to redistribute energy, generating a different solution from the one based on energy conservation alone. With this in mind we will show that the radial action J_r (to be defined below) is conserved in an average sense even during rapid potential changes.

This average conservation does not appear to have been derived elsewhere in the literature, although equivalent effects have been noted in numerical work by Dalal, Lithwick & Kuhlen (2010, see also Lithwick & Dalal 2011). We will first show how it can be derived from equations in Pontzen & Governato (2012) when the potential changes instantaneously, maintaining the sphericity of the halo. Then a more general (but more abstract) argument will be given which additionally shows that the other two actions (the z-component of the angular momentum j_z and the scalar angular momentum j) are also conserved in the same average sense. The second approach encompasses perturbations to the potential which have variations on arbitrary time-scales and may break spherical symmetry. However the first has a more intuitive content and therefore forms our starting point.

In a spherical system, the radial action J_r is defined by

\begin{equation} J_r = \frac{1}{\pi } \int _{r_{\mathrm{min}}}^{r_{\mathrm{max}}} \sqrt{2E-2\Phi (r; t) - j^2/r^2}\, \mathrm{d}r. \end{equation}

(2)

Here E is the specific energy, j is the specific angular momentum and Φ is the potential at a given radius r and time t; the r integral is taken over the region where the integrand is real.

The radial action J_r has the same units as specific angular momentum j. This reflects the similar conservation roles these two quantities play for the radial and angular components of the motion. In particular J_r is exactly conserved if any changes in the potential occur sufficiently slowly (‘adiabatically’) in time (e.g. Binney & Tremaine 1987).

On the other hand in the rapid, impulsive limit under an instantaneous change in energy E → E + ΔE and potential Φ(r) → Φ(r) + ΔΦ(r), the action of the particle is changed at first order:

\begin{equation} \Delta J_r = \left.\frac{\mathrm{\partial} J_r}{\mathrm{\partial} E}\right|_{\Phi } \Delta E + \int _0^{\infty } \mathrm{d}r\, \Delta \Phi \left. \frac{\delta J_r}{\delta \Phi } \right|_{E}. \end{equation}

(3)

In Pontzen & Governato (2012) we showed that the energy shift ΔE induced by the change of potential, averaged over possible orbital phases of the particle, is

\begin{equation} \langle \Delta E \rangle = - \frac{\int _0^{\infty } \mathrm{d}r \,\Delta \Phi \left. \delta J_r / \delta \Phi \right|_{E}}{\left. \mathrm{\partial} J_r / \mathrm{\partial} E \right|_{\Phi }}, \end{equation}

(4)

an exact result (see equation 12 of Pontzen & Governato 2012). Here angular brackets denote averaging over all possible phases of the orbit. Considering the probability distribution of radial actions after this change, one has 〈ΔJ_r〉 = 0 at first order, by substituting equation (4) in equation (3). Even though a specific particle will change its radial action, the ensemble average is conserved.

This result connects closely with the standard adiabatic argument that ΔJ_r = 0 if any changes to the potential occur on long time-scales. In our case, however, the necessary ‘phase averaging’ does not occur over time for an individual particle but instead via a statistical consideration of an ensemble of particles spread evenly through all possible phases.

We can generalize as follows. Adopting the complete set of action-angle coordinates for phase space (e.g. Binney & Tremaine 1987; Park 1990), the momenta are |$\boldsymbol {J}=(J_r, j, j_z)$|⁠, where j is the total angular momentum and j_z is its component in the z direction (so −j < j_z < j). The conjugate coordinates are |$\boldsymbol {\Theta }=(\psi _r, \phi , \chi )$|⁠, taken to be periodic with interval 2π, and the Hamiltonian is

\begin{equation} H(\boldsymbol {J}, \boldsymbol {\Theta }) = H_0(\boldsymbol {J}, t) + h(\boldsymbol {J}, \boldsymbol {\Theta }, t). \end{equation}

(5)

Here h is an arbitrary perturbation. It may consist of a long-lived term (perhaps a departure from spherical symmetry) and fluctuations on arbitrary time-scales. In the background, the |$\boldsymbol {\Theta }$| coordinates change at a constant rate,

\begin{equation} \boldsymbol {\Theta } = \boldsymbol {\Theta }_0 + \boldsymbol {\Omega } \, t, \end{equation}

(6)

and the frequencies for the radial and azimuthal motion Ω_r and Ω_j obey

\begin{equation} \Omega _r = \frac{\mathrm{\partial} {H_0}}{\mathrm{\partial} J_r}, \Omega _j = \frac{\mathrm{\partial} H_0}{\mathrm{\partial} j}, \end{equation}

(7)

by Hamilton’s equations. The frequencies may change slowly with time (dΩ/dt < Ω²). |$\boldsymbol {J}$| is conserved in the background but when h is non-zero, the equations of motion read

\begin{equation} \frac{\mathrm{d}\boldsymbol {J}}{\mathrm{d}t} = - \frac{\mathrm{\partial} {h}}{\mathrm{\partial} \boldsymbol {\Theta }}. \end{equation}

(8)

Because equation (6) shows that particles in the background move at a uniform rate in the |$\boldsymbol {\Theta }$| coordinates, an equilibrium distribution f has no |$\boldsymbol {\Theta }$| dependence, i.e. the density f is a function of |$\boldsymbol {J}$| alone. From this it follows that

\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t} \langle \boldsymbol {J} \rangle = -\int \mathrm{d}^3 J \mathrm{d}^3 \Theta f(\boldsymbol {J}) \frac{\mathrm{\partial} h}{\mathrm{\partial} \boldsymbol {\Theta }} = 0, \end{equation}

(9)

where the result is obtained via integration by parts. This means an individual particle’s |$\boldsymbol {J}$| can ‘diffuse’ (Binney & Lacey 1988) over large distances σ in action space,

\begin{equation} \sigma \equiv \langle (\delta J)^2 \rangle ^{1/2} = \mathcal {O}(\epsilon ), \end{equation}

(10)

compared to the variation in the mean μ of the population,

\begin{equation} \mu \equiv \langle \delta J \rangle = \mathcal {O}(\epsilon ^2). \end{equation}

(11)

We can inspect the diffusion in a simulation by calculating the relevant actions at two time-steps. As an example, Fig. 1 shows a histogram of changes in the radial action of tightly bound particles in the forming ‘dwarf’ halo (see Section 3) between z = 3.1 and 1.4, a time interval of approximately 2.7 Gyr. We define ‘tightly bound’ by selecting particles interior to 10 kpc at the earlier time-step, and calculate the J_r values of these particles in both outputs according to the numerical recipes given later (Section 3.1).

Figure 1.

An illustration of the diffusion of particles in action space. Here particles in the inner 10 kpc of a simulation have been selected and their radial actions J_r numerically calculated at two time-steps separated by Δt = 2.7 Gyr. The change in the population mean action is small (μ = −12.1 kpc km s⁻¹) compared against the magnitude of the random diffusion (σ = 287 kpc km s⁻¹).

Open in new tab Download slide

As expected from the linear analysis, Fig. 1 shows that individual simulated particles change their actions more rapidly than the population mean. Quantitatively, the change in the population mean μ is −12.1 kpc km s⁻¹ whereas a typical particle has moved by σ = 287 kpc km s⁻¹ ≫ |μ| from its original J_r value. Note also that the mean J_r value for these particles in the final time-step is 〈J_r〉 = 187 kpc km s⁻¹ < σ, meaning particles really do cover significant distances in action space. Similar results apply for j and j_z, respectively, the change in the population mean for these quantities μ is −14.9 and 9.9 kpc km s⁻¹ whereas the dispersion σ is 140 and 161 kpc km s⁻¹.

This analysis confirms that |$\langle \boldsymbol {J} \rangle$| evolves slowly and, although it is not exactly conserved, it forms a constraint of motion that cannot be ignored on finite time-scales. It is not yet clear to us whether or how this conservation relates to the phenomenologically conserved quantities identified by Stiavelli & Bertin (1987, see also Trenti & Bertin 2005; Trenti, Bertin & van Albada 2005). Indeed a complete description would require investigation of different moments of the distribution at higher order in perturbation theory. An exact link between phase mixing and invariance of mean action is suggested by the Poincaré invariant |$\oint \boldsymbol {J} \cdot \mathrm{d}\boldsymbol {\Theta }$|⁠, for instance, but the chaotic deformation of the contour makes it hard to formulate a precise statement from this starting point. For now the linear analysis and simulation statistics motivate a picture in which |$\langle \boldsymbol {J} \rangle$| evolves sufficiently slowly compared to the diffusion of an individual particle that it must be considered fixed in analysis of the distribution.

2.2 An aside on angular momentum

The third component of equation (9) seems to imply conservation of one component of the total angular momentum vector (proportional to 〈j_z〉). But the expected conservation of the other two components, proportional to 〈j_x, j_y〉, does not appear. In fact the normal angular momentum conservation laws arise only when the potential is generated self-consistently (through ∇²Φ = 4πGρ, where ρ is the real-space density corresponding to f), whereas the approximate conservation discussed in the section above requires only the weaker condition that interactions be collisionless. Nonetheless, it seems desirable to encode conservation of all components of the vector angular momentum in our final solution.

In the action-angle coordinate system j_x and j_y are indirectly represented by j, j_z and χ (the coordinate conjugate to j_z):

\begin{eqnarray} \left(\begin{array}{l} j_x \\ j_y \end{array}\right) = \sqrt{j^2 - j_z^2}\left(\begin{array}{l}\cos \chi \\ \sin \chi \end{array}\right). \end{eqnarray}

(12)

(Note that for particles orbiting in the unperturbed system H₀, χ is a constant of motion as well as the three actions since ∂H₀/∂j_z = 0.) Conservation of the x and y components of vector angular momentum is then expressed by a non-linear equation involving the j, j_z and χ distribution. For our current purposes we align the total angular momentum along the z-axis of our coordinate system; the requirement that 〈j_x〉 = 〈j_y〉 = 0 is satisfied if the particles are assumed phase mixed in the χ direction since then 〈sin χ〉 = 〈cos χ〉 = 0.

2.3 The new canonical ensemble

In the Introduction we explained that, to obtain a phase-space distribution function, we will maximize the entropy (1) subject to constraints on the particle population (further discussion is given in Appendix A). As well as energy conservation, we apply the three-vector of constraints on |$\langle \boldsymbol {J} \rangle$| discussed above. This gives rise to a total of four Lagrange multipliers in the resulting distribution function:

\begin{equation} f(\boldsymbol {J}) \propto \exp \left(-\boldsymbol {\beta } \cdot \boldsymbol {J} -\beta _E E(\boldsymbol {J})\right), \end{equation}

(13)

the Lagrange multipliers are |$\boldsymbol {\beta } = (\beta _r, \beta _j, \beta _z)$| and β_E; in the above expression |$E(\boldsymbol {J})=H_0(\boldsymbol {J})$| is the specific energy associated with a given orbit. In the absence of the new constraints, |$\boldsymbol {\beta }=0$| and β_E is identified with m/kT (where m is the particle mass and T is the thermodynamical temperature; note that the mass of the particle will not appear again and enters here only because of the conventional definition of temperature, and is not actually a consideration in the equilibrium achieved by a collisionless system).

All four constants can be determined in a variety of ways depending on the situation; for a complete account of structure formation one would like to be able to derive them from the initial conditions, but this lies beyond the scope of the current paper (although see Section 3.4 for further comments). The lack of any reference to |$\boldsymbol {\Theta }$| in equation (13) indicates that the solution is phase mixed, consistent with the initial assumptions.

Equation (13) is the essential prediction of the present work. As with any prediction derived from a maximum entropy argument, it will be able to fit the actual ensemble only if we have encapsulated enough of the dynamics within the constraints (Jaynes 1979b). The rest of this paper is concerned with testing to what extent that is the case.

First consider the probability of finding a particle with J_r in a given interval (ignoring j and j_z coordinates). This is given by

\begin{equation} p_r(J_r) = \int _0^{2\pi } \mathrm{d}^3 \Theta \int _0^{\infty } \mathrm{d}j \int _{-j}^j \mathrm{d}j_z\, f(\boldsymbol {J}), \end{equation}

(14)

because the action-angle coordinates are canonical, so the phase-space measure is constant. In the limit that the energy of the system becomes large at fixed action, equation (14) can be solved:

\begin{equation} p_r(J_r) \propto \exp - \beta _r J_r (\beta _E=0), \end{equation}

(15)

but it is not immediately clear whether we will be operating in this regime. More generally a closed form for p_r(J_r) is hard to obtain, but we can at least show that

\begin{equation} \frac{\mathrm{d}\ln p_r}{\mathrm{d}J_r} = -\beta _E \langle \Omega _r \rangle _{J_r} - \beta _r, \end{equation}

(16)

which can then be integrated numerically for a given case to give a concrete comparison between equation (13) and simulations. Here we have defined

\begin{equation} \langle \Omega _r \rangle _{J_r} \equiv \frac{1}{p_r(J_r)}\int f(\boldsymbol {J}) \Omega _r(\boldsymbol {J})\, \mathrm{d}j\, \mathrm{d}j_z, \end{equation}

(17)

which is the mean of the radial frequency for particles with a fixed J_r. With this definition, relation (16) can be derived from equation (14), recalling that the radial frequency Ω_r of the particle’s orbit obeys equation (7). We will investigate and explain the distribution of J_r values predicted by equation (16) in Section 3.2.

Now consider the distribution of total angular momentum j. We will follow exactly the same series of manipulations as for the radial action; however, in this case the marginalization over j_z introduces a non-trivial term:

\begin{eqnarray} p_j(j) &=& \int _0^{2 \pi } \mathrm{d}^3 \Theta \int _{-j}^{j} \mathrm{d}j_z\, \int _0^{\infty } \mathrm{d}J_r\, f(\boldsymbol {J}) \nonumber \\ && \propto \; \sinh \left( \beta _z j \right) \, \exp \left(-\beta _j j \right) \int _0^{\infty } \mathrm{d}J_r \exp \left( - \beta _E E \right). \end{eqnarray}

(18)

Once again, in the case β_E → 0, we have a fully analytic expression for p_j(j),

\begin{equation} p_j(j) \propto \sinh \left(\beta _z j \right) \exp \left( - \beta _j j \right) (\beta _E=0), \end{equation}

(19)

which will serve as a useful point of comparison. More generally we can differentiate equation (18) to obtain

\begin{equation} \frac{\mathrm{d}\ln p_j}{\mathrm{d}j} = \beta _z \coth \left(\beta _z j\right) -\beta _E \langle \Omega _j \rangle _j - \beta _j, \end{equation}

(20)

where Ω_j is the angular frequency of the orbit and

\begin{equation} \langle \Omega _j \rangle _j \equiv \frac{1}{p_j(j)} \int f(\boldsymbol {J}) \Omega _j(\boldsymbol {J})\, \mathrm{d}J_r \, \mathrm{d}j_z \end{equation}

(21)

is the mean angular frequency of particles at fixed j. Equation (20) for the angular momentum distribution (ignoring all other coordinates) is the equivalent of equation (16) for the radial action distribution. Once again we will investigate and explain the shape it predicts in Section 3.3. First, however, we will explain the simulations which serve as a point of comparison for the later discussions.

3 COMPARISON TO SIMULATIONS

3.1 Overview of the simulations

In the previous section we applied maximum entropy reasoning to conservation of energy and approximate conservation of action to derive an expected equilibrium phase-space distribution. We will now compare that expectation against simulated dark matter haloes. Our strategy is to integrate equations (16) and (20) numerically for these simulations and compare to the actual distribution of particles binned by J_r and j, respectively.

We will present results from three simulated dark matter haloes (shown in Fig. 2), chosen to span a wide range of masses with an approximately constant number of particles per halo (several million in each case). We also compared our results against the GHALO multibillion particle phase space (Stadel et al. 2009), finding good agreement similar to that described for our Milky Way (MW) halo here. This gives confidence that the mechanisms and results discussed in the paper are not sensitive to numerical resolution.

Images of the three dark matter simulations, accompanied by numerical properties. Respectively r200, M200, c, Npart and ϵ denote the radius at which the density exceeds the critical density by a factor 200, the mass within this radius, the ‘concentration’ c = r200/rs where rs is the NFW scale radius as described in the text, the number of particles within r200 and the gravitational softening length in physical units at z = 0. The images are scaled to show the virial sphere of the main halo. The brightness represents the column density of dark matter (scaled logarithmically to give a dynamic range of 3000 in each case); the colour corresponds to a density-weighted potential along the line of sight.

Figure 2.

Images of the three dark matter simulations, accompanied by numerical properties. Respectively r₂₀₀, M₂₀₀, c, N_part and ϵ denote the radius at which the density exceeds the critical density by a factor 200, the mass within this radius, the ‘concentration’ c = r₂₀₀/r_s where r_s is the NFW scale radius as described in the text, the number of particles within r₂₀₀ and the gravitational softening length in physical units at z = 0. The images are scaled to show the virial sphere of the main halo. The brightness represents the column density of dark matter (scaled logarithmically to give a dynamic range of 3000 in each case); the colour corresponds to a density-weighted potential along the line of sight.

Open in new tab Download slide

Our simulations are run from cosmological initial conditions at z ≃ 100 in a ‘zoom’ configuration (Navarro & White 1993), i.e. with high resolution for the main halo and its immediate surroundings and lower resolution for the cosmological environment. The softening lengths ϵ for the high-resolution region are listed in Fig. 2 and are fixed in physical units from z = 9, prior to which they scale linearly with cosmological scalefactor, a compromise motivated by numerical convergence studies (Diemand et al. 2004). We verified at the final output (z = 0) that the high-resolution regions have not been contaminated by low-resolution particles, and that the halo real-space density profiles are well described by a slowly rolling power law, in accordance with all recent simulations (e.g. Diemand et al. 2008; Stadel et al. 2009; Navarro et al. 2010, and references therein).

Each simulation output contains full Cartesian phase-space coordinates (⁠|$\boldsymbol {x}$|⁠, |$\boldsymbol {v}$|⁠). The position space is re-centred on the central density peak of the halo using the ‘shrinking sphere’ method of Power et al. (2003). The velocities are re-centred such that a central sphere of radius r₂₀₀/30 has zero net velocity, where r₂₀₀ is the radius at which the mean halo density is 200 times the critical density. Henceforth we only consider particles inside r₂₀₀.

From left to right in Fig. 2 the simulated haloes become more massive. The width of each panel is equal to 2r₂₀₀ and the luminosity is scaled to represent the column density over a dynamic range of 3000. The most conspicuous aspect of Fig. 2 is that the haloes become less centrally concentrated. We verified this by fitting a classic Navarro–Frenk–White (NFW; Navarro et al. 1996b) formula to the density profile. The NFW fit,

\begin{equation} \rho (r) = \frac{\rho _0}{(1+r/r_{\rm s})^2\,(r/r_{\rm s})}, \end{equation}

(22)

yields ρ₀, a characteristic density, and r_s, a scale radius. The latter is often expressed in a scale-free manner as a concentration value c = r₂₀₀/r_s; we have recorded the value for each halo in Fig. 2. As expected the concentration decreases with increasing mass, in agreement with previously known trends (e.g. Bullock et al. 2001b; Macciò et al. 2007). We thus have a sample of cosmological haloes which span a wide range in both mass and concentration. These different concentrations are thought to arise from different mean densities in the Universe at the epoch of collapse (Navarro et al. 1997; Bullock et al. 2001b).

All haloes in Fig. 2 exhibit large amounts of substructure; we will present results with this substructure subtracted, although we have verified that including the substructure does not have a qualitative impact on our results. The substructure is identified and removed using the ‘Amiga Halo Finder’ (Knollmann & Knebe 2009). For each remaining particle inside r₂₀₀, the specific scalar angular momentum is given by |$j=\left|\boldsymbol {v} \times \boldsymbol {x}\right|$|⁠. We calculate J_r by evaluating equation (2) numerically, using a spherically averaged potential Φ defined by

\begin{equation} \Phi (r) = \int _0^r \mathrm{d}r^{\prime } \frac{GM({<}r^{\prime })}{r^{\prime 2}}, \end{equation}

(23)

where M(<r^′) is the total mass enclosed by a sphere of radius r^′, and the specific energy E of each particle is defined as |$E=\boldsymbol {v}^2/2 + \Phi (r)$|⁠. J_r is evaluated using the true spherical potential out to r_term = 3r₂₀₀, beyond which (for reasons of numerical speed) the calculation is truncated and an analytic completion assuming a Keplerian (vacuum) potential is taken. We verified that changing r_term to 4r₂₀₀ had little impact on the results.

Before proceeding to a comparison, we need to derive appropriate β values. We calculate these using a Monte Carlo Markov chain (MCMC) to maximize the likelihood:

\begin{equation} \mathcal {L}(\boldsymbol {\beta }, \beta _E) = \prod _i f(\boldsymbol {J}_i; \boldsymbol {\beta }, \beta _E), \end{equation}

(24)

where f is the one-particle distribution function (13) normalized such that |$\int \mathrm{d}^3 \boldsymbol {J}\, \mathrm{d}^3 \boldsymbol {\Theta } \,f(\boldsymbol {J}) = 1$|⁠. This normalization must be accomplished numerically on a grid of J_r, j values; at each grid point E(J_r, j) is calculated by operating a bisection search on equation (2). This need only be done once, and then the evaluation of each link in the Markov chain is rapid.

The operation gives us maximum likelihood (i.e. ‘best-fitting’) parameters¹ (β_j, β_z, β_r, β_E) for a given simulation, optionally subject to constraints (such as β_E = 0 or β_j = β_z = β_r = 0). We are now fitting up to four parameters (excluding mass normalization), more than the one or two parameters normally used by simulators to describe their haloes (e.g. Navarro et al. 2004; Stadel et al. 2009). However the fitted real-space density profiles are purely phenomenological constructs; conversely here we are starting with a functional form derivable from physical considerations. As we have commented in Section 2.3 and will expand upon in Section 3.4, the βs should ultimately therefore be derived from initial conditions. For the present, however, the objective is to see whether our physical argument can correctly describe the phase-space distribution at all, for which fitting βs is the most pragmatic approach.

3.2 Comparison with MW: |$\boldsymbol {J_r}$| distribution

We will now start to test how closely equation (13) represents the distribution of particles in our simulations. We will investigate MW in some detail, before showing results for the other two simulations to which essentially the same discussion can be applied.

The distribution of J_r values in MW is shown by the histogram in the left-hand panel of Fig. 3. This can be compared with the thick solid curve which shows the distribution of J_r values according to expression (16); the parameters are β^{− 1}_E = 1.6 × 10⁴ km² s^{− 2} and β^{− 1}_r = 4.4 × 10³ kpc km s^{− 1}. The agreement is excellent over several orders of magnitude in probability density, spanning the values 0 < J_r < J_break, where J_break ≃ 3 × 10⁴ kpc km s⁻¹. Particles with J_r > J_break account for less than 0.1 per cent of the mass and are on long-period orbits, probably reflecting new material falling into the potential well. We will not consider them further.

Figure 3.

The distribution of particles’ radial action J_r (left-hand panel) and scalar angular momentum j (right-hand panel) in MW (both panels show the distribution of particles as a histogram). Also shown are maximum entropy solutions based on energy conservation alone (dotted curve), action conservation alone (dashed curve) and our advocated solution using both constraints (thick solid curve). The last of these provides a good reproduction of the distribution for J_r < J_{r, break} = 3 × 10⁴ kpc km s⁻¹ and j < j_break = 4 × 10⁴ kpc km s⁻¹, respectively, while extending over orders of magnitude in probability density. Less than 0.1 per cent of particles lie at J_r > J_{r, break}, approximately 0.2 per cent lie at j > j_break. Despite the good overall agreement, problems become apparent at very small j; a blow-up of the indicated range j < 3500 kpc km s⁻¹ is given in Fig. 4.

Open in new tab Download slide

The shape of the J_r solution can be understood as follows. We have already remarked that, in the limit β_E → 0, one recovers equation (15), an exact exponential (i.e. a straight line on the linear-log axes of Fig. 3). For comparison we have plotted the best-fitting distribution of this form (with β^{− 1}_r = 3.5 × 10³ kpc km s^{− 1}) as a dashed line. Since the period of an orbit increases with its energy (or radial action), the mean frequency |$\langle \Omega _r \rangle _{J_r}$| decreases for increasing J_r. So, inspecting equation (16), there will always be a J_r value above which |$\beta _E \langle \Omega _r \rangle _{J_r}$| becomes much smaller than β_r. Looking again at the thick solid curve in Fig. 3, the limiting solution at high J_r is indeed a pure exponential as this reasoning would suggest. At small J_r, however, the gradient of the solution is steeper because of the energy term.

Comparing the histogram, the thick solid curve and the dashed line in the left-hand panel of Fig. 3 thus leads us to the conclusion that J_r conservation (dashed line) accounts rather well for the qualitative form of the distribution, with an important correction from E conservation at low J_r. Finally, the dotted curve shows the best-fitting case with β_r = 0 – i.e. the normal statistical mechanical result in the absence of other constraints – and provides a poor fit at all J_r. In summary, the identification of the J_r constraint has resulted in dramatic improvements in the match to simulations.

3.3 Comparison with MW: |$\boldsymbol {j}$| distribution

Now consider the right-hand panel of Fig. 3 which shows the distribution of scalar angular momentum for the particles in our MW simulation. Once again the simulated particles are shown by the histogram; the best-fitting maximum entropy solution (β^{− 1}_j, β_z^{− 1} = 4.4 × 10³, 1.1 × 10⁴ kpc km s^{− 1}, with β_E as quoted above) is shown by the thick solid curve. It again reproduces the correct qualitative behaviour up to j_break = 4 × 10⁴ kpc km s⁻¹, with only 0.2 per cent of the mass at j > j_break. Although the angular momentum distribution has some fluctuations away from the predicted behaviour, the predictions remain nearly correct over two orders of magnitude in probability density. With the exception of a problem described below, we do not believe these fluctuations to be of particular importance beyond indicating the structure is not completely relaxed. In particular we will show later (Section 3.5) that these inhomogeneities can be ignored when reconstructing a density profile in real space. Certainly compared against a solution based on E conservation alone, again shown by a dotted curve, our solution can be counted a success.

The basic shape of the predicted j distribution can be understood in a similar way to the J_r distribution explained above. Consider again the case where β_E = 0 (so in effect the total energy is unconstrained), then the exact solution is given by equation (19). We can also take the isotropic limit, β_z → 0, giving

\begin{equation} p_j(j) \propto j\, \exp (-\beta _j j) (\beta _z=0,\textrm { }\beta _E=0). \end{equation}

(25)

This is analogous to the radial action case (15), but with a degeneracy factor j reflecting the increasing density of available states available as the angular momentum vector grows in size. The result is that the abundance of particles grows linearly with j for j < β^{− 1}_j and decays exponentially for j > β^{− 1}_j.

In light of the above discussion, it is notable that the turnover from growth to decay in p_j(j) occurs at j values much smaller than β^{− 1}_j. There are two ways to accomplish this. The first is to create a highly anisotropic set-up, β^{− 1}_z ≪ j₀, where j₀ is the smallest j value of interest. This packs orbits as much as possible into a single plane, generating a large net angular momentum and destroying the approximate spherical symmetry,² but effectively removing the degeneracy in j altogether:

\begin{equation} p_j(j) \propto \exp \left[(\beta _z-\beta _j) j\right] (\beta _z^{-1}\ll j, \textrm { }\beta _E=0). \end{equation}

(26)

Because it is maximally anisotropic, this solution cannot reflect the simulations; however, if we temporarily fit only j values using the functional form (19), we are pushed towards this unphysical limit (dashed line, Fig. 3, right-hand panel; (β^{− 1}_j, β_z^{− 1}) = (5.9, 6.3) × 10² kpc km s^{− 1}).

Luckily this is not the only way to overcome the shrinking phase space at low j. Equation (20) shows that if 〈Ω_j〉_j increases fast enough as j → 0 it can overcome the |$\coth (\beta _z j)$| degeneracy term. We have verified that in the full solution (thick solid line in Fig. 3, right-hand panel), this is the mechanism by which the turnover is pushed to low j.

$Despite good agreement over the majority of j space (see right-hand panel of Fig. 3), the fraction of simulated orbits (histogram) at very low angular momentum is substantially underestimated by the simplest maximum entropy argument (thin dotted curve). One fix discussed in the text is to postulate a second population at low energies (dashed curve). This yields a much better low-j fit (solid line) without affecting the high-j fit (except through a minor renormalization).$

Figure 4.

Despite good agreement over the majority of j space (see right-hand panel of Fig. 3), the fraction of simulated orbits (histogram) at very low angular momentum is substantially underestimated by the simplest maximum entropy argument (thin dotted curve). One fix discussed in the text is to postulate a second population at low energies (dashed curve). This yields a much better low-j fit (solid line) without affecting the high-j fit (except through a minor renormalization).

Open in new tab Download slide

Focusing attention on the low-j part of the distribution does, however, reveal a deficiency in our predictions. Fig. 4 shows the distribution of orbits with j < 3500 kpc km s⁻¹. The dotted line shows the same maximum entropy fit depicted by the solid line in Fig. 3. When the horizontal scale is expanded in this way, it becomes clear that the global fit undershoots the simulation values significantly at low j. This appears to be a systematic feature of all simulations we have inspected (the three detailed here, GHALO and various other lower resolution simulations which we used for testing purposes). It is possible to force a better fit by restricting the likelihood analysis to this region, but the global agreement is then considerably worse.

This suggests that the behaviour at low j is marginally decoupled from that in the rest of phase space. This could arise from the wide range of orbital periods: particles with small actions also have periods much shorter than the rest of the halo. The coupling between particles will necessarily be weak if their time-scales are very different (since particles on short orbits react adiabatically to fluctuations on long time-scales). This can substantially suppress redistribution of scalar angular momentum and is consistent with, although not reliant on, the early formation of a stable central cusp in simulations (e.g. Moore et al. 1998; Lu et al. 2006; Wang et al. 2011).

In principle this weakness of coupling between orbits in different regions could be expressed as a further constraint in the maximum entropy formalism. Further investigation awaits future work, but for now we will use this as a motivation to study a two-population system. We are not suggesting that there really are two sharply defined populations, but that this should anticipate the features of incomplete equilibrium.

Our maximum likelihood analysis is able to find a dramatically better fit in this case, placing 3.5 per cent of the mass in a second population at substantially lower temperature (β^{− 1}_E = 2.8 × 10³ km² s^{− 2}). The summed distribution is shown by the thick solid line in Fig. 4, with the contribution from the subdominant population indicated by the dashed line. Because the second distribution is so peaked near j = 0, the only difference at high j is a marginal renormalization. We also verified that the J_r distribution is barely affected.

It is undeniably disappointing that our solution does not automatically accommodate the behaviour at very low j, but we expect that future development of the ideas above can quantitatively account for the discrepancy. We consider other possible explanations in Section 4. However after focusing so much on one corner of phase space we should re-emphasize the major conclusion: the distribution over both J_r and j for 96 per cent of the particles is remarkably well described by the maximum entropy expression (13).

3.4 Other simulations

We confirmed that these basic conclusions persist in other simulations. The top row of Fig. 5 shows results from the ‘dwarf’ simulation. The J_r distribution (top left-hand panel) again shows excellent agreement for J_r < J_break, where J_break = 4 × 10³ kpc km s⁻¹. Only a tiny fraction of mass (<0.1 per cent) lies beyond this point of breakdown. As with MW, the dwarf’s angular momentum distribution (top right-hand panel) has more conspicuous fluctuations, but still roughly adheres to the maximum entropy solution up to j_break = 3 × 10³ kpc km s⁻¹, with around 1.3 per cent of mass lying beyond this point. We verified that at very low angular momenta j < 100 kpc km s⁻¹ there is again an overabundance of particles in the simulation, although in this case it accounts for less than 2 per cent of particles compared against the 3.5 per cent in MW. The parameters of the dwarf fit are (β^{− 1}_r, β_j^{− 1}, β^{− 1}_z) = (4.2, 1.9, 2.1) × 10² kpc km s^{− 1} and β^{− 1}_E = 2.2 × 10³ km s^{− 1}.

$As Fig. 3, but for the remaining two simulations. Once again, the maximum entropy distribution subject to $\boldsymbol {J}$ and E constraints (thick solid lines) predicts the simulations (histogram) accurately. The simulated distribution in Jr for the cluster (lower left-hand panel) has some noticeable fluctuations over large scales; this is likely because it is dynamically young.$

Figure 5.

As Fig. 3, but for the remaining two simulations. Once again, the maximum entropy distribution subject to |$\boldsymbol {J}$| and E constraints (thick solid lines) predicts the simulations (histogram) accurately. The simulated distribution in J_r for the cluster (lower left-hand panel) has some noticeable fluctuations over large scales; this is likely because it is dynamically young.

Open in new tab Download slide

Considering the cluster simulation (lower row of Fig. 5) gives similar results once again. This time the J_r distribution as well as the j distribution shows some notable fluctuations around the maximum entropy description. This may be because clusters assemble later (as we discussed in Section 3.1, this is reflected in the lower concentration value), so the system is dynamically young; however, we have not explicitly looked at time dependence of these distributions. The parameters of the cluster fit are (β^{− 1}_r, β_j^{− 1}, β^{− 1}_z) = (1.5, 0.8, 1.2) × 10⁵ kpc km s^{− 1} and β^{− 1}_E = 3.5 × 10⁵ km s^{− 1} with 0.1 and 1.1 per cent of the mass in the unrelaxed components beyond J_break = 7 × 10⁵ kpc km s⁻¹ and j = 6 × 10⁵ kpc km s⁻¹, respectively.

Naively one would expect β^{− 1}_E to scale approximately as

\begin{equation} \beta _E^{-1} \propto \frac{GM_{200}}{r_{200}} \propto M_{200}^{2/3}, \end{equation}

(27)

since M₂₀₀ and r₂₀₀ are by definition related through the fixed-mean-density condition M₂₀₀∝r³₂₀₀. Similarly the actions |$\boldsymbol {\beta }^{-1}$| should scale as

\begin{equation} \boldsymbol {\beta }^{-1} \propto r_{200} \sqrt{\frac{GM_{200}}{r_{200}}} \propto M_{200}^{2/3}. \end{equation}

(28)

Fig. 6 compares these expectations with the actual values, although we immediately caution against taking the scaling of three haloes too seriously. The upper panel shows the radial action, β^{− 1}_r, as a function of mass (crosses) with dotted lines indicating the scaling (28). The lower panel shows the same for the energy scales. Both panels show good agreement with the expected trends.

Figure 6.

The best-fitting scales for radial action (upper panel) and energy (lower panel) as a function of mass. The crosses show the values from the three simulations, while dotted lines give the expected scalings (27) and (28), which agree well with the simulations.

Open in new tab Download slide

For clarity we did not overplot the β_j values in Fig. 6, but these can be seen to be comparable to β_r. Because cosmological haloes are formed from near-cold collapse, their initial angular momentum will be small. The final dispersion of angular momentum is likely generated through a weak form of the radial orbit instability (e.g. Saha 1991; MacMillan, Widrow & Henriksen 2006; Bellovary et al. 2008; Barnes, Lanzel & Williams 2009). Thus the scales of the angular momentum distribution and the radial action distribution are likely to be intimately linked. This is one example of a dynamical consideration which should ultimately be used to link |$\boldsymbol {\beta }$| values to the initial conditions.

3.5 Real-space radial density profiles

We have shown that a first principles maximum entropy argument is capable of describing the phase-space distribution of particles in dark matter haloes, up to a small correction at low angular momentum. The natural next step is to ask what kind of real-space radial density profiles are implied by this phase-space distribution and whether these match the classic rolling-power-law shape given by simulations.

Calculating the density distribution corresponding to the phase-space distribution (13) is technically involved; a description is given in Appendix B. There we also explain how the same computer code can be used to calculate equilibrium density profiles from simulations (as opposed to analytic distributions). These profiles are generated subject to our simplifying assumptions of phase mixing and spherical symmetry. They agree well with traditional ‘binned’ estimates of the density, validating the assumptions. Furthermore in the new method, each simulated particle is smeared out over its orbit, resulting in considerably smaller Poisson noise than from traditional binned estimates.

Applying the algorithm to the maximum entropy solution, we find that the radial density profile implied by equation (13) follows a shallow power law in the centre and steepens with increasing radius, in qualitative agreement with the behaviour seen in numerical simulations (Navarro et al. 1997, 2004; Stadel et al. 2009, and references therein). However, when using the single population phase-space distribution fits, the central slope is too shallow (dotted line, Fig. 7). One can obtain higher central densities and inner slopes by changing the parameters, but then the outer slope becomes too steep. On the other hand, if one adopts the incomplete relaxation fit advocated in Section 3.3, a vastly improved real-space density profile is recovered (thick solid line, Fig. 7). This confirms that the ∼3.5 per cent population at low j is responsible for controlling the cusp. In the discussion below we will recap our current understanding of this issue and give directions for future investigation.

Figure 7.

The real-space density distribution (upper panel) of the one- and two-component maximum entropy solutions (dotted and solid lines, respectively) compared to the MW simulation binned density profiles (dots). The softening length in MW is 170 pc, so the profile should be reliable exterior to ∼700 pc (e.g. Power et al. 2003). The generic maximum entropy result is a density profile with slowly steepening power law to increasing radii, in agreement with the simulations. The one-component fit misses the central density cusp, showing that the ∼3.5 per cent correction to the low angular momentum orbits (Section 3.3) is required to reproduce this quintessential feature of simulated dark matter haloes. The lower panel shows the cumulative mass as a function of radius.

Open in new tab Download slide

4 DISCUSSION

We have shown that maximizing the entropy of a distribution function subject to constraints on total action and energy reproduces the phase-space density of particles in simulated dark matter haloes. Crucially, there is a clear physical motivation behind this choice of constraints. We started by explaining that, since any equilibrium distribution must be phase mixed, the late stages of relaxation approach this phase-mixed state. As a consequence |$\langle \boldsymbol {J} \rangle$| becomes a conserved quantity as equilibrium is approached (Section 2.1). This constitutes a dynamical barrier to continued evolution, preventing energy from being further redistributed.

The resulting canonical ensemble (i.e. the maximum entropy solution) is given by equation (13). From it we derived two key relationships which can be used to test the phase space of simulated haloes, respectively, equations (16) and (20). These were used to demonstrate a close agreement between simulations and theory (Sections 3.2 and 3.3) over orders of magnitude in probability density, and over a wide range of halo masses from dwarf galaxies to clusters (Section 3.4). We compared to the Lynden-Bell (1967) distribution which is obtained when energy can be arbitrarily redistributed between particles, finding that our new canonical ensemble offers a vastly improved fit (see dotted lines in Fig. 3). This strongly suggests that (a) maximum entropy with suitable dynamical constraints (representing incomplete violent relaxation) is a plausible route to understanding the 6D phase space of dark matter haloes; (b) the newly constrained quantities need not be conserved in general, but must be conserved whenever the system is close to equilibrium, so that their value becomes fixed as the dynamics settle down and (c) we have identified a physical argument leading to an important example of these constraints.

However, we found an overabundance of low angular momentum orbits in the simulations relative to the analytic predictions (Section 3.3). This implies that there is at least one more important constraint that we have not fully reflected in our analysis. Constructing the radial density profile (Section 3.5) confirms that, although a small fraction (≲4 per cent) of particles are causing the discrepancy, their existence is essential to understanding the origin of the central density cusps seen in numerical simulations.

The correspondingly large density of particles as j → 0 has been found by previous work, notably Bullock et al. (2001a), who offered a fitting formula which implies continually increasing particle numbers towards j = 0. With our higher resolution simulations, we can see that p_j(j) does eventually decrease at sufficiently low j, but slower than expected given the shrinking available phase space (Fig. 4). Furthermore, the large number of particles at low angular momentum can be linked directly to various discrepancies between Λ cold dark matter (ΛCDM) theory and observed galaxies (van den Bosch, Burkert & Swaters 2001; Dutton & van den Bosch 2009). In particular there must be mechanism to remove the low angular momentum baryons (Dekel & Silk 1986; Binney, Gerhard & Silk 2001; Governato et al. 2010; Brook et al. 2011). Understanding what causes the accumulation of low angular momentum material in the first place is now added to the list of puzzles in this area.

Our maximum entropy picture gives an interesting framework in which to interpret the situation. In Section 3.3 we gave an extensive analysis of equation (20) which suggests two routes to adding material at low j. The first option is to appeal to anisotropy (first term on the right-hand side); the second is to use a population of particles at low energy (high β_E in the second term on the right-hand side). We currently prefer the second explanation for the following reason. Particles near the centre have very short orbital periods, which make them decouple from fluctuations on the dynamical time-scale of the remainder of the halo. In numerical simulations, the cusps are indeed the first part of the halo to form, and they do not change much at late times (Moore et al. 1998; Syer & White 1998; Wang et al. 2011). Lu et al. (2006) construct an explicit two-phase model of the formation of haloes reflecting this differentiation, emphasizing the lack of equilibriation between the inner and outer parts of the halo (see also Lapi & Cavaliere 2011). Accordingly a time-scale constraint could be incorporated from the outset of the maximum entropy argument; we expect this would give similar results to our current approach of fitting a second population. This will be tackled explicitly in future work.

The alternative view is that the behaviour at low j may be sensitive to effects of asphericity. This could modify the effective degeneracy. But as we commented in Section 3.3, the only obvious method available is to pack orbits tightly into a plane, so making the phase space available uniform with j, rather than linearly increasing. Numerical results do show haloes become more anisotropic towards their centre (Jing & Suto 2002). On the other hand, when given a second population to fit, our code does not select this as a viable explanation for the existence of the cusp (Section 3.3).

If a full description of the physics generating low angular momentum orbits can be reached, the work in this paper lays the foundation for a complete description of the collisionless equilibria of dark matter haloes. Further questions of interest will include the following:

whether and how the constraint vector |$\boldsymbol {\beta }$| can be derived from initial conditions (which will likely lean heavily on an understanding of the radial orbit instability; e.g. Saha 1991; MacMillan et al. 2006; Bellovary et al. 2008; Barnes et al. 2009);
whether and how maximum entropy arguments can explain the power-law behaviour of the pseudo-phase-space density ρ(r)/σ³(r) (Taylor & Navarro 2001);
whether and how the phase mixing is maintained to sufficient accuracy to make the |$\langle \boldsymbol {J} \rangle$| conservation effective over periods of major disturbance (perhaps through chaotic mixing; e.g. Merritt & Valluri 1996; Henriksen & Widrow 1997);
how various moments of the distribution function evolve at higher order in perturbation theory (which is closely related to the previous question);
how the arguments are changed by adopting an explicitly aspherical background Hamiltonian (assuming this is not already answered when attacking the low-j question) and
how maximum entropy arguments can be adjoined to microphysical descriptions of cusp destruction (Navarro, Eke & Frenk 1996a; El-Zant, Shlosman & Hoffman 2001; Weinberg & Katz 2002; Read & Gilmore 2005; Mashchenko, Couchman & Wadsley 2006; Pontzen & Governato 2012) to shed further light on this essential area of galaxy formation.

Our substantial step forward should give confidence that a full statistical account of the distribution of particles in simulated dark matter haloes is achievable without any ad hoc assumptions or modifications to the well-established principle of maximum entropy. Such an account would be extremely powerful for practical and pedagogical aspects of understanding the behaviour of dark matter in the Universe.

AP gratefully acknowledges helpful conversations with Steven Gratton, James Binney, Justin Read, Simon White, Carlos Frenk, Julien Devriendt, James Wadsley, Phil Marshall, Julianne Dalcanton, Jorge Peñarrubia and John Magorrian, and thanks Kieran Finn for development of computer code for a related project. We thank the anonymous referee for helpful comments. The MW simulation was run by Alyson Brooks. The GHALO simulation was kindly made available by Joachim Stadel and Doug Potter. FG was funded by NSF grant AST-0908499, NSF grant AST-0607819 and NASA ATP NNX08AG84G. Simulations were run on NASA Advanced Supercomputing facilities. Simulation analysis was performed with the pynbody package (http://code.google.com/p/pynbody) on the DiRAC facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Oxford. This work was supported by the Oxford Martin School and the Beecroft Institute of Particle Astrophysics and Cosmology.

The MCMC technique also yields uncertainties on the β values, but these will not be considered further in the present work.

We note in passing that, technically, distribution functions with net angular momentum can nonetheless generate spherical potentials (Lynden-Bell 1960).

Lynden-Bell (1967) discusses an exclusion principle which arises from Liouville’s theorem and can modify the classical expression for entropy; the deviations will be significant if the initial phase-space density is comparable to the density in any regions of a final ‘coarse-grained’ view of phase space. This does not seem likely to apply in cosmological settings (Shu 1978), although see Barnes & Williams (2012) for a different view.

For a unique answer we also need to demand that entropy be invariant under coordinate transformations of the phase space, which reflects our initial ignorance of the distribution of particles before the dynamics are specified. For an alternative view, see Hjorth & Williams (2010) and Williams, Hjorth & Wojtak (2010) who suggested that the measure should be uniform in energy space. This is equivalent to imposing an a priori preference for some regions of phase space over others, one that should ultimately be derived from the equations of motion (and could then be re-expressed as a constraint).

REFERENCES

Barnes

E. I.

Williams

L. L. R.

. ,

ApJ

2012

, vol.

748

pg.

144

Month:	Total Views:
January 2017	1
February 2017	4
March 2017	3
April 2017	1
May 2017	1
June 2017	2
July 2017	2
August 2017	1
September 2017	6
October 2017	6
November 2017	5
December 2017	7
January 2018	11
February 2018	8
March 2018	15
April 2018	9
May 2018	8
June 2018	6
July 2018	14
August 2018	7
September 2018	8
October 2018	3
November 2018	6
December 2018	4
January 2019	6
February 2019	11
March 2019	15
April 2019	18
May 2019	17
June 2019	9
July 2019	12
August 2019	16
September 2019	6
October 2019	13
November 2019	18
December 2019	25
January 2020	9
February 2020	6
March 2020	15
April 2020	10

Article Contents

Abstract

1 INTRODUCTION

2 THE ANALYTIC PHASE-SPACE DISTRIBUTION

2.1 Conservation of action

2.2 An aside on angular momentum

2.3 The new canonical ensemble

3 COMPARISON TO SIMULATIONS

3.1 Overview of the simulations

3.2 Comparison with MW: |$\boldsymbol {J_r}$| distribution

3.3 Comparison with MW: |$\boldsymbol {j}$| distribution

3.4 Other simulations

3.5 Real-space radial density profiles

4 DISCUSSION

REFERENCES

APPENDIX A: WHY MAXIMIZE ENTROPY?

APPENDIX B: DYNAMICAL DENSITY ESTIMATES

APPENDIX C: CORRECTIONS AT APOCENTRE AND PERICENTRE

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only