{'fraction_non_alphanumeric': 0.038393930983439725, 'fraction_numerical': 0.02331390507910075, 'mean_word_length': 4.220719632938904, 'pattern_counts': {'":': 1, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Merger of Binary Supermassive Black Holes in Merging Galaxies

In spherical galaxies, binary supermassive black holes (SMBHs) have difficulty reaching sub-parsec separations due to depletion of stars on orbits that intersect the massive binary - the final-parsec problem. Galaxies that form via major mergers are substantially nonspherical, and it has been argued that the centrophilic orbits in triaxial galaxies might provide stars to the massive binary at a high enough rate to avoid stalling. Here we test that idea by carrying out fully self-consistent merger simulations of galaxies containing central SMBHs. We find hardening rates of the massive binaries that are indeed much higher than in spherical models, and essentially independent of the number of particles used in the simulations. Binary eccentricities remain high throughout the simulations. Our results constitute a fully stellar-dynamical solution to the final-parsec problem and imply a potentially high rate of events for low-frequency gravitational wave detectors like LISA.


INTRODUCTION
The prospect that low-frequency gravitational radiation will be detected by LISA, the laser-interferometer space antenna (Hughes 2003;Barack & Cutler 2004), has motivated theoretical studies into the formation and evolution of binary supermassive black holes (SMBHs). Such binaries would constitute the highest signal-tonoise ratio sources of low-frequency gravitational waves, but as is the case for virtually all potential LISA sources, the event rate is still poorly known, with estimates ranging from a few to a few thousand events per year (e.g. Wyithe & Loeb 2003;Rhook & Wyithe 2005;Sesana 2010). While astronomical evidence for the existence of SMBHs with masses M • > 10 6 M ⊙ in galaxy spheroids beyond the Local Group is increasingly compelling (Ferrarese & Ford 2005), the evidence for SMBHs with masses ≤ 10 6 M ⊙ , which lie more nearly in the LISA band, is still largely circumstantial (van der Marel 2004), and the evidence for binary SMBHs is weaker still (Komossa 2006). When two galaxies merge, a binary SMBH forms at the center of the new galaxy (Begelman et al. 1980;Roos 1981). A common practice when estimating event rates for LISA is to equate the binary SMBH coalescence rate with the galaxy merger rate, the latter derived from models of structure formation in which galaxies merge hierarchically (Haehnelt 1994;Menou, Haiman & Narayanan 2001;Volonteri, Haardt & Madau 2003;Wyithe & Loeb 2003;Jaffe & Backer 2003). However it is also possible that coalescence is delayed, perhaps indefinitely, if the massive binary is limited in its ability to exchange angular momentum with surrounding stars and gas. This is sometimes referred to as the "final parsec problem" (e.g. Milosavljevic & Merritt 2003) since at this approximate separation, a massive binary in a spherical galaxy would efficiently eject all stars on intersecting orbits Berczik et al. 2005). Further evolution of the binary is possible, even in the absence of gas, if these orbits are replenished, or if the galaxy geometry allows for a much larger population of centrophilic orbits than in the spherical case (Merritt & Poon 2004;Berczik et al. 2005).
In this paper, we argue that the latter situation may be generic. We carry out direct N -body merger simulations of galaxies containing central SMBHs, and show that the binary SMBHs formed in such simulations continue to shrink, via interaction with stars, well beyond the point at which they would be expected to stall in spherical galaxies. Of course, this mechanism can act in addition to torques from gas, if present, (e.g. Escala et al. 2005;Dotti et al. 2007;Cuadra et al. 2009), or to loss-cone repopulation due to "massive perturbers" (Perets & Alexander 2008), etc. However it appears that collisionless stellar dynamics is sufficient.
In §2 we describe the numerical methods and initial conditions. Simulations of isolated (spherical) models are presented in §3, while §4 describes the results of galaxy merger simulations. §5 estimates time scales for coalescence of massive binaries, and §6 sums up.

NUMERICAL METHODS AND INITIAL CONDITIONS
Our initial conditions are based on spherical galaxy models following Dehnen's (1993) density law: and with cumulative mass profile where M gal is the total mass of the galaxy, a is a scale radius, and γ defines the inner density profile slope. We adopted γ = 1, corresponding to a Hernquist model (Hernquist 1990), for all initial models. A massive particle representing a SMBH was placed at the center of each galaxy model. The isotropic distribution function reproducing ρ(r) in the combined potential of the stars and the SMBH in dynamical equilibrium was computed numerically and used to generate Monte-Carlo positions and velocities for the stars. In what follows, we adopt "model units" such that M gal = G = r 0 = 1. The models can be scaled to physical units by fixing galaxy mass M gal and scale radius r 0 combined with the relations [V ] = GM gal r 0 In a Hernquist model the effective (projected half-mass) radius r e is ∼ 1.81r 0 . For a typical, luminous (M B ≈ −20) elliptical galaxy or bulge, r e ≈ 1.5 kpc, M gal ≈ 10 11 M ⊙ and [T ] ≈ 1.1 Myr. Effective radii scale very weakly with galaxy luminosity (though with considerable scatter; e.g. Ferrarese et al. (2006)) and in the case of 10 9 M ⊙ galaxy, whose SMBH would have a mass more suited to detection by LISA, the physical unit of time would be closer to 6Myr. In the latter case, the length of our longest integrations amounts to ∼ 1.5 Gyr.
Two sets of simulations were carried out. In the first set, a single spherical galaxy was created, with a central SMBH of mass M • = 0.001 (Model A) or 0.01 (Model B). (In what follows, M • is always used to denote the mass of a single SMBH.) A second SMBH, of the same mass, was then placed at a distance of 0.5 from the center on a circular orbit. Integrations were continued until a time t max = 100.
In the second set, two identical spherical models were created, each containing a single central SMBH, and the two galaxies were merged. Here again, two values for M • were used: M • = 0.001 (Model C) or 0.01 (Model D) and the integrations were continued until a time t max = 250. In order to test the dependence of the results on particle number N , all simulations were carried out with N = (32, 64, 128, 256, 512) × 10 3 . Tables 1 and 2 summarize the parameters of the two sets of runs.
The N -body integrations were carried out using φGRAPE (Harfst et al. 2007), a parallel, directsummation N -body code that uses special-purpose hardware to compute pairwise forces between all particles: Here m i and r i are the mass and position of the ith particle respectively and ǫ is a force softening parameter. We set ǫ = 10 −5 in all the runs. φGRAPE integrates the equations of motion using a fourth-order Hermite integrator with individual block time steps similar to the well established Aarseth direct N-body codes (Makino & Aarseth 1992). The time step of particle i at time t was computed from the standard formula: where a is the acceleration of ith particle and the superscript k indicates the order of the time derivative. The time-step parameter η was set to 0.01 for all our runs. The choice of parameters was motivated by Berczik et al. (2005) who found that integration accuracy is more sensitive to η than to ǫ. Those authors adopted η = 0.01 and ǫ = 1 × 10 −4 for their study of SMBH binary evolution. In order to resove the orbit of the massive binary in our simulations down to a semi-major axis a ≈ 10 −4 we adopted the smaller softening length ǫ = 10 −5 . With these choices of η and ǫ, the relative energy error in our simulations was ∼ 10 −8 . For a more detailed discussion of these issues we refer the reader to section 2 of Merritt et al. (2007) . The N -body integrations were carried out on three, special-purpose computer clusters.
Titan, at the Astronomisches Rechen-Institut in Heidelberg, and gravitySimulator, at the Rochester Institute of Technology are both 32-node clusters employing GRAPE accelerator boards (Harfst et al. 2007). Some calculations were also carried out on the GPU-enhanced cluster "Kolob" at the University of Heidelberg. For the latter runs we used the modified version of the φGRAPE code including the SAPPORO library (Gaburov et al. 2009).

ISOLATED MODELS
Evolution of a binary SMBH in a gas-free galaxy can be divided into three phases (Begelman et al. 1980). In first phase, the SMBHs are unbound to each other and move independently in the galaxy potential. The evolution of the individual SMBH orbits during this phase is governed by dynamical friction against the background stars. The second phase begins when the two SMBHs are close enough together to form a bound pair. In this phase, the separation between the SMBHs decreases very rapidly due to the combined effects of dynamical friction and ejection of stars by the gravitational slingshot. (For a detailed discussion of this phase, we refer the reader to the careful analysis in Milosavljević & Merritt (2001).) At the end of the second phase, the two SMBHs form a hard binary, defined as a binary that ejects passing stars with positive (unbound) energies. Phase three consists of slow evolution of the binary as depleted orbits are repopulated, via two-body scattering or some other process. In a spherical galaxy with large N , the timescale for orbital re-population is essentially the two-body relaxation time which scales as ∼ N/ ln N . This is the origin of the "final-parsec problem" (Merritt & Milosavljevic 2005). Figure 1 plots the separation R between the two SMBH particles versus time in the isolated galaxies, Models A and B. Here we can see the three phases of binary evolution described above.
(1) In the first phase, the separation between the two SMBHs decreases due to dynamical friction. This phase ends when R ≈ r h , the gravitational influence radius. r h can be defined as the radius of a sphere around the binary SMBH that encloses a stellar mass equal to twice the binary mass. For Models A, M • = 0.001 and r h ≈ 0.07. From Figure 1, we can see that the first phase ends indeed when R ≈ 0.07. Here we want to emphasize that this estimate does not take into account the changes in the galaxy that are induced by the presence of the second SMBH.
For Models B, the radius of influence is ∼ 0.25. Judging from the bottom panel of Figure 1, it appears that the first phase ends approximately when R ≈ 0.25, although the transition is not as sharp as in model A, presumably because the initial separation is only ∼ 2r h .
(2) In the second phase, the separation between the two SMBHs decreases very rapidly. Dynamical friction, and the ejection of stars by the gravitational slingshot, act together to efficiently extract angular momentum from the massive binary. These processes are not well defined in the regime where the massive binary is neither very hard nor very soft; equations (13) and (14) of Milosavljević & Merritt (2001) give approximate expressions for the rate at which energy is transferred to stars by the two mechanisms. For Models A this process continues until t ≈ 40 while for Models B, this phase ends at t ≈ 20. The motion of two SMBHs in this phase is approximately Keplerian. We define the semi-major axis a and eccentricity e via the standard relations: where µ = G(M •1 + M •2) , v is the relative speed, and h is the specific angular momentum of the relative motion.
(3) The rapid phase of binary hardening comes to an end when a ≈ a h , where a h is the semi-major axis of a "hard binary": (e.g. Merritt et al. 2007). Here q = M •2 /M •1 ) is the mass ratio of two SMBHs; for both Models A and B, q = 1. For Models A, a h ≈ 0.004 with a −1 h ≈ 250 and for Models B, a h ≈ 0.016 with a −1 h ≈ 65. From Figure 2, we see that there is no clear dependence of the binary hardening rate on N until the binary becomes hard, a < a h . For Models A this happens around t = 25 and for Models B around t = 30. Beyond these times the binary hardening rate exhibits a clear N dependence in the sense that hardening is slower for larger N . Similar N -dependence has been seen in other studies of binary evolution in spherical galaxies (e.g. Berczik et al. 2006;Merritt et al. 2007). For real spherical galaxies, with much larger N , the binary would stop evolving beyond this point.
We estimated the hardening rate s in the N -body integrations, where by fitting straight lines to a −1 (t) in an interval ∆t = 50 from t = 50 to t = 100. Figure 3 shows the Ndependence of s. We see that for N > 50k, the hardening rate is a decreasing function of N ; scaling as ∝ N −0.6 ; this is consistent with the scalings found by other authors in the same N -range, e.g. Berczik et al. (2005); Merritt et al. (2007). Comparison of Figures 3a and b also shows lower hardening rates for Model B compared with Model A, by roughly an order of magnitude. For large N , the two-body relaxation time is long compared with orbital periods of stars near the massive binary, and the loss cone around the binary is nearly empty. In this regime, and in the asymptotic (large-N ) limit, the rate of binary evolution is predicted to scale as ∼ [(M •1 + M •2 )] −1 (e.g. equation (32) of Merritt et al. (2007)). Given that a is roughly a factor ten smaller, at a given time, in Models A as compared with Models B, there is an additional delay in the early phase of the hard binary evolution. However this delay is negligible for the decay time t 0 until which the semimajor axis reaches a 0 , when energy loss by GWs start to dominate. In section 5 a more detailed discussion of the total time to merge the binary SMBH is given. Figure 5 shows the density profile for model A5 in the final phase. The central logarithmic slope of the density profile has dropped from −1 to ∼ −0.6.

