6 N-Body Studies of Binary Black Hole Evolution

6 $N$ -Body Studies of Binary Black Hole Evolution

The interaction of a massive binary with point perturbers at the center of a galaxy is a straightforward problem for $N$ -body simulation. In principle, $N$ -body studies can reveal both the long-term evolution of the binary, as well as the effect of the binary on its stellar surroundings. The latter can be compared with observed nuclear density profiles as a test of the theory (Section 7).

Unfortunately, unless great care is taken, $N$ -body studies are likely to give misleading results. This follows from the result (Section 4.3) that time scales for two-body scattering of stars into the binary’s loss cone are of order $1010 yr$ or somewhat longer in real galaxies. In $N$ -body simulations, relaxation times are shorter by factors of $~ N/1011$ than in real galaxies, hence the long-term evolution of the binary is likely to be dominated by spurious loss cone refilling, wandering of the binary, and other noise-driven effects.

$N$ -body studies are most useful at characterizing the early stages of binary formation and decay, or simulating the disruptive tidal effects of a SBH on the nucleus of an infalling galaxy. Due to algorithmic limitations - primarily the difficulty of integrating galaxy models with high central concentrations - most such studies [35 , 75 , 126 , 125 , 158 , 157 , 86 , 25 , 128 ] have been based on galaxy models with unrealistically large cores.

Figure 7 is from the first [136 ] $N$ -body simulation of galaxy mergers in which the pre-merger galaxies contained power-law nuclear cusps as well as massive particles representing the SBHs. These simulations were run using GADGET [202 ], a tree code with inter-particle softening, and were not able to accurately follow the formation and decay of the massive binary. The SBH in the larger galaxy was found to tidally disrupt the steep cusp in the infalling galaxy, producing a remnant with only slightly higher central density than that of the giant galaxy initially. This result helps to explain the absence of dense cusps in bright galaxies [55 ], and suggests that the central structure of galaxies can only be understood by taking into account the destructive influence of SBHs on the stellar distribution during mergers. Additional results, using a similar $N$ -body code and a variety of mass ratios for the merging galaxies, were reported in [141 ].

Figure 7:

Final density profiles from a set of 10:1 merger simulations in which each galaxy contained a black hole (a-d) and in which neither galaxy contained a black hole (e-g) [136 ]. The four thin curves in each frame correspond to four different pre-merger orbits. (a), (e) Space density of stars initially associated with the secondary galaxy; thick curves are the initial density profile. (b), (f) Space density of stars initially associated with the primary galaxy; thick curves are the initial density profile. (c), (g) Space density of all stars. Lower thick curves are the initial density profile of the primary galaxy, and upper thick curves are the superposition of the initial density profiles of the primary and secondary galaxies. Lines of logarithmic slope -1 and -2 are also shown. (d), (h) Logarithmic slope of the surface density profiles of the merger remnants. Thick curves correspond to the initial primary galaxy.

Had these simulations been extended to longer times using a more accurate $N$ -body code, the massive binary would have ejected stars via the gravitational slingshot and lowered the central density still more. This was first demonstrated [150

] in an $N$ -body study that used a tree code for the early stages of the merger, and NBODY6, a high-precision, direct-summation code [1 ], for the later stages, when the binary separation fell below the tree code’s softening length. The pre-merger galaxies had steep, $-2 r ~ r$ density cusps and the mass ratio was 1:1. These simulations were continued until the binary separation had decayed by a factor of $~$ 10 below $ah$ . The initially steep nuclear cusps were converted to shallower, $r ~ r-1$ profiles shortly after the SBH particles had formed a hard binary; thereafter the nuclear profile evolved slowly toward even shallower slopes as the massive binary ejected stars. As Figure 7

illustrates, the stellar density around the binary drops very quickly after the binary reaches a separation $a ~~ ah$ .

Figure 8:

Lagrangian radii around each of the two SBH particles in an equal-mass merger simulation [150

]. From bottom to top, the radii enclose $-4 10$ , $-3.5 10$ , $-3 10$ , $- 2.5 10$ , $-2 10$ , $-1.5 10$ and $10- 1$ in units of the mass of one galaxy before the merger. The binary becomes “hard” at $t ~~ 11$ , and very rapidly heats the surrounding stellar fluid, lowering the local density.

