1 Introduction

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


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)


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 being a slow diffusive process, it is only in low-luminosity galaxies harboring MBHs of mass that central relaxation times are short enough to drive the binary to coalescence in less than a Hubble time (Merritt et al. 2007).

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
Table 1-body integrations. column: mass of the MBH binary; Other columns: particle number ; First two lines: simulations of spherical nuclei; Last four lines: Simulations of merging nuclei; measures the initial orbital angular momentum of the MBH binary ( for spherical nuclei), or otherwise it measures the initial orbital angular momentum of the merging nuclei ( for near-radial merger and for less radial merger). All nuclei have ; all binaries have equal mass .

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

Figure 1.— Binary hardening. Upper panel: in a spherical nucleus, decreases with . Middle panel: in a merging nucleus, is N-independent. Lower panel: hardening rates as a function of for different . Being much smaller, of the binary has been multiplied by to better fit in the plot. Labels ’s’ for spherical and ’m’ for merger.

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)


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


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


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)


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

Figure 2.— Triaxiality and mass flattening of merging nuclei. Shown are merging binaries of total mass with (upper panel), and with (almost radial mergers, in the middle and lower panels). and are measured in five mass shells between and , each of width . Triaxiality decreases over time, the faster the heavier the binary is. Mass flattening is constant.

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.

Figure 3.— Long term eccentricity evolution. Red and green lines represent NB and semi-analytic evolution without radiation reaction. Blue and magenta lines correspond to semi-analytic solution, including radiation reaction, for and , respectively.

Following Merritt & Poon (2004), we measure the triaxiality of the nucleus with . 12 Figure 2 depicts the evolution of and of the flattening for several mass shells of merging nuclei. The value of of each remnant, immediately after the merger, depends on the initial . In the case of a near radial merger, , the remnant is prolate and evolves over time towards an oblate spheroidal shape; for the remnant is an oblate spheroid from the very beginning. The triaxiality decreases over time, and the rate at which it changes is faster the larger is. The triaxiality remains significant, in the inner mass shells, until the binary reaches the relativistic phase in all models with the smallest (and more realistic) values of , and also in most of the other cases. The flattening is constant throughout in all cases, so the asymptotic shape of the merger is that of an oblate spheroid. We conclude that the rather mild triaxiality created during the merger supports a family of centrophilic orbits that keep the loss cone full () at all times until the binary reaches relativistic separations .

4. Eccentricity evolution and time scales for coalescence

Figure 4.— Upper panel: Range of coalescence time for binaries with and . Middle panel: distribution of eccentricities, for , at . Lower panel: distribution of eccentricities, for , when for PTAs.

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


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


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—most binaries will not be fully circularized by then. We expect therefore that the eccentricity in the LISA band will be non-negligible. Finally, the lower panel of Figure 4 depicts the eccentricity distribution at for the PTA band. We see that eccentricities are quite high—peaking at . The results presented here concerning the coalescence times and eccentricity growth corroborate recent three-body scattering studies—which had to treat the average hardening rate as a free parameter (Sesana 2010).

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.

