Dark Matter Halos: The Dynamical Basis of Effective Empirical Models

We investigate the dynamical basis of the classic empirical models (specifically, Sersic-Einasto and generalized NFW) that are widely used to describe the distributions of collisionless matter in galaxies. We submit that such a basis is provided by our \alpha-profiles, shown to constitute solutions of the Jeans dynamical equilibrium with physical boundary conditions. We show how to set the parameters of the empirical in terms of the dynamical models; we find the empirical models, and specifically Sersic-Einasto, to constitute a simple and close approximation to the dynamical models. Finally, we discuss how these provide an useful baseline for assessing the impact of the small-scale dynamics that may modulate the density slope in the central galaxy regions.


Introduction
The classic Sérsic (1963) models met a wide and lasting success as empirical representations of the projected (2-dimensional) light distributions in spheroidal galaxies (for a review, see Kormendy et al. 2009). Einasto (1965) developed and used a similar shape to describe in simple terms 3dimensional stellar mass profiles.
On the other hand, recent extensive N -body simulations (e.g., Navarro et al. 2004;Merritt et al. 2005;Gao et al. 2008;Stadel et al. 2009;Navarro et al. 2010) indicate that the Sérsic and Einasto functional forms also provide good patterns to represent the spherically-averaged mass distributions in dark matter (DM) halos ranging from galaxies to galaxy clusters. These apply at levels comparable to, or even better than the popular NFW formula (Navarro, Frenk & White 1997).
Still, no agreed understanding is available to explain the value in both the real and the virtual world of the Sérsic and Einasto representations (see discussions by Graham et al. 2006;Kormendy et al. 2009). Can we identify the underlying as-trophysical basis?

Density runs
The density runs of the SE family may be represented in the form ρ(r) =r −τ e −u (r η −1) , u = 2 − τ η . (1) Here, quantities are normalized to their value at r −2 , the reference radius where the logarithmic slope γ ≡ −d log ρ/d log r takes on the value 2; typically, in nearby elliptical galaxies r −2 corresponds to sizes of order 10 kpc, a few times the half-light radius R e . The parameters τ and η describe the inner slope and the middle curvature of the density run, re-spectively. The original Einasto profile belongs to this family, and is obtained when τ = 0. Note, however, that by deprojecting from the plane of the sky a Sérsic 2-dimensional run e −s 1/n with index n ≈ 3−4 (suited for normal ellipticals, see Kormendy et al. 2009) produces a cuspy inner run as in Eq. (1) with τ ≃ 1 − 1.19/2n + 0.22/4n 2 ≈ 0.8 significantly different from 0 and less than 1, as shown by Prugniel & Simien (1997).
On the DM side, recent simulations (see Gao et al. 2008;Stadel et al. 2009;Navarro et al. 2010) only provide an upper bound τ < 0.9 for the inner slope. When the original Einasto profile (with τ = 0) is adopted, the best-fit to simulated DM halos obtains for η ≈ 0.2; we will come back to this value later on.

Toward a single family
The main apparent difference between SE and gNFW is constituted by the former's exponential decline vs. the latter's powerlaw falloff ρ ∝ r −(τ +ηξ) for large r.
On the other hand, Eq.
(2) is to be considered for large values of ξ anyway, since a steep density run in the halo outskirts is indicated by observations of light distribution in spheroidal galaxies (other than cDs, see Kormendy et al. 2009), and of DM distributions from weak lensing in galaxies and galaxy clusters (e.g., Broadhurst et al. 2008;Oguri et al. 2009;Newman et al. 2009).
The circumstance is easily translated into the formal statement that the gNFW family converges 1 In the literature these runs are sometimes referred to as αβγ-models, and equivalently defined via the parameters γ = τ , α = 1/η, β = τ + ηξ.
to the SE for large ξ. This is seen on recastingρr τ from Eq. (2) in exponential form, to read (3) for approximating the middle and last terms we have used the circumstance that ξ ≫ 1 implies w ≫ 1 and so ξ w ≃ (2 − τ )/η ≡ u applies. Thus the two families in Eqs. (1) and (2) actually become one in this limit.
Thus in the following we focus mainly on the SE family, and proceed to discuss its dynamical basis in terms of the Jeans equation.