GALAXY MERGERS
We studied the evolution of binary SMBHs in merging galaxies by creating two identical galaxy models, each with a central SMBH, and placing the two galaxies on  bound relative orbits. Table 2 summarizes the parameters of the galaxy merger models. The initial separation of the two galaxy centers was 20r 0 . The initial relative velocity of the two galaxies was chosen such that the SMBH separation at first pericenter passage was ∼ r 0 ; in other words, the initial orbit of the binary galaxy was substantially eccentric. This was done in order to reduce the computational time required to bring the two SMBHs close together; in fact, a galaxy merger from such nearly "head-on" initial conditions is unlikely. Figure 4 illustrates the merging of two galaxies. In the initial phase (top panels) the two SMBH particles remain strongly associated with their respective density cusps. The two cusps merge into one at t ∼ 100, following which the cusp is gradually destroyed by the gravitational slingshot ejection of stars (Fig. 5).  Figure 6 shows the separation between the two SMBHs in the merger models. The first pericenter passage occurs around t ≈ 80. As the galaxies merge, the two SMBHs remain centrally located in their respective density cusps, until coming close enough together that a binary SMBH forms, at t ≈ 100, or slightly earlier in the case of Model D (phase one). The separation between the two SMBHs then decreases very rapidly due to the combined effects of dynamical friction and slingshot ejection of stars (phase two). Once the hard binary is formed, the hardening rate of the binary decreases (phase three). This behavior is qualitatively similar to what was seen in the isolated galaxy models.
However there is one important difference between the binary evolution in the isolated and merging galaxies. Figure 7 shows the evolution of the binary semi-major axis, and Figure 8 shows the binary hardening rates in the galaxy merger simulations; the latter were computed in the same way as in the isolated galaxy models, by linear fit of 1/a in the range 150 ≤ t ≤ 200. In the galaxy merger models, there is essentially no dependence of the binary hardening rate on particle number. Furthermore, in Models C and D, the hardening rates are much higher -more than five times higher than in the Models A and B, respectively, when N = 512k.
A likely explanation for this difference is that the nonspherical shapes of the merger remnants provide an additional source of torque on the stellar orbits, allowing them to interact with the central binary on a timescale shorter than the two-body relaxation time. It has been argued (Merritt & Poon 2004;Berczik et al. 2006) that allocating even a small fraction of the stellar orbits to a "centrophilic" family could lead to binary hardening rates that are much larger than those that result from two-body relaxation alone.
We evaluated the shapes of the merged galaxies in a number of ways. The central density contours for Model C3 are shown in Figure 9. Departures from axisymmetry are evident. Figure 10 shows the principle axis ratios of the merger remnant in Model C3 as a function of time; these were defined as the axis ratios of a homogeneous ellipsoid with the same inertia tensor. The departures from spherical symmetry are modest, but definite, and they appear to be nearly independent of time toward the end of the simulation.
The evolution of binary eccentricity in the merger models is presented in Figure 11. We notice high eccentricities for both models as soon as the SMBHs become bound, approaching unity for models D. The high eccentricities are due in part to the high eccentricity (e = 0.95) of the galaxies' relative orbit prior to the merger, and, in some runs, to post-merger evolution.