MP acknowledges support by DLR (Deutsches Zentrum für Luft- und Raumfahrt). We acknowledge support by the Chinese Academy of Sciences Visiting Professorship for Senior International Scientists, Grant Number 2009S1-5 (The Silk Road Project) (RS and PB). The special supercomputer Laohu at the High Performance Computing Center at National Astronomical Observatories, funded by Ministry of Finance under the grant ZDYZ2008-2, has been used. Simulations were also performed on the GRACE supercomputer (grants I/80 041-043 and I/84 678-680 of the Volkswagen Foundation and 823.219-439/30 and /36 of the Ministry of Science, Research and the Arts of Baden-Württemberg). We thank the DEISA Consortium (http://www.deisa.eu), cofunded through EU FP6 projects RI-508830 and RI-031513, for support within the DEISA Extreme Computing Initiative. The Kolob cluster is funded by the excellence funds of the University of Heidelberg in the Frontier scheme.


  1. affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, University of Heidelberg, D-69120 Heidelberg, Germany
  2. affiliation: Max Planck Institut für Gravitationsphysik (Albert-Einstein-Institut), D-14476 Potsdam, Germany
  3. affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, University of Heidelberg, D-69120 Heidelberg, Germany
  4. affiliation: Institut für Theoretische Astrophysik, Zentrum für Astronomie der University of Heidelberg, Albert-Ueberle-Str. 2, D-69120 Heidelberg, Germany
  5. affiliation: National Astronomical Observatories of China, Chinese Academy of Sciences NAOC/CAS, 20 A Datun Rd., Chaoyang District, Beijing 100012, China
  6. affiliation: Main Astronomical Observatory (MAO), National Academy of Sciences of Ukraine (NASU), Akademika Zabolotnoho 27, 03680 Kyiv, Ukraine)
  7. affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, University of Heidelberg, D-69120 Heidelberg, Germany
  8. affiliation: National Astronomical Observatories of China, Chinese Academy of Sciences NAOC/CAS, 20 A Datun Rd., Chaoyang District, Beijing 100012, China
  9. affiliation: Astronomisches Rechen-Institut, Zentrum für Astronomie, University of Heidelberg, D-69120 Heidelberg, Germany
  10. affiliation: The Kavli Institute for Astronomy and Astrophysics at Peking University
  11. 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).
  12. Models with and correspond to moderately oblate and prolate shapes, respectively.
  13. , where is the speed of light, is the Schwarzschild radius.


  1. Aarseth, S. Gravitational N-Body Simulations (Cambridge, UK: Cambridge University Press, November 2003.)
  2. Armitage, P. & Natarajan, P. 2005, ApJ, 634, 921
  3. Begelman, M. C., Blandford, R. D. & Rees, M. J. 1980, Nature, 287, 307
  4. Berczik, P., Merritt, D., & Spurzem, R. 2005, ApJ, 633, 680
  5. Berczik, P., Merritt, D., Spurzem, R. & Bischoff, H. P. 2006, ApJ, 642, 21
  6. Berczik, P., Berentzen, I., Nitadori, K. Preto, M. & Spurzem, R. 2011, to be submitted to ApJ
  7. Berczik, P., Nitadori, K., Hamada, T. & Spurzem, R., 2011, in preparation
  8. Berentzen, I., Preto, M., Berczik, P., Merritt, D. & Spurzem, R. 2009, ApJ, 695, 455 (Paper I)
  9. Callegari, S., Kazantzidis, S., Mayer, L., Colpi, M., Bellovary, J. M., Quinn, T. & Wadsley, J., 2011, ApJ, 729, 85
  10. 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
  11. Dehnen, W. 1993, MNRAS, 265, 250
  12. Ferrarese, L. & Ford, H. 2005, Space Sci. Rev., 116, 523
  13. Harfst, S. and Gualandris, A. and Merritt, D. and Spurzem, R. and Portegies Zwart, S. & Berczik, P. 2000, New Astronomy, 12, 357
  14. Lightman, A. P. & Shapiro, S. L. 1977, ApJ, 211, 244
  15. Kochfar, S. & Burkert, A. 2006, A&A, 445, 403
  16. Merritt, D. & Poon, M.Y. 2004, ApJ, 606, 788
  17. Merritt, D., Mikkola, S. & Szell, A. 2007, ApJ, 671, 53
  18. Milosavljević, M. & Merritt, D. 2003, ApJ, 596, 860
  19. Quinlan, G. D. 1996, New Astronomy, 1, 35
  20. Quinlan, G. D. & Hernquist, L. 1997, New Astronomy, 2, 533
  21. Perets, H. & Alexander, T. 2008, ApJ, 677, 146
  22. Peters, P.C. 1964, Phys Rev, 136, 1224
  23. Porter, E. & Sesana, A. 2010, arXiv:1005.5296
  24. Preto, M., Berentzen, I., Berczik, P., Merritt, D. & Spurzem, R. 2009, Journ of Phys Conf Ser, 154, 012049
  25. Preto, M. & Amaro-Seoane, P. 2010, ApJ, 708, 42
  26. Preto, M., Berentzen, I., Berczik, P. & Spurzem, R. 2011, to be submitted to ApJ
  27. Schödel, R., Merritt, D., & Eckart, A. 2009, A&A, 502, 91
  28. Sesana, A. 2010, ApJ, 719, 851
  29. Sesana, A.,Volonteri, M. & Haardt, F. 2010, MNRAS, 377, 1711
  30. Sesana, A.,Volonteri, M. & Haardt, F. 2009, Classical and Quantum Gravity, 26, 4033
  31. Spitzer, L. 1987, Dynamical evolution of globular clusters (Princeton, NJ, Princeton University Press, 1987, 191 p.)
  32. Tremaine, S., Richstone, D. O., Byun, Y.-I., Dressler, A., Faber, S. M., Grillmair, C., Kormendy, J., & Lauer, T. R. 1994, AJ, 107, 634
  33. Volonteri, M. 2010, A&A Rev., 18, 279
  34. Yu, Q. 2002, MNRAS, 331, 935
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description