text:
# Chaotic Loss Cones, Black Hole Fueling and the 𝑀_∙-𝜎 Relation

## 1 Introduction

Fueling of active galactic nuclei (AGNs) or quasars requires accretion at rates of $`\dot{M}1.5ϵ_{0.1}^1L_{46}_{}\mathrm{yr}^1`$ onto the central supermassive black hole (SBH), where $`L_{46}`$ is the energy output in units of $`10^{46}`$ ergs s<sup>-1</sup> and $`ϵ_{0.1}`$ is the mass conversion efficiency in units of its canonical value $`0.1`$. The fueling problem is usually broken into three parts: what is the fuel; how is it channeled into the SBH; and how does it radiate a substantial fraction of its energy before disappearing down the hole? This paper addresses the first two questions. Galactic spheroids have ample supplies of both stars and gas, but it is difficult to come up with mechanisms that can extract almost all of a mass element’s orbital angular momentum in a few crossing times, as required if the matter is to find its way into the event horizon of the black hole. Gas release through tidal disruption of stars is a possible fueling mechanism (Hills, 1975; Frank & Rees, 1976), but the rate at which stars are scattered onto low angular momentum orbits by gravitational encounters is too low to reproduce observed luminosities (Young, Shields & Wheeler 1977; Frank 1978). The fact that most AGNs are in spiral galaxies suggests that interstellar gas provides the bulk of the fuel. Gas can be driven into the SBH by torques from non-axisymmetric potential perturbations, due to stellar bars or to transient distortions of the potential during mergers or accretion events (Shlosman, Begelman & Frank 1990).

As a number of authors have pointed out, stellar feeding rates are also enhanced by the presence of non-axisymmetric perturbations (Norman & Silk, 1983; Gerhard & Binney, 1985). A barlike or triaxial potential is populated mainly by box orbits, and a star on a box orbit passes near the center once per crossing time. If a SBH is added, most of these “centrophilic” orbits become chaotic due to large-angle deflections by the SBH (Valluri & Merritt, 1998). Estimates of feeding rates due to centrophilic orbits in triaxial potentials are several times larger than rates due to scattering onto eccentric orbits in the spherical or axisymmetric geometries (Gerhard & Binney, 1985). While such an enhancement is significant, it is still not enough to explain the high luminosities of AGNs and quasars.

A number of factors motivated us to re-open the question of stellar fueling rates in non-axisymmetric nuclei. The $`M_{}\sigma `$ relation demonstrates a tight link between SBH masses and the kinematics of their stellar spheroids. While the origin of the relation is still uncertain, some scenarios postulate a role for stellar feeding (e. g. Zhao, Haehnelt & Rees 2002). Capture of stellar-mass objects by SBHs may provide an important source of signals for gravitational wave detectors like LISA (Hughes et al., 2001). Imaging of the centers of galaxies on scales of $`10`$ pc reveals a wealth of features in the stellar distribution that are not consistent with axisymmetry, including bars, bars-within-bars, nuclear spirals, and other misaligned structures (Wozniak et al., 1995; Rest et al., 2001; Peng et al., 2002; Erwin & Sparke, 2002). The steep power-law dependence of stellar density on radius near the centers of many galaxies revealed by HST (Crane et al., 1993; Gebhardt et al., 1996) implies higher stellar feeding rates than in older models (Norman & Silk, 1983; Gerhard & Binney, 1985) which postulated constant-density cores. Triaxiality has recently been shown to be sustainable in numerical models for black-hole nuclei (Poon & Merritt 2001, 2002, 2003, hereafter Papers I-III). Furthermore such models can contain a large population of centrophilic – typically chaotic – orbits, as high as $`75\%`$.

This paper, the fourth in a series on the dynamics of triaxial black-hole nuclei, examines the behavior of centrophilic orbits and the implications for black hole feeding. Test-particle integrations in fixed potentials are combined with knowledge of the orbital population in the self-consistent models to infer the rate at which stars would be supplied to the central SBH. We confirm and extend the results of earlier authors, who investigated stellar feeding rates in triaxial nuclei with constant-density cores. We find that stellar capture rates in steep power-law nuclei can be orders of magnitude greater than in nuclei with cores; in fact, feeding rates can meet or even exceed the so-called “full loss cone” rate, i.e. the capture rate in a spherical or axisymmetric galaxy in which the loss cone is continuously repopulated. The accretion rate in a dense, $`\rho r^2`$ nucleus scales as $`\sigma ^5`$ and the mass accumulated in $`10^{10}`$ yr is of the right order to reproduce the $`M_{}\sigma `$ relation.

A number of consequences follow from such high rates of stellar feeding. Tidal disruption events at the present epoch could be significantly more frequent than in models based on collisional loss-cone repopulation. The decay rate of a binary SBH could be signficantly enhanced compared with the rate in a spherical or axisymmetric nucleus, allowing binary SBHs to overcome the “final parsec problem” (Milosavljevic & Merritt 2003a, ) and coalesce by emission of gravitational radiation. Even if long-lived triaxiality should turn out to be rare, we show that transient departures from axisymmetry during mergers or galaxy interactions could induce feeding at rates approaching those inferred in AGNs.

## 2 Models and Units

The models of Papers II and III were based on the density law

$`\rho _{}`$ $`=`$ $`\rho _0m^\gamma ,`$ (1a)
$`m^2`$ $`=`$ $`{\displaystyle \frac{x^2}{a^2}}+{\displaystyle \frac{y^2}{b^2}}+{\displaystyle \frac{z^2}{c^2}}`$ (1b)
with $`\gamma =(1,2)`$ and $`c/a=0.5`$. The outer boundary (an equipotential surface) was chosen to contain roughly $`100(20)`$ times the black hole mass for $`\gamma =1(2)`$. Non-evolving solutions for both values of $`\gamma `$ were found for $`T=0.25`$ (nearly oblate) and $`T=0.50`$ (maximal triaxiality); $`T=(a^2b^2)/(a^2c^2)`$ is the triaxiality index. Models with $`T=0.75`$ (nearly prolate) evolved rapidly into precisely axisymmetric shapes and are not considered here. The black hole was represented by a central point with unit mass; the long-axis scale length $`a`$ and the constant of gravitation $`G`$ were also set to unity. The scale-free nature of the mass distribution allowed us to set $`\rho _0=1`$ without loss of generality.

Our focus is on models which include chaotic orbits. Four such models were presented in Paper III. Table 1 gives the axis ratios and the mass fractions on chaotic, tube and pyramid (box) orbits for these four models. (For definitions of the orbit families, see Papers I-III.) Chaotic orbits and tube orbits contributed roughly equally to the total mass in the self-consistent solutions, with regular box orbits a distant third. We therefore do not treat the box orbits separately from the chaotic orbits in what follows. The mass fraction in “centrophilic” orbits (chaotic orbits and pyramids) was $`50\%`$ in all models.

