Formation of galactic nuclei with multiple supermassive black holes at high redshifts
aa r X i v : . [ a s t r o - ph . C O ] J u l Mon. Not. R. Astron. Soc. , 000–000 (2011) Printed 14 November 2018 (MN L A TEX style file v2.2)
Formation of galactic nuclei with multiple supermassive black holesat high redshifts
Girish Kulkarni , ⋆ , Abraham Loeb † Institute for Theory & Computation, Harvard University, 60 Garden Street, Cambridge, MA 02138 USA Harish-Chandra Research Institute, Chhatnag Road, Jhunsi, Allahabad 211019 India
ABSTRACT
We examine the formation of groups of multiple supermassive black holes (SMBHs) in gas-poor galactic nuclei due to the high merger rate of galaxies at high redshifts. We calculatethe relative likelihood of binary, triple, and quadruple SMBH systems, by considering thetimescales for relevant processes and combining merger trees with N-body simulations for thedynamics of stars and SMBHs in galactic nuclei. Typical haloes today with mass M ≈ M ⊙ have an average mass M z = = × M ⊙ at z ∼
6, while rare haloes with current mass M & M ⊙ have an average mass M z = = × M ⊙ at that redshift. These cluster-sizehaloes are expected to host single galaxies at z ∼
6. We expect about 30% galaxies withinhaloes with present-day mass M ≈ M ⊙ to contain more than two SMBHs at redshifts2 . z .
6. For larger present-day haloes, with M & M ⊙ , this fraction is almost 60%. Theexistence of multiple SMBHs at high redshifts can potentially explain the mass deficienciesobserved in the cores of massive elliptical galaxies, which are up to 5 times the mass of theircentral BHs. Multiple SMBHs would also lead to an enhanced rate of tidal disruption of stars,modified gravitational wave signals compared to isolated BH binaries, and slingshot ejectionof SMBHs from galaxies at high speeds in excess of 2000 km s − . Key words: black hole physics – galaxies: nuclei – galaxies:evolution – galaxies:high redshift– galaxies:kinematics and dynamics
Most local galaxies host supermassive black holes (SMBHs)at their centres (Richstone et al. 1998; Ferrarese & Ford 2005).The SMBH mass M bh is correlated with properties of thespheroidal nucleus of the host galaxy, such as velocity disper-sion (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Ferrarese2002; G¨ultekin et al. 2009) and luminosity (Magorrian et al. 1998;McLure & Dunlop 2002; Marconi & Hunt 2003; G¨ultekin et al.2009). Detection of bright quasars at redshifts z & ∼ × M ⊙ already existed at z ∼
7. In the stan-dard Λ CDM cosmological model, growth of galaxies is hierar-chical and galaxy mergers are expected to be particularly fre-quent at redshifts z ∼ ff mann & Haehnelt 2000;Menou et al. 2001; Bromm & Loeb 2003; Volonteri et al. 2003;Hopkins et al. 2006; Tanaka & Haiman 2009).Existing merger tree models are based on the assumption ⋆ Email: [email protected] † Email: [email protected] that any binary black hole system, which inevitably forms in agalaxy’s merger history, coalesces on a short time-scale. How-ever, the evolution of SMBH binaries is a complex open prob-lem and it is unclear if a binary can merge within a Hubble time(Merritt & Milosavljevi´c 2005). One expects that during a mergerevent of two galaxies, the dynamics of their constituents SMBHswould proceed in three stages (Begelman et al. 1980). In the firststage, the SMBHs sink to the centre of the gravitational potentialof the merger remnant by dynamical friction and form a gravita-tionally bound binary. The newly-formed binary continues to loseenergy and angular momentum through its global gravitational in-teraction with many stars until the separation between the SMBHsreduces to a value at which the dominant mechanism of energy lossis the 3-body interaction between the binary and individual stars.This is the second stage of the binary’s evolution, and is knownas the ‘hard stage.’ The precise definition of a hard SMBH binaryvaries in the literature, but it is commonly assumed that the binarybecomes hard when its semi-major axis a reaches a value given by(Yu 2002) a ≈ a h ≡ Gm σ = . m M ⊙ !
200 km s − σ ! pc , (1)where stars in the galactic nucleus are assumed to have a one-dimensional velocity dispersion σ , and m denotes the mass of thelighter SMBH. Finally, once the SMBH separation decreases to a c (cid:13) Kulkarni & Loeb small-enough value, gravitational wave emission becomes the dom-inant mode of energy loss and the SMBHs coalesce rapidly. This isthe third stage of the SMBH binary evolution. The value of semi-major axis a at which the coalescence time scale due to gravita-tional wave emission alone is t is given by (Peters 1964; Loeb 2010) a ( t ) ≡ a gw ( t ) = . × − t yr ! / M × M ⊙ ! / pc , (2)where M is the total mass of the binary, and we have considered twoSMBHs with mass 10 M ⊙ each on a circular orbit (with shortertime scale at increasing eccentricity). Gravitational wave emissiontakes over as the dominant mode of energy loss when a = a gw ( t h ),where t h is the hardening time scale.Among these three stages of evolution of an SMBH binary,the largest uncertainty in the binary’s lifetime originates from thehard stage, which can be the slowest of the three stages sincethe binary quickly ejects all low angular momentum stars in itsvicinity, thus cutting o ff its supply of stars. This is known asthe “final parsec problem” (Milosavljevi´c & Merritt 2003b). Forexample, Yu (2002) studied coalescence of SMBH binaries ina sample of galaxies observed by Faber et al. (1997) and foundthat spherical, axisymmetric or weakly triaxial galaxies can allhave long-lived binary SMBHs that fail to coalesce. Similarly,Merritt & Milosavljevi´c (2005) found that the time spent by a bi-nary is less than 10 yr only for binaries with very low mass ratios( . − ). Furthermore, Merritt & Milosavljevi´c (2005) showedthat a binary may not be able to interact with all the stars in its losscone, thereby increasing the time spent in the hard stage even fur-ther; they found that in a nucleus with a singular isothermal spherestellar density profile, an equal-mass binary will stall at a separationof a ≈ a h / .
5, where we have defined a h in Equation (1). The finalseparation is expected to be even higher for galaxies with shallowerdensity profiles.Several ways have been discussed in the literature to e ffi -ciently extract energy and angular momentum from a hard SMBHbinary and overcome the final parsec problem. An example is thework by Armitage & Natarajan (2002), who suggested that gas cancatalyse the coalescence of a hard SMBH binary by serving asan e ff ective sink for the binary’s angular momentum. In particu-lar, they found that a binary with a separation of 0 . years with-out significant enhancement in the gas accretion rate. Similarly,Escala et al. (2004, 2005) found that in SPH simulations, cloudsof hot gas ( T gas ≈ T virial ) can induce decay of orbits of embeddedbinary point masses due to gravitational drag. A caveat to thesestudies is that feedback from gas accretion onto the SMBHs can re-move the rest of the gas from the merger remnant before the binarycoalesces. However, stellar dynamical processes could also acceler-ate binary coalescence, without gas. For example, Merritt & Poon(2004) considered the e ff ect of chaotic orbits in steep triaxial po-tentials. They found that stars are supplied to the central black holeat a rate proportional to the fifth power of the stellar velocity dis-persion and that the decay rate of a central black hole binary wouldbe enhanced even if only a few percent of the stars are on chaoticorbits, thus solving the final parsec problem. As another example,it was suggested that a third SMBH closely interacting with a hard However, for such low mass ratios the time taken by the lighter black holeto reach the galactic nucleus due to dynamical friction is itself expected toexceed the Hubble time.
SMBH binary can reduce the binary separation to a small value ei-ther due to the eccentricity oscillations induced in the binary via theKozai-Lidov mechanism (Blaes et al. 2002) or due to repopulationof the binary’s loss cone due to the perturbation in the large-scalepotential caused by the third black hole (Ho ff man & Loeb 2007).Blaes et al. (2002) found that the merger time scale of an inner cir-cular binary can be shortened by as much as an order of magnitude,and that general relativistic precession does not destroy the Kozai-Lidov e ff ect for hierarchical triples that are compact enough.In summary, there is substantial uncertainty in the currentunderstanding of the evolution of binary SMBHs. Clearly, if theSMBH binary coalescence time is longer than the typical time be-tween successive major mergers of the galaxy, then more than twoSMBHs may exist in the nucleus of a merger remnant. We studythis possibility in this paper. We calculate the relative likelihoodof binary, triple, and quadruple SMBH systems, by consideringthe timescales for relevant processes and combining galaxy mergertrees with direct-summation N-body simulations for the dynamicsof stars and SMBHs in galactic nuclei. An obvious question regard-ing galactic nuclei with multiple SMBHs is whether such systemscan be long-lived. We consider this question here. Finally, systemswith multiple SMBHs are likely to be interesting because of obser-vational e ff ects involving their e ff ect on the properties of the hostbulge, the enhancement in the rate of tidal disruption of stars, theirassociated gravitational wave and electromagnetic signals, and theslingshot ejection of SMBHs at high speeds. We study some ofthese e ff ects.In § § §
5, with their re-sults shown in §
6. We consider the observational signatures of ourfindings in §
7. Finally, we discuss and summarise our primary find-ings in § Galactic nuclei with multiple SMBHs were first studied bySaslaw et al. (1974), who computed orbits of three and four SMBHsystems by sampling the parameter space of the problem. Theyshowed that if an infalling SMBH is lighter than the components ofthe pre-existing binary, then the most probable outcome is a sling-shot ejection in which the infalling SMBH escapes at a velocitythat is about a third of the orbital velocity of the binary. Valtonen(1976) further showed that the ejection velocity can be significantlyenhanced if drag forces due to gravitational radiation are accountedfor in the three-body dynamics. The formation of systems withthree or four SMBHs in a hierarchical merger of smooth galac-tic potentials was first studied by Mikkola & Valtonen (1990) andValtonen et al. (1994) with the objective of understanding the struc-ture of extragalactic radio sources. This line of work was extendedto binary-binary scattering of SMBHs by Hein¨am¨aki (2001), andby Ho ff man & Loeb (2007), who studied repeated triple interac-tions in galactic nuclei. Both of these studies used cosmologicallyconsistent initial conditions based on the extended Press-Schechtertheory. Systems with a larger number of black holes were stud-ied by Hut & Rees (1992) and Xu & Ostriker (1994) using simpleanalytical models and numerical calculations of massive particlesin smooth galactic potentials. Xu & Ostriker (1994) concluded thatthe most-likely outcome in these cases is one in which most blackholes are ejected and the galactic center is left with zero, or one, c (cid:13) , 000–000 ultiple Black Holes or two black holes. Finally, full N-body simulations of galacticnuclei with constituent SMBHs were performed for the case oftwo successive mergers by Makino & Ebisuzaki (1996), Makino(1997), and Iwasawa et al. (2006). Much of this work on SMBHswas based on earlier studies of stellar-mass black holes in glob-ular clusters. Sigurdsson & Hernquist (1993) and Kulkarni et al.(1993) considered the evolution of ∼
100 stellar mass black holesin globular clusters. They concluded that after mass segregation,most of these black holes are ejected out on a short time scale,and the globular cluster is left with none or a few black holes.Mass segregation and associated e ff ects of stellar-mass black holesin a galactic nucleus with a central SMBH was also considered(Miralda-Escud´e & Gould 2000; Freitag et al. 2006).The possible formation of systems with multiple SMBHs dueto successive galactic mergers arises naturally in any model de-scribing the hierarchical assembly of galaxies. One approach tomodeling SMBH growth involves constructing semi-analytic pre-scriptions of various characteristic processes, like mergers of galax-ies, formation of spheroids, star formation, and gas thermodynam-ics, coupled with merger trees of dark matter haloes. This ap-proach has been adopted, for example, by Kau ff mann & Haehnelt(2000), who also extended it to study possible formation ofmultiple SMBH systems and implications for the M bh – σ rela-tion and density profiles observed in luminous elliptical galaxies(Haehnelt & Kau ff mann 2002). Another study by Volonteri et al.(2003) followed merger trees of dark matter haloes and their com-ponent SMBHs using Monte Carlo realizations of hierarchicalstructure formation in the Λ CDM cosmology. They modeled darkmatter haloes as singular isothermal spheres and calculated the in-spiral of less massive halos in more massive ones by using theChandrasekhar formula for dynamical friction. Gas accretion tothe SMBHs was modeled so as to reproduce the empirical M bh – σ relation and the SMBH dynamics was described with analyticprescriptions. In particular, the coalescence time of hard SMBHbinaries was calculated from a set of coupled di ff erential equa-tions based on scattering experiments involving the ejection of stel-lar mass from the loss cone due to the hard SMBH binary andthe resultant change in the hardening rate (Quinlan 1996; Merritt2000). For galaxies that underwent another major merger beforetheir constituent binary SMBH coalesced, a three-body interac-tion was implemented between the binary and the intruder SMBH.They found that the smallest SMBH was kicked out of the galaxyin 99% of cases, while the binary escapes the galaxy in 8 % ofcases. Thus, a significant fraction of galactic nuclei could endup with no SMBHs or o ff set SMBHs with mass lower than thatexpected from the M bh – σ relation. These results were later ex-tended to incorporate recoil in the SMBH merger remnant due toasymmetric emission of gravitational waves, which mainly a ff ectedthe M bh – σ relation for low-mass haloes by increasing the scatter(Volonteri & Rees 2006; Volonteri 2007; Blecha et al. 2011). Sim-ilar semi-analytic models were studied by several other authorsto understand the assembly of z ∼ Unless they coalesce rapidly, or get kicked out of the host galacticnucleus, we expect multi-SMBH systems to form in galactic nucleiat high redshift due to mergers of galaxies if the typical black holecoalescence timescale is longer than the feeding timescale of newincoming black holes. In this section, we establish a simple theoret-ical framework for this formation path using analytical estimates ofits relevant timescales: (i) the major merger time scale of galaxies; (ii) the time scale on which a satellite galaxy sinks to the center of ahost galaxy so that a close interaction between SMBHs can occur;and (iii) the time scale of SMBH coalescence.
Fakhouri et al. (2010) have quantified the average merger rate ofdark matter haloes per halo per unit redshift per unit mass ratio fora wide range of halo mass, progenitor mass ratios and redshift. Theresult is given by a fitting formula derived from the Millennium(Springel et al. 2005) and Millennium-II (Boylan-Kolchin et al.2009) simulations: dNd ξ dz ( M , ξ, z ) = A M M ⊙ ! α ξ β exp " ξ ˜ ξ ! γ (1 + z ) η . (3)Here, M is the halo mass at redshift z , and ξ is the mass ratio of pro-genitors. Mergers with ξ > . α = . β = − . γ = . η = . A = . ξ = . × − . Theaverage major merger rate per unit time is then given by dN m dt ( M , z ) = Z . d ξ dNd ξ dz ( M , ξ, z ) dzdt . (4)Fakhouri et al. (2010) also provide a fitting formula for averagemass growth rate of halos that can be used to calculate the halomass at redshift z for use in equation (3),˙ M ( z ) = . ⊙ yr (1 + . z ) p Ω m (1 + z ) + Ω Λ M M ⊙ ! . . (5)Using equation (4) we can now define the time scale of major merg-ers for a halo as t mrg = " dN m dt − . (6)The behavior of this quantity is shown in Figure 1 for three halomasses that discussed here: a Milky Way-like halo that has a mass M = M ⊙ at z =
0, the typical halo today that has mass M = M ⊙ at z =
0, and rare haloes with mass M = M ⊙ at z =
0. (In this paper, M always denotes the halo mass at redshift z =
0. We also refer to the average mass of such haloes at otherredshifts, by e. g. M z = and M z = . A halo with M = M ⊙ will c (cid:13) , 000–000 Kulkarni & Loeb
Figure 1.
Halo major merger time scale (mass ratio > . M = M ⊙ (blue solid line), 10 M ⊙ (blue dashed line) and 10 M ⊙ (blue dot-dashed line). The Hubbletime is shown by the solid red curve. Major mergers are more frequent athigher redshifts. On average, Milky Way-sized haloes are not expected toundergo a major merger for z .
1. Galaxy major merger time scale is alwayslonger than the subsequent dynamical friction time scale. have M z = = × M ⊙ . A halo with M = M ⊙ will have M z = = × M ⊙ . A cluster-size halo, with mass M = M ⊙ will have M z = = × M ⊙ and is expected to hold a singlegalaxy at that redshift.) This is the time scale at which we expectnew (satellite) haloes to enter the halo. As expected, halo mergersare more frequent at higher redshift. At redshift z . ff erent that the time scale for major mergers of dark matterhaloes calculated in Equation (6).The dynamical friction time scale is often estimated us-ing Chandrasekhar’s formula (Chandrasekhar 1943; Lacey & Cole1993; Binney & Tremaine 2008): t df = f df Θ orb ln Λ M host M sat t dyn , (7)where M host and M sat are the masses for the host and satellitehaloes respectively, ln Λ is the coulomb logarithm, Θ orb is a func-tion of the orbital energy and angular momentum of the satellite, f df is an adjustable parameter of order unity and t dyn is the halodynamical time scale calculated at the virial radius. Equation (7)is valid only in the limit of small satellite mass in an infinite,isotropic and homogeneous collisionless medium. Still, it has beenused in the literature even for large satellite masses by modifyingthe Coulomb logarithm. In recent years, deviations from predic-tions by equation (7) have been reported in both the M sat ≪ M host and M sat . M host regimes (Ta ff oni et al. 2003; Monaco et al. 2007;Boylan-Kolchin et al. 2008; Jiang et al. 2008; Wetzel et al. 2009).To correct the problems associated with Chandrasekhar’s for-mula, several groups have developed full dynamical models ofevolution of merging haloes (Taylor & Babul 2001; Gnedin 2003; Ta ff oni et al. 2003; Zentner et al. 2005). For example, one of theapproaches to overcome the limits of Chandrasekhar’s formula isthe theory of linear response (TLR; Colpi et al. 1999). TLR cap-tures the backreaction of the stellar distribution to the intrudingsatellite by correlating the instantaneous drag force on it with thedrag force at an earlier time via the fluctuation-dissipation theorem.Tidal stripping of a satellite halo is an important ingredient in thisformulation. In a singular isothermal sphere with 1D velocity dis-persion σ and density profile ρ ( r ) = σ / [2 π Gr ], TLR predicts asinking time t df = . r V cir GM sat ln Λ ǫ α , (8)where ǫ is the circularity (defined as the ratio between the angularmomentum of the current orbit relative to that of a circular orbit ofequal energy), r cir and V cir are the initial radius and velocity of thecircular orbit with the same energy of the actual orbit, and M S is themass of the incoming satellite halo. Numerical simulations suggesta value of 0 . − . α (van den Bosch et al. 1999;Colpi et al. 1999; Volonteri et al. 2003).Given the limitations of analytical treatments, we turn to re-sults of numerical simulations to understand the dynamical fric-tion time scale. Using N-body simulations, Boylan-Kolchin et al.(2008) give a fitting formula that accurately predicts the time-scalefor an extended satellite to sink from the virial radius of a host halodown to the halo’s centre for a wide range of mass ratios and orbits(including a central bulge in each galaxy changes the merging timescale by .
10 %). Their fitting formula is given by t df t dyn = A ξ − b ln(1 + /ξ ) exp " c jj cir ( E ) r cir ( E ) r vir d , (9)where A = . b = . c = . d = .
0. Here ξ is the massratio M sat / M host , j is the specific angular momentum of the satellitehalo, and j cir is the specific angular momentum of a circular orbitwith the same energy E . This formula is expected to be valid for0 . ξ .
0, and for circularities 0 . η ≡ j / j cir ( E ) . η ≈ . − . r cir ( E ) / r vir .
0. This covers the peak value of distribution seen in cosmologicalN-body simulations. We fix r cir ( E ) / r vir = . η = .
5, whichare the typical values found in simulations.We can now obtain the instantaneous merger rate of galax-ies by combining the halo merger rate and dynamical friction timescale. We closely follow the method of Shen (2009) and write B gal ( M , ξ, z ) = B [ M , ξ, z e ( z , ξ )] dz e dz , (10)where B ( M , ξ, z ) (per unit volume per unit mass per unit redshiftper unit mass ratio) is the instantaneous merger rate of halos withmass M , progenitors with mass ratio ξ at redshift z , B gal is the samequantity for galaxies. The redshift z e ( z , ξ ) is a function of z and ξ ,and is given implicitly by t ( z ) − t ( z e ) = t mrg ( ξ, z e ) , (11)where t ( z ) is the cosmic time at redshift z . Shen (2009) finds that dz e / dz is almost constant at all redshifts for ξ = . − dz e dz ≈ + . ξ . ln(1 + /ξ )] − , (12)for the fitting formula in equation (9). We assume this form in our c (cid:13) , 000–000 ultiple Black Holes Figure 2.
A comparison between the feeding time scale of incoming blackholes t in (black solid line; Eq. 13) and the time scale of black hole coales-cence t coal (black dot-dashed line; Eq. 23), for a halo mass M = M ⊙ and considering only mergers with a mass ratio ξ = .
4. The coalescencetime t coal has only a weak dependence on redshift because its dependenceon M bh and σ cancel out due to the M bh – σ relation. This figure shows thatat high redshift new black holes would arrive to the center of a galaxy fasterthan they could merge via dynamical processes. calculations. Once we have calculated B gal ( M , ξ, z ), we normalize itby n ( M , z ), the abundance of haloes of mass M at redshift z . We usethe Sheth-Tormen mass function (Sheth & Tormen 1999) to calcu-late n ( M , z ). This gives us the galaxy merger rate per halo per unit ξ per unit redshift, which is the galaxy’s counterpart of equation (3),and which we denote by dN gal / dz . The rate of mergers of galaxiesis the rate at which new black holes are added to the host halo’snucleus. Thus, the time scale of incoming black holes is t in = " dN gal dz dzdt − . (13)The result is shown by the solid black line in Figure 2 for a massratio of ξ = . M ⊙ at z = In order to find whether there is a generic possibility of formationof systems with multiple SMBHs, we compare the time scale onwhich new black holes are added to the galactic nucleus at a certainredshift with the coalescence time scale of a binary SMBH at thatredshift.As described in §
1, the formation and coalescence of a blackhole binary is expected to take place in three stages. We define thecoalescence time as the time that the binary spends in the secondof these stages, that is the time from when the binary separationis a = a h , defined in equation (1), up to when the separation is a = a gr at which point the binary enters the third stage of evolution,and gravitational waves become the dominant mechanism of en-ergy loss. For a hard binary, the dominant channel through whichenergy is lost is three-body interactions in which stars passing inclose proximity to the binary are ejected at a much higher velocity v ej = [ GM tot / a ] / , where M tot is the total mass of the binary. Thehardening time scale was quantified for a fixed stellar distribution by Quinlan (1996), who found a time scale of t h ( a ) ≡ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) a ˙ a (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = σ G ρ aH , (14)where a is the binary separation, ρ and σ are the density and one-dimensional velocity dispersion of the stellar background, and H is a dimensional parameter whose value was found from scatter-ing experiments to be 16 for a hard, equal-mass binary. In practice,however, the above expression for t h is valid only during the initialstages of the binary’s evolution. As the binary shrinks further, itejects stellar mass from the central regions and modifies the stellardensity ρ that appears in equation (14). This feedback can be quan-tified using a simple analytical model given by Merritt (2000), inwhich the binary evolution is described by two coupled equations,the first describing the binary’s hardening due to the presence ofstars, ddt a ! = H G ρσ , (15)and the second describing the change in stellar density due to ejec-tion of mass by the hard SMBH binary, dM ej d ln(1 / a ) = JM tot , (16)where M ej is the ejected mass, and J is another dimensionless pa-rameter that was measured by Quinlan (1996) to be close to unityand nearly independent of a .By assuming a singular isothermal sphere profile for the stel-lar density and assuming that the ejected stellar mass causes aconstant-density core to form at the center of this profile, Merritt(2000) finds that evolution of the binary separation can be describedas t − t init t = a h a " ln (cid:18) a h a (cid:19) − (cid:18) a h a (cid:19) + − aa h ! , (17)where a h is as defined in Equation (1), a ( t init ) = a h , and t is givenby t = π J H M tot m ! (cid:18) GM tot σ (cid:19) . (18)This result is found to closely match with the evolution observed inN-body simulation.On the other hand, the timescale for emission of gravitationalwaves is given by t gr = c a G m m M tot . (19)As a result, the binary will continue to harden only up to the timewhen hardening time t h = t gr , after which it will coalesce rapidlydue to gravitational wave emission. Using equation (17), it can beshown that this occurs when a = a gr where (Merritt 2000), a gr a h ≈ A | ln A | . , (20)and A = . m m ! . M tot m ! . (cid:18) σ c (cid:19) . (21)Here m and m are masses of the components of the SMBH binary.Finally, we can again use equation (17) to calculate the time it takesfor the binary to shrink from a = a h to a = a gr (Merritt 2000): t coal ≈ t A − | ln A | / , (22) c (cid:13) , 000–000 Kulkarni & Loeb which can be simplified as t coal ≈ . × yr m m ! . M tot m ! . M tot M ⊙ ! σ / s ! − . (23)Clearly, there is a possibility for the formation of multiple-SMBH system if t in < t coal . These two time scales are comparedin Figure 2 for a halo that has a mass of M = M ⊙ at z = ξ = .
4. At each redshift, we calculate t in from equation (13).In order to estimate t coal at a given redshift using equation (23),we first infer the mass of the halo at that redshift from the fittingfunction for the halo’s assembly history from equation (5). We thenassume that a galaxy belonging to a satellite halo with mass ratio ξ has merged with this host halo at this redshift.In order to estimate the mass of black holes in the nuclei ofthese galaxies, we follow the approach of Ho ff man & Loeb (2007)in employing the M bh – σ relation. The virial velocity (defined as thecircular velocity at virial radius) for a halo of mass M at redshift z is given by v vir = . M h − M ⊙ ! / " Ω m Ω zm ∆ c π / + z ! / km / s , (24)where Ω zm = Ω m (1 + z ) Ω m (1 + z ) + Ω Λ + Ω k (1 + z ) , (25)and ∆ c is the overdensity of the halo relative to the critical density,given for the Λ CDM cosmology by ∆ c = π + d − d , (26)where d = Ω zm − v c of its constituentspheroid and obtain the velocity dispersion of the spheroid usingthe relation (Ferrarese 2002) v c ≈ " σ / s . km / s . (27)This combined with the M bh – σ relation (Tremaine et al. 2002) σ / s ≈ M bh . × M ⊙ / . , (28)gives M halo M ⊙ ! = . M bh M ⊙ ! " Ω m Ω zm ∆ c π − / (1 + z ) − / . (29)We obtain the black hole masses in the host and the satellitehaloes using equation (29) and use the spheroid velocity dispersionfrom equation (27) to estimate the coalescence time from equation(23). The result is shown by the dashed line in Figure 2.At high redshift, early on in the assembly history of a halo,the galaxy merger rate is higher than the SMBH binary coalescencerate and systems with multiple SMBHs can form. Note that the timescale t coal obtained above will change if e ff ect of loss-cone replen-ishment and gas are taken into account. However, Yu (2002) findsthat in realistic spheroidal galaxies, even loss-cone replenishmentis insu ffi cient to cause early coalescence. We have described the literature on systems with more than twoSMBHs in §
2. If the infalling SMBH is less masssive than either
Figure 3.
An example merger tree form the Millennium simulation of ahalo that has a mass of M ∼ M ⊙ . This plot shows major mergers(mass ratio > . in all branches of the halo’s merger tree. of the components of a pre-existing binary then we expect the ulti-mate outcome to be ejection of the smaller SMBH and recoil of thebinary. Ho ff man & Loeb (2007) studied the statistics of close tripleSMBH encounters in galactic nuclei by computing a series of three-body orbits with physically motivated initial conditions appropriatefor giant elliptical galaxies. Their simulations included a smoothbackground potential consisting of a stellar bulge and a dark matterhalo, and also accounted for the e ff ect of dynamical friction due tostars and dark matter. They found that in most cases the intruderhelped the binary SMBH to coalesce via the Kozai-Lidov mecha-nism and by scattering stars into the binary’s loss cone. In this case,the intruder itself was left wandering in the galactic halo, or evenkicked out of the galaxy altogether. It was also found that escape ofall three black holes is exceedingly rare.Dynamical evolution of multiple massive black holes in glob-ular clusters has received much attention (Kulkarni et al. 1993;Sigurdsson & Hernquist 1993). From these studies, it is expectedthat systems with more than two SMBHs will last for about a cross-ing time. In order to accurately calculate the formation and evolution ofgalactic nuclei with multiple black holes, we perform direct-summation N-body simulations of galactic nuclei merging in a cos-mological context. This essentially involves generating physicallyconsistent initial conditions for galactic nuclei with SMBHs at highredshift and evolving them while taking into account the mergersof such nuclei and the resultant close interaction of their SMBHs.We obtain merger histories of galactic nuclei by extractingmerger trees of gravitationally bound subhaloes from the Millen-nium Simulation Database , which stores results of the MillenniumSimulation (Springel et al. 2005). The Millennium Simulation is apure dark matter simulation with a Λ CDM model with 2160 par-ticles in a periodic cube 500 h − Mpc on a side. This corresponds http: // / millennium / c (cid:13) , 000–000 ultiple Black Holes Simulation Mass of halo at z = ⊙ ) Max. BH no. SMBH Coalescences SMBH EscapesL1 1 . × . × . × . × . × . × . × . × Table 1.
Summary of simulations and results for haloes that have a mass of M ∼ M ⊙ . The maximum BH number denotes the number of black holes inthe biggest BH group found in a simulation. The last two columns show number of BH coalescences and escapes in the simulation. A halo with M = M ⊙ has average mass M z = = × M ⊙ .Simulation Mass of halo at z = ⊙ ) Max. BH no. SMBH Coalescences SMBH EscapesH1 1 . × . × . × . × . × . × . × . × . × . × . × . × . × . × . × . × . × Table 2.
Summary of simulation runs with haloes that have mass M & M ⊙ at z =
0. Various columns are same as Table 1. A halo with M = M ⊙ has average mass M z = = × M ⊙ . Figure 4.
Evolution of single and binary SMBHs in our simulations. (a) The left hand panel shows evolution of the x -component of the position of a 9 . × M ⊙ back hole near the centre of a Hernquist bulge of mass 5 . × M ⊙ and scale length of 0 . . × M ⊙ . The secular motionis due to that of the cusp. (b) The right hand panel shows evolution of the separation between SMBHs in a binary with initial separation 2 kpc and eccentricity0 .
5. The black hole masses were 8 . × M ⊙ and the binary evolved near the center of a Hernquist halo with mass 5 . × M ⊙ and scale length of 10 . . × M ⊙ .c (cid:13) , 000–000 Kulkarni & Loeb to a particle mass of 8 . × h − M ⊙ . The output of this simula-tion is stored in 64 snapshots between z =
127 and z =
0. Particlesin each snapshot are grouped into friends-of-friends (FOF) clus-ters that are expected to correspond to virialised structures. EachFOF halo contains substructure of gravitationally bound subhaloesthat can be related to each other across snapshots as progenitorsand descendants. Because a halo can contain multiple galaxies, weexpect the subhalo merger tree to reflect the merger history of thegalaxies within a halo. Since the goal of this paper is to under-stand formation and evolution of systems of multiple black holesdue to the hierarchical merger history of a galaxy, we extract sub-halo merger trees from the Millennium Simulation Database. Eachsuch merger tree typically shows growth of a subhalo via accretionof dark matter particles and via mergers. We process these mergertrees to keep only major mergers, which we define to be mergershaving mass ratio larger than 0 .
1. To identify the mass ratio of twosubhaloes, we use the masses of the distinct FOF haloes that thesesubhaloes were a part of before the FOF haloes merged. This is toaccount for the mass loss of the satellite subhalo due to tidal strip-ping after it enters the FOF group of the host subhalo, but beforethe eventual merger of the two subhaloes. (See discussion in § ff erence being our use of direct-summation N-body simula-tions instead of SPH simulations.Once we have a galaxy merger tree, we set up the initial condi-tions of our simulation in the “leaves” of the tree, that is, in haloesthat do not have a progenitor, and follow the evolution using an N-body calculation. The initial conditions of our simulation consist ofa stellar spheroid with a Hernquist density profile, ρ ( r ) = M π ar ( r + a ) , (30)where M is the total mass of the spheroid and the scale length a isrelated to the half mass radius r / of the spheroid by a = . r / .Values for the parameters M and a were obtained from the halomass as follows (Ho ff man & Loeb 2007). We first obtain the blackhole mass M bh from the halo mass M halo using Equation (29).We then use the empirical relation between the SMBH mass andthe spheroid’s virial mass (Magorrian et al. 1998; Marconi & Hunt2003; Peng et al. 2006) to obtain the latter as M sph = . × M ⊙ " M bh M ⊙ . . (31)The virial mass of the spheroid is related to its velocity dispersion σ e and half light radius R e by M sph = kR e σ e G . (32)We follow Marconi & Hunt (2003) and set k = R e / R e . These two methods are in essential agreement, as argued by Tremaine et al. (2002). As-suming a constant mass-to-light ratio for the Hernquist profile, wehave R e = . a and the velocity dispersion at radius R e / σ e = . GMa . (33)This lets us obtain the value of the parameter M of the Hernquistprofile as M = . M sph . The scale length a is readily obtained as a = GM sph κ σ , (34)where σ bh is obtained using the M − σ relation of equation (28).Having obtained a density profile for the bulge, we place a blackhole at its center and set the black hole mass to be ten times thatobtained from equation (29). This factor of ten is introduced tokeep the ratio between the black hole mass and the particle masshigh enough (Milosavljevi´c & Merritt 2001; Makino & Ebisuzaki1996). We confirm that the radius of influence r inf = Gm bh /σ ofthis black hole is still much smaller than the a . Velocities of the starsin the spheroid are then generated from the unique, isotropic ve-locity distribution that corresponds to the gravitational potential ofthe density profile in Equation (30) and the SMBH (Tremaine et al.2002). These initial conditions are then scaled to standard N-bodyunits of G = M = E = − .
25, where M is the total massof the system and E is its total energy (Heggie & Mathieu 1986;Aarseth 2003). In these units, in virial equilibrium, the mean squarevelocity h v i = / t cr = √
2, in-dependent of the number of particles. The conversion factors fromphysical units to these N-body units can be easily obtained via di-mensional analysis.Note that we ignore presence of gas in this set-up. Simu-lations of binary BHs in gaseous environment have not reachedsu ffi cient resolution to establish the role played by gas in evolu-tion of SMBHs in galactic nuclei (Merritt & Milosavljevi´c 2005;Colpi & Dotti 2009). Moreover, we expect that at high redshifts,quasar activity triggered by galaxy mergers could e ffi ciently drivegas away from the shallow potential wells of the galaxies.To perform the actual dynamical evolution of this system,we use the direct-summation code NBODY6 written by SverreAarseth (Aarseth 1999, 2003). This code has been well-tested forvarious applications since around 1992. Its purpose is to performan exact integration, without particle softening, of a large num-ber of particles. It integrates equations of motion of individualparticles using a fourth-order Hermite method with block timesteps (Makino & Aarseth 1992). This integrator is coupled with theAhmed-Cohen neighbour scheme (Ahmad & Cohen 1973), whichselects a subset of neighbours of a particle whose forces on it arecalculated at a higher time resolution that other, more distant, parti-cles. This scheme reduces the computational cost from O ( N ) toabout O ( N . ). Close two-body encounters are treated using theKustaanheimo-Stiefel (KS) regularization method that eliminatesthe r = η I , and the time-stepparameter for regular force polynomial, η R are set to 0.02. The en-ergy tolerance is set to Q E = × − and the regularized time-stepparameter is set to η U = . c (cid:13) , 000–000 ultiple Black Holes We check the stability of our initial conditions by evolvingstandalone realizations of the Hernquist bulge with a central BHand then traverse the merger tree of a given halo using NBODY6,starting from the initial conditions as described above. We scale thephysical time between two successive nodes of the tree to N-bodyunits and run NBODY6 for that duration. If a merger happens ata certain node, we place the two galactic nuclei at a distance of 2kpc apart and evolve in an head-on approach. Although such head-on mergers would be unlikely, we choose it to reduce the compu-tational time while still retaining some realism. When two galax-ies, that are in equilibrium separately, merge we expect some tran-sient response in the resulting dynamics. However, as discussed byMilosavljevi´c & Merritt (2001), any such e ff ects in the dynamicsof the central regions of the merger remnant of these galaxies areessentially negligible.Under these conditions, the component black holes approachafter a merger event and the remnant galactic nucleus is left withtwo black holes, which gradually harden due to dynamical frictionand three-body interactions with stars in their vicinity. Black holecoalescence is implemented in our simulation by monitoring theseparation of hard black hole binaries. Once members of a SMBHbinary get closer than a fixed distance d crit , we replace them witha single black hole with mass equal to the sum of the masses ofcomponent black holes. In all the runs reported in this paper, we set d crit = . M − σ relation, but the later growth of theseSMBHs occurs only via coalescence.Recoil due to anisotropic emission of gravitational waves isa natural consequence of asymmetric merger of black holes, ei-ther due to unequal masses or due to unequal spins (Peres 1962;Bekenstein 1973). Until recently, it was unclear whether this re-coil is large enough to be astrophysicaly relevant. However, recentresults from numerical relativity have revealed the resultant kickvelocities in a variety of merger configurations (Pretorius 2005;Baker et al. 2006). When the black hole spins are aligned with eachother and with the orbital spin, these simulations find revoil ve-locity of v recoil .
200 km s − (Baker et al. 2006; Gonz´alez et al.2007; Herrmann et al. 2007; Lousto & Zlochower 2009). In the ab-sence of spins, this recoil velocity is only a function of the ratio ofblack hole masses. For random orientations of spins, recoil veloc-ities as high as 2000 km s − have been obtained (Campanelli et al.2007a,b). Bogdanovi´c et al. (2007) argue that a circumbinary gasdisk can align the binary spins with the orbital axis thereby reduc-ing v recoil to about 200 km s − . In our simulations we assume a con-stant kick velocity of 200 km s − , which we impart to the remnantof every unequal-mass binary SMBH coalescence.We follow the approach of Makino & Aarseth (1992) and keepthe particle number fixed at N = throughout the simula-tion. Thus, at every merger, we combine particles in each merg-ing galactic nucleus and double the particle mass. This lets uskeep the particle number high throughout the merger tree of thehalo. The ratio of black hole mass to the stellar mass is typicallya few hundred, which is also roughly the ratio of the spheroid’stotal mass to the black hole’s mass. These values are compara-ble to other simulations of this kind (Makino & Ebisuzaki 1996;Milosavljevi´c & Merritt 2001).In summary, the unique features of our simulations are: (i) kinematically consistent initial conditions with black holes; (ii) cal-culation of mergers of galactic nuclei in a cosmological setting us-ing merger trees extracted from cosmological N-body simulations; (iii) calculation of merger of galactic nuclei resulting in a formation Figure 5.
Number of black holes as a function of redshift in a simulationwith M = . × M ⊙ . of SMBH binaries starting from the results of each nucleus havingevolved in isolation; and (iv) accurate calculation of SMBH-starand SMBH-SMBH dynamics throughout the assembly history of agalactic nucleus and its constituent SMBH with the e ff ect of gravi-tational wave recoil taken into account. We perform some basic checks on our code, such as ensuring en-ergy conservation and stable evolution of equilibrium systems. Inall of our runs, the relative error in the total energy is maintained at | ∆ E / E | < × − . The treatment of BH-BH and BH-star interac-tion is handled by the original nbody nbody ff ect. These three runs are excluded from the re-sults presented below. In a stellar environment, a single SMBH exhibits a random fluctuat-ing motion arising due to discrete interactions with individual stars.As a result, the e ff ect of the stellar environment on the SMBH canbe decomposed into two distinct components: (1) a smooth compo-nent arising due to the large scale distribution of the whole system,and (2) a stochastic fluctuating part coming form the interation withindividual stars (Chatterjee et al. 2002). This random motion is il-lustrated in the left hand panel of Figure 5, which shows evolutionof the x -component of the position of a 9 . × M ⊙ back holenear the centre of a Hernquist bulge of mass 5 . × M ⊙ and scalelength of 0 . . × M ⊙ . As expected,the SMBH wanders around due to stochastic interactions with thestars in its vicinity. The mean square amplitude of these fluctuations c (cid:13) , 000–000 Kulkarni & Loeb
Figure 6.
Histograms of ejection velocities of BHs. Left: Velocities of ejected black holes in all of our high mass runs. Note that this does not include ejectedblack holes with the highest velocities ( > − ). Right: number of ejections as a function of redshift in our high mass runs. is expected to be (Chatterjee et al. 2002; Milosavljevi´c & Merritt2003a) h x i ≈ m ∗ m BH r , (35)where r core is the radius within which the stellar distribution flat-tens out. The Hernquist distribution that we have used here doesnot have a well-defined core, since the density keeps rising as r − near the origin. Milosavljevi´c & Merritt (2003a) argue that the ef-fective core radius for such distribution can be taken as the ra-dius of influence of the black hole. The resultant mean squarevalue of fluctuations is somewhat smaller that that for Figure 5 byroughly a factor of 2 as is known to happen in N-body simulations(Quinlan & Hernquist 1997; Milosavljevi´c & Merritt 2003a).As described in §
1, the evolution of a binary black hole ina gas-poor galaxy takes place in three stages. Right hand panelof Figure 5 shows evolution of the separation between SMBHs ina binary with initial separation 2 kpc and eccentricity 0 . . × M ⊙ and the binaryevolved near the center of a Hernquist halo with mass 5 . × M ⊙ and scale length of 10 . . × M ⊙ . In the first stage of evolution, the SMBHs sink to the cen-tre of the galactic nucleus by losing energy via dynamical frictionand become bound to each other. This stage ends when the sep-aration between the SMBHs is equal to the radius of influence ofthe binary (Merritt & Milosavljevi´c 2005). In the second evolution-ary stage, the binary loses energy predominantly ejection of nearbystars via three-body interaction. The binary loses energy rapidly inthis stage, which continues until t ≈
200 Myr for the case depictedin Figure 5. The final stage of the SMBH binary evolution beginswhen the rapid hardening of the second stage stops. This happenswhen the binary semi-major axis takes the value given by Equation(1). The binary semi-major axis is related to the separation r by1 a = r − v µ , (36)where v is the relative velocity of the BHs and µ is the reduced mass(Makino & Funato 2004; Berczik et al. 2006; Merritt et al. 2007;Khan et al. 2011). In N -body simulations, the last stage is known tohave a dependence on the number of particles N such that the hard-ening rate decreases with increasing N (Makino & Funato 2004). Figure 7.
Escape velocities from the bulges of haloes in our three cate-gories of present-day masses of haloes. Solid line: M ≈ M ⊙ , Dashedline: M ≈ M ⊙ , Dot-dashed line: M & M ⊙ . Note that these areaverage values computed from the fitting functions to the Millennium sim-ulation. Therefore, case by case comparison with our runs is not straight-forward. For real spherical galaxies, the binary separation would stop evolv-ing after this point because the loss cone is empty.
We now run the simulation along merger trees of haloes drawn fromthe Millennium simulation as described in Section 5. These simula-tions are described in Tables 1 and 2. We randomly select 8 haloeswith mass M around 10 M ⊙ at z =
0. These correspond to thetypical haloes ( M ≈ M ∗ ) in the present epoch. We also randomlyselect 17 haloes whose present-day mass M is in excess of 10 M ⊙ . These are rare, high mass haloes that are expected to host theredshift 6 SDSS quasars (Li et al. 2007). Additionally, we have also c (cid:13) , 000–000 ultiple Black Holes Figure 8.
Number of black holes as a function of redshift in a few of our simulation runs.
Figure 9.
Projected stellar density contours in the presence of a binary in the simulation H5. Each panel is 400 pc on a side. Clockwise from top left to bottomright, the redshifts are z = . .
54, 7 .
27, 6 .
19, 5 .
28, and 4 .
52. The total time span is about 800 Myr. Core-SMBH oscillations are clearly visible.c (cid:13) , 000–000 Kulkarni & Loeb
Figure 10.
Projected stellar density contours in the presence of multiple BHs in the simulation H4. Each panel is 100 pc on a side. The total time span,clockwise from top left to bottom right, is about 1 Gyr. Most BHs are stripped of their cusps in nuclei with multiple BHs.
Figure 11.
The fraction of runs with multiple SMBHs at di ff erent redshift bins for haloes with a mass M ∼ M ⊙ at z =
0. The results of these runs aresummarised in Table 1. The three panels from left to right describe the occurrence of systems with more than 2, 3 and 4 black holes respectively. At eachredshift, this number can be interpreted as the likelihood of finding such systems in haloes of mass M ∼ M ⊙ at z =
0. It is seen that systems with multipleSMBHs are rare at redshift z .
2. Note that a halo with M = M ⊙ will have M z = = × M ⊙ . simulated 11 haloes with present-day mass similar to the MilkyWay halo ( M ∼ M ⊙ ).Using the prescriptions described in theprevious section, and using the N-body integrator, these simulationstell us about the e ff ect of multiple mergers of galactic nuclei withSMBHs.Figure 5 shows results from a typical simulation run, for ahalo of mass 1 . × M ⊙ . We plot here the number of BHs inthe bulge in the main branch of the galaxy’s merger tree at vari-ous redshifts. It is seen that the central bulge has more than oneSMBH for a wide redshift range (2 . z .
6; about 2 . . z . z = z &
6) there are no BHs in the central bulge. This is sim-ply an artifact of the limited numerical resolution of the Millenniumsimulation, because of which the halo merger tree is not resolvedat these redshifts. To ensure that this does not a ff ect our results for z .
6, we set up initial conditions at z ∼ M − σ relation, and by using a Hernquist bulge with inner slope c (cid:13) , 000–000 ultiple Black Holes Figure 12.
The fraction of runs with multiple SMBHs at di ff erent redshift bins for halo masses M & M ⊙ at z =
0. The results of these runs are summarisedin Table 2. The three panels from left to right describe the occurrence of systems with more than 2, 3 and 4 black holes respectively. At each redshift, thisnumber can be interpreted as the likelihood of finding such systems in haloes of mass M & M ⊙ at z =
0. It is seen that systems with multiple SMBHsare rare at redshift z .
2. These results can be compared with those in figure 11. Nuclei with multiple SMBHs are more likely in high mass haloes because ofhigher merger rate. Note that a halo with M = M ⊙ will have M z = ∼ M ⊙ . −
1. In the absence of gas, the systems with multiple SMBHs formgenerically, in high mass haloes with frequent of major mergers. Itis evident than such systems are usually short-lived and most oftenthese nuclei contain a single SMBH at z =
0. Most SMBHs es-cape into the halo, where they join a population of wandering blackholes or escape the halo completely.Similar results from a few other simulation runs for haloeswith mass M ∼ M ⊙ at z = − z = − , which is the GW recoil kick inour simulations. Note that this plot does not show kicks with veryhigh velocities, which we describe below.With the prescription that we have adopted in this paper, wefind that SMBH coalescence happens in each one of our simu-lations. Tables 1 and 2 give the number of BH coalescences oc-curring in our simulations. Due to the limitation on the particle number, our simulations implement BH coalescence by replacinga bound binary BH by a single BH whose mass is equal to the to-tal mass of the binary. As an example, Figure 9 shows the mergerof two bulges beginning from initial conditions at redshift 6.7 inthe run H5. In Figure 9, the hardening radius is a h = . t h =
500 Myr. We find the the BHs remain associated with theirhost cusps until cusp coalescence. It is known that by increasingthe e ff ective mass of the BHs, this increases the rate of coales-cence of the BHs by as much as ∼ ff erence from previous works is in the evolution of the densityprofile in the later stages of the merger. In our simulations, each co-alescence event is followed by recoil of the remnant at 200 km s − ,which at high redshift, usually results in the escape of the SMBHfrom the galaxy. At relatively low redshifts, the recoiled SMBH re-turns to the nucleus in few hundreds Myr. Because of this recoil,the remnant BH is detached from its cusp immediately. At the re-coil speed implemented here, this happens at a much smaller timescale that the local crossing time scale. As a result, the only e ff ectof the remnant on the cusp is due to subsequent core passages.Usually, most coalescences are assumed to take place due toBH hardening via BH-star encounters. In gas-free systems, thisleads to the final parsec problem. In our simulations, we find that inhigh mass haloes, roughly half of the SMBH coalescences are dueto three-body scattering with intruder SMBHs. This is expected,since in spite of higher major merger rate, high mass galaxies in ourmodel are still left with at most two SMBHs at z =
0. The dominantmechanism of coalescence is then three body interactions. Figure8 shows an example of the evolution of a multiple BH system thatundergoes three coalescences due to BH-BH dynamics. We findviolent oscillations of the cusp-BH system as shown in Figure 9.This has significant impact on the density distribution of the core,and also results in o ff -centre BHs, which slowly return to the centreof the cusp due to dynamical friction.About 10% of SMBH ejections in our simulations occur atvery high speeds of & − . In haloes with M ≈ M ⊙ these SMBHs will linger in the outskirts of the halo for 2 − c (cid:13) , 000–000 Kulkarni & Loeb
Figure 13.
Evolution of density profile for simulations H3 and H5 in N-body units. The solid line is the original Hernquist profile with an inner logarithmicslope of γ ≈ −
1. Dashed line shows the profile after one SMBH binary coalescence, dot-dashed line after the second coalesence and the dotted line after thethird coalescence. These plots are shown in N-body units to scale out the doubling of the half-mass radius. See text for details.
Gyr as can be seen by comparing with the bulge escape speedsin Figure 7. The SMBHs in the wandering phase that are intro-duced via this mechanism have markedly di ff erent properties thanthe BHs introduced due to galaxies that have not yet reached thehost galaxy’s center so as to have a close encounter (Volonteri et al.2003). The main di ff erence is that our ejected black holes are muchmore massive than those in the other category. Moreover, the veloc-ity of ejected SMBHs will typically be higher that black holes in theother category, which have already experience significant dynam-ical friction. Three of the 30 BH ejections in our runs are ejectedbinaries. From the results of our simulations, we can estimate the likelihoodof galactic nuclei with multiple black holes at high redshifts. Thehistograms in Figures 11 and 12 show fraction of runs with multipleSMBHs at each redshift for haloes with present-day masses of ∼ M ⊙ and ∼ M ⊙ , respectively. The three panels from leftto right describe the occurrence of systems with more than 2, 3and 4 black holes respectively. At each redshift, this number can beinterpreted as the likelihood of occurrence of such systems at thatredshift.Systems with more than 2 SMBHs are generically expected inthe central galaxies of haloes with M & M ⊙ at around z & z . ffi cient time to coalescence. This is consistentwith the expectation from our heuristic analysis of Section 3. Inother words, multiple black hole systems are numerous at aroundredshifts of 6, when there are many major mergers in the system.Our numerical simulations show that such systems can exist in suf-ficiently long-lived configurations of SMBHs separated on pc–kpcscale. Note that these histograms show the likelihood of such sys-tems to be zero at redshifts z &
10. However, this is simply becausethe Millennium simulation merger trees do not resolve progenitorsat these redshifts. As mentioned before, we have minimized the ef- fect of this shortcoming on our results by requiring that the SMBHsalways follow the M − σ relation initially.High mass galaxies ( M ≈ M ⊙ ) are more likely to havemultiple BHs in their nuclei at higher redshift. About 60% of thesegalaxies have more than 2 BHs between redshifts z ≈ M ≈ M ⊙ )The likelihood of occurrence of more than 3 and 4 BHs is similar,about 30%, in the two categories of simulation. However, for thehigh mass galaxies this likelihood is spread out over a wider rangein redshift, again due to the higher rate of major mergers.It is extremely rare for Milky Way-sized galaxies (halo mass M ≈ M ⊙ ) to have more than three SMBHs in their nuclei atany moment in their assembly history. Indeed, in our simulationsof these galaxies, only one run shows a triple BH system. The mainreason behind this is the smaller number of major mergers for thesegalaxies. Moreover, it is easier for SMBHs to escape the nuclei ofpredominantly small mass progenitors of these galaxies. ff ects on the stellar distribution Most bulges and early-type galaxies have a shallow cusp near theircentre. The mass distribution in this region can be described asa power law ρ ∝ r − γ . Most galaxies have slope 0 . . γ . . ff ect the mass distributionwithin its radius of influence. Only two galaxies, the Milky Way(Genzel et al. 2003) and M32 (Lauer et al. 1998), have been re-solved at these small distances. Both these galaxies have γ ≈ . M def , which is the di ff erencebetween the mass of the initial and final density distribution in a re-gion around the centre, is roughly M bh , the total mass of the SMBHbinary. The possibility of enhanced mass deficit because of re-peated core passages of recoiled black holes (Gualandris & Merritt2008) and due to repeated mergers (Merritt 2006) has been consid- c (cid:13) , 000–000 ultiple Black Holes Figure 14.
Mass deficiency versus number of coalescences averaged overten simulation runs. The presence of multiple SMBHs generally leads tolarger mass deficiency compared to a single hard SMBH binary. ered in the literature. Our simulations allow us to understand thee ff ect of both of these factors in addition to the mass deficit pro-duced by simultaneous presence of multiple SMBH in the galacticbulge.Figure 6.2 shows the cusp evolution in two of our simulations,each of which has four SMBHs and three coalescences. Densityprofiles after each coalescence is shown. Strong core formation isclearly seen. We calculate M def / M bh for ten such runs and showthe average result in Figure. Clearly M def / M bh is much larger whenmultiple SMBHs are present. Values of M def / M bh ≈ ∼ yr, we can expect themto carry the signature of core formation at high redshift due to mul-tiple SMBHs. At lower redshift our simulation are applicable tospiral bulges, which have much lower relaxation time scale ( ∼ yr). Indeed in the runs where a single black hole is left for z .
2, wefind the formation of a Bahcall-Wolf cusp. This is consistent withthe observed structure of the Milky Way bulge.The above considerations regarding cores in galaxy luminosityprofile are also applicable to dark matter cores. The ejection of darkmatter particles by the black holes will produce a core similar insize to the stellar core.
From the results of our simulations described above, we expectabout 30% of the galaxies within haloes with a present-day massof M ≈ M ⊙ to contain more than two SMBHs at redshifts2 . z .
6. For more massive haloes with M & M ⊙ , this frac-tion is almost 60%. However, since few such systems have beenunambiguously observed so far, we consider some observationalsignatures that would indicate their existence . Apart from their ef-fect on the stellar mass distribution, multiple SMBH systems lead Some systems with triple active galactic nuclei (AGNs) were reported sofar. Examples are NGC 6166 and 7720 (Tonry 1984) and SDSSJ1027 + to an enhanced rate of tidal disruption of stars, modified gravita-tional wave signals compared to isolated BH binaries, and slingshotejection of SMBHs from galaxies at high speeds.From the results of scattering experiments, Chen et al. (2009)found that the stellar tidal disruption rates due to three-body in-teractions between a hard, unequal-mass SMBH binary with fixedseparation and a bound stellar cusp is higher by several orders ofmagnitude than the corresponding rates for a single SMBH. In par-ticular, they find that the stellar tidal disruption rate is about 1 yr − for an isothermal stellar cusp with σ =
100 km s − containing anSMBH binary of total mass 10 M ⊙ . In comparison, the correspond-ing rate for a single 10 M ⊙ black hole is about 10 − yr − . The dura-tion of the tidal disruption phase is about 10 yr. This enhancementin the tidal disruption is due to the Kozai-Lidov e ff ect and due tochaotic resonant scattering (Chen et al. 2011). Tidal disruption of astar results in about half of the stellar mass being inserted in boundelliptical orbits. When it falls back in the black hole, this mass givesrise to a bright UV / X-ray emission (“tidal flare”) lasting for a fewyears. One such event may have already been recently observed inthe form of high-energy transients that can be modeled as suddenaccretion events onto an SMBH (Levan et al. 2011; Bloom et al.2011; Zauderer et al. 2011).We expect similar enhancement in the rate of stellar tidal dis-ruption in systems with multiple black holes. Firstly, the presenceof multiple SMBHs increases the combined tidal disruption crosssection of the black holes. (Although this will only enhance thetidal disruption rate by a factor of a few.) Secondly, even beforethey closely interact, the presence of a third SMBH a ff ects the tidaldisruption event rate onto an SMBH binary by scattering stars intothe binary’s loss cone at a rate that increases as inverse square of itsseparation from the binary (Ho ff man & Loeb 2007). Thirdly, as wesaw above, multiple SMBH systems are likely to contain recoiledblack holes, which have been kicked either due to anisotropic grav-itational wave emission after coalescence, or due to the gravita-tional slingshot. Sudden recoil promptly fills the loss cone of theseblack holes. The resultant enhancement in the tidal disruption eventrate can be substantial, increasing it up to 0.1 yr − (Stone & Loeb2011). Furthermore, if their recoil velocity is not too high, theserecoiled SMBHs oscillate around the stellar core with decreasingamplitude due to dynamical friction. This motion results in theirrepeated passages through the stellar core, thereby increasing thestellar tidal disruption event rate.Another observational signature of systems with multipleSMBHs is gravitational waves (GWs). The GW emission frombinary and triple SMBHs has been studied in the literature(Wyithe & Loeb 2003a; Sesana et al. 2004; Amaro-Seoane et al.2010). Space-based detectors like the Laser Interferometer SpaceAntenna (LISA) are expected to be sensitive in the frequencyrange ∼ − –10 − Hz. This corresponds to the inspiral of SMBHsystems with total mass ∼ − M ⊙ . Pulsar timing arrays(PTAs) like the Parkes PTA (Manchester 2008) and the Euro-pean PTA (Janssen et al. 2008) and ground-based detectors like theNorth American Nanohertz Observatory for Gravitational Waves(Jenet et al. 2009) are sensitive to even lower frequencies of ∼ − –10 − Hz.Yunes et al. (2011) studied modifications due to the presence (Liu et al. 2011). The first two objects are cD galaxies at z ≈ .
03 andthe latter is at z ≈ .
06. All three are kpc-scale triples. It is possible thatNGC 6166 is simply a superposition of a central cD galaxy and two low-luminosity elliptical galaxies (Lauer et al. 1998).c (cid:13) , 000–000 Kulkarni & Loeb of a secondary SMBH in the waveform of an extreme mass-ratioinspiral (EMRI) of a stellar mass objects into an SMBH. Theyfind that a 10 M ⊙ SMBH will produce detectable modificationsif it is within a few tenths of a parsec from the EMRI system, al-though this distance increases for higher mass SMBHs. In this pa-per, we have quantified the presence of such ‘massive perturbers.’The resultant modifications to gravitational waveforms will be adistinct signature of multiple-SMBH systems. Futhermore, suchsystems often contain binaries that have phases of very high ec-centricities, created via mechanisms like the Kozai-Lidov e ff ect(Ho ff man & Loeb 2007). Such binaries are expected to to emit in-tense bursts of high-frequency gravitational waves at the orbital pe-riapsis (Amaro-Seoane et al. 2010). As a result, sources that wouldnormally emit outside of the frequency windows of planned grav-itational wave searches may be shifted into observable range. Forexample, Amaro-Seoane et al. (2010) find that a few to a hundredgravitational wave bursts could be produced at a detectable (1 ns)level within the PTA frequency range if the fraction of SMBHtriplets is > . ff man & Loeb 2007). We have shown that about 10% of theSMBHs are ejected at velocities > − due to the sling-shot mechanism. This high-speed black holes will spend 1 −
10 Gyrin the outskirts of the halo. However, it is not clear whether detect-ing this population of wandering black holes will be possible.
In this work, we have addressed the formation of galactic nucleiwith mutiple SMBHs. We performed accurate N-body simulationsof mergers of galactic nuclei with SMBHs in a cosmological set-ting. Our calculation uniquely incorporated cosmological mergersof galaxies with an accurate treatment of dynamical interactions be-tween SMBHs and stars, which we achieved using the direct sum-mation N-body code, NBODY6. The need for such simulations hasbeen recognized in the literature (Merritt & Milosavljevi´c 2005).Our main conclusions are as follows: • In the absence of gas, high mass galaxies ( M & M ⊙ at z =
0) are generically expected to have had multiple SMBHsin their nuclei during their assembly history. Our simulations sug-gest that ∼
30% galaxies within haloes with a present-day mass of M ≈ M ⊙ ( M z = ≈ M ⊙ ) contain more than two SMBHs atredshifts 2 . z .
6. For more massive haloes, with M & M ⊙ ( M z = ≈ M ⊙ ), this fraction is almost 60%. This is in contrastto lower-mass galaxies ( M ≈ M ⊙ ; M z = ≈ M ⊙ ), whichrarely host more than two SMBHs in their nuclei at any moment intheir assembly history. • High mass galaxies as well as their low mass counterpartsare rarely expected to retain more than two SMBHs in their nucleiat the present epoch. SMBH coalescence and ejection reduces thenumber of SMBHs on the time scale of a Gyr. Furthermore, majormergers are rare at lower redshift. We also find that the number ofSMBHs in galactic nuclei is rarely reduced to zero at z =
0. Lessthan 5% of our high-mass runs resulted in such galaxies. • SMBH coalescence is common at high redshifts. Subsequentrecoil due to anisotropic gravitational wave emission often results in escaping SMBHs. Some of these SMBHs add to the wander-ing population of black holes in the galactic halo. In a few cases,this process also results in galactic nuclei with no SMBH near theircentres. BH-BH interaction also leads to ejected SMBHs via theslingshot mechanism. While most of ejected SMBHs have veloc-ities .
500 km s − , about 10% SMBHs are ejected at very highvelocities exceeding 2000 km s − . We also find binary SMBH ejec-tion in .
10% of the cases. • Multiple SMBHs have a strong e ff ect on the stellar distribu-tion due to three-body interactions and core passages. The resultingmass deficit is usually much larger than that due to a single SMBHbinary because of resonant BH-BH interactions and GW recoil ofthe BH remnant. We observe long-term oscillations of the BH-coresystem that could explain observations of o ff set AGNs. This hasimplications for recent observations by Civano et al. (2010) of a z = .
359 system that potentially contains a recoiled BH. • The presence of multiple SMBHs will have important e ff ectson the rate of tidal disruption of stars in galactic nuclei due to en-hanced tidal disruption cross section, scattering of stars by otherBHs, prompt loss cone refilling due to GW recoil and gravitationalslingshot. Similarly, the presence of more than two BHs in a hier-archical triple is expected to leave a signature in the GW emissionfrom the inner binary. This signature could be observed with futureGW observatories, such as LISA. Finally, we also expect such sys-tems to give rise to a distinct population of wandering SMBHs thatcould travel in large haloes over long time scales of a few Gyrs.The presence of gas could alter the above picture to some ex-tent. However, simulations of binary BHs in gaseous environmenthave not reached su ffi cient resolution to confirm this. Moreover,we expect that at high redshifts, AGN activity triggered by galaxymergers could e ffi ciently drive gas away from the shallow potentialwells of the galaxy. Our work can also be extended by calculatinglate stages of binary SMBH evolution more consistently. New reg-ularization techniques to do this are now available (Aarseth 2003);we defer their use to future work. Furthermore, multiple SMBHsystems can also form in additional ways, for example by fragmen-tation of disks (Goodman & Tan 2004). However, these systemswould evolve by migration (Kocsis et al. 2011) on a much shortertime scale than considered here. ACKNOWLEDGEMENTS
We acknowledge advice on various aspects of our simulations fromSverre Aarseth and would like to thank him for making his codesavailable. GK also acknowledges discussion with Jasjeet S. Bagla,Hagai Perets and Yue Shen, and the Institute for Theory and Com-putation for hospitality. Our simulations were run on the Odysseycluster supported by the FAS Sciences Division Research Com-puting Group at Harvard University. The Millennium Simulationdatabases used in this paper and the web application providingonline access to them were constructed as part of the activitiesof the German Astrophysical Virtual Observatory. This researchwas supported in part by a Fulbright-Nehru Professional and Pre-doctoral Fellowship from the US-India Educational Foundation, byNSF grant AST-0907890 and NASA grants NNX08AL43G andNNA09DB30A.
REFERENCES
Aarseth S. J., 1999, PASP, 111, 1333 c (cid:13) , 000–000 ultiple Black Holes Aarseth S. J., 2003, Gravitational N-Body SimulationsAhmad A., Cohen L., 1973, Journal of Computational Physics,12, 389Amaro-Seoane P., Freitag M. D., 2011, MNRAS, 412, 551Amaro-Seoane P., Sesana A., Ho ff man L., Benacquista M., Eich-horn C., Makino J., Spurzem R., 2010, MNRAS, 402, 2308Armitage P. J., Natarajan P., 2002, ApJ, 567, L9Baker J. G., Centrella J., Choi D.-I., Koppitz M., van Meter J.,2006, Physical Review Letters, 96, 111102Baker J. G., Centrella J., Choi D.-I., Koppitz M., van Meter J. R.,Miller M. C., 2006, ApJ, 653, L93Barkana R., Loeb A., 2001, Phys. Rep., 349, 125Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287,307Bekenstein J. D., 1973, ApJ, 183, 657Benson A. J., 2005, MNRAS, 358, 551Berczik P., Merritt D., Spurzem R., Bischof H.-P., 2006, ApJ, 642,L21Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edi-tion. Princeton University PressBlaes O., Lee M. H., Socrates A., 2002, ApJ, 578, 775Blecha L., Cox T. J., Loeb A., Hernquist L., 2011, MNRAS, 412,2154Bloom J. S., Giannios D., Metzger B. D., Cenko S. B., PerleyD. A., Butler N. R., Tanvir N. R., Levan A. J., O’ Brien P. T.,Strubbe L. E., De Colle F., Ramirez-Ruiz E., Lee W. H., Nayak-shin S., Quataert E., et al., 2011, ArXiv e-printsBogdanovi´c T., Reynolds C. S., Miller M. C., 2007, ApJ, 661,L147Boylan-Kolchin M., Ma C., Quataert E., 2008, MNRAS, 383, 93Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A.,Lemson G., 2009, MNRAS, 398, 1150Bromm V., Loeb A., 2003, ApJ, 596, 34Bundy K., Treu T., Ellis R. S., 2007, ApJ, 665, L5Campanelli M., Lousto C., Zlochower Y., Merritt D., 2007a, ApJ,659, L5Campanelli M., Lousto C. O., Zlochower Y., Merritt D., 2007b,Physical Review Letters, 98, 231102Chandrasekhar S., 1943, ApJ, 97, 255Chatterjee P., Hernquist L., Loeb A., 2002, ApJ, 572, 371Chen X., Madau P., Sesana A., Liu F. K., 2009, ApJ, 697, L149Chen X., Sesana A., Madau P., Liu F. K., 2011, ApJ, 729, 13Civano F., Elvis M., Lanzuisi G., Jahnke K., Zamorani G., BlechaL., Bongiorno A., Brusa M., Comastri A., Hao H., LeauthaudA., Loeb A., Mainieri V., Piconcelli E., Salvato M., Scoville N.,Trump J., Vignali C., et al., 2010, ApJ, 717, 209Colpi M., Dotti M., 2009, ArXiv e-printsColpi M., Mayer L., Governato F., 1999, ApJ, 525, 720Eisenstein D. J., Loeb A., 1995, ApJ, 443, 11Escala A., Larson R. B., Coppi P. S., Mardones D., 2004, ApJ,607, 765Escala A., Larson R. B., Coppi P. S., Mardones D., 2005, ApJ,630, 152Faber S. M., Tremaine S., Ajhar E. A., Byun Y.-I., Dressler A.,Gebhardt K., Grillmair C., Kormendy J., Lauer T. R., RichstoneD., 1997, AJ, 114, 1771Fakhouri O., Ma C., Boylan-Kolchin M., 2010, MNRAS, 406,2267Fan X., Narayanan V. K., Lupton R. H., Strauss M. A., KnappG. R., Becker R. H., et al., 2001, AJ, 122, 2833Ferrarese L., 2002, ApJ, 578, 90Ferrarese L., Cˆot´e P., Jord´an A., Peng E. W., Blakeslee J. P., Piatek S., Mei S., Merritt D., Milosavljevi´c M., Tonry J. L., West M. J.,2006, ApJS, 164, 334Ferrarese L., Ford H., 2005, Space Sci. Rev., 116, 523Ferrarese L., Merritt D., 2000, ApJ, 539, L9Freitag M., Amaro-Seoane P., Kalogera V., 2006, ApJ, 649, 91Gebhardt K., Bender R., Bower G., Dressler A., Faber S. M., Fil-ippenko A. V., Green R., Grillmair C., Ho L. C., Kormendy J.,Lauer T. R., Magorrian J., Pinkney J., Richstone D., TremaineS., 2000, ApJ, 539, L13Genzel R., Sch¨odel R., Ott T., Eisenhauer F., Hofmann R., LehnertM., Eckart A., Alexander T., Sternberg A., Lenzen R., Cl´enet Y.,Lacombe F., Rouan D., Renzini A., Tacconi-Garman L. E., 2003,ApJ, 594, 812Gnedin O. Y., 2003, ApJ, 589, 752Gonz´alez J. A., Hannam M., Sperhake U., Br¨ugmann B., Husa S.,2007, Physical Review Letters, 98, 231101Goodman J., Tan J. C., 2004, ApJ, 608, 108Graham A. W., 2004, ApJ, 613, L33Gualandris A., Merritt D., 2008, ApJ, 678, 780G¨ultekin K., Richstone D. O., Gebhardt K., Lauer T. R., TremaineS., Aller M. C., Bender R., Dressler A., Faber S. M., FilippenkoA. V., Green R., Ho L. C., Kormendy J., Magorrian J., PinkneyJ., Siopis C., 2009, ApJ, 698, 198Haehnelt M. G., Kau ff mann G., 2002, MNRAS, 336, L61Haiman Z., Loeb A., 1999, ApJ, 521, L9Heggie D. C., Mathieu R. D., 1986, in P. Hut & S. L. W. McMillaned., The Use of Supercomputers in Stellar Dynamics Vol. 267 ofLecture Notes in Physics, Berlin Springer Verlag, StandardisedUnits and Time Scales. pp 233– + Hein¨am¨aki P., 2001, A&A, 371, 795Herrmann F., Hinder I., Shoemaker D., Laguna P., 2007, Classicaland Quantum Gravity, 24, 33Ho ff man L., Loeb A., 2007, MNRAS, 377, 957Hopkins P. F., Hernquist L., 2010, MNRAS, 407, 447Hopkins P. F., Hernquist L., Cox T. J., Di Matteo T., RobertsonB., Springel V., 2006, ApJS, 163, 1Hopkins P. F., Hernquist L., Cox T. J., Robertson B., Krause E.,2007, ApJ, 669, 45Hut P., Rees M. J., 1992, MNRAS, 259, 27PIwasawa M., Funato Y., Makino J., 2006, ApJ, 651, 1059Janssen G. H., Stappers B. W., Kramer M., Purver M., JessnerA., Cognard I., 2008, in C. Bassa, Z. Wang, A. Cumming, &V. M. Kaspi ed., 40 Years of Pulsars: Millisecond Pulsars, Mag-netars and More Vol. 983 of American Institute of Physics Con-ference Series, European Pulsar Timing Array. pp 633–635Jenet F., Finn L. S., Lazio J., Lommen A., McLaughlin M., StairsI., Stinebring D., Verbiest J., Archibald A., Arzoumanian Z.,Backer D., Cordes J., Demorest P., Ferdman R., Freire P., other2009, ArXiv e-printsJiang C. Y., Jing Y. P., Faltenbacher A., Lin W. P., Li C., 2008,ApJ, 675, 1095Kau ff mann G., Haehnelt M., 2000, MNRAS, 311, 576Khan F., Just A., Merritt D., 2011, ArXiv e-printsKhochfar S., Burkert A., 2006, A&A, 445, 403Kocsis B., Yunes N., Loeb A., 2011, ArXiv e-printsKulkarni S. R., Hut P., McMillan S., 1993, Nature, 364, 421Lacey C., Cole S., 1993, MNRAS, 262, 627Lauer T. R., Faber S. M., Ajhar E. A., Grillmair C. J., ScowenP. A., 1998, AJ, 116, 2263Levan A. J., Tanvir N. R., Cenko S. B., Perley D. A., WiersemaK., Bloom J. S., Fruchter A. S., de Ugarte Postigo A., O’BrienP. T., Butler N., van der Horst A. J., Leloudas G., Morgan A. N., c (cid:13) , 000–000 Kulkarni & Loeb
Misra K., Bower G., Farihi J., et al., 2011, ArXiv e-printsLi Y., Hernquist L., Robertson B., Cox T. J., Hopkins P. F.,Springel V., Gao L., Di Matteo T., Zentner A. R., Jenkins A.,Yoshida N., 2007, ApJ, 665, 187Liu X., Shen Y., Strauss M. A., 2011, ArXiv e-printsLoeb A., 2010, Phys. Rev. D, 81, 047503Loeb A., Rasio F. A., 1994, ApJ, 432, 52Lousto C. O., Zlochower Y., 2009, Phys. Rev. D, 79, 064018Magorrian J., Tremaine S., Richstone D., Bender R., Bower G.,Dressler A., Faber S. M., Gebhardt K., Green R., Grillmair C.,Kormendy J., Lauer T., 1998, AJ, 115, 2285Makino J., 1997, ApJ, 478, 58Makino J., Aarseth S. J., 1992, PASJ, 44, 141Makino J., Ebisuzaki T., 1996, ApJ, 465, 527Makino J., Funato Y., 2004, ApJ, 602, 93Manchester R. N., 2008, in C. Bassa, Z. Wang, A. Cumming, &V. M. Kaspi ed., 40 Years of Pulsars: Millisecond Pulsars, Mag-netars and More Vol. 983 of American Institute of Physics Con-ference Series, The Parkes Pulsar Timing Array Project. pp 584–592Marconi A., Hunt L. K., 2003, ApJ, 589, L21McLure R. J., Dunlop J. S., 2002, MNRAS, 331, 795Menou K., Haiman Z., Narayanan V. K., 2001, ApJ, 558, 535Merritt D., 2000, in F. Combes, G. A. Mamon, & V. Charman-daris ed., Dynamics of Galaxies: from the Early Universe to thePresent Vol. 197 of Astronomical Society of the Pacific Confer-ence Series, Black Holes and Galaxy Evolution. pp 221– + Merritt D., 2006, ApJ, 648, 976Merritt D., Mikkola S., Szell A., 2007, ApJ, 671, 53Merritt D., Milosavljevi´c M., 2005, Living Reviews in Relativity,8, 8Merritt D., Poon M. Y., 2004, ApJ, 606, 788Merritt D., Szell A., 2006, ApJ, 648, 890Mikkola S., Aarseth S. J., 1990, Celestial Mechanics and Dynam-ical Astronomy, 47, 375Mikkola S., Valtonen M. J., 1990, ApJ, 348, 412Milosavljevi´c M., Merritt D., 2001, ApJ, 563, 34Milosavljevi´c M., Merritt D., 2003a, ApJ, 596, 860Milosavljevi´c M., Merritt D., 2003b, in J. M. Centrella ed., TheAstrophysics of Gravitational Wave Sources Vol. 686 of Amer-ican Institute of Physics Conference Series, The Final ParsecProblem. pp 201–210Miralda-Escud´e J., Gould A., 2000, ApJ, 545, 847Monaco P., Fontanot F., Ta ff oni G., 2007, MNRAS, 375, 1189Mortlock D. J., Warren S. J., Venemans B. P., Patel M., HewettP. C., McMahon R. G., Simpson C., Theuns T., Gonz´ales-SolaresE. A., Adamson A., Dye S., Hambly N. C., Hirst P., Irwin M. J.,Kuiper E., Lawrence A., R¨ottgering H. J. A., 2011, Nature, 474,616Peng C. Y., Impey C. D., Ho L. C., Barton E. J., Rix H.-W., 2006,ApJ, 640, 114Peres A., 1962, Physical Review, 128, 2471Peters P. C., 1964, Phys. Rev., 136, 1224Pretorius F., 2005, Physical Review Letters, 95, 121101Quinlan G. D., 1996, New Astron., 1, 35Quinlan G. D., Hernquist L., 1997, New Astron., 2, 533Richstone D., Ajhar E. A., Bender R., Bower G., Dressler A.,Faber S. M., Filippenko A. V., Gebhardt K., Green R., Ho L. C.,Kormendy J., Lauer T. R., Magorrian J., Tremaine S., 1998, Na-ture, 395, A14 + Saslaw W. C., Valtonen M. J., Aarseth S. J., 1974, ApJ, 190, 253Sesana A., Haardt F., Madau P., Volonteri M., 2004, ApJ, 611, 623 Shen Y., 2009, ApJ, 704, 89Sheth R. K., Tormen G., 1999, MNRAS, 308, 119Sigurdsson S., Hernquist L., 1993, Nature, 364, 423Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MN-RAS, 380, 877Springel V., White S. D. M., Jenkins A., Frenk C. S., Yoshida N.,Gao L., Navarro J., Thacker R., Croton D., Helly J., PeacockJ. A., Cole S., Thomas P., Couchman H., Evrard A., Colberg J.,Pearce F., 2005, Nature, 435, 629Stone N., Loeb A., 2011, MNRAS, 412, 75Ta ff oni G., Mayer L., Colpi M., Governato F., 2003, MNRAS,341, 434Tanaka T., Haiman Z., 2009, ApJ, 696, 1798Taylor J. E., Babul A., 2001, ApJ, 559, 716Tonry J. L., 1984, ApJ, 279, 13Tremaine S., Gebhardt K., Bender R., Bower G., Dressler A.,Faber S. M., Filippenko A. V., Green R., Grillmair C., Ho L. C.,Kormendy J., Lauer T. R., Magorrian J., Pinkney J., RichstoneD., 2002, ApJ, 574, 740Valtonen M. J., 1976, A&A, 46, 435Valtonen M. J., Mikkola S., Heinamaki P., Valtonen H., 1994,ApJS, 95, 69van den Bosch F. C., Lewis G. F., Lake G., Stadel J., 1999, ApJ,515, 50Volonteri M., 2007, ApJ, 663, L5Volonteri M., Haardt F., Madau P., 2003, ApJ, 582, 559Volonteri M., Rees M. J., 2006, ApJ, 650, 669Wetzel A. R., Cohn J. D., White M., 2009, MNRAS, 395, 1376Wyithe J. S. B., Loeb A., 2003a, ApJ, 590, 691Wyithe J. S. B., Loeb A., 2003b, ApJ, 595, 614Xu G., Ostriker J. P., 1994, ApJ, 437, 184Yoo J., Miralda-Escud´e J., 2004, ApJ, 614, L25Yu Q., 2002, MNRAS, 331, 935Yunes N., Coleman Miller M., Thornburg J., 2011, Phys. Rev. D,83, 044030Zauderer B. A., Berger E., Soderberg A. M., Loeb A., Narayan R.,Frail D. A., Petitpas G. R., Brunthaler A., Chornock R., Carpen-ter J. M., Pooley G. G., Mooley K., Kulkarni S. R., et al., 2011,ArXiv e-printsZentner A. R., Berlind A. A., Bullock J. S., Kravtsov A. V., Wech-sler R. H., 2005, ApJ, 624, 505 c (cid:13)000