Free Access
Volume 522, November 2010
Article Number A28
Number of page(s) 9
Section Cosmology (including clusters of galaxies)
Published online 28 October 2010

© ESO, 2010

1. Introduction

In a previous work (Henriksen et al., referred to here as Paper I, and references hereafter) we discussed therelation between the formation of black holes (hereafter BH) and of galaxies. We presented distribution functions (DF) for the co-eval formation of spherically symmetric cusps and bulges, as formed by radially infalling, collisionless matter (e.g. dark matter or stars). We also discussed the influence of a centrally dominant mass, including a point mass BH.

Observations (Kormendy & Richstone 1995; Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000) have established a strong correlation between the BH mass and the surrounding stellar bulge mass (or velocity dispersion), which we take to be an indication of co-eval growth. Such growth occurs in part during the dissipative baryon accretion by BH “seeds” in the AGN (Active Galactic Nuclei) phase, but there is as yet no generally accepted scenario for the origin of the seeds. Moreover, recent suggestions of very early very supermassive BHs (e.g. Kurk et al. 2007), together with changes in the normalization of the BH mass-bulge mass proportionality (relatively larger black holes at high red shift) (e.g. Maiolino et al. 2007), suggest an alternate early growth mechanism. Recently Peirani & de Freitas Pacheo (2008) have studied the possible size of the dark matter component in BH masses. They deduce that between 1% and 10% of the black hole mass could be due to dark matter.

We explored this possibility in the case of spherically symmetric radial infall (Paper I) and in this paper we will extend it to anisotropic infall. We use the same technique of inferring reasonable distribution functions for collisionless matter from limits of the time dependent Collisionless Boltzmann (CBE) and Poisson set, and are aided by some high resolution simulations of the evolution without a dominant central mass. We add a dominant central mass either analytically or by iteration. Just as in the radial case, the loss cones are not empty for substantial growth. The relaxation of collisionless matter is due to the temporal evolution (including the radial orbit instability) and, in addition, to possible “clump-clump” (two clump) interactions (Henriksen 2009; MacMillan & Henriksen 2003).

We use, in this paper as in the previous Paper (I), the Carter-Henriksen (Carter & Henriksen 1991) procedure to obtain a quasi-self-similar system of coordinates (Henriksen 2006a,b, hereafter H2006a, H2006b). This allows writing the CBE-Poisson set with explicit reference to possible transient self-similar dynamical relaxation. In this way we can remain “close” to self-similarity just as the numerical simulations appear to do.

Studies of BH-density cusps originated with the problem of BH feeding (Peebles 1972; Bahcall & Wolf 1976), and with the notion of adiabatic growth (Young 1980; Quinlan et al. 1995; MacMillan & Henriksen 2002). Observations of the central Milky Way have detected a mainly isotropic density cusp with logarithmic power in the range  −1.1  ±  0.3 (Gillessen et al. 2009). This is flatter than the adiabatic limit (MacMillan & Henriksen 2002) and in any case the adiabatic growth scenario does not produce the BH bulge correlations (ibid). All this has spurred the investigation of co-eval dynamical growth instead of adiabatic growth. Central cusps flatter than  −1.5 can be created by tight binary BH systems formed in mergers that “scour” the stellar environment (e.g. Merritt & Szell 2006; Nakano & Makino 1999). This process can produce log slopes as flat as  −1 or even  −0.5, and it is supported by a strong correlation between nuclear BH mass and central luminosity deficit (Kormendy & Bender 2009).

However the correlation in itself only implicates the influence of the BH. It does not necessarily require the merger history, which in any case is unlikely to be the same for different galaxies. Consequently we explore in this paper whether cusps as flat as those resulting from scouring might also be produced during the dynamical formation of the black hole.

We begin the next section with a summary of the general formulation in spherical symmetry. Subsequently we study a system comprised of anisotropic non-radial orbits in spherical symmetry. Finally we give our conclusions.

2. Dynamical equations and previous results

We will use the formulation of H2006a, wherein we transform to infall variables the collisionless Boltzmann and Poisson equations for a spherically symmetric anisotropic system. We begin with the “Fujiwara” form (e.g. Fujiwara 1983), namely



Here f is the phase-space mass density, Φ is the “mean” field gravitational potential, j2 is the square of the specific angular momentum and other notation is more or less standard.

The transformation to infall variables takes the form (e.g. H2006a)


The passage to the self-similar limit requires taking T = 0 when acting on the transformed variables. Thus the self-similar limit is a stationary system in these variables, which is a state that we refer to as “self-similar virialisation” (Henriksen & Widrow 1999, hereafter HW 1999; Le Delliou 2001). The virial ratio 2K/ | W |  is a constant in this state (although greater than one; K is kinetic energy and W is potential), but the system is not steady in physical variables as infall continues.

The single quantity a is the constant that determines the dynamical similarity, called the self-similar index. It is composed of two separate reciprocal scalings, α in time and δ in space, in the form a ≡ α/δ. As it varies it reflects all dominant physical constants with dimensions of mass as well as of length and time. This is because the mass reciprocal scaling μ has been reduced to 3δ −2α in order to maintain Newton’s constant G scale invariant (e.g. H2006a).

