text:
# Long Term Evolution of Massive Black Hole Binaries

## 1. Introduction

Larger galaxies grow through the agglomeration of smaller galaxies and protogalactic fragments. If more than one of the fragments contains a massive black hole (MBH), the MBHs will form a bound system in the merger product (Begelman, Blandford & Rees, 1980; Roos, 1981; Valtaoja, Valtonen & Byrd, 1989). This scenario has received considerable attention because the ultimate coalescence of such a pair would generate an observable outburst of gravity waves (Thorne & Braginskii, 1976). Begelman, Blandford & Rees (1980) showed that the evolution of a MBH binary can be divided into three distinct phases: 1. As the galaxies merge, the MBHs sink toward the center of the new galaxy via dynamical friction (Chandrasekhar, 1943) where they form a binary. 2. The binary continues to decay via gravitational slingshot interactions (Saslaw, Valtonen & Aarseth, 1974) in which stars on orbits intersecting the binary’s are ejected at velocities comparable to the binary’s orbital velocity, while the binary’s binding energy increases. 3. If the binary’s separation decreases to the point where the emission of gravity waves becomes efficient at carrying away the last remaining angular momentum, the MBHs coalesce rapidly. In this paper we explore Phase 2 and its transition to Phase 3; this transition is understood to be the bottleneck of a MBH binary’s path to coalescence.

The picture of Begelman, Blandford & Rees (1980) has remained essentially valid in the face of several important developments in intervening years. Early theoretical studies envisioned the centers of galaxies as constant-density cores, with small, embedded stellar cusps that develop around central MBHs via collisional relaxation (Peebles, 1972; Bahcall & Wolf, 1976; Young, 1977; Ipser, 1978). Observations at the time lacked the resolution to verify this picture. In the past decade, space-based observations have revealed that, in the majority of early-type galaxies, the central luminosity density increases continually inward as an approximate power law $`\rho (r)r^\gamma `$ down to regions where the MBH dominates the gravitational force (Ferrarese et al., 1994; Gebhardt et al., 1996; Rest et al., 2001). The typical luminosity profile of a faint elliptical galaxy is a nearly-featureless power law with $`1.5\gamma 2.5`$ at the resolution limit. The times required to build cusps via two-body relaxation are longer than the age of the Universe in most nuclei (e.g. Faber et al. 1997) and so relaxation is not a viable explanation for how the cusps formed. In this work we assume that galaxies inherited their steep central density profiles from an early epoch in which the MBH grew. Cold dark matter (CDM) halos produced by hierarchical structure formation that presently host galaxies and MBHs were cuspy at the time of virialization (e.g. Ghigna et al. 2000; Power et al. 2003). On smaller scales, the stellar cusps may have resulted from gaseous dissipation of star-forming medium inside pre-existing CDM cusps, from adiabatic response of the stars to growth of the MBH, or both. The same processes that led to cusp formation were likely to have fueled the early growth of MBHs and, thus, set the MBHs and galaxies on a course of co-evolution. We therefore assume that the loci of MBHs have coincided with the centers of cosmic overdensities (dark matter and stellar cusps) since the time MBHs first ascended into the massive range, $`M_{}10^6M_{}`$.

The cusps that accompany the MBHs endow them with an enlarged effective dynamical mass. Milosavljević & Merritt (2001) found that the dynamical friction time scale (Phase 1) is a function of this effective mass and much shorter than the orbital decay time implied by the masses of the MBHs alone. This is of critical importance for the capacity of MBHs originating in different, merging galaxies to capture each other and form a hard binary in the resulting galaxy. A MBH binary is customarily called “hard” if its binding energy per unit mass exceeds the specific kinetic energy of ambient stars, $`G\mu /4a\sigma ^2`$, where $`\mu =M_1M_2/(M_1+M_2)`$ is the reduced mass, $`a`$ is the binary’s semi-major axis and $`\sigma `$ is the 1D stellar velocity dispersion. Milosavljević & Merritt (2001) also found that in mergers of equal-mass $`\rho (r)r^2`$ galaxies, the true separations at which the binaries settle at the beginning of Phase 2 are $`10`$ times smaller than implied by arguments based on the equipartition of energy.

Masses of MBH’s correlate exceptionally well with the central line-of-sight velocity dispersions of their host galaxies (Ferrarese & Merritt, 2000; Gebhardt et al., 2000). It is likely that the masses are also related to the depth of the potential wells associated with the dark matter halos in which the host galaxies are embedded (Ferrarese, 2002a), which justifies attempts to characterize the evolution of the MBH population using the statistics of hierarchical formation of structure in dark matter cosmologies (Kauffmann & Haehnelt, 2000; Menou, Haiman & Narayanan, 2001; Volonteri, Haardt & Madau, 2002). As a corollary, one can estimate the frequency of MBH binary coalescence at redshifts that are beyond the reach of existing telescopes but within reach of a future space-based gravity wave detector such as LISA.111Laser Interferometer Space Antenna, http://lisa.jpl.nasa.gov The coalescence of two $`10^510^6M_{}`$ MBHs occurring practically anywhere in the observable Universe would generate a significant signal for LISA.

Early calculations (e.g. Haehnelt 1994) identified the MBH coalescence event rate with the galaxy merger rate. It was subsequently noted that stellar-dynamical processes are inefficient at facilitating the decay of MBH binaries when their separations shrink below about one parsec and that the ultimate coalescence should not be taken for granted (Valtonen, 1996a; Gould & Rix, 2000; Milosavljević & Merritt, 2001; Yu, 2002). Circumstantial evidence, however, suggests that long-lived MBH binaries are in fact rare. The coincidence of mean black hole mass densities derived from kinematical studies of local galaxies and distant quasars (Ferrarese (2002b)) suggests that MBHs are rarely ejected from galaxy centers by the strong interactions that would accompany the infall of a third MBH into a nucleus containing an uncoalesced binary. Furthermore there is evidence for MBH coalescences in the morphology of radio galaxies, at a rate approximately equal to the galaxy merger rate (Merritt & Ekers, 2002).

There are many processes that could drive the MBH binaries that form in galaxy mergers to coalescence. In this paper we present some new ideas regarding one such process, the gravitational slingshot ejection of stars.

Ejection of stars from the galaxy nucleus will result in a depletion of stars in the “loss cone.” We define the loss cone as the phase space domain in which orbits approach the binary closely enough to be influenced by the near field of the individual binary components. The orbits in the loss cone are subject to gravitational capture and slingshot ejection. The loss cone has thus traditionally been viewed as a sink in the phase space, surrounded by a distribution of stars in equilibrium with respect to two-body gravitational encounters. The latter assumption allows a steady-state solution for the phase space density near the loss cone boundary to be derived (e.g. Cohn & Kulsrud 1978); this steady state reflects a balance between the supply of stars to the loss cone by gravitational scattering and the loss due to capture or ejection. In this standard picture, stars that scatter into the loss cone are ejected from the system with net energy loss, resulting in a more tightly bound binary. After all of the stars originally inside the loss cone have been removed, continued decay of the binary hinges on the flux of stars scattering into the loss cone.

In this paper we study the physics of the “worst case scenario,” involving nearly spherical galaxies, in which the MBH binaries are least likely to coalesce in their lifetime. In § 2 we review the standard solution for the flux. In § 3 we question the applicability of some of the assumptions that were made to derive that solution. In § 4 and § 5 we introduce some novel aspects of the dynamics of stars in the loss cone and outline an iterative procedure for calculating the long-term evolution of a massive MBH binary in a spherical galaxy. The details emphasized here tend to shorten the timescale for the decay of MBH binaries, sometimes significantly. In § 6 we evaluate the effectiveness of $`N`$-body simulations to model the long-term evolution of MBH binaries. In § 7 we address the implications of the Brownian motion of MBH binaries for the long-term evolution. In § 8 we review the likelihood for the coalescence of MBH binaries. Extensions to axisymmetric and triaxial galaxies, as well as applications to specific galaxy models, are deferred for future work.

## 2. Review of Standard Loss Cone Theory

The theory of loss cone structure was originally developed to model the rates of tidal disruptions of stars by MBHs (Frank & Rees, 1976). Results reviewed here were derived in the context of globular cluster-like stellar systems containing MBHs (Lightman & Shapiro, 1977; Cohn & Kulsrud, 1978; Ipser, 1978) where the stars are tidally disrupted after coming within the tidal radius of the MBH:

$$r_{\mathrm{tidal}}=\left(\eta ^2\frac{M_{}}{m_{}}\right)^{1/3}r_{},$$
(1)
where $`m_{}`$ and $`r_{}`$ are, respectively, the stellar mass and the stellar radius, and $`\eta `$ is a form factor of order unity. Recently, the loss cone paradigm has been utilized to estimate rates of disruptions of stars by MBHs residing in galaxy nuclei (Sigurdsson & Rees, 1997; Syer & Ulmer, 1999; Magorrian & Tremaine, 1999), as well as to calculate the rates of scattering of stars into the capture zone of MBH binaries (Yu, 2002). At the heart of these studies lies the assumption that stellar systems, such as globular clusters or small galaxies, are old compared with local two-body relaxation times. This justifies the description of these systems in terms of the time-independent Fokker-Planck equation (Rosenbluth, MacDonald & Judd, 1957; Cohn & Kulsrud, 1978).

For simplicity we restrict attention to spherical galaxies; there the phase space distribution is a function of two orbital integrals—the energy $`E`$ and the angular momentum $`J`$ of a star. The Boltzmann equation reads:

$$\frac{f}{r}=\frac{1}{v_r}\frac{}{R}\left[\frac{(\mathrm{\Delta }R)^2}{2R}R\frac{f}{R}\right],$$
(2)
(e.g. Cohn & Kulsrud 1978, eq. 30), where $`f(E,R,r)d^3rd^3v`$ is the probability of finding a star in the phase space volume element $`d^3rd^3v`$, $`v_r`$ is the radial velocity of the star, $`R=J^2/J_c^2(E)`$ is the square angular momentum in units of the angular momentum of a circular orbit at energy $`E`$, and $`(\mathrm{\Delta }R)^2`$ is the diffusion coefficient associated with the orbital integral $`R`$. At no risk of confusion we will variably refer to $`R`$ as “angular momentum.” Note that all terms representing diffusion in energy have been discarded; this is justified because the gradients in $`E`$ are minuscule compared with the gradients in $`R`$ in the vicinity of the loss-cone boundary. The diffusion coefficient scales inversely with the number of stars in the galaxy:

$$(\mathrm{\Delta }R)^2N^1\frac{m_{}}{M_{\mathrm{galaxy}}}.$$
(3)
Conditions for the validity of the time-independent, orbit-averaged Fokker-Planck equation can be summarized as follows:

$$P(E)T_{\mathrm{relax}}(E)T_{\mathrm{age}},$$
(4)
where $`P(E)`$ is a typical orbital period (or the crossing time) of a star at energy $`E`$, $`T_{\mathrm{relax}}`$ is the local two-body relaxation time, and $`T_{\mathrm{age}}`$ is the total lifetime of the system. If the first inequality fails, perturbative expansion of the Boltzmann equation is no longer valid, large deflections dominate, and the system is in a gaseous phase. If the second inequality fails, the system need not be collisionally relaxed and can keep evolving; then the time-dependent Fokker-Planck equation is necessary. We argue in the next section that the second inequality is almost never satisfied in galaxies with MBH binaries.

Let $`r_{\mathrm{bin}}`$ denote the distance from the binary’s center of mass within which stellar orbits are strongly perturbed by the rotating quadrupole gravitational field of the binary. This distance is similar to the semi-major axis of the binary, $`r_{\mathrm{bin}}a`$. Stars that transgress a region around the binary $`rr_{\mathrm{bin}}`$ are subject to gravitational slingshot and are ejected from the system. This imposes a boundary condition on equation (2) of the form:

$$f(E,R,r)=0(rr_{\mathrm{bin}}).$$
(5)
Lightman & Shapiro (1977) derived an orbit averaged solution $`f(E,R)`$ to equation (2) subject to the boundary condition expressed by equation (5) to find:

$$f(E,R)P^1(E)\frac{dr}{v_r}f(E,R,r)=C\mathrm{ln}\left(\frac{R}{R_0}\right),$$
(6)
where $`C`$ is a normalization factor to be determined, and $`R_0(E)`$ is an energy-dependent angular momentum cutoff.

