LONG-TERM EVOLUTION OF MASSIVE BLACK HOLE BINARIES. IV. MERGERS OF GALAXIES WITH COLLISIONALLY RELAXED NUCLEI

Figures 3 and 4 illustrate the large-scale evolution of the bulge models described in Table 1 during the galaxy merger phase. The two bulge models start out either on circular (Model A) or eccentric (Model B) orbits about their common center of mass, with initial separations roughly one-half the radius of the larger bulge. Due to its more eccentric initial orbit, Model B evolves more quickly. The trajectories of the MBHs reflect the initial orbits of the parent systems, as can be seen in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Evolution of a binary MBH can be divided roughly into four phases. (1) Formation of a bound pair. Merger of two galaxies with central MBHs brings the two MBHs together, in a time comparable with the galaxy merger time, i.e., a few galaxy crossing times. (2) Formation of a hard binary. The separation between the two MBHs decreases very rapidly, due at first to dynamical friction against the stars, and then to the gravitational slingshot: near stars are ejected after gaining energy from the massive binary. A core is formed at this stage, with size roughly equal to the initial separation of the bound pair. (3) Binary hardening. If orbits of stars ejected by the gravitational slingshot are replenished, the massive binary will continue to shrink. Hypervelocity stars can be produced in this process. The stellar core will continue to grow. (4) Coalescence. If replenishment of stellar orbits continues, the binary separation decreases to a value such that emission of gravitational waves dominates its evolution, and the two MBHs coalesce.

Dynamical friction drives the evolution down to a separation a_f at which the stellar mass enclosed in the binary is of order twice the mass of the smaller MBH:

$$\begin{equation} M ({<} a_f) \approx 2\,M_2\,. \end{equation} \tag{ 8 }$$

This separation is smaller by a factor ∼M₂/M₁ than the radius of influence of the larger MBH.

At separations smaller than a_f, the binary hardens due to gravitational slingshot interactions with passing stars. A binary is defined as "hard" when it reaches a separation

$$\begin{equation} a_h \approx \frac{{\it GM}_2}{4\,\sigma ^2}, \end{equation} \tag{ 9 }$$

called the hard-binary separation; a binary is "hard" when its binding energy per unit mass, |E|/(M₁ + M₂), exceeds σ².

The binary enters the gravitational wave (GW) regime when the timescale for coalescence due to emission of gravity waves,

$$\begin{eqnarray} T_{\rm GW} & = & \frac{5}{256 F(e)} \frac{c^5}{G^3} \frac{a^4}{\mu \left(M_1 + M_2 \right)^2} \nonumber \\ & \approx & \frac{5.8 \times 10^{11} {\, \rm yr}}{F(e)} \left(\frac{a}{0.1{\, \rm pc}}\right)^4 \left(\frac{10^7\, M_\odot }{\mu }\right) \left(\frac{10^8\, M_\odot }{M_1+M_2}\right)^2,\nonumber \\ \end{eqnarray} \tag{ 10 }$$

becomes shorter than the time for hardening due to stellar interactions. Here

$$\begin{eqnarray*} &&F(e) = (1-e^2)^{-7/2} \left(1 + \frac{73}{24}e^2 + \frac{37}{96} e^4\right) \end{eqnarray*}$$

and μ ≡ M₁ M₂/(M₁ + M₂) is the reduced mass.

The evolution of the separation between the MBHs is shown in Figure 6. The first part of this evolution, which is driven by dynamical friction, ends when the separation reaches ∼a_f, Equation (8); a_f ≈ 0.9 in model units. The time to reach this separation was t_f ≈ 440 and ∼160, respectively, for models A and B. At about the same time, gravitational slingshot encounters with the stars begin to dominate the binaries' evolution and the effect on the galaxy structure starts to become apparent. Figure 7 clearly shows a decrease in the central densities of both the MS and the stellar BHs at a time t ≈ t_f. Similar expansions are observed in the WD and NS distributions. However, the stellar BHs are most affected, due to their higher, initial central concentration.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The orbital semi-major axis and eccentricity of the massive binary are shown as functions of time in Figure 8. Even in the circular-orbit merger (Model A), the massive binary forms with a slightly nonzero eccentricity. Subsequent evolution of a and e is driven by interactions with stars on intersecting orbits.

Zoom In Zoom Out Reset image size

We computed the time-dependent binary hardening rate (Quinlan 1996)

$$\begin{equation} s \equiv \frac{d}{dt} \left(\frac{1}{a}\right) \end{equation} \tag{ 11 }$$

by fitting a straight line to a⁻¹(t) in small time intervals between the time of binary formation and the end of the integration. The results are shown in the bottom panels of Figure 9 for both models. The top panels show the evolution of 1/a. The hardening rate appears roughly constant in time, which suggests that the binaries enter the hard phase rather quickly, and then continue hardening at an approximately constant rate. This behavior is a defining property of hard binaries, if the stellar background is unchanging, i.e., unaffected by the binary.

Zoom In Zoom Out Reset image size

Mikkola & Valtonen (1992), Quinlan (1996), and Sesana et al. (2006) performed three-body scattering experiments of circular binaries of varying mass ratio and hardness to derive estimates of the hardening rate. They provided fitting formulae for the dimensionless parameter H which is related to the hardening rate via

$$\begin{equation} H(a) = \frac{\sigma }{G\,\rho } \frac{d}{dt}\left(\frac{1}{a}\right)\,, \end{equation} \tag{ 12 }$$

where ρ is the stellar mass density and σ is the one-dimensional velocity dispersion, both assumed constant in space and time and unaffected by the presence of the binary. The H parameter obtained from these scattering experiments is a function of binary mass ratio and hardness; the latter is defined in terms of V_bin/σ, the ratio of binary circular speed to field-star velocity dispersion.

