Power-law Pseudo-phase-space Density Profiles of Dark Matter Halos: A Fluke of Physics?

Dark matter, the dominant mass component of the universe, is the scaffolding that provides the structure for galaxies and clusters of galaxies. Therefore, understanding the structure of dark matter halos is one of the most important goals of modern cosmology. While the equilibrium structure of stars has been known for about a century, the structure of equilibrium dark matter halos has proven harder to establish. Using the equation of hydrostatic equilibrium and the equation of state for gaseous material allows one to solve for the internal structure of stars. In the case of collisionless dark matter halos, it is not possible to solve the hydrostatic equilibrium equation—also known as the Jeans equation—because the equation of state of dark matter is unknown.¹

In 2001, an interesting observation was published by Taylor & Navarro (2001). Using N-body simulations, the authors measured mass density, ρ(r), and velocity dispersion σ(r), profiles of equilibrium dark matter halos. The quantity ρ/σ³ turned out to be a power law in radius over about three decades, despite the fact that neither density nor the velocity dispersion are power laws, but in fact significantly deviate from a scale-free form. Since then, it has been confirmed by a number of studies that dark matter halos formed in cosmological N-body simulations (Peirani et al. 2006; Hoffman et al. 2007; Knollmann et al. 2008; Ma et al. 2009; Ludlow et al. 2010; Navarro et al. 2010; Gao et al. 2012; Butsky et al. 2016; Nolting et al. 2016; Drakos et al. 2017), halos computed through an iterative collisionless spherical collapse (Austin et al. 2005), and even galaxies and clusters formed in the real universe (Chae 2014; Munari et al. 2014) are well characterized by a power law,

$$\begin{eqnarray}&&\displaystyle \frac{\rho (r)}{{\sigma }^{3}(r)}\propto {r}^{-\alpha }.\end{eqnarray} \tag{ 1 }$$

Because this quantity has the dimensions of phase-space density, it has been nicknamed pseudo-phase-space density.

Several papers made attempts to shed light on its origin (Austin et al. 2005; Henriksen 2006; Ludlow et al. 2011; Alard 2013). In the meantime, others continued to address the more general question of how to understand the structure of dark matter halos that develops so robustly in simulations (He & Kang 2010; Kang & He 2011; Salvador-Solé et al. 2012; Pontzen & Governato 2013; Beraldo e Silva et al. 2019). Hjorth & Williams (2010) proposed a theoretical derivation for the differential energy distribution of self-gravitating relaxed collisionless matter. Based on the principles of statistical mechanics, they proposed the most likely steady-state configuration of these systems. Their result, DARKexp, forms a one shape parameter family, with ϕ₀ characterizing the dimensionless depth of the central gravitational potential. DARKexp gives very good fits to the density profiles (Hjorth et al. 2015; Nolting et al. 2016) and, more importantly, to the differential energy distributions of simulated dark matter halos (Williams et al. 2010; Nolting et al. 2016). It also fits well the density profiles of observed equilibrium galaxy clusters (Beraldo e Silva et al. 2013). It was shown in Williams & Hjorth (2010) that the ρ/σ³ profiles of the DARKexp family are close to, but not exactly power laws for many values of ϕ₀, suggesting that it may not be a universal feature of relaxed systems. Hints of nonuniversality of this quantity have also been noted in other papers (Ma et al. 2009; Del Popolo 2011, 2015).

Recently, Nadler et al. (2017) challenged the physical origin of the pseudo-phase-space density. They consider 1D self-similar fluid collapse, following closely an earlier study by Bertschinger (1985). The authors follow the evolution of gas entropy, whose definition is effectively the same as that of the pseudo-phase-space density, ρ/σ³. Because they are dealing with gas, their treatment cannot incorporate velocity anisotropy, which is measured to be nonzero in numerical dark matter simulations (Hansen & Moore 2006; Hansen & Stadel 2006; Lemze et al. 2012), as well as some observations of galaxies and galaxy clusters (Hansen & Piffaretti 2007; Hansen et al. 2011; Longobardi et al. 2018), and stellar and globular cluster populations in the Milky Way (Vasiliev 2019).