It is convenient to define equivalent spherical models for the stellar mass distribution; as shown in Paper III, these spherical models have energy distributions $`N(E)`$ similar to those of the triaxial models and are useful when scaling the results found here to real galaxies. The equivalent spherical models have $`\rho _{}(r)=(r/\delta )^\gamma `$ with $`\delta =(abc)^{1/3}`$, or $`\delta =0.734(0.767),T=0.5(0.25)`$. Their gravitational potentials in model units are

$`\mathrm{\Phi }_{}(r)`$ $`=`$ $`2\pi \delta r,\gamma =1`$ (2a)
$`=`$ $`4\pi \delta ^2\left[\mathrm{ln}\left({\displaystyle \frac{r}{\delta }}\right)1\right],\gamma =2.`$ (2b)
The additive constants in these expressions for the potential have been chosen in the same way as for the triaxial models (see Paper I); we note that, for $`\gamma =2`$, the zero point of the potential occurs at a radius containing $`4\pi e\delta ^3`$ times the black hole mass $`M_{}`$, or $`13.51(15.41)M_{}`$ for $`T=0.5(0.25)`$. Note that the black hole has been omitted in these expressions for the potential.

We define $`r_h`$ to be the radius in the spherical model containing a mass in stars equal to twice the black hole mass (cf. Merritt 2003). This definition is equivalent to the standard one, $`r_h=GM_{}/\sigma ^2`$, when $`\gamma =2`$, the singular isothermal sphere. For $`\gamma =1`$, $`r_h=(\pi \delta )^{1/2}=0.659(0.643)`$ for $`T=0.5(0.25)`$, while for $`\gamma =2`$, $`r_h=(2\pi \delta ^2)^1=0.296(0.270)`$. The 1D stellar velocity dispersion $`\sigma `$ is equal to $`\sqrt{2\pi }\delta =1.84(1.92)`$ in model units for $`\gamma =2`$. The velocity dispersion in a nucleus with $`\gamma =1`$ depends on radius and is zero at the center in the absence of a black hole (Dehnen, 1993); hence we do not quote a value for $`\sigma `$ when $`\gamma =1`$.

We define $`E_h\mathrm{\Phi }(r_h)`$. In model units, $`E_h=3.04(3.10)`$ for $`\gamma =1`$ and $`12.93(15.01)`$ for $`\gamma =2`$.

The results presented below can be related to real galaxies using the following scale factors for mass, length and time, derived from the spherical models just defined. Numerical constants, when given, are for $`T=0.5`$ and differ only slightly for $`T=0.25`$.

$$\left[M\right]=M_{}$$
(3)
$`\left[L\right]`$ $`=`$ $`\left(\pi \delta \right)^{1/2}r_h152\mathrm{pc}\left({\displaystyle \frac{r_h}{100\mathrm{pc}}}\right),\gamma =1`$ (4a)
$`=`$ $`\left(2\pi \delta ^2\right)r_h37.4\mathrm{pc}\left({\displaystyle \frac{\sigma }{200\mathrm{km}\mathrm{s}^1}}\right)^2\left({\displaystyle \frac{M_{}}{10^8_{}}}\right),\gamma =2`$ (4b)
$`\left[T\right]`$ $`=`$ $`\left(\pi \delta \right)^{3/4}\sqrt{{\displaystyle \frac{r_h^3}{GM_{}}}}2.79\times 10^6\mathrm{yr}\left({\displaystyle \frac{r_h}{100\mathrm{pc}}}\right)^{3/2}\left({\displaystyle \frac{M_{}}{10^8_{}}}\right)^{1/2},\gamma =1`$ (5a)
$`=`$ $`\left(2\pi \delta ^2\right)^{3/2}\sqrt{{\displaystyle \frac{r_h^3}{GM_{}}}}=3.30\times 10^5\mathrm{yr}\left({\displaystyle \frac{\sigma }{200\mathrm{kms}^1}}\right)^3\left({\displaystyle \frac{M_{}}{10^8_{}}}\right),\gamma =2.`$ (5b)
In equation (5b), the relation $`\sigma ^2=GM_{}/r_h`$ was used.

The tidal disruption radius $`r_t`$ of a $`10^8_{}`$ black hole is roughly equal to its Schwarzschild radius $`r_s=2GM_{}/c^2`$. In model units, the Schwarzschild radius is