We wish to compare the N-body results for s with the predictions from the scattering experiments. Good agreement would imply a "full loss cone," i.e., that the rate of interaction of the binary with stars is unaffected by slingshot ejections. In spherical galaxy models, binary hardening rates fall much below the predictions from the scattering experiments, since orbital repopulation is driven by two-body scattering, which is a slow process for large N (Makino & Funato 2004; Berczik et al. 2005; Merritt et al. 2007b).

We defined the theoretically expected hardening rate to be

$$\begin{equation} s_{\rm D}(a) \equiv H(a)\left(\frac{\rho }{\sigma }\right)_D \end{equation} \tag{ 13 }$$

where H is taken from the published scattering experiments, and ρ and σ are measured in our N-body models, at a distance D from the MBH. (G = 1 in N-body units.) For H(a) we adopted the fitting formula of Sesana et al. (2006):

$$\begin{equation} H = A (1+a/a_0)^{\gamma }, \quad a_0 = 1.05 \frac{{\it GM}_2}{\sigma ^2}, \end{equation} \tag{ 14 }$$

with parameters for a 3:1 mass ratio: A = 15.82, γ = −0.95. Because ρ and σ vary with position (and time) in the N-body models, our computed values of s will depend on D. We found this dependence to be weak, as long as D is not too different (i.e., not more than a factor of two) from r_m: both ρ and σ are weakly dependent on radius for r ≈ r_m. In Figure 9 we plot values of s_D for D = 0.5, slightly larger than estimated influence radii in both models. We find that the N-body hardening rates, s, are quite consistent with the analytic predictions s_D. The contributions to the binary hardening from the different mass groups (assumed to scale in proportion to their mass densities) are shown in the bottom panels of Figure 9. The MS stars appear to be responsible for most of the binary hardening, followed by the WDs and the stellar BHs.

These results suggest that the massive binaries in our simulations are roughly in the "full-loss-cone" regime: they harden at a rate that is consistent with the expected hardening rate of a binary in an undepleted field of stars. This is in agreement with the results of Khan et al. (2011) and Preto et al. (2011), who also found that the stalling that occurs in spherical models is absent in simulations that start from a stage preceding merger of the two galaxies. Apparently, the nonspherical shapes of the merger remnants result in a large population of stars on "centrophilic" orbits: saucer orbits in the axisymmetric geometry (Sridhar & Touma 1999), pyramid orbits in the triaxial geometry (Merritt & Vasiliev 2010), etc. We discuss the shapes of our models in more detail in Section 6.

Figure 8 shows that the eccentricity of the massive binary remains roughly constant in both models. In a nonrotating galaxy, binary eccentricity is expected to increase gradually with time:

$$\begin{equation} \frac{de}{dt} =K \frac{d\ln (1/a)}{dt}, \end{equation} \tag{ 15 }$$

where K > 0 is a second dimensionless rate, also derivable from scattering experiments; for instance, Sesana et al. (2006) give

$$\begin{equation} K = A (1+a/a_0)^{\gamma } + B, \end{equation} \tag{ 16 }$$

and their Table 3 gives values of A and B for a 1:3 binary mass ratio. Stars that encounter the binary in a prograde (corotating) sense tend to circularize it; Sesana et al. (2011) found that when the fraction of corotating stars exceeded ∼0.7 (as opposed to 0.5 for a nonrotating galaxy), binaries tend to circularize. We found that our models have roughly this fraction (∼0.7) of corotating stars, consistent with the mild eccentricity growth that we see.

In our simulations, binary hardening rates s are essentially constant with time. If we assume that this remains true beyond the end of our simulations, we can use Equations (11) and (15) to extrapolate the binary elements (a, e) to arbitrarily later times. At some point, gravitational wave emission will dominate the evolution; calculating when this happens requires assigning physical units to the models. We considered two representative scalings, based on assumed MBH masses of 4.0 × 10⁶ M_☉ and 10⁸ M_☉, respectively. Scaling factors for length were set at 10 pc and 40 pc, respectively, yielding physical values of the influence radii of ∼3 pc and ∼12 pc. As discussed in Section 3, the unit of time is determined in this (full-loss-cone) regime by orbital periods, not relaxation times; hence the scaling factor for time, [T], is related to the scaling factors for mass and length by $[T] = \sqrt{[L]^3/(G[M])}$ and the binary hardening rate in physical units is [L]⁻¹[T]⁻¹ times the N-body hardening rate.

Given a value for s, and scale factors for mass and length, the evolution of the binary semi-major axis at late times is determined by

$$\begin{equation} \frac{da}{dt} = \frac{da}{dt}\bigg |_{\rm SI} + \frac{da}{dt}\bigg |_{\rm GW} = -s\,a^2(t) + \frac{da}{dt}\bigg |_{\rm GW}, \end{equation} \tag{ 17 }$$

where the first term on the right-hand side represents hardening due to interactions with stars, and the second term represents energy lost to gravitational waves. The rate of the latter process depends on e as well as a. In extrapolating e beyond the end of the N-body integrations, we ignored the effect of the galaxies' rotation on the binaries' eccentricity growth and simply applied Equations (15) and (16). The rate of change of the orbital elements due to GW emission is (Peters 1964)

$$\begin{eqnarray} &&\frac{da}{dt}\bigg |_{\rm GW} = -\frac{64}{5} \beta \frac{F(e)}{a^3} \end{eqnarray} \tag{ 18a }$$