The Dynamical model
The dynamical model of DM halos hinges upon the radial Jeans equation that expresses the selfgravitating, equilibrium of collisionless matter (see Binney & Tremaine 2008). The Jeans equation reads in terms of the density ρ(r), the related cumulative mass M (< r) ≡ 4π r 0 dx x 2 ρ(x), and the radial velocity dispersion σ 2 r (r). The last term on the r.h.s. describes the effects of anisotropic random velocities via the standard Binney (1978) parameter β ≡ 1 − σ 2 θ /σ 2 r . Note that the Jeans equation is designed to describe a (quasi-)static equilibrium, away from extreme major merger events like is the case with the Bullet Cluster (see Clowe et al. 2006). But even in relaxed conditions, solving Jeans requires an 'equation of state', i.e., a functional relation expressing the DM pressure ρ σ 2 r in terms of density (and possibly radius) only.

Equation of state
In seeking for such a relation, one can make contact with the classic theory of the non-linear collapse for DM halos in an expanding Universe; here self-similar arguments play the role of a pivotal pattern (see Fillmore & Goldreich 1983;Bertschinger 1985;Taylor & Navarro 2001). This still applies to modern views of the halo development (e.g., Mo & Mao 2004;Lu et al. 2006;Li et al. 2007;Lapi & Cavaliere 2009a,b;Fakhouri et al. 2010), that comprise two stages: an early collapse builds up the halo main body via a few major merger events and sets its phase-space structure by dynamical relaxation of DM particle orbits; this tails off into a secular development of the outskirts by smooth accretion and minor mergers.
The essence of the macroscopic equilibrium is conveyed by the self-similar scaling σ 2 r ∝ M/r adding to the geometric relation ρ ∝ M/r 3 . The macroscopic import of the halo phase-space structure is conveyed by combining these two quantities into the 'phase-space density' ρ/σ 3 r , or equivalently into the functional K(r) ≡ σ 2 r /ρ 2/3 often referred to as DM 'entropy' (see Bertschinger 1985;Taylor & Navarro 2001). For the latter quantity, one easily derives the scaling K(r) ∝ r M 1/3 implying K(r) ∝ r α ; whence one expects a slope α slightly exceeding unity.
To focus the values of α, in Lapi & Cavaliere (2009a,b) we have developed a full semianalytic treatment of the halo growth in the standard accelerating Universe (see Komatsu et al. 2010). We found constant values of α, that fall within the narrow range 1.25 − 1.3; on average, such values grow weakly with the mass of the halo body, from galaxies to rich clusters.
The halo development process has been probed, and the two-stage view confirmed by intensive Nbody simulations (e.g., Zhao et al. 2003;Wechsler et al. 2006;Diemand et al. 2007;Hoffman et al. 2007;Schmidt et al. 2008;Ascasibar & Gottlöber 2008;Vass et al. 2009;Genel et al. 2010;Wang et al. 2010). These also confirm that: (i) a (quasi-)static macroscopic equilibrium is attained at the end of the fast collapse, and is retained during the subsequent stage of secular, smooth mass addition; (ii) a persistent feature of such an equilibrium is constituted by powerlaw correlations holding in the form σ 2 r /ρ 2/3 ∝ r α , although it is still widely debated whether the radial or the total velocity dispersion best applies (see also the discussions by Schmidt et al. 2008 andby Navarro et al. 2010).
In building up our dynamical models we focus on the quantity K ≡ σ 2 r /ρ 2/3 ∝ r α that involves the radial dispersion σ 2 r (see also Dehnen & McLaughlin 2005). Operationally, this provides a direct expression for the radial pressure term ρ σ 2 r = K ρ 5/3 ∝ r α ρ 5/3 in the Jeans Eq. (4); anisotropies are accounted for by the last term on the r.h.s., as discussed in § 3.3 below.