The hardening rate of the binary in these simulations was found not to be strongly dependent on the number of particles. This result was subsequently shown [152

] to be due to the small $N$ : Stars were resupplied to the loss cone via collisions at a higher rate than they were being kicked out by the binary, ensuring a continuous supply of stars and allowing the binary to continue to shrink. While a qualitatively similar evolution may take place in some galaxies - for instance, loss cones in non-axisymmetric potentials can be continuously repopulated by stars on centrophilic orbits (Section 4.4) - collisional loss cone refilling is very unlikely to achieve anything like a full loss cone except in very small, dense galaxies. The long-term evolution of the binary in almost all published $N$ -body simulations are therefore not representative of what one would expect in real galaxies.

Figure 9:

Evolution of the binary semi-major axis (a) and hardening rate (b) in a set of high accuracy $N$ -body simulations; the initial galaxy model was a low-central-density Plummer sphere [20

]. Units are $G = Mgal = 1$ , $E = - 1/4$ , with $E$ the total energy. (a) Dashed lines are simulations with binary mass $M1 = M2 = 0.005$ and solid lines are for $M1 = M2 = 0.02$ , in units where the total galaxy mass is one. (b) Filled (open) circles are for $M = M = 0.005 1 2$ $(0.02)$ . Crosses indicate the hardening rate predicted by a simple model in which the supply of stars to the binary is limited by the rate at which they can be scattered into the binary’s influence sphere by gravitational encounters. The simulations with largest ( $M1, M2$ ) exhibit the nearly $N - 1$ dependence expected in the “empty loss cone” regime that is characteristic of real galaxies.

The relevant dimensionless parameter is

where $dJ$ is the change over one radial period in the angular momentum of a star on a low- $J$ orbit, and $Jlc$ is the angular momentum of an orbit at the edge of the binary’s loss cone. A value $q(E) » 1$ implies that the loss cone orbits at energy $E$ are re-populated at a much higher rate than they are de-populated by the binary, and the loss cone remains nearly full. A value $q« 1$ implies that the loss cone is essentially empty, and repopulation must take place diffusively, as stars scatter in from $J >~~ Jlc$ . In real galaxies, $N$ is large and $dJ$ is small, implying $q « 1$ . Achieving $q « 1$ in $N$ -body simulations requires large particle numbers, and/or a model for the galaxy that has an unrealistically low central density, so that the star-star relaxation time is long. Figure 9

shows a recent set of $N$ -body simulations that does both [20 ]. The massive binary was embedded in a Plummer [168 ] galaxy model, which has a core radius comparable to its half-mass radius; this model is very different from real galaxies but its very low degree of central concentration implies a long relaxation time and low rate of collisional loss-cone refilling. Large particle numbers were achieved, without sacrificing accuracy, by running the simulation on a parallel GRAPE cluster. The $N$ -dependence of the binary’s hardening rate is clear; in the simulations with binary mass $M1 = M2 = 0.02Mgal$ , the $N$ -dependence of the hardening rate is $s =_ (d/dt)(1/a) ~~ N -0.8$ , almost as steep as the $N -1$ dependence predicted for an “empty” loss cone [152

]. A similar study [128

], based on King-model galaxies, found a similar result.

The Plummer-model initial conditions used in the simulations of Figure 9 were identical to those adopted in two other $N$ -body studies based on a more approximate $N$ -body code [178 , 25 ]. Contrary to Figure 9 , Chatterjee, Hernquist & Loeb (2003) found that the binary hardening rate “saturated” at values of $5 N >~~ 2 × 10$ , remaining constant up to $5 N ~~ 4× 10$ . They speculated that this was due to a kind of Brownian-motion-mediated feedback, in which the binary maintains a constant supply rate by modulating the local density of stars. However no supporting evidence for this model was presented; for instance, it was not demonstrated that the central density was actually regulated by the binary, or that the amplitude of the Brownian wandering increased with $N$ in the manner postulated. Furthermore these authors provided no plots showing the claimed $N$ -dependence of the hardening rate. Chatterjee et al.’s conclusion, that “a substantial fraction of all massive binaries in galaxies can coalesce within a Hubble time”, is not substantiated by the more accurate $N$ -body simulations shown in Figure 9 .