This paper is a further attempt to understand ρ/σ³: is there some physical principle behind its power-law nature, or is it a mere coincidence? Demonstrating the existence of an underlying physical principle will have important implications for our understanding of self-gravitating collisionless systems. Here, we do not address the possible physical meaning of ρ/σ³, but instead assume that if one exists, the final equilibrium state of halos will satisfy Equation (1).

Our approach differs significantly from that of Nadler et al. (2017). While these authors dealt with evolution of isotropic gaseous material, the present work does not consider halo evolution, and instead concentrates on the equilibrium state of halos, where velocity anisotropy plays an important role. Though our methods are very different, our conclusions are essentially the same as theirs: pseudo-phase-space density is unlikely to have a fundamental physical interpretation, and hence cannot help in the understanding of dark matter halos.

Because we are dealing with equilibrium dark matter halos, the starting point of our analysis is the Jeans equation, a statement of hydrostatic (or, mechanical) equilibrium for collisionless matter. Based on the results of N-body simulations, all the quantities characterizing the spherically averaged equilibrium halos, namely (i) the density profile, (ii) velocity anisotropy profile, and (iii) pseudo-phase-space density profile can be modeled, to a good approximation, as simple analytic relations, with a total of four or six parameters, depending on how the density profile is represented.

One can make even stronger statements regarding (i) and (iii), going beyond fitting functions to simulations.

The density profiles of relaxed systems are given by a theoretically derived DARKexp, whose radial profile shape is known exactly. For ϕ₀ ≈ 4.5, this shape is very closely matched by Einasto profiles, which is why Einasto fit simulated halos very well. In Section 3.2.1 we represent density profiles by Einasto profiles, giving us a total of four model parameters for (i)–(iii). In Section 3.2.2 we assume the density profiles can come from a wider range of DARKexp models, and represent them with three Einasto segments, parameterized by three parameters, giving us a total of six model parameters for (i)–(iii).

If Equation (1) is a robust property of collisionless Newtonian gravity, then one expects this form to be very closely adhered to by equilibrium dark matter halos. Put differently, it would be pointless to try to explain the radial dependence of ρ/σ³ if it is not a nearly exact power law. In this paper we do not study dynamical evolution of halos; only their final equilibrium configurations. We explain the principle behind our analysis later in this section.

Of the relations (i)–(iii), the only one that does not have a nearly exact form is the velocity anisotropy profile, (ii). Here we have to rely on simulations. Fortunately, as we show in Section 3.2 the conclusions of our analysis do not depend on the exact shape of this relation. Hansen & Moore (2006) and Hansen & Stadel (2006) have shown based on a variety of initial conditions, that equilibrium halos have a tight linear relation between velocity anisotropy β, and the double logarithmic density profile slope, γ. Hansen & Moore (2006) give a range of parameters characterizing that relation, while Hansen & Stadel (2006) present a single set of best-fitting parameters. We use both of these results in our analysis.

At the end of their dynamical evolution, simulated halos attain a state where all three relations, (i)–(iii), have nearly exact parametric forms, represented with four or six parameters. The values of these parameters have been measured from simulations; let us call them collectively parameter set A. In addition to A, one can define another parameter set, B, which for the same set of parametric relations, solves the Jeans equation with the smallest residuals.

We can explain the meaning of set B as follows. Suppose you were told that the equilibrium halos had Einasto-like density profiles, pseudo-phase-space density profiles had power-law shape as a result of some physical principle, and anisotropy–density slope relation was approximately linear. You were then asked to obtain the parameter values characterizing these relations. The obvious way for you to proceed would be to solve the Jeans equation incorporating these relations and obtain the best-fitting set of parameters; that would be set B.

It is not a foregone conclusion that A and B are the same set. First, let us consider a situation—possibly hypothetical—where Equation (1) is a property of Newtonian dynamical evolution and characterizes all equilibrium halos. The premise of our analysis is that in this case, the two sets of parameters will be the same, i.e., evolution will find the parameter set for which ρ/σ³ is as close to a power law as it can be. On the other hand, if Equation (1) does not have a physical origin, parameter sets A and B need not be the same, and the fact that Equation (1) appears to be satisfied in simulated halos is a coincidence.