The DM α-profiles
In terms of K(r) ∝ r α , the Jeans equation may be recast into the compact form with γ ≡ −d log ρ/d log r representing the logarithmic density slope and v 2 c ≡ GM (< r)/r the circular velocity. Remarkably, by double differentiation this integro-differential equation for ρ(r) reduces to a handy 2 nd order differential equation for γ (see Austin et al. 2005;Dehnen & McLaughlin 2005).
Tackling first the isotropic case β = 0, we recall that the solution space of Eq. (6) spans the range α ≤ 35/27 = 1.296: the specific solution for the upper bound, and the behaviors of the others have been analytically investigated by Taylor & Navarro (2001), Austin et al. (2005) and Dehnen & McLaughlin (2005). In Lapi & Cavaliere (2009a), we explicitly derived the solutions in the full range α = 1.25−1.296, that are marked by a monotonically decreasing run and satisfy physical boundary conditions: a finite central pressure or energy density (equivalent to a round minimum of the gravitational potential); a steep outer run implying a finite and rapidly converging (hence a definite) overall mass. We dubbed α-profiles these physical solutions.
We shall use the following basic features of the latter. In the halo body at the point r 0 the αprofile is tangent to the pure powerlaw solution ρ ∝ r −γ0 of the Jeans equation; there v 2 c /σ 2 r ∝ r 2−γ0/3−α applies, to imply from Eq. (6) This is consistent for α = 1.25 with the self-similar slope; as such, it qualifies to provide a universal middle-range slope. Note that the point r 0 lies in the neighborhood of the radius r −2 (see § 2.1), specifically r 0 ≈ 1.74 − 1.51 r −2 holds for α ≈ 1.25 − 1.3. On the other hand, a monotonic density run implies the term v 2 c /σ 2 r ∝ r 2−γ/3−α to vanish at the center; this results in an inner powerlaw ρ ∝ r −γa with the slope γ a = 3 5 α .
This differs from zero as long as the entropy run grows from the center following K ∝ r α with α > 0. Finally, a finite mass implies v 2 c /σ 2 r ∝ r −1+2γ/3−α to hold in the outskirts, so as to yield a typical outer decline ρ ∝ r −γ b with slope This exceeds the value 3, and so constitutes the hallmark of a rapidly saturating mass; the circumstance makes less compelling here the role of a virial boundary.
Thus, compared to NFW the inner slope of the dynamical model is considerably flatter and the outer slope steeper; compared to the original Einasto profile, the main difference occurs in the inner regions where the dynamical model is (moderately) steeper.

Anisotropy
It is clear from Eq. (6) that anisotropies will steepen the density run for positive β, and flatten it for negative β. The latter condition is expected to prevail in the inner region, where tangential components develop from the angular momentum barrier (Nusser 2001;Lu et al. 2006). Moving outwards, radial motions are expected to prevail, so raising β up to values around 0.5 at r ≈ r 0 ; outwards of this, β is expected to saturate or even decrease, as one enters a region increasingly populated by DM particles on eccentric orbits with vanishing radial dispersions at their apocenters (see Bertschinger 1985). This view is supported by numerical simulations (see Austin et al. 2005;Dehnen & McLaughlin 2005;Hansen & Moore 2006;Navarro et al. 2010), which in detail suggest the average anisotropy-density slope relation to hold with parameters β(0) ≈ −0.1, β ′ ≈ 0.2, and the constraint β(r) 0.5. Note that at all radii the inequality γ(r) ≥ 2 β(r) is satisfied; this has been conjectured to constitute a necessary In Lapi & Cavaliere (2009b) we extended the dynamical model to such anisotropic conditions in the full range α ≈ 1.25 − 1.3, inspired by the analysis of Dehnen & McLaughlin (2005) for the upper bound of α. We note that the latter is now slightly modified to 35/27 − 4 β(0)/27 ≈ 1.31; likewise, the point r 0 where γ = γ 0 = 6 − 3 α applies moves slightly inwards, so that r 0 ≈ 1.58 − 1.38 r −2 now holds.
The main outcome, however, is that the density profile is somewhat flattened at the center relative to the isotropic case; the inner slope now reads γ a = 3 5 α + 6 5 β(0) .
In particular, even a limited central anisotropy (corresponding to values β(0) ≈ −0.1) causes an appreciable flattening down to γ a ≈ 0.63 − 0.66 for α ≈ 1.25 − 1.3. On the other hand, we stress that such small phenomenological anisotropies near the center imply the radial σ 2 r and the total dispersions σ 2 = σ 2 r [1 − 2 β/3] to be very close.

From Dynamical to Empirical Models
Here we discuss how the parameters of the empirical profiles (see § 2) can be set based on our dynamical model (see § 3); in such conditions, it will turn out that such profiles constitute close  Table 1 and 2. The lower panels highlight the corresponding logarithmic density slopes.

