STATISTICAL MECHANICS OF COLLISIONLESS ORBITS. I. ORIGIN OF CENTRAL CUSPS IN DARK-MATTER HALOS

The first major modification consists in noting that in an equilibrium collisionless system all particles retain their energies. This is not the case in an equilibrium collisional system, such as a classical gas, or a self-gravitating system dominated by two-body encounters, such as a globular cluster. The fundamental property of a collisionless system is that particles are distributed on orbits. In a relaxed system, once an energy is assigned to a particle it stays in a restricted portion of phase space, an energy cell. For an isotropic system, energy is the only isolating integral. Hence, we argue that in partitioning state space, the fundamental partition is energy space, not phase space (see also Efthymiopoulos et al. 2007).

This implies that the occupation number n, Equation (6), should be interpreted as an indicator of the number of particles with a given energy, i.e.,

$$\begin{equation} N(E) \propto \exp (-\beta E) \end{equation} \tag{ 8 }$$

and not as the number of particles in a parcel of classical phase space, as is usually assumed.

Binney (1982; see also, e.g., Hernquist 1990; Spergel & Hernquist 1992; Binney & Tremaine 1987) found that elliptical galaxies obeying the R^1/4 law have energy distributions very similar to Equation (8), however these authors' motivation for using this form was largely empirical.

Given N(E) one can, of course, recover the classical distribution function. Stiavelli & Bertin (1987) point out that the occupation numbers in the classical phase space can be related to those in the energy space if the former are assigned non-equal a priori weights, inversely proportional to the volume of phase space with a given energy, i.e.,

$$\begin{equation} f(E) \propto g(E)^{-1} \exp (-\beta E), \end{equation} \tag{ 9 }$$

from which Equation (8) is retrieved since N(E) ≡ f(E)g(E) for isotropic systems. The relation between the phase-space density f and the differential energy distribution N is obtained assuming isotropy from f(x, v)d³vd³x = N(E)dE and the density of states $g(E) = d^3{\bf v} d^3{\bf x}/{\it dE} = 16 \pi ^2 \int _0^{r_{\rm max}(E)} \sqrt{2(E-\Phi)} r^2 dr$ (Binney & Tremaine 1987).

Another interesting aspect of Equation (8), as pointed out by Binney (1982), is that for systems with less mass at the centers than in the outer parts, β, and hence the effective temperature, T = β⁻¹, need to be negative (see also Merritt et al. 1989). This is consistent with the self-gravitating nature of the systems, as these are characterized by negative heat capacity (Binney & Tremaine 1987).

2.2.1. Integer Approach

The second major modification is necessary because in systems with a finite potential depth the occupation numbers in energy cells with very bound energies can become small (Equation (8)). A similar problem is encountered when the classical phase-space distribution function, f(E), takes on an exponential form (Equation (7)), but in this case, because the temperature is positive, the small occupation numbers are found at the near escape energies, i.e., in the outer regions of systems. Madsen (1996) pointed out that in this case the Stirling approximation, Equation (4), breaks down. Following Simons (1994), he argued that the appropriate form for Equation (6) is

$$\begin{equation} n_i = [g_i \exp (-\alpha -\beta E_i)], \end{equation} \tag{ 10 }$$

where [·] means rounding down to the nearest integer.

Because the majority of particles in this latter case have energies near escape energies, the small occupation number modification has a dramatic effect on the resulting structures. Madsen (1996) showed that Equation (10) introduces a cutoff which results in finite-mass systems, similar to the King (1966) models of globular clusters. In effect, Madsen (1996) analytically derived the well-known King (1966) models, which were originally obtained as a simple heuristic modification of the isothermal sphere's distribution function, Equation (7).⁴

For collisionless systems, where the temperature is negative, the cutoff occurs at energies close to that of the central potential value, i.e., close to the center of the system. Hence, the effect on the structures is not so dramatic in terms of total mass. But as we show below, it determines the inner density profile of the equilibrium systems.

2.2.2. Continuous Approach

Because a physical system is not expected to have a step-like differential energy distribution N(E), a continuous version of Equation (10) is needed. We do this by introducing a superior approximation to Stirling's formula which, unlike Equation (5), is not limited to large occupation numbers.

We start by noting that when n = 0, one has the exact result ψ(1) = −γ, where γ ≈ 0.57721566... is Euler's constant. Combining this with the large n limit, Equation (5) leads to the approximation

$$\begin{equation} \psi (n+1)\approx \ln (n + \zeta), \end{equation} \tag{ 11 }$$