To assess the similarity of the two parameter sets, A and B, we evaluate the quality of solutions from a wide and finely sampled region of the model parameter space. Since all three relations, (i)–(iii), are of equal importance, one should treat them equally. This is not possible if one assumes exact forms for (i) and (ii), and then integrates the Jeans equation to get (iii). Fortunately, there is a way to place (i)–(iii), and their associated four or six model parameters, on the same footing. In Section 3.1 we show that the second derivative of the Jeans equation, combined with an Einasto profile, can be converted to a cubic algebraic equation. We then calculate how well each set of parameters satisfies this algebraic equation. In contrast to the integration of the Jeans equation, the algebraic equation does not single out ρ/σ³.

Since (i)–(iii) have been observed in simulations, we know that a rough agreement between sets A and B is guaranteed. Therefore, in order to provide support for Equation (1) having a physical meaning, we are looking for better than a rough agreement. In that case, we also expect set B to form a well defined and isolated trough in the parameter space and to be stable against small changes in parameterization, like changing from four to six parameters.

We note that our entire analysis is done with the assumption of spherical symmetry of halos, and is one of the limitations of our modeling. Real dark matter halos, even geometrically spherically symmetric equilibrium halos can have velocity structures that are not fully described by radial velocity anisotropy (Wojtak et al. 2013). However, most (if not all) analysis of equilibrium halos in N-body simulations is done after spherically averaging their properties, like density, radial and tangential velocity dispersions, and pseudo-phase-space density. Because the results of these analyses form the starting point of our analysis, we are necessarily confined to the case of spherical symmetry.

In this section we will work with the second derivative of the Jeans equation, and parameterized forms of the density profile, velocity anisotropy profile, and power-law profile of pseudo-phase-space density.

We start with the anisotropic Jeans equation,

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dr}}\left\{\displaystyle \frac{\rho (r){\sigma }^{2}(r)}{3-2\beta (r)}\right\}+\displaystyle \frac{2\beta (r)}{3-2\beta (r)}\displaystyle \frac{\rho (r){\sigma }^{2}(r)}{r}=-G\rho (r)\displaystyle \frac{M(r)}{{r}^{2}},\,\end{eqnarray} \tag{ 2 }$$

where M(r) is the mass enclosed within radius r, ρ is the density at that radius, σ is the total velocity dispersion,² and β is the anisotropy, defined using the tangential, ${\sigma }_{t}^{2}=\tfrac{1}{2}({\sigma }_{\phi }^{2}+{\sigma }_{\theta }^{2})$, and radial velocity dispersions of a dark matter halo: $\beta (r)\equiv 1-\tfrac{{\sigma }_{t}^{2}}{{\sigma }_{r}^{2}}$. We define dimensionless variables x ≡ r/r₀ and y ≡ ρ/ρ₀, and reduce the number of functions in Equation (2) by invoking the power-law nature of the pseudo-phase-space density profile, $\rho /{\sigma }^{3}\,={({\rho }_{0}/{\sigma }_{0}^{3})(r/{r}_{0})}^{-\alpha }$:

$$\begin{eqnarray}&&\displaystyle \frac{-{x}^{2}}{y}\left\{\displaystyle \frac{d}{{dr}}\left[\displaystyle \frac{{y}^{5/3}{x}^{2\alpha /3}}{3-2\beta (r)}\right]+\displaystyle \frac{2\beta (x)}{3-2\beta (x)}{y}^{\tfrac{5}{3}}{x}^{\tfrac{2\alpha }{3}-1}\right\}={BM}(x)\end{eqnarray} \tag{ 3 }$$

where $B=G/{r}_{0}{v}_{o}^{2}$. Using the assumption of Equation (1) leads to Jeans equation constrained. From now on, ACJ refers to ACJ equation.

3.1. Second Derivative of the ACJ Equation

We differentiate Equation (3) with respect to x, as was done in Taylor & Navarro (2001). Then following Williams et al. (2004) and Barnes et al. (2007), we differentiate it again with respect to x, arriving at

