# Fast coalescence of massive black hole binaries from mergers of galactic nuclei: implications for low-frequency gravitational-wave astrophysics

## Abstract

We investigate a purely stellar dynamical solution to the Final Parsec Problem. Galactic nuclei resulting from major mergers are not spherical, but show some degree of triaxiality. With -body simulations, we show that massive black hole binaries (MBHB) hosted by them will continuously interact with stars on centrophilic orbits and will thus inspiral—in much less than a Hubble time—down to separations at which gravitational wave (GW) emission is strong enough to drive them to coalescence. Such coalescences will be important sources of GWs for future space-borne detectors such as the Laser Interferometer Space Antenna (LISA). Based on our results, we expect that LISA will see between to such events every year, depending on the particular MBH seed model as obtained in recent studies of merger trees of galaxy and MBH co-evolution. Orbital eccentricities in the LISA band will be clearly distinguishable from zero with .

###### Subject headings:

black hole physics — galaxies: nuclei — stellar dynamics — gravitational waves## 1. Introduction

Massive black hole binaries (MBHBs) are one of the most interesting sources of gravitational waves (GWs) for future space-borne detectors such as the Laser Interferometer Space Antenna (LISA). They are expected to coalesce under the strong emission of GWs, after stellar- and/or gas-dynamical processes bring them to separations small enough ( pc) that GW emission is efficient in making them coalesce in less than a Hubble time (Milosavljevíc & Merritt 2003; Armitage & Natarajan 2005). It is still an open problem whether MBHB coalescences are generic and prompt, or whether long-lived binaries are the norm.

The paradigm for MBH binary evolution, after a merger of gas-poor galaxies, consists of three distinct phases (Begelman, Blandford & Rees 1980). First, the two MBHs sink towards the center due to the dynamical friction exerted by the stars. This process continues after they form a bound pair at a semimajor axis separation , where is the binary’s influence radius defined to be the radius which encloses twice the mass of the binary in stars. It stops when the binary reaches the hard binary separation (Quinlan 1996; Yu 2002)

(1) |

where is the binary’s reduced mass, is the local D velocity dispersion, is the binary’s mass ratio. Secondly, for , as dynamical friction becomes inefficient in further driving the inspiral, it is instead the slingshot ejection of stars, following three-body scattering with the binary, that dominates. Thirdly, the binary eventually reaches a separation at which the loss of orbital energy to GW emission drives the final coalescence. The transition from the first to the second phase is prompt provided that the mass ratio of the remnants is not too small (Colpi & Dotti 2009; Callegari et al. 2011). In contrast, the subsequent transition from the second to the third phase could constitute a bottleneck for the binary evolution towards final coalescence. This is the so-called Final Parsec Problem.

In quasi-steady spherical stellar environments, the binary’s hardening rate slows down significantly once it reaches separations a few times below
(Quinlan & Hernquist 1997; Milosavljevíc & Merritt 2003; Berczik et al. 2005). In these spherical and
gas-poor nuclei, two-body relaxation is the only mechanism for populating the binary’s loss
cone ^{11}

But spherical models are a worst case scenario—and not a very realistic one at that! Merger remnants will generally be irregular with some degree of triaxiality and, even if triaxiality would only be a rather mild and transient phenomenon, it may suffice to bring the binary down to (Merritt & Poon 2004). Berczik et al. (2006) and Berentzen et al. (2009) studied triaxial, rotating models of galactic nuclei using -body simulations—including the full post-Newtonian corrections to the MBHB. They have shown that MBHBs in such models do indeed coalesce in much less than a Hubble time. The next logical step is to study mergers of galactic nuclei to investigate if the latter results still hold true under more realistic models and initial conditions.

In this Letter, we use N-body simulations to show that: (1) in merging nuclei, the hardening rate is -independent—allowing the extrapolation of -body results to real galaxies; (2) the triaxiality depends on the orbital parameters of the progenitor galaxies: prolate shapes occur when the merger is almost radial, while an oblate morphology is the result of a less radial merger; (3) MBHs become bound with high eccentricities (up to ); (4) the eccentricity tends, on average, to increase in good agreement—often quantitative—with Quinlan (1996) predictions; (5) high eccentricities assist the MBHB into promptly coalescencing in much less than a Hubble time; (6) eccentricities in the LISA band are likely to be distinguishable from zero () even though GW circularizes the orbits, and will also be quite large () in the Pulsar Timing Array (PTA) band.

## 2. Models and Initial Conditions

We have performed two sets of -body experiments. In both sets, galactic nuclei are represented by spherically symmetric Dehnen models (Dehnen 1993; Tremaine et al. 1994). These models have a central power law density profile, , with logarithmic slope and a break radius which are both set equal to one. The total mass of each nucleus is set , and we adopt units where . The total mass of the binary MBH , and we take . We study unequal mass MBH coalescences in parallel papers (Berczik et al. 2011; Preto et al. 2011).