Parameters from dynamics
First we consider the original Einasto profile (τ = 0 in Eq. 1), since this has been widely used in the context of DM halo simulations. Here, τ is fixed to 0, and the only free parameter is the curvature η. This we set by requiring the logarithmic density slope γ(r) = 2r η to equal γ 0 at the pointr 0 . So we find the expression that takes on values η ≈ 0.15 − 0.2, see Table 1 and 2; remarkably, these turn out to agree with those derived from fits of state-of-the-art N -body simulations in terms of the same Einasto density run, as performed by Navarro et al. (2010). On the other hand, the flat central slope of the Einasto profile is at variance with the value γ a = 3 α/5 given by our dynamical models based on Jeans; to wit, consistency between pure Einasto and Jeans would require at the center a flat entropy distribution and a vanishing pressure.
Actually, the simulations quoted in § 3.1 within their finite mass resolution provide only an upper  Fig. 1 for the profiles of circular velocity, and for the corresponding logarithmic slopes. limit τ < 0.9 to the central slope. This grants scope to the full SE family of Eq. (1).
The latter features two parameters, the inner slope τ and the middle curvature η. These we set by requiring the logarithmic density slope to equal γ a forr → 0, and γ 0 atr 0 ; so we find Thus we predict the central slope to take on values τ = 0.6 − 0.8 and the corresponding curvature parameter to take on values η = 0.2−0.3, see Table  1 and 2. It will be worth fitting the outcomes of N -body simulations based on these extended SE profiles with τ > 0. Finally, we report the corresponding results for the empirical gNFW family. This features three parameters: inner slope τ , middle curvature η, and strength of the outer decline ξ; these we set by requiring the logarithmic density slope to equal γ a forr → 0, γ 0 atr 0 , and γ b for r → ∞. So we find The parameters so determined are listed in Table 1 and 2 for both the isotropic and anisotropic conditions.