$$\begin{eqnarray}&&\begin{array}{l}(2\alpha +\gamma -6)\left[\displaystyle \frac{2}{3}(\alpha -\gamma )+1\right](2\alpha -5\gamma )\\ =\,15\gamma ^{\prime\prime} +3\gamma ^{\prime} (8\alpha -5\gamma +4\beta +12\beta \theta {b}_{1}-5)\\ -\,3\theta [{b}_{1}(4{\alpha }^{2}+{\gamma }^{2}-8\alpha \gamma +8\alpha +7\gamma -15)]\\ -\,3{\theta }^{2}[6{b}_{1}{b}_{2}(\alpha -\gamma +1)]\\ -\,3{\theta }^{3}[{b}_{3}(54\beta +144{\beta }^{2}+24{\beta }^{3})]\\ -\,3\theta ^{\prime} [6{b}_{1}(\alpha -\gamma +1)+9{b}_{1}{b}_{2}\delta ]-3\theta ^{\prime\prime} (3{b}_{1}),\end{array}\end{eqnarray} \tag{ 4 }$$

where $\gamma =-d\mathrm{ln}y/d\mathrm{ln}x$ is the double logarithmic density slope, $\theta =d\mathrm{ln}\beta /d\mathrm{ln}x$, ${b}_{1}=2\beta /(3-2\beta )$, ${b}_{2}=(3+2\beta )/(3-2\beta )$, ${b}_{3}={(3-2\beta )}^{-3}$, and the primes indicate logarithmic derivatives. Using a range of initial conditions for collisionless N-body simulations, Hansen & Moore (2006) and Hansen & Stadel (2006) find that after equilibrium is achieved, all the halos end up having very similar shapes of the velocity anisotropy. All are well described by a linear relation³ between β and γ:

$$\begin{eqnarray}&&\beta (r)={\eta }_{1}+{\eta }_{2}\gamma (r).\end{eqnarray} \tag{ 5 }$$

This allows us to compute the logarithmic derivatives of θ:

$$\begin{eqnarray}&&\theta =-\displaystyle \frac{{\eta }_{2}}{\beta }\gamma ^{\prime} \end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&\theta ^{\prime} =-\displaystyle \frac{{\eta }_{2}}{\beta }\left\{\displaystyle \frac{{\eta }_{2}}{\beta }{(\gamma ^{\prime} )}^{2}+\gamma ^{\prime\prime} \right\}\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}&&\theta ^{\prime\prime} =-\displaystyle \frac{{\eta }_{2}}{\beta }\left\{\displaystyle \frac{2{(\gamma ^{\prime} )}^{2}\,{\eta }_{2}^{2}\,\gamma }{{\beta }^{3}}+2\gamma ^{\prime} \gamma ^{\prime\prime} +\displaystyle \frac{\gamma ^{\prime\prime} }{{\beta }^{2}}{\eta }_{2}\gamma +\displaystyle \frac{\gamma \prime\prime\prime }{\beta }\right\}.\end{eqnarray} \tag{ 6c }$$

Equation (4) is a second-order differential equation in γ and θ. It can be further expanded by using logarithmic derivatives in Equation (6), thereby converting the Equation (4) into a third-order differential equation in one variable only, γ. Due to the complexity of Equations (4) and (6), we are not showing the combined expression obtained after assembling all the parts together.

3.1.1. Single Einasto Algebraic (SEA) Equation

The logarithmic density γ can take many forms; a simple analytical expression that fits N-body density profiles well and is commonly used in the literature is the Einasto profile (Einasto 1965; Navarro et al. 2004):

$$\begin{eqnarray}&&\gamma ={{Ax}}^{p},\end{eqnarray} \tag{ 7 }$$

where A is a normalization constant. It has an interesting, and useful property that all its logarithmic derivatives have a linear dependence on γ:

$$\begin{eqnarray}&&\gamma ^{\prime} ={{Apx}}^{p}=p\gamma \end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&\gamma ^{\prime\prime} ={{Ap}}^{2}{x}^{p}={p}^{2}\gamma \end{eqnarray} \tag{ 8b }$$

$$\begin{eqnarray}&&\gamma \prime\prime\prime ={{Ap}}^{3}{x}^{p}={p}^{3}\gamma .\end{eqnarray} \tag{ 8c }$$