$`r_s`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{\pi \delta }}}{\displaystyle \frac{2GM_{}}{c^2r_h}}6.31\times 10^8\left({\displaystyle \frac{r_h}{100\mathrm{pc}}}\right)^1\left({\displaystyle \frac{M_{}}{10^8_{}}}\right),\gamma =1`$ (6a)
$`=`$ $`{\displaystyle \frac{1}{\pi \delta ^2}}{\displaystyle \frac{GM_{}}{c^2r_h}}0.591\left({\displaystyle \frac{\sigma }{c}}\right)^22.63\times 10^7\left({\displaystyle \frac{\sigma }{200\mathrm{km}\mathrm{s}^1}}\right)^2,\gamma =2.`$ (6b)
In order to keep the notation as simple as possible, the same symbols will be used for both dimensional and dimensionless quantities in what follows. When not noted explicitly, the distinction will be clear from the form of the equations.

## 3 Pericenter Distributions

Chaotic orbits pass near the central singularity once per crossing time, although the pericenter distance – the distance of closest approach to the black hole – varies quasi-randomly from passage to passage. We computed the statistics of pericenter passages by carrying out integrations of chaotic orbits for $`10^5T_D`$ in each of the four triaxial potentials; $`T_D`$ is the energy-dependent dynamical time defined in Paper I. To speed the integrations, the potential and forces due to the stars (equation 1b) were expressed in terms of a truncated basis set (Hernquist & Ostriker, 1992); Ten radial and angular functions were used. The routine RADAU (Hairer & Wanner, 1996) was used for the integrations. Pericenter distances were computed by locating two times $`(t_1,t_2)`$ between which the distance to the black hole reached a minimum and interpolating the solution on a fine grid in $`t_1<t<t_2`$ to find the distance of closest approach.

Figure 1 shows the cumulative distribution of pericenter distances for two chaotic orbits integrated in the triaxial potentials with $`\gamma =(1,2)`$ and $`T=0.5`$. The number of central encounters per unit time with pericenter distances less than $`d`$, $`N(r_p<d)`$, is nearly linear with $`d`$ over its entire range. The linear dependence was observed to extend down to pericenter distances of $`10^6`$ or smaller, of order the tidal disruption radius in model units. The approximately linear dependence of $`N_E`$ on $`d`$, combined with the gravitational focussing equation

$$r_b^2=d^2\left(1+\frac{2GM_{}}{V^2d}\right)\frac{2GM_{}}{V^2}d,$$
(7)
with $`r_b`$ the impact parameter, implies $`N(r_p<r_b)r_b^2`$, i. e. a nearly uniform filling of the two-dimensional cross section that defines the “throat” of the chaotic orbit (cf. Gerhard & Binney 1985).

Chaotic orbits were observed to quickly fill the configuration-space region accessible to them, and one therefore expects that the statistical properties of one chaotic orbit at a given energy will be the same as those of any other chaotic orbit at the same energy. We verified for a number of chaotic orbits that $`N(r_p<d)`$ is indeed only weakly dependent on the starting point for fixed energy $`E`$. Hence we can write $`N(r_p<d)=N_E(r_p<d)`$, the rate of pericenter passages for a star on a chaotic orbit of energy $`E`$.

We computed $`N_E(r_p<d)`$ at the energies associated with each of the mass shells defined in Paper II. Let $`A(E)d`$ be the rate at which a single star on a chaotic orbit of energy $`E`$ experiences pericenter passages with $`r_p<d`$. The dependence of $`A`$ on $`E`$ in the four triaxial potentials is shown in Figure 2. For $`\gamma =2`$, an exponential fits the data well:

$$\mathrm{ln}Aa+bE,$$
(8)
with

$`a=0.603,b=0.290(T=0.5)`$ (9a)
$`a=0.655,b=0.275(T=0.25).`$ (9b)
We note that $`b1/2\pi \delta ^2`$. Henceforth we set $`b1/2\pi \delta ^2`$; thus $`Ae^{E/\sigma ^2}`$ and $`A`$ scales with radius as $`r^2`$. This simple scaling will be useful in what follows. In physical units,

$$A(E)1.2\frac{\sigma ^5}{G^2M_{}^2}e^{(EE_h)/\sigma ^2}.$$
(10)
Figure 2 shows that this relation holds even at energies $`E\begin{array}{c}<\hfill \\ \hfill \end{array}E_h`$.

For $`\gamma =1`$, a power law provides a good fit in the energy range $`E\begin{array}{c}>\hfill \\ \hfill \end{array}E_h`$:

$$\mathrm{ln}Aa+b\mathrm{ln}E,$$
(11)
with

$`a=1.734,b=1.388(T=0.5)`$ (12a)
$`a=1.869,b=1.357(T=0.25).`$ (12b)
In physical units,

$$A(E)0.7\sqrt{\frac{GM_{}}{r_h^5}}\left(\frac{E}{E_h}\right)^{1.4}.$$
(13)
Thus $`Ar^{1.4}`$, $`r\begin{array}{c}>\hfill \\ \hfill \end{array}r_h`$.

The pericenter distribution at any $`E`$ is also characterized by a second quantity, the maximum pericenter distance $`r_{p,max}(E)`$ reached by chaotic orbits of energy $`E`$. This distance corresponds roughly to the width of the “throat” of a box orbit in an integrable triaxial potential. We find, for $`\gamma =2`$,

$$\mathrm{ln}r_{p,max}c+dE,$$
(14)
with

$`c=0.43,d=0.146(T=0.5)`$ (15a)
$`c=0.36,d=0.141(T=0.25).`$ (15b)
We note that $`d=1/4\pi \delta ^2`$ to within its uncertainties; thus

$$r_{p,max}0.3r_he^{(EE_h)/2\sigma ^2}$$
(16)
and $`r_{p,max}`$ scales approximately linearly with apocenter distance, $`r_{p,max}0.2r_{apo}`$.

For $`\gamma =1`$, we find, for $`E\begin{array}{c}>\hfill \\ \hfill \end{array}E_h`$,

$$\mathrm{ln}r_{p,max}c+d\mathrm{ln}E,$$
(17)
with

$`c=1.86,d=0.86(T=0.5)`$ (18a)
$`c=2.03,d=0.90(T=0.25).`$ (18b)
This implies

$$r_{p,max}0.60r_h\left(\frac{E}{E_h}\right)^{0.88}$$
(19)
and again $`r_{p,max}`$ scales roughly linearly with $`r_{apo}`$.

## 4 Feeding Rates

The feeding rates implied by Figure 2 are high compared with those in diffusive loss cone models; in fact they are comparable to the maximum possible capture rates in spherical models, the so-called “full loss cone” rate. We illustrate this by considering the rate at which stars on a single orbit pass within a distance $`r_t`$ of the center. In the triaxial models, this rate is $`A(E)r_t`$. Setting $`E`$ to the energy of a circular orbit at $`r_h`$, taking $`r_t2\times 10^7`$ in model units (equation 6b for $`\gamma =2`$), and using the expressions given above for $`A(E)`$, we find a rate of $`2\times 10^6T_D^1`$. In other words, roughly $`10^6`$ orbital periods are required for a star on the edge of the black hole’s sphere of influence to pass within its tidal radius. In a spherical model with a full loss cone, a fraction $`r_tr_h/r^2`$ of the stars at $`r`$ come within $`r_t`$ each orbital period; hence the mean capture rate from stars near $`r_h`$ is $`(r_t/r_h)T_D^1`$ or $`1.5\times 10^6T_D^1`$. In the spherical geometry, the net rate of pericenter passages is due to a small population of stars (those with sufficiently low angular momenta) being lost over a short period of time (an orbital period). In the triaxial geometry, a large population of stars (all stars on chaotic orbits) are supplied to the center over a longer period of time ($`10^5`$ periods). Of course, in the absence of loss cone repopulation, capture rates in the spherical geometry would drop to zero after a single orbital period.

Next we compute the full feeding rates in the triaxial models. Let $`_c(E,t)dE`$ be the mass in stars on chaotic orbits with energies from $`E`$ to $`E+dE`$; the time dependence reflects a possible loss of stars due to tidal disruption or capture. If the capture or disruption radius is $`r_t`$, the energy-dependent loss rate is

$`\dot{}_c(E,t)dE`$ $`=`$ $`r_tA(E)_c(E,t)dE`$ (20a)
$`=`$ $`r_tA(E)_c(E,0)e^{A(E)r_tt}dE`$ (20b)
and the total capture rate from chaotic orbits at all energies is

$$\dot{M}(t)=r_tA(E)_c(E,0)e^{A(E)r_tt}𝑑E.$$
(21)
We evaluated these expressions in two ways. $`_c(E,0)`$ can be computed directly from the orbital weights in the Schwarzschild solutions (cf. Figure 3 of Paper III). Alternatively, smooth approximations to $`_c(E,0)`$ can be computed using the equivalent spherical models defined above. In the latter case, we define $`_c(E,0)f_c(E)(E)`$ with $`(E)`$ the energy distribution in the equivalent spherical model and $`f_c(E)`$ the fraction of stars at energy $`E`$ assumed to be on chaotic orbits. For $`\gamma =1`$, the energy distribution in the spherical geometry is

$$(E)=\frac{8}{35}\frac{r_h^2}{G^2M_{}}E$$
(22a)
or, in model units,
$$(E)=\frac{8}{35\pi \delta }E.$$
(22b)
For $`\gamma =2`$, we have

$$(E)=\frac{2\sqrt{6}}{9}\frac{r_h}{G}e^{(EE_h)/2\sigma ^2}$$
(23a)
or, in model units,
$`(E)={\displaystyle \frac{2\sqrt{6}}{9}}e\delta e^{E/4\pi \delta ^2}.`$ (23b)
Since these approximations to $`(E)`$ were derived from expressions for the potential that exclude the central point mass, they are only strictly valid outside the black hole’s sphere of influence, $`E\begin{array}{c}>\hfill \\ \hfill \end{array}E_h`$. However the expression for $`\gamma =2`$ turns out to be reasonably correct even at lower energies (cf. Figure 3 of Paper III), a result that will be used below.

Figure 3 shows $`\dot{}_c(E,0)`$ computed in both ways. The chaotic orbit fraction $`f_c(E)`$ was set to a constant, $`\overline{f_c}`$, in each model; no attempt was made to match the detailed energy dependence of the chaotic orbital populations in the Schwarzschild solutions. Nevertheless the fits are reasonably good, particularly for $`\gamma =2`$. For $`\gamma =1`$, the orbital distributions in the Schwarzschild solutions are “noisy” but the analytic expression does a reasonable job of reproducing the mean dependence. The values of $`\overline{f_c}`$ used in Figure 3 were $`\overline{f_c}=0.25`$ ($`\gamma =2,T=0.5`$ and $`0.25`$), $`\overline{f_c}=0.5`$ ($`\gamma =1,T=0.5`$) and $`\overline{f_c}=0.4`$ ($`\gamma =1,T=0.25`$). These values are slightly smaller than the overall chaotic mass fractions in the numerical solutions (Table 1); the reason is that the Schwarzschild solutions put a large number of orbits in the outermost shells to compensate for the hard outer edge, and are correspondingly depleted at low energies (Figure 3 of Paper III).

Also shown in Figure 3 are the feeding rates predicted by a spherical, “full loss cone” model. The full-loss-cone feeding rate is

$`\dot{M}_{full}(E)dE`$ $`=`$ $`{\displaystyle \frac{_{lc}(E)}{P(E)}}dE`$ (24a)
$`=`$ $`4\pi ^2f(E)J_{lc}^2(E)dE`$ (24b)
where $`_{lc}(E)`$ is the number of stars in the spherical model with pericenters below $`r_t`$, $`f(E)`$ is the isotropic distribution function, $`P(E)`$ is the radial period, and $`J_{lc}`$ is the angular momentum of a star with pericenter $`r_t`$, $`J_{lc}^2=2r_t^2(E\mathrm{\Phi }(r_t))2GM_{}r_t`$. Equation (24b) describes the capture rate in a spherical model assuming that all of the orbits remain fully populated. Figure 3 shows that the capture rates in the triaxial models are similar to the full-loss-cone rates in the equivalent spherical models, and can even exceed them at high energies.

In a real spherical galaxy, the high feeding rates corresponding to a full loss cone would persist for just a single orbital period at each energy. In the triaxial models, by contrast, the capture rate from chaotic orbits decreases only by a factor $`\mathrm{exp}\left[A(E)P(E)r_t\right]`$ in one orbital period. At energies $`\begin{array}{c}>\hfill \\ \hfill \end{array}E_h`$, this factor is negligible meaning that feeding rates will remain large for many orbital periods.

The evolution of $`\dot{}_c(E,t)`$ is shown in Figure 4. Here we have set the capture radius $`r_t`$ to its smallest possible value, $`r_t=r_s`$, using equation (6b). For $`\gamma =1`$, this requires a choice for $`M_{}`$ and $`r_h`$; we took $`M_{}=10^8_{}`$ and $`r_h=100`$ pc. For $`\gamma =2`$, only $`\sigma `$ needs to be specified; we took $`\sigma =200`$ km s<sup>-1</sup>. The analytic approximations, equations (22b) and (23b), were used for $`M_c(E,0)`$ and we set $`\overline{f_c}=0.5`$. Figure 4 shows that, for $`\gamma =1`$, the capture rate is low enough that no significant changes occur in $`_c(E,t)`$ even at times as great as $`t=10^6`$, corresponding to $`10^{12}`$ yr using the scaling of equation (5b). For $`\gamma =2`$, the capture rate is higher and changes in $`_c(E,t)`$ begin to occur at energies above $`E_h`$ for $`t\begin{array}{c}>\hfill \\ \hfill \end{array}10^5`$, corresponding to $`10^{10}`$ yr.

The total capture rate from orbits at all energies is $`\dot{M}(t)=\dot{}_c(E,t)𝑑E`$. For $`\gamma =1`$ and $`r_tr_s`$, the captured mass is a small fraction of $`M_{}`$ and $`\dot{m}`$ is nearly independent of time. Thus we can write

$`\dot{M}`$ $``$ $`{\displaystyle _{E_1}^{E_2}}\dot{}_c(E,0)𝑑E=r_t{\displaystyle _{E_1}^{E_2}}A(E)f_c(E)(E)𝑑E`$ (25a)
$``$ $`{\displaystyle \frac{8}{35\pi \delta }}e^ar_t{\displaystyle _{E_1}^{E_2}}f_c(E)E^{1+b}𝑑E`$ (25b)
in model units. As lower integration limit we take $`E_1=0`$; this is reasonable since the value of the integral is not strongly dependent on the lower cutoff. However since $`1+b0.4`$, the integral diverges for large $`E_2`$. This corresponds physically to the fact that a $`\rho r^1`$ nucleus can not extend indefinitely. Nuclear density profiles in real, weak-cusp galaxies exhibit a break at a radius $`r_b10r_h`$. We define $`E_b\mathrm{\Phi }(r_b)`$ and set $`E_2=E_b`$, giving, in model units,

$$\dot{M}=\frac{8}{35(2+b)\pi \delta }e^ar_t\overline{f_c}E_b^{2+b}.$$
(26)
Setting $`r_t=r_s`$, this becomes, in physical units,

$$\dot{M}=(1.77,1.97)\times 10^5_{}\mathrm{yr}^1\overline{f_c}\left(\frac{r_h}{100\mathrm{pc}}\right)^{5/2}\left(\frac{M_{}}{10^8_{}}\right)^{5/2}\left(\frac{r_b}{10r_h}\right)^{(0.61,0.64)}$$
(27)
where the number pairs in parentheses refer to $`T=(0.5,0.25)`$ respectively. We note that the feeding rate is almost the same for the two values of the triaxiality index, suggesting a weak dependence of $`\dot{}`$ on the degree of departure from axisymmetry, for a given chaotic mass fraction. The feeding rates of equation (27) imply a captured mass of $`\begin{array}{c}<\hfill \\ \hfill \end{array}0.01M_{}`$ over a Hubble time.

For $`\gamma =2`$, the feeding rates are much higher, and the captured mass can be a large fraction of $`M_{}`$ even when $`r_t`$ is as small as $`r_s`$. Hence the population of chaotic orbits is significantly depleted and we can not ignore the time dependence of $`\dot{m}`$. We have

$`\dot{M}(t)`$ $`=`$ $`r_t{\displaystyle A(E)_c(E,0)e^{A(E)r_tt}𝑑E}`$ (28a)
$`=`$ $`{\displaystyle \frac{2\sqrt{6}e\delta }{9}}e^ar_t{\displaystyle f_c(E)e^{E/4\pi \delta ^2}\mathrm{exp}\left(r_tte^{a+bE}\right)𝑑E}`$ (28b)
in model units. It is reasonable to set the lower integration limit to $`\mathrm{}`$, since our expressions for $`A(E)`$ and $`(E,0)`$ are valid for $`EE_h`$ and the contribution to the integral from low energies falls off quickly with time as the most-bound particles are eaten; hence any errors due to the forms of $`A(E)`$ or $`(E)`$ at low energies are transient. Setting the upper integration limit to $`+\mathrm{}`$ is also reasonable since the product $`A(E)(E,0)`$ drops rapidly with $`E`$. (Furthermore it is possible that the combined, luminous plus dark matter density profiles in early-type galaxies are well described as $`\rho r^2`$ even far outside of the nucleus.) Making the substitution $`yr_tt\mathrm{exp}(a+bE)`$, we find

$`\dot{M}(t)`$ $`=`$ $`{\displaystyle \frac{4\pi \sqrt{6}e\delta ^3}{9}}e^{a/2}\overline{f_c}\left({\displaystyle \frac{r_t}{t}}\right)^{1/2}{\displaystyle _0^{\mathrm{}}}y^{1/2}e^y𝑑y`$ (29a)
$``$ $`(4.82,5.36)\overline{f_c}\left({\displaystyle \frac{r_t}{t}}\right)^{1/2}`$ (29b)
where the paired numbers refer to $`T=(0.5,0.25)`$. We note again the small dependence of the feeding rate on the degree of triaxiality for fixed $`\overline{f}_c`$. The total mass accreted after time $`t_{acc}`$ is

$$\mathrm{\Delta }M=(9.64,10.72)\overline{f_c}\left(r_tt_{acc}\right)^{1/2}.$$
(30)
In physical units, these expressions become

$`\dot{M}(t)`$ $`=`$ $`(1.48,1.54)\overline{f_c}{\displaystyle \frac{\sigma ^3}{G}}{\displaystyle \frac{\sigma }{c}}\left({\displaystyle \frac{r_t}{r_s}}\right)^{1/2}\left({\displaystyle \frac{t}{GM_{}/\sigma ^3}}\right)^{1/2}`$ (31b)
$``$ $`(4.28,4.46)\times 10^3_{}\mathrm{yr}^1\overline{f_c}\left({\displaystyle \frac{r_t}{r_s}}\right)^{1/2}\left({\displaystyle \frac{\sigma }{200\mathrm{km}\mathrm{s}^1}}\right)^{5/2}\times `$
$`\left({\displaystyle \frac{M_{}}{10^8_{}}}\right)^{1/2}\left({\displaystyle \frac{t}{10^{10}\mathrm{yr}}}\right)^{1/2}`$
and

$`\mathrm{\Delta }M`$ $`=`$ $`(2.97,3.09)\overline{f_c}M_{}{\displaystyle \frac{\sigma }{c}}\left({\displaystyle \frac{r_t}{r_s}}\right)^{1/2}\left({\displaystyle \frac{t_{acc}}{GM_{}/c^2}}\right)^{1/2}`$ (32b)
$``$ $`(8.56,8.92)\times 10^7_{}\overline{f_c}\left({\displaystyle \frac{r_t}{r_s}}\right)^{1/2}\left({\displaystyle \frac{\sigma }{200\mathrm{kms}^1}}\right)^{5/2}\times `$
$`\left({\displaystyle \frac{M_{}}{10^8_{}}}\right)^{1/2}\left({\displaystyle \frac{t_{acc}}{10^{10}\mathrm{yr}}}\right)^{1/2}`$
The $`r_t^{1/2}`$ dependence reflects the more rapid destruction of the cusp when $`r_t`$ is large. Equation (32b) implies a captured mass of order $`M_{}`$ in a Hubble time.

## 5 Tidal Disruption Rates

The feeding rate in a steep-cusp nucleus is given by equation (31b), as a function of the capture radius $`r_t`$. If $`M_{}\begin{array}{c}<\hfill \\ \hfill \end{array}10^8_{}`$, solar type stars are tidally disrupted at a radius

$$r_t\left(\frac{M_{}}{10^8_{}}\right)^{2/3}r_s$$
(33)
(Hills, 1975). ($`10^8_{}`$ is also roughly the mass of the largest black holes that sit in steep-cusp nuclei.) Combining equations (31b) and (33), the tidal disruption rate in a steep-cusp nucleus is

$$\dot{M}_{collisionless}(t)4\times 10^3_{}\mathrm{yr}^1\overline{f_c}\sigma _{200}^{5/2}M_{,8}^{1/6}t_{10}^{1/2}$$
(34)
with $`t_{10}`$ the time since cusp formation in units of $`10^{10}`$ yr, $`\sigma _{200}`$ the velocity dispersion in units of $`200`$ km s<sup>-1</sup> and $`M_{,8}`$ the black hole mass in units of $`10^8_{}`$. We refer to this as the “collisionless” tidal disruption rate to highlight that the supply of stars to the black hole is not being driven by gravitational scattering, as in the standard model (Frank & Rees 1976; Lightman & Shapiro 1977).

Supermassive black holes in galaxies at the current epoch satisfy the $`M_{}\sigma `$ relation,

$$\left(\frac{M_{}}{10^8_{}}\right)1.48\left(\frac{\sigma }{200\mathrm{km}\mathrm{s}^1}\right)^{4.65}$$
(35)
(e.g. Merritt & Ferrarese 2001b), which allows us to write the tidal disruption rate in nearby galaxies in terms of either $`M_{}`$ or $`\sigma `$ alone:

$`\dot{M}_{collisionless}(t)`$ $``$ $`5\times 10^3_{}\mathrm{yr}^1\overline{f_c}\sigma _{200}^{3.28}t_{10}^{1/2}`$ (36a)
$``$ $`4\times 10^3_{}\mathrm{yr}^1\overline{f_c}M_{,8}^{0.70}t_{10}^{1/2}.`$ (36b)
We compare this to the rate at which collisional (encounter-driven) loss-cone refilling supplies stars to the central black hole. In a spherical, $`\rho r^2`$ nucleus,

$$\dot{M}_{collisional}=2.8\times 10^4_{}\mathrm{yr}^1\sigma _{200}^{7/2}M_{,8}^1$$
(37)
(Wang & Merritt 2003); this expression again assumes that the disrupted stars have the mass and radius of the sun. Equation (37) can likewise be rewritten using the $`M_{}\sigma `$ relation:

$`\dot{M}_{collisional}`$ $``$ $`2\times 10^4_{}\mathrm{yr}^1\sigma _{200}^{1.15}`$ (38a)
$``$ $`2\times 10^4_{}\mathrm{yr}^1M_{,8}^{0.25}.`$ (38b)
Figure 5 plots $`\dot{M}(M_{})`$ in the collisionless (triaxial) and collisional (spherical) cases for a $`\rho r^2`$ nucleus. The two rates scale in the opposite sense with $`M_{}`$ and collisional feeding would dominate in sufficiently small galaxies. The rates are equal when

$$M_{}5\times 10^6_{}\overline{f_c}^1t_{10}^{0.5}.$$
(39)
For $`t_{10}^{1/2}/\overline{f_c}1`$, this mass is close to that of the smallest black holes with reliably determined masses, in the Milky Way and M32. Between this mass and $`M_{}10^8_{}`$, Figure 5 suggests that loss cone feeding driven by triaxiality can easily dominate collisional loss cone feeding, as long as the chaotic mass fraction $`f_c`$ is not too much smaller than one. Flaring rates should peak in the brightest galaxies with steep nuclear density profiles, at values of a few times $`10^3\mathrm{yr}^1`$. Since the disruption rate in the triaxial geometry varies as $`t^{0.5}(1+z)^{0.5}`$, flaring could be even more important at intermediate redshifts.

In a $`\rho r^1`$ nucleus, and again assuming $`M_{}\begin{array}{c}<\hfill \\ \hfill \end{array}10^8_{}`$, equation (27) gives a disruption rate for solar-type stars of

$$\dot{M}_{collisionless}2\times 10^5_{}\mathrm{yr}^1\overline{f_c}\left(\frac{r_h}{100\mathrm{pc}}\right)^{5/2}\left(\frac{M_{}}{10^8_{}}\right)^{11/6}\left(\frac{r_b}{10r_h}\right)^{5/8},$$
(40)
approximately independent of time. The corresponding collisional expression is

$$\dot{M}_{collisional}1\times 10^6_{}\mathrm{yr}^1\overline{f_c}\left(\frac{r_h}{100\mathrm{pc}}\right)^{7/3}\left(\frac{M_{}}{10^8_{}}\right)^{29/18}.$$
(41)
The functional dependence in equation (41) is taken from Wang & Merritt (2003), and the normalizing factor is based on those authors’ calculation of $`\dot{N}`$ for the galaxy NGC 3379 ($`\gamma 1.1`$). For a $`\gamma =1`$ nucleus, the collisional and collisionless disruption rates scale in almost the same way with $`M_{}`$ and $`r_h`$, and the collisionless rate exceeds the collisional rate as long as $`\overline{f_c}\begin{array}{c}>\hfill \\ \hfill \end{array}0.05`$.

In summary: even small fractional populations of chaotic orbits can produce tidal disruption rates that exceed those predicted by the standard model of collisional loss cone repopulation in a spherical nucleus (Syer & Ulmer 1999; Magorrian & Tremaine 1999; Wang & Merritt 2003). In triaxial nuclei with $`\rho r^2`$ and $`M_{}10^8_{}`$, tidal flaring rates can plausibly reach values as high as several times $`10^3`$ yr<sup>-1</sup> for solar type stars. Rates could be even higher if the cusps were recently formed.

## 6 Black Hole Growth and the $`M_{}\sigma `$ Relation

The fate of gas liberated by the tidal disruption of a star at $`r_t>r_s`$ is uncertain. But for $`M_{}\begin{array}{c}>\hfill \\ \hfill \end{array}10^8_{}`$, stars are directly captured by the black hole without being disrupted. Using equation (32b), setting $`r_t=r_s`$ and requiring $`\mathrm{\Delta }M=M_{}`$, we find the following relation between $`M_{}`$ and $`\sigma `$ in a steep triaxial cusp:

$$\frac{M_{}}{10^8_{}}0.8\overline{f_c}\left(\frac{t_{acc}}{10^{10}\mathrm{yr}}\right)\left(\frac{\sigma }{200\mathrm{km}\mathrm{s}^1}\right)^5.$$
(42)
This is remarkably similar to the $`M_{}\sigma `$ relation: the exponent on $`\sigma `$ is consistent with measured values, $`4.5\pm 0.5`$ (Merritt & Ferrarese 2001a, ) and even the normalization is of the right order if $`\overline{f_c}1`$.

The $`\dot{M}\sigma ^5`$ dependence can be understood in several ways. Writing $`\dot{}(E,0)=r_sA(E)_c(E,0)`$ and integrating over energies from $`E_h`$ to infinity, we find

$$\dot{M}(EE_h)2.5f_c\frac{\sigma ^3}{G}\left(\frac{\sigma }{c}\right)^2.$$
(43)
Alternatively, the total rate at which stars pass inward at radius $`r`$ is $`4\pi r^2\rho (r)\sigma =2\sigma ^3/G`$, and a fraction $`F=T_DAf_cr_s`$ of stars at energy $`E`$ pass within $`r_s`$ each crossing time. Equation (10) gives

$$A(r)\frac{\sigma ^5}{G^2M_{}^2}\left(\frac{r}{r_h}\right)^2$$
(44)
and $`T_D(r_h/\sigma )(r/r_h)`$, so that

$$\dot{M}(r)\overline{f_c}\frac{\sigma ^3}{G}\left(\frac{\sigma }{c}\right)^2\frac{r_h}{r}.$$
(45)
These relations show that the feeding rate due to stars originating outside of the SBH’s sphere of influence is of order

$$\frac{\sigma ^5}{Gc^2}10^3_{}\mathrm{yr}^1\left(\frac{\sigma }{200\mathrm{km}\mathrm{s}^1}\right)^5,$$
(46)
independent of $`M_{}`$ – large enough to contribute substantially to the masses of SBHs if accretion continues for $`10^9`$ yr or more.

For $`M_{}\begin{array}{c}<\hfill \\ \hfill \end{array}10^8_{}`$, the tidal disruption radius exceeds $`r_s`$ and scales as $`M_{}^{1/3}`$ (equation 33). If the growth of the black hole is determined by the rate of tidal disruptions, $`\mathrm{\Delta }M`$ scales more weakly with $`\sigma `$ in this regime, roughly as $`\sigma ^3`$. This would imply a flattening of the $`M_{}\sigma `$ relation at masses below $`10^8_{}`$. However it is not clear what the fate of tidally-liberated gas would be; some would fall into the black hole but some would escape (Frank, 1979). For masses below $`10^6_{}`$, Figure 5 suggests that collisional loss cone repopulation would dominate over collisionless feeding.

Equation (46) is similar to equation (8) of Zhao, Haehnelt & Rees (2002), who considered the consequences for black hole growth of a continuously resupplied loss cone in the spherical geometry. The agreement is reasonable since the feeding rates in the triaxial models are comparable to spherical, full-loss-cone rates for $`E\begin{array}{c}>\hfill \\ \hfill \end{array}E_h`$ (Figure 3). The triaxial models are different in one respect. In the spherical geometry, the contribution per orbital period to $`\mathrm{\Delta }M`$ from stars at radius $`r`$ scales as $`r^2`$. In the triaxial models, every centrophilic orbit is capable of visiting the center and the probability that a single star will do so in one orbital period varies as $`r^1`$. Thus while the feeding rates from stars at $`rr_h`$ are comparable in the two models, large-radius orbits contribute relatively more in the triaxial models. One consequence is that the growth of the black hole is not so strongly dependent on the stellar distribution near $`r_h`$.

Accretion of stars, or of gas liberated from stars, have long been discussed as mechanisms for growing SBHs in galactic nuclei (Hills 1975; Young, Shields & Wheeler 1977; Frank 1978; McMillan, Lightman & Cohn 1981). Hills pointed out already in 1975 that feeding at the full-loss-cone rate could plausibly grow a $`10^8_{}`$ black hole in a Hubble time. Subsequent authors (e.g. Young, Shields & Wheeler) tended criticize this claim on the grounds that gravitational encounters occur too infrequently to maintain a fully-populated loss cone. Our result – that sustained feeding rates in the triaxial geometry can be comparable to full-loss-cone rates in the spherical geometry – should help to revive Hills’ model. However, chaotic loss cones still fail to reproduce in a natural way the observed time dependence of quasar luminosities. The $`t^{1/2}`$ time dependence found here for stellar feeding in a $`\rho r^2`$ density cusp (equation 31b) – which results entirely from the depopulation of chaotic orbits – is much more gradual than the observed, $`(1+z)^3`$ falloff of quasar luminosities (e.g. Boyle, Shanks & Peterson 1988), and would require that much of the growth took place long after the quasar epoch. Slow, optically faint growth of SBHs has been proposed (Richstone et al. 1998), but only before it was discovered that SBH masses in nearby galaxies had been substantially overestimated (Merritt & Ferrarese 2001b, ). SBH masses are now consistent with accretion rates inferred from quasar statistics (Merritt & Ferrarese 2001c, ) leaving less room for growth at late times.

On the other hand, the time dependence derived here assumes that the steep density cusp forms only once and that it is slowly depleted by accretion of centrophilic orbits. Cusps might form episodically via star formation during mergers leading to repeated bursts of stellar fueling and larger average rates of accretion (Milosavljevic & Merritt 2003b, ). Other considerations could boost the mean capture rate as well. Growth of the SBH via accretion steepens the stellar density profile, to $`\rho r^{2.5}`$ near the SBH (Merritt, 2003), an effect that was ignored here. Depopulation of the chaotic orbits would be counteracted to some extent by scattering of stars onto these orbits. Dark matter could also be accreted (Zhao, Haehnelt & Rees 2002; Read & Gilmore 2003). Once the supply of star-forming gas was depleted (Kauffmann & Haehnelt, 2000), subsequent mergers would turn off the stellar feeding through the formation of a binary SBH which ejects stars from the nucleus and lowers its density (Milosavljevic & Merritt, 2001). This might be a natural way to explain the rapid observed falloff in quasar luminosities. In summary, while capture of stars is probably not the only way in which black holes grow, there would seem to be no reason to rule it out as an important contributor.

## 7 Transient Changes in the Potential

Deviations from axial symmetry are likely to occur at least temporarily following a merger or accretion event, yielding chaotic orbits near the black hole and resulting in enhanced feeding rates (Norman & Silk, 1983). Assuming a steep cusp (perhaps formed dissipatively during the merger), equation (32b) and the $`M_{}\sigma `$ relation imply an accreted mass of

$$\mathrm{\Delta }M2\times 10^7_{}\overline{f_c}\left(\frac{r_t}{r_s}\right)^{0.5}\left(\frac{\sigma }{200\mathrm{km}\mathrm{s}^1}\right)^5\left(\frac{\mathrm{\Delta }t}{10^8\mathrm{yr}}\right)^{0.5}$$
(47)
in time $`\mathrm{\Delta }t`$, or a mean feeding rate over $`\mathrm{\Delta }t`$ of

$$\frac{\mathrm{\Delta }M}{\mathrm{\Delta }t}0.2_{}\mathrm{yr}^1\overline{f_c}\left(\frac{r_t}{r_s}\right)^{0.5}\left(\frac{\sigma }{200\mathrm{km}\mathrm{s}^1}\right)^5\left(\frac{\mathrm{\Delta }t}{10^8\mathrm{yr}}\right)^{0.5}.$$
(48)
This is comparable to feeding rates inferred in AGNs for $`\sigma \begin{array}{c}<\hfill \\ \hfill \end{array}250`$ km s<sup>-1</sup> and in quasars for $`\sigma \begin{array}{c}>\hfill \\ \hfill \end{array}350`$ km s<sup>-1</sup>.

## 8 Decay of a Black-Hole Binary

Following the merger of two galaxies each containing a SBH, a binary forms with semi-major axis $`a=a_{hard}`$, where

$`a_{hard}`$ $``$ $`{\displaystyle \frac{G\mu }{4\sigma ^2}}2.7\mathrm{pc}(1+p)^1\left({\displaystyle \frac{m_2}{10^7_{}}}\right)\left({\displaystyle \frac{\sigma }{200\mathrm{km}\mathrm{s}^1}}\right)^2`$ (49a)
$`=`$ $`{\displaystyle \frac{p}{4(1+p)^2}}r_h;`$ (49b)
$`m_2`$ is the mass of the smaller black hole, $`pm_2/m_1`$ and $`\mu m_1m_2/M_{}`$ (Merritt, 2003). In a spherical or axisymmetric galaxy, the binary quickly ejects stars with pericenters below $`a_{hard}`$, after which changes in the binary separation take place on the much longer time scale associated with collisional refilling of the binary’s loss cone Yu (2002). In all but the densest nuclei, decay stalls at a separation $`a\begin{array}{c}>\hfill \\ \hfill \end{array}0.1a_{hard}`$, too large for the efficient emission of gravitational waves (Merritt, 2003). In a triaxial galaxy, decay can continue due to the much larger supply of centrophilic stars. The decay rate is determined by the supply of stars into a region of radius $`a`$ around the center of mass of the binary (Quinlan, 1996). In a steep triaxial cusp, the feeding rate from stars of energy $`E>\mathrm{\Phi }(a_{hard})`$ into a region of radius $`a_{hard}`$ is given by an equation similar to (43):

$`\dot{M}(E>E_{hard})`$ $``$ $`a_{hard}{\displaystyle _{E_{hard}}^{\mathrm{}}}A(E)M_c(E)𝑑E`$ (50a)
$``$ $`1.3\overline{f_c}{\displaystyle \frac{\sigma ^3}{G}}`$ (50b)
$``$ $`2500_{}\mathrm{yr}^1\overline{f_c}\left({\displaystyle \frac{\sigma }{200\mathrm{km}\mathrm{s}^3}}\right)^3,`$ (50c)
independent of $`M_{}`$ and $`p`$. If this rate were maintained, the binary would interact with its own mass in stars in a time of only $`10^5`$ yr. In fact, the feeding rate will decline with time as the centrophilic orbits are depleted. Equating the energy carried away by stars with the change in the binary’s binding energy gives

$$\frac{3}{2}\frac{G\mu }{a}dM\frac{Gm_1m_2}{2}d\left(\frac{1}{a}\right)$$
(51)
(Merritt, 2003). The coupled equations describing the change in the binary separation and the evolution of the stellar distribution are then

$`{\displaystyle \frac{d(1/a)}{dt}}`$ $`=`$ $`{\displaystyle \frac{3}{M_{}}}{\displaystyle \frac{r_t}{a}}{\displaystyle A(E)_c(E,t)𝑑E},`$ (52a)
$`{\displaystyle \frac{d_c}{dt}}`$ $`=`$ $`r_tA(E)_c(E,t)`$ (52b)
where $`r_t`$ is the pericenter distance below which stars can efficiently exchange energy with the binary. Note that we are ignoring any contribution to shrinkage of the binary from low-angular-momentum regular orbits, or from collisional repopulation of the loss cone.

We set $`r_t`$ to a (small) multiple of $`a`$, $`r_t(t)=a(t)`$, $`1`$, reflecting the fact that only stars which come to within roughly a binary separation of the black holes will interact strongly enough with them to be ejected. As lower integration limit on the integral in equation (52a), we take $`E_{hard}`$; this choice reflects the fact that the stellar distribution immediately following the formation of a hard binary is similar to that in the pre-merger galaxies outside of $`ra_{hard}`$ (Milosavljevic & Merritt, 2001). The evolution equations become

$`{\displaystyle \frac{d(1/a)}{dt}}`$ $`=`$ $`{\displaystyle \frac{3}{M_{}}}{\displaystyle _{E_{hard}}^{\mathrm{}}}A(E)_c(E,t)𝑑E,`$ (53a)
$`{\displaystyle \frac{d_c}{dt}}`$ $`=`$ $`a(t)A(E)_c(E,t).`$ (53b)
As initial conditions we take $`a(0)=a_{hard}`$ and $`_c(E,0)=\overline{f_c}(E)`$ with $`(E)`$ given by equation (23b).

Figure 6 shows solutions for various choices of $`\overline{f_c}`$ and for $`pm_2/m_1=(1,0.1)`$. At late times, the binary separation varies approximately as

$$\frac{a_{hard}}{a}\overline{f_c}^2\times t$$
(54)
in model units, for both values of $`p`$. In physical units,

$$\frac{a_{hard}}{a}3\times 10^4\overline{f_c}^2\left(\frac{\sigma }{200\mathrm{km}\mathrm{s}^1}\right)^3\left(\frac{M_{}}{10^8_{}}\right)^1\left(\frac{t}{10^{10}\mathrm{yr}}\right).$$
(55)
A $`1/at`$ dependence is expected in “full loss cone” situations like this one (Milosavljevic & Merritt 2003b). The nonlinear dependence of $`a(t)`$ on $`\overline{f_c}`$ reflects the fact that, when $`\overline{f_c}`$ is small, the decay rate is also small and the binary tends to stall at a large separation, causing it to interact with and deplete the supply of stars, further slowing the decay etc.

Coalescence due to gravitational wave emission in a time $`t_{gr}`$ occurs when $`a=a_{gr}`$, where

$$\frac{a_{hard}}{a_{gr}}75\frac{p^{3/4}}{\left(1+p\right)^{3/2}}\left(\frac{\sigma }{200\mathrm{km}\mathrm{s}^1}\right)^{7/8}\left(\frac{t_{gr}}{10^9\mathrm{yr}}\right)^{1/4}$$
(56)
for a circular-orbit binary (Merritt, 2003). For $`p=1(0.1)`$, equation (56) implies that decay by a factor of $`10^2(10^1)`$ is required in order for gravity-wave coalescence to occur in $`10^9`$ yr. Combining equations (55) and (56), this can be achieved in $`10^{10}`$ yr for chaotic mass fractions:

$$\overline{f_c}0.05\left(\frac{\sigma }{200\mathrm{km}\mathrm{s}^1}\right)^{9/16}$$
(57)
with a weak dependence on $`p`$. Thus, placing just a few percent of a galaxy’s mass on chaotic orbits is sufficient to overcome the “final parsec problem” (Milosavljevic & Merritt 2003a, ) and induce coalesence, if (as assumed in Figure 6) the chaotic orbits are present initially at energies much greater than $`E_h`$. Since the stars supplied to the binary come from well outside of $`r_h`$ in this case, the effect on the stellar density profile would be negligible even after the binary had ejected several times its own mass in stars.

## 9 Summary

In a $`\rho r^2`$ triaxial nucleus containing chaotic (centrophilic) orbits and a central black hole, stars are supplied to the black hole’s sphere of influence at a rate $`f_c\sigma ^3/G`$, and to the black hole itself at a rate $`f_c\sigma ^5/Gc^2`$, where $`f_c`$ is the fraction of the mass on centrophilic orbits and $`\sigma `$ is the stellar velocity dispersion. If $`f_c`$ is of order unity, as in recently published self-consistent models (Papers II, III), the mass accumulated in $`10^{10}`$ yr is of the same order as observed black hole masses. Feeding rates in such an environment fall off as $`t^{1/2}`$ as the chaotic orbits are depleted. Tidal disruption rates at the current epoch due to stars on chaotic orbits are as much as 1-2 orders of magnitude greater than in standard, spherical or axisymmetric models in which loss-cone refilling occurs via gravitational scattering onto eccentric orbits. Transient changes in the shape of a nucleus during mergers or accretion events could lead to episodic stellar feeding at rates of $`0.22_{}`$ yr<sup>-1</sup>, comparable to the rates inferred in active galactic nuclei and quasars. Decay of a black-hole binary at the center of a triaxial nucleus could be greatly enhanced compared with spherical or axisymmetric nuclei, even if only a few percent of the mass is on chaotic orbits, thus solving the “final parsec problem.”

This work was supported by NSF grants AST 00-71099 and AST 02-0631, and by NASA grants NAG5-6037 and NAG5-9046. M.Y. Poon is grateful to the Croucher Foundation for a postdoctoral fellowship.