Stellar-mass black holes in star clusters: implications for gravitational wave radiation
aa r X i v : . [ a s t r o - ph . S R ] O c t Mon. Not. R. Astron. Soc. , 1–11 (2009) Printed 12 October 2018 (MN L A TEX style file v2.2)
Stellar-mass black holes in star clusters: implications forgravitational wave radiation
Sambaran Banerjee , ⋆ , Holger Baumgardt and Pavel Kroupa Argelander-Institut f¨ur Astronomie, Auf dem H¨ugel 71, 53121, Bonn, Germany Alexander von Humboldt fellow
Accepted — — —. Received — — —; in original from — — —
ABSTRACT
We study the dynamics of stellar-mass black holes (BH) in star clusters with particu-lar attention to the formation of BH-BH binaries, which are interesting as sources ofgravitational waves (GW). In the present study, we examine the properties of theseBH-BH binaries through direct N-body simulations of star clusters using the NBODY6code on Graphical Processing Unit (GPU) platforms. We perform simulations for starclusters with low-mass stars starting from Plummer models with an initialpopulation of BHs, varying the cluster-mass and BH-retention fraction. Additionally,we do several calculations of star clusters confined within a reflective boundary mim-icking only the core of a massive star cluster which can be performed much fasterthan the corresponding full cluster integration. We find that stellar-mass BHs withmasses ∼ M ⊙ segregate rapidly ( ∼
100 Myr timescale) into the cluster core andform a dense sub-cluster of BHs within typically 0 . − . N & × , typically 1 - 2 BH-BHmergers occur per cluster within the first ∼ Key words: gravitational waves – stellar dynamics – scattering – methods: N-bodysimulations – galaxies: star clusters – black hole physics
Star clusters, e.g. , globular clusters (henceforth GC), youngand intermediate-age massive clusters and open clusters har-bor a large overdensity of compact stellar remnants com-pared to that in the field by virtue of their high density andthe mass segregation of the compact remnants towards thecluster core (Hut & Verbunt 1983). These compact stars, ⋆ E-mail: [email protected] (SB); [email protected] (BG); [email protected] (PK) which are neutron stars (henceforth NS) and black holes(henceforth BH), are produced by the stellar evolution ofthe most massive stars within the first ∼
50 Myr after clus-ter formation. These compact stars, being generally heav-ier than the remaining low-mass stars of the cluster, segre-gate quickly (within one or a few half-mass relaxation times)to the cluster core, forming a dense sub-cluster of compactstars. Such a compact-star sub-cluster is of broad interestas it efficiently produces compact-star binaries through dy-namical encounters (Hut & Verbunt 1983; Hut et al. 1992), e.g. , X-ray binaries and NS-NS and BH-BH binaries, which c (cid:13) S. Banerjee, H. Baumgardt and P. Kroupa are of interest for a wide range of physical phenomena. Forexample, X-ray binaries are the primary sources of GC X-rayflux, while mergers of tight NS-NS binaries (through emis-sion of gravitational waves) is the most likely scenario forthe production of short duration GRBs and both NS-NS andBH-BH mergers are very important sources of gravitationalwaves (henceforth GW) detectable by future gravitationalwave observatories like “Advanced LIGO” (AdLIGO) and“LISA” (Amaro-Seoane et al. 2007). Double-NS systems arealso interesting because they could be observable as double-pulsar systems (Ransom 2008), allowing important tests ofGeneral Relativity (Kramer 2008). In the present work, weinvestigate the dynamics of stellar-mass BHs in star clus-ters, with particular emphasis on dynamically formed BH-BH binaries. Such binaries are strong sources of GWs asthey spiral-in through GW radiation, a process detectableout to several thousand Mpc distances.Tight BH binaries that can merge within a Hub-ble time can also be produced in the galactic field fromtight stellar binaries, which result from stellar evolu-tion of the components ( e.g. , Bulik & Belczynski 2003;O’Shaughnessy et al. 2007, 2008; Belczynski et al. 2007).However, Belczynski et al. (2007) have shown with revisedbinary-evolution models that the majority of potential BH-BH binary progenitors actually merge because of commonenvelope (CE) evolution which occurs when any of the bi-nary members crosses the Hertzsprung gap. They found thatthis reduces the merger rate by a factor as large as ∼ ∼ − , subject to the uncertainties of the CE evolution modelthat has been incorporated. Note that the corresponding NS-NS merger rate, as the above authors predict, is considerablyhigher, ∼
20 yr − . In view of the above result, the majorityof merging BH binaries are possibly those that are formeddynamically in star clusters.As studied by several authors earlier (Merritt et al.2006; Mackey et al. 2007), black holes, formed through stel-lar evolution, segregate into the cluster core within ∼ . e.g. , Merritt et al. (2006); Mackey et al.(2007). Mackey et al. (2007) found notable agreement be-tween the core expansion as obtained from their N-body simulations with the observed age-core radius correlationfor star clusters in the Magellanic Clouds.The dynamics of stellar mass BHs and formationof BH binaries in star clusters and the resulting rateof GW emission from close enough BH binaries hasbeen studied by several authors using N-body integrations(Portegies Zwart & McMillan 2000), numerical 3-body scat-tering experiments (G¨ultekin et.al. 2004) or Monte-Carlomethods (O’Leary et.al 2006). Portegies Zwart & McMillan(2000) considered the rate of mergers of escaping BH bina-ries from various stellar systems, e.g. , massive GCs, youngpopulous clusters and galactic nuclei. They estimated themerger rates within 15 Gyr (their adopted age of the Uni-verse) by assuming the binding-energy ( E b ) distribution ofthe escaping BH binaries to be uniformly distributed inlog E b , as inferred from simulations of N ≈ i.e. , contribution from all types of starclusters) detection rate is negligible, it can be as high as ∼ − for AdLIGO. G¨ultekin et.al. (2004) performedsequential numerical integrations of BH binary-single BHclose encounters in a uniform stellar background, where, inbetween successive encounters, the BH-binaries were evolveddue to GW emission. From such simulations, these authorsstudied the growth of BHs through successive BH-binarymergers for the first time. In a more self-consistent studyof the growth of BHs and the possibility of formation of in-termediate mass black holes through successive BH-binarymergers, O’Leary et.al (2006) used the Monte-Carlo ap-proach and considered “pure” BH clusters that are dynam-ically detached from their parent clusters which can be ex-pected to form due to mass stratification instability (Spitzer1987). These authors considered dynamically formed BHbinaries from 3-body encounters in such BH clusters uti-lizing theoretical cross-sections of 3-body binary formation(see O’Leary et.al 2006 and references therein). Further-more, they included both binary-single and binary-binaryencounters. Considering the sub-set of BH binaries formedin the BH clusters that merge within a Hubble time, theydetermined typically a few AdLIGO detections per year forold GCs.While such results are remarkable and promising, amore detailed study of the dynamics of stellar-mass BHsin star clusters with a realistic number of stars is essen-tial. As the cross-sections of different processes governingthe dynamics of BH binaries, viz. , 3-body binary formation,binary-single star encounters and ionization have differentdependencies on the number of stars N of the cluster, an ex-trapolation to much different N can be problematic. Henceit is important to consider clusters with values of N appro-priate for massive clusters or GCs. Also, exact treatmentsof the various dynamical processes are crucial for realisticpredictions of BH-BH merger rates. Finally, a more care-ful study of the different dynamical processes leading to theformation and evolution of BH binaries in star clusters isneeded to understand better under which conditions tightinspiralling BH binaries can be formed dynamically from astar cluster.In the present work, we make a detailed study of thedynamics of BH-BH binaries formed in a BH sub-cluster, c (cid:13) , 1–11 tellar-mass black holes in star clusters Q t (Myr) Figure 1.
Typical evolution of virial coefficient Q of a star clusterwith a reflective sphere. Shown is the evolution of the run R3K180of Table 1. A rapid initial heating is observed, caused by super-elastic encounters between BHs, followed by a saturation of theheating curve, caused by the enhanced escape rate of stars due tothe heating (see text). as introduced above. In particular, we investigate whetherhard enough BH binaries that can merge via gravitationalradiation in a Hubble time within the cluster or after gettingejected from the cluster, can be formed in such a sub-cluster.To that end, we perform direct N-body integrations of con-centrated star clusters (half-mass radius r h N low-mass stars in which a certain number ofstellar-mass BHs is added, representing a star cluster withan evolved stellar population.The present paper is organized as follows: In Sec. 2 wedescribe our simulations in detail. We discuss the various ele-ments and assumptions of the simulations and summarize allthe runs that we perform (Sec. 2.2). We also discuss the useof a reflective boundary to simulate only the core of a starcluster (Sec. 2.3). In Sec. 3 we discuss our results with partic-ular emphasis on the dynamical BH binaries formed duringthe simulations. We discuss the BH-BH mergers that occurwithin the clusters and also the merger timescales of the BHbinaries that escape from the clusters (Sec. 3.2) and obtainthe distributions of merger-times for both cases (Sec. 3.3).Finally, in Sec. 4 we interpret our results in the context ofdifferent types of star clusters that are observed and provideestimates of BH-BH merger detection rates. To study the dynamics of BHs in star clusters, we performdirect N-body simulations using NBODY6 with star clustersin which a certain number of BHs are added initially. The ini-tial density distribution follows a Plummer model with half-mass radius r h N low-mass main-sequence stars in the mass-range 0 . M ⊙ m . M ⊙ . Ob-served half mass radii for GCs and open clusters are usuallylarger, but as we expect the clusters to expand considerablydue to the heating caused by the encounters in the BH core(see Sec. 1), we begin with more concentrated clusters.Black holes formed through stellar evolution are typ-ically within the mass-range 8 M ⊙ . M BH . M ⊙ for stellar progenitors with solar-like metallicity (Casares 2007;Belczynski et al. 2009). The exact form of the BH mass-function is still debated and depends on the metallicity ofthe parent stellar population and the stellar wind mass-lossmodel (Belczynski et al. 2009). In a dense stellar environ-ment it is further affected by the frequent stellar merg-ers and binary coalescence (Belczynski et al. 2008). In viewof the uncertainty of the BH mass-function, we consideronly equal-mass BHs, with M BH = 10 M ⊙ in this work.Such a value was also assumed by earlier authors, e.g. ,Portegies Zwart & McMillan (2000), Benacquista (2002).A specified number of BHs are added to each cluster,and we assume that the BHs follow the same spatial distribu-tion as the stars initially. The number of BHs added dependson the BH retention fraction. We explore both full retentionand the case where half of the BHs are ejected from the clus-ter by natal kicks. With such a cluster, we mimic the epochat which massive stars have already evolved and produceBHs. While such a cluster is not completely representativeof the more realistic case in which the BHs are producedfrom stellar evolution during the very early phases, it stillserves the purpose of studying the dynamics of the BHs sincethis becomes important only later after segregation of theBHs into the cluster core (see Sec. 3).We perform our simulations using NBODY6 code(Aarseth 2003) enabled for use with Graphical ProcessingUnits (GPU). We have made arrangements in NBODY6 toput a specified number N BH of BHs of a given mass M BH bypicking stars randomly from the cluster and replacing themby the BHs. Necessary rescalings have been done to compen-sate for the excess mass gained by the cluster. These runs areperformed on NVIDIA 9800 GX2/GTX 280/GTX 285 GPUplatforms, located at the Argelander-Institut f¨ur Astronomie(AIfA), University of Bonn, Germany, with the GPU en-abled version of NBODY6 (see Sverre Aarseth’s Homepage:sverre.com). To evolve the BH-BH binaries due to GW emission, Pe-ters’ formula (Peters 1964) is utilized in NBODY6 (also seeBaumgardt et al. 2006), which provides approximate orbit-averaged rates of evolution of binary semi-major axis a andeccentricity e due to GW emission. According to this for-mula, the merger time T mrg of an equal-mass BH-BH binarydue to GW emission is given by T mrg = 150Myr (cid:16) M ⊙ M BH (cid:17) (cid:18) aR ⊙ (cid:19) (1 − e ) / . (1)Peters’ formula is limited up to the mass quadrapoleterms of the radiating system. Close to the merger, higherorder PN terms become important, which modifies the abovevalue of the merger time (Blanchet 2006). In the presentwork however, we are primarily interested in an overallstatistics of the merger events rather than the detailed or-bit of a single BH-BH merger, so that the higher order PNcorrections are not crucial. Hence we restrict ourselves to Pe-ters’ formula. Moreover, the orbit of a tight BH-BH binarygets modified by dynamical encounters, which will anywaymodify its merger time to a much larger extent compared tothat for the unperturbed binary with the higher PN terms. c (cid:13) , 1–11 S. Banerjee, H. Baumgardt and P. Kroupa
In NBODY6, the orbital evolution of compact binariesis also considered when the binary is inside a hierarchy. Thusa tight BH binary will continue to shrink even if it acquiresan outer member forming a hierarchical triple, which canoften happen due to the strong focussing effect of BH bi-naries. Also, energy removal due to GW (bursts) during aclose hyperbolic passage between two BHs is considered inthe code.Numerical simulations of BH-BH mergers (see Hughes2009 for an excellent review) indicate that for unequal-massBHs or even for equal-mass BHs with unequal spins, themerged BH product acquires a velocity-kick of typically 100km s − or more due to an asymmetry in momentum out-flow from the system, associated with the GW emission.Although we consider equal-mass BHs, the merger kicks as-sociated with the inequality of the spins of the merging BHswould be generally sufficient to eject merged BHs from thecluster. Therefore, in our simulations we provide an arbi-trarily large kick of 150 km s − immediately after a BH-BHmerger, to make sure that the merged BH escapes. To study the rate of BH-BH mergers coming from a star clus-ter, we perform simulations of isolated star clusters with sin-gle low-mass stars and BHs as mentioned above. The forma-tion of the BH-core through mass segregation and its dynam-ics remain largely unaffected by the presence of a tidal field,which mainly affects stars near the tidal boundary. Whilethe enhanced removal of stars accelerates the core-collapse ofthe cluster (see e.g. , Spitzer 1987, Ch. 3), the latter is morestrongly enhanced by the collapse of the BHs themselves( ∼
100 Myr timescale, see below), so that the effect of anytidal field is only second-order. Hence, isolated clusters aregood enough for our purposes. Further, for simplicity, we donot take into account primordial binaries in this initial study.The primordial binary fraction in GCs and their perioddistribution is still widely debated (Ashman & Zepf 1998;Bellazzini et al. 2002; Sommariva et al. 2009). The presenceof a primordial binary population can however significantlyinfluence the dynamics of stellar-mass BHs which we deferfor a future more detailed study.For solar-like metallicity, Eggleton’s stellar evolu-tion model (Eggleton, Fitchett & Tout (1989), adapted inNBODY6) gives about N BH ≈
200 BHs for a cluster with N = 10 stars following a Kroupa IMF (Kroupa 2001). Theabove N BH (or its proportion with N ) is thus an upper limitto the number of BHs in a GC that corresponds to a full re-tention ( i.e. , no or low natal kicks for all BHs).We perform 2 simulations with N = 4 . × and N BH = 80, 2 runs with N = 6 . × , N BH = 110, i.e. ,about full BH-retention. Two of the above runs are repeatedwith half the above N BH s. Also, one run with N = 5 × and excess N BH = 200, appropriate for a top-heavy IMF hasbeen performed. Finally, we do 2 runs with N = 10 with N BH = 80 (about 50% retention fraction) and 200 (full re-tention). All the clusters consist of low-mass stars between0 . M ⊙ . m . . M ⊙ with a Kroupa IMF and the BHs have M BH = 10 M ⊙ , as discussed in detail in the beginning of thesection.In addition to these systems, we also study stellar-massBHs in clusters with smaller N representing open clusters, R ( p c ) t (Myr) 0 10 20 30 40 50 60 70 80 90 100 1 10 100 1000 N B H ( t ) t (Myr) Figure 2.
Typical example of the mass segregation of BHs.Shown is the radial position R vs. time t (top panel) of all BHs formodel C50K80 of Table 1. Each of the points represents a BH.The BHs segregate within 50 Myr during which N BH remainsunchanged (bottom panel). As the BH sub-cluster becomes denseenough that BH-BH binary formation takes place through 3-bodyencounters, BHs and BH binaries are ejected from the BH-coreand N BH starts to decrease. in order to estimate the lower-limit in cluster mass for theoccurrence of BH-BH mergers. We perform 10 runs with N = 5 × , N = 10 and N = 2 . × each with full BH-retention ( i.e. , N BH = 12, 20 and 50 respectively). Resultsof all our runs are summarized in Table 1. We also perform simulations with a smaller number of starsand BHs that are confined within a reflecting sphericalboundary. With such a dynamical system one can mimicthe core of a massive cluster, where the BHs are concen-trated after mass segregation. The advantage of this ap-proach is that one can simulate the evolution of a massivecluster with much fewer stars. We simulate N = 3000 − . ∼ M ⊙ pc − , appropriate for the core-density of amassive cluster. The initial BH population is taken to be N BH ≈
100 or 200, representing half or full BH-retentionrespectively, of a N = 10 star cluster. In this way, the sim- c (cid:13) , 1–11 tellar-mass black holes in star clusters ulation of the BH-core of a massive GC can be performedmuch faster than the simulation of the whole cluster.To implement the reflective boundary, we simply con-sider all outgoing stars beyond the reflective sphere R s within each block and reflect them elastically. This is donein the main integration loop. For these stars, we computethe force polynomials (Aarseth 2003) separately (using theCPU) which makes the code slower by factors of 2-3 depend-ing on the stellar density. While the behavior of a N-bodycode in presence of a reflective boundary needs to be stud-ied in more detail, the present implementation seems to befairly stable with single-stars and BHs, with energy errorsmarginally larger than that for a free cluster. (In fact, in theNBODY6 package, some routines are already available foroptionally implementing the reflective boundary). For starsand BHs beyond a pre-assigned speed v esc , representing theescape speed of the host cluster, we do allow them to escapethrough the reflective boundary.It is of course important to note that a cluster sim-ulation within a reflective boundary cannot represent theevolution of a free cluster, in particular, since the cluster ex-pands due to the BH-heating which continually reduces thestellar density in the core. With the reflective boundary, thestellar density still decreases because of the escape of starswith speed v > v esc , but it generally does not mimic thecluster expansion. The cooling effect of these escaping starsalso helps to inhibit the runaway heating of the cluster dueto the binding energy released by binary-star/binary-binaryencounters. The escaping single BHs and BH binaries areparticularly efficient in this mass-loss cooling due to theirlarger masses. Fig. 1 shows a typical example from one ofour runs (R3K180 of Table 1, see below) in which the be-havior of the virial coefficient Q ≡ T /V ( T = total kineticenergy of the system, V = total potential energy) is shown.While Q initially increases quickly indicating rapid heatingof the cluster, the escape rate of stars and BHs also increasesdue to the increased velocity dispersion. The latter effect inturn increases the cooling rate so that Q finally tends toflatten to a value for which the dynamical heating is bal-anced by the escape-cooling of the cluster. For the reflectiveboundary clusters used in our simulations (see above), theinitial velocity-dispersion of the stars within the reflectiveboundary is chosen to be v ≈ − and escape speed v esc ≈
24 km s − . Typically, during each run, the virial co-efficient Q increases by a factor of about 6.5 and most of thegrowth takes place during the first few hundred Myr, as inFig. 1. As v ∼ Q / , it increases to ≈
17 km s − , which istypical for the central velocity dispersion of a massive clus-ter. In a reflective boundary cluster, the BH binaries startforming immediately after the start of the integration sincethe system already initiates with a BH-core. We perform aset of 6 runs with reflective boundary clusters (see Table 1).We shall discuss results of our reflective boundary simula-tions in Sec. 3. In this section we discuss the results of the simulations in-troduced in Sec. 2.2. As already mentioned, we focus on tight BH-BH binaries both within the cluster and the es-caped ones that can merge within a few Gyr. Table 1 pro-vides an overview of the results of our simulations for bothisolated and reflective clusters. It shows the number of BH-BH mergers within the cluster for each computation andthe corresponding merger times. Also, for each computation,the numbers of escaped BH-BH binaries that merge within3 Gyr are shown. For convenience of discussion, we utilizeone of these models, viz. , the one with identity C50K80 (seeTable 1) — all others generally posses similar properties.Fig. 2 (top panel) demonstrates the mass segregationof the BHs in the cluster which takes about 50 Myr. Beforethe BHs segregate to within about 0 . N BH remains constant. Asthe BHs segregate within about 0 . N BH during thisphase is also shown in Fig. 2 (bottom panel).The time t seg taken by the BHs to segregate to thecluster core to form the BH sub-cluster is essentially thecore-collapse time for the initial BH cluster. t seg is given by(see e.g. , Spitzer (1987), Ch. 3): t seg = h m i M BH t cc , (2)where, t cc is the core-collapse time of the host star cluster it-self ( i.e. , without the BHs) and h m i is its mean stellar mass.According to numerical experiments by several authors ( e.g. ,Baumgardt, Hut & Heggie (2002)), a Plummer cluster takesabout t cc ≈ t rh (0) to reach the core-collapse stage, where t rh (0) is the initial half-mass relaxation time of the Plum-mer cluster. For the above example, t rh (0) ≈
73 Myr (seeEqn. (2-63) of Spitzer (1987)) and h m i ≈ . t seg ≈
65 Myr from Eqn. (2). This roughly agrees with theformation time of the central dense BH-cluster as observedin the N-body integration (see Fig. 2).
A particular dynamically formed BH-BH binary can initiallybe very wide, sometimes with semi-major-axis a more than5000 R ⊙ . However, it then shrinks very rapidly due to re-peated super-elastic encounters with the surrounding BHs.The collisional hardening rate is given by ˙ a C ∝ a (seeBanerjee & Ghosh 2006 and references therein). The vari-ation in a and e is random due to the stochastic natureof the encounters, but the binary hardens on average (see e.g. , Heggie & Hut 2003, Ch. 19 & Ch. 21) and the overallhardening rate decreases with a . These encounters includeboth flybys and exchanges with single BHs. When the binaryis wide, and the recoils it receives are low, so that it usu-ally remains in the cluster. However, as it becomes harder,the recoils become stronger, and the binary is ejected from c (cid:13) , 1–11 S. Banerjee, H. Baumgardt and P. Kroupa
100 1000 0 100 200 300 400 500 600 700 a ( R O · ) t (Myr) Figure 3.
Semi-major-axis a vs. evolution time t for all BH-BH pairs formed for the model C50K80, where a different symbol is usedfor each newly formed BH-BH binary. The sequences of the points indicate that each BH binary initially hardens quickly to smallervalues of a where the hardening rate decreases. The estimated upper-limit a esc for which a BH-BH pair can escape due to recoil froman encounter with a single BH is indicated by the horizontal line. Hardening rates of all BH pairs terminate below this line. Note thatdifferent traces of the same symbol represent different hardening sequences following encounters. the BH-core, but returns to the BH-core due to dynami-cal friction (see Sec. 1). Finally, the recoil becomes strongenough that the BH binary escapes from the cluster (see alsoBenacquista 2002).For equal mass binary-single star encounters, 40%of the binding energy of the binary is released perclose encounter on average (Spitzer 1987, Eqn. (6-25);Hut, McMillan & Romani 1992). This leads to the follow-ing expression for the average recoil velocity of a BH binarydue to encounters with single BHs: h v rec i = 6 . × − GM BH a . (3)Hence, the BH binary can be ejected due to recoil from acluster with escape velocity v esc for a < a esc , where a esc isobtained by setting the left hand side of Eqn. (3) to v esc .For the above Plummer cluster, we then get from the valueof its central potential (see Heggie & Hut 2003, Table 8.1) a esc ≈ R ⊙ , which is shown as solid line in Fig. 3. Theescaped binaries in our calculations are generally found tobe harder than a esc except for a few systems, which escapewithin triple-BHs.Fig. 3 provides an overall impression of the hardeningof the BH binaries where the values of a for all BH-BHbinaries formed during the computation are plotted. Eachnewly formed BH pair is represented by a different symbol.From the sequences of points, it can be seen that each of theBH-BH pairs exhibits an overall tendency of rapid hardeningand the steep a ( t ) dependence tends to flatten as a decreases.The termination of a particular curve generally indicatesescape of the corresponding BH binary from the cluster (ormore rarely a merger, see Sec. 3.2). Typically, a BH binaryleaves the cluster with a < R ⊙ (see Fig. 5). Can a BH-pair be hardened to small enough a (and eccentricenough at the same time) for it to merge via GW radiation?To investigate this question, we consider the positions of theBH-BH binaries within the cluster in a a vs . (1 − e ) plane(Fig. 4), where each newly formed BH pair is representedby a different symbol as in Fig. 3. The points are color-coded with the evolution time in Myr, the color-scale beingdisplayed on the right. Although for each BH-pair a and e fluctuate over the plane, the changes occur on a collisiontime-scale which is ∼ Myr. Since the orbital period of thebinaries corresponding to these points is much shorter (from ∼ days to years), these points generally represent binarieswhich are stable over many orbits. Overplotted on Fig. 4are lines of constant GW merger time, T mrg , as given byEqn. (1). While most of the points have very large mergertime, a few of them do lie around the T mrg = 10 Myr line.This indicates that these binaries are indeed hardened up tosmall enough a and/or acquire sufficient eccentricity that ifthey are left unperturbed, they can merge via GW emissionwithin several Myr. This feature is found to be generally truefor all the computations reported here (see Table 1), ie, eachof them does produce BH-pairs that are capable of mergingwithin several Myr. It is important to note that the shortmerger timescale of these binaries is due to their high eccen-tricities since they are generally far too wide ( ∼ R ⊙ ) tomerge unless they are highly eccentric (see Eqn. (1)). Thisis also true for the escaping BH binaries (see below).However, these merging BH-pairs, particularly due totheir large eccentricity and focussing effect, can still be per-turbed by further encounters on timescales ∼ Myr whichcan often prevent them from merging. In the particular ex- c (cid:13) , 1–11 tellar-mass black holes in star clusters O · )1-e T mrg = 10 MyrT mrg = 13 Gyr 10 100 1000 0.001 0.01 0.1 1 Figure 4.
Positions of all the BH binaries within the cluster in a 1 − e vs. a plane where different symbols are used to distinguishbetween different BH pairs (for C50K80). The color-coding of the points, as indicated by the color-scale, represents the evolution time(in Myr) of the cluster at which they appear at a particular location in the above plane. See text for details. ample shown in Fig. 4, only one of the “merging candi-dates” ( T mrg
10 Myr) could actually merge within thecluster (inverted triangles). For all the other candidates, in-cluding the one which escaped (filled dots), the eccentricitybecame much smaller due to encounters so that the mergertimescales increased considerably.As for the escapers, since they remain unperturbed af-terwards, all with GW merger times smaller than a Hub-ble time are of interest. However, we find that typically foreach cluster there are relatively few BH binaries with GWmerger times T mrg & e.g. , Fig. 5 and Fig. 6). In the aboveexample, i.e. , model C50K80, there is one BH-BH escaperwith T mrg ∼ ≈ ≈ <
10 Myr.The results of the present computations, as summarizedin Table 1, indicate that a medium-mass to massive clusterwith sufficient BH-retention is likely to have at least one BH-pair that can merge within the cluster within a time-span ofa few Gyr. Each of these clusters also typically eject a fewBH binaries that can merge within a Hubble time. In one ofthe models (C50K200), one BH-BH hyperbolic collision alsooccurred at t mrg ≈
820 Myr. As such hyperbolic mergersdid not show up in any of the other models, they must bemuch rarer than GW driven mergers.
As a representation of the merger time distribution of merg-ers within the cluster, a combined distribution of t mrg of allthe BH binaries that merged within the clusters is shown inFig. 6 (top panel). This is justified, as the values of BH-BH binary merger times t mrg within the cluster, shown in Ta-ble 1, apparently do not indicate any correlation of t mrg withthe cluster parameters, in particular the number of stars N .Both the isolated clusters and the reflective boundary clus-ters are included. The t mrg distribution in Fig. 6 indicatesthat N mrg decreases with t mrg . This can be expected, as atlater times, N BH in a cluster decreases, so the hardeningrate of BH binaries also decreases, making their merger lessprobable. Comparing our results for clusters with differentvalues of N , we find that the rate of depletion of BHs fromthe BH core during the first few Gyr does not depend appre-ciably on N . In our computations, BHs are depleted withina few Gyr for most of the models. Clusters with large N BH s, e.g. , C100K200, C50K200 do indicate longer BH-retention,but this is due to the expansion of the cluster (and thereforethe BH sub-cluster itself) because of the heating by the BH-core (see Sec. 1) and hence it is unlikely that tight BH-BHbinaries can be formed during the later phases.Fig. 6 (bottom panel) shows the distribution of themerger times for the escaping BH binaries from all the com-putations of Table 1. Note that the number of tight esca-pers in Table 1 also does not indicate any dependence on N . For each escaper t mrg includes its time of escape t esc , i.e. , t mrg = t esc + T mrg , where T mrg is calculated fromEqn. (1) using the values of a and e with which it escaped.The merging rate among escaped binaries also shows a de-crease with time. Such a decrease is expected from Eqn. (1).For example, if we consider that the escapers that mergewithin a Hubble time have a mean radius a ∼ R ⊙ and athermal eccentricity distribution dN mrg /de ∝ e , it is easyto show from Eqn. (1) that dN mrg /dT mrg ∝ T − / mrg . Theabove merger-time distribution indicate that significantlymore BH-BH mergers occur outside the clusters, i.e. , amongthe escapers, than within the clusters. We shall return to this c (cid:13) , 1–11 S. Banerjee, H. Baumgardt and P. Kroupa
Table 1.
Summary of the computations performed for isolated clusters and those with reflective boundary, see Sec. 2.2 & Sec. 2.3. Themeaning of different columns is as follows: Col. (1): Identity of the particular model — similar values with different names (ending withA,B etc) imply computations repeated with different random seeds, Col. (2): Total number of stars N , Col. (3): Number of simulations N sim with the particular cluster, Col. (4): Initial half-mass radius of the cluster r h (0) (isolated cluster) or radius of reflective sphere R s , Col. (5): Initial number of BHs N BH (0), Col. (6): Total number of BH-BH binary mergers within the cluster N mrg , Col. (7): Thetimes t mrg corresponding to the mergers, Col. (8): Number of escaped BH-pairs N esc — the three values of N esc are those with T mrg . R AdLIGO detected by AdLIGO assuming that the correspondingmodel cluster has a space density of ρ cl = 3 . h Mpc − (see Sec. 4.1).Model name N N sim r h (0) or R s (pc) N BH (0) N mrg t mrg (Myr) N esc R AdLIGO
Isolated clustersC5K12 5000 10 1.0 12 0 — — — — — —C10K20 10000 10 1.0 20 0 — — — — — —C25K50 25000 10 1.0 50 0 — — 3 1 1 — —C50K80 45000 1 1.0 80 1 698.3 3 1 0 28( ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± - e a (R O · )P (days) T mrg = 10 MyrT mrg = 100 MyrT mrg = 1 GyrT mrg = 3 GyrT mrg = 13 Gyr Figure 5.
Positions of all escaped BH binaries in a 1 − e vs. a plane for the model C65K110. Lines of constant merger times areplotted as in Fig. 4. point in Sec. 4. The median merger time for BH-mergerswithin the clusters is t mrg ≈ ≈ Our present study indicates that centrally concentrated starclusters, with N & . × are capable of dynamicallyproducing BH binaries that can merge within a few Gyr,provided a significant number of BHs are retained in theclusters after their birth. The results of our simulations (seeTable 1) imply that most of the BH-BH mergers occur withinthe first few Gyr of cluster evolution for both mergers withinthe cluster and mergers of escaped BH binaries.The above results imply that an important class of can-didates for dynamically forming BH binaries that mergeat the present epoch are star clusters with initial mass M cl & × M ⊙ , which are less than few Gyr old.Such clusters represent intermediate-age massive clusters(hereafter IMC) with initial masses close to the upper-limit of the initial cluster mass function (ICMF) in spiral(Wiedner, Kroupa & Larsen 2004; Larsen 2009a) and star-burst (Gieles et al. 2006) galaxies. While it is not impossi-ble to obtain BH-BH mergers within a Hubble time fromlower-mass clusters, the overall BH-BH merger and escape To correlate with the observed clusters, we consider the mass M cl of the parent cluster, i.e. , the cluster mass before the BHsare formed through stellar evolution, which are heavier than theclusters of low mass stars that we model, for the same value of N . For a given simulated cluster, we estimate the correspondingparent mass M cl by weighting the mean stellar mass of h m i cl ≈ . M ⊙ of a Kroupa IMF with star from 0 . M ⊙ upto 100 M ⊙ withthe total number of stars N for that cluster.c (cid:13) , 1–11 tellar-mass black holes in star clusters N m r g t mrg (Myr) N m r g t mrg (Myr) Figure 6. Top:
Distribution of the merger times t mrg for BHbinary mergers within the cluster for the models of Table 1. Bot-tom:
Distribution of the merger times t mrg for escaped BH bi-naries for the models of Table 1 (see text). rates strongly decrease with cluster-mass as Table 1 indi-cates. For M cl . . × M ⊙ , mergers already becomemuch rarer (see Table 1). Due to the statistical nature ofmerger or ejection events it is ambiguous to set any well-defined limit on the cluster-mass beyond which these eventsbecome appreciable (such an estimate would also require amuch larger number of N-body integrations). In view of ourresults, M cl ≈ × M ⊙ is a representative lower limit be-yond which an appreciable number of mergers and escapersmerging within a Hubble time can be obtained.Old globular clusters, which can be about 10 times ormore massive, are expected to produce mergers or escapersmore efficiently. As the timescale of depletion of BHs fromthe BH cluster is nearly independent of the parent clustermass (see Sec. 3.3), GCs can also be expected to produceBH-BH mergers over similar time-span as the IMCs, i.e. ,within the first few Gyr of evolution. Since GCs are typ-ically much older ( ∼
10 Gyr), they do not contribute sig-nificantly to the present-day merger rate, since most of themergers from them would have occurred earlier. Consider-ing the light-travel time of ≈ . D ≈ We now make an estimate of the BH-BH merger detec-tion rate from IMCs by ground-based GW observatorieslike LIGO and AdLIGO. In estimating the overall BH-BHmerger rate using the results of our model clusters, one needsto consider the distribution of the cluster parameters thatare varied over the models, viz. , cluster mass, half-mass ra-dius, and BH retention fraction. Such distributions are farfrom being well determined, except for the mass distributionfor young clusters in spiral and starburst galaxies (Bik et al.2003; Bastian & Lamers 2003; Gieles 2004). Therefore, de-termination of an overall merger rate considering the distri-bution of our computed clusters can be ambiguous. Hence,as a useful alternative, we determine the BH binary mergerdetection rates for each of the cluster models in Table 1 thatgives an appreciable number of mergers, for a representativedensity of IMCs. Such an approach has been considered byearlier authors, e.g. , O’Leary et.al (2006) and can provide areasonable idea of the rate of detection of BH-BH mergersfrom IMCs.As an estimate of the space density of IMCs,we adopt that for young populous clusters inPortegies Zwart & McMillan (2000), which has a simi-lar mass-range as the IMCs: ρ cl = 3 . h Mpc − , (4)where h is the Hubble parameter, defined as H /
100 km s − , H being the Hubble constant (Peebles 1993). The abovespace density has been derived from the space densitiesof spiral, blue elliptical and starburst galaxies (Heyl et al.1997) assuming that young populous clusters have the samespecific frequencies ( S N ) as old GCs (van den Bergh 1995;McLaughlin 1999), but in absence of any firm determinationof the S N s of the former. We compute the detection rate foreach model cluster assuming that it has a space density ofthe above value.The LIGO/AdLIGO detection rate of BH-BH merg-ers from a particular model cluster can be estimated from(Belczynski et al. 2007 and references therein) R LIGO = 43 πD ρ cl R mrg , (5)where R mrg is the compact binary merger rate from a clusterand D is the maximum distance from which the emitted GWfrom a compact-binary inspiral can be detected. D is givenby D = D (cid:18) M ch M ch,nsns (cid:19) / , (6)where D = 18 . M ch is the “chirp mass” of the com- c (cid:13) , 1–11 S. Banerjee, H. Baumgardt and P. Kroupa pact binary with component masses m and m , which isgiven by M ch = ( m m ) / ( m + m ) / , (7)and M ch,nsns = 1 . M ⊙ is that for a binary with two 1 . M ⊙ neutron stars.For a BH binary with m = m = 10 M ⊙ , M ch =8 . M ⊙ which gives D ≈ R AdLIGO (mean over 3 Gyr takinginto account the time of escape t esc of the escaped binaries)for the model clusters are shown in Table 1, where the cur-rently accepted value of the Hubble parameter h = 0 .
73 isassumed. The error in each detection rate is simply obtainedfrom the Poisson dispersion of the total number of mergersfor the corresponding cluster. Note that these detection ratesare for clusters with solar-like metallicity which is impliedby our assumption of 10 M ⊙ BHs for all the clusters.To obtain a basic estimate of the overall detection rateof BH-BH mergers from IMCs, we consider the subset ofour computed models that are isolated clusters with fullBH retention (see Table 1). We take the mass function ofthe IMCs to be a power law with index α = − e.g. , Gieles 2004; Larsen 2009a). Then theweighted average of the corresponding AdLIGO detectionrates is R AdLIGO ≈ ±
7) yr − , which estimates the totalpresent-day detection rate of BH-BH mergers from IMCs ex-pected for AdLIGO. The corresponding LIGO detection rateis negligible, R LIGO ≈ . × − yr − . Note that these BH-BH detection rates are only lower limits. First, the observedpopulation of star clusters can be an underestimation by afactor of 2 owing to their dissolution in the tidal field of theirhost galaxies (see Portegies Zwart & McMillan 2000 and ref-erences therein). Second, the above detection rates are onlyfrom IMCs and there can be additional contributions fromGCs (see above).In comparing the AdLIGO detection rates from ourcomputations with those from earlier works, we note thatour rates are typically an order of magnitude smaller thanthose of Portegies Zwart & McMillan (2000), but aboutone order of magnitude larger than those obtained byO’Leary et.al (2006). The principal origin of the former dif-ference is due to the fact that while we considered onlyIMCs distinguishing them as the most appropriate can-didates for producing present-day BH-BH mergers (seeabove), Portegies Zwart & McMillan (2000) also includedGCs, which have considerably larger spatial density and alsolarger fraction of BH-BH binaries merging within the Hubbletime, as obtained from their analytic extrapolations. Notethat the number of escapers as obtained by them for youngpopulous clusters and the fraction of them merging withina Hubble time (see Table 1 of Portegies Zwart & McMillan2000) is similar to that obtained from the present compu-tations, implying qualitative agreement. The above authorsapparently did not consider the time scales of formation anddepletion of the BH sub-systems in their preliminary study.On the other hand, although O’Leary et.al (2006) consid-ered clusters significantly more massive than our’s in theirMonte-Carlo approach, so that larger merger detection ratescan be expected, their clusters were much older (8 Gyr and 13 Gyr) than IMCs which, in accordance with our results,accounts for their much lower detection rate.It is interesting to note that the dynamical BH-BHmerger detection rates obtained by us are typically an orderof magnitude higher than that from primordial stellar bina-ries as predicted by Belczynski et al. (2007) based on theirrevised binary evolution model, and is similar to that for theisolated NS-NS binaries derived by them. Hence our resultsimply that dynamical BH-BH binaries constitute the domi-nant contribution to the BH-BH merger detection. Thus, thedynamical BH-BH inspirals from star clusters seem to be apromising channel for GW detection by the future AdLIGO,although their estimated detection rate with the presentLIGO detector is negligible, in agreement with the hithertonon-detection of GW. The work presented here is a first step towards a detailedstudy of the dynamics of stellar-mass BHs in star clustersand the consequences for GW-driven BH mergers, and im-provements in several directions are possible. First, we donot consider the initial phase of the cluster here when BHsform through stellar evolution and a more consistent ap-proach would be to begin with a star cluster with a full stel-lar spectrum and produce the BHs from stellar evolution.Note however that in the present study we have insertednumbers of BHs in our clusters of low-mass stars similar towhat would have been formed from the stellar evolution of acluster following a Kroupa IMF (see Sec. 2.2). Stellar evolu-tion also produces NSs (about twice in number as the BHs)which also segregate to the central region of the cluster. Itis interesting to study the dynamics of the NS cluster andhow it is affected by the (more concentrated) BH cluster —in particular the formation of tight inspiralling NS-NS andNS-BH binaries, which are important for both GW detectionand GRBs.Another aspect that we do not consider in our presentstudy is the effect of primordial stellar binaries. Tight stel-lar binaries aid the formation of compact binaries throughdouble-exchanges (Grindlay et.al. 2006), in addition to the3-body mechanism (see Sec. 1) which can increase the num-ber BH-BH (also BH-NS and NS-NS) binaries formed andhence the merger rate. Thus the study of BH-BH binaries instar clusters with primordial binaries is an important nextstep. Such studies are in progress and will be presented infuture papers.Finally, the number of N-body computations in this ini-tial study is not enough to obtain the BH-BH merger rate asa function of cluster mass and BH retention fraction with areasonable accuracy. There are typically one integration percluster model. To obtain merger rate dependencies with clus-ter parameters, e.g. , mass, half mass radius, binary fractionand BH retention, many computations are needed withinsmaller intervals, which involves much larger time and com-puting capacity than that utilized for the present project.Such results, combined with improved knowledge of clusterparameter distributions and BH formation through super-nova explosions of massive stars, would provide a robustestimate of the BH-BH merger detection rate. Conversely,with such a merger rate function, the detection of BH-BHmergers by GW detectors like AdLIGO in the near future c (cid:13) , 1–11 tellar-mass black holes in star clusters will shed light on the above mentioned long-standing ques-tions. ACKNOWLEDGMENTS
We are thankful to Sverre Aarseth of the Institute of As-tronomy, Cambridge, UK, for many discussions and sugges-tions and his timely modifications of the NBODY6 code,which have been very helpful for our computations. Wethank Keigo Nitadori for his development and regular up-dates of the GPU libraries utilized in NBODY6. We thankthe anonymous referee for very helpful comments and sug-gestions which have improved several parts of this paper sig-nificantly. This work has been supported by the Alexandervon Humboldt Foundation.
REFERENCES c (cid:13)000