$$\begin{eqnarray} &&\frac{de}{dt}\bigg |_{\rm GW} = -\frac{304}{15}\beta \frac{e G(e)}{a^4}, \end{eqnarray} \tag{ 18b }$$

where F(e) is given in Equation (11),

$$\begin{equation} G(e) = (1-e^2)^{-5/2} \left(1 + \frac{121}{304}e^2 \right), \end{equation} \tag{ 19 }$$

and

$$\begin{equation*} \beta = \frac{G^3}{c^5} M_1 M_2 \left(M_1+M_2\right).\nonumber \end{equation*}$$

We numerically solved the coupled equations for the evolution of the semi-major axis and eccentricity, starting from values at the end of the N-body phase or somewhat earlier, for each of the two sets of scaling factors. The results are shown in Figure 10 and Table 3. In this table, the time T_NB represents the time spent in the N-body hardening phase while T_HD represents the time from the end of the N-body integration until coalescence; the latter time includes both a stellar-interaction-driven and a GW-dominated regime. The total evolution time from the beginning of the galaxy merger to full coalescence is given by the sum of the two times. For a Milky Way type galaxy undergoing a 3:1 merger, this time is of the order of 10⁸ yr for an initial circular orbit and 6 × 10⁷ yr for an initially moderately eccentric orbit. The table also lists the values of the semi-major axis a_f and eccentricity e_f at the transition between the simulated N-body phase and the semi-analytical phase, as well as the value of the eccentricity when the separation reaches 10r_g = 10 GM₁/c². At smaller separations, equations like (18) are not valid.

Zoom In Zoom Out Reset image size

The pre-merger galaxies had central density profiles that were well approximated as power laws with respect to radius. Following the merger, the density on scales r ≲ r_h is strongly modified (lowered) by the action of the massive binary. The process of "core scouring" is illustrated in Figure 11, which compares the spatial density profiles of the different species in the larger galaxy at four different times. For this plot, only particles that were originally associated with the larger galaxy were used. By the time the binary has become hard, a mass deficit (Milosavljević et al. 2002) is clearly visible in all components on scales significantly larger than r_m. Model A shows substantially more evolution of the central density at early times; at late times the cores in the two models are more similar.

Zoom In Zoom Out Reset image size

We estimated core sizes in two ways. First, we fit core-Sérsic profiles to our models at the end of the merger phase and adopted the resulting values of the break radius R_b as estimates of the core radius. The projected density profiles of galaxies are globally well fit by the Sérsic (1968) law:

$$\begin{equation} I(R) = I(0) {\rm exp}\lbrace -b (R/R_e)^{1/n} \rbrace, \end{equation} \tag{ 20 }$$

where I(0) is the central intensity, R_e is the effective half-light radius, and n is a parameter controlling the curvature in a log–log plot. The term b is not a parameter but a function of n (see, e.g., Terzić & Graham 2005). Deviations from this law appear close to the center, where bright galaxies typically show central light deficits with respect to the inward extrapolation of the Sérsic law while faint galaxies show central light excesses (e.g., Graham & Guzmán 2003; Ferrarese et al. 2006b; Côté et al. 2007). Adding an inner power law with slope γ (Graham et al. 2003) yields

$$\begin{equation} I(R) = I^{\prime }\, [1+(R_b/R)^{\alpha }]^{\gamma /\alpha }{\rm exp}\left\lbrace -b\, \left(R^\alpha + R_b^\alpha \right)\big/R_e^\alpha \right\rbrace ^{1/\left(n \alpha \right)} \end{equation} \tag{ 21 }$$

with

$$\begin{equation} I^{\prime } = I_b\,2^{-\gamma /\alpha } {\rm exp}[ b\, 2^{1/\left(\alpha n\right)} \left(R_b/R_e\right)^{1/n} ] \,, \end{equation} \tag{ 22 }$$

the so-called core-Sérsic law. Here, α is an additional parameter regulating the transition from the inner power law to the outer Sérsic law and R_b, the break radius, marks the distance where the profile changes from one regime to the other.

We also considered a second definition of the core radius (King 1962), as the projected radius r_c at which the surface density falls to one-half its central value. Here, "central" was taken to be the value at a projected radius of 0.1r_m.

Both estimates, computed at the time when a ∼ 0.1 a_h, are listed in Table 4, separately for the MS stars and the stellar-mass BHs. These radii were computed using all the particles from the respective mass groups, without regard to the galaxy with which they were originally associated. We also give r_m, computed using Equation (6); we set M_• = M₁ + M₂ in that equation, and combined together all the stellar species from both galaxies when computing M_⋆.

In the case of the dominant (MS) mass component, Table 4 shows that r_c ≈ 2r_m for both Models A and B. However, R_b is somewhat larger than r_c. The reason is that, due to the flatness of the profile, the radius at which the inner profile transitions to the outer profile (i.e., the break radius) is much larger than the radius at which the density falls to half its central value (i.e., the core radius). This is illustrated in Figure 12 for Model B.

Zoom In Zoom Out Reset image size

In what follows, we will adopt r_c as the more robust measure of the core radius. We emphasize that the size of the core in the MS stars is larger than the influence radius of the massive binary in both of our models, whichever definition of "core radius" or "influence radius" is used. One important consequence, discussed in more detail below, is that regrowth of a Bahcall–Wolf cusp at radii ≲ r_m following coalescence of the two MBHs leaves the core structure essentially unchanged.

Stars can become unbound during a galaxy merger due to both tidal stripping in the early phases of the merger and gravitational slingshot interactions with the binary MBHs when this is hard. Unbound stars at the end of the merger phase are a result of both mechanisms. We determined the number of unbound stars in the two merger simulations by selecting stars with velocity in excess of the local escape speed from the system. At any given time, the escape velocity at a distance r from the binary center of mass is $V_{\rm esc}(r) = \sqrt{-2\Phi (r)}$, where the gravitational potential Φ can be expressed as (Binney & Tremaine 1987, Equations (2)–(22))