We assume that time, radius, velocity and density are measured in fiducial units ro/vo, ro, vo and ρo respectively. The unit of the DF is fo and that of the potential is . We remove constants from the transformed equations by taking


These transformations convert Eqs. (1), (2) to the respective forms




This integro-differential system is closed by


This completes the formalism that we will use to obtain the results below. The “cusps” we describe there will generally end in what is the central “bulge” surrounding the black hole, rather than in the black hole itself.

3. Spherically symmetric steady anisotropic bulges and cusps

We follow the same pattern of discussion in this section that we used previously (Paper I) for radial orbits. We begin in this section with steady bulge-black hole systems that retain a memory of the prior nearly self-similar relaxation. The relaxation proceeds by way of the relation


(with the appropriate total energy and potential energy). This includes clump-clump interactions and the radial orbit instability and produces finally the coarse grained (i.e. steady H2006a) system.

The steady self-similar cusps were found in (Henriksen & Widrow 1995, hereafter HW95) except for the case where a = 1. However we can recover them here by using the same procedure that we used for radial orbits (see Paper I). We employ the characteristics of (5), together with the identity . A combination of the characteristics leads to an expression for the scaled energy on a characteristic, namely




We impose the steady state in Eq. (9) by setting Ψ/∂s = 0 and requiring as in Paper I that Ψ ∝ Rp with p = 2(1 − a) so that


Consequently the characteristics of Eq. (5) are seen to imply


These combine to give the general self-similar steady DF in the form (a ≠ 1 and a ≠ 3/2)




The quantity Zo is the value of Z at s = αT = 0, that is j2. However in the ab-initio steady self-similar analysis, this dependence does not appear (HW95). It is however present in general in this limit from time dependent phase, and we use it in the next paper in this series (Paper III), on non-self-similar cusps.

The physical form of the self-similar DF becomes, on using the scaling relations (3),


where . This result has also been shown (e.g. H2006a) to be the zeroth order DF in the coarse graining approach that becomes the exact steady state in the long time limit (α → ∞). It was also used by Kulessa & Lynden-Bell (Kulessa & Lynden-Bell 1992) in their study of the mass of the Galaxy.

An earlier important paper by Stiavelli & Bertin (Stiavelli & Bertin 1985) also proposed equilibrium distribution functions for elliptical galaxies. These were based in general on three integral models, but were reduced to our two familiar integrals in the case of spherical symmetry. In our case the motivation for the equilibrium DF is different, because it derives from the self-similar dynamical growth. There are nevertheless some significant similarities, particularly as modified to an inverted energy distribution by Merritt et al. (1989, hereafter MTJ). We discuss these points further below and in the conclusions to this paper.

The density profile that accompanies this DF for a particular choice of  can be shown by direct integration over phase space to be ρ ∝ r − 2a, and subsequently the consistent potential is Φ ∝ r2(1 − a) so long as a ≠ 1,3/2.

The case a = 3/2 is excluded simply because it represents a point mass in a massless halo. However this is of interest precisely when describing the environment of a dominant central mass, whether this be the halo of a central black hole or the halo of a coarse-grained collisionless bulge of stars.

The general form of the DF becomes from Eq. (15)


This gives a number density of massless particles  ∝ r-3 for a potential Φ =  − M * /r. The rms radial velocity (the same as the radial velocity dispersion since ) will generally be  ∝ Φ for any power law choice of .

One can come close to imitating the Bahcall & Wolf (Bahcall & Wolf 1976) zero flux solution for a black hole cusp by choosing . This yields the DF


where F is an arbitrary function. In fact, they use for argument of their arbitrary function λ = (1 − e2)/2, where e is the orbital eccentricity, and for Keplerian orbits around a point mass M this becomes  | E | j2/(GM)2. In our units GM = 1.

If this were exactly the Bahcall-Wolf solution it would be remarkable, since their solution is the result of collisional diffusion into the loss cone. Here however we are obliged to have the additional dependence (j2) −5/4, which implies a divergent loss cone. Thus the necessary diffusion is simply assumed in this example in order to duplicate the energy and eccentricity dependence (arbitrary because of spherical symmetry).

One can obtain the Bahcall-Wolf energy dependence from the general form (15) only by setting a = 7/3, which is quite unrealistic and gives the wrong arbitrary dependence. We shall see in Paper III that it is possible to choose a simple non-self-similar DF that does fall between these two cases and imitates the Bahcall-Wolf result more precisely.

There are two simple limits of the form (18) for which there is some evidence in the numerical simulations of isolated dark matter bulges (MacMillan 2006). Just as in the radial case the following comparison with simulations suggest that this family of DF’s may actually be realized. These limits are found by taking first  to be a constant K, and then by taking it to be proportional to κ − q. This yields respectively



where we have defined w = (3 − a)/(4 − 2a).