Equation (6) assumes that $`f(E,R,r)`$ vanishes identically when $`RR_0`$. In the absence of two body relaxation, $`R_0=R_{\mathrm{lc}}`$, where $`R_{\mathrm{lc}}(E)`$ is the angular momentum of a particle with pericenter distance at the very edge of the sphere of capture, $`r_{}(E,R_{\mathrm{lc}})=r_{\mathrm{bin}}`$. In the presence of relaxation, some stars will enter and exit the consumption zone $`R<R_{\mathrm{lc}}`$ all within one orbital period, thereby evading consumption. Cohn & Kulsrud (1978) carried out a boundary-layer analysis to find the following relation between $`R_0`$ and $`R_{\mathrm{lc}}`$:

$$R_0(E)=R_{\mathrm{lc}}(E)\times \{\begin{array}{cc}\mathrm{exp}(q),\hfill & q(E)>1\hfill \\ \mathrm{exp}(0.186q0.824\sqrt{q}),\hfill & q(E)<1\hfill \end{array},$$
(7)
where $`q(E)`$ is the ratio of the orbital period at energy $`E`$ to the timescale for diffusional refilling of the consumption zone at this energy:

$$q(E)\frac{1}{R_{\mathrm{lc}}(E)}\frac{dr}{v_r}\underset{R0}{lim}\frac{(\mathrm{\Delta }R)^2}{2R}.$$
(8)
When the period is much shorter than the refilling time scale ($`q1`$), the system is in a “diffusive” regime, the loss cone is largely empty and $`R_0R_{\mathrm{lc}}`$. When the period is much longer than the refilling time scale ($`q1`$), the system is in a “pinhole” regime, the loss cone is largely full and $`R_0R_{\mathrm{lc}}`$. The energy at which $`q=1`$ is called the critical energy. If the loss cone is as narrow as in the tidal disruption problem, parts of the phase space at the largest energies $`E2\sigma ^2`$ will be in the diffusive regime, while the less bound regions will be in the pinhole regime. Most of the lost cone flux originates near the critical energy.

The normalization factor $`C`$ can be expressed in terms of the isotropized distribution function $`\overline{f}(E)=_0^1f(E,R)𝑑R`$ to obtain:

$$C(E)=\left(R_0\mathrm{ln}R_01\right)^1\overline{f}(E),$$
(9)
implying that the orbit-averaged solution for the time-independent Fokker-Planck equation reads (using $`R_01`$):

$$f(E,R)=\frac{\mathrm{ln}(R/R_0)}{\mathrm{ln}(1/R_0)1}\overline{f}(E).$$
(10)
In this case, the number flux of stars into the loss cone is given by:

$`(E)dE`$ $`=`$ $`4\pi ^2J_c^2\left\{{\displaystyle \frac{dr}{v_r}\underset{R0}{lim}\frac{(\mathrm{\Delta }R)^2}{2R}}\right\}`$ (11)
$`\times {\displaystyle \frac{\overline{f}}{\mathrm{ln}(1/R_0)1}}dE,`$
where all quantities on the right hand side depend on energy, and $`R_0`$ also depends on the geometric size of the BH binary, $`R_0G(M_1+M_2)a/J_c^2`$. Note that the product $`(\mathrm{\Delta }R)^2\overline{f}`$ is independent of the number of stars in the galaxy $`NM_{\mathrm{galaxy}}/m_{}`$. The total mass flow $`m_{}(E)𝑑E`$, however, scales as $`N^1`$.

## 3. Does the standard theory apply to MBH binaries?

The orbital separation $`a`$ of MBHs in a binary is much larger than the tidal disruption radius around an isolated MBH. Using equation (1):

$`{\displaystyle \frac{r_{\mathrm{bin}}}{r_{\mathrm{tidal}}}}`$ $``$ $`10^5\times \eta ^{2/3}\left({\displaystyle \frac{M_{}}{10^8M_{}}}\right)^{1/3}\left({\displaystyle \frac{m_{}}{M_{}}}\right)^{1/3}`$ (12)
$`\times \left({\displaystyle \frac{r_{}}{R_{}}}\right)^1\left({\displaystyle \frac{a}{1\text{ pc}}}\right).`$
Henceforth, $`M_{}=M_1+M_2`$ denotes the total mass in MBHs. In the context of a loss cone associated with a binary, most of the galaxy is in the strongly diffusive regime and the boundary layer is of negligible thickness, $`R_0R_{\mathrm{lc}}`$. Figure 3a shows the dependence of the quantity $`q`$ on energy for a model described by the $`\rho r^2(1+r^2)^2`$ stellar density profile and scaled to the galaxy M32 with $`0.1\%`$ of the galaxy mass in the central MBH. Note that $`q(E)1`$ when the MBHs first form a bound pair and remains $`1`$ even when the binary separation shrinks to $`0.01`$ pc. Scaling to other galaxies, assuming a nuclear density profile of $`\rho r^2`$ and a fixed ratio of MBH mass to the galaxy mass, we find:

$`q(E)`$ $``$ $`0.025\times \left({\displaystyle \frac{m_{}}{M_{}}}\right)\left({\displaystyle \frac{M_{}}{3\times 10^6M_{}}}\right)^1`$ (13)
$`\times \left({\displaystyle \frac{a}{GM_{}/8\sigma ^2}}\right)^1e^{E/2\sigma ^2}`$
with $`\sigma `$ the 1D stellar velocity dispersion. Since $`q(E)P(E)/T_{\mathrm{relax}}(E)`$, the typical deflection of a star due to the presence of other stars over one orbital period is small compared with the size of the loss cone. Therefore, $`f`$ is nearly uniform around the orbit of the star; the stellar system satisfies Liouville’s theorem; and the orbit-averaged Fokker-Planck equation (Lightman & Shapiro, 1977) is a valid description of the loss cone dynamics. We will rely on this conclusion throughout this work.

Another aspect of the standard theory fares less well in the case of MBH binaries. We have so far discussed the time-independent Fokker-Planck equation; strictly speaking, however, time-independent solutions do not exist as stars that diffuse into the loss cone are moved from one orbit to another, or removed from the system. The flux of stars into the loss cone given by the standard theory is shown in Figure 3b. As the stellar orbital integrals diffuse, the potential changes accordingly, and the isotropic distribution $`\overline{f}(E)`$ adjusts to the changing potential. Milosavljević & Merritt (2001) demonstrated that the change of the overall density profile of the galaxy can amount to a destruction of the steep central cusp. One may attempt to account for the changing galaxy profile by adjusting $`\overline{f}(E)`$ as dictated by the isotropic Jeans equation, while keeping the formal solution in equation (10) unchanged. Similarly, as the MBH binary hardens, the loss cone boundary $`R_{\mathrm{lc}}`$ decreases. One could try scaling the formal solution to accommodate a changing $`R_{\mathrm{lc}}`$.

This strategy meets a serious obstacle when considering the time required for the product of a galaxy merger to converge to the quasi-steady state described in the previous section. The merger and the subsequent formation of a hard binary proceed on the local dynamical time scale, which can be much shorter than the collisional relaxation time scale. Therefore the distribution function $`f(E,R)`$ immediately following the formation of a hard binary can be far from the solution in equation (10). The mass in stars near the loss cone in excess of the equilibrium solution can be orders of magnitute larger than the mass fed into the loss cone over a Hubble time assuming the equilibrium solution. Sudden draining of the loss cone during formation of the hard binary produces steep phase space gradients that are closer to the step function:

$$f(E,R)\{\begin{array}{cc}\overline{f}(E),\hfill & R>R_{\mathrm{lc}}\hfill \\ 0,\hfill & R<R_{\mathrm{lc}}\hfill \end{array}.$$
(14)
Since the collisional transport rate in phase space is proportional to the gradient of $`f`$ with respect to $`R`$, steep gradients imply an enhanced flux into the loss cone and an accelerated evolution toward the equilibrium form.

Figure 3c shows that in low mass ellipticals like M32, the time needed to achieve equilibrium depends on the energy, but still measures well in excess of the Hubble time for the relevant range of energies. In intermedate-mass and massive galaxies, the time is longer than the Hubble time regardless of the orbital energy. In § 4 we derive equations describing the evolution of the stellar distribution starting from arbitrary post-merger initial conditions and quantify the discrepancies with the steady state solution. We find that the steady-state solution is never a good description of a stellar system containing a binary MBH.

We have so far ignored the fate of the stars that are ejected by the binary. Unlike the case of tidal disruption of stars by a MBH, a star ejected by a binary MBH via the gravitational slingshot remains inside the galaxy and its host dark matter halo. Ejected stars may therefore return to interact with the binary for a second and subsequent times. Every ejected star can, in principle, interact with the binary infinitely many times. Its energy can increase or decrease at every encounter. This process has been ignored in the published Fokker-Planck treatments of tidal disruption; nonetheless it warrants consideration because of its capacity to shorten the lifetime of a MBH binary even after all of the stars inside the loss cone have been ejected at least once. We discuss the efficiency of reejection in § 5.

The standard treatment predicts that the decay time of MBH binaries scales with the inverse stellar mass $`a/\dot{a}m_{}^1N`$, where $`N`$ is the number of stars in the galaxy. Evidence that self-consistent $`N`$-body simulations reproduce this trend has so far been equivocal. Simulations with $`M_{}/m_{}10^3`$ ($`N10^5`$), such as those reported in Milosavljević & Merritt (2001), appear to have taken place in the “pinhole” regime, where the loss cone is almost full and the $`N`$-dependence is weak. In simulations where the MBH binary was allowed to wander in space, its Brownian motion was a likely source of additional, spurious diffusion of stars in the phase space. Simulations with $`M_{}/m_{}10^4`$ and a binary that does not wander are able to recover the trends characteristic of the diffusive regime. We discuss $`N`$-body simulations of MBH binary dynamics in § 6 and § 7.

## 4. Time-Dependent Evolution of the Loss Cone

### 4.1. The Loss Cone as an Initial Value Problem

In this section we study the diffusion of stars into the loss cone of a MBH binary following a merger of two galaxies. We assume that stars are removed from the final galaxy the first time they are ejected by the MBH binay; the continued role of these stars is addressed in § 5. The loss cones of galaxies harboring MBHs are in the diffusive regime. A large fraction of stars scattering into the geometric loss cone are kicked out by the binary in one orbital period. This suggests that $`R=R_{\mathrm{lc}}(E)`$ can be thought of as an absorbing boundary condition where the distribution function vanishes. Since $`R_{\mathrm{lc}}(E)`$ has no explicit dependence on the orbital phase, one can seek solutions of the Fokker-Planck equation that have no explicit dependence on $`r`$. The only remaining dependence on $`r`$ enters through the diffusion coefficients in the collision term of the Boltzmann equation and can be integrated out. As a result, one obtains the “orbit averaged” Fokker-Planck equation (Lightman & Shapiro 1977).

Instead of the distribution function $`f`$, we work in terms of the number density of stars in the $`(E,R)`$ plane, which is related to the distribution function via:

$$𝒩(E,R,t)dEdR=4\pi ^2P(E,R)J_c(E)^2f(E,R,t)dEdR.$$
(15)
The time-dependent Fokker-Planck equation including terms that describe diffusion in the angular momentum direction reads:

$$\frac{𝒩}{t}=\frac{}{R}\left(\mathrm{\Delta }R𝒩\right)+\frac{1}{2}\frac{^2}{R^2}\left[(\mathrm{\Delta }R)^2𝒩\right].$$
(16)
Using a relation (Binney & Lacey, 1988) between the the first- and the second-order diffusion coefficients,

$$\mathrm{\Delta }R=\frac{1}{2}\frac{}{R}\left(\mathrm{\Delta }R\right)^2,$$
(17)
equation (16) can be simplified to obtain:

$$\frac{𝒩}{t}=\frac{1}{2}\frac{}{R}\left[\left(\mathrm{\Delta }R\right)^2\frac{𝒩}{R}\right].$$
(18)
Next we expand the diffusion coefficient $`(\mathrm{\Delta }R)^2`$ in the limit $`R0`$:

$$\frac{1}{2}(\mathrm{\Delta }R)^2=\left[\underset{R0}{lim}\frac{(\mathrm{\Delta }R)^2}{2R}\right]R+𝒪(R^2),$$
(19)
where the coefficient in brackets depends on energy and radius. We then average over one orbital period $`P^1𝑑r/v_r`$ so that the Fokker-Planck equation ($`R1`$) becomes:

$$\frac{𝒩}{t}=4\mu \frac{}{R}\left(R\frac{𝒩}{R}\right),$$
(20)
where

$$\mu (E)\frac{1}{4P}\frac{dr}{v_r}\underset{R0}{lim}\frac{(\mathrm{\Delta }R)^2}{2R}.$$
(21)
Note that this definition of $`\mu (E)`$ differs from that in Magorrian & Tremaine (1999) by a factor of $`1/4`$. We finally carry out a change of variables $`R=j^2`$ to rewrite the equation as:

$$\frac{𝒩}{t}=\frac{\mu }{j}\frac{}{j}\left(j\frac{𝒩}{j}\right).$$
(22)
This is the heat equation in cylindrical coordinates with radial variable $`j`$ and diffusivity $`\mu `$. At every energy $`E`$, boundary conditions need to be specified at $`0<j_{\mathrm{lc}}<1`$ and $`j=1`$, where $`j_{\mathrm{lc}}\sqrt{R_{\mathrm{lc}}}`$. The boundary condition at $`j=1`$ is of the Neumann type:

$$\frac{𝒩}{j}|_{j=1}=0.$$
(23)
The boundary condition at $`j=j_{\mathrm{lc}}`$ is a perfectly absorbing, Dirichlet boundary condition:

$$𝒩(E,j)=0(jj_{\mathrm{lc}})$$
(24)
The geometrical meaning of the boundary conditions is illustrated in Figure 4.1.

The general solution of equation (22) subject to boundary conditions (23) and (24) can be obtained by means of Fourier-Bessel synthesis (e.g. Özisik 1980):

$`𝒩(E,j,t)`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{2}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left[\beta _mJ_0(\beta _mj_{\mathrm{lc}})\right]^2}{\left[J_0(\beta _mj_{\mathrm{lc}})\right]^2\left[J_1(\beta _m)\right]^2}}`$ (25)
$`\times A(\beta _m;j)e^{\mu \beta _m^2t}`$
$`\times {\displaystyle _{j_{\mathrm{lc}}}^1}j^{}A(\beta _m;j^{})𝒩(E,j^{},0)dj^{},`$
where $`J_n`$ and $`Y_n`$ ($`n=0,1`$) are Bessel functions of the first and second kind, $`A(x;y)`$ is a combination of the Bessel functions defined via:

$$A(x;y)J_0(xy)Y_1(x)J_1(x)Y_0(xy),$$
(26)
while the $`\beta _m`$ are consecutive solutions of the equation:

$$A(\beta _m;j_{\mathrm{lc}})=0.$$
(27)
We also evaluate the time-dependent version of the loss cone flux:

$`(E,t)dE`$ $`=`$ $`{\displaystyle \frac{d}{dt}}\left[2{\displaystyle _{j_{\mathrm{lc}}}^1}𝒩(E,j,t)j𝑑j\right]`$ (28)
$`=`$ $`\pi ^2{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\mu \beta _{m}^{}{}_{}{}^{3}j_{\mathrm{lc}}\left[J_0(\beta _mj_{\mathrm{lc}})\right]^2}{\left[J_0(\beta _mj_{\mathrm{lc}})\right]^2\left[J_1(\beta _m)\right]^2}}`$
$`\times B(\beta _m;j_{\mathrm{lc}})e^{\mu \beta _m^2t}`$
$`\times {\displaystyle _{j_{\mathrm{lc}}}^1}j^{}A(\beta _m;j^{})𝒩(E,j^{},0)dj^{}.`$
where $`B(x;y)`$ is another combination of the Bessel functions

$$B(x;y)J_1(xy)Y_1(x)J_1(x)Y_1(xy).$$
(29)
### 4.2. History of the Loss Cone Profile

In Figure 4.2a we show the evolution of $`𝒩(E,R,t)`$ at one, arbitrarily chosen $`E`$, assuming that the loss cone is empty within $`R_{\mathrm{lc}}`$ and that $`𝒩`$ is a constant function of $`R`$ outside of $`R_{\mathrm{lc}}`$ at the beginning ($`t=0`$). Here we have assumed that the loss cone boundary is static, $`R_{\mathrm{lc}}0.02`$. For comparison, we also show the equilibrium solution of equation (10) which can be expressed:

$$𝒩_{\mathrm{equi}}(E,R,t)=\frac{\mathrm{ln}(R/R_{\mathrm{lc}})}{\mathrm{ln}(1/R_{\mathrm{lc}})1}_{R_{\mathrm{lc}}}^1𝒩(E,R,t)𝑑R.$$
(30)
The phase-space gradients $`𝒩/R`$ decay rapidly at first and then more gradually as they approach the equilibrium solution; the former converges to, but never becomes equal to the latter. It is evident that the total population $`𝒩𝑑R`$ incurs a decrement of order unity before it has had time to reach the state of collisional equilibrium. A decrease in $`𝒩`$ implies a decrease in the spatial density gradients in the galaxy, or “cusp destruction.”

In Figure 4.2b we present the total mass consumed by the binary, which is given by equation (32). The time-dependent and the equilibrium masses are shown for comparison. The equilibrium solution (dot-dashed curve) yields a loss cone flux smaller by at least a factor of two within a Hubble time (Fig. 4.2b). The static loss cone boundary approximation employed here is valid when $`M_{\mathrm{lost}}M_{}`$ (§ 4.3), or below the dashed line in the figure; we consider the general case in § 4.4. Our model here is that of a galaxy with a steep density profile such as M32. The enhancement over the equilibrium fluxes is the strongest in the galaxies with shallower central profiles or “cores;” there, however, $`M_{\mathrm{lost}}`$ remains smaller than $`M_{}`$ during a Hubble time assuming perfect sphericity of the galactic potential.

The following remarks illustrate how the time dependence of the loss cone profile complicates the inferences that can be made about MBH binaries in real galaxies. Given $`\rho (r)`$, the Jeans equation can be solved for the isotropic distribution function $`f(E)`$ and the orbital population $`𝒩(E)𝒩(E,R)𝑑R`$. However from equation (20), the flux of stars scattering into the loss cone is proportional to the phase space gradient in $`R`$ at the boundary of the loss region, or:

$$(E,t)=4\mu R_{\mathrm{lc}}\frac{𝒩}{R}|_{R_{\mathrm{lc}}}.$$
(31)
If the nucleus is not collisionally relaxed, the angular momentum dependence of $`𝒩(E,R)`$ is unknown. Since the gradient $`𝒩/R`$ evolves in time (e.g. equation 28), the loss cone flux is likely to be a strong function of the age of the galaxy since the most recent merger event, as well as of other details like the mass ratio of the merger.

Attempts to estimate the separations of MBH binaries in presently observed galaxies (e.g. Yu 2002) are therefore problematic. As a first step, we would need to relate the total mass in stars ejected by the MBH binary (following the initial emptying of the loss cone) to the evolution of its semi-major axis. The most commonly used form of this relation (Quinlan, 1996), derived via scattering experiments and confirmed in $`N`$-body simulations in the case of equal-mass binaries (Milosavljević & Merritt, 2001), reads:

$`M_{\mathrm{lost}}(t)`$ $``$ $`m_{}{\displaystyle _0^t}{\displaystyle (E,t^{})𝑑E𝑑t^{}}`$ (32)
$``$ $`𝒥M_{}\mathrm{ln}{\displaystyle \frac{a(0)}{a(t)}},`$
where, again, $`M_{}=M_1+M_2`$ is the binary mass, and $`𝒥`$ is a numerical factor approximately equal to one when the binary is hard and $`M_1M_2`$. Milosavljević & Merritt (2001) found that $`𝒥0.5`$ correctly described the binary decay in their $`N`$-body simulations; in a general case $`𝒥`$ is a function both of the hardness of the binary and the MBH mass ratio $`M_1/M_2`$ (but see § 5.4). Therefore the number of $`e`$-foldings in the binary separation is proportional to the mass ejected. Since the evolution of the binary separation depends on the history of the loss cone flux from the binary formation onward, the assumption that the equilibrium flux obtained from the presently observed $`\rho (r)`$ was the same at all times in the past (Yu, 2002) is too conservative. Not only would the phase-space gradients $`𝒩/R`$ have been steeper in the past, but the density profile would have been steeper as well, providing significantly enhanced fluxes early in the life of the binary. Present-day MBH binary separations calculated assuming steady-state galaxies (Yu, 2002) should therefore be interpreted as upper limits.

### 4.3. Approximations

Several approximations have been made in deriving the solution in equations (25) and (28). These approximations relate to the assumption that the isotropized distribution function $`f(E)`$ changes very little over the time interval $`t[0,t_{\mathrm{max}}]`$ when the solution is expected to be valid. We summarize the approximations as follows.

1. The local crossing time is much shorter than the local relaxation time, $`P(E)T_{\mathrm{relax}}(E)`$. This is generally required of the perturbative expansion on which the derivation of the Fokker-Planck is based. Since we have calculated the time dependence explicitly, it is no longer necessary that the system be old compared with the relaxation time. It is however imperative that the physical system modeled by equation (22) have reached dynamical equilibrium prior to the time the model can be considered applicable.

2. The shallowing of the overall potential due to the ejection of stars from the loss cone is small. This guarantees that the diffusivity $`\mu (E)`$ also stays constant. In the tidal disruption problem this is always the case, as the consumed stellar mass is at most a small fraction of the MBH mass. In the massive binary problem, the potential is never constant, especially following mergers of galaxies with steep cusps when the binaries impart significant damage to the merged cusp. Indeed, if the ejected mass becomes comparable to the MBH mass—which is a likely outcome of mergers of galaxies with steep density cusps—the potential near the MBH will incur a relative change of order unity. This notwithstanding, if the potential changes slowly enough, adiabatic invariance suggests that the solution of equation (25) could be very close to the true solution modulo an appropriate reparametrization of the energy dependence.

3. The geometric loss cone boundary $`R_{\mathrm{lc}}(E)=j_{\mathrm{lc}}^{}{}_{}{}^{2}(E)`$ changes little in time. Again, this is true for the tidal disruption loss cone, but almost never true for the binary loss cone. As the binary siphons energy into stars that are captured and then ejected, the semi-major axis $`a`$ of the binary decreases. Since $`R_{\mathrm{lc}}a^2`$, $`j_{\mathrm{lc}}`$ shrinks proportionally to the binary separation. The solution in equation (25) will be valid only so long as $`\mathrm{\Delta }aa`$. The domain of validity of the static loss cone boundary approximation can be estimated by considering $`M_{\mathrm{lost}}(t)`$. Recall that the loss cone boundary scales with the binary separation, $`R_{\mathrm{lc}}a`$. Therefore the approximation in which the geometric loss cone boundary is static is valid only if $`M_{\mathrm{lost}}(t)𝒥M_{}`$.

### 4.4. Evolution of the Semi-Major Axis

A self-consistent solution accounting for the evolution of the loss-cone boundary in response to the binary’s orbital decay can be constructed by iterating over the appropriately chosen time steps. The steps are chosen such that the binary separation changes little within each step, $`\mathrm{\Delta }t𝒥M_{}/m_{}(E,t)𝑑E`$. At the end of each time step, one corrects the binary semi-major axis $`a`$ for the energy that the binary has exchanged with the ejected stars as dictated by equation (32). One then recalculates the Fourier-Bessel coefficients in equation (25) using the new, advanced $`𝒩(E,j,t+\mathrm{\Delta }t_i)`$ and the corrected loss cone boundary $`j_{\mathrm{lc}}(E,t+\mathrm{\Delta }t_i)`$ and thus extends the time-dependent solution beyond $`t+\mathrm{\Delta }t_i`$. Iteration of this kind overcomes the requirement for a static background and presents a computationally feasible numerical procedure for calculating loss cone dynamics over cosmological times.

We carried out this procedure in a model of the galaxy M32 and present the results in Figure 4.4a–c. The character of $`a(t)`$ does not vary with the choice of the initial separation $`a(0)`$ as long as the initial width of loss cone is consistent with the separation chosen; here we have set the geometric loss cone boundary to equal $`R_{\mathrm{lc}}(E)=GM_{}a/J_c^2(E)`$. The initial separation used in the middle panel of Figure 4.4, $`a(0)=0.1`$ pc, is the closest to what we expect to find in the merger of two equal galaxies with $`10^6M_{}`$ MBHs. Note that while $`|da/dt|`$ is bounded in the equilibrium solution (dashed line), in the time-dependent solution (solid line) $`|da/dt|`$ is subtantially larger for small $`t`$ and diverges in the limit $`t0`$ due to the very large initial gradients $`𝒩/R`$.

We can express the difference using the quantity (denoting the equilibrium solution with $`a_{\mathrm{equi}}`$):

$$\mathrm{\Theta }(t)\frac{a(0)a(t)}{a(0)a_{\mathrm{equi}}(t)}.$$
(33)
In Figure 4.4d–f we show the evolution of $`\mathrm{\Theta }(t)`$; the enhancement is well described by:

$`\mathrm{\Theta }(t)(1015)\times \left({\displaystyle \frac{3\times 10^6M_{}}{M_{}}}{\displaystyle \frac{m_{}}{M_{}}}{\displaystyle \frac{t}{10\text{ Myr}}}\right)^{(0.330.39)}.`$ (34)
This relation is a good fit for $`tt_{\mathrm{destr}}`$ where $`t_{\mathrm{destr}}`$ is the time scale on which the removal of stars from the stellar cusp (the cusp destruction) starts affecting the pool of stars available for diffusion into the loss cone. The time-scale can be identified with the diffusion of the stellar mass of about $`𝒥M_{}0.5M_{}`$ into the loss cone, which in turn corresponds to the change of the binary separation by the factor $`e^{1/2}0.6`$. From Figure 4.4a-b, $`a/a(0)=0.6`$ at approximately at 1 Gyr for a $`3\times 10^6M_{}`$ MBH. Scaling to an arbitrary binary mass gives:

$$t_{\mathrm{destr}}\frac{M_{}}{3\times 10^6M_{}}\left(\frac{m_{}}{M_{}}\right)^1\times 1\text{ Gyr}.$$
(35)
The difference between the two solutions decreases with time, to roughly a factor of two after a Hubble time.

There is, however, a plausible scenario in which the mechanism discussed in this section could lead to a much greater enhancement in the binary’s mean rate of decay. Galactic nuclei may episodically return to a strongly non-equilibrium state due to events like the capture of a dwarf galaxy, infall of a giant molecular cloud, or star formation from ambient gas. A rejuvenated loss cone will once again result in strong phase-space gradients as the binary rapidly ejects stars with $`R<R_{lc}`$. Then the the non-equilibrium enhancement averaged over the entire lifetime of the binary would be higher than suggested by Figure 4.4. For example, if the episodic replenishment of a MBH binary loss cone in a galaxy like M32 occurs every 10, 100, or 1,000 Myrs, equation (34) implies that the decay of the binary in that galaxy will occur, respectively $`\mathrm{\Theta }`$ 10, 5, or 3 times faster than what the equilibrium theory would have predicted.

In conclusion: the equilibrium solution in equation (10)—originally developed to describe the much smaller loss cones representing tidal disruption of stellar mass objects by massive ones—is never a good approximation to the structure of binary MBH nuclei. Instead, an epoch in which the flux of stars into the loss cone is enhanced compared to the equilibrium solution, is succeeded by one in which the cusp is significantly modified due to ejection of stars.

## 5. Secondary Slingshot

### 5.1. General Considerations

Until now we have ignored the possibility that the stars ejected by the binary once are ejected again as they return to the nucleus on radial orbits. The stationary solution for the loss cone dynamics (§ 2) and the time-dependent solution (§ 4) are based on the “sinkhole” paradigm modelled after the tidal disruption of stars by a MBH. In the tidal disruption picture stars are destroyed as soon as they transgress the tidal disruption zone of the MBH. This is not true for MBH binaries. The capture cross section of binaries is larger than that that of single MBHs by factors of $`10^510^{11}`$ and only a negligible number of stars are disrupted; the vast majority simply receive a kick $`(\mathrm{\Delta }E,\mathrm{\Delta }R)`$ that transports them to another orbit. Following an event of slingshot, the final angular momentum is comparable to the binary’s orbital angular momentum, $`J\sqrt{GM_{}a}`$. If the binary orbit decays slower than the orbital period, $`(d\mathrm{ln}a/dt)^1P`$, most stars inside the loss cone remain inside the loss cone, encounter the binary at their next pericenter passage, and are “reejected.” If the decrement in the binary separation is substantial, some stars that are ejected very near the loss cone boundary finish just outside boundary as the boundary shifts inward. These stars are lost at a rate proportional to the rate of decrease of the loss cone angular momentum $`R_{\mathrm{lc}}a`$.

In this section we explore the possibility that the reejection can prolong the decay of a MBH binary beyond the stalling point when the loss cone would have been depleted in the sinkhole paradigm. To simplify the forthcoming derivation, we consider an orbit-averaged system in which probability for slingshot of star of energy $`E`$ located inside the loss cone is uniform in time and totals to unity within one orbital period. Some stars can become quasi-permanently captured by the MBH binary; in reality such captures are possible but statistically negiligible. Furthermore, we assume that the orbital population function $`𝒩(E,R,t)`$ is constant in $`R`$ inside the loss cone. We define the loss-cone orbital population function:

$$𝒩(E,t)_0^{R_{\mathrm{lc}}}𝒩(E,R,t)𝑑R.$$
(36)
To derive the time-dependence of $`𝒩(E,t)`$, we consider the evolution over a small interval $`\mathrm{\Delta }t`$. The evolution has to account for the redistribution of stars due to the kicks imparted by the MBH binary, as well as for the gradual influx of stars due to the two-body relaxation. Therefore we have:

$$𝒩(E,t+\mathrm{\Delta }t)=𝒩(E^{},t)\zeta (E^{},E,\mathrm{\Delta }t)𝑑E^{}+(E,t)\mathrm{\Delta }t.$$
(37)
Here, $`\zeta (E^{},E,\mathrm{\Delta }t)`$ is the transition probability between the energy $`E^{}`$ at time $`t`$ and the energy $`E`$ at time $`t+\mathrm{\Delta }t`$, and $`(E,t)`$ is the relaxation flux of stars as derived in § 2 and § 4. Note that both depend implicitly on the binary’s semi-major axis. The transition probability can be understood as a combination of the kicks from $`E`$ into some other energy, the kicks from another energy $`E^{}`$ into $`E`$, and the attrition related to the decrease in the loss cone size $`R_{\mathrm{lc}}`$:

$`\zeta (E^{},E,\mathrm{\Delta }t)`$ $`=`$ $`\left[1{\displaystyle \frac{\mathrm{\Delta }t}{P(E^{})}}\right]\delta (EE^{})`$ (38)
$`+{\displaystyle \frac{\mathrm{\Delta }t}{P(E^{})}}\zeta _1(E^{},E)`$
$`+\left[{\displaystyle \frac{R_{\mathrm{lc}}(t+\mathrm{\Delta }t)}{R_{\mathrm{lc}}(t)}}1\right].`$
Here, $`\zeta _1(E^{},E)`$ is the probability that a star ejected from energy $`E^{}`$ will end up on an orbit with energy $`E`$. The probability is normalized to unity $`\zeta _1(E^{},E)𝑑E^{}=1`$.

Equation (38) is written under the assumption that the orbital phases of stars at a given energy are randomly distributed. In reality, this may not be the case as stars retain the memory of their initial orbital phases. Variance in the ejection velocities will generate dispersion in the distribution of stellar phases, which will destroy the phase coherence of the stars after they are ejected a few times. We do not believe that an initial non-uniform distribution of the phases will have a significant effect on the decay of black hole binaries.

Passing to the infinitesimal limit in equation (37) and using $`d/dt\mathrm{ln}R_{\mathrm{lc}}(t)=d/dt\mathrm{ln}a(t)`$ we obtain an equation of motion for the orbital population of the loss cone:

$`{\displaystyle \frac{𝒩(E,t)}{t}}`$ $`=`$ $`{\displaystyle \frac{𝒩(E,t)}{P(E)}}+{\displaystyle \frac{𝒩(E^{},t)}{P(E^{})}\zeta _1(E^{},E)𝑑E^{}}`$ (39)
$`+𝒩(E,t){\displaystyle \frac{d\mathrm{ln}a(t)}{dt}}+(E,t).`$
We proceed to derive an equation describing the evolution of the total number of stars inside the loss cone:

$`{\displaystyle \frac{d𝒩(t)}{dt}}`$ $``$ $`{\displaystyle \frac{d}{dt}}{\displaystyle 𝒩(E,t)𝑑E}`$ (40)
$`=`$ $`𝒩(t){\displaystyle \frac{d\mathrm{ln}a(t)}{dt}}+{\displaystyle (E,t)𝑑E},`$
and of the total energy budget of the loss cone:

$`{\displaystyle \frac{d(t)}{dt}}`$ $``$ $`{\displaystyle \frac{d}{dt}}{\displaystyle 𝒩(E,t)E𝑑E}`$ (41)
$`=`$ $`{\displaystyle \frac{𝒩(E,t)}{P(E)}\left[E+\zeta _1(E,E^{})E^{}𝑑E^{}\right]𝑑E}`$
$`+{\displaystyle 𝒩(E,t)\frac{d\mathrm{ln}a(t)}{dt}E𝑑E}`$
$`+{\displaystyle (E,t)E𝑑E}.`$
The physical meaning of the factor in brackets in the last row of equation (41) is immediately apparent—it is the average energy change that a particle originally at energy $`E`$ experiences as a consequence of the gravitational slingshot:

$$\mathrm{\Delta }EE+\zeta _1(E,E^{})E^{}𝑑E^{}.$$
(42)
The energy gained in the first term of equation (41) has to be compensated by the change in the binary’s binding energy

$$\frac{d}{dt}\left(\frac{GM_1M_2}{2a}\right)=m_{}\frac{𝒩(E,t)}{P(E)}\mathrm{\Delta }E𝑑E.$$
(43)
The hardening of a MBH binary coupled to an evolving stellar population inside the loss cone is described by equations (39) and (43). These equations are rendered in a form conducive to numerical integration.

### 5.2. Example: The Singular Isothermal Sphere

To obtain a sense of how the secondary slingshot contributes to the decay of the binary, we consider the following example: the galaxy is a singular isothermal sphere (SIS) with the density law:

$$\rho (r)=\frac{M}{4\pi r_0r^2},$$
(44)
and the potential:

$$\mathrm{\Phi }(r)=\frac{GM}{r_0}\mathrm{ln}\left(\frac{r}{r_0}\right),$$
(45)
where $`M`$ and $`r_0`$ are constants. Note that the potential due to this configuration diverges at $`r\mathrm{}`$; we have been free to fix $`\mathrm{\Phi }(r_0)=0`$.

The SIS is probably a better approximation to the potential of a generic early-type galaxy than it may seem. Indeed, while the luminous matter (the stars) in a galaxy follows the SIS profile only within a fraction of the effective radius, the combined stellar and CDM contributions to the potential appear to approximate the SIS over and beyond the extent of the luminous distribution. For simplicity we also assume that relaxation is absent in the model: $`(E,t)=0`$. The radial period inside the SIS potential is:

$$P(E)=\frac{dr}{v_r(E)}=\frac{\sqrt{\pi }r_0}{\sigma }e^{E/2\sigma ^2},$$
(46)
where $`\sigma \sqrt{GM/2r_0}`$ is the velocity dispersion outside the radius of influence of the MBH. Substituting this into equation (43) gives:

$$\frac{d}{dt}\left(\frac{GM_1M_2}{2a}\right)=\frac{m_{}}{P(0)}𝒩(E,t)e^{E/2\sigma ^2}\mathrm{\Delta }E𝑑E.$$
(47)
Naively it may appear that the integral in equation (47) diverges. This not the case because large (very bound) energies are forbidden inside the loss cone, $`𝒩(E,t)=0`$ for $`E>E_{\mathrm{max}}`$, since the coexistence of a star with large binding energy and a binary is very improbable. The distribution of stars inside the loss cone peaks at the energy corresponding to the slingshot ejection velocity (see Figure 5.2):

$$E_{\mathrm{eject}}\mathrm{\Phi }_{\mathrm{eject}}\frac{1}{2}\frac{G(M_1+M_2)}{a},$$
(48)
where $`\mathrm{\Phi }_{\mathrm{eject}}`$ is the fiducial potential energy of a star at the edge of the MBH binary’s sphere of influence and $`\sqrt{G(M_1+M_2)/a}`$ is the orbital velocity the MBHs. To be able to solve equation (47) analytically, we pretend that all stars inside the loss cone have the same energy equal to $`E_{\mathrm{eject}}`$ at time $`t`$:

$`𝒩(E,t)𝒩(t)\delta (EE_{\mathrm{eject}}).`$ (49)
Therefore:

$$\frac{d}{dt}\left(\frac{GM_1M_2}{2a}\right)\frac{m_{}}{P(0)}𝒩(t)\mathrm{\Delta }E_{E_{\mathrm{eject}}}e^{E_{\mathrm{eject}}/2\sigma ^2}$$
(50)
Equation (50) can be simplified further after analyzing the evolution of the product $`𝒩(t)\mathrm{\Delta }E_{E_{\mathrm{eject}}}`$. Equation (40) implies that in the absence of relaxation, the total number of stars inside the loss cone decays proportionally to the binary separation, $`𝒩(t)a(t)`$. We further argue that $`\mathrm{\Delta }E_{E_{\mathrm{eject}}}`$ exhibits scaling with $`a^1(t)`$, by considering two extreme cases:

1. Just after two MBH form a hard binary, the stars in the inner stellar cusp inside the loss cone interact with the MBH binary for the first time and receive a kick comparable to, if not larger than, than their kinetic energy:

$$\mathrm{\Delta }E\frac{v^2}{2},\left(\frac{v^2}{2}+\mathrm{\Delta }E\right)\frac{1}{2}\frac{G(M_1+M_2)}{a}.$$
(51)
Thus $`\mathrm{\Delta }Ea^1(t)`$ is indeed a good approximation.

2. When the binary is old and evolving slowly, one may wonder whether a steady state ($`\mathrm{\Delta }E=0`$) can be achieved, in which the stars in the loss cone exchange no net energy with the binary.222We thank C. Pryor for raising this question. A steady state might be expected when the majority of stars inside the loss cone have been ejected to large energies at which the energy exchange with the binary is very inefficient. Assuming that we have indeed succeeded preparing such a steady state, we recognize that the second moment $`(\mathrm{\Delta }E)^2`$ must be finite. Defining $`\delta E(\mathrm{\Delta }E)^2^{1/2}>0`$, we consider the radial periods of stars ejected slower than the average ($`E_1=\overline{E}+\delta E`$) and faster than the average ($`E_2=\overline{E}\delta E`$). Since the former return to re-encounter the binary sooner than the latter, $`P(E_1)<P(E_2)`$, at time $`t+P(\overline{E})`$ the distribution of stars impinging on the binary will have acquired a finite first moment:

$$\overline{E}\overline{E}+\left(\delta E\right)^2\frac{d\mathrm{ln}P^1}{dE}|_{\overline{E}}=\overline{E}+\frac{\left(\delta E\right)^2}{2\sigma ^2}.$$
(52)
Thus the average effective energy of stars as experienced by the binary is shifted by $`(\delta E)^2/2\sigma ^2`$, and the binary will respond with a decrease of the semi-major axis by $`\mathrm{\Delta }a=a^2(\delta E)^2𝒩m_{}/GM_1M_2\sigma ^2`$. This in turn will generate a shift in $`E_{\mathrm{eject}}`$ by the amount:

$`E_{\mathrm{eject}}`$ $``$ $`E_{\mathrm{eject}}{\displaystyle \frac{G(M_1+M_2)}{2a^2}}\mathrm{\Delta }a`$ (53)
$`=`$ $`E_{\mathrm{eject}}{\displaystyle \frac{𝒩m_{}}{\mu }}{\displaystyle \frac{(\delta E)^2}{2\sigma ^2}},`$
where $`\mu M_1M_2/(M_1+M_2)`$ is the reduced mass. Finally, the true distribution of stellar energies shifts in response to the incremental binary hardening by the amount $`\mathrm{\Delta }E=\mathrm{\Delta }E_{\mathrm{eject}}`$. Recall that $`𝒩a`$. From the scaling of the restricted three-body problem we know that $`\delta EG(M_1+M_2)/2a`$. Therefore, again:

$$\mathrm{\Delta }E\frac{𝒩(0)m_{}}{\mu }\frac{a(t)}{a(0)}\frac{G^2(M_1+M_2)^2}{8\sigma ^2a^2(t)}a^1(t).$$
(54)
We have shown that when the binary just becomes hard and evolves rapidly, as well as when it is very hard and evolves very gradually, the product $`𝒩(t)\mathrm{\Delta }Ea^1a^1a^0`$ is approximately constant (although the constant of proportionality could be different in the above two regimes). We proceed to integrate equation (50) where we substitute $`E_{\mathrm{eject}}`$ from equation (48):

$$\frac{1}{a(t)}=\frac{1}{a(0)}+\frac{4\sigma ^2}{G(M_1+M_2)}\mathrm{ln}\left[1+\frac{m_{}𝒩\mathrm{\Delta }E}{2\mu \sigma ^2}\frac{t}{P(E_0)}\right]$$
(55)
for $`t>0`$, where $`a(0)`$ is the initial separation and $`E_0`$ is the energy of ejected stars at the outset:

$$E_0=E_{\mathrm{eject}}|_{t=0}=\mathrm{\Phi }_{\mathrm{eject}}\frac{G(M_1+M_2)}{2a(0)}.$$
(56)
Equation (55) describes the evolution of a MBH binary in the SIS potential ignoring relaxation. In the secondary slingshot mechanism, the inverse semi-major axis and, hence, the binding energy, are logarithmic functions time. This result has been derived assuming that relaxation-driven loss-cone diffusion was negligible. The logarithmic behaviour can be contrasted with the effect of diffusion in which the binding energy grows linearly. In a realistic galaxy both mechanisms are at work; the reejection amplifies the effect of diffusion on the hardening rate of the MBH binary.

It is encouraging that the logarithmic behavior seems to hold in the simulations where the interstellar interactions have been replaced by a smooth potential to prevent relaxation. Figure 5.2 shows the evolution of the inverse semi-major axis in a simulation with $`𝒩18,000`$ stars in the loss cone at $`t=0`$. The observed rate of decay $`da^1/d\mathrm{ln}t1.7\times 10^4`$ comes close to the prediction of equation (55): $`4\sigma ^2/G(M_1+M_2)=2\times 10^4`$.

### 5.3. Reejection in Galaxies

An upper limit on the effectiveness of the mechanism presented here can be estimated as follows. In a merger of two galaxies, the MBHs form a hard binary when $`a(0)a_{\mathrm{hard}}=G\mu /4\sigma ^2`$. The total mass of stars inside the loss cone just after the hard binary forms is $`m_{}𝒩(0)10.0\mu `$. We can also assume that $`M_1=M_2`$ and $`\mathrm{\Delta }E\text{few}\times 2\sigma ^2`$. In a Hubble time, $`a^1`$ will have increased by the factor:

$$\frac{a(0)}{a(t_{\mathrm{Hubble}})}1+0.25\mathrm{ln}\left[\text{few}\times 10\times \frac{t_{\mathrm{Hubble}}}{P(E_0)}\right].$$
(57)
For $`P(E_0)=10^{3,5,7}`$ years, we obtain $`a(0)/a(t_{\mathrm{Hubble}})5`$, $`4`$, and $`3`$, respectively.

This conclusion may be optimistic. In a circular binary, orbital velocities of the MBHs can be written as:

$$V_1=2\sigma \sqrt{\frac{a_{\mathrm{hard}}}{a}\frac{M_2}{M_1}},V_2=2\sigma \sqrt{\frac{a_{\mathrm{hard}}}{a}\frac{M_1}{M_2}}$$
(58)
In the isothermal sphere potential, the stars ejected with velocities $`V_1vV_2`$ will leave the binary’s sphere of influence $`r_\mathrm{h}=GM_{}/2\sigma ^2`$ and travel to a radius $`r_{\mathrm{max}}`$:

$$\mathrm{exp}\left(\frac{M_2}{M_1}\frac{a_{\mathrm{hard}}}{a}\right)\frac{r_{\mathrm{max}}}{r_\mathrm{h}}\mathrm{exp}\left(\frac{M_1}{M_2}\frac{a_{\mathrm{hard}}}{a}\right).$$
(59)
Realistically, reejection becomes ineffective if $`r_{\mathrm{max}}`$ is too large: first because the isothermal sphere potential may not extend indefinitely; and second because a star with large $`r_{\mathrm{max}}`$ is easily perturbed from its nearly-radial orbit on the way in or out. If we suppose that $`r_{\mathrm{max}}100r_\mathrm{h}`$ is an effective upper limit, equation (59) suggests that $`a_{\mathrm{min}}0.2a_{\mathrm{hard}}`$ is an effective lower limit for reejection when $`M_1M_2`$.

However, we note that many stars are ejected with velocities much less than the binary’s orbital velocity (the velocity change varies between 0 and $`\text{few}\times V_2`$, and can also be negative, depending on the star’s orbital parameters and the binary’s phase). This is particularly true when $`M_2M_1`$; the average velocity change of a star encountering the binary is $`\mathrm{\Delta }vV_2`$, although every star inside the loss cone can be pumped to an ejection at $`V_2`$ after some number of encounters with the binary. If $`\mathrm{\Delta }vV_1`$, which is statistically favored when $`M_2M_1`$, the lower limit on the semi-major axis is weakened, $`a_{\mathrm{min}}0.2a_{\mathrm{hard}}`$. We observe that the transfer of energy from the binary to the stars is an incremental, stochastic process involving multiple encounters; the formalism developed in this section is applicable whenever the form of the galaxy’s gravitational potential permits repeated stellar encounters with the MBH binary.

### 5.4. A New Mass-Ejection Formula

The reejection paradigm suggests a new and more general expression for the relation between the ejected mass and the shrinkage factor of a binary. The standard expression, equation (32), was motivated by scattering experiments in which stars are assumed lost if they exit the binary’s sphere of influence with enough velocity to escape the binary. In the case of a binary embedded in a galactic potential, the critical quantity is the energy gained by a star between the time it enters and exits the loss cone. Equating this with the change in the binary’s absolute binding energy, we find:

$$\frac{G(M_1+M_2)\mu }{2}\left(\frac{1}{a_{\mathrm{final}}}\frac{1}{a_{\mathrm{initial}}}\right)=M_{\mathrm{lost}}(\overline{E}_{\mathrm{enter}}\overline{E}_{\mathrm{exit}}).$$
(60)
In an SIS potential, the greatest contribution to the loss cone comes from stars near $`r_\mathrm{h}`$. Defining $`(\overline{E}_{\mathrm{enter}}\overline{E}_{\mathrm{exit}})\mathrm{\Delta }\mathrm{\Phi }`$ and passing to the infinitesimal limit $`a_{\mathrm{initial}}a_{\mathrm{final}}a`$ in equation (60), we obtain:

$$\frac{1}{(M_1+M_2)}\frac{dM_{\mathrm{lost}}}{d\mathrm{ln}a^1}=\frac{1}{(\mathrm{\Delta }\mathrm{\Phi }/2\sigma ^2)(a/a_{\mathrm{hard}})}$$
(61)
Comparing this relation to equation (32) we identify the effective value of the mass-ejection parameter $`𝒥`$:

$$𝒥\frac{1}{(\mathrm{\Delta }\mathrm{\Phi }/2\sigma ^2)(a/a_{\mathrm{hard}})}.$$
(62)
To shrink the binary by one $`e`$-folding, a stellar mass of $`𝒥M_{}`$ must be transported from an energy marginally bound to the black hole, to the galactic escape velocity.

### 5.5. Reejection in Other Geometries

Reejection would occur differently in nonspherical (axisymmetric or triaxial) potentials, since angular momentum would not be conserved and ejected stars could miss the binary on their second passage. However there will generally exist a subset of orbits defined by a maximum pericenter distance that is $`a`$ and stars on such orbits will interact with the binary once per crossing time, just as in the spherical case. As an example, consider centrophilic (chaotic) orbits in a triaxial nucleus (Poon & Merritt, 2002). There is one such orbit per energy (in a time-averaged sense), and the cumulative distribution of pericenter distances for a single orbit is found to be linear in $`r_\mathrm{p}`$ up to a maximum value, $`r_{\mathrm{p},\mathrm{max}}(E)`$ (Merritt & Poon, 2003). Thus the star passes within a distance $`r_{\mathrm{p},\mathrm{max}}`$ of the center each crossing time. (Merritt & Poon, 2003) evaluated $`r_{\mathrm{p},\mathrm{max}}(E)`$ for chaotic orbits in the self-consistent triaxial models of Poon & Merritt (2002), which have density profiles $`\rho r^2`$, and found:

$$r_{\mathrm{p},\mathrm{max}}(E)0.3r_\mathrm{h}e^{\left(\mathrm{\Phi }_\mathrm{h}E\right)/2\sigma ^2}$$
(63)
with $`\mathrm{\Phi }_\mathrm{h}\mathrm{\Phi }(r_\mathrm{h})`$ and $`r_\mathrm{h}=GM_{}/2\sigma ^2`$. The probability of a star passing within $`a`$ during a single traversal of the galaxy is $`a/r_{\mathrm{p},\mathrm{max}}`$ for $`a<r_{\mathrm{p},\mathrm{max}}`$ and one for $`a>r_{\mathrm{p},\mathrm{max}}`$. It follows that all stars on chaotic orbits with $`E\mathrm{\Phi }_\mathrm{h}`$ pass inside of $`a_{\mathrm{hard}}`$ on each orbit, and a substantial fraction of stars would continue to engage in reejection even after the binary had shrunk by a factor of several below $`a_{\mathrm{hard}}`$. Hence reejection in the triaxial geometry might occur in roughly the same way we have described for spherical systems. The situation would be more complex in axisymmetric (nonspherical) galaxies but there would always exist some subset of orbits defined by $`r_{\mathrm{p},\mathrm{max}}<a`$ particularly if the deviations from spherical symmetry were modest.

## 6. $`N`$-Body Simulations and Real Galaxies

In this section, we ask whether the long-term evolution of MBH binaries seen in $`N`$-body simulations is consistent with the theory derived above, and whether the simulations are capable of reproducing the evolution of binaries in real galaxies.

Does the standard theory (§ 2) reproduce the correct time-dependence for the evolution of the binary semi-major axis? To solve for $`a(t)`$, we start from the differential form of equation (32):

$$\frac{d}{dt}\mathrm{ln}\left(\frac{1}{a}\right)=\frac{1}{𝒥M_{}}\frac{dM_{\mathrm{lost}}}{dt}=\frac{m_{}}{𝒥M_{}}(E,a,t)𝑑E.$$
(64)
Note that both $`𝒥`$ and $``$ are functions of $`a`$. In the equilibrium theory (equation 11), $`\mathrm{ln}^1(a_1/a)`$ where $`a_1=J_c^2(E)/GM_{}e`$ is independent of $`a`$. Near the radius of influence of the MBH where the peak loss cone flux originates, $`a_1GM_{}/2\sigma ^2`$. Figure 5a of Quinlan (1996) implies that for hard binaries,

$$𝒥0.1\times \mathrm{ln}\left(\frac{4a_{\mathrm{hard}}}{a}\right),$$
(65)
where as before $`a_{\mathrm{hard}}=G\mu /4\sigma ^2`$. Therefore the evolution of the semi-major axis in the standard model has the approximate form:

$$\frac{d}{dt}\mathrm{ln}\left(\frac{1}{a}\right)\frac{\mathrm{\Gamma }/N_{}}{\mathrm{ln}(a_1/a)\mathrm{ln}(a_2/a)}.$$
(66)
Here, $`\mathrm{\Gamma }/N_{}`$ is a constant proportional to the angular momentum diffusion coefficient. We have also explicitly factored out the $`N`$-dependence of this coefficient via $`N_{}M_{}/m_{}`$, the number of stars whose mass equals that of the MBH. This equation can be easily solved assuming $`a_1a_24a_{\mathrm{hard}}`$, valid for approximately equal-mass MBHs ($`M_{}=4\mu `$), yielding:

$$\frac{a_{\mathrm{hard}}}{a(t)}=\frac{1}{4}e^{(3\mathrm{\Gamma }t/N_{})^{1/3}}.$$
(67)
For $`1<\mathrm{\Gamma }t/N_{}<10`$, the solution (67) is accurately approximated by the linear function:

$$\frac{a_{\mathrm{hard}}}{a(t)}\frac{1}{2}\left(1+\frac{\mathrm{\Gamma }t}{N_{}}\right)(1<\mathrm{\Gamma }t/N_{}<10).$$
(68)
For larger values of $`\mathrm{\Gamma }t/N_{}`$, we are in the regime where the density profile of the galaxy incurs significant damage and many assumptions of the standard model do not apply.

In $`N`$-body experiments with $`N10^5`$ and $`N_{}10^3`$ (Milosavljević & Merritt, 2001) we observed a linear growth of $`1/a(t)`$. However, the binary decay in our simulations went beyond the range $`a(0)/a(t)1`$ where the assumptions inherent in the standard model are satisfied. Additional (unpublished) simulations with $`N10^6`$ also produced a linear growth of $`1/a(t)`$.

Strikingly, however, these simulations did not yield the linear dependence of $`1/a(t)`$ on $`N_{}`$ predicted by equation (68) (see Milosavljević & Merritt (2001), Fig. 8b). Fearing that the Brownian motion of the MBH binary might have contributed to the loss-cone diffusion, we repeated the simulations in a symmetric mode that eliminated the binary’s Brownian motion.333Brownian motion will be identically nil in mergers of equal galaxies that are symmetric with respect to reflection $`(𝐱,𝐯)(𝐱,𝐯)`$. The repeated experiments did yield a decrease in the hardening rate with $`N_{}`$, but it was still much weaker than linear. Similar behavior had been reported by Makino (1997) whose simulated MBHs were more massive relative to the galaxy than ours. Makino observed that the decay timescale was sub-linear in the number of bodies, $`da^1/dtN_{}^{}{}_{}{}^{1/3}`$, which is still strongly at odds with the prediction of equation (68).

The explanation for this discrepancy can be found in Figure 6, showing the quantity $`q(E)`$ in these simulations. Recall that $`q(E)1`$ means that the loss cone is almost empty (“diffusion regime”) while $`q(E)1`$ implies a nearly-full loss cone (“pinhole regime”). The figure shows that the simulations of Milosavljević & Merritt (2001) took place almost entirely in the pinhole regime – the loss cone was nearly full at all times. In § 3 we argued that the MBH binary loss cones in real galaxies are in the diffusive regime throughout the relevant range of energies. In the full loss cone regime, the flux of stars into the loss cone is given by the number of stars inside the loss cone boundary assuming no depletion, divided by the radial period on which these stars arrive to within the binary’s sphere of influence, or:

$$^{(\mathrm{pinhole})}(E)dE\frac{𝒩(E)R_{\mathrm{lc}}(E)}{P(E)}dE$$
(69)
which is valid for $`E<E_{\mathrm{crit}}9\sigma ^2`$ (see Fig. 6), $`q(E_{\mathrm{crit}})1`$. The mass flow into the maw of the binary $`m_{}^{(\mathrm{pinhole})}(E)𝑑E`$ is independent of $`N_{}`$ and always much larger than in the diffusive case. Ignoring reejection and noting that $`R_{\mathrm{lc}}a`$, we find that the semi-major axis is predicted to be an exactly linear function of the time:

$$\frac{1}{a(t)}=\frac{1}{a(0)}+t_0^{E_{\mathrm{crit}}}\frac{Gm_{}𝒩(E)dE}{𝒥P(E)J_c^2(E)},$$
(70)
as observed in the $`N`$-body simulations (Makino, 1997; Milosavljević & Merritt, 2001).

What value of $`N`$ is required in order that a numerical simulation be in the appropriate, diffusive loss-cone regime? The scaling of $`E_{\mathrm{crit}}`$ with $`N`$ can be deduced from equation (13). Assuming a density profile $`\rho r^2`$ and a potential of the form $`2\sigma ^2\mathrm{ln}(r/r_0)`$ such that $`r_0=10^3GM_{}/2\sigma ^2`$ we find:

$$E_{\mathrm{crit}}2\sigma ^2\mathrm{ln}\left(\frac{7.5\times 10^4}{N_{}}\frac{GM_{}/8\sigma ^2}{a}\right).$$
(71)
The transition from a pinhole-like loss cone to a diffusive loss cone occurs when $`N_{}10^410^5`$. Since a typical MBH contains 0.1% of its host galaxy’s mass (Merritt & Ferrarese, 2001), and thus $`N10^3N_{}`$, an $`N`$-body simulation would have to contain $`10^{45}\times 10^3=10^{78}`$ bodies to reproduce the correct, diffusive behavior of a real galaxy. This requirement is a severe one for direct-summation $`N`$-body codes, which rarely exceed particle numbers of $`10^6`$ even on parallel hardware (e.g. Dorband, Hemsendorf & Merritt 2003). One route might be to couple the special purpose GRAPE hardware444http://grape.astron.s.u-tokyo.ac.jp/grape/, which is limited to $`N10^6`$, to algorithms that can handle large particle numbers by swapping with a fast front end.

Finally, we can ask how reejection would affect the predicted time-dependence of the semi-major axis in $`N`$-body simulations, assuming a loss cone always full at energies $`E<E_{\mathrm{crit}}`$. For simplicity we assume that the Brownian motion has been suppressed; the role of Brownian motion at enhancing binary decay is discussed in detail in the next section. Return and reejection of stars is a function of the depth of the potential well in which the model MBH binary has been placed. The appropriate mass-ejection coefficient $`𝒥`$ is the one derived above, equation (62), which takes into account the total energy exchange with a star before it finally escapes from the galaxy. Note that $`\mathrm{\Delta }\mathrm{\Phi }`$—the potential barrier between the energy at which stars are fed into the loss cone, and the energy at which they escape the galaxy—is a function of loss cone entrance energy $`E`$. Then:

$`{\displaystyle \frac{d}{dt}}\mathrm{ln}\left({\displaystyle \frac{1}{a}}\right)`$ $`=`$ $`{\displaystyle \frac{m_{}}{M_{}}}{\displaystyle _0^{E_{\mathrm{crit}}}}{\displaystyle \frac{^{(\mathrm{pinhole})}(E)}{𝒥(E)}}𝑑E`$ (72)
$`=`$ $`{\displaystyle \frac{𝒟}{2}}\left({\displaystyle \frac{a}{a_{\mathrm{hard}}}}\right)^2,`$
where, ignoring the implicit dependence of $`E_{\mathrm{crit}}`$ on $`a`$ (equation 71), $`𝒟`$ is a constant:

$$𝒟=\frac{G^2m_{}\mu }{4\sigma ^4}_0^{E_{\mathrm{crit}}}\frac{𝒩(E)\mathrm{\Delta }\mathrm{\Phi }(E)}{P(E)J_c^2(E)}𝑑E.$$
(73)
The integrand of $`𝒟`$ is independent of $`N_{}`$, however $`E_{\mathrm{crit}}`$ depends weakly on $`N_{}`$. Integration gives:

$$\frac{a_{\mathrm{hard}}}{a(t)}=\sqrt{1+𝒟t}.$$
(74)
This result and expression (55) were both derived taking reejection into account; their difference stems from the assumptions made about the rate at which the loss cone is being refilled. In equation (55), the refilling is absent and the binary decays by reejection only. In equation (74), the refilling is at its maximum (as expected in $`N`$-body simulations) and the reejection serves to amplify the effect of refilling on the binary’s decay, converting a $`\mathrm{log}t`$ dependence to a $`t^{1/2}`$ dependence.

## 7. Brownian Motion

Another potentially important source of flux into the loss cones of simulations in Milosavljević & Merritt (2001) was the Brownian motion of the MBH binary. Single and binary MBHs alike acquire and maintain a random, quasi-periodic drift coming from the discrete encounters with stars at the center of the galaxy. The drift can be understood in terms of the equipartition of kinetic energy between the MBH(s) and the ambient stars. The Brownian motion was recently studied by Merritt (2001), Chatterjee, Hernquist & Loeb (2002) and Dorband, Hemsendorf & Merritt (2003). These studies concluded that the rms excursions of the MBH’s possition and velocity are given by (ignoring factors of order unity):

$$x^2\frac{m_{}}{M_{}}r_{\mathrm{core}}^{}{}_{}{}^{2},\dot{x}^2\frac{m_{}}{M_{}}\sigma ^2,$$
(75)
where $`r_{\mathrm{core}}`$ is the “core radius” of the stellar distribution within which the stellar density approaches a finite, central value. In many real galaxies, however, the core radius is not well-defined as the density increases near the MBH as an approximate power-law, $`\rho r^\gamma `$ with $`0\gamma 2`$.

$`N`$-body simulations (Milosavljević & Merritt, 2001) reveal that the MBH binary remained centered on the stellar density cusp as it executes the Brownian motion; the binary “carries the cusp with it.” The cusp endows the MBHs with a larger effective dynamical mass. We distinguish the stars closely following the MBH (“satellites”), from the stars that orbit within the stellar cusp and come close to the MBH during a fraction of their orbital period (“visitors”). The satellites respond to the binary’s movement and follow the binary as it wanders; their angular momentum with respect to the binary remains approximately constant. The visitors, however, stay embedded in the static galactic potential and experience the binary as moving target. As the binary moves, the visitors drift along nearly constant-energy trajectories in the $`E`$$`R`$ plane that can penetrate the loss cone boundary. The drift can contribute significantly to the loss cone refilling flux at energy $`E`$ if $`\mathrm{\Delta }R_{\mathrm{brown}}R_{\mathrm{lc}}`$, where $`\mathrm{\Delta }R_{\mathrm{brown}}GM_{}r_{\mathrm{brown}}/J_c^2(E)`$.

The meaning of the “core radius” in galaxies without a core in the density profile can be understood in the light of the satellite-visitor dichotomy. The satellites drift along with the MBH binary in the background of the visitors and other stars not interacting with the binary. The satellites, with average velocities larger than the ambient stellar velocity dispersion, form a bound subsystem in combination with the MBHs; the subsystem executes the Brownian motion in the field of other stars. We can therefore think of the background potential in which the binary-satellite subsystem moves as the stellar potential with the contribution of the satellites subtracted. The characteristic radial scale separating the satellites from other stars is the the dynamical radius of influence of the black holes; therefore the effective value for the “core radius” in coreless density profiles such as $`\rho r^2`$ reads:

$$r_{\mathrm{core}}\frac{GM_{}}{\sigma ^2}.$$
(76)
Assume the binary executes a small acceleration $`𝒜`$, which could be the associated with the Brownian motion or some other astrophysical process. We think of $`𝒜`$ as a small perturbation to an otherwise orbital integral-preserving motion of a test star. In the presence of Brownian motion, the acceleration can be expressed using the quantities presented in equation (75):

$$𝒜^2\frac{\dot{x}^2^2}{x^2}\frac{m_{}}{M_{}}\frac{\sigma ^4}{r_{\mathrm{core}}^{}{}_{}{}^{2}}\frac{m_{}\sigma ^8}{G^2M_{}^3}.$$
(77)
Now let $`𝐫`$ and $`𝐉`$ denote, respectively, the position and the angular momentum of the star relative to the binary. Then:

$$\frac{dJ^2}{dt}=2𝐉\frac{d}{dt}\left(𝐫\times 𝐯\right)=2𝐉(𝐫\times 𝒜).$$
(78)
The total change in $`J^2`$ over one orbital period equals:

$$\mathrm{\Delta }J^2=2_0^P𝐉(𝐫\times 𝒜)𝑑t.$$
(79)
Note that $`\mathrm{\Delta }J^2(E,𝐉,𝒜)=\mathrm{\Delta }J^2(E,J,|𝒜|,\widehat{\mathrm{\Omega }})`$ is a function of the angles $`\widehat{\mathrm{\Omega }}`$ defining the inclination of the unperturbed orbit relative to the acceleration vector. To calculate the second-order Fokker-Planck diffusion coefficient for the diffusion in the $`R=J^2/J_c^2(E)`$ (see § 2), we square, divide by the orbital period, and average over $`|𝒜|`$ and $`\widehat{\mathrm{\Omega }}`$:

$$(\mathrm{\Delta }R)^2=\frac{\left[\mathrm{\Delta }J^2(E,J,|𝒜|,\widehat{\mathrm{\Omega }})\right]^2_{|𝒜|,\widehat{\mathrm{\Omega }}}}{P(E)J_c^4(E)}$$
(80)
The diffusion time scale near the loss cone boundary at $`J^2=J_{\mathrm{lc}}^2(E)GM_{}r_{\mathrm{bin}}`$, where $`r_{\mathrm{bin}}`$ is the radius from the center of the binary within which gravitational slingshot can occur (see § 2), will then be:

$`t_{\mathrm{brown}}`$ $`=`$ $`{\displaystyle \frac{R_{\mathrm{lc}}^{}{}_{}{}^{2}}{(\mathrm{\Delta }R)^2_{\mathrm{lc}}}}`$ (81)
$``$ $`{\displaystyle \frac{G^2M_{}^2r_{\mathrm{bin}}^2P(E)}{\left[\mathrm{\Delta }J^2(E,J_{\mathrm{lc}},|𝒜|,\widehat{\mathrm{\Omega }})\right]^2_{|𝒜|,\widehat{\mathrm{\Omega }}}}}`$
To estimate $`t_{\mathrm{brown}}`$ we ignore the stellar self-gravity; the test star is then orbiting in the approximately-Keplerian potential of the MBH binary. Consider an elliptical orbit with the closest approach to the binary at $`r_{}>r_{\mathrm{bin}}`$. Let $`\phi `$ denote the angle between $`𝒜`$ and $`𝐉`$ (see Fig. 7). Since $`dt=(r^2/J)d\theta `$ where $`\theta `$ is the azimuthal angle in the orbital plane, $`\mathrm{\Delta }J^2`$ becomes:

$$\mathrm{\Delta }J_{\mathrm{Kepler}}^2=2|𝒜|\mathrm{sin}\phi _0^{2\pi }r^3\mathrm{sin}(\theta )𝑑\theta .$$
(82)
The perturbation in the orbit is assumed linear in $`𝒜`$; therefore it is sufficient to substitute in equation (82) the Keplerian orbit:

$$r(\theta )=\frac{GM_{}}{2E}\frac{(1e^2)}{1+e\mathrm{cos}(\theta \psi )},$$
(83)
where $`\psi `$ is the inclination of the orbit with respect to the projection of $`𝒜`$ on the orbital plane. Integration yields:

$$\mathrm{\Delta }J_{\mathrm{Kepler}}^2=6\pi |𝒜|\left(\frac{GM_{}}{2E}\right)^3e\sqrt{1e^2}\mathrm{sin}\phi \mathrm{sin}\psi .$$
(84)
Note that in Kepler’s potential, $`e^2=1R`$; $`J_c(E)=GM_{}/\sqrt{2E}`$; and $`P(E)=2\pi GM_{}(2E)^{3/2}`$. Therefore, knowing that $`R_{\mathrm{lc}}1`$ at the loss cone boundary, we have:

$$e^2(1e^2)_{\mathrm{lc}}R_{\mathrm{lc}}r_{\mathrm{bin}}\left(\frac{GM_{}}{2E}\right)^1.$$
(85)
Substituting equations (84) and (85) in equation (81) gives the time scale on which Brownian motion refills the loss cone:

$$t_{\mathrm{brown},\mathrm{Kepler}}\frac{(GM_{})^{3/2}r_{\mathrm{bin}}}{6\pi 𝒜^2}\left(\frac{GM_{}}{2E}\right)^{7/2}.$$
(86)
The time scale was derived ignoring the contribution to the mean potential from the stars in the galaxy. If the galactic potential had been taken into account and the total potential had not been Keplerian, the functional dependence on $`E`$ in the time scale (86) would have had a different form. The scaling in the other parameters ($`𝒜`$, $`r_{\mathrm{bin}}`$, $`M_{}`$), however, will be the same in the generic case. Substituting the expression for the acceleration from equation (77) yields:

$`t_{\mathrm{brown}}`$ $``$ $`400\text{ Myr}\times {\displaystyle \frac{r_{\mathrm{bin}}}{a_{\mathrm{hard}}}}\left({\displaystyle \frac{m_{}}{M_{}}}\right)^1`$ (87)
$`\times \left({\displaystyle \frac{M_{}}{10^6M_{}}}\right)^{23/\alpha }e^{14/\alpha }K(x),`$
where $`a_{\mathrm{hard}}GM_{}/8\sigma ^2`$ is the hard binary separation, $`\alpha d\mathrm{ln}M_{}/d\mathrm{ln}\sigma `$ is the slope of the relation between the black hole mass and the velocity dispersion, and $`K(x)`$ is a dimensionless factor encapsulating the dependence on the energy $`E2\sigma ^2x`$. Using $`\alpha 4.5`$ (Ferrarese & Merritt, 2000) and $`r_{\mathrm{bin}}/a_{\mathrm{hard}}0.2`$ in a galaxy like M32, we obtain $`t_{\mathrm{brown}}5`$ Gyr at energies corresponding to the radius of influence of the black hole where $`K(x)1`$. The total mass inside the loss cone at these energies is small; thus the Brownian motion-driven diffusion of stars into the loss cone does not appear to be efficient in galaxies.

The effect of the Brownian motion in galaxies could be enhanced if clumpy inhomogeneities are present in the galactic potential. Backer & Sramek (1999) point out that the mass of inhomogeneities at the center of the Milky Way varies with the projected distance from position of Sgr A, assuming the distance of 8.5 kpc to the Galactic center, as:

$$m(r)10^3M_{}\left(\frac{r}{\text{pc}}\right)^{1.5}.$$
(88)
The accelerations implied by these inhomogeneities, $`𝒜Gm(r)/r^210^6\text{ cm s}^1`$ at $`r=1`$ pc, are comparable to the accelerations due to the stars that are given by equation (77) and amount to $`𝒜_{}=0.53\times 10^6\text{ cm s}^1`$ using $`m_{}=M_{}`$, $`M_{}=2.6\times 10^6M_{}`$, and $`\sigma =75110\text{ km s}^1`$. The presence of compact inhomogeneities that do not suffer tidal truncation in the galactic potential, such as intermediate-mass black holes, could in principle enhance the amplitude of the Brownian motion beyond the predictions of this section.

Scaling the results of this section to the simulations of Milosavljević & Merritt (2001), we find that the time scale on which the Brownian motion fed stars into the loss cone amounted to $`t_{\mathrm{brown}}1\text{ Myr}\times K(x)`$, which is of the order of the dynamical time. While we argued above (§ 6) that gravitational scattering into the loss cone was sufficient to keep it full in the $`N`$-body simulations, Brownian motion probably also contributed substantially to the maintenance of a full loss cone in these simulations. This underscores a critical difference between the published $`N`$-body simulations and real galaxies; the $`N`$-body simulations are dominated by small-$`N`$ effects and real galaxies are not.

## 8. Coalescence

When the MBHs approach each other sufficiently closely, emission of gravity waves becomes the dominant mechanism extracting the binding energy. The binary coalesces within a time $`t_{\mathrm{gr}}`$ when its semimajor axis $`a_{\mathrm{gr}}`$ and the eccentricity $`e_{\mathrm{gr}}`$ satisfy the relation derived by Peters (1964):

$$a_{\mathrm{gr}}^4=\frac{64}{5}\frac{G^3M_1M_2(M_1+M_2)F(e_{\mathrm{gr}})}{c^5}t_{\mathrm{gr}},$$
(89)
where $`F(e)`$ has the form:

$$F(e)=\left(1e^2\right)^{7/2}\left(1+\frac{73}{24}e^2+\frac{37}{96}e^4\right).$$
(90)
The ratio between the hard binary separation $`a_{\mathrm{hard}}=G\mu /4\sigma ^2`$ at which MBH binaries first form, and the coalescence separation $`a_{\mathrm{gr}}`$ for circular binaries reads:

$`{\displaystyle \frac{a_{\mathrm{hard}}}{a_{\mathrm{gr}}}}`$ $``$ $`34\times e^{9.2/\alpha }\left({\displaystyle \frac{M_{}}{10^6M_{}}}\right)^{1/42/\alpha }`$ (91)
$`\times p^{3/4}(1+p)^{3/2}\left({\displaystyle \frac{t_{\mathrm{gr}}}{10^9\text{ yr}}}\right)^{1/4}`$
where, again, $`\alpha `$ is the slope of the $`M_{}`$$`\sigma `$ relation, and $`pM_2/M_11`$ is the ratio of the MBH masses. The factor of $`10100`$ that the MBH binary needs to decay on its path from initial formation to coalescence is sometimes called the “final parsec problem.” The required decay factor decreases with decreasing mass ratio; however when $`p10^3`$, the smaller MBH’s host galaxy is likely to be severely tidally disrupted during the inspiral and thus deposit the black hole on an orbit larger than $`a_{\mathrm{hard}}`$. The time scale on which the smaller MBH spirals inward due to dynamical friction could then be too long to permit the formation of a hard binary (Merritt, 2000; Yu, 2002) and the MBHs would not coalesce even if the mechanisms discussed in this paper imply that they should.

We proceed to a qualitative discussion of the scenarios for the evolution of MBH binaries with masses lying within several ranges of values. Bright, early type galaxies ($`M_V<21`$) typically have shallow nuclear density profiles, $`|d\mathrm{log}\nu /d\mathrm{log}r|1`$, within a break radius $`r_b`$ (e.g., Gebhardt et al. 1996), although some exceptions to this rule have been found (e.g., Laine et al. 2003). The break radius is usually located well outside the dynamical radius of influence of the MBH. Therefore MBH binaries with initial black hole masses $`M_{}10^8M_{}`$ are from the outset situated in low-density stellar cusps and the mass initially inside the loss cone is small compared with $`M_{}`$. In addition, the time scales for diffusion into the loss cone are very long ($`t_{\mathrm{diff}}10^{11}`$ yr). Therefore, binaries that form in the mergers of galaxies with black hole masses above $`10^8M_{}`$ have difficulty decaying much beyond the orbital separation of a hard binary. In this mass range, $`a_{\mathrm{hard}}/a_{\mathrm{gr}}10`$, implying that these binaries would be long-lived. Exceptions might be merger events with extreme mass ratios, $`p10^2`$, for which the required decay factors are smaller. If collisional loss cone refilling is very inefficient, the dominant loss cone replenishment mechanism might be via chaotic or centrophilic orbits if the potential is appreciably non-axisymmetric (Merritt & Poon, 2003). Even in axisymmetric, non-spherical potentials, larger stellar masses are available inside the loss cone (Yu, 2002) than suggested by the numerical estimates of the preceding sections (generalization of the results §4–§7 to axisymmetric potentials would be straightforward but tedious). Furthermore, the longer a binary persists in an uncoalesced state, the more likely that a third MBH will fall in during a subsequent merger, facilitating decay by three-body interactions (Makino & Ebisuzaki, 1994; Valtonen, 1996b). Thus we would still expect coalescence to sometimes occur for high-mass binaries.

At the other end of the mass range, $`M_{}\text{ few}\times 10^6M_{}`$, the black holes are expected to be embedded in approximately isothermal stellar cusps, $`|d\mathrm{log}\nu /d\mathrm{log}r|2`$. The formidable factor that the binaries need to cover from formation to coalescence, $`a_{\mathrm{hard}}/a_{\mathrm{gr}}100`$ is slightly larger than for more massive binaries due to the $`M_{}`$-dependence in equation (91). But a combination of processes contribute to effective orbital decay in these galaxies. Foremost, $`N`$-body simulations (Milosavljević & Merritt, 2001) have shown that the binaries shrink by a factor of $`510`$ beyond the hard binary separation immediately following the formation of a binary. This takes place within a dynamical time and is a consequence of the abundance of cold (large binding energy) stars in the nascent loss cone to which the kinetic energy of the MBHs can be effectively transferred in slingshot ejections (see §5.4). The ejection of most bound stars leads to a reduction of the density profile close to the sphere of influence of the MBH. However there are good reasons to believe that galaxies with $`M_{}10^6M_{}`$ are in position to repair cusp damage on a time scale comparable to the binary hardening time scale, since both time scales are proportional to the relaxation time. Indeed, the phase space transport due to two-body relaxation will tend to rebuild a cusp profile $`\rho r^{7/4}`$ around the MBH (Bahcall & Wolf, 1976). Moreover, the nearest and the best resolved MBH nucleus—Sgr A cluster at the center of the Milky Way—exhibits evidence for young, massive stars of spectral type O-B deep inside the radius of influence of the MBH (e.g., Figer et al. 2000; Ghez et al. 2003). This suggests that even in galaxy nuclei with $`M_{\mathrm{SgrA}^{}}4\times 10^6M_{}`$ black holes, star formation can repopulate the central cusp in as short a time as $`10^7`$ yr.

After the loss cone is empty of its original content, two-body relaxation lets stars diffuse into the loss region. Note that the binary is still a factor $`1020`$ wider than the separation conducive to gravitational coalescence in a Hubble time. From equation (11), the equilibrium diffusion supplies an influx of stars at the rate:

$$\frac{dM_{\mathrm{lost}}}{dt}\frac{3M_{}}{10^{10}\text{ yr}}\left(\frac{M_{}}{10^6M_{}}\right)^1\left(\frac{m_{}}{M_{}}\right),$$
(92)
and the non-equilibrium diffusion could in principle supply stars at a rate larger by a factor of $`_{\mathrm{noneq}}/_{\mathrm{eq}}210`$. Combining equations (92) and (61), setting $`m_{}=M_{}`$, and using $`\mathrm{\Delta }\mathrm{\Phi }/2\sigma ^210`$ which is reasonable for isothermal galaxies with $`M_{}10^3M_{\mathrm{gal}}`$, we find:

$`{\displaystyle \frac{d\mathrm{ln}a^1}{dt}}`$ $``$ $`3\times 10^{10}\text{ yr}^1\left({\displaystyle \frac{a}{0.1a_{\mathrm{hard}}}}\right)`$ (93)
$`\times \left({\displaystyle \frac{M_{}}{10^6M_{}}}\right)^1\left({\displaystyle \frac{_{\mathrm{noneq}}}{_{\mathrm{eq}}}}\right),`$
where we have explictly indicated that diffusion becomes the dominant driver of binary hardening only after the immediate post-merger shrinkage to $`a0.1a_{\mathrm{hard}}`$ has been completed. The factor describing the enhancement of the loss cone flux out of equilibrium, $`_{\mathrm{noneq}}/_{\mathrm{eq}}`$, is greatly uncertain and large ($`10`$) when the phase-space is subject to episodic perturbations excited by the substructure inside the galaxy (orbiting giant molecular clouds, super-star clusters, etc.) which can restore the large phase-space gradients at the loss cone boundary (see § 4). We integrate equation (93) to obtain:

$$\frac{0.1a_{\mathrm{hard}}}{a}\frac{t}{3\times 10^9\text{ yr}}\left(\frac{M_{}}{10^6M_{}}\right)^1\left(\frac{_{\mathrm{noneq}}}{_{\mathrm{eq}}}\right).$$
(94)
Although this result suffers from a variety of uncertainties, it suggests that for $`M_{}\text{ few}\times 10^6M_{}`$, the binary separation can decay by a factor of at least a few out of the factor of $`\times 1020`$ that remain toward gravitational coalescence. Hence we expect that stellar-dynamical processes might often serve to bring MBH binaries in this mass range to complete coalescence.

Finally, MBHs in the intermediate range of $`10^{6.5}M_{}M_{}10^8M_{}`$ live in nuclei with a range of cusp slopes, $`0.5<|d\mathrm{log}\nu /d\mathrm{log}r|<2`$ (Gebhardt et al., 1996). The mechanisms discussed here would facilitate coalescence in galaxies with steeper cusps, while in galaxies with shallower cusps, binaries might be shrunk to separations within a factor of a few of $`a_{\mathrm{gr}}`$, where other mechanisms (triaxiality, gas dissipation etc.) could plausibly bridge the remaining gap.

## 9. Summary and Discussion

We conclude with a tabulated summary of physical regimes that have been discussed so far (Table 1).

Clearly the long-term evolution of a massive black hole binary can be very different in different environments. We distinguish three characteristic regimes.

1. Collisional. The relaxation time is shorter than the lifetime of the system and the phase-space gradients at the edge of the loss cone are given by steady-state solutions to the Fokker-Planck equation (§ 2). The densest galactic nuclei may be in this regime. Resupply of the loss cone takes place on the time scale associated with scattering of stars onto eccentric orbits. Most of this scattering occurs in the “diffusive” (as opposed to “pinhole”) regime and the decay time of a MBH binary scales as $`|a/\dot{a}|m_{}^1N`$. In the densest galactic nuclei, collisional loss cone refilling may just be able to drive a MBH binary to coalescence in a Hubble time (e.g. Yu 2002). $`N`$-body studies with any feasible $`N`$ are guaranteed to be in the collisional regime, and may also be in the pinhole regime corresponding to an always-full loss cone, leading to artifically high decay rates, $`a^1t`$ (§ 6). In $`N`$-body studies with $`N10^6`$, Brownian motion of the binary further enhances the decay (§ 7).

2. Collisionless. The relaxation time is longer than the system lifetime and gravitational encounters between stars can be ignored. Examples are the low-density nuclei of bright elliptical galaxies; dark matter density cusps are also in this regime if the dark matter particles are much less massive than stars. The MBH binary quickly interacts with stars whose pericenters lie within its sphere of influence; in a low-density nucleus, the associated mass is less than that of the MBH binary and the decay tends to stall at a separation too large for gravity wave emission to be effective. However evolution can continue due to reejection of stars that lie within the binary’s loss cone but have not yet escaped from the system. In the spherical geometry, reejection implies $`|a/\dot{a}|(1+t/t_0)/a`$, leading to a logarithmic dependence of binary hardness on time (§ 5). Reejection in galactic nuclei may contribute a factor of $``$ a few to the change in $`a`$ in MBH binaries over a Hubble time.

3. Intermediate. The relaxation time is of order the age of the system. While gravitational encounters contribute to the re-population of the loss cone, not enough time has elapsed for the phase space distribution to have reached a collisional steady state. We argued that most galactic nuclei are likely to be in this regime. The flux of stars into the loss cone can be substantially higher than predicted by the steady-state theory, due to strong gradients in the phase space density near the loss cone boundary produced when the MBH initially formed (§ 4). This transitory enhancement would be most important in a nucleus that continues to experience mergers or infall, in such a way that the loss cone repeatedly returns to an unrelaxed state with its associated steep gradients.

The evolution of a real MBH binary may reflect a combination of the different regimes summarized above, as well as other mechanisms that we have discussed only briefly or not at all. We note a close parallel between the “final parsec problem” and the problem of quasar fueling: both requre that of order $`10^8M_{}`$ be supplied to the inner parsec of a galaxy in a time shorter than the age of the universe. Nature clearly accomplishes this in the case of quasars, probably through gas flows driven by torques from stellar bars (Shlosman, Begelman & Frank, 1990). The same inflow of gas could contribute to the decay of a MBH binary in a number of ways: by leading to the renewed formation of stars which subsequently interact with the binary; by inducing torques which extract angular momentum from the binary (Armitage & Natarajan, 2002); through accretion, increasing the masses of one or both of the MBHs and reducing their separation; etc.

In massive elliptical galaxies or spheroids, the fraction of mass in the form of gas during the last major merger event is likely to have been small (Kauffmann & Haehnelt, 2000). In addition to the mechanisms discussed in this paper, decay of a MBH binary could be enhanced by the presence of a significant population of “centrophilic” orbits, orbits that pass near the center of the potential once per crossing time. Such orbits may constitute a large fraction of the mass of a triaxial nucleus (Poon & Merritt, 2002). Even in galaxies where none of these mechanisms is effective and the decay stalls, infall of a third MBH following a merger or accretion event could accelerate the decay via the Kozai mechanism (Blaes, Lee & Socrates, 2002), or the strong three-body gravitational interactions (Makino & Ebisuzaki, 1994; Valtonen, 1996b), which induce large changes in the binary’s eccentricity and enhanced rates of gravity wave emission at pericenter.

We have not attempted in this paper to make quantitative estimates of the decay rates or stalling radii of MBH binaries in real galaxies. Such a program would be problematic for a number of reasons, most importantly the uncertain history of the loss cone and the unknown elapsed time and character of the most recent merger (§ 4). However we argued that a recent study (Yu, 2002) could have overestimated stalling radii, due both to the neglect of the various new mechanisms discussed here, and also by ignoring the influence that MBH formation must have had on the form of the nuclear density profile: this profile must have been steeper before the MBH binary formed, leading to substantially enhanced loss cone fluxes and more rapid decay early in the life of the binary. Nevertheless, none of our results would lead us to rule out the existence of uncoalesced MBH binaries in at least some galaxies, particularly galaxies with large, low-density nuclei and a little gas (§ 8).

A striking and disappointing conclusion of this study is the difficulty of using numerical $`N`$-body simulations to follow the long-term evolution of MBH binaries (§ 6). We argued that a variety of collisional effects associated with small $`N`$ would lead to a more rapid decay of the binary than in real galaxies. While there is still a useful role for $`N`$-body simulations in following the initial formation and early evolution of a MBH binary, when much of the decay and associated cusp destruction takes place, the rapid decay of $`a(t)`$ seen in most published simulations is due to loss-cone repopulation occuring at rates much higher than expected in real galaxies. We believe that the future role of $`N`$-body simulations in this field will be limited to predicting the form of the phase-space gradients produced during the formation of a MBH binary, which can then be adapted as initial conditions for a collisionless or Fokker-Planck treatment. A carefully designed set of $`N`$-body experiments might also be used to understand the $`N`$-dependence of collisional loss cone effects which could then be extrapolated to the larger-$`N`$ regime characteristic of real galaxies.

We are grateful to A. Kosowsky, C. Pryor, and R. Spurzem for their detailed comments on an early version of the manuscript. We acknowledge helpful discussions with J. Binney, S. Phinney, F. Rasio, and M. Rees. This work was supported by NSF grants AST 96-17088 and AST 00-71099 and by NASA grants NAG5-6037 and NAG5-9046 and STScI grant HST-AR-08759. M. M. thanks the Sherman Fairchild Foundation for support.