$$\begin{equation} \Phi (r) = \Phi (r^{\prime } < r) + \Phi (r^{\prime } > r) \end{equation} \tag{ 23 }$$

with

$$\begin{equation} \Phi (r^{\prime } < r) = -4\pi G \, \frac{1}{r} \int _0^r \rho (r^{\prime }) r^{\prime 2} dr^{\prime } \approx \frac{-{\it GM}({<}r)}{r} \end{equation} \tag{ 24 }$$

having defined M(< r) the mass enclosed within a sphere of radius r, and

$$\begin{equation} \Phi (r^{\prime } > r) = -4\pi G \int _r^{\infty } \! \rho (r^{\prime }) r^{\prime } dr^{\prime } \approx -G \sum _{i=1}^N \frac{m_i}{r_i}, \mbox{ for } r_i > r {.} \end{equation} \tag{ 25 }$$

Table 5 reports the total fraction of unbound stars f_e at the time when a ∼ 0.1a_h; the fraction of unbound stars originally belonging to the first and second galaxy f_e1, f_e2, and the fraction of unbound MS, WD, NS, and BHs f_MS, f_WD, f_NS, andf_BH, normalized to the total number of stars in their mass group. We find that about 2% of all stars are ejected by the end of the merger phase. Model A, which has a more gradual evolution, produces more unbound stars than Model B. Stars from the different mass groups show roughly equal probabilities of being ejected.

We distinguished stars ejected via tidal stripping (TS) versus gravitational slingshot (GS) based on the time they become unbound: if a star becomes unbound before the binary has become hard we consider it ejected by TS whereas if it becomes unbound after the binary has become hard we consider it ejected by GS. Based on this selection, we find that TS is responsible for the ejection of about 90% of the escapers. The smaller galaxy is the most susceptible to TS, as shown by the fact that f_e2 ≫ f_e1. Table 5 also lists the mass in stars unbound by TS in units of the total stellar mass M_gal and by the mass in stars unbound by GS in units of the binary mass. High-velocity escapers, with velocities as high as several times the central escape speed, are mainly produced by GS. These stars are analogs of the hypervelocity stars detected in the halo of the Milky Way (e.g., Brown et al. 2005, 2006; Gualandris et al. 2005; Gualandris & Portegies Zwart 2007; Gvaramadze et al. 2009). Stars unbound by TS tend to have lower velocities, and would be less numerous in a real galaxy, due to the deeper potential well.

We carried out integrations of the isotropic Fokker–Planck equation describing galaxies containing two stellar mass groups and an MBH. We used these integrations to address two questions. (1) What is the nature of the "accelerated cusp growth" that Preto & Amaro-Seoane observed in their Fokker–Planck integrations? (2) Is the evolution that we observe in our multi-component N-body models consistent with the predictions of Fokker–Planck models?

We begin by summarizing the evolution equations.³

Let f_i(E, t) be the phase-space number density of the ith species at time t and (binding) energy E, where E = −v²/2 + ψ and ψ = −Φ with Φ the gravitational potential. Let i = 1 denote MS, of mass m₁, while i = 2 denotes stellar-mass BHs, of mass m₂; we set m₂/m₁ = 10 as in the N-body integrations. The evolution equation for the ith component is

$$\begin{eqnarray} &&4\pi ^2 p(E)\frac{\partial f_i}{\partial t} = -\frac{\partial F_i}{\partial E} \end{eqnarray} \tag{ 26a }$$

$$\begin{eqnarray} &&F_i = \sum _{j=1,2}\left(-{D_{EE}}_{ij}\frac{\partial f_i}{\partial E} -{D_E}_{ij}f_i\right) \end{eqnarray} \tag{ 26b }$$

$$\begin{eqnarray} {D_{EE}}_{ij} &=& 16\pi ^3\Gamma m_j^2\left [q(E)\int _0^Ef_j(E^{\prime },t){\it dE}^{\prime }\right.\nonumber\\ &&\left. +\, \int _E^{\infty }f_j(E^{\prime },t)q(E^{\prime }){\it dE}^{\prime }\right ] \end{eqnarray} \tag{ 26c }$$

$$\begin{eqnarray} &&{D_E}_{ij} = -16\pi ^3\Gamma m_im_j\int _E^\infty f_j(E^{\prime },t)p(E^{\prime }){\it dE}^{\prime } \end{eqnarray} \tag{ 26d }$$

(Merritt 2012). Here, p(E) and q(E) are given by

$$\begin{eqnarray} p(E) &=& 4\int _0^{\psi ^{-1}(E)} r^2 v(E,r) dr, \end{eqnarray} \tag{ 27a }$$

$$\begin{eqnarray} q(E) &=& \frac{4}{3}\int _0^{\psi ^{-1}(E)} r^2 v^3(E,r) dr, \end{eqnarray} \tag{ 27b }$$

with ψ⁻¹(E) the inverse of the potential function, $v(E,r) = \sqrt{2\left[\psi (r)-E\right]}$, and Γ = 4πG²ln Λ. To simplify the calculations, we made the assumption (as in numerous earlier studies, e.g., Preto et al. 2004; Merritt et al. 2007a; Preto & Amaro-Seoane 2010) that the gravitational potential,

$$\begin{equation} \psi (r) = \frac{GM_\bullet }{r} + \psi _\star (r), \end{equation} \tag{ 28 }$$