We expect (see discussion to follow) that the DF (18) will apply near the centre of a system where a < 1. The self similar potential is then increasing with radius and so must be taken positive to obtain an attractive force. This allows us to calculate the corresponding density explicitly as


Here B(x,y) is the standard Beta function, which may be expressed in terms of the gamma function Γ(x) as Γ(x)Γ(y)/Γ(x + y).

The Poisson equation shows that any power law solution for the potential with this density goes as Φ ∝ r2(1 − a) and hence also ρ ∝ r−2a, both as they must for self-similarity. The only reason for giving this explicit form here is that it may be used in an iterative calculation of the transition from halo to black hole dominance. In such a calculation a point mass potential plus the halo potential is inserted to give a new density, which in turn may be used in the Poisson equation to give a new potential and so on. The halo potential must dominate however so as to keep the net potential positive.

As an example we assume that the potential has the form Φ =  − C1/r + C2r2(1 − a) (C1 and C2 arbitrary positive constants) so that by (20) the density becomes ρ ∝ (C2 − C1/r3 − 2a)r − 2a. From the Poisson equation the new potential now follows as


where C3 should be positive to maintain the point mass contribution. The inner boundary of this expression is at the radius where it equals zero and the outer boundary would be inside the NFW (Navarro et al. 1996, hereafter NFW 1996) scale radius (which we may set equal to our scale radius ro). The radial mean square velocity  (also the radial squared dispersion) is proportional to Φ and so we predict a minimum value just outside the radius rm where the central mass becomes dominant (i.e., Φ(rm) = 0), followed by a steady rise that tends to r2(1 − a).

The density corresponding to this potential follows from (20) as


so that it has for a < 1 a rapid rise near rm followed by the self-similar regime wherein ρ ∝ r − 2a. The radius rm would be between 1 and 10 pc in the halos of typical galaxies.

Of course this iteration has not fully converged but subsequent cycles will produce higher order terms. The result (21) seems to indicate that the black hole mass is enhanced by an r-3 cusp that peaks just outside rm.

The radial velocity dispersion is generally equal to the root mean square radial velocity for the DF (18) since the mean radial velocity is zero. This is proportional to Φ ∝ r2(1 − a). The squared tangential velocity dispersion however becomes


where q <  − 2 or 1 > a > 1/3. The factor F(a) in this equation that multiplies  is shown in Fig. 1 as a function of a.

thumbnail Fig. 1

The curve shows the dependence of the factor that multiplies Φ in the tangential velocity dispersion as a function of a. The adiabatic variation of a with r is approximately a = r0.2 (H2007), where r is measured in units of the NFW scale radius rs.

Unlike the radial velocity dispersion, which should steadily increase with radius proportionally to Φ (which flattens as a → 1), the tangential velocity dispersion is expected to show a more precipitous drop as r tends towards the scale radius and a → 1 as F(a) declines. The variation of a is only adiabatic with radius however, varying roughly as r0.2. We shall see evidence for this below.

We turn first to an examination of the other limit (19). We expect this to apply in the outer region to which much angular momentum has been transferred by the bar due to the radial orbit instability (MacMillan et al. 2006, hereafter MWH 2006). Hence we choose the regime where a > 1 and the energy/potential is negative. Iteration is not really relevant in this limit since we are far from the black hole, but because of its importance in calculating mean quantities we give the density as


Once again the Poisson equation has the consistent solution Φ ∝  − r2(1 − a) and ρ ∝ r − 2a. The integral I0(w) is given by


where the lower limit must be greater than zero for convergence since w ≥ 1 for a ≥ 1.

The mean square radial velocity becomes in this limit


where, since a > 1, this decreases with radius. The tangential velocity dispersion has once again a more complicated expression namely


Although w(a) is slowly varying for 1 < a < 3/2, the factor multiplying Φ declines to zero as w → 3/2.

It is time however to compare these DF’s chosen for their self-similar memory to the results of simulations, in order to justify our attention.

In references (H2006a) and (Henriksen 2007, hereafter H2007) it was concluded that a ≈ 0.72 near the outer boundary of the relaxed region. Moreover we know from those papers, and more generally from extensive cosmological simulations, that a ≈ 0.5 in the interior, well relaxed, region. This is referred to as adiabatic self-similarity (H2006a) since the index a dynamically evolves, but relatively slowly (a ∝ rαF, αF = O(0.2), H2007) as already employed above. To match the simulation results we require a ≈ 3/2 in the vicinity of the NFW scale radius, after which there may be a tidal truncation to r-4 (Henriksen 2004, hereafter H2004).

There is evidence in the simulations for the persistence of self-similarity, even in the case that most closely resembles an isolated cosmological halo (i.e. cosmological perturbations are included in the initial conditions). Thus in the cosmological-like halo simulation of (MacMillan 2006), the virial mass and virial radius are found to grow in a power law manner according to t2.16 and t1.30 (we do not quote numerical errors, which are at the level of a few percent) respectively. The logarithmic density slope in the outer relaxed region is approximately  − 1.4 while the pseudo-density logarithmic slope is  − 2.16.