with ζ = exp(−γ) ≈ 0.561459... . While simple, this expression is a remarkably good approximation as shown in Figure 1 along with the classical large n approximation (Equation (5); ζ = 0). The deviation from the exact ψ(n + 1) is always positive, with a single maximum difference of 0.0237 at n = 0.680. Using this in Equation (2), one obtains a modification to Equation (6),

$$\begin{equation} n = g \exp (-\alpha -\beta E) - \zeta. \end{equation} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

Similar modifications have previously been proposed (Landsberg 1954; Hjorth 1993; Menon et al. 2006; Dubey et al. 2008).

The Lagrange multiplier α is determined by requiring that n = 0 at some energy E'. Thus, we get

$$\begin{equation} n = \zeta (\exp (-\beta [E - E^{\prime }])-1), \end{equation} \tag{ 13 }$$

thereby eliminating the cell size, g, which does not have a unique, physically meaningful value in gravitational systems.

For a system which vanishes at the escape energy, E' = 0, and so the classical phase-space occupation number function is n = ζ(exp(−βE) − 1). The systems we are interested in vanish at a finite central potential, E' = Φ₀, and so the energy-space occupation number function becomes n = ζ(exp(−β(E − Φ₀)) − 1). This is the continuous version of Equation (10).

Equation (13) with E' = Φ₀ incorporates the two modifications we introduced in Sections 2.1 and 2.2. Identifying the occupation number n as being proportional to N, the differential energy distribution, one obtains N(E) = A(exp(−β[E − Φ₀]) − 1), where A is determined by the mass of the system. In this expression, E and Φ₀ have units of energy, and β is the inverse of temperature. We convert this to dimensionless form, using ε = βE and ϕ₀ = βΦ₀ (note that both ε and ϕ₀ are positive quantities for bound systems),

$$\begin{equation} N(\varepsilon)=A(\exp [\phi _0-\varepsilon ]-1). \end{equation} \tag{ 14 }$$

We dub this expression the "DARKexp". It represents our prediction for fully relaxed, collisionless, self-gravitating, isotropic systems. Because the arguments presented in this paper do not apply to non-isotropic systems, Equation (14) cannot be directly compared to the results of N-body simulations. Having said that, we note that N(ε) is not very sensitive to anisotropy, as it depends primarily on the density profile and not on the dynamics of the system.

The detailed structure of the DARKexp systems, including their density and velocity dispersion profiles, will be considered in an accompanying paper (Williams & Hjorth 2010). In the next section, we discuss the limiting behavior of DARKexp models at small and large radii.

In a general spherically symmetric structure, the limiting form for the central potential can be assumed to be

$$\begin{equation} \phi =\phi _0-\phi _\alpha r^\alpha +\cdots \end{equation} \tag{ 15 }$$

Possion's equation ∇²ϕ = 4πGρ yields a limiting power-law behavior of the central density,

$$\begin{equation} \rho (r)\propto r^{\alpha -2}. \end{equation} \tag{ 16 }$$

For ε → ϕ₀, the distribution function then becomes

$$\begin{equation} f(\varepsilon)\propto (\phi _0-\varepsilon)^{-(4+\alpha)/2\alpha } \end{equation} \tag{ 17 }$$

and the density of states becomes

$$\begin{equation} g(\varepsilon)\propto (\phi _0-\varepsilon)^{(6+\alpha)/2\alpha }. \end{equation} \tag{ 18 }$$

For similar expressions, see Henriksen & Widrow (1995), Widrow (2000), and Arad et al. (2004). The differential mass distribution N(ε) = f(ε)g(ε) then becomes

$$\begin{equation} N(\varepsilon)\propto (\phi _0-\varepsilon)^{1/\alpha }. \end{equation} \tag{ 19 }$$

In general, the value of α ranges from 0 for the singular isothermal sphere to α = 1 for Navarro–Frenk–White or Hernquist profiles, to shallower slopes, reaching a flat core for α = 2. Increasing α corresponds to increasingly steeper ln N(ε) versus ε curves. For α → ∞, the system develops a hole in the central density profile, which is unphysical.

In the DARKexp model, Equation (14), N(ε) ∝ (ϕ₀ − ε) as ε → ϕ₀. In other words, α = 1 and we retrieve the central density cusps

$$\begin{equation} \rho (r) \propto r^{-1}, \end{equation} \tag{ 20 }$$

known from numerical simulations. The corresponding distribution function is