The set (A) consists of a single spherical nucleus where two MBHs are placed symmetrically about the center, on an unbound orbit, with initial separation , initial angular momentum , where is the angular momentum of the local circular orbit. The set (B) consists on the equal-mass merger of two initially bound—but well-separated—spherical nuclei, each of which has a single MBH at the center with zero initial velocity with respect to its nucleus. For B, the initial separation refers to both nuclei taken as if they were point masses located at each center of mass. The half-mass radius of each nucleus is ; accordingly, and in order to have an initial configuration with two well separated nuclei, while minimizing the computing time, we set the initial separation equal to . For the initial orbital angular momentum of the binary nuclei, we have taken two values and given the nearly-parabolic encounters typical of major galaxy mergers seen in cosmological simulations (Khochfar & Burkert 2006). During the first pericenter passages, the MBH separations are and , respectively. Table 1 lists the runs and adopted parameters.

64K | 128K | 256K | 512K | 1M | |
---|---|---|---|---|---|

\toprule0.005 | |||||

0.1 | — | ||||

\toprule0.005 | |||||

0.01 | — | ||||

0.02 | — | ||||

0.1 | — |

We have performed the -body simulations using the parallel -GPU code. This is a yet unpublished variant of the parallel direct N-body code -GRAPE (Harfst et al. 2007), which uses GPU accelerator cards on parallel clusters (Berczik et al. 2011). It includes a fourth-order Hermite integration scheme, with block time steps, analogous to NBODY1 (Aarseth 2003).

The code does not include regularization of close encounters, and softening of the gravitational interaction is adopted instead. The softening length has to be chosen small enough that it reproduces the refilling of the binary’s loss cone by two-body relaxation. After some testing, we adopt a softening length in model units. We set the time step parameter (Aarseth 2003) to for the field stars and for the BHs. Furthermore, we force the MBHs to be advanced synchronously at all times with the smallest step. With the parallelized version of the -GPU code, one can study models with very large number of particles and the results agree with NBODY4 (Aarseth 2003) as far as single stars and distant encounters are concerned. For the high velocity dispersions present in nuclei with a MBH, the effect of close encounters between field stars is negligible for the bulk evolution of the stellar system (Preto & Amaro-Seoane 2010).

## 3. MBH evolution in spherical versus in merging nuclei

The stars that drive the orbital decay of a hard MBHB are those that enter the loss cone orbits. The MBHB’s hardening rate is thus determined by the product of the flux of stars entering the loss cone with the average kinetic energy they receive when ejected—at the expense of the MBHB’s orbital energy—through the slingshot mechanism. Denoting by the time-dependent flux into the loss cone and by the mean kinetic energy imparted to stars which are scattered off by the binary, the hardening rate is given by (Yu 2002)

(2) |

where , and is the gravitational potential due to the stars. The mean kinetic energy is given by

(3) |

where is a dimensionless quantity which was measured from three-body scattering experiments (Quinlan 1996). Therefore, the hardening rate can be rewritten as

(4) |

The time evolution of the flux depends on the symmetries of the gravitational potential—and on the orbit families it supports—of the nuclei in question. In principle, in the spherical case can be obtained from Fokker-Planck calculations that take into account the diffusion of stars in phase space (Merritt et al. 2007; Preto & Amaro-Seoane 2010). Here we derive simple scaling relations which are useful in interpreting the N-body results. For each energy , where is the number of stars of energy per unit energy and is the local two-body relaxation time; the latter scales as (Spitzer 1987). The flux of stars into the loss cone is expected to peak around (Perets & Alexander 2008), so we evaluate these quantities there. Hence, —where follows from the relation (Ferrarese & Ford 2005). Then, obtains. On the other hand, for a fixed galaxy mass, we have and therefore . Since in our N-body models, and are unchanged and only changes as is varied, we find that the hardening rate scales with and as . The case of a triaxial nucleus is different: for each star is not conserved, thus stars may precess into the loss cone on a time scale (Merritt & Poon 2004); and will depend only on the global gravitational potential of the galaxy. In this case, the mass flux into the loss cone , and also , will be independent of the number of stars.

In Figure 1, we see that is N-dependent in a spherical nucleus, while it is -independent in the merging one. In the former case, , with and for binaries of and . These results can be interpreted as follows. In the empty loss cone limit (), the stars repopulate the loss cone at a rate much lower than the that with which they are ejected by the binary, which is . In the full loss cone limit (), stars enter the loss cone at a rate which is similar to the rate at which they are ejected by the binary. A measure for the loss cone refilling rate is given by (Lightman & Shapiro 1977)