The numerical predictions follow by using a = 0.72 from (H2007). The self-similar mass growth inside any growing radius (fixed R) is (HW 1999)  ∝ t(3/a − 2) ∝ t2.16, while r ∝ t1/a ∝ t1.38 (see Eq. (3)). The logarithmic density slope is predicted to be  − 2a = 1.44, and the pseudo-density logarithmic power (e.g. H2007, Hansen 2004) is given by a − 3 =  − 2.28. These reasonable correspondences with the numerical simulation (only the predicted pseudo-density differs significantly from the simulated value) encourage us to adopt a DF as one of the above steady forms. The self-similar memory is incorporated either in a = 0.72 or in a = 0.5 in the relaxed region, which boundary is approximately (in fact the region where a < 1 extends beyond the scale radius by about a factor 2 in the simulation) coincident with the NFW scale radius.

The velocity dispersion in the relaxed region is based on the DF (18) and has been given in Eq. (23) for the tangential case, while the radial dispersion is proportional to . Hence, away from the transition to the black hole region, we can test our predicted dispersions against the simulation. Using the value a = 0.72 this gives the radial (σr) velocity dispersion going as r0.28 which should be compared with Fig. 2. The cosmological case is the right hand panel, which the phase space diagram (not shown) shows to be more relaxed than the initially unperturbed simulation on the left. The r0.28 rise is a reasonable description of the numerical behaviour, especially as some adiabatic evolution in a towards unity is to be expected. This will flatten the rate of rise Between 4 and 60 kpc, the radial dispersion rises by about a factor 2 while the predicted value using a = 0.72 is 2.13. Allowing a to increase only to 0.75 on average, improves the fit.

Eventually at large enough radius, according to the complete cosmological simulations, we require a → 3/2 and so we predict a decline in the radial dispersion proportional to r − 1/2 in this region. Figure 2 indicates that this begins at about twice the indicated scale radius. Before this point but still outside the scale radius the decline is slower. In fact the simulation gives a ≈ 5/4 in this region so that the expected decline would be only as r − 1/4.

The tangential dispersion velocity has the same behaviour at small r but rolls over more rapidly. This is indicated very clearly in the graphs of the anisotropy parameter on the panel below the velocity dispersion. This is consistent qualitatively with the action of the factor F(a) defined in Eq. (23), and displayed in Fig. 1 over the relevant range of a. The weak adiabatic dependence of a on r is not weak enough however to cover the range in r over which σ ⊥  stays flat. There is some uncertainty in this factor and for a ∝ r0.15 one arrives at nearly a factor two in radius as a varies from 0.72 to 0.8. However Fig. 1 shows that F has dropped by more than a factor two in this range, while  increases by only a factor of 1.2 at best. The implied factor two decline is not seen. The only explanation is that a does not vary significantly in this region, so that there is little adiabatic evolution.

At the end of this plateau we have entered the region where a ≈ 5/4 > 1. The decline in  is however only as r-0.25 there, whereas the observed decline in Fig. 2 is more like r-1. Once again the difference between this behaviour and that of the radial dispersion must lie in the factor multiplying Φ.

thumbnail Fig. 2

The figure shows two distinct simulations from (MacMillan 2006). The right hand panels display the result (with continuing infall) of an isolated halo simulation starting from a set of particles perturbed by a cosmological spectrum of density and velocity fluctuations. The panel on the left indicates the same stage as developed from a non-perturbed halo. The upper line on the upper panel in each case is for a similar halo, but allowing only purely radial collapse. The middle line is the radial velocity dispersion while the lower line is the tangential dispersion. The lower panels display the anisotropy parameter for non-radial collapse.

So we have found above that the implications of Eqs. (18) and (19) are reasonable, but do we have any evidence concerning the DF itself?

In the outer part of the relaxed region MacMillan (ibid, and Fig. 3) finds that the DF is mainly dependent on angular momentum. This is evident from the figure by noting that numerically g(E,j2) ∝ dM/dEin the range  − 100 to  − 200 energy units (MacMillan 2006). Hence from the definition


we infer that f ∝ dlndM/dE/dj2. This is independent of energy in the outer region. Note that in terms of j2 the figure is simply stretched by a factor of 2 in the j direction.

We therefore fit the self-similar DF (19) with a = 1.25 to obtain f ∝ (j2)-1.6. This is in fact a reasonable fit to the DF in the less tightly bound region with angular momentum greater than about 500 units (Fig. 3). Such a cut-off in angular momentum was contained in the Stiavelli & Bertin DF (Stiavelli & Bertin 1985) where it was exponential. MTJ gave a heuristic justification for a cut-off in j2 and also wondered whether an exponential or a power-law cut-off was superior. They also asked whether such behaviour arose naturally in the course of the dynamical formation of galaxies. Here our answer is in favour of a power-law cut-off following self-similar infall.

An intermediate region where the logarithmic density slope varies from just above 2 to about 3 does occur in all simulations of collisionless matter. Moreover this region tends to coincide with the passage from isotropy to radial orbits (e.g. 2). According to the preceding remarks (i.e. the DF (19) with a > 1) this region can be characterized by roughly constant numbers of particles over a broad range of energy, but with numbers rapidly increasing towards a minimum in angular momentum. Some particles would have too much angular momentum to enter the central region, while those below a critical angular momentum can populate the central regions.