$$\begin{equation} f(\varepsilon) \propto (\phi _0-\varepsilon)^{-5/2} \end{equation} \tag{ 21 }$$

for ε → ϕ₀, which is the same as that of the Hernquist (1990) model.

Thus, our proposed differential energy distribution, the DARKexp, naturally predicts that the central density slopes of collisionless self-gravitating systems should asymptote to −1 at the centers of structures. Note that this slope is not the result of any specific dynamical process operating during the formation of halos, but a generic consequence of full relaxation.

The Appendix discusses the limiting behavior for more general cutoffs to N(ε) ∝ exp(−ε).

While we are focusing on the inner parts of halos in this paper, we note that a full model for the entire system can be obtained by assuming that N(ε) is finite at the escape energy and zero above (as plotted in Figure 2). The rationale behind this is that during violent relaxation the escape energy is no special location and N(ε) is expected to be continuous at what will eventually become the escape energy (Jaffe 1987; Hjorth & Madsen 1991). For such a model f(ε) ∝ ε^5/2 and ρ(r) ∝ r⁻⁴ for ε → 0 and r → ∞, similar to the Hernquist model and broadly consistent with numerical simulations of dark-matter halos. However, if relaxation is not complete at radii where particles have near escape velocities, then the outer density profile slope may deviate from −4.

Zoom In Zoom Out Reset image size

In this paper, we used statistical mechanics to predict the energy distribution function of relaxed systems. A key feature of statistical mechanics approaches is that they deal only with the final equilibrium states and derive these by equating them to the most probable or maximum entropy states. We stress that our theory is therefore limited to the description of the final state of self-gravitating collisionless systems. Because the final state is an equilibrium state, it is the result of full relaxation and therefore does not rely on any additional assumptions.

Implicit in our derivation is the assumption of equal a priori probabilities in state space. While this can be accomplished through efficient mixing (ergodicity), our theory does not address specific physical mechanism for attaining full mixing. In a companion paper (Williams & Hjorth 2010), we use the Extended Secondary Infall Model (ESIM) to test the DARKexp prediction, but stress that ESIM is a restricted physical model and is not equivalent to N-body simulations.

A straightforward way to evaluate the applicability of physical mechanisms for relaxation is to compare their end states. For example, the scattering model for violent relaxation introduced by Spergel & Hernquist (1992) predicts N(E) ∝ [(E − Φ₀)⁻² + C(E − Φ₀)^−1/2]⁻¹exp(−βE) which is clearly inconsistent with the DARKexp final state. On the other hand, the final state of ESIM halos is well fit with DARKexp.

Our prediction for N(ε) applies to the end states of collisionless dynamics. Collisionlessness has several different consequences, many of which can be used in one's theory. Lynden-Bell (1967) used the collisionless Boltzmann equation to introduce an exclusion principle to account for the incompressibility of the phase-space fluid. We use a different property: that collisionless dynamics implies the constancy of the energy per unit mass for each particle. Therefore, the collisionless nature of the problem is automatically included and there is no need for an exclusion principle.

In a system driven toward equilibrium by two-body relaxation, the situation is quite different. In this case, efficient relaxation implies that any particle can in principle end up anywhere in phase space and the usual partition of μ space is appropriate. In this case, taking into account low occupation numbers, one retrieves the King (1966) f(E) which is rounded down at near escape energies, while DARKexp's collisionless N(E) is rounded down at highly bound energies.

The density profiles of the two are clearly different, but both are very good approximations to what one sees in globular clusters and simulated dark-matter haloes, respectively.

In this paper, we dealt exclusively with isotropic systems for which the differential energy distribution is a function of energy only. The full simulated structures, however, are definitely anisotropic and exhibit a correlation between the density slope and anisotropy (Hansen & Moore 2006).

In principle, our statistical mechanics formalism can be extended to spherical non-isotropic systems, and we plan to do so in a future publication. In a fully developed theory, the distribution function must depend explicitly on angular momentum, especially in the outer parts (e.g., Stiavelli & Bertin 1987). Whether the energy distribution of the full theory will be significantly affected by the introduction of angular momentum is yet unclear; numerical experiments of MacMillan et al. (2006, Figure 9) seem to indicate that the DARKexp form is still a good description of systems that underwent radial orbit instability.

For now we note that even though the DARKexp prediction may not be expected to be an excellent fit to simulations at all energies, it should apply to the central regions of dark-matter halos, as these are isotropic. Therefore, our explanation of the central density cusp of −1 applies to simulated halos.