the sum of the stellar and MBH potentials, was constant with time and given by its value at t = 0. This is a reasonable approximation at all radii: even if the stellar density evolves significantly at r ≲ r_m, the potential at these radii is dominated by the MBH and is nearly unchanging. At r ≳ r_m, there is no significant evolution in the density and the approximation is again valid.

The coupled equations (26) were solved by standard techniques, starting from initial conditions in which the two components had configuration-space densities

$$\begin{equation} \rho _i(r) = \rho _i(0)\left(\frac{r}{r_0}\right)^{-\gamma } \left(1+\frac{r}{r_0}\right)^{\gamma -4} \end{equation} \tag{ 29 }$$

with the same (r₀, γ). This is the same initial mass distribution adopted by Preto & Amaro-Seoane (2010), who set γ = (1/2, 1). Parameters defining the initial conditions were γ; M_•/M_gal, the ratio of MBH mass to the total mass in components 1 and 2; and the number density ratio

$$\begin{equation} {\cal R} \equiv \frac{N_{\rm BH}}{N_{\rm MS}}. \end{equation} \tag{ 30 }$$

As unit of time we adopted the relaxation time defined above, evaluated at the MBH influence radius r_m, with r_m defined as in the N-body models; note that the value of log Λ becomes irrelevant when the time is expressed in this way.

The right panels of Figure 15 show results from a Fokker–Planck integration of the same initial conditions used for the N-body model C2 described above. The agreement is very good.

Figure 16 shows integrations of two models both with γ = 0.5 and M_•/M_gal = 0.05. The models differ in the fraction of heavy objects: ${\cal R}=0$ and ${\cal R}=10^{-3}$. The model with ${\cal R}=10^{-3}$ is the same model plotted by Preto & Amaro-Seoane (2010) in their Figure 3. We have plotted only the density of the lighter (MS) component.

Zoom In Zoom Out Reset image size

The superficial appearance of these plots depends strongly on the radial range plotted. If the range is sufficiently large, extending to radii ≪r_m, an MS cusp with slope dlog ρ/dlog r ≈ −3/2 catches the eye in the model with BHs. This is the feature emphasized by Preto & Amaro-Seoane (2010). Viewed more closely, in the region 0.02 ≲ r/r_m ≲ 1, the plots tell a different story: the model including BHs exhibits a lower density in MS stars at all times. Except at very small radii, addition of the BHs has the expected effect of reducing the density of the MS stars.

Preto & Amaro-Seoane (2010) attributed the "accelerated cusp growth" in the lighter component to "mass segregation," without stating explicitly the connection between the two phenomena. In fact, as we now argue, the mechanism driving the evolution of the lighter component at small radii is not mass segregation; it is scattering of the light component by the heavy component (Merritt et al. 2007a; Alexander & Hopman 2009).

Consider a two-component system; as before, species no. 1 is the light (MS) component and species no. 2 is the heavy (BH) component. The four diffusion coefficients that appear in the Fokker–Planck equation for the light (MS) component scale with m_i and f_i as

$$\begin{eqnarray*} D_{EE_{11}} &\simeq & m_1^2 f_1\simeq m_1\rho _1,\ \ D_{E_{11}} \simeq m_1^2f_1 \simeq m_1\rho _1\nonumber \\ D_{EE_{12}} &\simeq & m_2^2 f_2\simeq m_2\rho _2,\ \ D_{E_{12}} \simeq m_1m_2f_2 \simeq m_1\rho _2\nonumber. \end{eqnarray*}$$

If m₂ρ₂ ≫ m₁ρ₁, $D_{EE_{12}} \gg D_{EE_{11}}$; in other words, self-scattering is negligible compared with scattering off of BHs. The first-order coefficients are smaller than $D_{EE_{12}}$ by factors of (m₁ρ₁)/(m₂ρ₂) and m₁/m₂ and can also be ignored. The evolution equation for the light component becomes in this limit

$$\begin{equation} \frac{\partial f_{\rm MS}}{\partial t} \approx \frac{1}{4\pi ^2p} \frac{\partial }{\partial E} \left(D_{EE}\frac{\partial f_{\rm MS}}{\partial E}\right), \end{equation} \tag{ 31 }$$

with D_EE given by Equation (26c) after setting m_j = m_BH. The steady-state solution is obtained by setting the term in parentheses (the flux) to zero, yielding

$$\begin{equation} \frac{\partial f_{\rm MS}}{\partial E} =0,\ \ \ \ f_{\rm MS}={\rm const,} \end{equation} \tag{ 32 }$$

independent of f_BH. A constant f_MS corresponds, in a point-mass potential, to a density

$$\begin{equation} \rho _{\rm MS} \propto r^{-3/2}, \end{equation} \tag{ 33 }$$

which describes the MS cusp that appears, at early times and at small radii, in the Fokker–Planck models.

The manner in which this steady state is reached will depend on the initial conditions. In the models considered here, the initial MS density is ρ_MS∝r^−1/2, which corresponds to an f(E) that tends to zero near the MBH—a "hole" in phase space at low energies. Scattering of the MS component by the BHs rapidly "fills in" this hole as it drives f toward a constant value. The result is a sharp increase in the MS density at early times at the smallest radii.

The condition that the evolution of the light (MS) component be dominated by scattering off the heavy (BH) component—as opposed to self-interactions, which tend to build a steeper, Bahcall–Wolf cusp—is ρ_BH ≳ (m_MS/m_BH)ρ_MS ≈ 0.1ρ_MS. In Figure 17, the first time at which this condition is satisfied is marked by open circles. This figure shows evolution of the local slope, |dlog ρ_MS/dlog r|, at two radii, 0.05r_m and 0.002r_m, for several different values of ${\cal R}$, as computed via the Fokker–Planck equation. At the larger radius, the BHs remain a small fraction of the total in most of these models, and the evolution of the MS component is not strongly affected. But because the Bahcall–Wolf cusp builds "from the outside in," its initial growth is hardly reflected at much smaller radii. Here, as the lower panel shows, the dominant effect is scattering by BHs, and the effect of the scattering on the MS density profile is strongly dependent on (roughly proportional to) ${\cal R}$.