Results and comments
With the parameters focused as discussed in the previous subsection, Fig. 1 illustrates how the empirical compare with our dynamical models. We plot the density run of the latter (specifically for the α-profile with α ≈ 1.25 suitable to galactic halos), compared to those of the Einasto, SE, and gNFW models. The left and right panels refer to isotropic and anisotropic conditions, respectively; the popular NFW profile is also shown for reference. To make comparisons easier, we plot in the lower panels the corresponding logarithmic density slopes.
It turns out that the closest approximation to the dynamical model is provided by SE, which shares with it not only the central slope by construction, but also the body and the outer behaviors. The original Einasto profile provides an acceptable approximation in the middle and outer ranges, but not at the center, because of its flatness. On the other hand, the gNFW family provides an acceptable approximation in the inner and middle ranges, but not in the outskirts where its slope is too flat. Finally, the NFW profile provides an acceptable approximation only in the middle range.
Similar conclusions concern the profiles of circular velocities v 2 c (r) ≡ G M (< r)/r, that are analytically dealt with in the Appendix and illustrated in Fig. 2.
We stress that the handy SE representation is convenient in analyzing data in several contexts, including: the DM particle annihilation signal expected from the Galactic Center (see Lapi et al. 2010a); rotation curves of dwarf and normal spiral galaxies (see Salucci et al. 2007); individual and statistical properties of elliptical and spiral galaxies (see Cook et al. 2009); strong and weak gravitational lensing (see Lapi et al. 2009b), currently observed in clusters (e.g., Zitrin et al. 2010) and soon in massive elliptical galaxies (see discussion by Bradač et al. 2009); X-ray emission from the intracluster plasma (see Cavaliere et al. 2009;Fusco-Femiano et al. 2009, Lapi et al. 2010b.

Discussion
We first stress that the dynamical model (as well as its approximations in terms of empirical models) is in keeping with the basic features of standard DM, i.e., its cold and collisionless nature. In fact, it implies σ 2 r (r) → 0 for large r ≫ r −2 , a behavior expected in the outskirts for cold matter dominating the potential well.
At the inner end, with decreasing r we expect σ 2 r (r) to increase toward a maximum, correspond-ing to effective conversion of inflow kinetic into random energy. In fact, toward the center Jeans requires d log σ 2 r /d log r = γ − GM (< r)/r 2 → γ a to hold as the gravitational force vanishes there, to the effect that σ 2 r (r) ∝ r γa → 0. Concerning the collisionless nature of the DM, the boundary conditions at the center imply a finite, non zero pressure (and energy density), while a long collisional mean free path allows the pressure gradient dp/dr to diverge. Conversely, with a short mean free path λ the pressure gradient cannot diverge on scales r λ, where a finite σ 2 and a flatter γ apply. In fact, weakly collisional conditions have been proposed to explain the cored light profiles observed in many spheroidal galaxies (see Ostriker 2000).
On approaching the center of a galactic halo, one expects the basic dynamical model from largescale Jeans equilibrium to be altered to an increasing degree by small-scale dynamics and/or energetics related to baryons. These processes are specifically related to following issues: transfer of energy/angular momentum from baryons to DM during galaxy formation; scouring baryons by the energy feedback from central active galactic nuclei; any 'adiabatic' contraction of the baryons. Such issues will be briefly discussed in turn, with a warning that they enter increasingly debated grounds.

Energy/Angular Momentum Transfers
Flattening of the inner density profile may be caused by transfer of energy and/or angular momentum from the baryons to the DM during the galaxy formation process (see El-Zant et al. 2001;Tonini et al. 2006, Romano-Díaz et al. 2008. In detail, upon transfer of tangential random motions from the baryons to an initially isotropic DM structure, the density in the inner region is expected to behave as (Tonini et al. 2006) Thus for β < 0 the profile is flattened relative to the original γ a , down to the point of developing a core for β −γ a /2 (2 − γ a ) ≈ −0.3. However, a reliable assessment of the amount of angular momentum transferred from the baryons to the DM is still wanting, and would require aimed numerical simulations of better resolution than presently achieved.

Other processes on inner scales
Less agreed processes may affect galactic scales r 10 2 pc. For example, at the formation of a spheroid, central starbursts and a supermassive black hole may easily discharge enough energy (∼ 10 62 erg for a black hole mass M • ∼ 10 9 M ⊙ ) with sufficient coupling ( 1%) to blow most of the gaseous baryonic mass m ∝ r 3−γa /(3 − γ a ) out of the gravitational potential well. This will cause an expansion of the DM and of the stellar distributions (see Fan et al. 2008), that flattens the central slope.
In addition, binary black hole dynamics following a substantial merger may eject on longer timescales formed stars from radii r ≈ 10 (M • / 10 8 M ⊙ ) 1/(3−γa) pc containing an overall mass of a few times the black hole's, and so may cause a light deficit in some galaxy cores (see Merritt 2004;Lauer et al. 2007;Kormendy et al. 2009). A full discussion of the issue concerning cored vs. cusped ellipticals is beyond the scope of the present paper.

Adiabatic Contraction?
On the other hand, some steepening of the inner density profile may be induced by any 'adiabatic' contraction of the diffuse star-forming baryons into a disc-like structure, as proposed by Blumenthal et al. (1986) and Mo et al. (1998) but currently under scrutiny, see Abadi et al. (2010).
On the basis of the standard treatments, it is easily shown that in the inner region an initial powerlaw ρ(r i ) ∝ r −γa i is modified into ρ ∝ r −3/(4−γa) ; this yields typical slopes around 0.9, steeper than the original γ a ≤ 0.78 but still significantly flatter than 1. However, recent numerical simulations (see discussion by Abadi et al. 2010) suggest that the treatment of adiabatic contraction leading to Eq. (19) is likely to be extreme; actually, in the inner region the contraction is ineffective and the density slope hardly modified. Again, highly resolved N -body experiments are needed to clarify the issue.

Conclusions
We have discussed the dynamical basis of the Sérsic-Einasto empirical models, in terms of wellbehaved solutions of the Jeans equation with physical boundary conditions comprising: a finite central energy density, a closely self-similar body, a finite (definite) overall mass.
We find the SE profile to be particularly suitable to represent the general run of the dynamical solution. Specifically, we have discussed how to tune the parameters of SE in terms of the dynamical model; in such conditions, we find the former to constitute a simple and close approximation to the latter.
The resulting SE profile shares with the dynamical model the following features: an outer steep decline, hence a definite overall mass; a closely selfsimilar body with slope γ 0 = 6−3α; an inner slope around γ a = 3α/5, hence flatter than −1. The latter slope provides an useful baseline for discussing alterations of the inner behavior caused by additional baryonic processes.
In conclusion, we submit that the dynamical models discussed here, namely the α-profiles, provide the astrophysical basis for understanding the empirical success of the SE profiles in fitting the real and the virtual observables, from galaxies to galaxy clusters.