TIME SCALES FOR COALESCENCE
SMBH binaries are a potentially important source of gravitational wave (GW) emission. Gravitational waves extract energy and angular momentum from the binary, hence changing its orbital elements. Peters (1964) gives approximate, orbit-averaged expressions for the rates of change of a binary's semimajor axis and eccentricity due to GW emission: We define t GW (a 0 , e 0 ) as the time required, according to the coupled equations (10a, 10b), for the binary semimajor axis to shrink to zero under the influence of GW emission, starting from the initial values (a 0 , e 0 ). This time is a strong function of e 0 for e 0 ≈ 1, as it is in some of our simulations. The full time to coalescence, t coal , includes also the time from the start of the simulation until the GW regime is entered. We estimated t coal as follows. In the N -body merger models, the post-merger hardening rate of the binary, s, is essentially independent of a and N (Figures 7,  8). The rate of change of the binary semi-major axis in this regime is approximately Given values for M gal and r 0 , this dimensionless rate can be converted into physical units. At some time in the simulation, a will have fallen to such a small value that (da/dt) NB as given by equation (11) will equal (da/dt) GW , where the latter quantity is computed from both a and e. We define t 0 as the time at which these two rates are equal. The full time to coalescence, t coal , is the time from the start of the simulation until t 0 , plus t GW (a 0 , e 0 ), where a 0 ≡ a(t = t 0 ), e 0 ≡ e(t = t 0 ). In most cases, t 0 exceeded the final time t final of the simulation. In these cases, the time ∆t between t final and t 0 was assumed to be In such cases, we also assumed that e(t) was constant, e(t) = e(t final ), for t > t final when computing t GW . Computation of the GW evolution rates required adoption of particular values for the physical units of length  (t = 60, 68, 76, 84, 92). The figure also shows the initial density profile (γ = 1) and a fit to the data using Eq. 1 with γ = 0.6 and r 0 = 0.8. The central density drops as the binary ejects mass from the core. Right: Spatial density profile for model C4 averaged over five different time steps (t = 140, 160, 180, 200, 220). The figure also shows the initial density profile (γ = 1) and a fit to the data using equation 1 with γ = 0.3 and r 0 = 0.6. and time, hence of M gal and r 0 . Table 3 gives results for two different assumed values of the galaxy mass: M gal = (10 9 , 10 11 )M ⊙ .
Because of the strong dependence of t GW on eccentricity, Table 4 gives a GW and ∆t, for run C3 only, under two different assumptions about e 0 : e 0 = 0.5 and e 0 = 0. The total decay time increases by a factor of 2-3 if the orbit is assumed circular at the start of the GW phase. For low mass-galaxies (10 9 M ⊙ ), the total time to merge the SMBHs exceeds ∼ 1 Gyr, whereas for highmass galaxies (10 11 M ⊙ ) the coalescence time is well below 0.5 Gyr. This is contrary to the usual statement that binary coalescensce is slower in more massive binaries. That statement reflects the much longer timescales for gravitational encounters in massive galaxies; but in our merger models, evolution prior to the GW regime occurs  much more rapidly than predicted based on collisional loss cone repopulation rates, and the dominant influence on the coalescence time is the timescale for GW energy loss, which is longer for smaller SMBH masses. We note that all of the coalescence times in Table 3 are short compared with those expected in spherical models (e.g. Figure 19 of Merritt, Mikkola & Szell 2007).
The lower mass scaling corresponds to sources detectable by LISA. We find coalescence times of a few Gyr  in these cases, which, while relatively short, are nevertheless long enough that SMBH binaries in these galaxies may not achieve full coalescence before a subsequent galaxy merger occurs.