While Brownian motion probably does affect the decay rate of binaries in $N$ -body simulations [152 ], it is doubtful that the effect is significant in real galaxies. The Brownian velocity of single black holes is found in $N$ -body integrations to be [114 ]

where $2 ~s$ is the 1D, mean square stellar velocity within a region $r <~ 0.5rh$ around the black hole (and includes the influence of the black hole on the stellar motions), and $m*$ is the stellar mass. In the case of the Milky Way black hole, Equation (42

) implies $Vrms ~~ 0.2km s-1$ (assuming $m = M * o .$ ) and an rms displacement of $<~ 0.1pc$ . Brownian motion of a massive binary is larger than that of a single black hole, but only by a modest factor [133

, 150

, 128 ]. The rms displacement of a binary from its otherwise central location would therefore be very small in a real galaxy, probably even less than the separation between the two components of the binary.

Figure 10:

Results from a set of $N$ -body integrations of a massive binary in a galaxy with a $r ~ r -0.5$ density cusp [205 ]. Each curve is the average of a set of integrations starting from different random realizations of the same initial conditions. (a) Evolution of the “mass deficit” (Equation 43

), i.e. the mass in stars ejected by the binary. For a given value of binary separation $a$ , the mass deficit is nearly independent of particle number $N$ , implying that one can draw conclusions from observed mass deficits about the binary that produced them. (b) Evolution of binary eccentricity. The eccentricity evolution is strongly $N$ -dependent and tends to decrease with increasing $N$ , suggesting that the eccentricity evolution in real binaries would be modest.

The goal of $N$ -body studies is to simulate binary evolution in galaxies with realistic density profiles, and with large enough particle numbers that re-population of the binary’s loss cone takes place diffusively, as in real galaxies. Two avenues are open for making further progress in this area.

One could combine a Monte-Carlo treatment of stellar encounters [57 , 58 ] with a lookup table, derived from scattering experiments, of energy and angular momentum changes experienced during close passages of stars to the binary [178 , 133 ]. Such a hybrid algorithm would allow one to adjust the degree of collisionality at will and record the effects on both the binary’s evolution, and the influence of the binary on the stellar distribution. This approach would be difficult to generalize to non-spherical geometries however.
A straightforward $N$ -body approach is also feasible, but particle numbers in excess of $7 ~ 10$ are required [152 ]. Such large particle numbers are just now becoming feasible for direct-summation $N$ -body codes, by combining the GRAPE accelerator boards [127 ] with a parallel architecture [33 ]. Indeed just such an approach was used for the integrations of Figure 9 , although the galaxy models in that study were rather unphysical.

Figure 10 shows a promising early step in this direction. The initial galaxy models had $-0.5 r ~ r$ density cusps; integrations were carried out on a GRAPE-6 computer, limiting the total particle number to $~$ 256 K, but the motion of the black hole binary (of mass $M1 = M2 = 0.005Mgal$ ) and nearby stars was carried out using the Mikkola-Aarseth chain regularization algorithm [145 , 146 , 3 ]. Because of the models’ higher central density and limited particle numbers, the binary’s loss cone was only partially empty, $q >~~ 1$ , and the $N$ -dependence of the hardening rate was shallower than expected in real galaxies, $- 0.4 (d/dt)(1/a) ~ N$ . Figure 10 (a) shows that the “damage” inflicted by the binary on the nucleus is not strongly dependent on $N$ , as expected (cf. Equation 18 ). This is an encouraging result since it implies that one can hope to learn something definite about pre-existing binaries by comparing $N$ -body simulations with observations of the centers of current-day galaxies (Section 7). On the other hand, Figure 10 (b) suggests that the evolution of the binary’s eccentricity is strongly $N$ -dependent. This may explain the rather disparate results on eccentricity evolution in past $N$ -body studies [150 , 86 , 2 ].

Much progress on this problem is expected in the next few years.

6 -Body Studies of Binary Black Hole Evolution

6 $N$ -Body Studies of Binary Black Hole Evolution