In the simulations this redistribution of angular momentum is attributed to the onset of the Radial Orbit Instability (ROI) in (MWH 2006) and the consequent formation of a bar. In the current spherically symmetric analysis, we are forced to patch these different regions together “by hand”, guided only by the different possibilities for a.

In the inner relaxed region we expect the DF (18) to apply with a ≈ 0.72. This gives an isotropic f ∝ E-4. The simulations (NFW 1996) suggest that a → 1/2 in the very centre of the system, which would give f ∝ E − 5/2. This is in accord with previous discussions (e.g. H2006a). The fact that our energy is positive, while the simulation energy is negative, reflects the different zero points for the potential. We use the centre of the system as the reference zero in our analysis, while the simulation uses infinity.

Numerically, g(E) ∝  | E | -2.5 in this region so that d2M/dEdj2 should vary as E-6.5 (E-5 if a = 0.5). There is very little dependence on j2 here so this prediction might be compared to the behaviour of dM/dE. This quantity is falling steeply beyond E =  − 300 (3), but the evidence is weak at the limit of resolution.

At slightly higher angular momentum there is a rising dependence on j. This can be imitated with a self-similar solution by taking in Eq. (15), which gives the DF f ∝ E-5(j2)1/5.

thumbnail Fig. 3

The figure shows a contour plot of d2M/dEdJ for the full simulation of an isolated halo (MacMillan 2006). The separate plots of dM/dE and dMdJ are also indicated. The corresponding density of states is almost independent of J and is the expected (e.g. Binney & Tremaine 1987)  | E | -2.5 law. The line on the upper left is J(E) for circular orbits.

The above prolonged discussion suggests that a DF with a self-similar memory can be used to parametrise the phase space of these collisionless simulations. At least this is true if we accept the adiabatic variability of a, which has been made plausible on the grounds of a local entropy maximum (e.g. H2007). However what has this to do with the black-hole bulge mass correlation that the co-eval growth was supposed to explain?

If one assumes that the halo around the seed black hole is accreted eventually through two-body and clump-clump relaxation (MacMillan & Henriksen 2003), then the mass of the black hole/halo should be proportional to the mass of the bulge simply because of the power law density (with nearly constant power inside rs). Letting M be the black hole mass, Ms be the bulge mass and rhrs, be the black hole halo and bulge radii respectively (rs is the scale radius); we find the relation


If the appropriate a = 0.72, then the power of the radial dependence is 1.56. For a true universal relation however, we must assume that galaxy growth is not only self-similar, but that there is a similarity between galaxies so that the scale ratio is similar. This is not exactly the case (e.g. Navarro et al. 2010), but it is nearly so. Such a relation might be tested by high resolution simulations containing a seed black hole.

3.1. More general self-similar distribution functions

The limits of Eq. (15) that we discussed above are quite arbitrary, except for the expectation of central isotropy and the dominance of angular momentum in the outer regions. In this section we explore more general possibilities, even though empirically the preceding limits seem adequate.

The distribution function (15) can be put in a form that seems to generalize the FPDF (i.e. Fridmann & Polyachenko DF). We choose to obtain


The density integral over this DF requires a lower cut-off in angular momentum for a > 1 and an upper cut-off in energy for a < 1 in order to avoid singularities. We have exhibited the upper energy cut-off in (30) explicitly. It appears as the arbitrary constant in the potential so that the difference Eo − E still has the self-similar behaviour (r2(1 − a)), as does the density (r − 2a). Moreover for a > 1 the minimum angular momentum should be a fixed fraction of the value 2r2(E − Φ), say 0 < k1 ≪ 1, in order for the density profile and potential associated with Eq. (15) to apply. The radial velocity dispersion is also  ∝ r2(1 − a). It is only when a = 1 that one obtains the r-2 profile of the radial FPDF, and this must be discussed separately below.

If one computes formally the density by integrating the DF (30) in a negative energy region (negative although a < 1, because the central mass dominates), one finds that for a < 1


As in the previous section, this can be used in an iterative fashion.

Thus in a region dominated by a central point mass (and taking Eo = 0) we see that the cusp would be, at lowest order,


The next order would be found by using this in the Poisson equation to obtain a new potential and then a new density from Eq. (31). Once again this low order result can not be flatter than r-1.5. The steepest limit is r-2 as a → 1 − . We remark that at a = 2/3 one obtains the Bahcall & Wolf zero flux cusp r − 7/4, due to the filled loss cone in the DF. For a = 0.72, the cusp is like r-1.78, very close to the Bahcall & Wolf cusp and still with a filled loss cone. In Paper III of this series we shall find an exact example of this type.

3.2. Singular isothermal sphere: the global inverse square law