(5) |

where is the mean change in , per orbital period, of a star on a low-.

In the limit when , the loss cone is said to be empty; while in the full loss cone limit. For a given nucleus, and for , we expect to be independent of . As a result, the weaker dependence of on for lighter binaries, placed in a spherical nucleus, follows from ; at the same , of the binary is larger than that of the one. We would need to use a larger for the lighter binaries, , to enter deep into the empty loss cone limit ; the heavier binary, , almost reaches this limit. The dependence of on is more straightforward to interpret. For the spherical case, the lighter binary is expected to harden at a rate higher than the heavier, which is indeed the case. In the merger case, is -independent and therefore . Since the mass ratio between the binaries is , also differs by a corresponding factor of two.

Following Merritt & Poon (2004), we measure the triaxiality of the nucleus with
. ^{12}

## 4. Eccentricity evolution and time scales for coalescence

The hardening rate, due to the slingshot ejection of stars, is found to be in our -body simulations essentially independent of the binary’s orbital elements, while that due to GWs is strongly dependent on them: (Peters 1964). As a result, the time a binary takes to coalesce depends strongly on its eccentricity. In paper I, we did follow this evolution self-consistently with N-body simulations of rotating King models. Since such calculations are extremely CPU-intensive, we estimate the full evolution using a semi-analytic approach (Quinlan 1996). The advantage is that we can calibrate the average hardening rate with our N-body simulations—which would remain a free parameter otherwise—, and thus make quantitative predictions on both the coalescence times and the long term eccentricity evolution.

The evolution of the MBHB orbital elements, including the effect due to orbital energy lost to GWs, is given by

(6) |

The GW terms are given in Peters (1964). The eccentricity evolution, driven by the stars, is obtained from three-body scattering experiments

(7) |

where and the constants are taken from Quinlan (1996). In order to assess the quality of the fits to the -body results, Figure 3 compares the -body evolution of the binary with that obtained from the semi-analytic model. We take as initial conditions for the integration of equations (4) an instant of time in the early hard binary phase. Given the differences between our -body models and the assumptions embodied by the semi-analytic description the agreement is quite remarkable.

We then include the GW terms due to radiation reaction to compute the time it takes for the binary to coalesce. To scale our models to binaries with and , we adopt the most recent observational values for the mass normalization of the Milky Way nucleus, (Schödel et al. 2009) and use the relation to extrapolate to different MBH masses. The results are shown in the upper panel of Figure 4. We see that coalescence times range between yrs and yrs. These times are not longer that the mean time between successive major mergers. In contrast, for a spherical nucleus, coalescence times for the lower mass would become , while binaries with would stall (Preto et al. 2011).

We also follow the long term evolution of the eccentricity. In the -body runs, the binaries
become bound with high eccentricities (up to ) on average—in agreement with
previous works (Berentzen et al. 2009; Preto et al. 2009). Since LISA will be sensitive to
the inspiral signal of binaries, it is important for data analysis purposes to
estimate whether they will enter the band with non-negligible eccentricity ()
(Porter & Sesana 2010). The middle panel of Figure 4 displays the distribution of
eccentricities at ^{13}

## 5. Summary

With our results, we are moving closer towards a consistent solution to the Final Parsec Problem, and thus of providing a dynamical substantiation to the cosmological scenario where prompt coalescences are assumed during the course of galaxy evolution (Sesana et al. 2007; Volonteri 2010). Our results suggest that the formation of eccentric binaries, followed by a quick orbital decay, could result from the expected development of global non-axisymmetries in galaxies after they merge. Our gas-poor merger models show only rather mild departures from axisymmetry and a small amount of rotation; we believe that stronger departures from axisymmetry—to be expected from higher amount of rotation—, and the presence of gas will only reinforce our conclusions. It seems, therefore, probable that prompt coalescences result from mergers of irregular galaxies expected to be common at high redshift. Based on our prompt MBHB coalescence results, we expect that LISA will see events per year depending on the MBH seed model (Sesana et al. 2009; Volonteri 2010). Moreover, even though GWs circularizes the MBHBs during the late relativistic phase of inspiral, they are likely to have some residual () eccenticity when entering the LISA band and a broad distribution () in the PTA band.

### Footnotes

- affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, University of Heidelberg, D-69120 Heidelberg, Germany
- affiliation: Max Planck Institut für Gravitationsphysik (Albert-Einstein-Institut), D-14476 Potsdam, Germany
- affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, University of Heidelberg, D-69120 Heidelberg, Germany
- affiliation: Institut für Theoretische Astrophysik, Zentrum für Astronomie der University of Heidelberg, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany
- affiliation: National Astronomical Observatories of China, Chinese Academy of Sciences NAOC/CAS, 20 A Datun Rd., Chaoyang District, Beijing 100012, China
- affiliation: Main Astronomical Observatory (MAO), National Academy of Sciences of Ukraine (NASU), Akademika Zabolotnoho 27, 03680 Kyiv, Ukraine)
- affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, University of Heidelberg, D-69120 Heidelberg, Germany
- affiliation: National Astronomical Observatories of China, Chinese Academy of Sciences NAOC/CAS, 20 A Datun Rd., Chaoyang District, Beijing 100012, China
- affiliation: The Kavli Institute for Astronomy and Astrophysics at Peking University
- The loss cone is the region of phase space corresponding, roughly speaking, to orbits that cross the binary, i.e. with angular momentum , where (Lightman & Shapiro 1977).
- Models with and correspond to moderately oblate and prolate shapes, respectively.
- , where is the speed of light, is the Schwarzschild radius.

### References

- Aarseth, S. Gravitational N-Body Simulations (Cambridge, UK: Cambridge University Press, November 2003.)
- Armitage, P. & Natarajan, P. 2005, ApJ, 634, 921
- Begelman, M. C., Blandford, R. D. & Rees, M. J. 1980, Nature, 287, 307
- Berczik, P., Merritt, D., & Spurzem, R. 2005, ApJ, 633, 680
- Berczik, P., Merritt, D., Spurzem, R. & Bischoff, H. P. 2006, ApJ, 642, 21
- Berczik, P., Berentzen, I., Nitadori, K. Preto, M. & Spurzem, R. 2011, to be submitted to ApJ
- Berczik, P., Nitadori, K., Hamada, T. & Spurzem, R., 2011, in preparation
- Berentzen, I., Preto, M., Berczik, P., Merritt, D. & Spurzem, R. 2009, ApJ, 695, 455 (Paper I)
- Callegari, S., Kazantzidis, S., Mayer, L., Colpi, M., Bellovary, J. M., Quinn, T. & Wadsley, J., 2011, ApJ, 729, 85
- Colpi, M. & Dotti, M. 2009, arXiv:0906.4339 , Invited Review to appear on Advanced Science Letters (ASL), Special Issue on Computational Astrophysics, edited by Lucio Mayer
- Dehnen, W. 1993, MNRAS, 265, 250
- Ferrarese, L. & Ford, H. 2005, Space Sci. Rev., 116, 523
- Harfst, S. and Gualandris, A. and Merritt, D. and Spurzem, R. and Portegies Zwart, S. & Berczik, P. 2000, New Astronomy, 12, 357
- Lightman, A. P. & Shapiro, S. L. 1977, ApJ, 211, 244
- Kochfar, S. & Burkert, A. 2006, A&A, 445, 403
- Merritt, D. & Poon, M.Y. 2004, ApJ, 606, 788
- Merritt, D., Mikkola, S. & Szell, A. 2007, ApJ, 671, 53
- Milosavljević, M. & Merritt, D. 2003, ApJ, 596, 860
- Quinlan, G. D. 1996, New Astronomy, 1, 35
- Quinlan, G. D. & Hernquist, L. 1997, New Astronomy, 2, 533
- Perets, H. & Alexander, T. 2008, ApJ, 677, 146
- Peters, P.C. 1964, Phys Rev, 136, 1224
- Porter, E. & Sesana, A. 2010, arXiv:1005.5296
- Preto, M., Berentzen, I., Berczik, P., Merritt, D. & Spurzem, R. 2009, Journ of Phys Conf Ser, 154, 012049
- Preto, M. & Amaro-Seoane, P. 2010, ApJ, 708, 42
- Preto, M., Berentzen, I., Berczik, P. & Spurzem, R. 2011, to be submitted to ApJ
- Schödel, R., Merritt, D., & Eckart, A. 2009, A&A, 502, 91
- Sesana, A. 2010, ApJ, 719, 851
- Sesana, A.,Volonteri, M. & Haardt, F. 2010, MNRAS, 377, 1711
- Sesana, A.,Volonteri, M. & Haardt, F. 2009, Classical and Quantum Gravity, 26, 4033
- Spitzer, L. 1987, Dynamical evolution of globular clusters (Princeton, NJ, Princeton University Press, 1987, 191 p.)
- Tremaine, S., Richstone, D. O., Byun, Y.-I., Dressler, A., Faber, S. M., Grillmair, C., Kormendy, J., & Lauer, T. R. 1994, AJ, 107, 634
- Volonteri, M. 2010, A&A Rev., 18, 279
- Yu, Q. 2002, MNRAS, 331, 935