We take advantage of this property of the Einasto profile. In Equations (4) and (6) we eliminate all derivatives of γ by replacing them with their counterparts in Equation (8). This converts differential Equation (4) into an algebraic equation of cubic order. From now on, we will refer to it as the SEA equation that consists of a collection of Equations (4)–(8). It depends on four parameters: the power-law slope of the pseudo-phase-space density α, the two velocity anisotropy parameters, η₁, η₂, and the Einasto parameter p.

3.1.2. Triple Einasto Algebraic (TEA) Equation

Though Einasto profiles fit N-body dark matter halo density profiles quite well, they are still only fitting functions. DARKexp models have the advantage of being theoretically derived from fundamental statistical mechanical principles. They also fit the density profiles, and the energy distributions of N-body halos very well. It is not surprising that they resemble Einasto profiles for a certain range of DARKexp shape parameter ϕ₀. DARKexp density profiles do not have an analytic expression, but since for the relevant range of ϕ₀ they are not too different from the Einasto shape, DARKexp can be approximated by three Einasto segments, with three different slopes p₁, p₂, and p₃, respectively. The profiles we use cover three decades in radius, from x = 10⁻² to x = 10¹, and x = 1 corresponds to the radius where γ = 2, i.e., has the isothermal slope. Note that the shape of the profile at the very center, at x < 10⁻² is irrelevant since it is excluded from the analysis.

