Issue 
A&A
Volume 522, November 2010



Article Number  A28  
Number of page(s)  9  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/200913648  
Published online  28 October 2010 
Black holes and galactic density cusps
Spherically symmetric anisotropic cusps
^{1}
Instituto de Física Teórica UAM/CSIC, Facultad de Ciencias, CXI,
Universidad Autónoma de Madrid,
Cantoblanco,
28049
Madrid,
Spain
email: Morgan.LeDelliou@uam.es
^{2}
Queen’s University,
Kingston, Ontario, Canada
email: henriksn@astro.queensu.ca
^{3}
Faculty of Science, University of Ontario Institute of Technology,
Oshawa,
Ontario, L1H 7K4, Canada
email: joseph.macmillan@gmail.com
Received:
11
November
2009
Accepted:
20
February
2010
Aims. In this paper we study density cusps that may contain central black holes. The actual coeval selfsimilar growth would not distinguish between the central object and the surroundings.
Methods. To study the environment of a growing black hole we seek descriptions of steady “cusps” that may contain a black hole and that retain at least a memory of selfsimilarity. We refer to the environment in brief as the “bulge” and on smaller scales, the “halo”.
Results. We find simple descriptions of the simulations of collisionless matter by comparing predicted densities, velocity dispersions and distribution functions with the simulations. In some cases central point masses may be included by iteration. We emphasize that the coeval selfsimilar growth allows an explanation of the black hole bulge mass correlation between approximately similar collisionless systems.
Conclusions. We have derived our results from first principles assuming adiabatic selfsimilarity and either selfsimilar virialisation or normal steady virialisation. We conclude that distribution functions that retain a memory of selfsimilar evolution provide an understanding of collisionless systems. The implied energy relaxation of the collisionless matter is due to the time dependence. Phase mixing relaxation may be enhanced by clumpclump interactions.
Key words: cosmology: theory / dark matter / galaxies: halos / galaxies: nuclei / black hole physics / gravitation
© 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 coeval 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 coeval 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 massbulge 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 “clumpclump” (two clump) interactions (Henriksen 2009; MacMillan & Henriksen 2003).
We use, in this paper as in the previous Paper (I), the CarterHenriksen (Carter & Henriksen 1991) procedure to obtain a quasiselfsimilar system of coordinates (Henriksen 2006a,b, hereafter H2006a, H2006b). This allows writing the CBEPoisson set with explicit reference to possible transient selfsimilar dynamical relaxation. In this way we can remain “close” to selfsimilarity just as the numerical simulations appear to do.
Studies of BHdensity 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 coeval 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 nonradial 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 phasespace mass density, Φ is the “mean” field gravitational potential, j^{2} 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 selfsimilar limit requires taking ∂_{T} = 0 when acting on the transformed variables. Thus the selfsimilar limit is a stationary system in these variables, which is a state that we refer to as “selfsimilar 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 selfsimilar 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 r_{o}/v_{o}, r_{o}, v_{o} and ρ_{o} respectively. The unit of the DF is f_{o} 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
and
This integrodifferential 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 bulgeblack hole systems that retain a memory of the prior nearly selfsimilar relaxation. The relaxation proceeds by way of the relation
(with the appropriate total energy and potential energy). This includes clumpclump interactions and the radial orbit instability and produces finally the coarse grained (i.e. steady H2006a) system.
The steady selfsimilar 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
where
We impose the steady state in Eq. (9) by setting ∂Ψ/∂s = 0 and requiring as in Paper I that Ψ ∝ R^{p} with p = 2(1 − a) so that
Consequently the characteristics of Eq. (5) are seen to imply
These combine to give the general selfsimilar steady DF in the form (a ≠ 1 and a ≠ 3/2)
where
The quantity Z_{o} is the value of Z at s = αT = 0, that is j^{2}. However in the abinitio steady selfsimilar 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 nonselfsimilar cusps.
The physical form of the selfsimilar 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 & LyndenBell (Kulessa & LyndenBell 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 selfsimilar 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 Φ ∝ r^{2(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 coarsegrained 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 − e^{2})/2, where e is the orbital eccentricity, and for Keplerian orbits around a point mass M_{•} this becomes  E  j^{2}/(GM_{•})^{2}. In our units GM_{•} = 1.
If this were exactly the BahcallWolf 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 (j^{2})^{ −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 BahcallWolf 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 nonselfsimilar DF that does fall between these two cases and imitates the BahcallWolf 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 Φ ∝ r^{2(1 − a)} and hence also ρ ∝ r^{−2a}, both as they must for selfsimilarity. 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 Φ = − C_{1}/r + C_{2}r^{2(1 − a)} (C_{1} and C_{2} arbitrary positive constants) so that by (20) the density becomes ρ ∝ (C_{2} − C_{1}/r^{3 − 2a})r^{ − 2a}. From the Poisson equation the new potential now follows as
where C_{3} 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 r_{o}). The radial mean square velocity (also the radial squared dispersion) is proportional to Φ and so we predict a minimum value just outside the radius r_{m} where the central mass becomes dominant (i.e., Φ(r_{m}) = 0), followed by a steady rise that tends to r^{2(1 − a)}.
The density corresponding to this potential follows from (20) as
so that it has for a < 1 a rapid rise near r_{m} followed by the selfsimilar regime wherein ρ ∝ r^{ − 2a}. The radius r_{m} 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 r_{m}.
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 Φ ∝ r^{2(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.
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 = r^{0.2} (H2007), where r is measured in units of the NFW scale radius r_{s}. 
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 r^{0.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 Φ ∝ − r^{2(1 − a)} and ρ ∝ r^{ − 2a}. The integral I_{0}(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 selfsimilar 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 selfsimilarity (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 selfsimilarity, 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 cosmologicallike halo simulation of (MacMillan 2006), the virial mass and virial radius are found to grow in a power law manner according to t^{2.16} and t^{1.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 pseudodensity logarithmic slope is − 2.16.
The numerical predictions follow by using a = 0.72 from (H2007). The selfsimilar mass growth inside any growing radius (fixed R) is (HW 1999) ∝ t^{(3/a − 2)} ∝ t^{2.16}, while r ∝ t^{1/a} ∝ t^{1.38} (see Eq. (3)). The logarithmic density slope is predicted to be − 2a = 1.44, and the pseudodensity logarithmic power (e.g. H2007, Hansen 2004) is given by a − 3 = − 2.28. These reasonable correspondences with the numerical simulation (only the predicted pseudodensity differs significantly from the simulated value) encourage us to adopt a DF as one of the above steady forms. The selfsimilar 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 r^{0.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 r^{0.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 ∝ r^{0.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 Φ.
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 nonperturbed 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 nonradial 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,j^{2}) ∝ dM/dEin the range − 100 to − 200 energy units (MacMillan 2006). Hence from the definition
we infer that f ∝ dlndM/dE/dj^{2}. This is independent of energy in the outer region. Note that in terms of j^{2} the figure is simply stretched by a factor of 2 in the j direction.
We therefore fit the selfsimilar DF (19) with a = 1.25 to obtain f ∝ (j^{2})^{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 cutoff in angular momentum was contained in the Stiavelli & Bertin DF (Stiavelli & Bertin 1985) where it was exponential. MTJ gave a heuristic justification for a cutoff in j^{2} and also wondered whether an exponential or a powerlaw cutoff 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 powerlaw cutoff following selfsimilar 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 d^{2}M/dEdj^{2} should vary as E^{6.5} (E^{5} if a = 0.5). There is very little dependence on j^{2} 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 selfsimilar solution by taking in Eq. (15), which gives the DF f ∝ E^{5}(j^{2})^{1/5}.
Fig. 3
The figure shows a contour plot of d^{2}M/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 selfsimilar 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 blackhole bulge mass correlation that the coeval growth was supposed to explain?
If one assumes that the halo around the seed black hole is accreted eventually through twobody and clumpclump 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 r_{s}). Letting M_{•} be the black hole mass, M_{s} be the bulge mass and r_{h}, r_{s}, be the black hole halo and bulge radii respectively (r_{s} 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 selfsimilar, 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 selfsimilar 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 cutoff in angular momentum for a > 1 and an upper cutoff in energy for a < 1 in order to avoid singularities. We have exhibited the upper energy cutoff in (30) explicitly. It appears as the arbitrary constant in the potential so that the difference E_{o} − E still has the selfsimilar behaviour (r^{2(1 − a)}), as does the density (r^{ − 2a}). Moreover for a > 1 the minimum angular momentum should be a fixed fraction of the value 2r^{2}(E − Φ), say 0 < k_{1} ≪ 1, in order for the density profile and potential associated with Eq. (15) to apply. The radial velocity dispersion is also ∝ r^{2(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 E_{o} = 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
where
and Z_{o} is again j^{2} at s = 0 on the characteristic.
By writing the selfsimilar physical form at T = 0 we find (for a = 1)
Once again so that κ = E − Ψ/2lnj^{2}.
This limit with a = 1 is an extension of the inverse square density law to nonradial orbits, but the function of κ is arbitrary. The DF is certainly nonunique for the same density profile.
If we calculate the density profile for κ > 0 by integrating over (36) we arrive eventually at
where u ≡ j^{2}/r^{2}. The limits u_{1}, u_{2} 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 j^{2} 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 j^{2} ≫ GM_{•}r, 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 cutoff and a possible powerlaw cutoff 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 selfsimilar 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 (v_{r}, v_{ ⊥ }) 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_{•}/r_{s}, 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_{•}/c^{2}. For a = 2/3 the growth is exponential, but the efolding time is longer than the age of the Universe (5 × 10^{19} s) even if the bulge has 10^{10} 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 (r_{h} and M_{h} are replaced by the black hole values to obtain the previous equation)
where M_{h} is the mass accreted inside the radius r_{h}. This radius would be an intermediate scale between the bulge scale r_{s} 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 r_{h}/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 selfsimilar 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 selfsimilar growth of a central bulge or black hole.
The final result concerning steady, selfsimilar 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 rederived the steady selfsimilar 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 selfsimilarity 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 coeval growth of blackhole/halo and bulge, we explain the blackhole bulge mass correlation simply (29). This relies on the approximate similarity between different halos.
The forms that we find follow from the assumed selfsimilar 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 cutoff 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 r_{a} in the earlier papers, so inner and outer have similar meanings.
The variation is a powerlaw 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 selfsimilar 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 powerlaw. 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 selfsimilar 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 selfsimilarity. 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 selfsimilar DFs.
In the next paper in this series (Paper III), we shall widen our scope to the study and production of nonselfsimilar cusps and DF.
Acknowledgments
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 FPA200605807, at the IFT, Universidad Autónoma de Madrid, Spain.
References
 Bahcall, J., & Wolf, R. A. 1976, ApJ, 209, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, New Jersey: Princeton University Press) [Google Scholar]
 Carter, B., & Henriksen, R. N. 1991, JMP, 32, 2580 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Fridman, A. M., & Polyachenko, V. L. 1984, Physics of Gravitating Systems (New York: Springer) [Google Scholar]
 Fujiwara, T. 1983, PASJ, 35, 547 [NASA ADS] [Google Scholar]
 Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13 [NASA ADS] [CrossRef] [Google Scholar]
 Gillessen, S., Eisenhauer, F., Trippe, S., et al. 2009, ApJ, 692, 1075 [NASA ADS] [CrossRef] [Google Scholar]
 Hansen, S. 2004, MNRAS, 352, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Henriksen, R. N. 2004, MNRAS, 355, 1217 [NASA ADS] [CrossRef] [Google Scholar]
 Henriksen, R. N. 2006a, ApJ, 653, 894 [NASA ADS] [CrossRef] [Google Scholar]
 Henriksen, R. N. 2006b, MNRAS, 366, 697 [NASA ADS] [CrossRef] [Google Scholar]
 Henriksen, R. N. 2007, ApJ, 671, 1147 [NASA ADS] [CrossRef] [Google Scholar]
 Henriksen, R. N. 2009, ApJ, 690, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Henriksen, R. N., & Widrow, L. M. 1995, MNRAS, 276, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Henriksen, R. N., & Widrow, L. M. 1999, MNRAS, 302, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Henriksen, R. N., & Le Delliou, M. 2002, MNRAS, 331, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Kulessa, A. S., & LyndenBell, D. 1992, MNRAS, 255, 105 [NASA ADS] [Google Scholar]
 Kurk, J. D., Walter, F., Fan, X., et al. 2007, ApJ, 669, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Kormendy, J., & Bender, R. 2009, ApJ, 691, L142 [NASA ADS] [CrossRef] [Google Scholar]
 Kormendy, J., & Richstone, D. 1995, ARA&A, 33, 581 [NASA ADS] [CrossRef] [Google Scholar]
 Le Delliou, M. 2001, Ph.D. Thesis, Queen’s University, Kingston, Canada [Google Scholar]
 Le Delliou, M., Henriksen, R. N., & MacMillan, J. D. 2009, submitted [arXiv:0911.2238] (Paper III) [Google Scholar]
 Le Delliou, M., Henriksen, R. N., & MacMillan, J. D. 2010 [arXiv:0911.2232] (Paper I) [Google Scholar]
 MacMillan, J. 2006, Ph.D. Thesis, Queen’s University at Kingston, Canada [Google Scholar]
 MacMillan, J. D., & Henriksen, R. N. 2002, ApJ, 569, 83 [NASA ADS] [CrossRef] [Google Scholar]
 MacMillan, J. D., & Henriksen, R. N. 2003, Carnegie Observatories Astrophysics Series, 1, Coevolution of Black Holes and Galaxies, L.C. Ho [Google Scholar]
 MacMillan, J. D., Widrow, L. M., & Henriksen, R. N. 2006, ApJ, 653, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285 [NASA ADS] [CrossRef] [Google Scholar]
 Maiolino, R., Neri, R., Beelen, A., et al. 2007, A&A, 472, L33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Merritt, D., & Szell, A. 2006, ApJ, 648, 890 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, D., Tremaine, S., & Johnstone, D. 1989, MNRAS, 236, 829 [NASA ADS] [Google Scholar]
 Mutka, P. 2009, Proceeding of Invisible Universe, Palais de l’UNESCO, Paris, ed. J.M. Alimi [Google Scholar]
 Nakano, T., & Makino, M. 1999, ApJ, 525, L77 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Ludlow, A., Springel, V., et al. 2010, MNRAS, 402, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Peirani, S., & de Freitas Pacheo, J. A. 2008, Phys. Rev. D, 77 (6), 064023 [Google Scholar]
 Peebles, P. J. E. 1972, Gen. Rel. Grav., 3, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Quinlan, G. D., Hernquist, L., & Sigurdsson, S. 1995, ApJ, 440, 554 [NASA ADS] [CrossRef] [Google Scholar]
 Stiavelli, M., & Bertin, G. 1985, MNRAS, 217, 735 [NASA ADS] [Google Scholar]
 Young, P. 1980, ApJ, 242, 1232 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
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 = r^{0.2} (H2007), where r is measured in units of the NFW scale radius r_{s}. 

In the text 
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 nonperturbed 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 nonradial collapse. 

In the text 
Fig. 3
The figure shows a contour plot of d^{2}M/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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.