So far we have ignored the case when a = 1, the singular isothermal limit. In fact the steady DF with a = 1 must be treated separately because of the logarithmic potential Φ = Ψlnr (Ψ constant) that accompanies the inverse square density law. We may still deduce it from the characteristics of (5) if we take T = 0 finally, since a steady state is the same at any time. The appropriate characteristics become, with a = 1,


These integrate to give




and Zo is again j2 at s = 0 on the characteristic.

By writing the self-similar physical form at T = 0 we find (for a = 1)



Once again so that κ = E − Ψ/2lnj2.

This limit with a = 1 is an extension of the inverse square density law to non-radial orbits, but the function of κ is arbitrary. The DF is certainly non-unique for the same density profile.

If we calculate the density profile for κ > 0 by integrating over (36) we arrive eventually at


where u ≡ j2/r2. The limits u1, u2 are the two roots found in terms of the Lambert function by setting the argument of the square root equal to zero. When κ/Ψ is large, the lower root becomes zero and the upper root approaches 2κ. There is slow convergence of the inner integral. Provided that  is such as to make the outer integral converge, the profile is r-2.

We note that for any finite κ, the possible range of u has a definite upper limit and a definite lower limit that may be close to zero. This means that at given r there is a finite range in j2 that is present, and that this range collapses on zero as r → 0. This is as expected with spherical symmetry.

It is possible to have solutions with κ < 0 whereupon


where again the inner limits are the roots found by setting the argument of the square root equal to zero. For a real range of values, Ψ/ | κ |  must be greater than a minimum value that depends on the value of Ψ. The same restriction on angular momentum as r → 0 applies. In all cases the velocity dispersion is constant.

Unlike the radial FPDF that has an inverse square cusp, a black hole is not readily incorporated into this distribution. However at large enough r (or basically where in phase space where ) it suffices to modify the integral by a point mass potential so that


The velocity dispersion is no longer strictly constant, unless j2 ≫ GMr, but the influence of the black hole is weak in this regime.

One special model of an inverse square, steady, anisotropic bulge is provided by Eq. (36) by taking


where b is any real number.

In velocity space the anisotropic DF becomes explicitly ()


This form is quite closely related to the DF suggested by Stiavelli & Bertin (Stiavelli & Bertin 1985). However it offers a combination of an exponential cut-off and a possible power-law cut-off in angular momentum. It does have the inverted distribution in energy (it increases outward) as suggested by MTJ. This is also true of the DF that we suggest for the inner region, namely Eq. (18) with a < 1. Thus we obtain this inversion naturally as a consequence of the self-similar evolution.

When b = 0 we have a pure Gaussian in the energy and once again the singular isothermal sphere. If b < 0 there is a filled loss cone and for b > 0 the loss cone is empty. Provided b >  − 1 this DF is elliptical in velocity space (vrv ⊥ ) at each r and becomes circular (isotropic) at b = 0. The isothermal form with b = 0 also appears in the fully asymmetric case of zeroth order coarse graining as discussed in the appendix of (H2004).

This example is only of interest in the context of the persistent r-2 mass distribution in galaxies and clusters referred to in Paper I. There is clearly a lack of uniqueness however, as several distribution functions may produce an inverse square density profile. Adding a black hole to this example is inconvenient (unless it is negligible) as the form of the DF is not suitable for iteration. However just as in the general case above, κ may be taken as that given in Eq. (40).

One can calculate the growth of a central black hole embedded in an inner bulge with the DF (18). We can expect to maintain the steady bulge DF until the black hole mass is a significant fraction of the bulge mass. We find after obvious but tedious calculations that


where x ≡ r/rs, the subscript s refers to characteristic bulge quantities, and a > 1/3 for convergence of the integrals. We recall that the last stable radius of a Schwarzschild black hole is r ≡ 6GM/c2. For a = 2/3 the growth is exponential, but the e-folding time is longer than the age of the Universe (5  ×  1019   s) even if the bulge has 1010   M ⊙  in one kiloparsec. The result is similar for other values of a. However growth of halo cloaking the actual black hole could be much faster. In general the mass growth from the DF (18) is (rh and Mh are replaced by the black hole values to obtain the previous equation)


where Mh is the mass accreted inside the radius rh. This radius would be an intermediate scale between the bulge scale rs and r, determined ultimately by the minimum angular momentum actually available. The accretion rate is proportional to the radius of the structure and is faster than the black hole rate by rh/r when a = 1. One assumes that all of this mass is eventually accreted by the black hole through relaxation processes (e.g. MacMillan & Henriksen 2003), although multiple steps might be required.Then the black hole growth is limited by the time to empty this reservoir onto the centre.

This calculation is done without any bias towards the loss cone (j ≈ 0), and with a self-similar DF of the form (30), one can actually grow the black hole directly. We shall save this calculation for Paper III.

4. Conclusions

In this paper, we have developed distribution functions that describe both dark matter bulges and a central black hole or at least a central mass concentration. We succeed mainly in describing the dark matter bulges.