$$\begin{eqnarray}\gamma (x)=\left\{\begin{array}{ll}2\times {10}^{1.3({p}_{1}-{p}_{2})}\,{x}^{{p}_{1}}, & \mathrm{log}x\leqslant -1.3\\ 2\,{x}^{{p}_{2}}, & -1.3\leqslant \mathrm{log}x\leqslant 0\\ 2\,{x}^{{p}_{3}}, & 0\leqslant \mathrm{log}x.\end{array}\right.\end{eqnarray} \tag{ 9 }$$

The continuity conditions are derived at $\mathrm{log}x=-1.3$ and log x = 0 by forcing the density slopes γ to be equal at these two points, with the help of a constant multiplier in the logarithmic units. We also require that p₂ ≤ p₁, in accordance with the shape of DARKexp density profiles (Williams et al. 2010).

When Equation (9) is substituted into Equation (4), the final equation becomes a function of six model parameters: the power-law slope of the pseudo-phase-space density α, the two velocity anisotropy parameters, η₁, η₂, and three Einasto parameters p₁, p₂, and p₃, which collectively approximate DARKexp. We call it the TEA equation, which consists of Equations (4)–(9).

3.1.3. Figure of Merit for Evaluating the Solutions

If the simple analytical assumptions—Equations (1), (5), and (7) or Equation (9)—are exactly correct, then the SEA and TEA equations will be exactly true at all radii, x. However, this is not the case for any set of four or six parameters. To estimate how well any given set of parameters solves the two equations over three decades in radius, we divide this entire range into 30 logarithmic intervals, calculate the normalized absolute difference between the left-hand side (LHS) and the right-hand side (RHS) of the equations, and then average these values. We chose log intervals because that is what is used by other researchers when fitting density and pseudo-phase-space profiles to N-body simulated halos. The resulting value, which we call δ, is defined as,

$$\begin{eqnarray}&&\delta ={\left\langle \displaystyle \frac{| \mathrm{LHS}-\mathrm{RHS}| }{\sqrt{({\mathrm{LHS}}^{2}+{\mathrm{RHS}}^{2})/2}}\right\rangle }_{\mathrm{log}x}.\end{eqnarray} \tag{ 10 }$$

After scanning a large portion of the parameter space, we will identify regions that represent best solutions of the two equations (global minima of δ) and compare the corresponding parameter values to those found in computer simulations.

3.2. Grid-based Search for Approximate Solutions

3.2.1. SEA Equation

Here, we consider the single Einasto profile of Equation (7), resulting in a four-dimensional parameter space. The parameter ranges covered are $1.8\leqslant \alpha \leqslant 2.2$, −0.60 ≤ η₁ ≤ 0.40, −0.5 ≤ η₂ ≤ 0.36, 0.02 ≤ p ≤ 0.3, with a step size of 0.02 for all parameters. Our parameter range for α values is wider than what is found in the literature, which generally span ∼1.85–1.96. The range of density profile slopes is also somewhat wider. The velocity anisotropy parameter ranges encompass those found in N-body simulations, as analyzed in Hansen & Moore (2006).

We start by calculating the figure of merit, δ, Equation (10), at every location in the four-dimensional parameter space. Given our step size and parameter limits, we consider a total of around 5 × 10⁵ parameter sets. The parameter sets with the lowest δ values correspond to the best solutions to the SEA equation. Figure of merit δ is always greater than zero, meaning that no solution is exact. As was pointed out in Section 2, the advantage of using the second derivative of the Jeans equation is that, by design, it contains no information about which one of the three analytic assumptions—density, anisotropy, or ρ/σ³ profile, or what combination of the three—is the reason for δ being nonzero. In the rest of this subsection we study the interdependence of these parameters, examine the properties of the best solutions, and compare them to those of dark matter halos found in simulations.

To display the solution space, we chose the density profile slope p corresponding to the lowest global δ value (p = 0.12), and plot four cuts through the remaining 3D parameter space, corresponding to four values of α: 1.8, 1.9, 2.0, and 2.1; see Figure 1. The color scale in these heat maps is the same for all four panels. The gray color indicates very large values of δ, outside our range. The black box outlines the anisotropy parameters given in Hansen & Moore (2006), while the black dot marks the single set presented in Hansen & Stadel (2006); both are based on simulated halos.

Zoom In Zoom Out Reset image size

The pattern revealed in these approximate solutions is complicated. There is a symmetry axis outlined by white pixels, where the denominator of δ is zero for some values of x. It goes roughly diagonally through all four panels, and the color pattern is inverted across the axis. The trough of minima is to the right of this axis for α = 1.8–1.9, and shifts to the left for large values of α.

The global δ minimum is at $(\alpha ,p,{\eta }_{1},{\eta }_{2})\,=(1.80,0.12,0.38,0.30)$ i.e., set B for the SEA equation, and does not correspond closely to the parameters of the simulated halos. Specifically, the anisotropies are too large, even exceeding β = 1 (which is not allowed because of the definition of β), the Einasto density profile slope, p, is too small, and the slope of the pseudo-phase-space density profile, α is a little shallower than the range seen in simulations.

To quantify the difference between set B parameters and those obtained in N-body simulations, set A, we use Ludlow et al. (2011), who find that the slope of the Einasto density profile slope, p, and the exponent of the pseudo-phase-space density, α, are linearly correlated; see their Figure 7. We estimate the uncertainty, σ, in these set A values from that figure: the dispersion around the linear relation is ${\sigma }_{p\alpha }\approx 0.015$. Here, pα is not a product, but is meant to represent that the two parameters are correlated. We also assume that the two anisotropy parameters are not correlated. Their uncertainty is taken to be half the range given in Hansen & Moore (2006), i.e., for the anisotropy parameter η₁, the uncertainty is ${\sigma }_{{\eta }_{1}}=0.27$, and for η₂, it is ${\sigma }_{{\eta }_{2}}=0.12$.

For case SEA we estimate the significance of the difference between set A and set B parameters by vectorizing

$$\begin{eqnarray*}&&{\sigma }_{\mathrm{SEA}}=\left(\displaystyle \frac{| p{\alpha }_{A}-p{\alpha }_{B}| }{{\sigma }_{p\alpha }},\displaystyle \frac{| {\eta }_{{1}_{A}}-{\eta }_{{1}_{B}}| }{{\sigma }_{{\eta }_{1}}},\displaystyle \frac{| {\eta }_{{2}_{A}}-{\eta }_{{2}_{B}}| }{{\sigma }_{{\eta }_{2}}}\right),\end{eqnarray*}$$

where $| p{\alpha }_{A}-p{\alpha }_{B}|$ is the distance in the plane of p versus α of our set B parameters and the closest set A parameters in Figure 7 of Ludlow et al. (2011). We find that ${\sigma }_{\mathrm{SEA}}=7.4$, hence the two sets are significantly different from each other.

The best 150 solutions in the space of the anisotropy–density relation are displayed in Figure 2(a). The global minimum is the thick purple line, and it is very different from the anisotropy–density relation found in simulations. The thin lines (150 solutions) are color coded by their α value. The range found in simulations, $1.85\leqslant \alpha \leqslant 1.96$ (blue) does not occur in this set.

Zoom In Zoom Out Reset image size

The fact that the parameter set corresponding to the global δ minimum (set B) does not coincide with the one observed in simulations (set A) is an indication that ρ/σ³ is not as close to a power law as it could have been, and hence does not have a physical significance. This conclusion is insensitive to the exact shape of the anisotropy–density relation, and remains valid even if the relation is only roughly linear. This is because the global δ minimum is outside of the range obtained from simulations (black box in Figure 1).

3.2.2. TEA Equation

In this section we extend the range of density profiles we consider, by including profiles that resemble DARKexp, which can be approximated by three joined Einasto-like segments, Equation (9).

The corresponding hexa-dimensional parameter space we search spans the same range for α, η₁, and η₂ as in Section 3.2.1. The values of the slope of the two Einasto segments of the density profile at smaller radii, p₁ and p₂, cover the same range as p, but the segment that applies to larger radii spans $0.16\leqslant {p}_{3}\leqslant 0.24$, because the density profile is expected to be steeper there. We impose an additional constraint that p₂ ≤ p₁, as indicated by the shape of DARKexp density profiles. The step size in each of the six dimensions is 0.02, so the total number of parameter sets we consider is nearly 2.5 × 10⁷.

Figure 3 is similar to Figure 1, but only shows η₁ and η₂ ranges that are indicated by black rectangles in the latter figure. This is the velocity anisotropy range found in simulations (Hansen & Moore 2006). The four panels correspond to the same four α values as in that figure, while the three values of the density profile slopes of triple Einasto are those that resemble the density profile of DARKexp with ϕ₀ = 4.5.

Zoom In Zoom Out Reset image size

Just like the δ maps in Figure 1, the ones in Figure 3 also show complicated patterns. (Note that the color scale is different in the four panels.) The global δ minimum, i.e., our set B, is at $(\alpha ,{p}_{1},{p}_{2},{p}_{3},{\eta }_{1},{\eta }_{2})=(1.94,0.2,0.18,0.16,0.06,-0.08)$. Its parameter values are not the same as those found in simulations. The α value is a little too large, the density profile is only marginally well approximated by either Einasto or DARKexp, and the anisotropy β is mildly tangentially anisotropic at large radii in contrast to the radially anisotropic velocity profiles of simulated halos. We quantify the difference between set A and set B parameters using the same standard deviation analysis described in Section 3.2.1. Since published papers on N-body simulations present only single, not triple, Einasto fits, we use the same result from Ludlow et al. (2011) for case TEA. This is acceptable because set B values for $({p}_{1},\,{p}_{2},\,{p}_{3})$ are very close to each other, (0.20, 0.18, 0.16), so we can use the average, p = 0.18. With this, set A and set B differ at the 3.3σ level.

Figure 2(b) shows the anisotropy–density relation for the 150 best solutions. The global minimum is shown as a thick purple line. The average β − γ relation from simulations (Hansen & Stadel 2006) is plotted for reference as the thick black line. Most, if not all of the best 150 solutions have properties different from those of simulated equilibrium halos. All green lines and the thick purple line have an anisotropy–density slope opposite to that seen in simulations, and α values that are larger than those seen in simulations; α ≥ 1.96. The solutions represented by red lines have α values that are too small, α ≤ 1.85. Furthermore, many of these have β > 1 in the relevant range of density slopes. The solutions that have α in the observed range, 1.85 ≤ α ≤ 1.96 (blue lines) have approximately isotropic velocity distributions at all radii, and are thus only marginally consistent with simulations.

Figure 4 shows the density profile slopes versus log radius, for the 150 best solutions, color coded by α as in Figures 2(b), and 4(b) plots the subset of these 150 solutions that are within the Hansen & Moore (2006) anisotropy range. Though we allowed a rather broad range of slopes for p₁, p₂, and p₃, the best solutions are all clustered around a narrow range of values, resulting in considerable overlap of profiles in the figure. These values are such that ${p}_{1}\approx {p}_{2}\approx 0.16\mbox{--}0.18$, and p₃ ≈ 0.18–0.20. In other words, even though the density profile had the option of deviating from an Einasto form, the best solutions still have p₁ ≈ p₂ ≈ p₃. The closest DARKexp has ϕ₀ = 4.5, and is represented as a light blue line. So the density profiles corresponding to best solutions are similar to, but not the same, as those found in N-body simulations.

Zoom In Zoom Out Reset image size

To sum up, the parameter set found in simulations is (i) somewhat similar to, but (ii) not very close to, the best solution to the ACJ equation. The first statement is an a posteriori conclusion, which was already shown to be the case by Taylor & Navarro (2001). The second statement suggests that the power-law nature of ρ/σ³ found in simulations is a coincidence. If Equation (1) had physical significance, nature would have found a different set of parameters—namely, our set B—to satisfy Equation (1), while keeping density profiles Einasto or DARKexp-like, and anisotropy–density slope relation approximately linear.

A further conclusion drawn from these results argues against Equation (1) having a physical significance. If it were, one would expect the parameters of the single Einasto (Section 3.2.1) and triple Einasto (Section 3.2.2) best solutions to be similar. This is not the case. While the density profiles in both cases are Einasto-like (for the triple Einasto ${p}_{1}\approx {p}_{2}\approx {p}_{3}$), the slopes are different, and the power-law exponents α, and the velocity anisotropies are very different. Thus, a relatively small change in the parameterization of the problem significantly changes the parameters of the global minimum. The unstable solution suggests that Equation (1) is not a universal feature of dynamical evolution.

It has been known for almost 20 yr that the pseudo-phase-space density profiles of equilibrium dark matter halos are well approximated by a power law in an over ∼3 decade radius. The main goal of this paper is to determine whether this scale-free behavior could have a physical origin, or if it is simply a curious coincidence. While we do not address the possible physical meaning of ρ/σ³ being a power law, we assume that if one exists, equilibrium halos will obey the Equation (1) relation.

We work with the final equilibrium halos, instead of analyzing their dynamical evolution. We use the Jeans equation of hydrostatic equilibrium, and the three simple parametric relations that describe the radial behavior of three profiles: density, velocity anisotropy profile, and ρ/σ³, which are parameterized by a total of four or six model parameters, depending on whether the density is described by a single or triple Einasto profile. We then search the parameter space, over a wide and finely sampled range, for the best solution to the Jeans equation (set B). Because all three relations are seen in simulations, we already know that parameters describing N-body halos (set A) will correspond to reasonably good solutions of the Jeans equation. However, if ρ/σ³ ∝ r^−α has a physical origin, we expect N-body halos to correspond to the best solutions (i.e., set A and B to be the same within uncertainties), and expect these to form a well defined, isolated, and stable global minimum trough.

What we actually find is quite different. The structure of the solution space is complicated, with no indication of an isolated trough. The parameter set found in N-body simulations (set A) and the one that corresponds to a global minimum we find (set B) are different at ∼3–7σ significance. Since numerical work like ours can explore only a limited range of parameters, and will necessarily leave large portions of parameter space unexplored, we need to consider the possibility that the true global minimum (set B) is outside of our range. That means that set B is even further away from set A, which only strengthens our conclusion that the two sets are not the same.

Furthermore, we find that if the parameterization of the density profile is changed somewhat, from a single to triple Einasto, the parameters of the best solution, and especially those of velocity anisotropy change significantly indicating that the global minimum in the solution space is not stable.

Because the parameters describing the density, velocity dispersion, and ρ/σ³ radial profiles observed in simulations are not the ones that would result in the ρ/σ³ being closest to a power law, we conclude that the approximate power-law nature of the pseudo-phase-space density seen in N-body simulations and semianalytical collapses does not have a physically meaningful origin, and so does not shed light on the effective equation of state of self-gravitating dark matter halos.