Zoom In Zoom Out Reset image size

We summarize our findings in this section as follows.

• 1.  
  One- and two-component Fokker–Planck models reproduce well the behavior seen in N-body integrations with the same initial parameters.
• 2.  
  In two-component (MS + BH) models, addition of the heavy (BH) component results in a slightly lower MS density at radii r ≳ 0.05r_m, due to heating of the stars by the BHs.
• 3.  
  The same heating more promptly modifies the MS density profile at small radii, r ≲ 0.05r_m, at least in models that start from a flat core. The rate of growth of this "mini-cusp" is strongly dependent on the BH fraction.

The fact that the effects of scattering of MS stars by BHs is restricted to small radii, r ≲ 0.05r_m, is consistent with the fact that we do not observe this phenomenon in the N-body simulations: for instance, the density profiles of Figure 14 extend down to only ∼0.05r_m. Only a handful of particles were ever present at smaller radii. At the radii resolvable by the N-body simulations, the Fokker–Planck models predict that the BHs should retard the growth of the Bahcall–Wolf cusp in the MS component (Figure 16), not accelerate it, and this is what we see in the N-body simulations.

We now return to a discussion of the post-merger N-body integrations. Figure 13 showed the evolution of the mass enclosed within r_m for each of the four species in each of the N-body integrations. In this section, we look more closely at how the number of BH particles evolves with time on smaller scales.

Figure 18 plots the cumulative radial distribution of BHs in Model A1 at different times. Time zero in this plot corresponds to the moment that the two MBH particles were combined into one. The smallest radii at which there are any BH particles in the N-body models decrease from ∼0.1r_m at early times to ∼0.01r_m at late times. In order to estimate BH numbers at smaller radii, we carried out regression fits of log N_BH to log r. The fits were performed in the radial interval [r₁, r₂] such that N_BH(< r₁) = 5 and N_BH(< r₂) = 25. The resulting slopes, α ≡ |dlog N_BH/dlog r|, are plotted in Figure 19 for all models as a function of time in units of the initial relaxation time. We find slopes 1 ≲ α ≲ 3 for N_BH(r), which imply slopes in the mass density profile in the range 0 ≲ −dlog ρ/dlog r ≲ 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 20 shows the inferred number of BHs at r < 0.01 pc as a function of time. In making this plot, we did a rough scaling of our models to the nucleus of the Milky Way, assuming an influence radius r_m of 3 pc. The factor

$$\begin{equation*} \frac{4\times 10^6\, M_\odot }{10\, M_\odot } \times \frac{m_{\rm BH}}{M_\bullet } = 3.3\times 10^3 \nonumber \end{equation*}$$

was used to convert the number of BH particles in the simulations to the actual number; the first of these factors contains masses in physical units, the second in the units of the N-body code. Since the number of BH particles at these radii is small, we extrapolated inward under two assumptions: a density profile of constant power-law index and a constant density inside r = 0.1r_m. In the former case, we used the slopes derived from the fits to N_BH(r). Since the fitted slopes are typically large, the former assumption results in much larger, inferred numbers of BHs.

Zoom In Zoom Out Reset image size

It is tempting to compare the numbers so obtained with estimates of N_BH obtained from steady-state Fokker–Planck models of the Milky Way nucleus (Hopman & Alexander 2006; Freitag et al. 2006). Here we note one ambiguity associated with such comparisons. In the Milky Way, the two influence radii r_h and r_m are similar. The first is

$$\begin{equation} r_{\rm h}\equiv \frac{G\,M_\bullet }{\sigma ^2}\approx 3.5{\, \rm pc}\left(\frac{M_\bullet }{4\times 10^6\, M_\odot }\right) \left(\frac{\sigma }{70\,{\rm km\ s}^{-1}}\right)^{-2}. \end{equation} \tag{ 34 }$$