In the discussion of cusps and bulges based on purely radial orbits (Paper I), we were able to distinguish the Distribution function of Fridmann & Polyachenkofrom that of Henriksen & Widrow. The FPDF was found to describe accurately the purely radial simulations of isolated collisionless halos carried out in (MacMillan 2006). Moreover a point mass could be added without changing the form of the DF, which allows a self-similar growth of a central bulge or black hole.

The final result concerning steady, self-similar radial orbits concerned the special case a = 1. The DF is a Gaussian that has been found previously in coarse graining (Henriksen & Le Delliou 2002). We included it there as a second example of a radial DF that produces an r-2 density profile (Mutka 2009).

The inclusion of angular momentum here leads to more realistic situations. We re-derived the steady self-similar DFs from first principles in Eq. (15). We showed that these can be used to describe the simulated collisionless halos calculated in (MacMillan 2006), if we use the value of a ≈ 0.72 identified in (H2007) and take two limits. In one (18) the DF is isotropic and describes approximately the most relaxed central region of the bulge. The other limit (19) describes the outer region and the transition across the NFW scale radius. The qualitative correspondence supports the relevancy of this family. The potential of a central mass can not be included exactly in these distribution functions, but iteration is possible. An example was given for the DF (18).

In particular, velocity dispersions have been calculated and compared qualitatively with the results of a high resolution numerical simulation of an isolated halo. These correspond to reasonable densities when adiabatic self-similarity is employed. Unfortunately these do not apply to the vicinity of the black hole itself. However the eventual dominance of the black hole through the r-1 potential becomes obvious, as seen in the iterated potential (21).

Simple forms for the DF have been suggested for the various regions of the halo. And given the dynamic co-eval growth of black-hole/halo and bulge, we explain the black-hole bulge mass correlation simply (29). This relies on the approximate similarity between different halos.

The forms that we find follow from the assumed self-similar path to equilibrium. Remarkably they are in accord with the general forms for the DF as presented in Stiavelli & Bertin (Stiavelli & Bertin 1985) and as modified in MTJ. That is, these distribution functions contain an inner inverted distribution in energy (referred to as negative temperature in MTJ) and a cut-off in the distribution in angular momentum farther from the centre. The scale of the inner region that is argued to coincide with NFW scale radius in this paper is the parameter ra in the earlier papers, so inner and outer have similar meanings.

The variation is a power-law generally in both energy and angular momentum, although in the very relevant case of an inverse square density profile, an exponential variation is possible. These conclusions follow from the assumed adiabatically self-similar evolution towards equilibrium. As such they tend to answer questions posed in MTJ such as whether these distributions arise naturally and if the variations are exponential or power-law. We now know however that part of the answer lies in the anisotropy produced by the radial orbit instability (MWH 2006).

We have also found a self-similar generalization of the radial FPDF in Eq. (30). Unfortunately this DF does not have the property of yielding a density that is independent of the potential and so a point mass potential is incompatible with self-similarity. It might describe a system at large radii with a large central mass.

As always a = 1 must be treated separately and we give a derivation from first principles. The result is new and it generalizes the FPDF in the sense that ρ ∝ r-2 always. Although we can not place a central mass inside this bulge exactly, it allows an anisotropic DF in the r-2 bulge region.

ONK7L 3N6, Generally we find that we can not describe elegantly anisotropic bulges containing black holes with self-similar DFs.

In the next paper in this series (Paper III), we shall widen our scope to the study and production of non-self-similar cusps and DF.


R.N.H. acknowledges the support of an operating grant from the canadian Natural Sciences and Research Council. The work of MLeD is supported by CSIC (Spain) under the contract JAEDoc072, with partial support from CICYT project FPA2006-05807, at the IFT, Universidad Autónoma de Madrid, Spain.


  1. Bahcall, J., & Wolf, R. A. 1976, ApJ, 209, 214 [NASA ADS] [CrossRef] [Google Scholar]
  2. Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, New Jersey: Princeton University Press) [Google Scholar]
  3. Carter, B., & Henriksen, R. N. 1991, JMP, 32, 2580 [NASA ADS] [CrossRef] [Google Scholar]
  4. Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9 [NASA ADS] [CrossRef] [Google Scholar]
  5. Fridman, A. M., & Polyachenko, V. L. 1984, Physics of Gravitating Systems (New York: Springer) [Google Scholar]
  6. Fujiwara, T. 1983, PASJ, 35, 547 [NASA ADS] [Google Scholar]
  7. Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13 [NASA ADS] [CrossRef] [Google Scholar]
  8. Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 1075 [NASA ADS] [CrossRef] [Google Scholar]
  9. Hansen, S. 2004, MNRAS, 352, L41 [NASA ADS] [CrossRef] [Google Scholar]
  10. Henriksen, R. N. 2004, MNRAS, 355, 1217 [NASA ADS] [CrossRef] [Google Scholar]
  11. Henriksen, R. N. 2006a, ApJ, 653, 894 [NASA ADS] [CrossRef] [Google Scholar]
  12. Henriksen, R. N. 2006b, MNRAS, 366, 697 [NASA ADS] [CrossRef] [Google Scholar]
  13. Henriksen, R. N. 2007, ApJ, 671, 1147 [NASA ADS] [CrossRef] [Google Scholar]
  14. Henriksen, R. N. 2009, ApJ, 690, 102 [NASA ADS] [CrossRef] [Google Scholar]
  15. Henriksen, R. N., & Widrow, L. M. 1995, MNRAS, 276, 679 [NASA ADS] [CrossRef] [Google Scholar]
  16. Henriksen, R. N., & Widrow, L. M. 1999, MNRAS, 302, 321 [NASA ADS] [CrossRef] [Google Scholar]
  17. Henriksen, R. N., & Le Delliou, M. 2002, MNRAS, 331, 423 [NASA ADS] [CrossRef] [Google Scholar]
  18. Kulessa, A. S., & Lynden-Bell, D. 1992, MNRAS, 255, 105 [NASA ADS] [Google Scholar]
  19. Kurk, J. D., Walter, F., Fan, X., et al. 2007, ApJ, 669, 32 [NASA ADS] [CrossRef] [Google Scholar]
  20. Kormendy, J., & Bender, R. 2009, ApJ, 691, L142 [NASA ADS] [CrossRef] [Google Scholar]
  21. Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581 [NASA ADS] [CrossRef] [Google Scholar]
  22. Le Delliou, M. 2001, Ph.D. Thesis, Queen’s University, Kingston, Canada [Google Scholar]
  23. Le Delliou, M., Henriksen, R. N., & MacMillan, J. D. 2009, submitted [arXiv:0911.2238] (Paper III) [Google Scholar]
  24. Le Delliou, M., Henriksen, R. N., & MacMillan, J. D. 2010 [arXiv:0911.2232] (Paper I) [Google Scholar]
  25. MacMillan, J. 2006, Ph.D. Thesis, Queen’s University at Kingston, Canada [Google Scholar]
  26. MacMillan, J. D., & Henriksen, R. N. 2002, ApJ, 569, 83 [NASA ADS] [CrossRef] [Google Scholar]
  27. MacMillan, J. D., & Henriksen, R. N. 2003, Carnegie Observatories Astrophysics Series, 1, Co-evolution of Black Holes and Galaxies, L.C. Ho [Google Scholar]
  28. MacMillan, J. D., Widrow, L. M., & Henriksen, R. N. 2006, ApJ, 653, 43 [NASA ADS] [CrossRef] [Google Scholar]
  29. Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285 [NASA ADS] [CrossRef] [Google Scholar]
  30. Maiolino, R., Neri, R., Beelen, A., et al. 2007, A&A, 472, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Merritt, D., & Szell, A. 2006, ApJ, 648, 890 [NASA ADS] [CrossRef] [Google Scholar]
  32. Merritt, D., Tremaine, S., & Johnstone, D. 1989, MNRAS, 236, 829 [NASA ADS] [Google Scholar]
  33. Mutka, P. 2009, Proceeding of Invisible Universe, Palais de l’UNESCO, Paris, ed. J.-M. Alimi [Google Scholar]
  34. Nakano, T., & Makino, M. 1999, ApJ, 525, L77 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  35. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [NASA ADS] [CrossRef] [Google Scholar]
  36. Navarro, J. F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 402, 21 [NASA ADS] [CrossRef] [Google Scholar]
  37. Peirani, S., & de Freitas Pacheo, J. A. 2008, Phys. Rev. D, 77 (6), 064023 [Google Scholar]
  38. Peebles, P. J. E. 1972, Gen. Rel. Grav., 3, 63 [NASA ADS] [CrossRef] [Google Scholar]
  39. Quinlan, G. D., Hernquist, L., & Sigurdsson, S. 1995, ApJ, 440, 554 [NASA ADS] [CrossRef] [Google Scholar]
  40. Stiavelli, M., & Bertin, G. 1985, MNRAS, 217, 735 [NASA ADS] [Google Scholar]
  41. Young, P. 1980, ApJ, 242, 1232 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail Fig. 1

The curve shows the dependence of the factor that multiplies Φ in the tangential velocity dispersion as a function of a. The adiabatic variation of a with r is approximately a = r0.2 (H2007), where r is measured in units of the NFW scale radius rs.

In the text
thumbnail Fig. 2

The figure shows two distinct simulations from (MacMillan 2006). The right hand panels display the result (with continuing infall) of an isolated halo simulation starting from a set of particles perturbed by a cosmological spectrum of density and velocity fluctuations. The panel on the left indicates the same stage as developed from a non-perturbed halo. The upper line on the upper panel in each case is for a similar halo, but allowing only purely radial collapse. The middle line is the radial velocity dispersion while the lower line is the tangential dispersion. The lower panels display the anisotropy parameter for non-radial collapse.

In the text
thumbnail Fig. 3

The figure shows a contour plot of d2M/dEdJ for the full simulation of an isolated halo (MacMillan 2006). The separate plots of dM/dE and dMdJ are also indicated. The corresponding density of states is almost independent of J and is the expected (e.g. Binney & Tremaine 1987)  | E | -2.5 law. The line on the upper left is J(E) for circular orbits.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.