The second is somewhat less certain, but dynamical estimates of the mass distribution in the inner few parsecs (e.g., Schödel et al. 2009; Oh et al. 2009) give r_m ≈ 2 pc, consistent within a factor of two with r_h. The near-equality of r_h and r_m in the Milky Way is due to the fact that the density profile is similar to that of the singular isothermal sphere, ρ ∼ r⁻²; the radius of the Milky Way's core is ∼0.5 pc, substantially smaller than both r_h and r_m. In our post-merger N-body models, on the other hand, core radii and r_m are both substantially larger than r_h. This fact precludes a unique scaling of our models to the Milky Way—at least in the Galaxy's current state. (The Milky Way's core may have been larger in the past (Merritt 2010).)

With this caveat in mind, we assume that r_m determines the scaling of our models to the Milky Way. Figure 20 is based on this scaling. In their collisionally relaxed models, Hopman & Alexander (2006) found

$$\begin{equation} N_{\rm BH}(r<0.01{\, \rm pc}) = 150. \end{equation} \tag{ 35 }$$

These authors assumed a mass function with the same four species as in our models, and with the same relative numbers as in our models A2 and B2 (Table 1). In Figure 20, the inferred number of BHs inside 0.01 pc is very uncertain at the relevant times, i.e., t ≈ 10¹⁰ yr, fluctuating between zero and maximum values of ∼1 (Model A2) and ∼10 (Model B2). These upper limits are factors of ∼10² and ∼10¹, respectively, smaller than in the Hopman & Alexander (2006) models. There is an independent way to reach a similar conclusion. After several relaxation times, N_BH(< 0.01 pc) is ∼10 (Model A2) and ∼100 (Model B2). These numbers, presumably representing the mass distribution in a near steady state, are ∼10 times larger than their values at t = 10¹⁰ yr.

It is interesting that the BH distribution takes so long to reach this (nearly) steady state—a time of at least twice the relaxation time as defined by the dominant (MS) population (Figure 14). This may be due in part to the fact that the MS distribution is also continuously evolving. Another reason is the persistence of a core in the dominant component. Chandrasekhar's dynamical friction coefficient, in its most widely used form (Binney & Tremaine 1987), predicts that the frictional force near an MBH drops essentially to zero if ρ(r) increases more slowly than r^−1/2 toward the center. The (isotropic) phase-space density f(E) corresponding to that density profile has no stars moving more slowly than the local circular speed, and Chandrasekhar's formula identifies the frictional force exclusively with field stars moving more slowly than the massive body. This property of the dynamical friction force was automatically incorporated into the Fokker–Planck calculations presented above, since that equation is based on Chandrasekhar's coefficients in their standard forms (Rosenbluth et al. 1957). In reality, some part of the frictional force in the N-body simulations will come from stars moving faster than the test star, and including the contribution from these stars keeps the frictional force from falling identically to zero (Antonini & Merritt 2011). Nevertheless, one expects dynamical friction to act more slowly in these models than expected based on the application of standard formulae that neglect the special properties of f.

In galaxies containing a single stellar population, a ρ∝r^−7/4 Bahcall–Wolf cusp is expected to appear at radii r ⩽ r_BW ≈ 0.2r_m around the MBH. While the growth time of the cusp is dependent on the initial conditions, simulations starting from a shallow cusp inside r_m show that the stellar density will have reached an approximate steady state after roughly one relaxation time at r_m. We presented an example of such evolution in Figure 15.

In nuclei with two mass groups, e.g., solar-mass stars (MS) and 10 M_☉ BHs, evolution toward a steady state near the MBH depends on the relative numbers and masses in the two groups, as well as on the initial conditions. The results of our N-body and Fokker–Planck integrations of two-component models were presented in Section 5.1. Figure 17 showed that addition of the BHs reduces slightly the rate of formation of a Bahcall–Wolf cusp in the MS (light) component, and also affects the final value of the density-profile slope, which varies from ρ_MS ∼ r^−7/4 when ρ_BH/ρ_MS is small to ∼r^−3/2 as ρ_BH/ρ_MS is increased (Bahcall & Wolf 1977). In addition, if the initial MS distribution is very flat near the MBH, scattering by the heavier BHs can dominate the evolution of the MS component at early times for r ≲ 0.05r_m, converting ρ_MS ∼ r^−1/2 to ρ_MS ∼ r^−3/2 at these small radii, even before the Bahcall–Wolf cusp has fully formed at larger radii.

The full galaxy merger simulations presented here allowed us, for the first time, to evaluate the evolution of multi-component nuclei starting from initial conditions that were motivated by a well-defined physical model. Our pre-merger galaxies contained mass-segregated nuclei with four mass components, representing an evolved stellar population. These initial distributions were modified both by the galaxy merger and by the formation of an MBH binary, which created a large core in each of the components (Figures 7 and 11). The core radius of the heaviest (BH) component was somewhat smaller than the cores in the three lighter components, a relic of the earlier mass segregation, and of the incomplete cusp destruction process during the binary MBH phase (Figure 11). Evolution during the merger phase set the "initial conditions" of the nucleus at the time when the two MBHs were combined into one. Unlike the rather ad hoc initial conditions used in many earlier studies of nuclear evolution (e.g., Freitag et al. 2006; Preto & Amaro-Seoane 2010), our post-merger models have cores that should be reasonable representations of the cores in real galaxies that formed via dissipationless mergers, with mass densities that decline toward the center and anisotropic kinematics that reflect the action of the massive binary on the stellar orbits (Figure 23).

By continuing the evolution of these multi-component models for a time greater than T_r(r_m) after coalescence of the two MBHs, we found that the cores characterizing the distribution of the dominant, MS component persisted; in fact, the MS density at radii ∼r_m decreased gradually with time (Figure 13)—a predictable consequence of the initial extent of the cores, several times r_m, implying T_r(r_c) > T_r(r_m), and of continued "heating" by the heavier BHs. Nevertheless, a Bahcall–Wolf cusp gradually reformed in the lighter components at radii r ≲ 0.2r_m ≪ r_c; growth times were found to be ≳ T_r(r_m), somewhat greater than in earlier models with more idealized initial conditions (e.g., Freitag et al. 2006).

Our pre-merger galaxies (models A1, B1) contained larger numbers of remnants than would be expected based on a standard IMF; this was done in order to better resolve the distributions of those components near the MBH. We tested the dependence of the post-merger evolution on the assumed mass function by carrying out a second set of integrations in which we decreased the relative numbers of remnants (NS, WD, BH) by factors of a few to values more consistent with standard IMFs (models A2, B2). Evolution of the dominant, MS component in the latter models differed only modestly from its evolution in the models with larger remnant fractions; the main difference was a lower rate of core expansion reflecting a lower rate of heating by the BHs (Figure 13). The rate of growth of the Bahcall–Wolf cusp in the MS component was essentially unchanged (Figure 14).

Our results on the regeneration of Bahcall–Wolf cusps following dissipationless mergers are consistent with those obtained in simulations of single-component galaxies (Merritt & Szell 2006). We can summarize these results by stating that regrowth of a cusp in the dominant stellar component requires a time comparable with, or somewhat longer than, the relaxation time of that component measured at the MBH influence radius. The new simulations presented here suggest that timescales for cusp regrowth are only weakly dependent on the number of heavy remnants (BHs), at least if the fraction of mass in the heavier population does not exceed a few percent of the total.

The Milky Way nucleus is near enough that a Bahcall–Wolf cusp could be resolved if present, and the dominant stellar population is believed to be old. Since r_h ≈ r_m ≈ 2–3 pc in the Milky Way, the expected outer radius of the Bahcall–Wolf cusp is r_BW ≈ 0.5 pc ≈ 10''. As is well known, number counts of the dominant, old stellar population show no evidence of a rise in density at this radius; instead the number counts are flat, or even falling, from ∼10'' into at least 1'' projected radius (Buchholz et al. 2009; Do et al. 2009; Bartko et al. 2010).

Assuming solar-mass stars, the relaxation time at the influence radius of Sgr A* is 20–30 Gyr (Merritt 2010), while the mean stellar age in the nuclear star cluster is estimated to be ∼5 Gyr (Figer et al. 2004). The time available for formation of a cusp is therefore (5/25)T_r(r_m) ≈ 0.2T_r(r_m). Our simulations (e.g., Figure 14) suggest that a Bahcall–Wolf cusp in the dominant component is unlikely to have formed in so short a time, and this is consistent with the lack of a Bahcall–Wolf cusp at the Galactic center.

The core in the Milky Way has a radius of ∼0.5 pc, somewhat smaller than r_h or r_m, while the cores formed in our merger models are somewhat larger than r_m (Table 4). It has been argued that the Milky Way core is small enough that gravitational encounters would cause it to shrink appreciably in 10 Gyr, as the stellar distribution evolves toward a Bahcall–Wolf cusp (Merritt 2010). The cores in our N-body models are so large that they do not evolve appreciably after the binary MBH has been replaced by a single MBH; the Bahcall–Wolf cusp forms at radii smaller than r_c, leaving the core structure essentially unchanged. Without necessarily advocating a merger model for the origin of the Milky Way core (it is unclear whether our Galaxy even contains a bulge; Martinez-Valpuesta & Gerhard 2011), we note that the sizes of cores formed by binary MBHs scale with the binary mass ratio. Presumably, we could have produced cores more similar in size to the Milky Way's if we had adjusted this ratio.

It has been suggested (Löckmann et al. 2010) that heating by BHs could be responsible for the lack of a Bahcall–Wolf cusp at the Galactic center. We do not find support for this hypothesis neither in our four-component N-body models nor in our two-component Fokker–Planck models. The latter showed that even a quite large BH population still allowed a Bahcall–Wolf cusp to form in the MS stars on a timescale of ∼T_r(r_m); the main effect of the BHs is to decrease the asymptotic slope of the cusp from ∼r^−7/4 to ∼r^−3/2.

A dense cluster of stellar-mass BHs has been invoked as a potential solution to a number of problems of collisional dynamics at the Galactic center. Examples include randomization of the orbits of young stars via gravitational scattering (e.g., Perets et al. 2009), production of hypervelocity stars through encounters with BHs (e.g., O'Leary & Loeb 2008), and warping of young stellar disks (Kocsis & Tremaine 2011). These treatments typically assume a collisionally relaxed state for the Galactic center. In the relaxed models, the mass in BHs inside 0.1 pc is ∼10⁴ M_☉, i.e., N_BH(r < 0.1 pc) ≈ 10³ and N_BH(r < 0.01 pc) ≈ 10² (Hopman & Alexander 2006; Freitag et al. 2006), assuming "standard" IMFs. A high density of stellar BHs at the centers of galaxies like the Milky Way is also commonly assumed in discussions of the EMRI problem (Amaro-Seoane et al. 2007).

Models like these are called into question by the lack of a Bahcall–Wolf cusp in the late-type stars at the center of the Milky Way (Buchholz et al. 2009; Do et al. 2009; Bartko et al. 2010). If the Galactic center is not collisionally relaxed, as these observations seem to suggest, then computing the distribution of the heavy remnants becomes a more difficult, time-dependent problem (Merritt 2010), and knowledge of the initial conditions is essential.

Dissipationless mergers imply "initial conditions" characterized by a core in the dominant stellar component. The cores formed in our (major) merger simulations are large enough that they do not evolve appreciably (i.e., shrink) even after ∼3 relaxation times at the MBH influence radius. As a result, evolution of the BH distribution takes place against a stellar background with a very different density profile than in the steady-state models. In a core around an MBH, dynamical friction is much weaker than one would estimate by plugging the local density into standard formulae for orbital decay, due to the absence of stars moving more slowly than the local circular velocity (Antonini & Merritt 2011). Scaling our N-body models to the Milky Way, we predicted numbers of BHs inside r_h that were substantially smaller than in the collisionally relaxed models; in the inner 10⁻² pc—the radii most relevant to the EMRI problem (Merritt et al. 2010)—the number of BHs was, at most, 10–100 times smaller than predicted by these models, even after 10 Gyr (Figure 20).

It is unclear how relevant merger models are to the center of the Milky Way. If the core observed at the Galactic center had some other origin, the distribution of stellar BHs might have little connection with the distribution of the giant stars. However, if cores of size r_c ≈ r_h are common features of galactic nuclei, and if at some early time both the BHs and the stars had a common core radius, our models imply that the distribution of BHs should be considered very uncertain, even in galaxies for which nuclear half-mass relaxation times are as short as the age of the universe.